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

HAL
open science

# Blazars at Very High Energies: Emission Modelling

Hélène Sol, Andreas Zech

► To cite this version:

Hélène Sol, Andreas Zech. Blazars at Very High Energies: Emission Modelling. Galaxies, 2022, 10 (6), pp.105. 10.3390/galaxies10060105. hal-03843642

HAL Id: hal-03843642

https://hal.science/hal-03843642

Submitted on 15 Dec 2022

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

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
Review

# Blazars at Very High Energies: Emission Modelling

Hélène Sol * and Andreas Zech * |

Laboratoire Univers et Théories, Observatoire de Paris, Université PSL, CNRS, Université Paris Cité,
5 Place Jules Janssen, F-92190 Meudon, France

* Correspondence: helene.sol@observatoiredeparis.psl.eu (H.S.); andreas.zech@observatoiredeparis.psl.eu (A.Z.)

**Abstract:** Blazars are very broadband cosmic sources with spectra spanning over twenty orders of magnitude in frequency, down to the 100 MHz regime in the radio range, up to VHE at several tens of TeV. The modelling of their spectral energy distribution at high energies currently considers two main classes of models, leptonic and lepto-hadronic, which both succeed fairly well in describing the observed spectra for the two populations of blazars, namely BL Lac objects (BL Lacs) and flat spectrum radio quasars (FSRQs). However they are both confronted with difficulties, in particular to reproduce flaring phenomena monitored with a good multi-spectral and temporal coverage, or to reproduce extreme sources which challenge the basic descriptions. Such a situation has led to a diversity of specific scenarios, the positioning of which in relation to the general context of the sources is generally not clearly fixed. The identification of the dominant particle acceleration mechanism at work and a better understanding of the location of the TeV emitting zone would make it possible to break the degeneracies between models. Multi-wavelength and multi-messenger studies should also help in this regard, with the perspective to elaborate a general reference scenario of blazars and AGNs.

**Keywords:** blazars; very high energy astrophysics; non-thermal emission

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

check for
updates

**Citation:** Sol, H.; Zech, A. Blazars at Very High Energies: Emission Modelling. *Galaxies* **2022**, *10*, 105. https://doi.org/10.3390/galaxies10060105

Academic Editor: Junhui Fan

Received: 01 August 2022

Accepted: 25 October 2022

Published: 3 November 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

The broadband spectral energy distribution (SED) from blazars can stretch over the full observationally accessible spectrum, from the radio-band up to very-high-energy (VHE) gamma-rays. This continuum emission is dominated by non-thermal radiation in the form of two spectral bumps, with the first peaking in the optical, ultraviolet or X-ray range and the second in the high-energy or VHE gamma-ray band. Depending on the source, non-thermal emission from an extended radio jet and thermal emission from the host galaxy, the accretion disk and the dusty torus can lead to a more complex SED.

The comparison of multi-wavelength (MWL) emission from blazars with numerical models that calculate radiative emission and transfer for given assumptions on the particle content and characteristics of the emission regions, represents the principal tool used to probe the microphysics inside the source and its physical parameters. Information from the VHE part of the SED is a necessity to constrain model parameters for high-frequency peaked (BL Lac) objects, which radiate a significant part of their overall emission in that energy range, but observational VHE constraints can also be important for softer sources, since they probe the very limit of particle acceleration and are also most sensitive to photon-photon absorption on internal or external photon fields.

VHE observations are also crucial for the characterisation of the variability time scales involved in variable emission. The short cooling time at these energies leads to the most rapid variability, imposing rigorous constraints on the size of the emission region, due to light travel time arguments.

We first provide a short review of the main emission scenarios, their applications to different blazar types and their limitations in Section 2. Then we turn our attention to the underlying particle acceleration mechanisms in Section 3. A critical question, which might
not have a unique answer, concerns the location of the emission zone(s) responsible for the VHE emission (Section 4). We finally discuss how MWL and multi-messenger (MM) observations constrain the large variety of proposed emission models in Section 5 and conclude in Section 6 with a brief outlook on the spin-offs of a better understanding of blazar physics.

## 2. Blazar Emission Modelling

While spectra observed in a limited wavelength range can be directly fitted with power-law and log-parabola functions or simple physical models, the interpretation of the full SED is usually based on a comparison with the synthetic MWL emission from radiative models of different degrees of complexity. Such modelling of the SED can exploit proper fitting in the case of simple models (e.g., [1]), but for models with a high number of free parameters and for the often sparse and not strictly simultaneous MWL data sets, an adjustment “by eye” is still the standard approach. Given the limited observational coverage and the intrinsic degeneracy of even the simplest emission models, the objective of SED modelling is rather to identify a range of physically acceptable models that can explain the data, than to find the best formal fit model, which might spuriously lead to a set of extreme parameter values.

### 2.1. Leptonic Emission Modelling

While the low-energy bump in the non-thermal emission is generally ascribed to synchrotron emission from electrons or pairs inside the jet, the high-energy bump can be interpreted with leptonic or hadronic processes. In the standard leptonic approach, electrons or pairs up-scatter the synchrotron radiation they have produced and photons from external fields to (very) high energies through the Inverse Compton process (e.g., [2–9]).

If external photon fields can be neglected, the stationary single-zone synchrotron-self-Compton (SSC) model provides the most basic approach to describe the steady MWL emission (e.g., [10,11]). For a phenomenological approach, the emission region is usually described as a spherical “blob” filled with a tangled magnetic field of homogeneous strength and a single particle distribution parameterised by a broken power law or a logparabola, which generally provide the necessary free parameters to represent the observed shape of the synchrotron and inverse-Compton bumps. Alternatively, a particle spectrum is injected and evolved, based on a differential (Fokker-Planck) equation, which can account for particle injection, acceleration, radiative and adiabatic cooling, and particle escape, into a steady-state particle spectrum. A general form is given by [12,13] as:

$$ \frac{\partial N_e(\gamma, t)}{\partial t} = \frac{\partial}{\partial \gamma} \left[ (b_c \gamma^2 - a\gamma - 2D_0\gamma) \cdot N_e(\gamma, t) \right] + \frac{\partial}{\partial \gamma} \left( D_0 \gamma^2 \frac{\partial N_e(\gamma, t)}{\partial \gamma} \right) - \frac{N_e(\gamma, t)}{t_{\text{esc}}} + Q_{\text{inj}}(\gamma, t) \quad (1) $$

Here $b_c\gamma^2$ is the radiative cooling rate, $a\gamma$ the Fermi-I acceleration rate, $D_0\gamma^2$ the energy diffusion coefficient, $t_{\text{esc}}$ the particle escape time and $Q_{\text{inj}}$ the injection function. This particular description is appropriate in the “hard-sphere” scattering approximation.

The whole emission region is moving towards the observer at relativistic speed, with bulk Lorentz factor $\Gamma$ under a small viewing angle $\theta$, leading to a Doppler boosting of the emission that is characterized by the bulk Doppler factor $\delta = 1/\Gamma(1 - \beta \cos \theta)$ Adiabatic expansion of the “blob” during the considered duration of its emission is mostly neglected in such models. Under certain conditions, the emission region can also be interpreted as a standing shock through which plasma is moving at relativistic bulk speed. The single-zone SSC scenario requires only 8 or 9 free parameters and, even with this simplistic approach, usually succeeds in reproducing the available data for most blazars detected at VHE, but does not provide a unique solution for a given MWL SED. It is however possible
to carry out a scan of the available parameter space and to identify the acceptable range of solutions for a given confidence interval [11,14].

Figure 1 shows the interpretation of an SED of a high-frequency peaked BL Lac object, where the synchrotron and SSC components of a compact emission region account for the X-ray and VHE emission, while additional synchrotron emission from the extended jet dominates the lower frequencies. Such an extended component is generally required to correctly describe at least the low-frequency radio emission. In this case, emission from the host galaxy (thin brown line) is completely dominated by the non-thermal emission. Absorption by the extragalactic background light (EBL) is noticeable at VHE although it remains weak at small redshifts.

<figure><img src="image_5.png" /><figcaption>**Figure 1.** A typical stationary SSC model for the high-frequency peaked BL Lac object PKS 2155−304. See text for details. Credit: [15], A&A, reproduced with permission @ESO.</figcaption></figure>

Leptonic models that account for interactions with external photon fields, from the accretion disk, broad line region, dusty torus, surrounding starlight (or photons from the EBL and the CMB in the case of extended jets or intergalactic particle cascades), become rapidly more complex, even when using a single-zone approach. The large number of additional free parameters to describe these photon fields can be reduced either through direct observational constraints or by using scaling relations that link the different contributions to the accretion disk luminosity or supposed mass accretion rate and radiative efficiency of the central engine [16]. The inclusion of external Inverse Compton (EIC) radiation makes it possible to reproduce the high Compton dominance observed in luminous blazar classes (e.g., [4]).

Figure 2 shows a stationary model including, in addition to synchrotron and SSC emission, EIC radiation from several external photon fields, while the SSC emission dominates at X-ray energies, the high Compton dominance in the gamma-ray range is due to these latter components. Direct emission from the accretion disk and the dust torus is strongly dominated by the non-thermal emission in this particular case.
<figure><img src="image_6.png" /><figcaption>**Figure 2.** A typical stationary model, including a range of EIC components, for the FSRQ 3C 279 (taken from [17] based on a code by [18]). Credit: [17] @AAS. Reproduced with permission.</figcaption></figure>

While such steady-state models can be used to broadly characterize different activity states of blazars, in order to truly model rapidly varying emission, as encountered during blazar flares, a time-dependent description of the electron distribution and/or of the source parameters is required. Flares can then be modelled through an interplay of particle injection or acceleration with particle cooling and escape, following, e.g., Equation (1) as a function of time, or through sudden changes in the magnetic field, Doppler factor or physical extension of the emission region (e.g., [19–21]).

The radiative processes involved in leptonic models, synchrotron and Inverse Compton radiation, as well as photon absorption through pair production, are well known and documented. Beyond the single-zone scenario, a description of the emission from the extended jet is usually achieved by representing the latter as a combination of homogeneous slices or cells and by calculating the radiative transfer along a given line of sight, while this is straightforward for the synchrotron and EIC emission, the calculation of SSC emission from such an extended region is more complicated, since one ideally needs to take into account the radiation emanating from all regions of the jet beyond the SSC emission region under study [22,23].

## 2.2. Lepto-Hadronic Emission Modelling

While leptonic models represent the most economic approach in terms of free parameters and energy requirements, there are several reasons to consider an additional hadronic component that might contribute to the observed broadband emission. One motivation is the open question on the origin of ultra-high-energy cosmic rays (UHECRs) and high-energy neutrinos. Observations with current detectors, such as the Pierre Auger Observatory and IceCube, seem to favour an extragalactic origin at least for a fraction of the observed particle fluxes [24,25], leaving blazars as one of the few potential emitters, together with other types of AGNs, starburst galaxies or GRBs. A direct link between astrophysical sources and detected particles is notoriously difficult to establish, but leptohadronic models offer the possibility of connecting a potential site of acceleration of UHECRs with the expected emission of neutrinos and photons up to the VHE range.
From a more fundamental point of view, there is strong evidence from observations and modelling that relativistic blazar jets contain an important hadronic component [26]. Depending on the assumed particle acceleration process (see Section 3), the existence of a highly relativistic hadron population may thus seem natural.

Even within the single-zone framework, different flavours of such models exist. Given the complexity of nuclear interactions and the presumed strong dominance of hydrogen nuclei, at least for a large range of energies, most models restrict themselves to the treatment of protons [27–29], although models with a more comprehensive treatment of nuclei exist [30]. In the extreme scenario of purely hadronic models [31], all emission is traced back to a relativistic population of hadrons. Leptonic synchrotron emission comes only from secondary electrons created in hadronic interactions. Such scenarios require usually extreme jet powers.

A more common approach is the “lepto-hadronic” framework, where the low-energy emission bump of the blazar SED is ascribed to electron synchrotron emission, as in leptonic models, while relativistic hadrons contribute to—and in some cases dominate—the high-energy emission.

Two distinct regimes in parameter space exist to explain a hadron dominated high-energy bump (e.g., [31,32]), (i) strong magnetic fields, of the order of 1 G to 100 G and large proton Lorentz factors of up to $\gamma_{p,max} \approx 10^{10}$ lead to a dominant proton-synchrotron emission at high energies, or (ii) smaller magnetic fields combined with a large particle and/or photon field density result in an emission dominated by synchrotron-pair or IC-pair cascades. For low magnetic field strengths, mixed solutions can be found where leptonic and hadronic processes contribute both at high energies [33]. Generally, to contribute significantly to high energies, the energy density of relativistic hadrons needs to largely exceed that of relativistic leptons; this could be ascribed to specificities of the acceleration process or to the much faster radiative cooling of electrons.

Different hadronic interaction processes contribute to the emission of very energetic photons and to the production of secondary leptons. Such models need to account for synchrotron emission from hadrons and secondary particles (electrons, muons, and to a lesser extent pions and other mesons), for pair production in the Bethe-Heitler process and, at the highest energies, in photo-pion production processes. A precise description of the latter requires Monte Carlo simulations based on measured cross-sections or descriptions using parameterisations of such simulations. In interactions of ultra-relativistic protons with target photons or protons, the emission region is opaque for first generations of secondary particles, including gamma-rays, leading to successive formations of the above mentioned internal particle cascades. If relativistic nuclei are included, spallation processes need to be taken into consideration as well. A more detailed description of all relevant interactions and emission processes can be found, e.g., in [34] and references therein.

To reduce the additional number of free parameters in lepto-hadronic models that are required to describe the energy distribution of the hadronic population, often the hypothesis of a co-acceleration of leptons and hadrons, resulting in a similar spectral shape for both populations, is invoked. When introducing a certain number of constraints, a parameter scan can be undertaken with lepto-hadronic single-zone models to explore the range of accessible solutions [32,33].

