###### Abstract

Supermassive black holes lying in the center of galaxies can launch
relativistic jets of plasma along their polar axis. The physics of black-hole
jets is a very active research topic in astrophysics, owing to the fact that
many questions remain open on the physical mechanisms of jet launching, of
particle acceleration in the jet, and on the radiative processes. In this work
I focus on the last item, and present a review of the current understanding of
radiative emission processes in supermassive-black-hole jets.

###### keywords:

relativistic Jets; active galactic nuclei; theoretical emission models

xx 1 5 Received: date; Accepted: date; Published: date yes Leptonic and
Hadronic Radiative Processes in Supermassive-Black-Hole Jets Matteo Cerruti
Matteo Cerruti

##  1 Introduction

###  1.1 Active Galactic Nuclei

The brightest, persistent objects in the Universe are active galactic nuclei
(AGNs), compact regions in the center of galaxies that can outshine the host
itself (composed of hundreds of billions of stars), and can be detected on
cosmological distances, up to a redshift $z=7.5$ (Bañados et al., 2018), which
means a comoving distance of about 9 Gpc in the current $\Lambda$CDM
cosmological model (Planck Collaboration et al., 2018). About 75 years of
observations at all wavelengths have collectively shaped the so-called AGN
unified model: the engine powering the system is a supermassive black hole
(SMBH) of $10^{8-9}\ M_{\odot}$ (Event Horizon Telescope Collaboration et al.,
2019), which accretes matter in the form of an accretion disk (Abramowicz and
Fragile, 2013) and which is surrounded at the parsec scale by an obscuring
dusty torus (Jaffe et al., 1993). The thermal emission from the disk ionizes
clouds of matter orbiting the SMBH, which then re-emit emission lines in
optical/UV (the broad and narrow line-regions, BLR/NLR, the first one close to
the black hole at the sub-parsec scale (Greene and Ho, 2005; Bentz et al.,
2009), the second one at much larger distances, hundreds of parsecs and more
(Bennert et al., 2002), beyond the torus, see a schematic representation in
Figure 1). The inclination of the system with respect to the observer is
responsible for the variety of AGN types, due to the fact that when the system
is seen face-on the observer has direct access to the SMBH, while when the
system is seen edge-on the inner regions are obscured by the torus (Netzer,
2015). Quasars (quasi-stellar objects) are now understood as luminous AGNs
seen face-on. The host galaxy of the first identified quasars (such as 3C 273
(Schmidt, 1963)) can now be resolved and studied with current instrumentation
(Bahcall et al., 1997).

A major dichotomy exists among AGNs: a minority of around 10$\%$ of AGNs are
bright in the radio band (Wilson and Colbert, 1995) and are hence dubbed
radio-loud AGNs, in contrast with the more common radio-quiet AGNs. The radio
non-thermal emission is associated with a pair of relativistic jets of plasma
which are launched along the polar axis of the SMBH and perpendicularly to the
accretion disk. These jets can travel up to the Mpc scale, exceeding the
(visible) size of the AGN host galaxy. For the nearest and brightest objects,
the emission from the relativistic jet has been resolved not only in the radio
band (Kellermann et al., 1998), but also in optical (Butcher et al., 1980),
X-rays (Marshall et al., 2005), and more recently in $\gamma$-rays (Abdo et
al., 2010; De Naurois, 2019). AGNs with extended radio emission from the jets
are called radio-galaxies. In this case the observer sees the jet from the
side. When the relativistic jet points instead in the direction of the
observer, relativistic effects boost the emission (Blandford and Rees, 1978;
Blandford and Königl, 1979), making these AGNs particularly luminous in the
Universe. These radio-loud AGNs observed down-the-jet are called blazars (Urry
and Padovani, 1995), and are characterized by high luminosity, rapid
variability, and high degree of polarization (Angel and Stockman, 1980).
Although blazars are a minority among AGNs their high luminosity makes them
ubiquitous in all extragalactic surveys at all wavelengths. It is remarkable
that in $\gamma$-rays, where the radio-quiet AGNs do not emit anymore, around
$70\%$ of all the known extragalactic sources are blazars (The Fermi-LAT
collaboration, 2019).

Figure 1: Schematic representation (not in scale) of the central regions of a
radio-loud AGN.

###  1.2 Blazars

Blazars come in two flavors: BL Lacertae objects, whose optical/UV spectrum is
characterized by a featureless non-thermal continuum (Shaw et al., 2013), and
Flat-Spectrum-Radio-Quasars (FSRQs), which show instead broad emission lines
from the BLR typical of quasars (Shaw et al., 2012). This dichotomy is also
observed in the radio-galaxies parent population, which are divided into
Fanaroff–Riley type I radio-galaxies, with strong core emission and relatively
fainter jets, and Fanaroff–Riley type II radio-galaxies, with much stronger
jets that terminate in a shock with the inter-galactic medium (Fanaroff and
Riley, 1974).

The fact that the relativistic jet of plasma moves towards the observer with
speed $v=\beta c$, Lorentz factor $\Gamma=(1-\beta^{2})^{-0.5}$, and viewing
angle $\vartheta_{view}$ has three important consequences coming directly from
special relativity:

  * •

plasma moving in the jet will be seen in the observer’s frame as showing a
projected superluminal speed $\beta_{project}>1$. It can be shown that
$\beta_{project}=(\beta\sin\vartheta_{view})/(1-\beta\cos\vartheta_{view})$,
which can be as high as $\beta\Gamma$. Indeed, the detection of superluminal
motion in blazar jets is one of the key observations supporting the AGN
unified model (Lister et al., 2013).

  * •

the emission from the jet is boosted in the direction of movement, which
translates into a flux density in the observer’s frame111Here and in the
following, primed quantities refer to the jet’s frame, while non-primed
quantities to the observer’s frame. which scales as
$F_{\nu}=\delta^{3}F^{\prime}_{\nu^{\prime}}$, where the quantity
$\delta~{}=~{}[\Gamma(1-~{}\beta\cos\vartheta_{view})]^{-1}$ is defined as the
Doppler factor. Typical values for the Doppler factor in blazars are of the
order of ten, which means that the flux density in the observer’s frame is
typically a factor of a thousand higher than in the jet’s frame. Similarly, it
means that if we have a radio-galaxy and a blazar at the same redshift, and
with identical flux densities in the jet’s frame, the blazar will be a
thousand times brighter than the radio-galaxy as observed from Earth.

  * •

the effect of time-compression translates into a flux variability in the
observer’s frame much faster than the one in the jet’s frame: the variability
timescale, conventionally defined as the flux-doubling timescale $\tau_{var}$,
is equal to $(1+z)\ \tau_{var}^{\prime}/\delta$. For a typical blazar Doppler
factor of ten, a flare with a timescale in the jet’s frame of one hour, is
compressed in the observer’s frame into a minutes-long flare.

The spectral energy distribution (SED) 222We define SED the energy density
representation $\nu F_{\nu}$, where $F_{\nu}$ is the flux density in units of
erg cm-2 s-1 Hz-1. It is easy to show that while $F_{\nu}$ is different from
$F_{\lambda}$ (the flux density per unit wavelength) and $F_{E}$ (the flux
density per unit energy), the SED is constant and $\nu F_{\nu}=\lambda
F_{\lambda}=EF_{E}$. of blazars is characterized by a non-thermal continuum
extending from radio to very-high-energy $\gamma$-rays, with only subdominant
thermal emission in optical/UV associated with the accretion disk (the big-
blue-bump) and emission lines from the BLR (only seen in FSRQs, and not in BL
Lacs). This non-thermal continuum is comprised of two distinct radiative
components: the first one peaks in infrared-to-X-rays, the second one peaks in
the $\gamma$-ray band, ranging from MeV to TeV. A remarkable discovery by
$\gamma$-ray telescopes is that in FSRQs the majority of the energy is emitted
in the $\gamma$-ray band (Hartman et al., 1992), and in BL Lacs the energy
emitted in the $\gamma$-ray band is as much as the one emitted in the radio-
to-X-ray band (Punch et al., 1992). The SED peak frequencies are not fixed,
not within the blazar population, nor in a single blazar. FSRQs show typically
a first SED peak in the infrared band, and a second SED peak in the MeV band;
BL Lacs, on the other hand, show a huge variety of peak frequencies, from
infrared to X-rays, and from MeV to TeV. The peak frequency is used to further
classify BL Lacs (Padovani and Giommi, 1995; Abdo et al., 2010) into low-
frequency-peaked BL Lacs (LBLs, with $\nu_{peak}<10^{14}$ Hz, in the infrared
band, similarly to FSRQs) and high-frequency-peaked BL Lacs (HBLs, with
$\nu_{peak}>10^{15}$ Hz, in the UV/X-ray band). The transition among the two
sub-classes is smooth, and there exist intermediate-frequency-peaked BL Lacs
(IBLs, with $10^{14}$ Hz $<\nu_{peak}<10^{15}$ Hz, in optical, (Laurent-
Muehleisen et al., 1999)). It is not clear how the peak frequency distribution
ends: observations by Cherenkov telescopes in the TeV band have identified a
population of extremely high-frequency-peaked blazars (EHBLs) with
$\nu_{peak}>10^{18}$ Hz, in the hard-X-rays (Costamante et al., 2001; Bonnoli
et al., 2015), and, correspondingly, a high-energy SED peak extending into the
TeV band (Aharonian et al., 2007). More recent studies seem to indicate
however that the EHBL population is not homogeneous, with all various
combinations of soft/hard X-rays and soft/hard TeV spectra possible (Foffano
et al., 2019; Biteau et al., 2020; Costamante, 2020). It seems that the peak
frequency is anti-correlated with the blazar luminosity: the brightest blazars
are the ones with the lowest $\nu_{peak}$, while the faintest ones are the
ones with the highest $\nu_{peak}$. This anti-correlation defines the blazar
sequence (Fossati et al., 1998; Meyer et al., 2011; Ghisellini et al., 2017).
The existence of this sequence is still an open question in the blazar
research field, and several outliers to the sequence have been identified
(Padovani et al., 2003; Nieppola et al., 2008; Padovani et al., 2012; Giommi
et al., 2012). If the blazar sequence represents a true physical property of
blazars, it will help understanding what physical processes shape the blazar
SEDs (Böttcher and Dermer, 2002; Ghisellini and Tavecchio, 2008; Finke, 2013).

The framework described up to here is pretty solid, and has been tested
against a variety of observations at all wavelengths. Nonetheless, there is
still many questions which remain open. Only limiting the discussion to radio-
loud AGNs, it is still not completely understood (i)  how relativistic jets
are formed and powered, and why only some AGNs can produce them; (ii) why
there are two intrinsically different jet types as seen in FRI/FRII and BL
Lacs/FSRQs; (iii) how particles are accelerated in jets; (iv) how photons are
produced in jets. The last item is the subject of this review. I will
summarize the state-of-the-art answer to the following question: how do
supermassive-black-hole jets shine? For a review of recent progress on
observational results see (Blandford et al., 2019; Hovatta and Lindfors,
2020). For a broader discussion of relativistic jets at all scales (AGNs and
micro-quasars) see (Romero et al., 2017).

The blazar non-thermal continuum emission implies the existence of a non-
thermal population of particles that can radiate photons over the whole
electromagnetic spectrum. Leptonic radiative processes are the ones associated
with electrons/positrons in the jet; hadronic radiative processes are the ones
associated with protons and nuclei. In the following I will discuss how
photons (and as it will be shown later, neutrinos) can be emitted from jets,
assuming that there exists an efficient particle accelerator in the jet. For
particle acceleration in AGN jets, see (Kirk et al., 2000; Sironi et al.,
2015)

##  2 Leptonic Radiative Processes

###  2.1 Electron Synchrotron Emission

Synchrotron radiation is produced by charged particles moving in a magnetic
field. For sake of simplicity, let’s start with the assumption that the
emitting region is a sphere of radius $R$ in the jet, moving with bulk Lorentz
factor $\Gamma$, and filled with a tangled, homogeneous magnetic field
$B^{\prime}$. If a particle with electric charge $e$ and Lorentz factor
$\gamma^{\prime}=E^{\prime}/mc^{2}$ is put in the region, it radiates a
synchrotron power (in units of erg s-1 Hz-1) per frequency $\nu^{\prime}$
equal to 333For further details on this and other radiative mechanisms see
(Blumenthal and Gould, 1970; Rybicki and Lightman, 1979; Longair, 1994; Dermer
and Menon, 2009).

$P^{\prime}_{syn,\nu^{\prime}}(\nu^{\prime},\gamma^{\prime},\varphi^{\prime})=\frac{\sqrt{3}e^{3}B^{\prime}\sin{\varphi^{\prime}}}{mc^{2}}\frac{\nu^{\prime}}{\nu^{\prime}_{c}}\int^{\infty}_{\frac{\nu^{\prime}}{\nu^{\prime}_{c}}}{dx\
K_{5/3}(x)}$ (1)

where $\varphi^{\prime}$ is the pitch angle between the particle momentum and
the magnetic field lines, the frequency $\nu^{\prime}_{c}$ is defined as
$\nu^{\prime}_{c}(\gamma^{\prime},\varphi^{\prime})=\frac{3eB}{4\pi
mc^{2}}\gamma^{\prime 2}\sin{\varphi^{\prime}}$, and $K_{5/3}(x)$ is the
modified Bessell function of the second kind of order $5/3$. The sign of the
charge does not matter: the synchrotron emission from electrons and positrons
is not distinguishable, and in the following I will simply talk about
electrons.

Let’s now assume that in the emitting region there is not one single electron,
but an isotropic population of electrons with a power-law distribution of
Lorentz factors $\gamma^{\prime}_{e}$, defined between $\gamma^{\prime}_{min}$
and $\gamma^{\prime}_{Max}$:

$N^{\prime}_{e}(\gamma^{\prime}_{e})=N^{\prime}_{0}\gamma_{e}^{\prime-n}$ (2)

where $N^{\prime}_{0}$ is the normalization of the electron distribution in
units of cm-3, and $n$ the electron spectral index. The collective synchrotron
emissivity (in units of erg s-1 cm-3 Hz-1 sr-1) from this population of
electrons is then computed by integrating the synchrotron power over
$\gamma^{\prime}_{e}$ and $\varphi^{\prime}$ as

$j^{\prime}_{syn,\nu^{\prime}}(\nu^{\prime})=\frac{1}{8\pi}\int_{\gamma^{\prime}_{min}}^{\gamma^{\prime}_{Max}}{d\gamma^{\prime}_{e}\
N^{\prime}_{e}(\gamma^{\prime}_{e})}\ \int_{0}^{\pi}{d\varphi^{\prime}\
\sin{\varphi^{\prime}}\
P^{\prime}_{syn,\nu^{\prime}}(\nu^{\prime},\gamma^{\prime},\varphi^{\prime})}$
(3)

Given the emissivity, the $\nu F_{\nu}$ flux in the observer’s frame (in units
of erg cm-2 s-1) is computed as

$\nu
F_{\nu}(\nu)=\frac{4\pi}{3}\frac{R^{3}}{d_{L}^{2}}\delta^{4}\nu^{\prime}j^{\prime}_{syn,\nu^{\prime}}(\nu^{\prime})$
(4)

where $d_{L}$ is the luminosity-distance of the source. It can be shown by
performing analytically the integrals that if the electron distribution is a
power-law distribution with index $n$ as in Equation (2), then the
differential photon flux $dN/dE$ (which is a quantity commonly used in high-
energy astrophysics, equal to $F_{E}(E)/E$) is also a power-law distribution
with index $\Gamma_{dN/dE}=(n+1)/2$, the flux density $F_{\nu}$ has index
$\alpha=(n-1)/2$, and the energy flux $\nu F_{\nu}(\nu)$ has index
$p=(n-3)/2$. Given that $\nu=\delta\nu^{\prime}/(1+z)$, and that
$j^{\prime}_{\nu}\propto\nu^{\prime-\alpha}$, the energy flux is proportional
to $\delta^{3+\alpha}$.

At low energies, an important physical process is synchrotron self-absorption:
the photons that the electrons radiate can be re-absorbed by the electrons
themselves. Following Einstein’s approach, it is possible to calculate the
self-absorption coefficient (in units of cm-1) as

$\mu^{\prime}_{syn}(\nu^{\prime})=-\frac{1}{16\pi
m\nu^{2}}\int_{\gamma^{\prime}_{min}}^{\gamma^{\prime}_{Max}}{d\gamma^{\prime}\
\gamma^{\prime
2}\frac{d}{d\gamma}\left(\frac{N^{\prime}_{e}(\gamma^{\prime}_{e})}{\gamma^{\prime
2}}\right)}\ \int_{0}^{\pi}{d\varphi^{\prime}\ \sin{\varphi^{\prime}}\
P^{\prime}_{syn,\nu^{\prime}}(\nu^{\prime},\gamma^{\prime},\varphi^{\prime})}$
(5)

In the presence of this absorption term, the intensity of the synchrotron
radiation (in units of erg s-1 cm-2 Hz-1 sr-1) from this spherical emitting
region in the jet can be computed as (Bloom and Marscher, 1996)

$I^{\prime}_{syn,\nu^{\prime}}(\nu^{\prime})=\frac{j^{\prime}_{syn,\nu^{\prime}}(\nu^{\prime})}{\mu^{\prime}_{syn}(\nu^{\prime})}\left\\{1-\frac{2}{\tau(\nu^{\prime})^{2}}\left[1-\mathrm{e}^{-\tau(\nu^{\prime})}(\tau(\nu^{\prime})+1)\right]\right\\}$
(6)

where $\tau(\nu^{\prime})=2R\mu^{\prime}_{syn}(\nu^{\prime})$. The effect of
self-absorption on the radiated synchrotron spectrum is a change in the power-
law index: at frequencies below $\nu^{\prime}_{self-abs}$ the radiated flux
density is proportional to $\nu^{5/2}$. Said otherwise, at frequencies below
$\nu^{\prime}_{self-abs}$ the emitting region is optically thick to
synchrotron radiation. Self-absorption in the emitting region has one
important consequence: a single emitting region in the relativistic jet cannot
reproduce blazar SEDs in the radio band. The emission in the radio band must
come from other regions which are necessarily optically thin to synchrotron
radiation, as will be discussed in Section 2.1.2.