An example for a stationary lepto-hadronic model is shown in Figure 3. In this particular configuration, proton-synchrotron emission dominates in the X-ray band, while the highest energies are reproduced as the sum of synchrotron radiation from cascades that are triggered by proton-photon interactions (Bethe-Heitler and photo-pion processes). The significant contribution from the latter leads also to the emergence of high-energy neutrinos (see Section 5.3).
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Example of a complex lepto-hadronic model, applied to the SED of the BL Lac object 1H 1914-194, including synchrotron and SSC emission from primary leptons, proton-synchrotron emission, and radiation from secondary leptons generated in Bethe-Heitler pair production ("pe") and photo-pion production ("pπ"). In addition, the neutrino spectrum emerging from hadronic interactions is shown at the highest energies. (taken from [35]).</figcaption></figure>

Proton-proton interactions are usually negligible in lepto-hadronic blazar models, given the high particle density and the extreme jet powers that a significant emission from such processes would require [36]. A possible exception is given by models where the jet interacts with an obstacle, which will be discussed shortly below.

## 2.3. Application of Single-Zone Models to Specific Blazar Types

### 2.3.1. Applications to FSRQs and Luminous BL Lac Objects

Steady-state emission from luminous flat-spectrum radio quasars (FSRQs) is generally well described by single-zone models when taking into account external photon fields, while both the broad-line region (BLR) and the dusty torus may provide the necessary photon density, detection of VHE gamma rays from such sources puts the emission region outside or at most at the outer edge of the BLR, otherwise photon-photon absorption would be too significant. Such types of scenarios lead to rather large emitting zones, which have difficulty explaining ultrafast-varying active states (see Section 4.4). As a result, flare emission from the well-known FSRQ 3C 279 proves very challenging for single-zone models, be they leptonic or lepto-hadronic (e.g., [37]).

Intermediate blazar types (LBLs and IBLs, i.e. low-frequency peaked and intermediate-frequency peaked BL Lac objects) may also include a significant EIC component (e.g., [38]). A scenario for the LBL AP Librae, taking into account information from VLBI data, proposes for instance a contribution from the photon field generated by the extended jet surrounding the high-energy emission region [39]. Another model for the same source invokes upscattering of photons from the cosmic microwave background (CMB) by electrons or pairs in a large-scale jet [40]. A lepto-hadronic scenario is also a possible alternative [41].

However a general study opposing leptonic and lepto-hadronic single-zone interpretations of steady-state SEDs from FSRQs, IBLs and LBLs [42], concludes that, unless we modify our current understanding of accretion and jet launching mechanisms, the lepto-hadronic interpretation is disfavoured due to the very large required jet powers, largely exceeding the Eddington limit, that are in contradiction with estimates of jet powers by other methods, requiring at the same time an exceedingly low radiative accretion efficiency.
### 2.3.2. Applications to HBLs and Extreme Blazars

Blazars of the HBL type are generally very well described with single-zone SSC models, as long as the data sets are limited to steady states. These models require usually small magnetic fields, in the range from 0.01 to 0.1 G, far away from equipartition between magnetic and kinetic energy density in the emitting gamma-ray zone. A representation of the synchrotron bump with a single power-law electron distribution, including a spectral break with a change of index $\Delta n = 1$ as expected from a simple radiative cooling, is not possible in most cases. A broken power-law with a much steeper second spectral part or possibly a logparabolic distribution are required. This seems to indicate that the basic homogeneous single-zone picture is a simplification, where a free parametrization of the particle distribution is necessary to model a more complex combination of particle populations and inhomogeneities in the emission region, considering additional effects such as particle injection, acceleration, and escape. Lepto-hadronic models always provide an alternative explanation, with the same caveats as for the luminous blazars. In the case of HBLs, it is possible to find solutions with sub-Eddington luminosities, but the large difference between required jet power and very low disk luminosity still remains puzzling. As for the more luminous blazars, modelling of flares pushes the single-zone approach, whether leptonic or lepto-hadronic, to its limit also for HBLs. The steep flux ratio between the VHE and X-ray band during the gigantic 2006 flares of PKS 2155-304 for example, seems to require at minimum a combination of a region dominating the steady-state emission and a more compact region responsible for the flare emission [15]. To reproduce the MWL light curve and spectral evolution during the 2010 flare of Mrk 421 also requires contributions from several emission or acceleration regions [21], to cite another example among several others.

The validity of the single-zone SSC model is also put into question, when trying to apply it to the most extreme HBLs. Such “extreme-TeV” blazars [43] exhibit very narrow spectral bumps, peaking in the X-ray and TeV range. To reproduce their extremely hard spectra requires a very low magnetic field in the $\mu$ G range, and either very hard particle spectra, beyond what seems acceptable for standard particle acceleration mechanisms, or the artificial introduction of a very high value of the minimum electron Lorentz factor, on the order of $10^3$ or $10^4$ (e.g., [44,45]). Acceleration through magnetic reconnection, while able to explain hard particle spectra, requires relatively high magnetization to be efficient (e.g., [46]), contrary to what is derived from the spectral distributions of these objects. A physical explanation for such unusual parameters for extreme blazars has been proposed in the form of electron-proton co-acceleration on relativistic shocks [47]. Electrons are seen in PIC simulations to be boosted to higher Lorentz factors, as interactions with protons lead to an energy transfer to electrons, until the energies of the two particle species are roughly in equipartition. Low magnetisation is a prerequisite for efficient shock acceleration and arises naturally in this model. While technically a lepto-hadronic scenario, the hadronic emission component is completely negligible, when assuming an equal number density of electrons and protons. To explain the hardest spectra in this scenario requires re-acceleration on a second shock or on turbulences within the jet [48]. Alternative scenarios evoke significant adiabatic losses and specific stochastic acceleration models [49], different types of lepto-hadronic models [33], external cascades in the intergalactic medium [50] or EIC on CMB photons along the jet [51].

## 2.4. Strengths and Weaknesses of Emission Scenarios

As presented above, current emission scenarios can usually reproduce the SED observed during quiescent states of blazars. However, two general problems arise in the lepto-hadronic framework. First, if hadronic processes contribute significantly to the high-energy emission, jet powers are always very large, either dominated by the magnetic energy density in proton-synchrotron dominated scenarios or by the kinetic energy density of the hadrons in scenarios where proton-photon interactions dominate. Sub-Eddington solutions can be found for certain sources and scenarios, but it is generally difficult to reconcile the
large required jet power with the accretion power inferred from known processes and from constraints on the luminosity of the accretion disk. Second, most lepto-hadronic models face difficulties in accounting for rapid variability observed at very high energies, due to the much longer acceleration and cooling time scales for hadrons, as compared to electrons. The characteristic synchrotron cooling time $t_{syn}$ of a particle, for example, depends on its mass $m$ as:

$$t_{syn} = \frac{3m^2c^3}{4\sigma_T U_B \beta^2 E} \qquad (2)$$

with $\sigma_T$ the Thomson cross section, $U_B$ the magnetic energy density, $\beta$ and $E$ the initial velocity (in units of $c$) and energy of the particle. To reach short variability time scales requires scenarios invoking very high magnetic field strengths or possibly jets colliding with matter.

Leptonic emission models can more easily reproduce the short time scales seen in VHE blazar flares, and explain the often observed simultaneous variations of the low-energy and high-energy bumps of the SED. However, the interpretation of very rapid variability down to flux doubling times of only a few minutes, which has been observed in a few exceptional flares [52,53], requires extreme parameters or models that differ greatly from the standard view of an emission region embedded inside the jet and with an extension comparable to the jet radius. Other limitations of current emission models come from certain emission signatures, such as “orphan” high-energy flares without counterparts at lower energies (e.g., [54]) or “extreme” blazars with very hard spectra in the high-energy band (discussed above) (cf. [43] and references therein), which cannot be described with single-zone leptonic models with standard parameters and generally require lepto-hadronic solutions or more sophisticated multi-zone models. It is also interesting to note that, at least in the absence of strong external photon fields, both leptonic and lepto-hadronic single-zone models generally provide solutions that are far out of equipartition between the particle energy densities and the magnetic energy density, which can be debatable in particular for the modeling of the stationary states.

## 2.5. Beyond the Single-Zone Model

The present lack of strong observational constraints on key parameters such as the magnetic field strength and particle densities of the emission regions or their physical extension leads to a degeneracy of even the simplest single-zone descriptions of the steady emission. In the high-energy range, where AGN jets are generally not spatially resolved, only high-quality multi-wavelength data of flaring events can break this degeneracy and, in some cases, reject the simplistic single-zone approach. Conceptually, the single zone is most commonly understood as a plasma “blob” moving along the jet or as an approximation for a stationary region inside the jet with a plasma flowing through it. The emission from particles leaving the single zone is neglected, evoking a strong difference in the physical conditions within and outside of the emission region, such as the magnetic field strength or particle density. An observational justification for the single-zone approach can be found in the compact emission regions seen in VLBI images of moving and standing “radio knots”, which are also seen in the optical and X-ray bands for the most nearby radio galaxies [55,56]. The implicit assumption is that one compact emission region, amplified by Doppler boosting in beamed jets, dominates largely the emission from intermediate to very high energies, while the large-scale emission from the jet is ignored. Simplistic as they may be, single-zone models have been very successful in describing the steady multi-wavelength emission from most BL Lac type blazars and, with the addition of external photons fields, for the more luminous FSRQs. They can also constitute the basic building block for more complex models that may be required to describe emission from flares or extended regions. Describing the physical source with only a few basic parameters—the redshift, radius of the emission region, magnetic field strength and Doppler factor, in most cases—permits to study the impact of the particle distribution and its evolution with time (e.g., [21,57]).
Inhomogeneous jet models propose a more realistic treatment of an extended emission region, which might be identified with the whole of the relativistic jet or a large part of it. The more complete description of the source geometry and its physical parameters is often done at the expense of a less detailed characterisation of the particle distribution. Models beyond the single zone are clearly needed when emission is thought to be extended along the jet and when one wants to describe the different acceleration and cooling time scales of electrons injected into the jet at different energies and propagating along it, leading to non-trivial multi-frequency variability and polarization patterns (e.g., [58–61]).

The simulation of polarization in the radio and optical band in multi-zone or MHD models clearly provides a very useful additional observable that can be verified against observations to probe the structure and degree of order of the magnetic field in different regions of the jet (for a recent comprehensive review see [62]). Temporal changes in the direction of the optical polarization angle have also been identified with gamma-ray flares, providing clues to the structure and location of the emission regions of flares, thus further constraining emission scenarios (e.g., [63,64]). Apart from probing the magnetic field structure, future polarization studies in the X-ray and gamma-ray band, may also show some promise in distinguishing between leptonic and (lepto-)hadronic emission models [65]. This is due to the fact that the high-energy emission in leptonic models is dominated by Inverse Compton emission, leading to a lower degree of polarization than the synchrotron emission expected from proton primaries or from secondary cascade particles in (lepto-)hadronic models.

It should be noted that the treatment of the internal light-travel time effects becomes important when dealing with an extended emission region [66]. This more complex multi-zone approach, when extended to lepto-hadronic modelling, requires in addition the follow-up of several particle species over a large range of simulated cells (e.g., [67]).

A particular example of inhomogeneous scenarios is given by the spine-in-sheath model, where the jet is modelled as a fast moving spine, often consisting of an electron-positron pair plasma, within an often baryonic, less relativistic, sheath layer (e.g., [68–71]). This scenario has been initially proposed to reconcile the high bulk Lorentz factors of a few 10 s inferred from radiative modelling of gamma-loud blazars, with the significantly lower values derived from direct observations of the movement of VLBI radio knots. It can also reconcile very fast variability from a thin spine with a more slowly varying emission from the sheath. Low-energy emission from the surrounding sheath creates also an additional photon field that can add to the EIC component produced by electrons inside the spine. Even if jets from AGNs seem to carry an important fraction of hadrons, the exact lepton to hadron ratio is an open question. To reconcile jet power estimates from radiative modelling with much lower observational estimates from studies of X-ray cavities and the energy input into radio lobes, a contribution from electron-positron pairs seems to be required [72]. This is another indication for a structured jet that contains a pair plasma in addition to a baryonic component.

Multi-zone models can also make it possible to consider a separate (purely leptonic) origin for rapid flares, while still allowing for a lepto-hadronic solution for the steady emission. However, the case of very rapid variability during exceptional flares, with flux doubling time scales down to a few minutes, most strikingly observed in the HBL PKS 2155-304 [53] in 2006, still remains a general problem for current models. Among several propositions to explain these observations are “minijet” [46,73] and “needle-in-jet” [74] scenarios, where the existence of multiple, very compact emission regions inside the jet, induced by magneto-centrifugal acceleration of electron beams or magnetic reconnection in a Poynting flux dominated jet, try to explain the inferred small size and very large Doppler factor of the emission region behind the flares. One way to incorporate very rapid variability in lepto-hadronic models of extended jets is by assuming collisions with clouds, possibly from the broad-line region, or stellar envelopes (e.g., [75,76]).

Taking the description of the jet physics one step further, several models include emission processes in magneto-hydrodynamic (MHD) and particle-in-cell (PIC) simulations
of relativistic jets, while such an approach—computationally expensive of course—leads to a large extension of the parameter space, it carries the promise of mutually constraining the macrophysics of the jet evolution inside the ambient medium and the microphysics of the naturally arising acceleration and emission regions inside the jet. Such models try to directly link the observed emission with the assumed acceleration processes inside the jet, such as acceleration on moving and recollimation shocks—often limited to a link with synchrotron emission in the radio domain—e.g., [77–82], shear acceleration (cf. [83] and references therein), or magnetic reconnection (e.g., [84]). Figure 4 illustrates for example a synthetic radio-synchrotron map of a simulated jet, where the appearance of MWL flares is attributed to shock acceleration of electrons or pairs in collisions of moving and stationary shocks, as observed on kpc-scale in the jet of 3C264 [85]. Such a combination of different approaches seems best suited to exploit at the same time the kinematic and dynamic information on the jet one gains from time-resolved VLBI observations and the information on acceleration and emission processes that are most strongly constrained by the highly variable gamma-ray emission.

<figure><img src="image_8.png" /><figcaption>**Figure 4.** (Left): Synchrotron emission map at $10^{10}$ Hz of an MHD simulation showing the impact of a fast ejecta propagating through a series of recollimation shocks in a relativistic jet, seen under 90°. The x- and y-axes are given in units of the jet-radius $R_{jet}$. (Right): Light curve of the emission in the X-ray range showing the contribution from the ejecta (blue) and the remaining jet (red). The units on the X-axis are years in the observer frame (taken from [86]). Credit: [86], A&A, reproduced with permission @ESO.</figcaption></figure>

**3. Brief Overview of Particle Acceleration Mechanisms**

All current emission models require the presence of a population of charged particles accelerated up to very high energies, which are the true origin of the VHE radiation. However the dominant acceleration mechanisms at work to generate such populations in VHE blazars are not yet well identified. Several possibilities are open as illustrated in Section 2. They have been extensively studied in the literature and are still the subject of in-depth and constantly improving analyzes due to the exponential growth in the power of supercomputers. We only recall here the main ones, supposed to play a major role in AGN jets.

**3.1. Fermi Processes and Shear Acceleration**

Turbulence and shocks are naturally expected to develop in the collisionless plasma of astrophysical jets, providing a wide range of possibilities for activating Fermi-type acceleration mechanisms where particles can reach very high energies through their successive interactions with moving magnetized scattering centers.
The growth of Kelvin-Helmholtz instabilities for instance at the interface of the jet with its surrounding ambient medium, or at internal interfaces between various jet layers or knots with different bulk velocities in stratified jets, can generate turbulence and induce type II Fermi acceleration. Relativistic particles gain energy with an average gain of $<\Delta E> \propto (u_c/c)^2 E$ per interaction, where E is the energy of the particle and $u_c$ the velocity of the scattering center. The characteristic acceleration time $\tau_{acc}$ is independent of E and a population of particles with a differential distribution $dN(E) \propto E^{-p}dE$ can be generated with a power law index $p = 1 + (\tau_{acc}/\tau_{esc})$ if $\tau_{esc}$, the characteristic escape time from the acceleration zone, is also energy-independent. Such type of stochastic acceleration has been shown to generate very high energy particles, both in the quasilinear theory and in the regime of strong nonlinear turbulence [87–93], and is expected to play a role in quiescent and flaring blazars [13, 21, 94].

The presence of shearing flows inside the jet and at its periphery can even enhance the efficiency of such processes by the so-called “shear acceleration”. It was analyzed that instabilities and turbulence induced by shearing in an outer layer of thickness of the order of $0.1R_{jet}$ around a jet of radius $R_{jet}$ can accelerate electrons up to PeV energies and protons up to EeV energies, without destroying the jet by instabilities [95]. The average gain of energy per interaction by a particle can be expressed as $<\Delta E> \propto (u/c)^2 E \propto (l/c)^2 (\partial v_z/\partial x)^2 E$ for a flow with speed $v_z$ along the z-axis sheared in the transverse direction $x$, where $u = l(\partial v_z/\partial x)$ is the effective speed change of the scattering centers as seen by a particle crossing the flow along the x-axis on one mean free path $l$. The acceleration time being inversely proportional to the mean free path, $\tau_{acc} \propto c/l(\partial v_z/\partial x)^2$, shearing acceleration appears faster for particles already accelerated at high energies, and more efficient for protons with large $l$ than for electrons. The distribution of accelerated particles tends to a power law with index depending on the properties of the underlying turbulence [83].

A wide variety of shocks are also expected to develop in supersonic AGN jets, such as reconfinement shocks, stationary or moving internal shocks along the jet, or large scale bow shocks related to hot spots at the extremity of Fanaroff-Riley type II radiosources. All these shocks are privileged sites to launch the activation of type I Fermi acceleration mechanisms [96]. Fast particles can be partially trapped in the vicinity of the shock by their scattering on the magnetic turbulence that they themselves amplify, and cross non-relativistic shocks several times with an important average gain of energy for each round trip $<\Delta E> \propto (\Delta u/c)E \propto [(r-1)/r](u_{shock}/c)$, where $\Delta u$ is the relative velocity between the upstream and downstream medium, $u_{shock}$ the shock speed in the frame of the upstream medium, and $r$ the shock compression ratio, of the order of 4 for strong very supersonic shocks. Here the ratio $\tau_{acc}/\tau_{esc}$ is energy-independent because the same scattering mechanisms are responsible for both particle acceleration and escape, and the differential distribution of particles is then a power law with index $p = (r+2)/(r-1)$ which is of the order of 2 for strong shocks, and presents a rather universal character often observed in the cosmos and obtained in numerical simulations.

The case of ultra-relativistic shocks which can develop in AGN jets with a bulk Lorentz factor $\Gamma_{shock} \gg 1$ is quite special because only energetic particles with individual Lorentz factor $\gamma \gg \Gamma_{shock}$ can successfully cross the shock from downstream to upstream. Moreover in the case of magnetized shocks such crossing from downstream to upstream is possible only for parallel or quasi-parallel shocks, with a small angle between the magnetic field and the flow direction. Type I Fermi processes are then much less efficient, especially in perpendicular and quasi-perpendicular shocks, and typical spectral indices are then softer, of the order of 2.2 to 2.3 as estimated by semi-analytical methods, and Monte Carlo (MC) or Particle in Cell (PIC) simulations [96–103]. However, shock surfing and shock drift acceleration can still develop in quasi-perpendicular shocks because of the non-zero electric field in $V \times B/c$ induced by particle motion towards the magnetized shock. Such mechanisms sometimes described as “fast Fermi processes” allow particles to increase their energy by an order of magnitude in a single encounter with the shock. For weak levels of turbulence, shock drift acceleration can become dominant in oblique shocks and generate
flat-spectrum population of energetic particles up to the highest energies as obtained in MC simulations. Acceleration in relativistic shocks can then produce populations of VHE particles with a large range of power-law indices, from hard index $p \simeq 1$ to very steep ones, depending critically on the properties of turbulence and magnetic field, on the shock speed, and shock obliquity for mildly relativistic shocks [99]. The coupling of such MC simulations with radiation transfer models made it possible to reproduce quite well the typical flaring behavior of the two archetypal blazars 3C 279 and Mrk 501 [20]. One difficulty with the MC approach resides in the powerful test particle approximation and the fact that feedback mechanisms of the accelerated particles on the ambient plasma, shock and turbulence are mostly neglected. This can be questionable for substantial populations of non-thermal particles. Conversely first principles PIC simulations provide a self-consistent description of the microphysics which reduces the efficiency of the Fermi acceleration in relativistic shocks likely due to the growth of instabilities and micro-turbulence which reduces the trapping of particles near the shock. Multiple shock re-acceleration scenarios could then provide a way out to explain the hardest gamma-ray spectra observed in extreme TeV blazars [47]. However, due to its high computational cost, the PIC approach also has its own limitations. In particular it does not yet describe the long term evolution of the shock and large scale phenomena, and has difficulties to reach the high individual particle Lorentz factors needed to reproduce VHE spectra.

Interestingly it has been shown that the different aspects of Fermi-type acceleration processes can be described within a same formalism applicable to sub- or ultra-relativistic cases, turbulent, shear and shock acceleration, in flat or nonflat spacetime [104]. Such types of approach could make it possible to better identify the domains of predominance of the various Fermi acceleration regimes, including also the immediate vicinity of massive black holes.

## 3.2. Magnetic Reconnection

Breaking and reconnecting magnetic field lines can release substantial amount of magnetic energy transmitted to the plasma in various forms, including significant particle acceleration. The time evolution of the magnetic field written in the context of MHD illustrates this type of phenomena, $\partial B/\partial t = \nabla \times (u \times B) - \nabla \times \eta(\nabla \times B)$, where $\eta = c^2/4\pi\sigma$ is the plasma magnetic diffusivity, with $\sigma$ its conductivity. In the highly conducting astrophysical plasmas, the first term on the right hand side generally dominates the magnetic field behavior, the magnetic field is frozen in the plasma and no reconnection occurs. However the diffusive term can become dominant in various situations, for instance when the first term dwindles to zero like around stagnation points, lines or surfaces with $u \simeq 0$. The magnetic field then diffuses in a time $\tau_{diff} \simeq L_B^2/\eta$, $L_B$ being the spatial scale of magnetic field variations, possibly inducing strong currents and current sheets with non-zero electric field, resulting in plasma heating and particle acceleration. Growing instabilities and turbulences induced by the currents can then lead to the development of an anomalous resistivity which further amplifies the dissipation of the magnetic field. The study of faster and even explosive magnetic reconnection events, which can occur in hot and diffuse plasmas, requires describing the microphysics where ions and electrons decouple and new terms appear especially in the Ohm's law [105]. Numerical simulations based on fundamental kinetic Vlasov-Maxwell equations show that relativistic magnetic reconnection can then accelerate particles up to very high energies with a power law distribution with index $p \simeq 1.5$ to 1 for high magnetisation [106,107]. Despite the fact that the various aspects of magnetic reconnection are not yet fully understood, there have been great advances in the field thanks to numerous theoretical works and a large number of studies presenting simulations 2D and 3D PIC simulations of relativistic magnetic reconnection [108–116] (see [106,107] for reviews). Such phenomena should naturally occur in AGN jets where various instabilities can grow, with complex and turbulent velocity fields, eddies and sheared flows, and are also observed in simulations of black hole magnetospheres [117–119]. They attract more and more attention to explain AGN flares and variabilities [46,73,120–123].
### 3.3. Pulsar-like Acceleration around Black Hole

Like pulsar magnetospheres [124], it seems that the magnetosphere of rotating black holes embedded in a magnetic field can very efficiently accelerate particles, either by centrifugal forces present throughout the magnetosphere [125,126], or by strong electric fields which grow in gaps, voids of plasma which can develop around critical surfaces in the vicinity of the black hole [127,128]. Particles are expected to be present in the magnetospheres, possibly generated by pair creation or injected from the turbulent accretion flow and disk. For weak enough ambient radiation field and typical magnetic field of a few tens of Gauss, they can then reach high energies at the outer light surface, estimated to possibly reach 10 to 100 TeV for electrons and positrons, and up to $10^{18}$ to $10^{20}$ eV for protons which do not suffer high synchrotron losses [129–138]. Such effects have been proposed especially to explain ultra-fast varying AGN flares (see Section 4) and should at least contribute to inject a population of energetic particles at the base of jets and winds from the central AGN engine (e.g., [133,137,139]).

A combination of the different acceleration processes mentioned above is possible and even likely, pending a full understanding of the main accelerator. Different acceleration zones and mechanisms could contribute depending on the local plasma properties such as velocity fields, turbulence, magnetization. Magnetic reconnection requires high magnetization and could significantly contribute for instance at the jet base which is magnetically dominated for jets extracted by the Blandford-Znajek process, while efficient shock acceleration which requires lower magnetization could dominate further in the jet, when enough magnetic field has been dissipated. Indeed it is already often considered that pre-acceleration by magnetic reconnection is at work inside jets and provides particles energetic enough to efficiently interact with plasma waves and structures in Fermi-type acceleration zones. Finding the location of the emitting zone should help identify the dominant acceleration mechanisms and the relative importance of leptonic versus hadronic processes.

## 4. Locating the VHE Emitting Zone: A Critical Missing Link

### 4.1. *Black Hole Magnetospheres*

Very fast variability detected in some blazars as PKS2155-304 and Mrk501 and in the radiogalaxies M87 and IC310 at TeV energies suggests from the causality argument that the size of the VHE zone should be comparable to or even less than the size of their black hole horizon, even when assuming rather high jet bulk Lorentz factors. Compact black hole magnetospheres where particle acceleration can be very quick and efficient then appear as a natural site for the origin of highly variable VHE flares, depending on the relative importance of the γ-γ absorption phenomena. In the presence of bright accretion disks usually found in FRSQs, such absorption is expected to drastically limit the propagation of VHE photons. However absorption effects can be overcome in AGNs like BL Lac objects or radiogalaxies with sub-luminous disks in the ADAF regime and sufficiently small magnetic fields.

Advanced instruments such as the Event Horizon Telescope and GRAVITY are now able to detect very fine structures in the core of nearby AGNs, which recently motivated massive numerical simulations of black hole magnetospheres, using PIC codes in general relativity as illustrated in Figure 5. In a quasi force-free approximation, they show the launching of relativistic jets by the Blandford-Znajek process, as well as the formation of a significant current sheet with magnetic reconnection events in the equatorial plane of a black hole in extreme rotation [135]. Considering a plasma of electrons, positrons and photons, and taking into account radiation transfer with inverse-Compton scattering and pair creation, GRPIC simulations can explore the formation of gaps where fast particle acceleration occurs due to the unscreened electric field, together with the generation of high energy radiation in the magnetosphere [136,140]. An intermittent gap develops at the base of the jet near the inner light surface for a certain range of optical depth, and ejects electron-positron plasmons which induce fast-varying gamma-ray emission. Although the amplitude of these self-consistently generated bursts appears at first sight to be slightly
lower than the amplitudes of the gamma-ray flares currently observed from TeV-emitting AGNs, such a phenomenon deserves attention because it could allow us to probe at VHE the immediate vicinity of the supermassive black hole's horizons. Sudden changes in the magnetosphere parameters, such as the local accretion rate, could possibly enhance the observed variabilities [140]. Strong and rapid TeV flares dominated by curvature emission of particles accelerated in a spark gap developing at a few gravitational radius $r_g$ from the horizon could be launched for instance by abrupt changes of the disk emission or of magnetospheric currents [141].