####  2.1.1 Self-Consistent Electron Distribution

As electrons radiate synchrotron photons, they lose energy, or, using the
jargon of astrophysicists, cool. When calculating the synchrotron spectrum
from a stationary electron population in the jet, it is thus important to
first compute the self-consistent electron energy distribution at equilibrium.
The differential equation that describes this process is

$\frac{\partial}{\partial
t}N^{\prime}_{e}(\gamma^{\prime}_{e},t)=\frac{\partial}{\partial\gamma^{\prime}_{e}}\left[\gamma^{\prime}_{e}\frac{N^{\prime}_{e}(\gamma^{\prime}_{e},t)}{\tau_{syn}(\gamma^{\prime}_{e})}\right]-\frac{N^{\prime}_{e}(\gamma^{\prime}_{e},t)}{\tau_{ad}}+Q^{\prime}_{e}(\gamma^{\prime}_{e})$
(7)

where $Q^{\prime}_{e}(\gamma^{\prime}_{e})$ is the function describing
continuous (time-independent) injection of electrons,
$\tau_{syn}(\gamma^{\prime})$ is the synchrotron cooling timescale, and
$\tau_{ad}$ is the adiabatic timescale. The latter is the timescale associated
with the expansion of the emitting region as it travels along the jet, and it
is energy-independent. It is usually assumed to be $\tau_{ad}\simeq R/c$. From
a mathematical point of view, the term in $\tau_{ad}$ can also be interpreted
as an energy-independent particle escape from the emitting region. The
synchrotron cooling timescale is defined as
$\tau_{syn}(\gamma_{e}^{\prime})=\left(\frac{1}{\gamma_{e}^{\prime}}\frac{\partial\gamma_{e}^{\prime}}{\partial
t}\right)^{-1}$ and it is equal to

$\tau_{syn}(\gamma^{\prime}_{e})=\frac{3m_{e}c}{4\sigma_{T}u^{\prime}_{B}}\frac{1}{\gamma^{\prime}_{e}}$
(8)

where $u^{\prime}_{B}=B^{\prime 2}/8\pi$ is the magnetic energy density. The
particle distribution at equilibrium can be computed putting Equation (7)
equal to zero. The integral solution is (Inoue and Takahara, 1996)

$N^{\prime}_{e}(\gamma^{\prime}_{e})=e^{-\gamma^{\prime}_{e,break}/\gamma^{\prime}_{e}}\frac{\gamma^{\prime}_{e,break}\
\tau_{ad}}{\gamma^{\prime 2}_{e}}\int_{\gamma^{\prime}_{e}}^{\infty}d\zeta
Q^{\prime}_{e}(\zeta)e^{\gamma^{\prime}_{e,break}/\zeta}$ (9)

where $\gamma^{\prime}_{e,break}=\frac{3mc^{2}}{4\sigma_{T}u^{\prime}_{B}R}$
is the value of $\gamma^{\prime}_{e}$ that satisfies
$\tau_{syn}(\gamma^{\prime}_{e})=\tau_{ad}$. It can be shown that if the
injection term is
$Q^{\prime}_{e}(\gamma^{\prime}_{e})=Q^{\prime}_{0}\gamma^{\prime-q}_{e}$, the
effect of synchrotron cooling on the particle distribution is a change in the
spectral index, which softens and becomes equal to $q+1$ for particle Lorentz
factors above $\gamma^{\prime}_{e,break}$. From the point of view of the
energy flux distribution, it also softens above $\nu^{\prime}_{e,break}$, and
becomes equal to $p+1/2$. This is true only if $q\geq 2$. If $q<2$, the
particle distribution shows a pile-up at $\gamma^{\prime}_{e,break}$
(Kardashev, 1962; Saugé and Henri, 2004). Another specific case happens when
the injection function is also characterized by a minimum Lorentz factor
$\gamma^{\prime}_{min,inj}$. If
$\gamma^{\prime}_{e,break}>\gamma^{\prime}_{min,inj}$ (slow-cooling regime),
then the solution is the same as before. If however
$\gamma^{\prime}_{e,break}<\gamma^{\prime}_{min,inj}$ the steady-state
solution is completely cooled (fast-cooling regime), and characterized by an
index $q+1$ above $\gamma^{\prime}_{min,inj}$, and an index of $2$ between
$\gamma^{\prime}_{e,break}$ and $\gamma^{\prime}_{min,inj}$ (Sari and Esin,
2001).

We are now in the position to quantitatively compute a self-consistent
synchrotron spectrum from an electron distribution in a spherical emitting
region in the relativistic jet. The SED can be parametrized by a broken-power-
law distribution, proportional to $\nu^{7/2}$ below the self-absorption
frequency $\nu_{self-abs}$, to $\nu^{-p}$ between $\nu_{self-abs}$ and
$\nu_{break}$, and to $\nu^{-(p+1/2)}$ above $\nu_{break}$.

It would be surprising if such a simple scenario would be able to reproduce
blazar SEDs, and indeed, it does not. Electron synchrotron radiation from a
single emitting region in the jet can reproduce the first blazar SED
component, at one price: while the emission is consistent with a broken-power-
law distribution, the electron indexes $n_{1,2}$ below and above
$\gamma^{\prime}_{e,break}$ are not consistent with a simple break of 1 as
calculated above. What we usually observe are much stronger breaks, with
$n_{2}-n_{1}>1$, with $n_{1}=1.5-2.5$ and $n_{2}=3.5-4.5$. For this reason,
the most common approach is to consider $n_{1,2}$ and
$\gamma_{e,break}^{\prime}$ as free independent model parameters (Tavecchio et
al., 2010; Zhang et al., 2012), and discuss a posteriori their consistency
with a synchrotron cooling break. In a minority of cases, observations are
consistent with a broken power-law with a spectral break of one (Aharonian et
al., 2008; Aleksić et al., 2015; Archambault et al., 2015; H. E. S. S.
Collaboration et al., 2011). Alternatively, it is possible to assume that the
electron distribution at equilibrium is in the fast-cooling regime. In this
case, if $n_{1}$ is consistent with $2.0$, it is possible to fit the SED
assuming an injection index $q\simeq 2.5-3.5$ (Aliu et al., 2013; Baloković et
al., 2016). This approach has however two drawbacks: the first one is that it
needs an extra free parameter $\eta$ that is the scaling factor of the
adiabatic timescale that has to be increased to be in the fast-cooling regime
$\tau_{ad}=\eta R/c$; the second one is that an injection index $q\simeq 3.0$
is significantly softer than standard acceleration mechanisms. The reason we
observe stronger breaks may be that one of the assumptions we made are not
correct: there could be multiple acceleration/emitting regions (Marscher,
2014), non-linear cooling terms (relevant when adding the inverse-Compton
process, see Section 2.2), a non-homogeneous magnetic field, incomplete
cooling (i.e., the SED is from an electron population that is not at
equilibrium) (Chiang and Böttcher, 2002; Graff et al., 2008). Alternatively,
the injected particle distribution may not follow a simple power-law function,
and could show intrinsic curvature in the form of a log-parabolic function
(Massaro et al., 2004; Tramacere et al., 2009), or a power-law one with
exponential cut-off at $\gamma^{\prime}_{Max}$. The observational measurements
of $n_{1}$ and $n_{2}$ may thus be the asymptotic power-law fits over a narrow
energy band.

####  2.1.2 Low-Energy SED Modeling

As discussed above, a single emitting region in the jet radiating synchrotron
photons cannot explain the observed spectrum from radio up to optical-X-rays,
due to self-absorption in the radio band. Blazar SED modeling is thus,
inevitably, at least two-zones: a small, compact emitting region is
responsible for the high-energy part of the first SED component (optical
frequencies and above), the one that shows rapid variability; a larger, less
dense emitting region is where the radio photons are produced. It is natural
to consider that this radio emission comes from the extended radio-jet. The
exact geometry is not clear: the radio emission can be produced downstream
with respect to the optical-X-ray one, or be associated with a second flow
that encompasses the inner jet (Sol et al., 1989), or the whole emission could
be produced not in a small overdensity in the jet, but in a conical jet, with
varying magnetic field (Ghisellini et al., 1985; Potter and Cotter, 2012). For
most blazars, a two-zone model is enough to fit the first SED component
(Archer et al., 2018), although the transition between the two components
remains unclear, with potential contribution of the extended region up to the
optical/UV band (H. E. S. S. Collaboration et al., 2012; Ahnen et al., 2016).

A fit of the model to the observations can be used to put constraints on the
physics of the emitting region, or of the acceleration mechanisms. The number
of free parameters of the synchrotron model is nine: three for the emitting
region ($\delta$, $B$, and $R$), and six for the particle distribution
($N_{0}$, $\gamma^{\prime}_{min,break,Max}$, and $n_{1,2}$). If we force a
self-consistent electron distribution, $\gamma^{\prime}_{break}$, and $n_{2}$
can be expressed as a function of other parameters, and the number of free
parameters becomes seven; if we further assume that
$\gamma^{\prime}_{min,Max}$ are respectively low-enough ($\simeq$1) and high
enough ($\geq$106) that they cannot be constrained by observations, we are
left with five free parameters. Still, the true number of observables is much
lower: the synchrotron spectrum is fully characterized by its peak frequency
and luminosity, and by its index, i.e., we have only three observables. This
is true both for the small plasmoid in the jet, and for the large component
emitting in the radio band. We could constrain the model further if we could
see the self-absorption frequency, but the one for the plasmoid is hidden
below the extended-jet emission. Synchrotron modeling of the low-energy SED of
blazars is, by itself, degenerate , which means that although it is possible
to find a set of model parameters that can reproduce the SED, we cannot
constrain the parameter values.

####  2.1.3 Time-Dependent Modeling

Until now we focused on the synchrotron emission by a stationary distribution
of electrons in the jet, comparing theoretical expectations with blazar SEDs,
that can be seen as snapshots of the blazar emission at a given instant. One
of the key observing properties of blazars is however variability, and key
pieces of physical information can be extracted by studying blazar flares. It
is possible to simulate blazar light-curves by solving Equation (7) as a
function of time, following the evolution of the electron distribution with
time, and the associated synchrotron radiation.

As a first approximation, a synchrotron blazar flare can be modeled with an
instantaneous injection of particles $Q^{\prime}_{e}(t^{\prime})=Q_{0}\
\delta(t^{\prime}-t^{\prime}_{0})$. In this case the time-evolution is fully
determined by the cooling timescale $\tau^{\prime}_{syn}$ for Lorentz factors
higher than $\gamma^{\prime}_{break}$. The measure of $\tau_{var}$ can thus
provide a direct measurement of
$\tau^{\prime}_{syn}=\tau^{\prime}_{syn}(B^{\prime},\gamma^{\prime}_{e})$ and
thus help constrain the parameter space.

Increasing complexity, all plausible acceleration processes are energy-
dependent. In a more realistic scenario we can thus define an acceleration
timescale $\tau^{\prime}_{acc}=\tau^{\prime}_{acc}(\gamma^{\prime}_{e})$. The
flare is thus expected to be characterized by an asymmetric profile, due to
the different values of $\tau^{\prime}_{acc}$ and $\tau^{\prime}_{syn}$. The
energy dependency of both acceleration and cooling implies that the spectrum
changes during the flare. In a plot $\Gamma_{dN/dE}$ vs $F_{\nu}$ this
translates into a hysteresis cycle. Equivalently, it will appear as a time-lag
among light-curves computed at different energies (Kardashev, 1962;
Mastichiadis and Kirk, 1997; Kirk et al., 1998; Böttcher and Chiang, 2002).
Both hysteresis and time-lags have been observed in X-ray flare of blazars
(Takahashi et al., 1996; Kataoka et al., 2000; Zhang et al., 2002; Sokolov et
al., 2004), and can be used to simultaneously constrain the magnetic field
strength and the Doppler factor of the emitting region (with a degeneracy
between the two), resulting in $B\simeq 0.1$ for $\delta\simeq 10$. When
modeling blazar flares, there is another timescale that needs to be taken into
account that is the light-crossing time $R/c$: the variability will thus be
diluted as observed from Earth if $R/c>\tau^{\prime}_{syn}$ (Chiaberge and
Ghisellini, 1999).

###  2.2 Inverse-Compton Emission

Compton scattering is the interaction between a charged particle and a photon.
Following Section 2.1 on synchrotron radiation, let us start again with an
electron. In the direct Compton scattering a photon scatters over an electron
at rest in an atom, and loses energy in the interaction. The photon wavelength
after scattering changes as

$\Delta\lambda=\frac{h}{m_{e}c}(1-\cos\zeta)$ (10)

with $\zeta$ the scattering angle, and with $\frac{h}{m_{e}c}=\lambda_{C}$ the
Compton wavelength.

In the inverse-Compton scattering a high-energy electron interacts with a low-
energy photon, and transfers part of its energy to it. The inverse-Compton
scattering is a very efficient process to produce X-ray and $\gamma$-ray
photons in astrophysical sources starting from low-energy photon fields. Let
us mark with a star the reference frame at rest with the electron, and let’s
mark with the $1,2$ subscripts the values before and after scattering,
respectively. The photon a-dimensional energies are
$\varepsilon^{\star}_{1,2}~{}=~{}h\nu^{\star}_{1,2}/m_{e}c^{2}$. The electron
Lorentz factors are $\gamma^{\star}_{e;1,2}=E^{\star}_{e;1,2}/m_{e}c^{2}$,
with $\gamma^{\star}_{e;1}=1$. It can be shown with simple kinematic and
energy/momentum conservation equations that in the reference frame at rest
with the electron

$\varepsilon^{\star}_{2}=\frac{\varepsilon^{\star}_{1}}{1+\varepsilon^{\star}_{1}(1-\cos\xi^{\star})}$
(11)

where $\xi^{\star}$ is the angle between the incoming and scattered directions
of the photon. When $\varepsilon^{\star}_{1}\ll 1$, Equation (11) simplifies
as $\varepsilon^{\star}_{2}\simeq\varepsilon^{\star}_{1}$, which is the photon
energy in the electron rest frame is unchanged. The condition
$\varepsilon^{\star}_{1}\ll 1$ represents the Thomson regime, where the photon
energy is much lower than the electron rest mass. The equations describing the
change of reference frame from and to the observer (unstarred) are:

$\varepsilon^{\star}_{1}=\varepsilon_{1}\gamma_{e}(1-\beta_{e}\cos\vartheta_{1})$
(12)
$\varepsilon_{2}=\varepsilon^{\star}_{2}\gamma_{e}(1+\beta_{e}\cos\vartheta^{\star}_{2})$
(13)

where $\vartheta_{1}$ is defined as the angle between the electron and the
incoming photon, in the observer’s frame, and $\vartheta^{\star}_{2}$ the one
of the scattered photon, in the electron’s rest frame. In the Thomson regime,
we thus have

$\varepsilon_{2}\simeq\varepsilon^{\star}_{1}\gamma_{e}(1+\beta_{e}\cos\vartheta^{\star}_{2})\simeq\varepsilon_{1}\gamma^{2}_{e}(1-\beta_{e}\cos\vartheta_{1})(1+\beta_{e}\cos\vartheta^{\star}_{2})$
(14)

In the case of a relativistic electron with $\gamma_{e}\gg 1$ and
$\beta_{e}\simeq 1$, the photon energy in the observer’s frame after
scattering can be as high as $\varepsilon_{2}\simeq
4\varepsilon_{1}\gamma^{2}_{e}$. This occurs in the so-called head-on
approximation, which is for $\cos\vartheta_{1}=-1$ and
$\cos\vartheta^{\star}_{2}=1$.

The expression of the inverse-Compton cross-section is equal to $\sigma_{T}$
in the Thomson regime, while for the general case it has been calculated by
Klein and Nishina (1929), and is known as the Klein–Nishina formula. The
differential cross-section, in the reference frame of the electron, is

$\frac{d\sigma_{C}}{d\varepsilon_{2}^{\star}d\Omega^{\star}}=\frac{3\sigma_{T}}{16\pi}\left(\frac{\varepsilon_{2}^{\star}}{\varepsilon_{1}^{\star}}\right)^{2}\left(\frac{\varepsilon_{2}^{\star}}{\varepsilon_{1}^{\star}}+\frac{\varepsilon_{1}^{\star}}{\varepsilon_{2}^{\star}}-\sin^{2}\xi^{\star}\right)$
(15)

where $d\Omega^{\star}=d\varphi^{\star}d\cos\xi^{\star}$. The integral cross-
section is

$\sigma_{C}(\varepsilon_{1}^{\star})=\frac{3\sigma_{T}}{8\varepsilon_{1}^{\star
2}}\left(4+\frac{2\varepsilon_{1}^{\star
2}(1+\varepsilon_{1}^{\star})}{(1+2\varepsilon_{1}^{\star})^{2}}+\frac{\varepsilon_{1}^{\star
2}-2\varepsilon_{1}^{\star}-2}{\varepsilon_{1}^{\star}}\ln(1+2\varepsilon_{1}^{\star})\right)$
(16)

The cross-section drops steadily above $\varepsilon^{\star}_{1}\gg 1$, and
asymptotically behaves as
$\sigma_{C}(\varepsilon_{1}^{\star})~{}\simeq~{}\frac{3}{8}\frac{\sigma_{T}}{\varepsilon^{\star}_{1}}\left(\log(2\varepsilon^{\star}_{1})+1/2\right)$.
The regime $\varepsilon^{\star}_{1}\gg 1$ is known as the Klein–Nishina
regime. Solving Equations (11) and (13) for the Klein–Nishina case, it is easy
to show that in the best case $\varepsilon_{2}\simeq\gamma_{e}$.