<figure><img src="image_9.png" /><figcaption>**Figure 5.** 3D general relativistic particle-in-cell simulations (GRPIC) of a rotating black hole magnetosphere. Snapshot of the normalized total plasma density (in brown) with some magnetic field lines (in blue) at $t = 40 r_g/c$, showing a persistent fragmented equatorial current sheet where particles are accelerated by magnetic reconnection, very close to the black hole horizon (in black). Taken from [118].</figcaption></figure>

## 4.2. Nuclear VLBI Jets

The base and the inner parts of the jets observed as radio core and knots by VLBI on pc and sub-pc scales are often considered as the most natural place for the VHE emitting zone since it ensures that observed fluxes and variabilities can be directly enhanced by Doppler boosting if the VHE zone is embedded in the relativistic jet bulk outflow. A wide variety of emission models make this assumption over years, either explicit or implied. From our current views on hadronic and leptonic models, one can think that leptonic scenarios with magnetic field in the range of 0.01 to 0.1 G could apply to VHE emitting zones along the VLBI jets while hadronic scenarios which require much higher magnetic field could concern the base and very inner parts of the jets, possibly inside still unresolved radio cores.

Indeed, there is growing evidence for correlation between some HE and VHE flares and various phenomena detected in radio VLBI such as variations in the core, emergence of new knots at the jet base, or appearance of stationary features [142–147], which proves that the emission of at least a part of the flares is related to the nuclear jet structure (see Figure 6). However the causal link has not yet been fully established, and the statistical significance for the whole blazar population is still difficult to assess.
<figure><img src="image_10.png" /></figure>

**Figure 6.** Time evolution of the blazar OJ 287 on parsec-scale, showing the appearance of a new VLBI feature (knot K) in the jet following the first ever VHE flare detected from this source during the period 1–4 February 2017, and the occurrence of a high radio emission state around the TeV activity. (a) Total intensity VLBI images at five epochs from 86 GHz GMVA observation. Red circles represent the model fit components. The new feature K appeared in March 2017, at 0.2 mas from the radio core, between the two quasi-stationary components S1 and S2. (b) Multiwavelength radio light curves from 2015 to 2017. The vertical dashed line shows the time of the VHE flaring state and shaded areas the periods of enhanced activity. The TeV detection was also coincident with increased activity in the X-ray and gamma-ray ranges. Credit: [147], A&A, reproduced with permission @ESO.

Search for specific variability patterns in HE and VHE light curves is a promising way to further explore possible underlying signatures of the VLBI structure, especially in HBLs which often show a series of stationary VLBI knots possibly interpreted as shock diamonds where particle acceleration and enhanced radiation can occur. Some clues suggesting that a perturbation moving at constant speed along the jet induces a series of flares when crossing the successive shocks have been found for instance in the case of Mrk 421, where the main emission zone would be at about 10 pc from the core, with secondary zones further down in the jet [148]. This could be confirmed if a unique variability pattern was detected in X-rays or gamma-rays after each strong flare. Such a temporal coincidence between the occurrence of a bright flare and the interaction of a fast-moving VLBI knot with a stationary jet feature has been also recently observed in the FSRQ PKS1510-089 [149]. If that correlation is true, the VHE fast-varying flare comes from a zone located at about 50 pc from the core, outside the BLR, which is consistent with the lack of clear BLR absorption features in the spectrum but raises the question of the origin of the rapid variability.

### 4.3. Large Scale Jets

Recent gamma-ray observations of the radiogalaxy Cen A that resolved its radio and X-ray jet at VHE definitely prove that large scale AGN jets can radiate in the TeV domain and suggest that in situ acceleration of ultra-relativistic electrons can be at work
over several kpc along AGN jets [150]. Such extended VHE emission is often ascribed to comptonization process of various types of soft photon background [151,152]. In some cases the geometry of the inverse-Compton scattering process can even lead to large angle gamma-ray emission, for instance if an isotropic distribution of relativistic electrons from the decelerated kpc-scale jet interact with low energy photons from the inner relativistic jet. Beyond blazars, unbeamed large scale jets could therefore also produce significant VHE fluxes, reachable by present IACTs for the brightest ones, or by CTA.

Such large scale emitting jets could in principle explain or contribute to the stationary VHE fluxes detected from blazars the gamma-ray emission of which is not angularly resolved, but they are poorly adapted to explain any fast variability. This shows in filigram that at least two types of VHE emission zones are present in AGN, although blazars and generally all ultrafast AGN outbursts are likely dominated by compact Doppler boosted emitting zones.

## 4.4. Constraints from VHE Observations

The γ-γ absorption effect can impose severe constraints on the location of the VHE emitting zone in the presence of strong ambient photons with energy $\epsilon_{bkgd} \gtrsim m_e^2 c^4 / \epsilon_\gamma$ since gamma-rays with energy $\epsilon_\gamma$ can not escape when the optical depth due to pair creation on ambient photons becomes too important. Although possibly negligible for BL Lac objects, this effect should drastically affect FSRQs with bright disks, dusty tori and BLR, and impose some minimal distance of the VHE emitting zone from the central engine and the BLR. This suggests that the emitting gamma-ray zone is located beyond the BLR in FRSQs and appears coherent with studies searching for absorption effects in the gamma-ray spectra of blazars, which detect very few features of absorption by the BLR [153] but find some signatures of absorption by the more extended narrow line region (NLR) in the sub-TeV range [154]. Conversely, it is puzzling to note that correlation was found between the variations of BLR emission lines and gamma-ray flux in some AGN such as PKS1510-089 [155], which suggests on the contrary that part of their HE emission comes from Inverse-Compton scattering with ambient photons from the BLR.

Moreover, as mentioned in Section 2, ultrafast variability appears rather frequent in blazars and radiogalaxies and still remains hard to grasp despite the development of several scenarios (see for instance [156]). A strong tension therefore remains on the determination of the location of the VHE emitting zone especially for fast-varying FSRQ such as 3C 279 or PKS 1222+216 for which the causality argument provides a severe upper limit on the size of the emitting plasma while absorption effect constrains the emission to come from rather large spatial scales, far from the central engine. Two main alternatives can be considered for such sources. Specific geometry of the BLR could allow HE photons to escape from the compact core and jet base if the jet is shielded from external photons by a plasma sheath [157]. Another possibility is to assume that the fast-varying flares come from further down the outflow, from active small sub-structures in a more extended jet [46,74,158]. More exotic alternatives such as the existence of axion-like particles could also solve the puzzle by allowing VHE gamma-ray to escape from the black hole vicinity below the BLR due to a significant reduction in γ-γ absorption.

The new generation of VHE gamma-ray instruments such as the future imaging atmospheric Cherenkov Telescope Array (CTA), together with the VHE wide field detectors HAWC, LHAASO and SWGO, will detect much larger samples of BL Lacs and FRSQs, with better temporal and spectral coverage and resolution. They will search for specific spectral features and especially for additional components in the multi-TeV range which can be expected in the presence of hadronic cascades. They should also make it possible to gather many complete and detailed datasets on blazar gamma-ray flares, which are currently very few. Search for periodicities and quasi-periodicities, analysis of the power spectral density of VHE lightcurves, of possible lognormality and of the nature of the noise for different activity levels, studies of time delays, spectral features and hysteresis for various sources, should better constrain the location, geometry and kinematics of
the emitting zone(s), as well as the main leptonic or hadronic scenarios. For example, the detection or not of time delays in the arrival times of VHE photons as a function of their energy should strongly constrain the origin of flares in single-zone models, often interpreted as due to competing effects between the particle injection and acceleration and the radiative cooling, which generally induce intrinsic spectral lags, although there is still no firm detection confirmed by current instruments [159,160]. The persistent non-detection of such time delays could be an argument in favor of flare scenarios based on a simple change in the bulk Doppler factor of the emitting zone, which have been proposed in some sources for instance in the case of helical jets [161] and in the context of strong recollimation shocks [162]. The exploration of very fast sub-minute time scales accessible thanks to the expected jump in sensitivity promises to be particularly interesting and should drastically constrain the properties of the emitting zones and emission processes [163]. Deep VHE analysis and long-term monitoring of archetypal sources of HBL, IBL, LBL and FSRQ populations, of the extreme UHBL blazars which challenge the current scenarios, and of some peculiar sources such as the binary black hole system OJ 287 [147,164] or the gravitationally lensed QSO B0218+357 [165] should further characterize links between the VHE properties of the AGN and the global characteristics and fine structure of their central core and jet.

**5. Multi-Wavelength and Multi-Messenger Astrophysics**

Beyond analysis of unidentified VHE sources off the galactic plane, which may reveal as yet unknown types of extragalactic cosmic sources, the MWL and MM approach will be essential to characterize the general properties of high energy AGN, to generate a still-missing sample of tracked flares with good temporal and frequency coverage, and to firmly break the degeneracy between the different versions of plausible models of gamma-ray blazars, and in particular between the leptonic and hadronic emission scenarios. The current development of large simultaneous campaigns and monitoring, alert networks, ToO programs and data-driven detection appears mandatory to optimize the scientific return of many experiments and large infrastructures catching blazar data. We just mention a few examples below.

**5.1. The Infrared, Optical, X-ray and Low-Energy Gamma-Ray Domains**

Non-thermal emission is usually expected to dominate in blazars, especially in BL Lacs, and in the gamma-ray domain. However both non-thermal and thermal emission can be detected in the lower energy ranges. Although not always easy to disentangle, the two of them are important to monitor in order to better constrain and describe, respectively, (i) the populations of HE and VHE particles, (ii) the AGN environment with its ambient photon fields (disk, corona, torus, BLR, NLR, clouds, gaseous and stellar components, host galaxy), and to deduce the causal links between the different AGN components.

Several outstanding projects and infrastructures will explore the sky from IR to X-rays in the coming decades, such as the JWST, Vera Rubin Observatory (LSST), ELT, Athena and many others. They will offer plenty of opportunities to conduct deep and wide surveys as well as short and long term MWL monitoring of AGN, making possible the statistical analysis of correlations and time lags between different frequencies. This should characterize the populations of emitting particles and identify the sequence of processes, shedding critical new light on the long-standing question of the origin of blazar variability. High-quality spectra together with polarimetric monitoring will be crucial to explore the nature of the optical and X-ray emission, probe the magnetic fields, and disentangle its various components for the different types of sources [62]. Spectra will also be requisite to complete our knowledge on the still unknown redshifts for many BL Lacs [166]. The Imaging X-ray Polarimetry Explorer (IXPE) launched on 2021 December 9 will observe several AGN and blazars during its first years and already provided spectropolarimetric constraints on the radiogalaxy Cen A, suggesting that the X-rays from its core are due to Compton scattering of low-energy photons by non-thermal electrons accelerated in regions with highly disor-
dered magnetic fields within a few parsecs of the black hole, possibly around the innermost jet [167]. In blazars, X-rays explore either the high-frequency part of the synchrotron bump, or the low-frequency part of the (leptonic or hadronic) high-energy bump, and X-ray polarization will be a strong diagnostic of the different radiation and acceleration models. IXPE and other upcoming X-ray polarization missions, also combined with MWL data, should shed new light on many fundamental questions of blazar physics, such as the magnetic field properties and the nature of the dominant emission processes [168], the relevance of multiple emission zone models in turbulent jets [169,170], and the relative importance of shocks, turbulence and relativistic magnetic reconnection to accelerate particles during HBL flares [171,172].

Indeed IR and X-ray data have already drastically modified our standard view of the so-called "dusty tori" which now appear clumpy and dynamic across various spatial scales, and often present a polar component in addition to the equatorial toroidal structure [173]. IR interferometry of many local AGN has shown two nuclear components on scales of 0.1 to 10 pc, possibly described as an equatorial inflowing disk and a polar-extended feature due to a dusty outflowing wind [174]. If applied to FSRQ, such results would alter the properties of the expected ambient photon field around and between the broad and narrow line regions and would have a significant impact on γ-γ absorption and external inverse-Compton emission at the location of the nuclear jet.

The next generation of gamma-ray space telescopes will be crucial for simultaneously accessing the entire gamma-ray domain. In particular, projects such as Amego and Amego-X, the All-sky Medium Energy Gamma-ray Observatory, should fill the “MeV gap” in sensitivity [175]. Gamma-ray polarimetry in principle has the ability to distinguish between different radiation mechanisms and to probe hadronic signatures in blazars. The detection of a high MeV polarization would be a direct indication that the proton synchrotron dominates the high energy bump of the source

Other MWL missions, such as the Transient High Energy Sky and Early Universe Surveyor (THESEUS) with two high-energy instruments covering X-rays and soft gamma-rays and one IR telescope on board, selected for an assessment study with a possible launch in the 2030s, will be especially important to explore the time domain of AGN and blazars [176]. Due to its large field of view, THESEUS should for instance be able to catch and follow transient phenomena from a large population of AGN, discover new cases of Quasi-Periodic Eruptions (QPE), and study the peculiar “changing state AGN” which evolve from one AGN type to another one over various timescales. It should gather very accurate time lags between various high-energy bands, and cover the full X-ray duty cycle of HBL flares, directly probing particle acceleration and radiation processes.

## 5.2. *Synergy with Radio VLBI and Absolute Sub-Mas Astrometry*

As already illustrated in Section 4.2, the synergy between the very high angular resolution provided by radio VLBI and the very high temporal resolution explored by IACT at VHE appears very promising, already now (see Figure 7) but also with the first operations of the SKA, CTA, and other facilities expected for this decade, together with the concept of Global VLBI Alliance under development, which will change the landscape for the study of the non-thermal universe [177]. VLBI networks will have the power to resolve blazars down to their inner core-jet region. Fast and multi-frequency VLBI imaging in total and polarized light of the vicinity of supermassive black holes and nuclear jets for a sample of gamma-ray emitting AGN will explore the detailed properties and the evolution of the radio-emitting plasma and of the magnetic field, and should clarify the jet-launching, particle acceleration and energy release mechanisms. Coordinated with VHE monitoring, this might be the unique way to firmly and accurately identify the location of the VHE flare emitting zones, and to understand the respective role of the black hole magnetosphere and nuclear jet in the generation of quiescent and fast-varying gamma rays. The cases of binary black hole VHE blazars and of gravitationally lensed ones, although more complex, might in addition offer the opportunity to observe the same relativistic jet from different viewing
angles and configurations, providing a first insight into the three-dimensional properties of
the phenomena.

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

**Figure 7.** Innermost regions of an AGN: standard sketch of a relativistic jet, showing the launching, collimation, acceleration and beginning of propagation zones as a function of the BH gravitational radius $R_g$. A typical resolution is given here for an AGN at low redshift, for four global and space (RadioAstron) VLBI networks at different frequencies, illustrating the possibility to well explore the jet collimation and propagation at the multi-parsec scale, and to probe even into the jet launching region far below the parsec scale, depending on the opacity in the jet and surroundings. Taken from [177].

The International Celestial Reference Frame (ICRF3) and Gaia Data Releases now provide absolute sub-millisecond astrometry at three radio frequencies and in the optical range with a precision of the order of 0.1 mas. Positional differences for more than 3500 AGN cross-identified in ICRF3 and EDR3 reveal some significant offsets between the radio and optical centroids at the mas level. A part of the optical centroids appear coincident with VLBI knots, and a few of them even show high linear polarization [178,179]. Such results suggest that some of the Gaia centroids from blazars are dominated by a non-thermal synchrotron optical emission from an active zone along the nuclear jet, and could in such a case be a direct signature of the location of the gamma-ray emitting zone [180]. This is particularly promising for BL Lacs where other contributions to the optical emission, such as radiation from the accretion disk, are expected to be very low. Progress in this area will require coordination between VHE monitoring and sub-mas astrometry surveys in the radio and optical domains.

5.3. *The Link with Neutrinos and UHECRs*

Due to their deflection in intergalactic and Galactic magnetic fields, it is notoriously
difficult to establish the sources of ultra-high-energy cosmic rays (UHECRs) as presently
observed with cosmic ray detectors such as the Telescope Array and the Pierre AUGER
Observatory. Spatial correlation between the arrival directions of UHECR and various
populations of VHE cosmic sources such as AGN have been analyzed, but no significant
one was found up to now [181]. More specific answers are expected soon with the upgraded
detectors, and later with new-generation experiments such as POEMMA, GRAND and
GCOS. It is also interesting to note that, since AGN jets accelerate electrons, protons and
ions, very long term jet variability could be related to features of the UHECR spectrum [182].
For the time being, only the detection by IceCube, or by new projects like the KM3NeT in
the years to come, of high-energy neutrinos from astrophysical objects permits to directly
infer the acceleration of (U)HECRs in well identified cosmic sources. Lepto-hadronic
emission models then provide a direct link between a population of cosmic rays and the
expected MWL and neutrino emission for a given scenario.
Given the low detection rate of astrophysical neutrinos and the limited angular resolution of their arrival directions, temporal correlations with variable emitters help to significantly improve the identification of their sources. The recent first detection of a high-energy neutrino from the blazar TXS 0506+056 in a flaring state [25] provided thus an important opportunity to compare different emission scenarios to the MWL and multi-messenger data. A large number of lepto-hadronic scenarios have been confronted to this data set¹; a short review is given in [183]. The scenarios that most convincingly describe the observed MWL and neutrino emission during the flaring state imply a combination of an SSC model and a hadronic component that requires in addition external photon fields, from the surrounding jet or emission regions in the AGN, to reach the estimated neutrino flux. The subsequent detection of a neutrino flare from the same source, which was not accompanied by any significant flux increase in the MWL emission [184], presents clearly a challenge for single-zone lepto-hadronic models, and seems to indicate that neutrinos and gamma-rays may well be produced in separate emission regions.

Statistical studies suggest as well that jetted AGN with strong Doppler boosting could be important sources of neutrinos, as shown for instance by the fact that the sub-group of AGN associated with the highest energy neutrinos ($E > 200$ TeV) appears more core-dominated than the whole complete VLBI-flux-density limited sample analyzed in [185]. It has also been pointed out that in some blazars that may coincide with neutrinos recorded by IceCube, the gamma-ray emission is temporarily suppressed during efficient neutrino production due to increased $\gamma$-$\gamma$ opacity [186]. If its existence were confirmed, such a particular class of blazars could contribute to the isotropic cosmic flux of HE neutrinos discovered by IceCube in 2013 [187], the origin of which is still unknown.

## 5.4. Perspective with Gravitational Waves

Our knowledge of the gravitational wave (GW) universe will improve dramatically over the next few decades. Several events directly related to the massive BH of blazars and AGN, such as merger, capture of compact stars, spiralling phase, should be explored by Pulsar Timing Arrays and future gravitational wave space interferometers like LISA, whose launch is expected in 2034. As already observed in GRB and in TDE, the spiralling phases can lead to the formation of relativistic jets from which an emission of gamma-rays is expected and could be detected for “not too distant” transient jets directed towards the solar system. CTA’s wide field of view would be well suited to search for such gamma counterparts of GW events.

Currently, the origin and formation sites of the population of binary black hole (BBH) mergers detected by LIGO-Virgo remains uncertain as electromagnetic counterparts are very weak and difficult to detect in the absence of a dense gaseous environment. However there are theoretical arguments in favor of a possible AGN-disk origin for some BBH events. It could be the case of GW190521 for which an optical AGN flare was suggested as a possible electromagnetic counterpart (although debated) [188]. Standard AGN accretion disks are conducive to stellar mass BBH mergers and can lead to hierarchical formation of black holes, providing an explanation for the large BH masses currently observed. Such environment also induce an excess of eccentric mergers like GW190521 [189], and could produce significant electromagnetic counterpart due to the dense ambient plasma surrounding the merger. Indeed, non-zero eccentricity and non-zero spin-orbit tilt directly found in gravitational data could characterize the mergers occurring in AGN disks, in contrast to isolated mergers or mergers in stellar clusters. This might be difficult to achieve with LIGO-Virgo but will be possible with the Einstein telescope project. Prompt optical-UV counterpart due to the formation of shocks has been predicted from such stellar mass mergers, when the sphere of influence of the binary system is larger than the disk thickness, or depending on the disk opacity [190]. Transient broad band emission from radio to high energies, as well as HE neutrinos and cosmic rays, can also be expected from 10 to 100 solar mass mergers in AGN disks during episodes where the BBH accrete sufficient gas at high rate, exceeding the Eddington limit and generating fast outflows and transient relativistic
jets through the formation of mini-disks [191,192]. Gamma-rays could be detectable with present and future instruments in the case of energetic outflows, even before the merging event. Such MWL and MM events would correspond to a very specific type of blazar flares which have not yet been much considered, and would provide completely new information on AGN formation and accretion disks properties.

**6. Conclusions and Outlook**

Beyond blazar physics, blazar data and modelling at VHE are also mandatory for several questions of fundamental physics and cosmology, gamma-ray emitting blazars being commonly used as beacons to explore their line of sight and probe the nature of space-time through the analysis of the propagation of their VHE photons on cosmic scales. Any progress in the understanding and description of the intrinsic gamma-ray spectra and variability of blazars will significantly impact all results on EBL, intergalactic magnetic fields (IGMF), Lorentz invariance violation (LIV) studies, and axion-like particles searches. A better knowledge on blazar behavior would also usefully contribute improving the accuracy of celestial reference frames, which use a majority of blazars as cosmic reference points and are essential for many activities such as geophysics, GPS, space navigation. Blazar research clearly has a multidisciplinary interest. As illustrated in this review, fundamental difficulty in consolidating any emission modeling is the current absence of a general baseline scenario of blazars and active galactic nuclei sufficiently developed to consistently include the new information collected in recent years. This should be overcome in the coming decade where a wealth of high quality data is expected from several new MWL and MM infrastructures, with the advances in computing power and theoretical models.

**Author Contributions:** Conceptualization, H.S. and A.Z.; methodology, H.S. and A.Z.; software, H.S. and A.Z.; validation, H.S. and A.Z.; formal analysis, H.S. and A.Z.; investigation, H.S. and A.Z.; resources, H.S. and A.Z.; data curation, H.S. and A.Z.; writing—original draft preparation, H.S. and A.Z.; writing—review and editing, H.S. and A.Z.; visualization, H.S. and A.Z.; supervision, H.S. and A.Z.; project administration, H.S. and A.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

**Note**

¹ cf. http://tevcat.uchicago.edu/?mode=1&showsrc=309 (accessed on 31 July 2022).

**References**

1. Nigro, C.; Sitarek, J.; Gliwny, P.; Sanchez, D.; Tramacere, A.; Craig, M. agnpy: An open-source python package modelling the radiative processes of jetted active galactic nuclei. *Astron. Astrophys.* **2022**, *660*, 18. [CrossRef]

2. Konigl, A. Relativistic jets as X-ray and gamma-ray sources. *Astrophys. J.* **1981**, *243*, 700. [CrossRef]

3. Dermer, C.D.; Schlickeiser, R. Model for the High-Energy Emission from Blazars. *Astrophys. J.* **1993**, *416*, 458. [CrossRef]

4. 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. [CrossRef]

5. Levinson, A.; Blandford, R. Pair cascades in extragalactic jets. II. The beamed X-ray spectrum. *Astrophys. J.* **1995**, *449*, 86. [CrossRef]

6. Inoue, S.; Takahara, F. Electron Acceleration and Gamma-Ray Emission from Blazars. *Astrophys. J.* **1996**, *463*, 555. [CrossRef]

7. Ghisellini, G.; Madau, P. On the origin of the gamma-ray emission in blazars. *Mon. Not. R. Astron. Soc.* **1996**, *280*, 67. [CrossRef]

8. Ghisellini, G.; Maraschi, L. Bulk acceleration in relativistic jets and the spectral properties of blazars. *Astrophys. J.* **1989**, *340*, 181. [CrossRef]

9. Katarzyński, K.; Sol, H.; Kus, A. The multifrequency emission of Mrk 501, from radio to TeV gamma-rays. *Astron. Astrophys.* **2001**, *367*, 809–825. [CrossRef]

10. Bednarek, W.; Protheroe, R.J. Testing the homogeneous synchrotron self-Compton model for gamma-ray production in Mrk421. *Mon. Not. R. Astron. Soc.* **1997**, *292*, 646. [CrossRef]

11. Tavecchio, F.; Maraschi, L.; Ghisellini, G. Constraints on the Physical Parameters of TeV Blazars. *Astrophys. J.* **1998**, *509*, 608. [CrossRef]
12. Kardashev, N.S. Nonstationarity of Spectra of Young Sources of Nonthermal Radio Emission. Sov. Astron. **1962**, *6*, 317.

13. Tramacere, A.; Massaro, E.; Taylor, A.M. Stochastic Acceleration and the Evolution of Spectral Distributions in Synchro-Self-Compton Sources: A Self-consistent Modeling of Blazars' Flares. Astrophys. J. **2011**, *739*, 66. [CrossRef]

14. 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*, 47. [CrossRef]

15. Abramowski, A. et al. [H.E.S.S. Collaboration]. A multiwavelength view of the flaring state of PKS 2155-304 in 2006. Astron. Astrophys., **2012**, *539*, 149.

16. Ghisellini, G.; Tavecchio, F. Canonical high-power blazars. Mon. Not. R. Astron. Soc. **2009**, *397*, 985. [CrossRef]

17. Lewis, T.R.; Finke, J.D.; Becker, P.A. Electron Acceleration in Blazars: Application to the 3C 279 Flare on 2013 December 20. Astrophys. J. **2019**, *884*, 116. [CrossRef]

18. Finke, J.D. External Compton scattering in blazar jets and the location of the gamma-ray emitting region. Astrophys. J. **2016**, *830*, 94. [CrossRef]

19. Joshi, M.; Böttcher, M. Time-dependent Radiation Transfer in the Internal Shock Model Scenario for Blazar Jets. Astrophys. J. **2011**, *727*, 21. [CrossRef]

20. Böttcher, M.; Baring, M.G. Multi-wavelength variability signatures of relativistic shocks in blazar jets. Astrophys. J. **2019**, *887*, 133. [CrossRef]

21. Dmytriiev, A.; Sol, H.; Zech, A. Connecting steady emission and very high energy flaring states in blazars: The case of Mrk 421. Mon. Not. R. Astron. Soc. **2021**, *505*, 2712. [CrossRef]

22. Katarzyński, K.; Sol, H.; Kus, A. The multifrequency variability of Mrk 421. Astron. Astrophys. **2003**, *410*, 101. [CrossRef]

23. 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. [CrossRef]

24. Aab, A. et al. [The Pierre Auger Collaboration]. Observation of a large-scale anisotropy in the arrival directions of cosmic rays above $8 \times 10^{18}$ eV. Science **2007**, *357*, 1266–1270.

25. The IceCube Collaboration; Fermi-LAT; MAGIC; AGILE; ASAS-SN; HAWC; H.E.S.S.; INTEGRAL; Kanata; Kiso; et al. Multi-messenger observations of a flaring blazar coincident with high-energy neutrino IceCube-170922A. Science **2018**, *361*, eaat1378. [CrossRef]

26. Celotti, A.; Ghisellini, G. The power of blazar jets. Mon. Not. R. Astron. Soc. **2008**, *385*, 283–300. [CrossRef]

27. Mannheim, K. The proton blazar. Astron. Astrophys. **1993**, *269*, 67.

28. Aharonian, F.A. TeV gamma rays from BL Lac objects due to synchrotron radiation of extremely high energy protons. New Astron. **2000**, *5*, 377. [CrossRef]

29. Mücke, A.; Protheroe, R.J. A proton synchrotron blazar model for flaring in Markarian 501. Astropart. Phys. **2001**, *15*, 121. [CrossRef]