Let’s discuss now the inverse-Compton emission from a distribution of
electrons $N^{\prime}_{e}(\gamma^{\prime}_{e})$ scattering a distribution of
photons $n^{\prime}_{ph}(\varepsilon^{\prime})$. The computation of the
inverse-Compton emissivity is analytically complex, requiring transformation
of the Klein–Nishina cross-section from the electron rest frame to the jet
frame. It has been investigated under several approximations (e.g., Thomson
and Klein–Nishina approximation, head-on collisions, monochromatic photons).
It is best solved numerically performing integrals over the electron and
photon distributions times the cross-section (Equation (15)). For an isotropic
distribution of electrons and photons, the inverse-Compton emissivity can be
expressed as:

$j^{\prime}_{C}(\nu^{\prime})=\frac{h^{2}\nu^{\prime}}{4\pi m_{e}c^{2}}\int
d\varepsilon^{\prime}n^{\prime}_{ph}(\varepsilon^{\prime})\int
d\gamma^{\prime}N^{\prime}_{e}(\gamma^{\prime}_{e})\
C(\varepsilon^{\prime},\gamma^{\prime}_{e},\nu^{\prime})$ (17)

where $C(\varepsilon^{\prime},\gamma^{\prime}_{e},\nu^{\prime})$ is the Jones’
Compton kernel (Jones, 1968).

As for the synchrotron radiation, the fact that part of the electron energy is
transferred to the photons imply that electrons are cooling in the process. It
is thus important to include the inverse-Compton cooling losses when computing
the self-consistent electron distribution at equilibrium (see Section 2.1.1).
Let’s start first with the Thomson regime: in this case the electron losses
due to inverse-Compton scattering (for an isotropic photon field) take a form
which is very similar to the synchrotron one, and the associated timescale is

$\tau_{C,Th}(\gamma^{\prime}_{e})=\frac{3mc}{4\sigma_{T}u^{\prime}_{ph}}\frac{1}{\gamma^{\prime}_{e}}$
(18)

where $u^{\prime}_{ph}$ is the energy density of the scattered photon field.
When computing the electron distribution at equilibrium, it is just possible
to simply use a single cooling timescale
$\tau_{cooling}(\gamma^{\prime}_{e})=\frac{3mc}{4\sigma_{T}(u^{\prime}_{B}+u^{\prime}_{ph})}\frac{1}{\gamma^{\prime}_{e}}$,
and the effect on the electron distribution is again a broken-power-law with
indexes $n$ and $n+1$ below and above $\gamma^{\prime}_{break}$ respectively
(in the slow cooling regime).

In the Klein–Nishina regime, the inverse-Compton losses are also suppressed:
if the inverse-Compton losses dominate over the synchrotron ones for high
Lorentz factors (as could occur in the presence of bright external fields),
the result on the steady-state electron distribution is a hardening at the
highest energies (Moderski et al., 2005).

####  2.2.1 Synchrotron-Self-Compton

Going back to the spheroidal plasmoid travelling with Lorentz factor $\Gamma$
in the SMBH relativistic jet, filled with both an homogeneous magnetic field
$B$ and an electron population $N^{\prime}_{e}(\gamma_{e}^{\prime})$, it is
natural that the synchrotron photons will suffer inverse-Compton scattering
over the same electrons population that produced them. This kind of emission
is called synchrotron-self-Compton (SSC), and it is a very efficient process
to transfer energy from particles to photons (Gould, 1979).

In the Thomson regime, it can be shown analytically that if the electron
distribution $N^{\prime}_{e}(\gamma_{e}^{\prime})$ is a power-law with index
$n$, not only the synchrotron energy flux has index $p=(n-3)/2$, but the SSC
energy flux has also index $p=(n-3)/2$. The SSC emission is thus mirroring the
synchrotron radiation at higher energies. It is thus natural to investigate
SSC emission as the radiative process responsible for the high-energy SED
component in blazars. Now let us investigate two key aspects: how high are
these SSC energies, and what is the relative luminosity $L_{SSC}/L_{syn}$
(Tavecchio et al., 1998):

  * •

Let’s assume that the synchrotron peak frequency $\nu_{syn,peak}$ is directly
related to $\gamma^{\prime}_{break}$, i.e., let’s assume that $n_{1}<3$ and
$n_{2}>3$. In this case

$\nu_{syn,peak}=\frac{e}{2\pi m_{e}c}B^{\prime}\gamma^{\prime
2}_{break}\frac{\delta}{1+z}\ \textrm{Hz}$ (19)

The SSC peak frequency in the Thomson regime can be easily computed assuming
that the electrons at $\gamma^{\prime}_{break}$ are the ones primarily
scattering the photons at $\nu_{syn,peak}$, and thus
$\nu_{SSC,peak}~{}=~{}4/3\ \nu_{syn,peak}\gamma^{\prime 2}_{break}$. Assuming
typical values for HBLs, i.e., a $\nu_{syn,peak}\simeq 0.1$ keV, in the soft
X-rays, and $\gamma^{\prime}_{break}\simeq 10^{4}$ (see next Section),
$\nu_{SSC,peak}$ can thus reach 10 GeV. The ratio
$\nu_{SSC,peak}/\nu_{syn,peak}$ can be used to put a strong constraint on the
model parameters, given that

$\gamma^{\prime}_{break}=\left(\frac{3}{4}\frac{\nu_{SSC,peak}}{\nu_{syn,peak}}\right)^{1/2}$
(20)

  * •

the total luminosities of the synchrotron and SSC components $L_{syn,SSC}$ can
be expressed, in a first approximation, as
$L_{syn,SSC}\simeq\nu^{\prime}_{syn,SSC;peak}\
L_{syn,SSC}(\nu^{\prime}_{syn,SSC;peak})$, which is their peak SED
luminosities. The ratio of the two SED peak fluxes depends only on the
magnetic energy density and the synchrotron photons energy density
$u^{\prime}_{ph,syn}$:

$\frac{\nu_{SSC,peak}\ F_{SSC}(\nu_{SSC,peak})}{\nu_{syn,peak}\
F_{syn}(\nu_{syn,peak})}=\frac{\nu^{\prime}_{SSC,peak}\
L_{SSC}(\nu^{\prime}_{SSC,peak})}{\nu^{\prime}_{syn,peak}\
L_{syn}(\nu^{\prime}_{syn,peak})}=\frac{u^{\prime}_{ph,syn}}{u^{\prime}_{B}}$
(21)

In the Klein–Nishina regime on the other hand, the SSC emission is suppressed,
producing a major correction to the shape of the SSC component at higher
energies. As we increase the peak of the synchrotron emission, the SSC peak
frequency increases quadratically until the scattering enters into the
Klein–Nishina regime, and higher-energy photons do not contribute any more to
the emission. Even for the Klein–Nishina regime, simple relations between the
frequencies and luminosities of the synchrotron and SSC components can be
computed (Tavecchio et al., 1998). In is important to underline that the
computation of best-fit solutions for the SSC model relies on a good estimate
of the peak frequencies and fluxes. Given the wide variety of peak frequencies
in the blazar population, these may fall in relatively unexplored energy
bands.

When discussing the modeling of the low-energy SED component of blazars as
synchrotron radiation (see Section 2.1.2), we showed that the modeling of the
synchrotron component alone is degenerate. This is not the case when we add
the information from the SSC radiative component. We recall that the number of
free parameters for a spherical plasmoid in the jet, filled with an electron
population parametrized as broken-power-law function, is nine. Assuming that
$\gamma^{\prime}_{min,Max}$ are respectively low and high enough that they
cannot be constrained by the data, we are left with seven model parameters.
With the information from both SED components, the number of observables is
six: the frequencies and luminosities of the synchrotron and SSC peaks, and
the two indexes of the synchrotron distribution before and after the peak. We
are left with just one missing constraint, which can be the self-consistency
of the electron distribution, or can be the observed fastest variability
timescale $\tau_{var}$. It can be related to the size of the emitting region
$R$ and the Doppler factor $\delta$ via a causality argument

$R\leq\frac{\tau_{var}}{c}\frac{\delta}{(1+z)}$ (22)

A one-zone SSC model can describe well the SEDs of HBLs (assuming that the
radio emission comes from a larger, less dense region in the jet). An example
of an SSC modeling of the HBL Markarian 421 is shown in the top-left panel of
Figure 2. In the recent years, a number of algorithms have been published to
fully constrain the SSC model parameter space, and have been applied
successfully to HBLs (Finke et al., 2008; Mankuzhiyil et al., 2011; Zhang et
al., 2012; Cerruti et al., 2013; Ahnen et al., 2017). Typical best-fit
solutions are:

$\begin{array}[]{l}\delta\simeq 10\textrm{--}50\\\ R\simeq
10^{16\textrm{--}17}~{}\textrm{cm}\\\ B^{\prime}\simeq
10\textrm{--}100~{}\textrm{mG}\\\ \gamma^{\prime}_{min}\leq 10^{2-3}\\\
\gamma^{\prime}_{break}\simeq 10^{3\textrm{--}4}\\\ \gamma^{\prime}_{Max}\geq
10^{5\textrm{--}6}\\\ n_{1}\simeq 1.5\textrm{--}2.5\\\ n_{2}\simeq
3.5\textrm{--}4.5\\\ L_{e}=10^{43-44}\ \textrm{erg s}^{-1}\end{array}$ (23)

where $L_{e}=2\pi R^{2}c\Gamma^{2}u^{\prime}_{e}$ is the electron luminosity
in the observer’s frame, and $u^{\prime}_{e}~{}=~{}m_{e}c^{2}\int
d\gamma_{e}^{\prime}\gamma_{e}^{\prime}N(\gamma_{e}^{\prime})$ is the electron
energy density. The SSC scenario fails to reproduce the SED of other blazar
sub-classes: it cannot reproduce the SEDs of LBLs, nor FSRQs; it can reproduce
the SEDs of EHBLs at the expenses of a high Doppler factor ($\delta\geq 50$)
and of a high value of $\gamma^{\prime}_{min}\geq 10^{3}$ (Biteau et al.,
2020; Cerruti et al., 2015).

####  2.2.2 Synchrotron-Self-Compton: Time Signatures

Now that we have a model which is able to describe the blazar SEDs (at least
for HBLs), we can add the time information to further test it and constrain
the parameter values. As for the synchrotron case, we can compute the
evolution of the SED with time by solving Equation (7) (including inverse-
Compton losses) and calculating both synchrotron and SSC emission as a
function of time. One difficulty of SSC time-dependent modeling is that the
inverse-Compton cooling term depends on the density of the synchrotron
photons, which depends on the integral of the electron distribution itself as
a function of time (which is what we want to solve). In presence of
significant SSC losses the differential equation becomes thus non-linear and
the solution is significantly more complex (Zacharias and Schlickeiser, 2010;
Zacharias, 2014). Blazar light-curves in the $\gamma$-ray band have been
successfully fit with the SSC model (H. E. S. S. Collaboration et al., 2012;
Chen et al., 2011). Similarly to synchrotron radiation, time-lags and
hysteresis are also expected in SSC light-curves due to the energy dependency
of both the acceleration and cooling terms (Perennes et al., 2020). But
contrarily to what is observed in X-rays, there has not been any significant
detection of time-lags or of hysteresis in the $\gamma$-ray band (Albert et
al., 2007; Abeysekara et al., 2017). This negative result is likely due to the
intrinsic lower sensitivity of gamma-ray telescopes compared to X-ray ones,
and not to a failure of the model.