30. Rodrigues, X.; Fedynitch, A.; Gao, S.; Boncioli, D.; Winter, W. Neutrinos and Ultra-high-energy Cosmic-ray Nuclei from Blazars. Astrophys. J. **2018**, *854*, 54. [CrossRef]

31. Mastichiadis, A.; Petropoulou, M.; Dimitrakoudis, S. Mrk 421 as a case study for TeV and X-ray variability in leptohadronic models. Mon. Not. R. Astron. Soc. **2013**, *434*, 2684. [CrossRef]

32. Zech, A.; Cerruti, M.; Mazin, D. Expected signatures from hadronic emission processes in the TeV spectra of BL Lacertae objects. Astron. Astrophys. **2017**, *602*, 25. [CrossRef]

33. 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. [CrossRef]

34. Cerruti, M. Leptonic and Hadronic Radiative Processes in Supermassive-Black-Hole Jets. Galaxies **2020**, *8*, 72. [CrossRef]

35. Petropoulou, M.; Dimitrakoudis, S.; Padovani, P. Photohadronic origin of γ-ray BL Lac emission: Implications for IceCube neutrinos. Mon. Not. R. Astron. Soc. **2015**, *448*, 2412. [CrossRef]

36. Reynoso, M.M.; Medina, M.C.; Romero, G.E. A lepto-hadronic model for high-energy emission from FR I radiogalaxies. Astron. Astrophys. **2011**, *531*, 30. [CrossRef]

37. Abdalla, H.; Adam, R.; Aharonian, F.; Benkhali, F.A.; Angüner, E.O.; Arakawa, M.; Arcaro, C.; Armand, C.; Ashkar, H.; Backes, M.; et al. Constraints on the emission region of 3C 279 during strong flares in 2014 and 2015 through VHE gamma-ray observations with H.E.S.S. Astron. Astrophys. **2019**, *627*, 159.

38. Böttcher, M.; Reimer, A.; Sweeney, K.; Prakash, A. Leptonic and Hadronic Modeling of Fermi-detected Blazars. Astrophys. J. **2013**, *768*, 54. [CrossRef]

39. Hervet, O.; Boisson, C.; Sol, H. Linking radio and gamma-ray emission in Ap Librae. Astron. Astrophys. **2015**, *578*, 69. [CrossRef]

40. Zacharias, M.; Wagner, S.J. The extended jet of AP Librae: Origin of the very high-energy γ-ray emission?. Astron. Astrophys. **2016**, *588*, 110. [CrossRef]

41. Petropoulou, M.; Vasilopoulos, G.; Giannios, D. The TeV emission of Ap Librae: A hadronic interpretation and prospects for CTA. Mon. Not. R. Astron. Soc. **2017**, *464*, 2213. [CrossRef]

42. Zdziarski, A.; Boettcher, M. Hadronic models of blazars require a change of the accretion paradigm. Mon. Not. R. Astron. Soc. **2015**, *450*, L21–L25. [CrossRef]

43. 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. [CrossRef]
44. Katarzyński, K.; Ghisellini, G.; Tavecchio, F. Hard TeV spectra of blazars and the constraints to the infrared intergalactic background. *Mon. Not. R. Astron. Soc.* **2006**, *368*, 52. [CrossRef]

45. Costamante, L.; Bonnoli, G.; Tavecchio, F. The NuSTAR view on hard-TeV BL Lacs. *Mon. Not. R. Astron. Soc.* **2018**, *477*, 4257. [CrossRef]

46. Giannios, D.; Uzdensky, D.A.; Begelman, M.C. Fast TeV variability in blazars: Jets in a jet. *Mon. Not. R. Astron. Soc.* **2009**, *395*, 29. [CrossRef]

47. Zech, A.; Lemoine, M. Electron-proton co-acceleration on relativistic shocks in extreme-TeV blazars. *Astron. Astrophys.* **2021**, *654*, 96. [CrossRef]

48. Tavecchio, F.; Costa, A.; Sciaccaluga, A. Extreme blazars: The result of unstable recollimated jets? *Mon. Not. R. Astron. Soc.* **2022**. [CrossRef]

49. Lefa, E.; Rieger, F.M.; Aharonian, F. Formation of Very Hard Gamma-Ray Spectra of Blazars in Leptonic Models. Astrophys. J. **2011**, 740, 64. [CrossRef]

50. Tavecchio, F. On the hadronic cascade scenario for extreme BL Lacs. *Mon. Not. R. Astron. Soc.* **2014**, *438*, 3255. [CrossRef]

51. Böttcher, M.; Dermer, C.D.; Finke, J.D. The Hard VHE γ-Ray Emission in High-Redshift TeV Blazars: Comptonization of Cosmic Microwave Background Radiation in an Extended Jet? Astrophys. J. **2008**, *679*, 9. [CrossRef]

52. Gaidos, J.A.; Akerlof, C.W.; Biller, S.; Boyle, P.J.; Breslin, A.C.; Buckley, J.H.; Carter-Lewis, D.A.; Catanese, M.; Cawley, M.F.; Fegan, D.J.; et al. Extremely rapid bursts of TeV photons from the active galaxy Markarian 421. Nature **1996**, *383*, 219-320. [CrossRef]

53. Aharonian, F. et al. [H.E.S.S. Collaboration]. An Exceptional Very High Energy Gamma-Ray Flare of PKS 2155-304. Astrophys. J. **2007**, 664, L71 [CrossRef]

54. Krawczynski, H.; Hughes, S.B.; Horan, D.; Aharonian, F.; Aller, M.F.; Aller, H.; Boltwood, P.; Buckley, J.; Coppi, P.; Fossati, G.; et al. Multiwavelength Observations of Strong Flares from the TeV Blazar 1ES 1959+650. Astrophys. J. **2004**, *601*, 151. [CrossRef]

55. Snios, B.; Nulsen, P.E.; Kraft, R.P.; Cheung, C.C.; Meyer, E.T.; Forman, W.R.; Jones, C.; Murray, S.S. Detection of Superluminal Motion in the X-ray Jet of M87. Astrophys. J. **2019**, *879*, 8. [CrossRef]

56. Snios, B.; Wykes, S.; Nulsen, P.E.J. Variability and Proper Motion of X-Ray Knots in the Jet of Centaurus A. Astrophys. J. **2019**, *871*, 2, 248. [CrossRef]

57. Dimitrakoudis, S.; Mastichiadis, A.; Protheroe, R.J.; Reimer, A. The time-dependent one-zone hadronic model. Astron. Astrophys. **2012**, *546*, 120. [CrossRef]

58. Boutelier, T.; Henri, G.; Petrucci, P.-O. An inhomogeneous jet model for the rapid variability of TeV blazars. Mon. Not. R. Astron. Soc. **2008**, *390*, L73. [CrossRef]

59. Potter, W.J.; Cotter, G. Synchrotron and inverse-Compton emission from blazar jets—II. An accelerating jet model with a geometry set by observations of M87. Mon. Not. R. Astron. Soc. **2013**, *429*, 1189. [CrossRef]

60. Marscher, A.P. Turbulent Extreme Multi-Zone Model for simulating flux and polarization variability in Blazars. Astrophys. J. **2014**, *780*, 87. [CrossRef]

61. Lucchini, M.; Markoff, S.; Crumley, P. Breaking degeneracy in jet dynamics: Multi-epoch joint modelling of the BL Lac PKS 2155-304. Mon. Not. R. Astron. Soc. **2019**, *482*, 4798. [CrossRef]

62. Aller, M.F.; Gomez, J.L.; Perlman, E. (Eds.) *Polarimetry as a Probe of Magnetic Fields in AGN Jets; Galaxies Special Issue; MDPI: Basel, Switzerland, 2021.*

63. The Fermi-LAT Collaboration and Members of the 3C 279 Multi-Band Campaign. A change in the optical polarization associated with a γ-ray flare in the blazar 3C 279. Nature **2010**, *463*, 919-923. [CrossRef] [PubMed]

64. Blinov, D.; Pavlidou, V.; Papadakis, I.; Kiehlmann, S.; Liodakis, I.; Panopoulou, G. V.; Angelakis, E.; Baloković, M.; Hovatta, T.; King, O. G. et al. RoboPol: Connection between optical polarization plane rotations and gamma-ray flares in blazars. Mon. Not. R. Astron. Soc. **2018**, *474*, 1296. [CrossRef]

65. Zhang, H.; Böttcher M. X-Ray and Gamma-Ray Polarization in Leptonic and Hadronic Jet Models of Blazars. Astrophys. J. **2013**, *774*, 18 [CrossRef]

66. Chiaberge, M.; Ghisellini, G. Rapid variability in the synchrotron self-Compton model for blazars. Mon. Not. R. Astron. Soc. **1999**, *306*, 551. [CrossRef]

67. Zacharias, M.; Reimer, A.; Boisson, C.; Zech, A. EXHALE-JET: An extended hadro-leptonic jet model for blazars—I. Code description and initial results. Mon. Not. R. Astron. Soc. **2022**, *512*, 3948. [CrossRef]

68. Sol, H.; Pelletier, G.; Asseo, E. Two-flow model for extragalactic radio jets. Mon. Not. R. Astron. Soc. **1989**, *237*, 411. [CrossRef]

69. Tavecchio, F.; Ghisellini, G. Spine-sheath layer radiative interplay in subparsec-scale jets and the TeV emission from M87. Mon. Not. R. Astron. Soc. **2008**, *385*, 98. [CrossRef]

70. Sikora, M.; Rutkowski, M.; Begelman, M.C. A spine-sheath model for strong-line blazars. Mon. Not. R. Astron. Soc. **2016**, *457*, 1352. [CrossRef]

71. Vuillaume, T.; Henri, G.; Petrucci, P.-O. A stratified jet model for AGN emission in the two-flow paradigm. Astron. Astrophys. **2018**, *620*, 41. [CrossRef]

72. Sikora, M. Powers and Magnetization of Blazar Jets. Galaxies **2016**, *4*, 12. [CrossRef]

73. Nalewajko, K.; Giannios, D.S.; Begelman, M.C. Radiative properties of reconnection-powered minijets in blazars. Mon. Not. R. Astron. Soc. **2011**, *413*, 333. [CrossRef]
74. Ghisellini, G.; Tavecchio, F. Rapid variability in TeV blazars: The case of PKS2155-304. *Mon. Not. R. Astron. Soc.* **2008**, *386*, 28. [CrossRef]

75. Barkov, M.V.; Aharonian, F.A.; Bogovalov, S.V. Rapid TeV Variability in Blazars as a Result of Jet-Star Interaction. *Astrophys. J.* **2012**, *749*, 119. [CrossRef]

76. Bosch-Ramon, V.; Perucho, M.; Barkov, M.V. Clouds and red giants interacting with the base of AGN jets. *Astron. Astrophys.* **2012**, *539*, 69. [CrossRef]

77. Gomez, J.L.; Marti, J.M.A.; Marscher, A.P. Parsec-Scale Synchrotron Emission from Hydrodynamic Relativistic Jets in Active Galactic Nuclei. *Astrophys. J.* **1995**, *449*, 19. [CrossRef]

78. Agudo, I.; Gómez, J.-L.; Martí, J.-M. Jet Stability and the Generation of Superluminal and Stationary Components. *Astrophys. J.* **2001**, *549*, 183. [CrossRef]

79. Mimica, P.; Aloy, M.-A.; Agudo, I. Spectral Evolution of Superluminal Components in Parsec-Scale Jets. *Astrophys. J.* **2009**, *696*, 1142. [CrossRef]

80. Porth, O.; Fendt, C.; Meliani, Z. Synchrotron Radiation of Self-collimating Relativistic Magnetohydrodynamic Jets. *Astrophys. J.* **2011**, *737*, 42. [CrossRef]

81. Fromm, C.M.; Perucho, M.; Mimica, P. Spectral evolution of flaring blazars from numerical simulations. *Astron. Astrophys.* **2016**, *588*, 101. [CrossRef]

82. Fichet de Clairfontaine, G.; Meliani, Z.; Zech, A. Flux variability from ejecta in structured relativistic jets with large-scale magnetic fields. *Astron. Astrophys.* **2021**, *647*, 77. [CrossRef]

83. Rieger, F.M. An Introduction to Particle Acceleration in Shearing Flows. *Galaxies* **2019**, *7*, 78. [CrossRef]

84. Petropoulou, M.; Sironi, L.; Spitkovsky, A. Relativistic Magnetic Reconnection in Electron-Positron-Proton Plasmas: Implications for Jets of Active Galactic Nuclei. *Astrophys. J.* **2019**, *880*, 37. [CrossRef]

85. Meyer, E.T.; Georganopoulos, M.; Sparks, W.B.; Perlman, E.; van der Marel, R.P.; Anderson, J.; Sohn, S.T.; Biretta, J.; Norman, C.; Chiaberge, M. A kiloparsec-scale internal shock collision in the jet of a nearby radio galaxy. *Nature* **2015**, *521*, 495–497. [CrossRef] [PubMed]

86. Fichet de Clairfontaine, G.; Meliani, Z.; Zech, A. Flare echos from relaxation shocks in perturbed relativistic jets. *Astron. Astrophys.* **2022**, *661*, A54. [CrossRef]

87. Schlickeiser, R. *Cosmic Ray Astrophysics*; Springer: Berlin/Heidelberg, Germany, 2002.

88. Stawarz, L.; Petrosian, V. On the Momentum Diffusion of Radiating Ultrarelativistic Electrons in a Turbulent Magnetic Field. *Astrophys. J.* **2008**, *681*, 1725. [CrossRef]

89. Zhdankin, V.; Werner, G.R.; Uzdensky, D.A.; Begelman, M.C. Kinetic Turbulence in Relativistic Plasma: From Thermal Bath to Nonthermal Continuum. *Phys. Rev. Lett.* **2017**, *118*, 055103. [CrossRef]

90. Zhdankin, V.; Uzdensky, D.A.; Werner, G.R.; Begelman, M.C. System-size Convergence of Nonthermal Particle Acceleration in Relativistic Plasma Turbulence. *Astrophys. J.* **2018**, *867*, L18. [CrossRef]

91. Comisso, L.; Sironi, L. Particle Acceleration in Relativistic Plasma Turbulence. *Phys. Rev. Lett.* **2018**, *121*, 255101. [CrossRef]

92. Nättilä, J.; Beloborodov, A.M. Radiative Turbulent Flares in Magnetically Dominated Plasmas. *Astrophys. J.* **2021**, *921*, 87. [CrossRef]

93. Bresci, V.; Lemoine, M.; Gremillet, L.; Comisso, L.; Sironi, L.; Demidemet, C. Nonresonant particle acceleration in strong turbulence: Comparison to kinetic and MHD simulations. *Phys. Rev. D* **2022**, *106*, 023028. [CrossRef]

94. Asano, K.; Hayashida, M. The Most Intensive Gamma-Ray Flare of Quasar 3C 279 with the Second-order Fermi Acceleration. *Astrophys. J. Lett.* **2015**, *808*, L18. [CrossRef]

95. Rieger, F.M.; Duffy, P. Turbulence and particle acceleration in shearing flows. *Astrophys. J. Lett.* **2021**, *907*, L2. [CrossRef]

96. Spitkovsky, A. Particle Acceleration in Relativistic Collisionless Shocks: Fermi Process at Last? *Astrophys. J. Lett.* **2008**, *682*, L5. [CrossRef]

97. Nishikawa, K.-I.; Hardee, P.; Richardson, G.; Preece, R.; Sol, H.; Fishman, G. J. Particle Acceleration and Magnetic Field Generation in Electron-Positron Relativistic shocks. *Astrophys. J.* **2005**, *622*, 927. [CrossRef]

98. Sironi, L.; Spitkovsky, A. Particle Acceleration in Relativistic Magnetized Collisionless Electron-Ion Shocks. *Astrophys. J.* **2011**, *726*, 75. [CrossRef]

99. Summerlin, E.J.; Baring, M.G. Diffusive acceleration of particles at oblique, relativistic, magnetohydrodynamic shocks. *Astrophys. J.* **2012**, *745*, 63. [CrossRef]

100. Ellison, D.C.; Warren, D.C.; Bykov, A.M. Monte Carlo simulations of nonlinear particle acceleration in parallel trans-relativistic shocks. *Astrophys. J.* **2013**, *776*, 46. [CrossRef]

101. Guo, X.; Sironi, L.; Narayan, R. Non-thermal electron acceleration in low Mach number collisionless shocks. I. Particle energy spectra and acceleration mechanism. *Astrophys. J.* **2014**, *794*, 153. [CrossRef]

102. Sironi, L.; Keshet, U.; Lemoine, M. Relativistic shocks: Particle acceleration and magnetization. *Space Sci. Rev.* **2015**, *191*, 519. [CrossRef]

103. Groselj, D.; Sironi, L.; Beloborodov, A.M. Microphysics of relativistic collisionless electron-ion-positron shocks. *Astrophys. J.* **2022**, *933*, 74. [CrossRef]

104. Lemoine, M. Generalized Fermi acceleration. *Phys. Rev. D* **2019**, *99*, 083006. [CrossRef]
105. Treumann, R.A.; Baumjohann, W. Spontaneous magnetic reconnection—Collisionless reconnection and its potential astrophysical relevance. *Astron. Astrophys. Rev.* **2015**, *23*, 4. [CrossRef]

106. Marcowith, A.; Ferrand, G.; Grech, M.; Meliani, Z.; Plotnikov, I.; Walder, R. Multi-scale simulations of particle acceleration in astrophysical systems. *Living Rev. Comput. Astrophys.* **2020**, *6*, 1. [CrossRef]

107. Nishikawa, K.; Dutan, I.; Kohn, C.; Mizuno, Y. PIC methods in astrophysics: Simulations of relativistic jets and kinetic physics in astrophysical systems. *Living Rev. Comput. Astrophys.* **2021**, *7*, 1. [CrossRef] [PubMed]

108. Zenitani, S.; Hoshino, H. The role of the guide field in relativistic pair plasma reconnection. *Astrophys. J.* **2008**, *677*, 530. [CrossRef]

109. Sironi, L.; Spitkovsky, A. Relativistic reconnection: An efficient source of non-thermal particles. *Astrophys. J. Lett.* **2014**, *783*, L21. [CrossRef]

110. Sironi, L.; Giannios, D. Petropoulou, Plasmoids in relativistic reconnection, from birth to adulthood: First they grow, then they go. *Mon. Not. R. Astron. Soc.* **2016**, *462*, 48. [CrossRef]

111. Werner, G.R.; Uzdensky, D.A. Nonthermal particle acceleration in 3D relativistic magnetic reconnection in pair plasma. *Astrophys. J. Lett.* **2017**, *843*, L27. [CrossRef]

112. Christie, I.M.; Petropoulou, M.; Sironi, L.; Giannios, D. Interplasmoid Compton scattering and the Compton dominance of BL Lacs. *Mon. Not. R. Astron. Soc.* **2020**, *492*, 549. [CrossRef]

113. Nishikawa, K.; Mizuno, Y.; Gomez, J.L.; Ioana, D.; Jacek, N.; Oleh, K.; Nicholas, M.; Athina, M.; Martin, P.; Kouichi, H. Rapid particle acceleration due to recollimation shocks and turbulent magnetic fields in injected jets with helical magnetic fields. *Astrophys. J.* **2020**, *493*, 2652. [CrossRef]

114. Zhang, H.; Sironi, L.; Giannios, D. Fast particle acceleration in three-dimensional relativistic reconnection. *Astrophys. J.* **2021**, *922*, 261. [CrossRef]

115. Sironi, L.; Rowan, M.E.; Narayan, R. Reconnection-driven particle acceleration in relativistic shear flows. *Astrophys. J. Lett.* **2021**, *907*, L44. [CrossRef]

116. Meli, A.; Nishikawa, K.; Kohn, C.; Dutan, I.; Mizuno, Y.; Kobzar, O.; MacDonald, N.; Gomez, J. L.; Hirotani, K. 3D PIC simulations for relativistic jets with a toroidal magnetic field. *Mon. Not. R. Astron. Soc.* **2022**, submitted.

117. Bransgrove, A.; Ripperda, B.; Philippov, A. Magnetic Hair and Reconnection in Black Hole Magnetospheres. *Phys. Rev. Lett.* **2021**, *127*, 055101. [CrossRef] [PubMed]

118. Crinquand, B.; Cerutti, B.; Dubus, G.; Parfrey, K.; Philippov, A. Images of magnetospheric reconnection-powered radiation around supermassive black holes. *arXiv* **2022**, arXiv:2202.04472.

119. Nathanail, A.; Mpisketzis, V.; Porth, O.; Fromm, C.M.; Rezzolla, L. Magnetic reconnection and plasmoid formation in three-dimensional accretion flows around black holes. *Mon. Not. R. Astron. Soc.* **2022**, *513*, 4267. [CrossRef]

120. Medina-Torrejón, T.E.; de Gouveia Dal Pino, E.M.; Kadowaki, L.H.S.; Kowal, G.; Singh, C.B.; Mizuno, Y. Particle Acceleration by Relativistic Magnetic Reconnection Driven by Kink Instability Turbulence in Poynting Flux-Dominated Jets. *Astrophys. J.* **2021**, *908*, 193. [CrossRef]

121. Morris, P.J.; Potter, W.J.; Cotter, G. The feasibility of magnetic reconnection powered blazar flares from synchrotron self-Compton emission. *Mon. Not. R. Astron. Soc.* **2019**, *486*, 1548. [CrossRef]

122. Asano, K.; Hayashida, M. Synchrotron gamma-ray emission model of the giant outburst of the quasar 3C 279 in 2015: Fast reconnection or stochastic acceleration with electromagnetic cascade? *Astrophys. J.* **2020**, *890*, 56. [CrossRef]

123. Sobacchi, E.; Nattila, J.; Sironi, L. A fully kinetic model for orphan gamma-ray flares in blazars. *Mon. Not. R. Astron. Soc.* **2021**, *503*, 688. [CrossRef]

124. Goldreich, P.; Julian, W.H. Pulsar Electrodynamics. *Astrophys. J.* **1969**, *157*, 869. [CrossRef]

125. Gangadhara, R.T.; Lesch, H. On the nonthermal emission in active galactic nuclei. *Astron. Astrophys.* **1997**, *323*, L45.

126. Osmanov, Z.; Rogava, A.; Bodo, G. On the efficiency of particle acceleration by rotating magnetospheres in AGN. *Astron. Astrophys.* **2007**, *470*, 395. [CrossRef]

127. Beskin, V.S.; Istomin, Y.N.; Pariev, V.I. Filling the magnetosphere of a supermassive black hole with plasma. *Sov. Astron.* **1992**, *36*, 64.

128. Hirotani, K.; Okamoto, I. Pair plasma production in a force-free magnetosphere around a supermassive black hole. *Astrophys. J.* **1998**, *497*, 563. [CrossRef]

129. Neronov, A.; Aharonian, F.A. Production of TeV gamma radiation in the vicinity of the supermassive black hole in the giant radio galaxy M87. *Astrophys. J.* **2007**, *671*, 85. [CrossRef]

130. Rieger, F.M.; Aharonian, F.A. Variable VHE gamma-ray emission from non-blazar AGNs. *Astron. Astrophys.* **2008**, *479*, L5. [CrossRef]

131. Istomin, Y.; Sol, H. Acceleration of particles in the vicinity of a massive black hole. *Astrophys. Space Sci.* **2009**, *321*, 57. [CrossRef]

132. Osmanov, Z. Is very high energy emission from the BL Lac 1ES 0806+524 centrifugally driven? *New Astron.* **2010**, *15*, 351. [CrossRef]

133. Broderick, A.E.; Tchekhovskoy, A. Horizon-scale lepton acceleration in jets: Explaining the compact radio emission in M87. *Astrophys. J.* **2015**, *809*, 97. [CrossRef]

134. Katsoulakos, G.; Rieger, F.M. Magnetospheric Gamma-Ray Emission in Active Galactic Nuclei. *Astrophys. J.* **2018**, *852*, 112. [CrossRef]
135. Parfrey, K.; Philippov, A.; Cerruti, B. First-principles plasma simulations of black hole jet launching. *Phys. Rev. Lett.* **2019**, *122*, 03521. [CrossRef]

136. Crinquand, B.; Cerutti, B.; Philippov, A.; Parfrey, K.; Dubus, G. Multidimensional simulations of ergospheric pair discharges around black holes. *Phys. Rev. Lett.* **2020**, *124*, 145101. [CrossRef]

137. Katsoulakos, G.; Rieger, F.M. Gap-type Particle Acceleration in the Magnetospheres of Rotating Supermassive Black Holes. *Astrophys. J.* **2020**, *895*, 99. [CrossRef]

138. Istomin, Y.; Gunya, A. Acceleration of the high-energy protons in active galactic nuclei. *Astron. Nachr.* **2021**, *342*, 182. [CrossRef]

139. Aleksić, J.; Ansoldi, S.; Antonelli, L.A.; Antoranz, P.; Babic, A.; Bangale, P.; Barrio, J.A.; Becerra González, J.; Bednarek, W.; Bernardini, E.; et al. Black hole lightning due to particle acceleration at subhorizon scales. *Science* **2014**, *346*, 1080–1084. [CrossRef]

140. Crinquand, B.; Cerutti, B.; Dubus, G.; Parfrey, K.; Philippov, A. Synthetic gamma-ray lightcurves of Kerr black hole magnetospheric activity from particle-in-cell simulations. *Astron. Astrophys.* **2021**, *650*, A163. [CrossRef]

141. Kisaka, S.; Levinson, A.; Toma, K.; Niv, I. The Response of Black Hole Spark Gaps to External Changes: A Production Mechanism of Rapid TeV Flares? *Astrophys. J.* **2022**, *924*, 22. [CrossRef]

142. Akiyama, K.; Lu, R.; Fish, V.L.; Doeleman, S.S.; Broderick, A.E.; Dexter, J.; Hada, K.; Kino, M.; Nagai, H.; Honma, M. et al. 230 GHz VLBI observations of M87: Event-horizon-scale structure during an enhanced very-high-energy gamma-ray state in 2012. *Astrophys. J.* **2015**, *807*, 150. [CrossRef]

143. Jorstad, S.; Marscher, A. The VLBA-BU-BLAZAR multiwavelength monitoring program. *Galaxies* **2016**, *4*, 47. [CrossRef]

144. Rani, B.; Jorstad, S. G.; Marscher, A. P.; Agudo, I.; Sokolovsky, K. V.; Larionov, V. M.; Smith, P.; Mosunova, D. A.; Borman, G. A.; Grishina, T. S. et al. Exploring the connection between parsec-scale jet activity and broadband outbursts in 3C 279. *Astrophys. J.* **2018**, *858*, 80. [CrossRef]

145. Larionov, V.M.; Jorstad, S.G.; Marscher, A.P.; Villata, M.; Raiteri, C.M.; Smith, P.S.; Agudo, I.; Savchenko, S.S.; Morozova, D.A.; Acosta-Pulido, J.A. et al. Multiwavelength behaviour of the blazar 3C 279: Decade-long study from gamma-ray to radio. *Mon. Not. R. Astron. Soc.* **2020**, *492*, 3829. [CrossRef]

146. Kim, J.-Y.; Krichbaum, T.P.; Broderick, A.E.; Wielgus, M.; Blackburn, L.; Gómez, J.L.; Johnson, M.D.; Bouman, K.L.; Chael, A.; Akiyama, K.et al. Event Horizon Telescope imaging of the archetypal blazar 3C 279 at an extreme 20 microarcsecond resolution. *Astron. Astrophys.* **2020**, *640*, A69. [CrossRef]