Figure 2: Spectral energy distribution of the high-frequency-peaked BL
Lacertae object Markarian 421, showing in gray archival data from the SSDC SED
builder (https://tools.ssdc.asi.it/SED/). Top left: one-zone SSC modeling of
the SED, with in red synchrotron emission by electrons and in green
synchrotron-self-Compton emission; model parameters are $\delta=30$,
$R=1.6\times 10^{16}\ \textrm{cm}$, $B^{\prime}=40\ \textrm{mG}$,
$\gamma^{\prime}_{e,min}=100$, $\gamma^{\prime}_{e,break}=6\times 10^{4}$,
$\gamma^{\prime}_{e,Max}=10^{6}$, $n_{e,1}=1.8$, $n_{e,2}=3.5$,
$L_{e}=5.1\times 10^{43}\ \textrm{erg s}^{-1}$. Top right: one-zone lepto-
hadronic modeling of the SED, with in blue synchrotron emission by protons, in
dotted-blue synchrotron emission by secondary pairs from proton–synchrotron,
in dotted-red synchrotron emission by the Bethe–Heitler cascade, in dotted
yellow synchrotron emission by the $\pi^{0}$ cascade, and in dotted pink
synchrotron emission by the $\pi^{\pm}$ cascade; neutrino emission is in
violet; model parameters are same as the first plot, with in addition a proton
distribution with $\gamma^{\prime}_{p,min}=1$,
$\gamma^{\prime}_{p,Max}=2\times 10^{8}$, $n_{p}=1.8$, $L_{p}=1.3\times
10^{47}\ \textrm{erg s}^{-1}$ Bottom: one-zone hadronic modeling of the SED;
model parameters are $\delta=30$, $R=2.5\times 10^{15}\ \textrm{cm}$,
$B^{\prime}=50\ \textrm{G}$, $\gamma^{\prime}_{e,min}=300$,
$\gamma^{\prime}_{e,Max}=3\times 10^{4}$, $n_{e}=2.9$, $L_{e}=5.4\times
10^{39}\ \textrm{erg s}^{-1}$, $\gamma^{\prime}_{p,min}=1$,
$\gamma^{\prime}_{p,Max}=3\times 10^{9}$, $n_{p}=1.9$, $L_{p}=4.7\times
10^{43}\ \textrm{erg s}^{-1}$.

Another important prediction of the SSC scenario is the correlation between
the variability of the two SED components, which are expected to vary in
concert. If the variability is simply due to an achromatic change in the
particle density in the emitting region, the SSC component is expected to vary
as the square of the synchrotron component (due to both the change in the
particle distribution, and in the synchrotron photons). A correlation plot
between X-rays and $\gamma$-rays is thus expected to show a simple quadratic
relation. If on the other hand the variability is energy-dependent, the
expected correlations are more complex, and depend on the energy bands of the
observations, the SED peak frequencies, the energy-dependency of the injection
(Katarzyński et al., 2005). Correlation studies among different wavebands are
indeed a powerful tool to test radiative models when the statistics in an
individual band does not allow light-curve fitting (Fossati et al., 2008;
Acciari et al., 2011; Aleksić et al., 2015; Furniss et al., 2015; Baloković et
al., 2016; González et al., 2019).

####  2.2.3 External Inverse-Compton

Synchrotron photons are not the only low-energy photons that can be up-
scattered to $\gamma$-rays. The SMBH environment is very bright, and there is
a huge variety of photon fields that can go through inverse-Compton scattering
with electrons in the jet. This radiative process is usually called External-
inverse-Compton (EIC) and depends strongly on the properties of the external
field and on the exact location of the emitting region $r$ (measured from the
SMBH). The computation of the EIC component is similar to the SSC one, but
requires as first step the transformation of the external photon field
distribution into the jet’s frame. For blazar modeling, the relevant photon
fields, starting from the AGN core outwards, are:

  * •

the SMBH accretion disk is a reservoir of thermal photons. EIC emission over
disk photons is however hindered due to the Doppler deboosting (Dermer and
Schlickeiser, 1993): the plasma in the jet is travelling away from it with
Lorentz factor $\Gamma$, and thus the disk photon field energy density, in the
reference frame of the plasmoid, is strongly suppressed by a factor
$\Gamma^{-3}$. EIC scattering over disk photons can thus dominate the overall
$\gamma$-ray emission only if the emitting region is located close to the disk
itself. In this case the plasma can see part of the disk at angles large
enough to reduce the de-boosting.

  * •

the BLR produces bright emission lines (Ly$\alpha$ being the dominant one)
that can serve as target photon field for EIC emission (Sikora et al., 1994).
In this case the emission is boosted or deboosted in the reference frame of
the emitting region as a function of $r$ and $r_{BLR}$. It is important to
underline that the BLR is structured, and that different emission lines are
produced at different distances $r_{BLR,i}$ from the SMBH. The size of the BLR
can be expressed as a function of its luminosity (or other proxies such as the
luminosity of the continuum at $5100$ Å, or the disk luminosity (Greene and
Ho, 2005; Kaspi et al., 2000)). The EIC emission is thus a superposition of
several components coming from the different lines (Finke, 2016), with their
relative strength depending on $r$. A second radiative component coming from
the BLR is the photon field from the accretion disk which can be Thomson-
scattered by electrons in the BLR and thus seen boosted in the reference frame
of the blob. This component depends also on the optical depth $\tau_{BLR}$ of
the BLR. A parameter which is not so well known is the aperture angle
$\alpha_{BLR}$ of the BLR: it can have a spherical geometry ($\alpha_{BLR}=0$)
or be flattened over the disk ($\alpha_{BLR}\rightarrow\pi/2$) (Tavecchio and
Ghisellini, 2012; Lei and Wang, 2014). The assumed BLR geometry has a direct
impact on the photon field seen in the blob reference frame, and in the limit
of a flat BLR the computation of the associated EIC component becomes similar
to the one from the accretion disk.

  * •

the emission from the dusty torus is thermal and can be well described with a
single-temperature black-body distribution with temperature $T\simeq 1000$ K .
The torus location can also be parametrized as a function of the disk
luminosity $r_{torus}=2.5\times 10^{18}\sqrt{L_{disk}/10^{45}~{}\textrm{erg
s}^{-1}}$ cm (Sikora et al., 2009). Similarly to the BLR aperture angle, the
aperture angle of the torus is also not well known. It is commonly assumed to
be of around $\pi/4$. As with the BLR photon field, the blob of plasma sees a
constant photon density from the torus for $r\ll r_{torus}$, and then a
rapidly decreasing photon field while crossing $r_{torus}$ and leaving the
torus behind (Błażejowski et al., 2000). An example of an EIC modeling of the
FSRQ 3C279 is shown in Figure 3.

  * •

A ubiquitous photon field is the one from the Cosmic Microwave Background
(CMB): electrons in the jet can also scatter these photons to higher energies.
This EIC emission has the important properties of being the only one dependent
on the redshift of the source. EIC over the CMB is not so efficient when the
emitting region is located close to the SMBH, because there are always
stronger photon fields. But it can become an important radiative process to
explain X-ray and $\gamma$-ray emission from SMBH jets at larger scales
((Tavecchio et al., 2000), but see (Meyer et al., 2017)). lSince the CMB
intensity scales as $(1+z)^{4}$, this emission process can play a more
important role for high-redshift quasars (Schwartz, 2002).

  * •

as we discussed before, structured jets can explain the synchrotron emission
from radio up to optical-X-rays. The emission from the other regions of the
jet can thus serve as target photon field for EIC scattering. In
configurations such as the spine-sheath scenario (a fast inner jet surrounded
by a slower layer), the emission from the external jet can be boosted in the
$\gamma$-ray emitting region due to the relative motion of the two components
(Ghisellini et al., 2005; Hervet et al., 2015; MacDonald et al., 2015).

Figure 3: Spectral energy distribution of the flat spectrum radio quasar
3C279, showing in gray archival data from the SSDC SED builder
(https://tools.ssdc.asi.it/SED/). One-zone external-inverse-Compton modeling
of the SED, with in red synchrotron emission by electrons; in green
synchrotron-self-Compton emission; in blue EIC emission over dust photons; and
in pink EIC emission over Ly$\alpha$ photons. Parameters are $\delta=30$,
$R=1.3\times 10^{16}\ \textrm{cm}$, $B^{\prime}=1.5\ \textrm{G}$,
$\gamma^{\prime}_{min}=50$, $\gamma^{\prime}_{break}=350$,
$\gamma^{\prime}_{Max}=10^{4}$, $n_{1}=2.5$, $n_{2}=3.5$, $L_{e}=7\times
10^{44}\ \textrm{erg s}^{-1}$, $u^{\prime}_{Ly\alpha}=1\ \textrm{erg
cm}^{-3}$, $u^{\prime}_{torus}=0.7\ \textrm{erg cm}^{-3}$

####  2.2.4 Energy Budget of the Emitting Region

An important characteristic of all radiative models is its energy budget: how
much energy is required to reproduce the observations? and how is energy
distributed among the various components? The total jet power (in erg s-1) can
be expressed as an explicit function of the various energy densities as

$L_{jet}=2\pi R^{2}c\Gamma^{2}(u^{\prime}_{B}+u^{\prime}_{e}+u^{\prime}_{ph})$
(24)

where the factor of 2 takes into account that there are always two jets.

Equipartition between the various radiative components is always attractive
especially in physical systems at equilibrium, because it directly provides a
minimum power solution. It is thus interesting to note that SSC solutions for
HBLs are not near equipartition, the emitting region being always particle
dominated (Kino et al., 2002; Tavecchio and Ghisellini, 2016). On the other
hand, EIC solutions for FSRQs show near-equipartition conditions between
particles, magnetic field, and external photons (Dermer et al., 2014).

###  2.3 Electron-Positron Pair Production

The last leptonic process that is relevant to calculate the radiative output
of SMBH jets is pair-production via photon–photon interaction (Gould and
Schréder, 1967)

$\gamma+\gamma\rightarrow e^{+}+e^{-}$ (25)

This process is the inverse of pair annihilation, in which one electron and
one positron annihilate into 0.511 MeV photons. Pair-production can occur only
if the energy of the two photons is higher than the rest mass of the leptons
that are produced, which is
$\varepsilon_{1}\varepsilon_{2}=h\nu_{1}h\nu_{2}/m_{e}^{2}c^{4}>1$. The cross-
section of the interaction is (Coppi and Blandford, 1990)

$\sigma_{\gamma\textrm{-}\gamma}(s)=\frac{3\sigma_{T}}{16}(1-s^{2})\left[2s(s^{2}-2)+(3-s^{4})\ln\left(\frac{1+s}{1-s}\right)\right]$
(26)

where $s=\varepsilon_{1}\varepsilon_{2}(1-\cos{\vartheta})/2$, and $\vartheta$
the interaction angle. The cross-section can be reasonably well approximated
as

$\sigma_{\gamma\textrm{-}\gamma}(\varepsilon_{1},\varepsilon_{2})\simeq\frac{2}{3}\sigma_{T}\delta\left(\varepsilon_{2}-\frac{2}{\varepsilon_{1}}\right)\frac{1}{\varepsilon_{1}}$
(27)

This approximation is known as the $\delta$-approximation. The pair-production
process is more efficient when the energy of the two photons are inversely
proportional. The two photons disappear in the interaction, and while for low-
energy photon fields the disappearance of a single photon is negligible, for
the high-energy photon field a single photon can contribute significantly to
the overall energy flux. For $\gamma$-rays, pair-production can result in
significant absorption, which can be expressed via the following absorption
coefficient (similarly to Equation (5); note that the following equation is
correct only if an angle-averaged cross-section is used, otherwise an
additional integration over the angle is required):

$\mu^{\prime}_{\gamma\textrm{-}\gamma}(\nu^{\prime})=\int_{\varepsilon^{\prime}_{min}}^{\varepsilon^{\prime}_{Max}}d\varepsilon^{\prime}\
n_{ph}^{\prime}(\varepsilon^{\prime})\sigma_{\gamma\textrm{-}\gamma}(\varepsilon^{\prime},h\nu^{\prime}/m_{e}c^{2})$
(28)

where $n_{ph}^{\prime}(\varepsilon^{\prime})$ is the number density of the
low-energy photon field, which can be computed from the intensity

$n_{ph}^{\prime}(\varepsilon^{\prime})=\frac{4\pi}{hc}\frac{I^{\prime}_{\varepsilon^{\prime}}(\varepsilon^{\prime})}{\varepsilon^{\prime}}$
(29)

The most relevant $\gamma$-$\gamma$ absorption term for $\gamma$-ray blazars
is external to the source, and it is represented by the integrated emission
from all galaxies (stars and dust) that permeates the Universe, the
Extragalactic Background Light (EBL). The EBL acts as an absorber for TeV
photons in the Universe, and this absorption effect increases with the
distance from Earth and with the energy of the $\gamma$-ray photon (Salamon
and Stecker, 1998). The result is a softening of the spectrum of TeV blazars,
and an horizon for TeV astronomy. The TeV horizon is defined as the energy-
dependent redshift for which the opacity $\tau_{EBL}=1$. Although the EBL
significantly hinders TeV observations of distant blazars (the most distant
source observed in the TeV band with Cherenkov telescopes is at z = 0.944,
although detected only up to 0.175 TeV (Ahnen et al., 2016)), this absorption
effect can be used to measure the EBL itself, opening up the path to
$\gamma$-ray cosmology (Biteau and Williams, 2015). The pairs produced in the
interaction with the EBL are not lost, and can produce further emission along
the line of sight mainly via inverse-Compton scattering over CMB photons,
although energy losses from plasma instabilities are also investigated
(Broderick et al., 2012). This emission can be detected both as a new
radiative component in the GeV band (Plaga, 1995), or as extended gamma-ray
pair halos (Aharonian et al., 1994), both modulated by the strength of the
magnetic field seen by the pairs. None of these features has been detected as
per today (H. E. S. S. Collaboration et al., 2014; Archambault et al., 2017;
Ackermann et al., 2018), and these observations can be used to put constraints
on the inter-galactic magnetic field (Neronov and Vovk, 2010; Taylor et al.,
2011).

Going back to the inverse-Compton scattering in SMBH jets, all photon fields
which are inverse-Compton scattered to $\gamma$-ray energies also act as an
absorber to the same $\gamma$-ray photons via pair-production. Let us start by
looking at the SSC case. Under the $\delta$-function approximation for
$\sigma_{\gamma\textrm{-}\gamma}$, and by forcing an equality to the causality
relation (Equation (22)), it is possible to write
$\tau_{\gamma\textrm{-}\gamma}<1$, where the pair-production opacity depends
only on the Doppler factor $\delta_{SSC}$ (and physical observables) (Dondi
and Ghisellini, 1995), and it is thus possible to provide a lower limit on
$\delta_{SSC}$

$\delta_{SSC}\geq\left[\frac{\sigma_{T}}{5hc^{2}}d_{L}^{2}(1+z)^{2\alpha}\
\frac{F_{\nu}(\nu_{0})}{\tau_{var}}\right]^{1/(4+2\alpha)}$ (30)

where $\nu_{0}=\left(1.6\times 10^{40}/\nu_{\gamma}\right)$ Hz is the
frequency of the low-energy photons that act as primary absorber for photons
at $\nu_{\gamma}$, and $\alpha$ is the index of the flux density $F_{\nu}$ at
$\nu_{0}$. For typical blazars parameters, Equation (30) results in
$\delta_{SSC}\geq 5-10$. The detection of TeV photons from HBLs can thus be
used to put a strong limit on $\delta$, constraining further the parameter
space of the SSC model.

####  2.3.1 On the Location of the $\gamma$-ray Emitting Region in FSRQs

For the EIC case, the opacity to $\gamma$-$\gamma$ pair production depends on
the location $r$ of the emitting region with respect to the external fields.
The opacity is expressed as the integral from $r$ outwards of the absorption
coefficient:

$\tau_{\gamma\textrm{-}\gamma}(\nu^{\prime})=\int_{r}^{\infty}dr^{\prime}\
\mu^{\prime}_{\gamma\textrm{-}\gamma}(\nu^{\prime})$ (31)

The presence of the BLR is very relevant for opacity calculations: if
$r<r_{BLR}$, radiation from the BLR will significantly absorb $\gamma$-ray
photons. Assuming that the most relevant emission line is Ly$\alpha$ (13.6
eV), this will cause an absorption feature at $\simeq$100 GeV. Detailed
calculations including all relevant BLR lines and the thermal continuum from
the disk scattered by the BLR (Ghisellini and Madau, 1996; Donea and
Protheroe, 2003; Liu and Bai, 2006; Ghisellini and Tavecchio, 2009; Finke,
2016) show that the very detection of $\gamma$-ray photons with energies
greater than 100 GeV from FSRQs is enough to put the location of the emitting
region beyond the BLR. Observations of FSRQs at $E>100$ GeV confirmed this
scenario and proved that, at least for the sources that have been detected at
high energies, the emitting region is at or beyond $r_{BLR}$ (Aleksić et al.,
2011; Nalewajko et al., 2014; Ahnen et al., 2015; Abeysekara et al., 2015;
Costamante et al., 2018; H. E. S. S. Collaboration et al., 2019, 2020).

The dusty torus also leaves an absorption imprint on the $\gamma$-ray
spectrum. Given its typical temperature of $\simeq$1000 K, the pair-production
absorption on the torus photons is expected to be seen in the TeV band in
$\gamma$-ray FSRQs. Due to the simultaneous absorption on the EBL, and to the
fact that FSRQs have an SED peak in the 100-MeV range with very few TeV
photons, such a cut-off has never been observed as per today (Galanti et al.,
2020).

####  2.3.2 Pair Injection

The electron-positron pairs that have such an important role for $\gamma$-ray
absorption in SMBH jet are injected in the emitting region, and can be treated
as a secondary injection of leptons. They will behave exactly as the primary
ones, and will radiate synchrotron and inverse-Compton photons, reaching their
own equilibrium distribution. The big advantage is that while we do not know
the details of the primary accelerator and we simply parametrized it, we have
the information on the energy distribution of the pairs injected via pair-
production (Boettcher and Schlickeiser, 1997):

$\begin{array}[]{l}Q^{\prime}(\gamma^{\prime}_{e})=\frac{3\sigma_{T}}{32}\int_{\gamma^{\prime}_{e}}^{\infty}d\varepsilon_{1}\frac{n_{ph}^{\prime}(\varepsilon_{1})}{\varepsilon^{3}_{1}}\int_{\frac{\varepsilon_{1}}{4\gamma^{\prime}_{e}(\varepsilon_{1}-\gamma^{\prime}_{e})}}^{\infty}d\varepsilon_{2}\frac{n_{ph}^{\prime}(\varepsilon_{2})}{\varepsilon^{2}_{2}}\times\\\
\times\left[\frac{4\varepsilon^{2}_{1}}{\gamma^{\prime}(\varepsilon_{1}-\gamma^{\prime})}\log\left(\frac{4\gamma^{\prime}_{e}\varepsilon_{2}(\varepsilon_{1}-\gamma^{\prime})}{\varepsilon_{1}}\right)-8\varepsilon_{1}\varepsilon_{2}+\frac{2\varepsilon^{2}_{1}(\varepsilon_{1}\varepsilon_{2}-1)}{\gamma^{\prime}(\varepsilon_{1}-\gamma^{\prime})}-\left(1-\frac{1}{\varepsilon_{1}\varepsilon_{2}}\right)\left(\frac{\varepsilon^{2}_{1}}{\gamma^{\prime}(\varepsilon_{1}-\gamma^{\prime})}\right)^{2}\right]\end{array}$
(32)

where $n_{ph}^{\prime}(\varepsilon_{1})$ and
$n_{ph}^{\prime}(\varepsilon_{2})$ are the photon densities of the two photon
fields. The injection term can then be used to compute the steady-state
distribution and the associated synchrotron and inverse-Compton emission
following Equation (7). For single-zone blazar leptonic models, emission by
this secondary leptonic population is usually negligible, hidden by the
primary one.

##  3 Hadronic Radiative Processes

###  3.1 Proton–Synchrotron Emission

Let us now consider the case of a jet in which there are also protons (and
possibly higher-Z nuclei) together with leptons. As with electrons, the proton
energy distribution is defined between $\gamma^{\prime}_{p,min}$ and
$\gamma^{\prime}_{p,Max}$ as a power-law function:

$N^{\prime}_{p}(\gamma^{\prime}_{p})=N^{\prime}_{0,p}\gamma^{\prime-
n_{p}}_{p}$ (33)

Going back to our spherical plasmoid in the jet, the first radiative process
that needs to be studied is synchrotron radiation by protons.With respect to
electrons, the emission from protons is suppressed and needs to be compensated
by increasing either the magnetic field value, the particle density, or both.
It is interesting to test the possibility that proton–synchrotron radiation is
responsible for the high-energy SED component in blazars. Indeed,
proton–synchrotron radiation would have the same index
$\alpha_{p}=(n_{p}-1)/2$ as SSC radiation, assuming that the injection index
of protons is similar to the one of electrons. Following Equation (19)
(replacing the electron mass with the proton mass) it is possible to have a
proton–synchrotron peak frequency at around $10$-$100$ GeV, consistent with
HBL SEDs, assuming $\gamma_{p,Max}^{\prime}\simeq 10^{9}$, $B^{\prime}\simeq
10$–$100$ G, and $\delta=10$–$50$. The proton–synchrotron solution for HBLs
lies thus in a completely different part of the parameter space with respect
to the SSC solution, with a much stronger (a thousand times higher) magnetic
field (Mannheim, 1993; Aharonian, 2000; Mücke and Protheroe, 2001). The number
of free parameters of the proton–synchrotron scenario is much larger than the
SSC one, due to the fact that the proton distribution adds 5 additional free
parameters (assuming an additional $\gamma^{\prime}_{p,break}$ due to cooling,
and the two indexes $n_{p,1/2}$), for a total of 14 model parameters. Such a
model is clearly degenerate, although the number of free parameters can be
considerably reduced by assuming some physically motivated constraints, such
as: self-consistency of break Lorentz factors, co-acceleration of electrons
and protons (and thus equality of $n_{e,1}$ and $n_{p,1}$), constraint on
$\gamma_{p,Max}$ via equation of acceleration and cooling, constraint on $R$
via causality arguments. But even imposing all these constraints the model
remain degenerate, and and we can only constrain allowed regions in the
parameter space (Cerruti et al., 2015).

###  3.2 Proton-Photon Interactions

Protons in the jet interact with the low-energy photon fields, both internal
(i.e., the electron synchrotron radiation) and external. But contrarily to
electrons, for which the dominant interaction is inverse-Compton scattering,
the dominant processes for protons are photo-meson production and
Bethe–Heitler pair-production. Inverse-Compton scattering for protons does
happen, but it can be safely neglected for the modeling of AGN emission.

####  3.2.1 Photo-Meson Production

The photo-meson (or photo-pion) process is the production of pions in proton-
photon interactions:

$\begin{array}[]{l}p+\gamma\rightarrow p^{\prime}+\pi^{0}\\\
p+\gamma\rightarrow n+\pi^{+}\\\ p+\gamma\rightarrow
p^{\prime}+\pi^{+}+\pi^{-}\\\ \end{array}$ (34)

The photo-meson cross-section is dominated by resonances ($\Delta$ and $N$
baryons), direct pion production and multiple pion production. For
astrophysical applications, the relevant information is the production of
pions, both neutral and charged. In order to produce pions, the photon energy
in the proton frame has to be higher than
$m_{\pi}c^{2}+(m_{\pi}/2m_{p})c^{2}\simeq 145$ MeV. In the emitting-region
frame, $\epsilon^{\prime}\gtrsim(145/\gamma^{\prime}_{p})$ MeV, meaning that
protons with Lorentz factor $\gamma^{\prime}_{p}\geq 10^{7}$ can photo-meson-
produce over ultra-violet photons.

Pions then decay primarily ($\sim 100\%$ branching ratio) as

$\begin{array}[]{l}\pi^{0}\rightarrow 2\gamma\\\
\pi^{+}\rightarrow\mu^{+}+\nu_{\mu}\rightarrow
e^{+}+\nu_{e}+\bar{\nu}_{\mu}+\nu_{\mu}\\\
\pi^{-}\rightarrow\mu^{-}+\bar{\nu}_{\mu}\rightarrow
e^{-}+\bar{\nu}_{e}+\nu_{\mu}+\bar{\nu}_{\mu}\\\ \end{array}$ (35)

Photo-meson production has thus one key property: neutrinos are produced
together with photons, and they can escape the emitting region without
suffering from any additional absorption or energy losses. Detection of
neutrinos from an AGN is thus a smoking-gun signal for the presence of
relativistic protons in the jet. This implies that AGNs can accelerate protons
to high energies, and can be responsible for the flux of particles that we
call, once they reach Earth, cosmic rays. Joint neutrino and photon
observations of AGNs can be used to constrain three key parameters of the
model: how many protons are accelerated, what is their energy distribution,
and what is their maximum energy.

Together with the neutrinos, photo-meson production injects photons and
electron/positron pairs in the emitting region. The spectrum of the injected
pions has a maximum energy equal to the proton maximum energy but reduced by
the inelasticity $\kappa\simeq 0.2$. The maximum energy of photons from the
$\pi^{0}$ decay can be as high as $0.1\gamma^{\prime}_{p}$. For protons with
$\gamma^{\prime}_{p}=10^{7}$, it translates into photons with $E\simeq 10$
PeV. These photons do not reach us: they are absorbed via pair-production both
in the jet, and in the path to Earth. The absorption occurring in the source
produces an electron/positron pair that starts radiating (and cooling). Given
the high energy of the first photon, the synchrotron radiation by these
secondary leptons is at energies high enough that the synchrotron photons can
suffer a second pair-production, and so on. Photons from the $\pi^{0}$ decay
can thus trigger a synchrotron-supported pair-cascade in the emitting region
(or, if the soft photon field energy density is larger than the magnetic one,
an inverse-Compton supported pair-cascade). Emission from the cascade can be
computed calculating the pair distribution at equilibrium and their
synchrotron radiation. When $\tau_{\gamma\textrm{-}\gamma}\gg 1$ the cascade
is saturated, and the resulting synchrotron emission is a power-law with index
$\alpha=1$, or in terms of energy flux $p=0$: the energy transferred from
protons to pions and then to high-energy photons, is then redistributed to low
energies.

Electrons and positrons produced in the decay of charged pions go through a
process similar to the photons from $\pi^{0}$: they have extremely high-
energy, and their synchrotron radiation is also in the PeV band. They thus
trigger as well a pair-cascade, which behaves exactly like the one from
$\pi^{0}$.

Another important emission process from photo-meson interactions is the
synchrotron radiation by $\mu^{\pm}$. Although not stable leptons, muons in
the jet can survive long enough to radiate synchrotron photons before decaying
to electrons/positrons. To compute the muon energy distribution at
equilibrium, the muon decay term must be added to Equation (7) as follows:

$\frac{\partial}{\partial
t}N^{\prime}_{\mu}(\gamma^{\prime}_{\mu},t)=\frac{\partial}{\partial\gamma^{\prime}_{\mu}}\left[\gamma^{\prime}_{\mu}\frac{N^{\prime}_{\mu}(\gamma^{\prime}_{\mu},t)}{\tau_{syn}(\gamma^{\prime}_{\mu})}\right]-\frac{N^{\prime}_{\mu}(\gamma^{\prime}_{\mu},t)}{\tau_{ad}}-\frac{N^{\prime}_{\mu}(\gamma^{\prime}_{\mu},t)}{\gamma^{\prime}_{\mu}\tau_{dec}}+Q^{\prime}_{\mu}(\gamma^{\prime}_{\mu})$
(36)

where $\tau_{dec}=2.2\times 10^{-6}$ seconds. The corresponding integral
solution for the muon distribution at equilibrium is:

$N^{\prime}_{\mu}(\gamma^{\prime}_{\mu})=\exp\left[-\frac{\gamma^{\prime}_{\mu,break}}{\gamma^{\prime}_{\mu}}-\frac{\gamma^{\prime}_{\mu,break}\tau_{ad}}{2\gamma^{\prime
2}_{\mu}\tau_{dec}}\right]\frac{\gamma^{\prime}_{\mu,break}\tau_{ad}}{\gamma^{\prime
2}_{\mu}}\int_{\gamma^{\prime}_{\mu}}^{\infty}d\zeta
Q^{\prime}_{\mu}(\zeta)\exp\left[\frac{\gamma^{\prime}_{\mu,break}}{\zeta}+\frac{\gamma^{\prime}_{\mu,break}\tau_{ad}}{2\zeta^{2}\tau_{dec}}\right]$
(37)

Synchrotron radiation by muons emerges as a third radiative component, at
higher energies compared to the synchrotron radiation by the parent protons.
In some parts of the parameter space it can represent the dominant radiative
emission in the TeV band. Another important consequence is that every time the
muon synchrotron radiation is important, it means that muons can cool
efficiently before decaying into electrons/positrons, and this energy loss has
to be included when calculating their injection in the emitting region.

An interesting property of photo-meson interactions is the creation of
neutrons, which can escape the emitting region without interacting with
magnetic fields, nor producing Bethe–Heitler pairs (see Section 3.2.2). These
neutrons can transfer a significant amount of energy at much larger distances
downstream the jet (Begelman et al., 1990; Atoyan and Dermer, 2003). In
presence of dense photon fields they can trigger additional photo-meson
production, or they can decay into protons (with a life-time of
$\gamma_{n}\times 880$ seconds) and radiate again synchrotron photons in the
presence of magnetic fields. This “neutral beam” model has drawn attention for
its possibility to naturally produce two separate emitting regions in the jet,
physically separated but casually connected (Dermer et al., 2012; Zhang et
al., 2020).

The main problem we are facing with photo-meson production is that the exact
computation of the energy distribution of all secondary particles injected in
the emitting region is a complex numerical task. The best approach is to
compute them via Monte-Carlo simulations (Mücke et al., 2000). Several
approaches to fit the results of Monte-Carlo simulations and provide a
parametrization of the injected secondary particles over a wide parameter
space are particularly useful (Kelner and Aharonian, 2008; Hümmer et al.,
2010).

####  3.2.2 Bethe–Heitler Pair-Production

Bethe–Heitler pair-production is the following process

$p+\gamma\rightarrow p^{\prime}+e^{+}+e^{-}$ (38)

It is a process which is in competition with the photo-meson production,
although it happens at lower energies. The threshold for this interaction is
lower than the photo-meson one by a factor $m_{e}/m_{\pi}\simeq 0.004$. In
terms of proton Lorentz factors, protons with $\gamma^{\prime}_{p}\geq 10^{5}$
can pair-produce over $10$ eV photons.

For low-energy protons, the Bethe–Heitler pair production is the main proton-
photon interaction, but as soon as the proton energy becomes greater than the
energy threshold for photo-meson production, the latter takes over. As
expected, pairs injected via this process have a lower energy compared to
pairs from the photo-meson (Petropoulou and Mastichiadis, 2015). The energy
distribution of the injected Bethe–Heitler pairs can be computed analytically
(Blumenthal, 1970; Chodorowski et al., 1992; Kelner and Aharonian, 2008), or
using a Monte-Carlo approach (Protheroe and Johnson, 1996; Mastichiadis et
al., 2005). Their fate is the same as all other leptons in the emitting
region: they radiate synchrotron and inverse-Compton photons, and if their
energy is high enough (and if the low-energy photon field is dense enough)
they can trigger a pair-cascade.

###  3.3 Hadronic and Lepto-Hadronic Models

Blazar hadronic models are a large family of models which share the presence
of a parent population of protons that are accelerated in the SMBH jet. But as
discussed above, from the electromagnetic point-of-view the only true hadronic
radiative process is proton–synchrotron emission. Solutions in which the high-
energy emission is produced by secondary leptons produced in p-$\gamma$
interactions are usually referred to as lepto-hadronic models.

Pure hadronic (proton–synchrotron) models can fit blazar SEDs (Mannheim, 1993;
Aharonian, 2000; Mücke and Protheroe, 2001; Böttcher et al., 2013), but they
face a major problem: the proton density required to fit the data is rather
high, and when expressed in terms of $L_{p}$ it gets close or clearly much
higher than the Eddington luminosity of the SMBH (Sikora et al., 2009).
Although $L_{Edd}$ should be treated as order-of-magnitude estimate for the
available accretion power, if the power in the jet required to fit the data is
super-Eddington by several orders of magnitudes it means that either the
relation between accretion and jet needs to be revisited, or, more likely,
that this particular solution is not viable. This energy-crisis for hadronic
models is particularly true for FSRQs (Zdziarski and Bottcher, 2015). For low-
luminosity HBLs and EHBLs it is not the case, and hadronic solutions with
$L_{jet}<L_{Edd}$ can be found (Cerruti et al., 2015; Petropoulou et al.,
2015). Even if the high-energy SED component is dominated by
proton–synchrotron radiation, emission by secondary leptons from p-$\gamma$
interactions can emerge (namely for compact and dense solutions) in the X-rays
and in the TeV band (Mücke et al., 2003; Abdo et al., 2011; Zech et al.,
2017). In terms of energy budget of the emitting region, proton–synchrotron
solutions are out of equipartition, with $u^{\prime}_{B}\geq u^{\prime}_{p}$.
An example of a hadronic modeling of the HBL Mrk 421 is shown in the second
plot of Figure 2.

In the lepto-hadronic solutions, the proton–synchrotron radiative component is
suppressed, and this can be easily achieved by reducing the magnetic field to
$B^{\prime}\leq 1$ G. We are back in the part of the parameter space where the
SSC emission dominates and indeed the simplest lepto-hadronic model is an SSC
solution loaded with relativistic protons that produce secondary leptons via
p-$\gamma$ interactions over the electron synchrotron photon field (this
solution is called one-zone lepto-hadronic model, to highlight the fact that
all radiative mechanisms are coming from a single emitting region in the jet,
and there external photon fields are negligible). The emission by secondary
leptons emerges again in the X-rays (as Bethe–Heitler component) and in the
TeV band (as photo-meson component) (Petropoulou et al., 2015). This kind of
lepto-hadronic solutions are typically more promising in terms of $\nu$ output
than pure hadronic solutions. As with hadronic solutions, lepto-hadronic
solutions also face energetic issues: especially if the goal of the modeling
is to maximize the neutrino output, the required jet power can quickly become
very high. A possible way-out to this issue is to use as target photon field
not the electron synchrotron radiation, but external photons. By increasing
the density of the target photon field, it is possible to achieve the same
p-$\gamma$ output lowering the required proton power. The discovery of the
blazar TXS0506 + 056 as first neutrino blazar candidate (IceCube Collaboration
et al., 2018), indicates that a lepto-hadronic solution over an external field
may be the favored scenario (Ansoldi et al., 2018; Keivani et al., 2018;
Murase et al., 2018; Gao et al., 2019; Cerruti et al., 2019; Righi et al.,
2019). In terms of energy budget of the emitting region, mixed lepto-hadronic
are usually particle dominated, although the result depends on the exact shape
of the proton distribution (a soft/hard distribution has a significant impact
on $u^{\prime}_{p}$, but a marginal effect on the photon and neutrino
emission). An example of a mixed lepto-hadronic modeling of the HBL Mrk 421 is
shown in the third plot of Figure 2.

Hadronic radiative models are of particular interest any time leptonic ones
face difficulties. One of the most interesting cases is the unusual SED of the
nearby radio-galaxy Centaurus A, which shows a unique third radiative
component emerging at $\simeq$100 GeV (H. E. S. S. Collaboration et al.,
2018). This additional component can be naturally explained in a hadronic or
lepto-hadronic scenario (Petropoulou et al., 2014; Fraija, 2014; Sahu et al.,
2012).

Among blazar hadronic radiative models, it is worth mentioning another
potential contribution produced not in the jet, but rather in the path from
the source to the observer. If AGNs accelerate cosmic rays, a hadronic beam
can be launched emitting photons and neutrinos while travelling towards Earth
(Essey and Kusenko, 2010; Murase et al., 2012). So far no evidence for this
hadronic contribution has been detected. For a recent study of propagation
effects comparing both hadronic cascades and electromagnetic cascades (the
ones triggered by the absorption over the EBL) see (Dzhatdoev et al., 2017).

####  3.3.1 Hadronic Models: Time Signatures

From the spectral point of view both leptonic and hadronic radiative models
provide similarly good fit to current electromagnetic observations, but they
can in principle be distinguished from their temporal behavior. Time-dependent
hadronic models require the solution of a system of coupled differential
equations of the kind:

$\frac{\partial N^{\prime}_{X}(t,E)}{\partial
t}=Q^{\prime}_{X}(t,E)-L^{\prime}_{X}(t,E)$ (39)

for each species $X$ (protons, photons, neutrinos, leptons), where
$Q^{\prime}_{X}(t,E)$ and $L^{\prime}_{X}(t,E)$ represents the injection and
loss terms, respectively. From a numerical point of view, hadronic time-
dependent codes are significantly more complex than leptonic ones, and have
been extensively developed only in the recent years (Mastichiadis et al.,
2005; Dimitrakoudis et al., 2012; Weidinger and Spanier, 2015; Diltz et al.,
2015). As with leptonic models they can be tested on data in two main ways:
fitting light-curve profiles, and studying multi-wavelength correlations.

For proton–synchrotron models, although the two SED components are produced by
two distinct particle populations, it is still possible to reproduce the
observed correlations assuming that electrons and protons share the same
acceleration mechanism (Mastichiadis et al., 2013). On the other hand, absence
of correlations can also be reproduced and indeed, one of the best
applications for blazar hadronic models, are the so-called orphan flare,
bright $\gamma$-ray flares without multi-wavelength counterparts, which cannot
be reproduced in a leptonic one-zone scenario (Krawczynski et al., 2004;
Böttcher, 2005; Sahu et al., 2013).

Time-dependent models have also the advantage to allow the study of non-linear
effects in the development of the hadronic pair-cascades, in which photons
from the cascade become targets for pair-production and can result in a unique
pulsating signature (Petropoulou and Mastichiadis, 2012).

###  3.4 Proton–Proton Interactions

The final hadronic process that needs to be taken into account is pion
production in proton–proton interactions. As with the photo-meson production,
a highly relativistic proton interacts with a target low-energy proton. For
astrophysical applications, the only relevant process is the production of
neutral pions $\pi^{0}$ and $\eta$ mesons, which then decay into $\pi^{\pm}$.
The result of the interaction is then

$\begin{array}[]{l}p+p\rightarrow\pi^{0}+Y\\\
p+p\rightarrow\pi^{+}+\pi^{-}+Z\\\ \end{array}$ (40)

where $Y$ and $Z$ are any other particle produced in the interaction. The
decay of charged pions is responsible for the associated neutrino emission. As
with the photo-meson process, the outcome of proton–proton interactions is
better studied via Monte-Carlo simulations, although useful parametrizations
have been developed (Kelner et al., 2006; Kafexhiu et al., 2014). In blazar’s
jets, proton–proton interactions represent usually a subdominant contribution.
In order to have a significant $p$-$p$ output a very proton-loaded jet is
needed (Reynoso et al., 2011). An alternative scenario is represented by
interactions of the jet with obstacles, such as BLR clouds, or stars (Barkov
et al., 2012, 2012; Bosch-Ramon et al., 2012).

##  4 Polarization Signatures

An important observable that has the potential to uniquely determine the
radiation mechanism, as well as the physical properties of the jet, is
polarization. A specific signature of synchrotron radiation is indeed that the
emission is linearly polarized. For a power-law electron energy distribution
with index $n$, the maximum polarization fraction is

$\Pi(n)=\frac{n+1}{n+\frac{7}{3}}$ (41)

which is equal to $0.69$ for $n=2$. Detection of polarized radiation in
blazars is thus a key prove that what is observed from radio to optical is
synchrotron radiation (Pavlidou et al., 2014). The polarization angle can be
used to map the orientation of the magnetic field in the jet. Recent multi-
wavelength campaigns have shown not only that optical polarization is a
characteristic of blazars, but that both the polarization fraction and the
polarization angle show drastic variability, often correlated with flares at
other wavelengths (Blinov et al., 2016). For an updated review of optical
polarization in AGN jets, see (Zhang, 2019) in this special issue. Variability
in the optical polarization can also be used to put constraints on the
geometry of the emitting regions, and in particular on the relative role of
the small plasmoid (the single-zone that emits from optical to gamma-rays) and
the larger jet (Barres de Almeida et al., 2014; MAGIC Collaboration et al.,
2020).

Inverse-Compton scattering is also expected to be polarized, depending on the
polarization of the target photon field (Krawczynski, 2012). X-ray and gamma-
ray polarimetry have the potential to significantly constrain blazar radiative
models: proton–synchrotron radiation is expected to be significantly much more
polarized than SSC radiation or EIC radiation (which is no polarized at all if
the external field is also not polarized). Leptonic and hadronic radiative
processes can thus be distinguished using this unique observable (Zhang and
Böttcher, 2013; Paliya et al., 2018).

##  5 Summary and Perspectives

Relativistic jets from supermassive black holes emit photons over the whole
electromagnetic spectrum, from radio to gamma-rays, and since 2017 we have the
first evidence that they might be bright in neutrinos as well. The last
decades have seen a significant improvement in our knowledge of their gamma-
ray emission, and the high-energy peak of their SED is now constrained as well
as the low-energy one. These advances on the observational side have
significantly constrained our models for the photon emission. The state-of-
the-art view is that the most likely emission process in FSRQs is external-
inverse-Compton scattering over the bright photon fields around the SMBH
environment. If hadronic processes play a role, they should be subdominant in
FSRQs, due to the large amount of power they require. For HBLs the main
radiation mechanism is SSC, or EIC in a structured jet. In this case, pure
hadronic models (proton–synchrotron dominated) cannot be excluded looking at
the jet energy budget. Observations of the first blazar neutrino candidate,
TXS 0506 + 056, favor a mixed scenario with a leptonic dominated SED with
subdominant hadronic components in the form of pair-cascades emerging in the
hard-X-rays and the TeV band. X-ray and TeV observations helped identifying a
population of blazars with extreme peak frequencies (in the hard-X-rays and
the TeV band, respectively), now known as EHBLs: in this case as well it is
not possible to disentangle the different radiative mechanisms, and both
leptonic and hadronic solutions provide good description of the data. In the
upcoming future new observations will provide key constraints on radiative
models, and will significantly increase our understanding of the physics of
relativistic jets. Neutrino astronomy is just in its early years and
additional data from IceCube and its upgrades (Aartsen et al., 2019), ANTARES
(Ageron et al., 2011), and KM3NET (Adrián-Martínez et al., 2016) will provide
new constraints on hadronic emission processes in AGN jets, whether new
neutrino AGNs will be detected, or not. X-ray polarimetry will also impact
significantly the field, providing unique constraints on the radiative
mechanism responsible for the X-ray emission in AGNs. The Imaging X-ray
Polarimetry Explorer (IXPE) satellite (Weisskopf et al., 2016) is expected to
be launched in late 2021. Since they are the most common sources in the gamma-
ray sky, answers to the current open problems and new questions on SMBH jets
physics will come from the next generation gamma-ray observatories, CTA in the
TeV band (Acharya et al., 2013; Cherenkov Telescope Array Consortium et al.,
2019), and future proposed MeV observatories, such as AMEGO (McEnery et al.,
2019).

M. Cerruti has received financial support through the Postdoctoral Junior
Leader Fellowship Programme from la Caixa Banking Foundation, grant n.
LCF/BQ/LI18/11630012. Funding for this work was partially provided by the
Spanish MINECO under project MDM-2014-0369 of ICCUB (Unidad de Excelencia
’María de Maeztu’). Silvia Cerruti is acknowledged for the AGN scheme in
Figure 1.

The authors declare no conflict of interest. The funders had no role in the
design of the study; in the collection, analyses, or interpretation of data;
in the writing of the manuscript, or in the decision to publish the results.

The following abbreviations are used in this manuscript:

AGN | Active Galactic Nucleus
---|---
BLR | Broad Line Region
EBL | Extragalactic Background Light
EHBL | Extremely High-Frequency Peaked BL Lac objects
EIC | External-Inverse-Compton
FR | Fanaroff–Riley (I and II)
FSRQ | Flat Spectrum Radio Quasar
HBL | High-Frequency Peaked BL Lac objects
IBL | Intermediate-Frequency Peaked BL Lac objects
$\Lambda$CDM |  $\Lambda$-Cold-Dark-Matter
LBL | Low-Frequency Peaked BL Lac objects
NLR | Narrow Line Region
SED | Spectral Energy Distribution
SMBH | Supermassive Black Hole
SSC | Synchrotron-Self-Compton

References

## References

  * Bañados et al. (2018) Bañados, E.; Venemans, B.P.; Mazzucchelli, C.; Farina, E.P.; Walter, F.; Wang, F.; Decarli, R.; Stern, D.; Fan, X.; Davies, F.B.; et al.  An 800-million-solar-mass black hole in a significantly neutral Universe at a redshift of 7.5.  Nature 2018, 553, 473–476, doi:10.1038/nature25180.
  * Planck Collaboration et al. (2018) Aghanim, N.; Akrami, Y.; Ashdown, M.; Aumont, J.; Baccigalupi, C.; Ballardini, M.; Banday, A.J.; Barreiro, R.B.; Bartolo, N.; Basak, S.; et al.  Planck 2018 results. VI. Cosmological parameters.  arXiv 2018, arXiv:1807.06209.
  * Event Horizon Telescope Collaboration et al. (2019) Akiyama, K.; Alberdi, A.; Alef, W.; Asada, K.; Azulay, R.; Baczko, A.K.; Ball, D.; Baloković, M.; Barrett, J.; Bintley, D.; et al.  First M87 Event Horizon Telescope Results. I. The Shadow of the Supermassive Black Hole.  Astrophys. J. Lett. 2019, 875, L1, doi:10.3847/2041-8213/ab0ec7.
  * Abramowicz and Fragile (2013) Abramowicz, M.A.; Fragile, P.C.  Foundations of Black Hole Accretion Disk Theory.  Living Rev. Relativ. 2013, 16, 1, doi:10.12942/lrr-2013-1.
  * Jaffe et al. (1993) Jaffe, W.; Ford, H.C.; Ferrarese, L.; van den Bosch, F.; O’Connell, R.W.  A large nuclear accretion disk in the active galaxy NGC4261.  Nature 1993, 364, 213–215, doi:10.1038/364213a0.
  * Greene and Ho (2005) Greene, J.E.; Ho, L.C.  Estimating Black Hole Masses in Active Galaxies Using the H$\alpha$ Emission Line.  Astrophys. J. 2005, 630, 122–129, doi:10.1086/431897.
  * Bentz et al. (2009) Bentz, M.C.; Walsh, J.L.; Barth, A.J.; Baliber, N.; Bennert, V.N.; Canalizo, G.; Filippenko, A.V.; Ganeshalingam, M.; Gates, E.L.; Greene, J.E.; et al.  The Lick AGN Monitoring Project: Broad-line Region Radii and Black Hole Masses from Reverberation Mapping of H$\beta$.  Astrophys. J. 2009, 705, 199–217, doi:10.1088/0004-637X/705/1/199.
  * Bennert et al. (2002) Bennert, N.; Falcke, H.; Schulz, H.; Wilson, A.S.; Wills, B.J.  Size and Structure of the Narrow-Line Region of Quasars.  Astrophys. J. Lett. 2002, 574, L105–L109, doi:10.1086/342420.
  * Netzer (2015) Netzer, H.  Revisiting the Unified Model of Active Galactic Nuclei.  Annu. Rev. Astron. Astrophys. 2015, 53, 365–408, doi:10.1146/annurev-astro-082214-122302.
  * Schmidt (1963) Schmidt, M.  3C 273: A Star-Like Object with Large Red-Shift.  Nature 1963, 197, 1040, doi:10.1038/1971040a0.
  * Bahcall et al. (1997) Bahcall, J.N.; Kirhakos, S.; Saxe, D.H.; Schneider, D.P.  Hubble Space Telescope Images of a Sample of 20 Nearby Luminous Quasars.  Astrophys. J. 1997, 479, 642–658, doi:10.1086/303926.
  * Wilson and Colbert (1995) Wilson, A.S.; Colbert, E.J.M.  The Difference between Radio-loud and Radio-quiet Active Galaxies.  Astrophys. J. 1995, 438, 62, doi:10.1086/175054.
  * Kellermann et al. (1998) Kellermann, K.I.; Vermeulen, R.C.; Zensus, J.A.; Cohen, M.H.  Sub-Milliarcsecond Imaging of Quasars and Active Galactic Nuclei.  Astron. J. 1998, 115, 1295–1318, doi:10.1086/300308.
  * Butcher et al. (1980) Butcher, H.R.; van Breugel, W.; Miley, G.K.  Optical observations of radio jets.  Astrophys. J. 1980, 235, 749–754, doi:10.1086/157677.
  * Marshall et al. (2005) Marshall, H.L.; Schwartz, D.A.; Lovell, J.E.J.; Murphy, D.W.; Worrall, D.M.; Birkinshaw, M.; Gelbord, J.M.; Perlman, E.S.; Jauncey, D.L.  A Chandra Survey of Quasar Jets: First Results.  Astrophys. J. Suppl. 2005, 156, 13–33, doi:10.1086/425578.
  * Abdo et al. (2010) Abdo, A.A.; Ackermann, M.; Ajello, M.; Atwood, W.B.; Baldini, L.; Ballet, J.; Barbiellini, G.; Bastieri, D.; Baughman, B.M.; Bechtol, K.; et al.  Fermi Gamma-Ray Imaging of a Radio Galaxy.  Science 2010, 328, 725, doi:10.1126/science.1184656.
  * De Naurois (2019) De Naurois, M.  Resolving the kpc jet of Centaurus A in TeV gamma-rays.  In Proceedings of the 36th International Cosmic Ray Conference (ICRC2019), Madison, WI, USA, 1 August 2019; Volume 36, p. 657.
  * Blandford and Rees (1978) Blandford, R.D.; Rees, M.J.  Extended and compact extragalactic radio sources: Interpretation and theory.  Phys. Scr. 1978, 17, 265–274, doi:10.1088/0031-8949/17/3/020.
  * Blandford and Königl (1979) Blandford, R.D.; Königl, A.  Relativistic jets as compact radio sources.  Astrophys. J. 1979, 232, 34–48, doi:10.1086/157262.
  * Urry and Padovani (1995) Urry, C.M.; Padovani, P.  Unified Schemes for Radio-Loud Active Galactic Nuclei.  Publ. Astron. Soc. Pac. 1995, 107, 803, doi:10.1086/133630.
  * Angel and Stockman (1980) Angel, J.R.P.; Stockman, H.S.  Optical and infrared polarization of active extragalactic objects.  Annu. Rev. Astron. Astrophys. 1980, 18, 321–361, doi:10.1146/annurev.aa.18.090180.001541.
  * The Fermi-LAT collaboration (2019) The Fermi-LAT Collaboration.  Fermi Large Area Telescope Fourth Source Catalog.  arXiv 2019, arXiv:1902.10045.
  * Shaw et al. (2013) Shaw, M.S.; Romani, R.W.; Cotter, G.; Healey, S.E.; Michelson, P.F.; Readhead, A.C.S.; Richards, J.L.; Max-Moerbeck, W.; King, O.G.; Potter, W.J.  Spectroscopy of the Largest Ever $\gamma$-Ray-selected BL Lac Sample.  Astrophys. J. 2013, 764, 135, doi:10.1088/0004-637X/764/2/135.
  * Shaw et al. (2012) Shaw, M.S.; Romani, R.W.; Cotter, G.; Healey, S.E.; Michelson, P.F.; Readhead, A.C.S.; Richards, J.L.; Max-Moerbeck, W.; King, O.G.; Potter, W.J.  Spectroscopy of Broad-line Blazars from 1LAC.  Astrophys. J. 2012, 748, 49, doi:10.1088/0004-637X/748/1/49.
  * Fanaroff and Riley (1974) Fanaroff, B.L.; Riley, J.M.  The morphology of extragalactic radio sources of high and low luminosity.  Mon. Not. R. Astron. Soc. 1974, 167, 31P–36P, doi:10.1093/mnras/167.1.31P.
  * Lister et al. (2013) Lister, M.L.; Aller, M.F.; Aller, H.D.; Homan, D.C.; Kellermann, K.I.; Kovalev, Y.Y.; Pushkarev, A.B.; Richards, J.L.; Ros, E.; Savolainen, T.  MOJAVE. X. Parsec-scale Jet Orientation Variations and Superluminal Motion in Active Galactic Nuclei.  Astron. J. 2013, 146, 120, doi:10.1088/0004-6256/146/5/120.
  * Hartman et al. (1992) Hartman, R.C.; Bertsch, D.L.; Fichtel, C.E.; Hunter, S.D.; Kanbach, G.; Kniffen, D.A.; Kwok, P.W.; Lin, Y.C.; Mattox, J.R.; Mayer-Hasselwander, H.A.; et al.  Detection of High-Energy Gamma Radiation from Quasar 3C 279 by the EGRET Telescope on the Compton Gamma Ray Observatory.  Astrophys. J. Lett. 1992, 385, L1, doi:10.1086/186263.
  * Punch et al. (1992) Punch, M.; Akerlof, C.W.; Cawley, M.F.; Chantell, M.; Fegan, D.J.; Fennell, S.; Gaidos, J.A.; Hagan, J.; Hillas, A.M.; Jiang, Y.; et al.  Detection of TeV photons from the active galaxy Markarian 421.  Nature 1992, 358, 477–478, doi:10.1038/358477a0.
  * Padovani and Giommi (1995) Padovani, P.; Giommi, P.  The Connection between X-ray- and Radio-selected BL Lacertae Objects.  Astrophys. J. 1995, 444, 567, doi:10.1086/175631.
  * Abdo et al. (2010) Abdo, A.A.; Ackermann, M.; Agudo, I.; Ajello, M.; Aller, H.D.; Aller, M.F.; Angelakis, E.; Arkharov, A.A.; Axelsson, M.; Bach, U.; et al.  The Spectral Energy Distribution of Fermi Bright Blazars.  Astrophys. J. 2010, 716, 30–70, doi:10.1088/0004-637X/716/1/30.
  * Laurent-Muehleisen et al. (1999) Laurent-Muehleisen, S.A.; Kollgaard, R.I.; Feigelson, E.D.; Brinkmann, W.; Siebert, J.  The RGB Sample of Intermediate BL Lacertae Objects.  Astrophys. J. 1999, 525, 127–143, doi:10.1086/307881.
  * Costamante et al. (2001) Costamante, L.; Ghisellini, G.; Giommi, P.; Tagliaferri, G.; Celotti, A.; Chiaberge, M.; Fossati, G.; Maraschi, L.; Tavecchio, F.; Treves, A.; et al.  Extreme synchrotron BL Lac objects. Stretching the blazar sequence.  Astron. Astrophys. 2001, 371, 512–526, doi:10.1051/0004-6361:20010412.
  * Bonnoli et al. (2015) Bonnoli, G.; Tavecchio, F.; Ghisellini, G.; Sbarrato, T.  An emerging population of BL Lacs with extreme properties: Towards a class of EBL and cosmic magnetic field probes?  Mon. Not. R. Astron. Soc. 2015, 451, 611–621, doi:10.1093/mnras/stv953.
  * Aharonian et al. (2007) Aharonian, F.; Akhperjanian, A.G.; Barres de Almeida, U.; Bazer-Bachi, A.R.; Behera, B.; Beilicke, M.; Benbow, W.; Bernlöhr, K.; Boisson, C.; Bolz, O.; et al.  New constraints on the mid-IR EBL from the HESS discovery of VHE $\gamma$-rays from 1ES 0229+200.  Astron. Astrophys. 2007, 475, L9–L13, doi:10.1051/0004-6361:20078462.
  * Foffano et al. (2019) Foffano, L.; Prandini, E.; Franceschini, A.; Paiano, S.  A new hard X-ray-selected sample of extreme high-energy peaked BL Lac objects and their TeV gamma-ray properties.  Mon. Not. R. Astron. Soc. 2019, 486, 1741–1762, doi:10.1093/mnras/stz812.
  * Biteau et al. (2020) Biteau, J.; Prandini, E.; Costamante, L.; Lemoine, M.; Padovani, P.; Pueschel, E.; Resconi, E.; Tavecchio, F.; Taylor, A.; Zech, A.  Progress in unveiling extreme particle acceleration in persistent astrophysical jets.  Nat. Astron. 2020, 4, 124–131, doi:10.1038/s41550-019-0988-4.
  * Costamante (2020) Costamante, L.  TeV-peaked candidate BL Lac objects.  Mon. Not. R. Astron. Soc. 2020, 491, 2771–2778, doi:10.1093/mnras/stz3018.
  * Fossati et al. (1998) Fossati, G.; Maraschi, L.; Celotti, A.; Comastri, A.; Ghisellini, G.  A unifying view of the spectral energy distributions of blazars.  Mon. Not. R. Astron. Soc. 1998, 299, 433–448, doi:10.1046/j.1365-8711.1998.01828.x.
  * Meyer et al. (2011) Meyer, E.T.; Fossati, G.; Georganopoulos, M.; Lister, M.L.  From the Blazar Sequence to the Blazar Envelope: Revisiting the Relativistic Jet Dichotomy in Radio-loud Active Galactic Nuclei.  Astrophys. J. 2011, 740, 98, doi:10.1088/0004-637X/740/2/98.
  * Ghisellini et al. (2017) Ghisellini, G.; Righi, C.; Costamante, L.; Tavecchio, F.  The Fermi blazar sequence.  Mon. Not. R. Astron. Soc. 2017, 469, 255–266, doi:10.1093/mnras/stx806.
  * Padovani et al. (2003) Padovani, P.; Perlman, E.S.; Landt, H.; Giommi, P.; Perri, M.  What Types of Jets Does Nature Make? A New Population of Radio Quasars.  Astrophys. J. 2003, 588, 128–142, doi:10.1086/373899.
  * Nieppola et al. (2008) Nieppola, E.; Valtaoja, E.; Tornikoski, M.; Hovatta, T.; Kotiranta, M.  Blazar sequence—An artefact of Doppler boosting.  Astron. Astrophys. 2008, 488, 867–872, doi:10.1051/0004-6361:200809716.
  * Padovani et al. (2012) Padovani, P.; Giommi, P.; Rau, A.  The discovery of high-power high synchrotron peak blazars.  Mon. Not. R. Astron. Soc. 2012, 422, L48–L52, doi:10.1111/j.1745-3933.2012.01234.x.
  * Giommi et al. (2012) Giommi, P.; Padovani, P.; Polenta, G.; Turriziani, S.; D’Elia, V.; Piranomonte, S.  A simplified view of blazars: Clearing the fog around long-standing selection effects.  Mon. Not. R. Astron. Soc. 2012, 420, 2899–2911, doi:10.1111/j.1365-2966.2011.20044.x.
  * Böttcher and Dermer (2002) Böttcher, M.; Dermer, C.D.  An Evolutionary Scenario for Blazar Unification.  Astrophys. J. 2002, 564, 86–91, doi:10.1086/324134.
  * Ghisellini and Tavecchio (2008) Ghisellini, G.; Tavecchio, F.  The blazar sequence: A new perspective.  Mon. Not. R. Astron. Soc. 2008, 387, 1669–1680, doi:10.1111/j.1365-2966.2008.13360.x.
  * Finke (2013) Finke, J.D.  Compton Dominance and the Blazar Sequence.  Astrophys. J. 2013, 763, 134, doi:10.1088/0004-637X/763/2/134.
  * Blandford et al. (2019) Blandford, R.; Meier, D.; Readhead, A.  Relativistic Jets from Active Galactic Nuclei.  Annu. Rev. Astron. Astrophys. 2019, 57, 467–509, doi:10.1146/annurev-astro-081817-051948.
  * Hovatta and Lindfors (2020) Hovatta, T.; Lindfors, E.  Relativistic Jets of Blazars.  arXiv 2020, arXiv:2003.06322.
  * Romero et al. (2017) Romero, G.E.; Boettcher, M.; Markoff, S.; Tavecchio, F.  Relativistic Jets in Active Galactic Nuclei and Microquasars.  Space Sci. Rev. 2017, 207, 5–61, doi:10.1007/s11214-016-0328-2.
  * Kirk et al. (2000) Kirk, J.G.; Guthmann, A.W.; Gallant, Y.A.; Achterberg, A.  Particle Acceleration at Ultrarelativistic Shocks: An Eigenfunction Method.  Astrophys. J. 2000, 542, 235–242, doi:10.1086/309533.
  * Sironi et al. (2015) Sironi, L.; Petropoulou, M.; Giannios, D.  Relativistic jets shine through shocks or magnetic reconnection?  Mon. Not. R. Astron. Soc. 2015, 450, 183–191, doi:10.1093/mnras/stv641.
  * Blumenthal and Gould (1970) Blumenthal, G.R.; Gould, R.J.  Bremsstrahlung, Synchrotron Radiation, and Compton Scattering of High-Energy Electrons Traversing Dilute Gases.  Rev. Mod. Phys. 1970, 42, 237–271, doi:10.1103/RevModPhys.42.237.
  * Rybicki and Lightman (1979) Rybicki, G.B.; Lightman, A.P.  Radiative Processes in Astrophysics; Wiley: Hoboken, NJ, USA, 1979.
  * Longair (1994) Longair, M.S.  High Energy Astrophysics. Volume 2. Stars, the Galaxy and the Interstellar Medium ; Cambridge University Press: Cambridge, UK, 1994; Volume 2.
  * Dermer and Menon (2009) Dermer, C.D.; Menon, G.  High Energy Radiation from Black Holes: Gamma Rays, Cosmic Rays, and Neutrinos; Princeton University Press: Princeton, NJ, USA, 2009.
  * Bloom and Marscher (1996) Bloom, S.D.; Marscher, A.P.  An Analysis of the Synchrotron Self-Compton Model for the Multi—Wave Band Spectra of Blazars.  Astrophys. J. 1996, 461, 657, doi:10.1086/177092.
  * Inoue and Takahara (1996) Inoue, S.; Takahara, F.  Electron Acceleration and Gamma-Ray Emission from Blazars.  Astrophys. J. 1996, 463, 555, doi:10.1086/177270.
  * Kardashev (1962) Kardashev, N.S.  Nonstationarity of Spectra of Young Sources of Nonthermal Radio Emission.  Sov. Astron. 1962, 6, 317.
  * Saugé and Henri (2004) Saugé, L.; Henri, G.  TeV Blazar Gamma-Ray Emission Produced by a Cooling Pileup Particle Energy Distribution Function.  Astrophys. J. 2004, 616, 136–146, doi:10.1086/424905.
  * Sari and Esin (2001) Sari, R.; Esin, A.A.  On the Synchrotron Self-Compton Emission from Relativistic Shocks and Its Implications for Gamma-Ray Burst Afterglows.  Astrophys. J. 2001, 548, 787–799, doi:10.1086/319003.
  * Tavecchio et al. (2010) Tavecchio, F.; Ghisellini, G.; Ghirlanda, G.; Foschini, L.; Maraschi, L.  TeV BL Lac objects at the dawn of the Fermi era.  Mon. Not. R. Astron. Soc. 2010, 401, 1570–1586, doi:10.1111/j.1365-2966.2009.15784.x.
  * Zhang et al. (2012) Zhang, J.; Liang, E.W.; Zhang, S.N.; Bai, J.M.  Radiation Mechanisms and Physical Properties of GeV-TeV BL Lac Objects.  Astrophys. J. 2012, 752, 157, doi:10.1088/0004-637X/752/2/157.
  * Aharonian et al. (2008) Aharonian, F.; Akhperjanian, A.G.; Barres de Almeida, U.; Bazer-Bachi, A.R.; Behera, B.; Beilicke, M.; Benbow, W.; Bernlöhr, K.; Boisson, C.; Borrel, V.; et al.  Discovery of VHE $\gamma$-rays from the high-frequency-peaked BL Lacertae object RGB J0152+017.  Astron. Astrophys. 2008, 481, L103–L107, doi:10.1051/0004-6361:200809603.
  * Aleksić et al. (2015) Aleksić, J.; Ansoldi, S.; Antonelli, L.A.; Antoranz, P.; Babic, A.; Bangale, P.; Barres de Almeida, U.; Barrio, J.A.; Becerra González, J.; Bednarek, W.; et al.  Discovery of very high energy $\gamma$-ray emission from the blazar 1ES 0033+595 by the MAGIC telescopes.  Mon. Not. R. Astron. Soc. 2015, 446, 217–225, doi:10.1093/mnras/stu2024.
  * Archambault et al. (2015) Archambault, S.; Archer, A.; Beilicke, M.; Benbow, W.; Bird, R.; Biteau, J.; Bouvier, A.; Bugaev, V.; Cardenzana, J.V.; Cerruti, M.; et al.  VERITAS Detection of $\gamma$-Ray Flaring Activity From the BL Lac Object 1ES 1727+502 During Bright Moonlight Observations.  Astrophys. J. 2015, 808, 110, doi:10.1088/0004-637X/808/2/110.
  * H. E. S. S. Collaboration et al. (2011) Abramowski, A.; Acero, F.; Aharonian, F.; Akhperjanian, A.G.; Anton, G.; Barnacka, A.; Barres de Almeida, U.; Bazer-Bachi, A.R.; Becherini, Y.; Becker, J.; et al.  Simultaneous multi-wavelength campaign on PKS 2005-489 in a high state.  Astron. Astrophys. 2011, 533, A110, doi:10.1051/0004-6361/201016170.
  * Aliu et al. (2013) Aliu, E.; Archambault, S.; Arlen, T.; Aune, T.; Beilicke, M.; Benbow, W.; Bird, R.; Bouvier, A.; Buckley, J.H.; Bugaev, V.; et al.  Long Term Observations of B2 1215+30 with VERITAS.  Astrophys. J. 2013, 779, 92, doi:10.1088/0004-637X/779/2/92.
  * Baloković et al. (2016) Baloković, M.; Paneque, D.; Madejski, G.; Furniss, A.; Chiang, J.; Ajello, M.; Alexander, D.M.; Barret, D.; Blandford, R.D.; Boggs, S.E.; et al.  Multiwavelength Study of Quiescent States of Mrk 421 with Unprecedented Hard X-ray Coverage Provided by NuSTAR in 2013.  Astrophys. J. 2016, 819, 156, doi:10.3847/0004-637X/819/2/156.
  * Marscher (2014) Marscher, A.P.  Turbulent, Extreme Multi-zone Model for Simulating Flux and Polarization Variability in Blazars.  Astrophys. J. 2014, 780, 87, doi:10.1088/0004-637X/780/1/87.
  * Chiang and Böttcher (2002) Chiang, J.; Böttcher, M.  Synchrotron and Synchrotron Self-Compton Spectral Signatures and Blazar Emission Models.  Astrophys. J. 2002, 564, 92–96, doi:10.1086/324294.
  * Graff et al. (2008) Graff, P.B.; Georganopoulos, M.; Perlman, E.S.; Kazanas, D.  A Multizone Model for Simulating the High-Energy Variability of TeV Blazars.  Astrophys. J. 2008, 689, 68–78, doi:10.1086/592427.
  * Massaro et al. (2004) Massaro, E.; Perri, M.; Giommi, P.; Nesci, R.  Log-parabolic spectra and particle acceleration in the BL Lac object Mkn 421: Spectral analysis of the complete BeppoSAX wide band X-ray data set.  Astron. Astrophys. 2004, 413, 489–503, doi:10.1051/0004-6361:20031558.
  * Tramacere et al. (2009) Tramacere, A.; Giommi, P.; Perri, M.; Verrecchia, F.; Tosti, G.  Swift observations of the very intense flaring activity of Mrk 421 during 2006. I. Phenomenological picture of electron acceleration and predictions for MeV/GeV emission.  Astron. Astrophys. 2009, 501, 879–898, doi:10.1051/0004-6361/200810865.
  * Sol et al. (1989) Sol, H.; Pelletier, G.; Asseo, E.  Two-flow model for extragalactic radio jets.  Mon. Not. R. Astron. Soc. 1989, 237, 411–429, doi:10.1093/mnras/237.2.411.
  * Ghisellini et al. (1985) Ghisellini, G.; Maraschi, L.; Treves, A.  Inhomogeneous synchrotron-self-compton models and the problem of relativistic beaming of BL Lac objects.  Astron. Astrophys. 1985, 146, 204–212.
  * Potter and Cotter (2012) Potter, W.J.; Cotter, G.  Synchrotron and inverse-Compton emission from blazar jets—I. A uniform conical jet model.  Mon. Not. R. Astron. Soc. 2012, 423, 756–765, doi:10.1111/j.1365-2966.2012.20918.x.
  * Archer et al. (2018) Archer, A.; Benbow, W.; Bird, R.; Brose, R.; Buchovecky, M.; Bugaev, V.; Cui, W.; Daniel, M.K.; Falcone, A.; Feng, Q.; et al.  HESS J1943+213: An Extreme Blazar Shining through the Galactic Plane.  Astrophys. J. 2018, 862, 41, doi:10.3847/1538-4357/aacbd0.
  * H. E. S. S. Collaboration et al. (2012) Abramowski, A.; Acero, F.; Aharonian, F.; Akhperjanian, A.G.; Anton, G.; Balzer, A.; Barnacka, A.; Barres de Almeida, U.; Becherini, Y.; Becker, J.; et al.  A multiwavelength view of the flaring state of PKS 2155-304 in 2006\.  Astron. Astrophys. 2012, 539, A149, doi:10.1051/0004-6361/201117509.
  * Ahnen et al. (2016) Ahnen, M.L.; Ansoldi, S.; Antonelli, L.A.; Antoranz, P.; Babic, A.; Banerjee, B.; Bangale, P.; Barres de Almeida, U.; Barrio, J.A.; Becerra González, J.; et al.  Investigating the peculiar emission from the new VHE gamma-ray source H1722+119.  Mon. Not. R. Astron. Soc. 2016, 459, 3271–3281, doi:10.1093/mnras/stw689.
  * Mastichiadis and Kirk (1997) Mastichiadis, A.; Kirk, J.G.  Variability in the synchrotron self-Compton model of blazar emission.  Astron. Astrophys. 1997, 320, 19–25.
  * Kirk et al. (1998) Kirk, J.G.; Rieger, F.M.; Mastichiadis, A.  Particle acceleration and synchrotron emission in blazar jets.  Astron. Astrophys. 1998, 333, 452–458.
  * Böttcher and Chiang (2002) Böttcher, M.; Chiang, J.  X-ray Spectral Variability Signatures of Flares in BL Lacertae Objects.  Astrophys. J. 2002, 581, 127–142, doi:10.1086/344155.
  * Takahashi et al. (1996) Takahashi, T.; Tashiro, M.; Madejski, G.; Kubo, H.; Kamae, T.; Kataoka, J.; Kii, T.; Makino, F.; Makishima, K.; Yamasaki, N.  ASCA Observation of an X-ray/TeV Flare from the BL Lacertae Object Markarian 421.  Astrophys. J. Lett. 1996, 470, L89, doi:10.1086/310302.
  * Kataoka et al. (2000) Kataoka, J.; Takahashi, T.; Makino, F.; Inoue, S.; Madejski, G.M.; Tashiro, M.; Urry, C.M.; Kubo, H.  Variability Pattern and the Spectral Evolution of the BL Lacertae Object PKS 2155-304.  Astrophys. J. 2000, 528, 243–253, doi:10.1086/308154.
  * Zhang et al. (2002) Zhang, Y.H.; Treves, A.; Celotti, A.; Chiappetti, L.; Fossati, G.; Ghisellini, G.; Maraschi, L.; Pian, E.; Tagliaferri, G.; Tavecchio, F.  Four Years of Monitoring Blazar PKS 2155-304 with BeppoSAX: Probing the Dynamics of the Jet.  Astrophys. J. 2002, 572, 762–785, doi:10.1086/340349.
  * Sokolov et al. (2004) Sokolov, A.; Marscher, A.P.; McHardy, I.M.  Synchrotron Self-Compton Model for Rapid Nonthermal Flares in Blazars with Frequency-dependent Time Lags.  Astrophys. J. 2004, 613, 725–746, doi:10.1086/423165.
  * Chiaberge and Ghisellini (1999) Chiaberge, M.; Ghisellini, G.  Rapid variability in the synchrotron self-Compton model for blazars.  Mon. Not. R. Astron. Soc. 1999, 306, 551–560, doi:10.1046/j.1365-8711.1999.02538.x.
  * Klein and Nishina (1929) Klein, O.; Nishina, T.  Über die Streuung von Strahlung durch freie Elektronen nach der neuen relativistischen Quantendynamik von Dirac.  Z. Phys. 1929, 52, 853–868, doi:10.1007/BF01366453.
  * Jones (1968) Jones, F.C.  Calculated Spectrum of Inverse-Compton-Scattered Photons.  Phys. Rev. 1968, 167, 1159–1169, doi:10.1103/PhysRev.167.1159.
  * Moderski et al. (2005) Moderski, R.; Sikora, M.; Coppi, P.S.; Aharonian, F.  Klein-Nishina effects in the spectra of non-thermal sources immersed in external radiation fields.  Mon. Not. R. Astron. Soc. 2005, 363, 954–966, doi:10.1111/j.1365-2966.2005.09494.x.
  * Gould (1979) Gould, R.J.  Compton and synchrotron processes in spherically-symmetric non-thermal sources.  Astron. Astrophys. 1979, 76, 306–311.
  * Tavecchio et al. (1998) Tavecchio, F.; Maraschi, L.; Ghisellini, G.  Constraints on the Physical Parameters of TeV Blazars.  Astrophys. J. 1998, 509, 608–619, doi:10.1086/306526.
  * Finke et al. (2008) Finke, J.D.; Dermer, C.D.; Böttcher, M.  Synchrotron Self-Compton Analysis of TeV X-ray-Selected BL Lacertae Objects.  Astrophys. J. 2008, 686, 181–194, doi:10.1086/590900.
  * Mankuzhiyil et al. (2011) Mankuzhiyil, N.; Ansoldi, S.; Persic, M.; Tavecchio, F.  The Environment and Distribution of Emitting Electrons as a Function of Source Activity in Markarian 421.  Astrophys. J. 2011, 733, 14, doi:10.1088/0004-637X/733/1/14.
  * Cerruti et al. (2013) Cerruti, M.; Boisson, C.; Zech, A.  Constraining the parameter space of the one-zone synchrotron-self-Compton model for GeV-TeV detected BL Lacertae objects.  Astron. Astrophys. 2013, 558, A47, doi:10.1051/0004-6361/201220963.
  * Ahnen et al. (2017) Ahnen, M.L.; Ansoldi, S.; Antonelli, L.A.; Antoranz, P.; Babic, A.; Banerjee, B.; Bangale, P.; Barres de Almeida, U.; Barrio, J.A.; Becerra González, J.; et al.  Multiband variability studies and novel broadband SED modeling of Mrk 501 in 2009.  Astron. Astrophys. 2017, 603, A31, doi:10.1051/0004-6361/201629540.
  * Cerruti et al. (2015) Cerruti, M.; Zech, A.; Boisson, C.; Inoue, S.  A hadronic origin for ultra-high-frequency-peaked BL Lac objects.  Mon. Not. R. Astron. Soc. 2015, 448, 910–927, doi:10.1093/mnras/stu2691.
  * Zacharias and Schlickeiser (2010) Zacharias, M.; Schlickeiser, R.  Blazar synchrotron emission of instantaneously power-law injected electrons under linear synchrotron, non-linear SSC, and combined synchrotron-SSC cooling.  Astron. Astrophys. 2010, 524, A31, doi:10.1051/0004-6361/201015284.
  * Zacharias (2014) Zacharias, M.  Inverse Compton light curves of blazars under non-linear, time-dependent synchrotron-self Compton cooling.  Mon. Not. R. Astron. Soc. 2014, 443, 3001–3021, doi:10.1093/mnras/stu1372.
  * Chen et al. (2011) Chen, X.; Fossati, G.; Liang, E.P.; Böttcher, M.  Time-dependent simulations of multiwavelength variability of the blazar Mrk 421 with a Monte Carlo multizone code.  Mon. Not. R. Astron. Soc. 2011, 416, 2368–2387, doi:10.1111/j.1365-2966.2011.19215.x.
  * Perennes et al. (2020) Perennes, C.; Sol, H.; Bolmont, J.  Modeling spectral lags in active galactic nucleus flares in the context of Lorentz invariance violation searches.  Astron. Astrophys. 2020, 633, A143, doi:10.1051/0004-6361/201936430.
  * Albert et al. (2007) Albert, J.; Aliu, E.; Anderhub, H.; Antoranz, P.; Armada, A.; Baixeras, C.; Barrio, J.A.; Bartko, H.; Bastieri, D.; Becker, J.K.; et al.  Variable Very High Energy $\gamma$-Ray Emission from Markarian 501.  Astrophys. J. 2007, 669, 862–883, doi:10.1086/521382.
  * Abeysekara et al. (2017) Abeysekara, A.U.; Archambault, S.; Archer, A.; Benbow, W.; Bird, R.; Buchovecky, M.; Buckley, J.H.; Bugaev, V.; Cardenzana, J.V.; Cerruti, M.; et al.  A Search for Spectral Hysteresis and Energy-dependent Time Lags from X-ray and TeV Gamma-Ray Observations of Mrk 421.  Astrophys. J. 2017, 834, 2, doi:10.3847/1538-4357/834/1/2.
  * Katarzyński et al. (2005) Katarzyński, K.; Ghisellini, G.; Tavecchio, F.; Maraschi, L.; Fossati, G.; Mastichiadis, A.  Correlation between the TeV and X-ray emission in high-energy peaked BL Lac objects.  Astron. Astrophys. 2005, 433, 479–496, doi:10.1051/0004-6361:20041556.
  * Fossati et al. (2008) Fossati, G.; Buckley, J.H.; Bond, I.H.; Bradbury, S.M.; Carter-Lewis, D.A.; Chow, Y.C.K.; Cui, W.; Falcone, A.D.; Finley, J.P.; Gaidos, J.A.; et al.  Multiwavelength Observations of Markarian 421 in 2001 March: An Unprecedented View on the X-ray/TeV Correlated Variability.  Astrophys. J. 2008, 677, 906–925, doi:10.1086/527311.
  * Acciari et al. (2011) Acciari, V.A.; Aliu, E.; Arlen, T.; Aune, T.; Beilicke, M.; Benbow, W.; Boltuch, D.; Bradbury, S.M.; Buckley, J.H.; Bugaev, V.; et al.  TeV and Multi-wavelength Observations of Mrk 421 in 2006–2008.  Astrophys. J. 2011, 738, 25, doi:10.1088/0004-637X/738/1/25.
  * Aleksić et al. (2015) Aleksić, J.; Ansoldi, S.; Antonelli, L.A.; Antoranz, P.; Babic, A.; Bangale, P.; Barres de Almeida, U.; Barrio, J.A.; Becerra González, J.; Bednarek, W.; et al.  Unprecedented study of the broadband emission of Mrk 421 during flaring activity in March 2010.  Astron. Astrophys. 2015, 578, A22, doi:10.1051/0004-6361/201424811.
  * Furniss et al. (2015) Furniss, A.; Noda, K.; Boggs, S.; Chiang, J.; Christensen, F.; Craig, W.; Giommi, P.; Hailey, C.; Harisson, F.; Madejski, G.; et al.  First NuSTAR Observations of Mrk 501 within a Radio to TeV Multi-Instrument Campaign.  Astrophys. J. 2015, 812, 65, doi:10.1088/0004-637X/812/1/65.
  * González et al. (2019) González, M.M.; Patricelli, B.; Fraija, N.; García-González, J.A.  Reconcilement of VHE $\gamma$-ray/X-ray correlation studies in Mrk 421 and break-down at high fluxes.  Mon. Not. R. Astron. Soc. 2019, 484, 2944–2953, doi:10.1093/mnras/stz109.
  * Dermer and Schlickeiser (1993) Dermer, C.D.; Schlickeiser, R.  Model for the High-Energy Emission from Blazars.  Astrophys. J. 1993, 416, 458, doi:10.1086/173251.
  * Sikora et al. (1994) Sikora, M.; Begelman, M.C.; Rees, M.J.  Comptonization of Diffuse Ambient Radiation by a Relativistic Jet: The Source of Gamma Rays from Blazars?  Astrophys. J. 1994, 421, 153, doi:10.1086/173633.
  * Kaspi et al. (2000) Kaspi, S.; Smith, P.S.; Netzer, H.; Maoz, D.; Jannuzi, B.T.; Giveon, U.  Reverberation Measurements for 17 Quasars and the Size-Mass-Luminosity Relations in Active Galactic Nuclei.  Astrophys. J. 2000, 533, 631–649, doi:10.1086/308704.
  * Finke (2016) Finke, J.D.  External Compton Scattering in Blazar Jets and the Location of the Gamma-Ray Emitting Region.  Astrophys. J. 2016, 830, 94, doi:10.3847/0004-637X/830/2/94.
  * Tavecchio and Ghisellini (2012) Tavecchio, F.; Ghisellini, G.  “Flat” broad line region and gamma-ray absorption in blazars.  arXiv 2012, arXiv:1209.2291.
  * Lei and Wang (2014) Lei, M.; Wang, J.  The effect of the broad-line region with geometrical structures on gamma-ray absorption in blazars.  Publ. Astron. Soc. Jpn. 2014, 66, 7, doi:10.1093/pasj/pst008.
  * Sikora et al. (2009) Sikora, M.; Stawarz, Ł.; Moderski, R.; Nalewajko, K.; Madejski, G.M.  Constraining Emission Models of Luminous Blazar Sources.  Astrophys. J. 2009, 704, 38–50, doi:10.1088/0004-637X/704/1/38.
  * Błażejowski et al. (2000) Błażejowski, M.; Sikora, M.; Moderski, R.; Madejski, G.M.  Comptonization of Infrared Radiation from Hot Dust by Relativistic Jets in Quasars.  Astrophys. J. 2000, 545, 107–116, doi:10.1086/317791.
  * Tavecchio et al. (2000) Tavecchio, F.; Maraschi, L.; Sambruna, R.M.; Urry, C.M.  The X-ray Jet of PKS 0637-752: Inverse Compton Radiation from the Cosmic Microwave Background?  Astrophys. J. Lett. 2000, 544, L23–L26, doi:10.1086/317292.
  * Meyer et al. (2017) Meyer, E.T.; Breiding, P.; Georganopoulos, M.; Oteo, I.; Zwaan, M.A.; Laing, R.; Godfrey, L.; Ivison, R.J.  New ALMA and Fermi/LAT Observations of the Large-scale Jet of PKS 0637-752 Strengthen the Case Against the IC/CMB Model.  Astrophys. J. Lett. 2017, 835, L35, doi:10.3847/2041-8213/835/2/L35.
  * Schwartz (2002) Schwartz, D.A.  X-ray Jets as Cosmic Beacons.  Astrophys. J. Lett. 2002, 569, L23–L26, doi:10.1086/340482.
  * Ghisellini et al. (2005) Ghisellini, G.; Tavecchio, F.; Chiaberge, M.  Structured jets in TeV BL Lac objects and radiogalaxies. Implications for the observed properties.  Astron. Astrophys. 2005, 432, 401–410, doi:10.1051/0004-6361:20041404.
  * Hervet et al. (2015) Hervet, O.; Boisson, C.; Sol, H.  Linking radio and gamma-ray emission in Ap Librae.  Astron. Astrophys. 2015, 578, A69, doi:10.1051/0004-6361/201425330.
  * MacDonald et al. (2015) MacDonald, N.R.; Marscher, A.P.; Jorstad, S.G.; Joshi, M.  Through the Ring of Fire: Gamma-Ray Variability in Blazars by a Moving Plasmoid Passing a Local Source of Seed Photons.  Astrophys. J. 2015, 804, 111, doi:10.1088/0004-637X/804/2/111.
  * Kino et al. (2002) Kino, M.; Takahara, F.; Kusunose, M.  Energetics of TeV Blazars and Physical Constraints on Their Emission Regions.  Astrophys. J. 2002, 564, 97–107, doi:10.1086/323363.
  * Tavecchio and Ghisellini (2016) Tavecchio, F.; Ghisellini, G.  On the magnetization of BL Lac jets.  Mon. Not. R. Astron. Soc. 2016, 456, 2374–2382, doi:10.1093/mnras/stv2790.
  * Dermer et al. (2014) Dermer, C.D.; Cerruti, M.; Lott, B.; Boisson, C.; Zech, A.  Equipartition Gamma-Ray Blazars and the Location of the Gamma-Ray Emission Site in 3C 279.  Astrophys. J. 2014, 782, 82, doi:10.1088/0004-637X/782/2/82.
  * Gould and Schréder (1967) Gould, R.J.; Schréder, G.P.  Pair Production in Photon-Photon Collisions.  Phys. Rev. 1967, 155, 1404–1407, doi:10.1103/PhysRev.155.1404.
  * Coppi and Blandford (1990) Coppi, P.S.; Blandford, R.D.  Reaction rates and energy distributions for elementary processes in relativistic pair plasmas.  Mon. Not. R. Astron. Soc. 1990, 245, 453–507.
  * Salamon and Stecker (1998) Salamon, M.H.; Stecker, F.W.  Absorption of High-Energy Gamma Rays by Interactions with Extragalactic Starlight Photons at High Redshifts and the High-Energy Gamma-Ray Background.  Astrophys. J. 1998, 493, 547–554, doi:10.1086/305134.
  * Ahnen et al. (2016) Ahnen, M.L.; Ansoldi, S.; Antonelli, L.A.; Antoranz, P.; Arcaro, C.; Babic, A.; Banerjee, B.; Bangale, P.; Barres de Almeida, U.; Barrio, J.A.; et al.  Detection of very high energy gamma-ray emission from the gravitationally lensed blazar QSO B0218+357 with the MAGIC telescopes.  Astron. Astrophys. 2016, 595, A98, doi:10.1051/0004-6361/201629461.
  * Biteau and Williams (2015) Biteau, J.; Williams, D.A.  The Extragalactic Background Light, the Hubble Constant, and Anomalies: Conclusions from 20 Years of TeV Gamma-ray Observations.  Astrophys. J. 2015, 812, 60, doi:10.1088/0004-637X/812/1/60.
  * Broderick et al. (2012) Broderick, A.E.; Chang, P.; Pfrommer, C.  The Cosmological Impact of Luminous TeV Blazars. I. Implications of Plasma Instabilities for the Intergalactic Magnetic Field and Extragalactic Gamma-Ray Background.  Astrophys. J. 2012, 752, 22, doi:10.1088/0004-637X/752/1/22.
  * Plaga (1995) Plaga, R.  Detecting intergalactic magnetic fields using time delays in pulses of $\gamma$-rays.  Nature 1995, 374, 430–432, doi:10.1038/374430a0.
  * Aharonian et al. (1994) Aharonian, F.A.; Coppi, P.S.; Voelk, H.J.  Very High Energy Gamma Rays from Active Galactic Nuclei: Cascading on the Cosmic Background Radiation Fields and the Formation of Pair Halos.  Astrophys. J. Lett. 1994, 423, L5, doi:10.1086/187222.
  * H. E. S. S. Collaboration et al. (2014) Abramowski, A.; Aharonian, F.; Ait Benkhali, F.; Akhperjanian, A.G.; Angüner, E.; Anton, G.; Backes, M.; Balenderan, S.; Balzer, A.; Barnacka, A.; et al.  Search for extended $\gamma$-ray emission around AGN with H.E.S.S. and Fermi-LAT.  Astron. Astrophys. 2014, 562, A145, doi:10.1051/0004-6361/201322510.
  * Archambault et al. (2017) Archambault, S.; Archer, A.; Benbow, W.; Buchovecky, M.; Bugaev, V.; Cerruti, M.; Connolly, M.P.; Cui, W.; Falcone, A.; Fernández Alonso, M.; et al.  Search for Magnetically Broadened Cascade Emission from Blazars with VERITAS.  Astrophys. J. 2017, 835, 288, doi:10.3847/1538-4357/835/2/288.
  * Ackermann et al. (2018) Ackermann, M.; Ajello, M.; Baldini, L.; Ballet, J.; Barbiellini, G.; Bastieri, D.; Bellazzini, R.; Bissaldi, E.; Blandford, R.D.; Bloom, E.D.; et al.  The Search for Spatial Extension in High-latitude Sources Detected by the Fermi Large Area Telescope.  Astrophys. J. Suppl. Ser. 2018, 237, 32, doi:10.3847/1538-4365/aacdf7.
  * Neronov and Vovk (2010) Neronov, A.; Vovk, I.  Evidence for Strong Extragalactic Magnetic Fields from Fermi Observations of TeV Blazars.  Science 2010, 328, 73, doi:10.1126/science.1184192.
  * Taylor et al. (2011) Taylor, A.M.; Vovk, I.; Neronov, A.  Extragalactic magnetic fields constraints from simultaneous GeV-TeV observations of blazars.  Astron. Astrophys. 2011, 529, A144, doi:10.1051/0004-6361/201116441.
  * Dondi and Ghisellini (1995) Dondi, L.; Ghisellini, G.  Gamma-ray-loud blazars and beaming.  Mon. Not. R. Astron. Soc. 1995, 273, 583–595, doi:10.1093/mnras/273.3.583.
  * Ghisellini and Madau (1996) Ghisellini, G.; Madau, P.  On the origin of the gamma-ray emission in blazars.  Mon. Not. R. Astron. Soc. 1996, 280, 67–76, doi:10.1093/mnras/280.1.67.
  * Donea and Protheroe (2003) Donea, A.C.; Protheroe, R.J.  Radiation fields of disk, BLR and torus in quasars and blazars: implications for /$\gamma$-ray absorption.  Astropart. Phys. 2003, 18, 377–393, doi:10.1016/S0927-6505(02)00155-X.
  * Liu and Bai (2006) Liu, H.T.; Bai, J.M.  Absorption of 10-200 GeV Gamma Rays by Radiation from Broad-Line Regions in Blazars.  Astrophys. J. 2006, 653, 1089–1097, doi:10.1086/509097.
  * Ghisellini and Tavecchio (2009) Ghisellini, G.; Tavecchio, F.  Canonical high-power blazars.  Mon. Not. R. Astron. Soc. 2009, 397, 985–1002, doi:10.1111/j.1365-2966.2009.15007.x.
  * Aleksić et al. (2011) Aleksić, J.; Antonelli, L.A.; Antoranz, P.; Backes, M.; Barrio, J.A.; Bastieri, D.; Becerra González, J.; Bednarek, W.; Berdyugin, A.; Berger, K.; et al.  MAGIC Discovery of Very High Energy Emission from the FSRQ PKS 1222+21.  Astrophys. J. Lett. 2011, 730, L8, doi:10.1088/2041-8205/730/1/L8.
  * Nalewajko et al. (2014) Nalewajko, K.; Begelman, M.C.; Sikora, M.  Constraining the Location of Gamma-Ray Flares in Luminous Blazars.  Astrophys. J. 2014, 789, 161, doi:10.1088/0004-637X/789/2/161.
  * Ahnen et al. (2015) Ahnen, M.L.; Ansoldi, S.; Antonelli, L.A.; Antoranz, P.; Babic, A.; Banerjee, B.; Bangale, P.; Barres de Almeida, U.; Barrio, J.A.; Bednarek, W.; et al.  Very High Energy $\gamma$-Rays from the Universe’s Middle Age: Detection of the z = 0.940 Blazar PKS 1441+25 with MAGIC.  Astrophys. J. Lett. 2015, 815, L23, doi:10.1088/2041-8205/815/2/L23.
  * Abeysekara et al. (2015) Abeysekara, A.U.; Archambault, S.; Archer, A.; Aune, T.; Barnacka, A.; Benbow, W.; Bird, R.; Biteau, J.; Buckley, J.H.; Bugaev, V.; et al.  Gamma-Rays from the Quasar PKS 1441+25: Story of an Escape.  Astrophys. J. Lett. 2015, 815, L22, doi:10.1088/2041-8205/815/2/L22.
  * Costamante et al. (2018) Costamante, L.; Cutini, S.; Tosti, G.; Antolini, E.; Tramacere, A.  On the origin of gamma-rays in Fermi blazars: Beyondthe broad-line region.  Mon. Not. R. Astron. Soc. 2018, 477, 4749–4767, doi:10.1093/mnras/sty887.
  * H. E. S. S. Collaboration et al. (2019) Abdalla, H.; Adam, R.; Aharonian, F.; Ait Benkhali, F.; Angüner, E.O.; Arakawa, M.; Arcaro, C.; Armand, C.; Ashkar, H.; Backes, M.; et al.  Constraints on the emission region of <ASTROBJ>3C 279</ASTROBJ> during strong flares in 2014 and 2015 through VHE $\gamma$-ray observations with H.E.S.S.  Astron. Astrophys. 2019, 627, A159, doi:10.1051/0004-6361/201935704.
  * H. E. S. S. Collaboration et al. (2020) Abdalla, H.; Adam, R.; Aharonian, F.; Ait Benkhali, F.; Angüner, E.O.; Arakawa, M.; Arcaro, C.; Armand, C.; Ashkar, H.; Backes, M.; et al.  H.E.S.S. detection of very high-energy $\gamma$-ray emission from the quasar PKS 0736+017.  Astron. Astrophys. 2020, 633, A162, doi:10.1051/0004-6361/201935906.
  * Galanti et al. (2020) Galanti, G.; Landoni, M.; Tavecchio, F.; Covino, S.  Probing the absorption of gamma-rays by IR radiation from the dusty torus in FSRQs with the Cherenkov telescope array.  Mon. Not. R. Astron. Soc. 2020, 495, 3463–3473, doi:10.1093/mnras/staa1391.
  * Boettcher and Schlickeiser (1997) Boettcher, M.; Schlickeiser, R.  The pair production spectrum from photon-photon annihilation.  Astron. Astrophys. 1997, 325, 866–870.
  * Mannheim (1993) Mannheim, K.  The proton blazar.  Astron. Astrophys. 1993, 269, 67–76.
  * Aharonian (2000) Aharonian, F.A.  TeV gamma rays from BL Lac objects due to synchrotron radiation of extremely high energy protons.  New Astron. 2000, 5, 377–395, doi:10.1016/S1384-1076(00)00039-7.
  * Mücke and Protheroe (2001) Mücke, A.; Protheroe, R.J.  A proton synchrotron blazar model for flaring in Markarian 501.  Astropart. Phys. 2001, 15, 121–136, doi:10.1016/S0927-6505(00)00141-9.
  * Begelman et al. (1990) Begelman, M.C.; Rudak, B.; Sikora, M.  Consequences of Relativistic Proton Injection in Active Galactic Nuclei.  Astrophys. J. 1990, 362, 38, doi:10.1086/169241.
  * Atoyan and Dermer (2003) Atoyan, A.M.; Dermer, C.D.  Neutral Beams from Blazar Jets.  Astrophys. J. 2003, 586, 79–96, doi:10.1086/346261.
  * Dermer et al. (2012) Dermer, C.D.; Murase, K.; Takami, H.  Variable Gamma-Ray Emission Induced by Ultra-high Energy Neutral Beams: Application to 4C +21.35.  Astrophys. J. 2012, 755, 147, doi:10.1088/0004-637X/755/2/147.
  * Zhang et al. (2020) Zhang, B.T.; Petropoulou, M.; Murase, K.; Oikonomou, F.  A Neutral Beam Model for High-energy Neutrino Emission from the Blazar TXS 0506 + 056.  Astrophys. J. 2020, 889, 118, doi:10.3847/1538-4357/ab659a.