147. Lico, R.; Casadio, C.; Jorstad, S.G.; Gomez, J.L.; Marscher, A.P.; Traianou, E.; Kim, J.Y.; Zhao, G.Y.; Fuentes, A.; Cho, I.; et al. New jet feature in the parsec-scale jet of the blazar OJ 287 connected to the 2017 teraelectronvolt flaring activity. *Astron. Astrophys.* **2022**, *658*, L10. [CrossRef]

148. Hervet, O.; Williams, D.A.; Falcone, A.D.; Kaur, A. Probing an X-Ray Flare Pattern in Mrk 421 Induced by Multiple Stationary Shocks: A Solution to the Bulk Lorentz Factor Crisis. *Astrophys. J.* **2019**, *877*, 26. [CrossRef]

149. Abdalla, H.; Adam, R.; Aharonian, F.; Benkhali, F.A.; Angüner, E.O.; Arcaro, C.; Armand, C.; Armstrong, T.; Ashkar, H.; Backes, M.; et al. H.E.S.S. and MAGIC observations of a sudden cessation of a very-high-energy γ-ray flare in PKS 1510-089 in May 2016. *Astron. Astrophys.*, **2021**, *648*, A23.

150. The H.E.S.S. Collaboration. Resolving acceleration to very high energies along the jet of Centaurus A. *Nature* **2020**, *582*, 356–359.

151. Bednarek, W. GeV-TeV γ-rays produced by electrons in the kpc-scale jet as a result of Comptonization of the inner jet emission. *Mon. Not. R. Astron. Soc.* **2019**, *483*, 1003. [CrossRef]

152. Stawarz, L.; Sikora, M.; Ostrowski, M. High-Energy Gamma Rays from FR I Jets A. *Astrophys. J.* **2003**, *597*, 186. [CrossRef]

153. Costamante, L.; Cutini, S.; Tosti, G.; Antolini, E.; Tramacere, A. On the origin of gamma-rays in Fermi blazars: Beyond the broad-line region. *Mon. Not. R. Astron. Soc.* **2018**, *477*, 4749. [CrossRef]

154. Foffano, L.; Vittorini, V.; Tavani, M.; Menegoni, E. Absorption Features in Sub-TeV Gamma-Ray Spectra of BL Lac Objects. *Astrophys. J.* **2022**, *926*, 95. [CrossRef]

155. Isler, J.C.; Urry, C.M.; Bailyn, C. The SMARTS multi-epoch optical spectroscopy atlas (SaMOSA): An analysis of emission line variability in southern hemisphere Fermi blazars. *Astrophys. J.* **2015**, *804*, 7. [CrossRef]

156. Aharonian, F.A.; Barkov, M.V.; Khangulyan, D. Scenarios for Ultrafast Gamma-Ray Variability in AGN. *Astron. Astrophys.* **2012**, *539*, 149. [CrossRef]

157. Meyer, M.; Scargle, J.D.; Blandford, R.D. Characterizing the Gamma-Ray Variability of the Brightest Flat Spectrum Radio Quasars Observed with the Fermi LAT. *Astrophys. J.* **2019**, *877*, 39. [CrossRef]

158. Narayan, R.; Piran, T. Variability in blazars: Clues from PKS 2155-304. *Mon. Not. R. Astron. Soc.* **2012**, *420*, 604. [CrossRef]

159. 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. [CrossRef]

160. Levy, C.; Sol, H.; Bolmont, J. Modeling intrinsic time-lags in flaring blazars in the context of Lorentz Invariance Violation searches. In Proceedings of the 37th International Cosmic Ray Conference, Berlin, Germany, 12–23 July 2021; p. 907. Available online: https://pos.sissa.it/cgi-bin/reader/conf.cgi?confid=395 (accessed on 31 July 2022).

161. Roy, A.; Sarkar, A.; Chatterjee, A.; Gupta, A.C.; Chitnis, V.R.; Wiita, P.J. Transient quasi-periodic oscillations at γ-rays in the TeV blazar PKS 1510-089. *Mon. Not. R. Astron. Soc.* **2022**, *510*, 3641. [CrossRef]

162. Hervet, O. et al. [VERITAS Collaboration]. A multiwavelength look at the 2017 flare of OJ287. In Proceedings of the 7th Heidelberg International Symposium on High-Energy Gamma-Ray Astronomy, Barcelona, Spain, 4–8 July 2022. Available online: https://indico.icc.ub.edu/event/46/contributions/1257/ (accessed on 31 July 2022).
163. Zech, A.; Mazin, D.; Biteau, J.; Kishimoto, M.; Pott, J.-U.; Almeida, C.R.; Surdej, J.; Tristram, K.R.W. Key Science Project: Active Galactic Nuclei. In *Science with the Cherenkov Telescope Array*; World Scientific: Singapore, 2019; p. 231.

164. O'Brien, S. et al. [VERITAS collaboration]. VERITAS detection of VHE emission from the optically bright quasar OJ 287. In Proceedings of the 35th International Cosmic Ray Conference, Busan, Korea, 10–20 July 2017.

165. Acciari, V.A. et al. [MAGIC Collaboration]. Multiwavelength study of the gravitationally lensed blazar QSO B0218+357 between 2016 and 2020. *Mon. Not. R. Astron. Soc.* **2022**, *510*, 2344. [CrossRef]

166. Goldoni, P.; Pita, S.; Boisson, C.; Max-Moerbeck, W.; Kasai, E.; Williams, D. A.; D'Ammando, F.; Navarro-Aranguiz, V.; Backes, M.; Barres de Almeida, U. et al. Optical spectroscopy of blazars for the Cherenkov Telescope Array. *Astron. Astrophys.* **2021**, *650*, A106. [CrossRef]

167. Ehlert, S.E.; Ferrazzoli, R.; Marinucci, A.; Marshall, H.L.; Middei, R.; Pacciani, L.; Perri, M.; Petrucci, P.-O.; Puccetti, S.; Barnouin, T.; et al. Limits on X-ray polarization at the core of Centaurus A as observed with the imaging X-ray polarimeter explorer. *Astrophys. J.* **2022**, *935*, 116. [CrossRef]

168. Peirson, A.L.; Liodakis, I.; Romani, R.W. Testing high-energy emission models for blazars with X-ray polarimetry. *Astrophys. J.* **2022**, *931*, 59. [CrossRef]

169. Peirson, A.L.; Romani, R.W. The polarization behavior of relativistic synchrotron self-Compton jets. *Astrophys. J.* **2019**, *885*, 76. [CrossRef]

170. Marscher, A.P.; Jorstad, S.G. Frequency and time dependence of linear polarization in turbulent jets of blazars. *Galaxies* **2021**, *9*, 27. [CrossRef]

171. Zhang, H.; Li, X.; Giannios, D.; Guo, F. First-principles prediction of X-ray polarization from magnetic reconnection in high-frequency BL Lacertae objects. *Astrophys. J.* **2021**, *912*, 129. [CrossRef]

172. Di Gesu, L.; Tavecchio, F.; Donnarumma, I.; Marscher, A.; Pesce-Rollins, M.; Landoni, M. Testing particle acceleration models for BL Lac jets with the imaging X-ray polarimeter explorer. *Astron. Astrophys.* **2022**, *662*, A83. [CrossRef]

173. Ramos Almeida, C.; Ricci, C. Nuclear obscuration in active galactic nuclei. *Nat. Astron.* **2017**, *1*, 679. [CrossRef]

174. Hoenig, S.; Kishimoto, M. Dusty winds in active galactic nuclei: Reconciling observations with models. *Astrophys. J. Lett.* **2017**, *838*, L20. [CrossRef]

175. Engel, K.; Goodman, J.; Huentemeyer, P.; Kierans, C.; Lewis, T.R.; Negro, M.; Santander, M.; Williams, D.A.; Allen, A.; Aramaki, T.; et al. Gamma-Ray Experiments: The Future of Gamma-Ray Experiments in the MeV-EeV Range. In Proceedings of the US Community Study on the Future of Particle Physics (Snowmass 2021), Seattle, WA, USA, 17–27 July 2022. Available online: https://arxiv.org/pdf/2203.07360.pdf (accessed on 31 July 2022).

176. Mereghetti, S.; Balman, S.; Caballero-Garcia, M.; Del Santo, M.; Doroshenko, V.; Erkut, M. H.; Hanlon, L.; Hoeflich, P.; Markowitz, A.; Osborne, J.P.; et al. Time domain astronomy with the THESEUS satellite. *arXiv* **2021**, arXiv:2104.09533v1.

177. Venturi, T.; Paragi, Z.; Lindqvist, M. (Eds.) VLBI20-30: *A Scientific Roadmap for the Next Decade; The Future of the European VLBI Network*; JIVE: Dwingeloo, The Netherlands; EVN: Nordrhein-Westfalen, Germany, 2020. Available online: https://arxiv.org/ftp/arxiv/papers/2007/2007.02347.pdf (accessed on 28 October 2022).

178. Kovalev, Y.Y.; Zobnina, D.I.; Plavin Blinov, D. Optical polarization properties of AGNs with significant VLBI-Gaia offsets. *Mon. Not. R. Astron. Soc.* **2020**, *493*, L54 [CrossRef]

179. Lambert, S.; Liu, N.; Arias, E.F.; Barache, C.; Souchay, J.; Taris, F.; Liu, J. C.; Zhu, Z. Parsec-scale alignments of radio-optical offsets with jets in AGNs from multifrequency geodetic VLBI, Gaia EDR3, and the MOJAVE program. *Astron. Astrophys.* **2021**, *651*, A64. [CrossRef]

180. Lambert, S.; Pierron, A.; Sol, H. Locating the blazar gamma-ray zone from astrometric VLBI and Gaia data? In Proceedings of the GAMMA2022, 7th Heidelberg International Symposium on High Energy gamma-ray Astronomy, Barcelona, Spain, 4–8 July 2022. to be published. Available online: https://pos.sissa.it/417/ (accessed on 31 July 2022).

181. Coleman, A.; Eser, J.; Mayotte, E.; Sarazin, F.; Schröder, F. G.; Soldin, D.; Venters, T. M.; Aloisio, R.; Alvarez-Muñiz, J.; Batista, R.A.; et al. Ultra-High-Energy Cosmic Rays: The Intersection of the Cosmic and Energy Frontiers. In Proceedings of the US Community Study on the Future of Particle Physics (Snowmass 2021), Seattle, WA, USA, 17–27 July 2022. Available online: https://arxiv.org/pdf/2205.05845.pdf (accessed on 31 July 2022).

182. Matthews, J.H.; Taylor, A.M. Particle acceleration in radio galaxies with flickering jets: GeV electrons to ultrahigh energy cosmic rays. *Mon. Not. R. Astron. Soc.* **2021**, *503*, 5948. [CrossRef]

183. Cerruti, M. Neutrinos from blazars. *J. Pharm. Health Care Sci.* **2020**, *1468*, 2094. [CrossRef]

184. IceCube Collaboration; Aartsen, M.; Ackermann, M.; Adams, J.; Aguilar, J.A.; Ahlers, M.; Ahrens, M.; Samarai, I.A.; Altmann, D.; Andeen, K.; et al. Neutrino emission from the direction of the blazar TXS 0506+056 prior to the IceCube-170922A alert. *Science* **2018**, *361*, 147–151.

185. Plavin, A.; Kovalev, Y.Y.; Kovalev, Y.A.; Troitsky, S. Observational Evidence for the Origin of High-energy Neutrinos in Parsec-scale Nuclei of Radio-bright Active Galaxies. *Astrophys. J.* **2020**, *894*, 101. [CrossRef]

186. Kun, E.; Bartos, I.; Becker Tjus, J. Cosmic Neutrinos from Temporarily Gamma-suppressed Blazars. *Astrophys. J. Lett.* **2021**, *911*, L18. [CrossRef]

187. Aartsen, M.G.; Ackermann, M.; Adams, J.; Aguilar, J.A.; Ahlers, M.; Ahrens, M.; Alispach, C.; Andeen, K.; Anderson, T.; Ansseau, I.; et al. Time-integrated Neutrino Source Searches with 10 years of IceCube Data. *Phys. Rev. Lett.* **2020**, *124*, 051103. [CrossRef]
188. Graham, M.J.; Ford, K.E.S.; McKernan, B.; Ross, N. P.; Stern, D.; Burdge, K.; Coughlin, M.; Djorgovski, S. G.; Drake, A. J.; Duev, D.; et al. Candidate Electromagnetic Counterpart to the Binary Black Hole Merger Gravitational Wave Event S190521g. *Phys. Rev. Lett.* **2020**, *124*, 251102. [CrossRef]

189. Samsing, J.; Bartos, I.; D'Orazio, D.J.; Haiman, Z.; Kocsis, B.; Leigh, N. W. C.; Liu, B.; Pessah, M. E.; Tagawa, H. AGN as Potential Factories for Eccentric Black Hole Mergers. *Nature* **2022**, *603*, 237–240. [CrossRef]

190. McKernan, B.; Ford, K.E.S.; Bartos, I. Ram-pressure Stripping of a Kicked Hill Sphere: Prompt Electromagnetic Emission from the Merger of Stellar Mass Black Holes in an AGN Accretion Disk. *Astrophys. J. Lett.* **2019**, *884*, L50. [CrossRef]

191. Murase, K.; Kashiyama, K.; Mészáros, P. Ultrafast outflows from black hole mergers with a minidisk. *Astrophys. J. Lett.* **2016**, *822*, L9. [CrossRef]

192. Bartos, I.; Kocsis, B.; Haiman, Z.; Márka, S. Rapid and Bright Stellar-mass Binary Black Hole Mergers in Active Galactic Nuclei. *Astrophys. J.* **2017**, *835*, 165. [CrossRef]