Leptonic and Hadronic Radiative Processes in Supermassive-Black-Hole Jets

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


Active Galactic Nuclei
The brightest, persistent objects in the Universe are active galactic nuclei (AGNs), compact regions in the center of galaxies that can outshine the host itself (composed of hundreds of billions of stars), and can be detected on cosmological distances, up to a redshift z = 7.5 [1], which means a comoving distance of about 9 Gpc in the current ΛCDM cosmological model [2]. About 75 years of observations at all wavelengths have collectively shaped the so-called AGN unified model: the engine powering the system is a supermassive black hole (SMBH) of 10 8−9 M [3], which accretes matter in the form of an accretion disk [4] and which is surrounded at the parsec scale by an obscuring dusty torus [5]. The thermal emission from the disk ionizes clouds of matter orbiting the SMBH, which then re-emit emission lines in optical/UV (the broad and narrow line-regions, BLR/NLR, the first one close to the black hole at the sub-parsec scale [6,7], the second one at much larger distances, hundreds of parsecs and more [8], beyond the torus, see a schematic representation in Figure 1). The inclination of the system with respect to the observer is responsible for the variety of AGN types, due to the fact that when the system is seen face-on the observer has direct access to the SMBH, while when the system is seen edge-on the inner regions are obscured by the torus [9]. Quasars (quasi-stellar objects) are now understood as luminous AGNs seen face-on. The host galaxy of the first identified quasars (such as 3C 273 [10]) can now be resolved and studied with current instrumentation [11].
A major dichotomy exists among AGNs: a minority of around 10% of AGNs are bright in the radio band [12] and are hence dubbed radio-loud AGNs, in contrast with the more common radio-quiet AGNs. The radio non-thermal emission is associated with a pair of relativistic jets of plasma which are launched along the polar axis of the SMBH and perpendicularly to the accretion disk. These jets can travel up to the Mpc scale, exceeding the (visible) size of the AGN host galaxy. For the nearest and brightest objects, the emission from the relativistic jet has been resolved not only in the radio band [13], but also in optical [14], X-rays [15], and more recently in γ-rays [16,17]. AGNs with extended radio emission from the jets are called radio-galaxies. In this case the observer sees the jet from the side. When the relativistic jet points instead in the direction of the observer, relativistic effects boost the emission [18,19], making these AGNs particularly luminous in the Universe. These radio-loud AGNs observed down-the-jet are called blazars [20], and are characterized by high luminosity, rapid variability, and high degree of polarization [21]. Although blazars are a minority among AGNs their high luminosity makes them ubiquitous in all extragalactic surveys at all wavelengths. It is remarkable that in γ-rays, where the radio-quiet AGNs do not emit anymore, around 70% of all the known extragalactic sources are blazars [22].

Blazars
Blazars come in two flavors: BL Lacertae objects, whose optical/UV spectrum is characterized by a featureless non-thermal continuum [23], and Flat-Spectrum-Radio-Quasars (FSRQs), which show instead broad emission lines from the BLR typical of quasars [24]. This dichotomy is also observed in the radio-galaxies parent population, which are divided into Fanaroff-Riley type I radio-galaxies, with strong core emission and relatively fainter jets, and Fanaroff-Riley type II radio-galaxies, with much stronger jets that terminate in a shock with the inter-galactic medium [25].
The fact that the relativistic jet of plasma moves towards the observer with speed v = βc, Lorentz factor Γ = (1 − β 2 ) −0.5 , and viewing angle ϑ view has three important consequences coming directly from special relativity: • plasma moving in the jet will be seen in the observer's frame as showing a projected superluminal speed β project > 1. It can be shown that β project = (β sin ϑ view )/(1 − β cos ϑ view ), which can be as high as βΓ. Indeed, the detection of superluminal motion in blazar jets is one of the key observations supporting the AGN unified model [26]. • the emission from the jet is boosted in the direction of movement, which translates into a flux density in the observer's frame 1 which scales as F ν = δ 3 F ν , where the quantity δ = [Γ(1 − β cos ϑ view )] −1 is defined as the Doppler factor. Typical values for the Doppler factor in blazars are of the order of ten, which means that the flux density in the observer's frame is typically a factor of a thousand higher than in the jet's frame. Similarly, it means that if we have a radio-galaxy and a blazar at the same redshift, and with identical flux densities in the jet's frame, the blazar will be a thousand times brighter than the radio-galaxy as observed from Earth. • the effect of time-compression translates into a flux variability in the observer's frame much faster than the one in the jet's frame: the variability timescale, conventionally defined as the flux-doubling timescale τ var , is equal to (1 + z) τ var /δ. For a typical blazar Doppler factor of ten, a flare with a timescale in the jet's frame of one hour, is compressed in the observer's frame into a minutes-long flare.
The spectral energy distribution (SED) 2 of blazars is characterized by a non-thermal continuum extending from radio to very-high-energy γ-rays, with only subdominant thermal emission in optical/UV associated with the accretion disk (the big-blue-bump) and emission lines from the BLR (only seen in FSRQs, and not in BL Lacs). This non-thermal continuum is comprised of two distinct radiative components: the first one peaks in infrared-to-X-rays, the second one peaks in the γ-ray band, ranging from MeV to TeV. A remarkable discovery by γ-ray telescopes is that in FSRQs the majority of the energy is emitted in the γ-ray band [27], and in BL Lacs the energy emitted in the γ-ray band is as much as the one emitted in the radio-to-X-ray band [28]. The SED peak frequencies are not fixed, not within the blazar population, nor in a single blazar. FSRQs show typically a first SED peak in the infrared band, and a second SED peak in the MeV band; BL Lacs, on the other hand, show a huge variety of peak frequencies, from infrared to X-rays, and from MeV to TeV. The peak frequency is used to further classify BL Lacs [29,30] into low-frequency-peaked BL Lacs (LBLs, with ν peak < 10 14 Hz, in the infrared band, similarly to FSRQs) and high-frequency-peaked BL Lacs (HBLs, with ν peak > 10 15 Hz, in the UV/X-ray band). The transition among the two sub-classes is smooth, and there exist intermediate-frequency-peaked BL Lacs (IBLs, with 10 14 Hz < ν peak < 10 15 Hz, in optical, [31]). It is not clear how the peak frequency distribution ends: observations by Cherenkov telescopes in the TeV band have identified a population of extremely high-frequency-peaked blazars (EHBLs) with ν peak > 10 18 Hz, in the hard-X-rays [32,33], and, correspondingly, a high-energy SED peak extending into the TeV band [34]. More recent studies seem to indicate however that the EHBL population is not homogeneous, with all various combinations of soft/hard X-rays and soft/hard TeV spectra possible [35][36][37]. It seems that the peak frequency is anti-correlated with the blazar luminosity: the brightest blazars are the ones with the lowest ν peak , while the faintest ones are the ones with the highest ν peak . This anti-correlation defines the blazar sequence [38][39][40]. The existence of this sequence is still an open question in the blazar research field, and several outliers to the sequence have been identified [41][42][43][44]. If the blazar sequence represents a true physical property of blazars, it will help understanding what physical processes shape the blazar SEDs [45][46][47]. 1 Here and in the following, primed quantities refer to the jet's frame, while non-primed quantities to the observer's frame. 2 We define SED the energy density representation νF ν , where F ν is the flux density in units of erg cm −2 s −1 Hz −1 . It is easy to show that while F ν is different from F λ (the flux density per unit wavelength) and F E (the flux density per unit energy), the SED is constant and νF ν = λF λ = EF E .
The framework described up to here is pretty solid, and has been tested against a variety of observations at all wavelengths. Nonetheless, there is still many questions which remain open. Only limiting the discussion to radio-loud AGNs, it is still not completely understood (i) how relativistic jets are formed and powered, and why only some AGNs can produce them; (ii) why there are two intrinsically different jet types as seen in FRI/FRII and BL Lacs/FSRQs; (iii) how particles are accelerated in jets; (iv) how photons are produced in jets. The last item is the subject of this review. I will summarize the state-of-the-art answer to the following question: how do supermassive-black-hole jets shine? For a review of recent progress on observational results see [48,49]. For a broader discussion of relativistic jets at all scales (AGNs and micro-quasars) see [50].
The blazar non-thermal continuum emission implies the existence of a non-thermal population of particles that can radiate photons over the whole electromagnetic spectrum. Leptonic radiative processes are the ones associated with electrons/positrons in the jet; hadronic radiative processes are the ones associated with protons and nuclei. In the following I will discuss how photons (and as it will be shown later, neutrinos) can be emitted from jets, assuming that there exists an efficient particle accelerator in the jet. For particle acceleration in AGN jets, see [51,52]

Electron Synchrotron Emission
Synchrotron radiation is produced by charged particles moving in a magnetic field. For sake of simplicity, let's start with the assumption that the emitting region is a sphere of radius R in the jet, moving with bulk Lorentz factor Γ, and filled with a tangled, homogeneous magnetic field B . If a particle with electric charge e and Lorentz factor γ = E /mc 2 is put in the region, it radiates a synchrotron power (in units of erg s −1 Hz −1 ) per frequency ν equal to 3 where ϕ is the pitch angle between the particle momentum and the magnetic field lines, the frequency ν c is defined as ν c (γ , ϕ ) = 3eB 4πmc 2 γ 2 sin ϕ , and K 5/3 (x) is the modified Bessell function of the second kind of order 5/3. The sign of the charge does not matter: the synchrotron emission from electrons and positrons is not distinguishable, and in the following I will simply talk about electrons.
Let's now assume that in the emitting region there is not one single electron, but an isotropic population of electrons with a power-law distribution of Lorentz factors γ e , defined between γ min and γ Max : where N 0 is the normalization of the electron distribution in units of cm −3 , and n the electron spectral index. The collective synchrotron emissivity (in units of erg s −1 cm −3 Hz −1 sr −1 ) from this population of electrons is then computed by integrating the synchrotron power over γ e and ϕ as 3 For further details on this and other radiative mechanisms see [53][54][55][56].
Given the emissivity, the νF ν flux in the observer's frame (in units of erg cm −2 s −1 ) is computed as where d L is the luminosity-distance of the source. It can be shown by performing analytically the integrals that if the electron distribution is a power-law distribution with index n as in Equation (2), then the differential photon flux dN/dE (which is a quantity commonly used in high-energy astrophysics, equal to F E (E)/E) is also a power-law distribution with index Γ dN/dE = (n + 1)/2, the flux density F ν has index α = (n − 1)/2, and the energy flux νF ν (ν) has index p = (n − 3)/2. Given that ν = δν /(1 + z), and that j ν ∝ ν −α , the energy flux is proportional to δ 3+α .
At low energies, an important physical process is synchrotron self-absorption: the photons that the electrons radiate can be re-absorbed by the electrons themselves. Following Einstein's approach, it is possible to calculate the self-absorption coefficient (in units of cm −1 ) as In the presence of this absorption term, the intensity of the synchrotron radiation (in units of erg s −1 cm −2 Hz −1 sr −1 ) from this spherical emitting region in the jet can be computed as [57] I syn,ν (ν ) = j syn,ν (ν ) where τ(ν ) = 2Rµ syn (ν ). The effect of self-absorption on the radiated synchrotron spectrum is a change in the power-law index: at frequencies below ν sel f −abs the radiated flux density is proportional to ν 5/2 . Said otherwise, at frequencies below ν sel f −abs the emitting region is optically thick to synchrotron radiation. Self-absorption in the emitting region has one important consequence: a single emitting region in the relativistic jet cannot reproduce blazar SEDs in the radio band. The emission in the radio band must come from other regions which are necessarily optically thin to synchrotron radiation, as will be discussed in Section 2.1.2.

Self-Consistent Electron Distribution
As electrons radiate synchrotron photons, they lose energy, or, using the jargon of astrophysicists, cool. When calculating the synchrotron spectrum from a stationary electron population in the jet, it is thus important to first compute the self-consistent electron energy distribution at equilibrium. The differential equation that describes this process is where Q e (γ e ) is the function describing continuous (time-independent) injection of electrons, τ syn (γ ) is the synchrotron cooling timescale, and τ ad is the adiabatic timescale. The latter is the timescale associated with the expansion of the emitting region as it travels along the jet, and it is energy-independent. It is usually assumed to be τ ad R/c. From a mathematical point of view, the term in τ ad can also be interpreted as an energy-independent particle escape from the emitting region. The synchrotron cooling timescale is defined as τ syn (γ e ) = 1 γ e ∂γ e ∂t −1 and it is equal to where u B = B 2 /8π is the magnetic energy density. The particle distribution at equilibrium can be computed putting Equation (7) equal to zero. The integral solution is [58] N e (γ e ) = e −γ e,break /γ e γ e,break τ ad dζQ e (ζ)e γ e,break /ζ (9) where γ e,break = 3mc 2 4σ T u B R is the value of γ e that satisfies τ syn (γ e ) = τ ad . It can be shown that if the injection term is Q e (γ e ) = Q 0 γ −q e , the effect of synchrotron cooling on the particle distribution is a change in the spectral index, which softens and becomes equal to q + 1 for particle Lorentz factors above γ e,break . From the point of view of the energy flux distribution, it also softens above ν e,break , and becomes equal to p + 1/2. This is true only if q ≥ 2. If q < 2, the particle distribution shows a pile-up at γ e,break [59,60]. Another specific case happens when the injection function is also characterized by a minimum Lorentz factor γ min,inj . If γ e,break > γ min,inj (slow-cooling regime), then the solution is the same as before. If however γ e,break < γ min,inj the steady-state solution is completely cooled (fast-cooling regime), and characterized by an index q + 1 above γ min,inj , and an index of 2 between γ e,break and γ min,inj [61].
We are now in the position to quantitatively compute a self-consistent synchrotron spectrum from an electron distribution in a spherical emitting region in the relativistic jet. The SED can be parametrized by a broken-power-law distribution, proportional to ν 7/2 below the self-absorption frequency ν sel f −abs , to ν −p between ν sel f −abs and ν break , and to ν −(p+1/2) above ν break .
It would be surprising if such a simple scenario would be able to reproduce blazar SEDs, and indeed, it does not. Electron synchrotron radiation from a single emitting region in the jet can reproduce the first blazar SED component, at one price: while the emission is consistent with a broken-power-law distribution, the electron indexes n 1,2 below and above γ e,break are not consistent with a simple break of 1 as calculated above. What we usually observe are much stronger breaks, with n 2 − n 1 > 1, with n 1 = 1.5 − 2.5 and n 2 = 3.5 − 4.5. For this reason, the most common approach is to consider n 1,2 and γ e,break as free independent model parameters [62,63], and discuss a posteriori their consistency with a synchrotron cooling break. In a minority of cases, observations are consistent with a broken power-law with a spectral break of one [64][65][66][67]. Alternatively, it is possible to assume that the electron distribution at equilibrium is in the fast-cooling regime. In this case, if n 1 is consistent with 2.0, it is possible to fit the SED assuming an injection index q 2.5 − 3.5 [68,69]. This approach has however two drawbacks: the first one is that it needs an extra free parameter η that is the scaling factor of the adiabatic timescale that has to be increased to be in the fast-cooling regime τ ad = ηR/c; the second one is that an injection index q 3.0 is significantly softer than standard acceleration mechanisms. The reason we observe stronger breaks may be that one of the assumptions we made are not correct: there could be multiple acceleration/emitting regions [70], non-linear cooling terms (relevant when adding the inverse-Compton process, see Section 2.2), a non-homogeneous magnetic field, incomplete cooling (i.e., the SED is from an electron population that is not at equilibrium) [71,72]. Alternatively, the injected particle distribution may not follow a simple power-law function, and could show intrinsic curvature in the form of a log-parabolic function [73,74], or a power-law one with exponential cut-off at γ Max . The observational measurements of n 1 and n 2 may thus be the asymptotic power-law fits over a narrow energy band.

Low-Energy SED Modeling
As discussed above, a single emitting region in the jet radiating synchrotron photons cannot explain the observed spectrum from radio up to optical-X-rays, due to self-absorption in the radio band. Blazar SED modeling is thus, inevitably, at least two-zones: a small, compact emitting region is responsible for the high-energy part of the first SED component (optical frequencies and above), the one that shows rapid variability; a larger, less dense emitting region is where the radio photons are produced. It is natural to consider that this radio emission comes from the extended radio-jet. The exact geometry is not clear: the radio emission can be produced downstream with respect to the optical-X-ray one, or be associated with a second flow that encompasses the inner jet [75], or the whole emission could be produced not in a small overdensity in the jet, but in a conical jet, with varying magnetic field [76,77]. For most blazars, a two-zone model is enough to fit the first SED component [78], although the transition between the two components remains unclear, with potential contribution of the extended region up to the optical/UV band [79,80].
A fit of the model to the observations can be used to put constraints on the physics of the emitting region, or of the acceleration mechanisms. The number of free parameters of the synchrotron model is nine: three for the emitting region (δ, B, and R), and six for the particle distribution (N 0 , γ min,break,Max , and n 1,2 ). If we force a self-consistent electron distribution, γ break , and n 2 can be expressed as a function of other parameters, and the number of free parameters becomes seven; if we further assume that γ min,Max are respectively low-enough ( 1) and high enough (≥10 6 ) that they cannot be constrained by observations, we are left with five free parameters. Still, the true number of observables is much lower: the synchrotron spectrum is fully characterized by its peak frequency and luminosity, and by its index, i.e., we have only three observables. This is true both for the small plasmoid in the jet, and for the large component emitting in the radio band. We could constrain the model further if we could see the self-absorption frequency, but the one for the plasmoid is hidden below the extended-jet emission. Synchrotron modeling of the low-energy SED of blazars is, by itself, degenerate , which means that although it is possible to find a set of model parameters that can reproduce the SED, we cannot constrain the parameter values.

Time-Dependent Modeling
Until now we focused on the synchrotron emission by a stationary distribution of electrons in the jet, comparing theoretical expectations with blazar SEDs, that can be seen as snapshots of the blazar emission at a given instant. One of the key observing properties of blazars is however variability, and key pieces of physical information can be extracted by studying blazar flares. It is possible to simulate blazar light-curves by solving Equation (7) as a function of time, following the evolution of the electron distribution with time, and the associated synchrotron radiation.
As a first approximation, a synchrotron blazar flare can be modeled with an instantaneous injection of particles Q e (t ) = Q 0 δ(t − t 0 ). In this case the time-evolution is fully determined by the cooling timescale τ syn for Lorentz factors higher than γ break . The measure of τ var can thus provide a direct measurement of τ syn = τ syn (B , γ e ) and thus help constrain the parameter space.
Increasing complexity, all plausible acceleration processes are energy-dependent. In a more realistic scenario we can thus define an acceleration timescale τ acc = τ acc (γ e ). The flare is thus expected to be characterized by an asymmetric profile, due to the different values of τ acc and τ syn . The energy dependency of both acceleration and cooling implies that the spectrum changes during the flare. In a plot Γ dN/dE vs F ν this translates into a hysteresis cycle. Equivalently, it will appear as a time-lag among light-curves computed at different energies [59,[81][82][83]. Both hysteresis and time-lags have been observed in X-ray flare of blazars [84][85][86][87], and can be used to simultaneously constrain the magnetic field strength and the Doppler factor of the emitting region (with a degeneracy between the two), resulting in B 0.1 for δ 10.
When modeling blazar flares, there is another timescale that needs to be taken into account that is the light-crossing time R/c: the variability will thus be diluted as observed from Earth if R/c > τ syn [88].

Inverse-Compton Emission
Compton scattering is the interaction between a charged particle and a photon. Following Section 2.1 on synchrotron radiation, let us start again with an electron. In the direct Compton scattering a photon scatters over an electron at rest in an atom, and loses energy in the interaction. The photon wavelength after scattering changes as with ζ the scattering angle, and with h m e c = λ C the Compton wavelength. In the inverse-Compton scattering a high-energy electron interacts with a low-energy photon, and transfers part of its energy to it. The inverse-Compton scattering is a very efficient process to produce X-ray and γ-ray photons in astrophysical sources starting from low-energy photon fields. Let us mark with a star the reference frame at rest with the electron, and let's mark with the 1, 2 subscripts the values before and after scattering, respectively. The photon a-dimensional energies are ε 1,2 = hν 1,2 /m e c 2 . The electron Lorentz factors are γ e;1,2 = E e;1,2 /m e c 2 , with γ e;1 = 1. It can be shown with simple kinematic and energy/momentum conservation equations that in the reference frame at rest with the electron where ξ is the angle between the incoming and scattered directions of the photon. When ε 1 1, Equation (11) simplifies as ε 2 ε 1 , which is the photon energy in the electron rest frame is unchanged. The condition ε 1 1 represents the Thomson regime, where the photon energy is much lower than the electron rest mass. The equations describing the change of reference frame from and to the observer (unstarred) are: where ϑ 1 is defined as the angle between the electron and the incoming photon, in the observer's frame, and ϑ 2 the one of the scattered photon, in the electron's rest frame. In the Thomson regime, we thus have In the case of a relativistic electron with γ e 1 and β e 1, the photon energy in the observer's frame after scattering can be as high as ε 2 4ε 1 γ 2 e . This occurs in the so-called head-on approximation, which is for cos ϑ 1 = −1 and cos ϑ 2 = 1.
The expression of the inverse-Compton cross-section is equal to σ T in the Thomson regime, while for the general case it has been calculated by Klein and Nishina [89], and is known as the Klein-Nishina formula. The differential cross-section, in the reference frame of the electron, is where dΩ = dϕ d cos ξ . The integral cross-section is The cross-section drops steadily above ε 1 1, and asymptotically behaves as log(2ε 1 ) + 1/2 . The regime ε 1 1 is known as the Klein-Nishina regime. Solving Equations (11) and (13) for the Klein-Nishina case, it is easy to show that in the best case ε 2 γ e .
Let's discuss now the inverse-Compton emission from a distribution of electrons N e (γ e ) scattering a distribution of photons n ph (ε ). The computation of the inverse-Compton emissivity is analytically complex, requiring transformation of the Klein-Nishina cross-section from the electron rest frame to the jet frame. It has been investigated under several approximations (e.g., Thomson and Klein-Nishina approximation, head-on collisions, monochromatic photons). It is best solved numerically performing integrals over the electron and photon distributions times the cross-section (Equation (15)). For an isotropic distribution of electrons and photons, the inverse-Compton emissivity can be expressed as: where C(ε , γ e , ν ) is the Jones' Compton kernel [90].
As for the synchrotron radiation, the fact that part of the electron energy is transferred to the photons imply that electrons are cooling in the process. It is thus important to include the inverse-Compton cooling losses when computing the self-consistent electron distribution at equilibrium (see Section 2.1.1). Let's start first with the Thomson regime: in this case the electron losses due to inverse-Compton scattering (for an isotropic photon field) take a form which is very similar to the synchrotron one, and the associated timescale is where u ph is the energy density of the scattered photon field. When computing the electron distribution at equilibrium, it is just possible to simply use a single cooling timescale τ cooling (γ e ) = 3mc 4σ T (u B +u ph ) 1 γ e , and the effect on the electron distribution is again a broken-power-law with indexes n and n + 1 below and above γ break respectively (in the slow cooling regime).
In the Klein-Nishina regime, the inverse-Compton losses are also suppressed: if the inverse-Compton losses dominate over the synchrotron ones for high Lorentz factors (as could occur in the presence of bright external fields), the result on the steady-state electron distribution is a hardening at the highest energies [91].

Synchrotron-Self-Compton
Going back to the spheroidal plasmoid travelling with Lorentz factor Γ in the SMBH relativistic jet, filled with both an homogeneous magnetic field B and an electron population N e (γ e ), it is natural that the synchrotron photons will suffer inverse-Compton scattering over the same electrons population that produced them. This kind of emission is called synchrotron-self-Compton (SSC), and it is a very efficient process to transfer energy from particles to photons [92].
In the Thomson regime, it can be shown analytically that if the electron distribution N e (γ e ) is a power-law with index n, not only the synchrotron energy flux has index p = (n − 3)/2, but the SSC energy flux has also index p = (n − 3)/2. The SSC emission is thus mirroring the synchrotron radiation at higher energies. It is thus natural to investigate SSC emission as the radiative process responsible for the high-energy SED component in blazars. Now let us investigate two key aspects: how high are these SSC energies, and what is the relative luminosity L SSC /L syn [93]: • Let's assume that the synchrotron peak frequency ν syn,peak is directly related to γ break , i.e., let's assume that n 1 < 3 and n 2 > 3. In this case The SSC peak frequency in the Thomson regime can be easily computed assuming that the electrons at γ break are the ones primarily scattering the photons at ν syn,peak , and thus ν SSC,peak = 4/3 ν syn,peak γ 2 break . Assuming typical values for HBLs, i.e., a ν syn,peak 0.1 keV, in the soft X-rays, and γ break 10 4 (see next Section), ν SSC,peak can thus reach 10 GeV. The ratio ν SSC,peak /ν syn,peak can be used to put a strong constraint on the model parameters, given that • the total luminosities of the synchrotron and SSC components L syn,SSC can be expressed, in a first approximation, as L syn,SSC ν syn,SSC;peak L syn,SSC (ν syn,SSC;peak ), which is their peak SED luminosities. The ratio of the two SED peak fluxes depends only on the magnetic energy density and the synchrotron photons energy density u ph,syn : ν SSC,peak F SSC (ν SSC,peak ) ν syn,peak F syn (ν syn,peak ) = ν SSC,peak L SSC (ν SSC,peak ) ν syn,peak L syn (ν syn,peak ) In the Klein-Nishina regime on the other hand, the SSC emission is suppressed, producing a major correction to the shape of the SSC component at higher energies. As we increase the peak of the synchrotron emission, the SSC peak frequency increases quadratically until the scattering enters into the Klein-Nishina regime, and higher-energy photons do not contribute any more to the emission. Even for the Klein-Nishina regime, simple relations between the frequencies and luminosities of the synchrotron and SSC components can be computed [93]. In is important to underline that the computation of best-fit solutions for the SSC model relies on a good estimate of the peak frequencies and fluxes. Given the wide variety of peak frequencies in the blazar population, these may fall in relatively unexplored energy bands.
When discussing the modeling of the low-energy SED component of blazars as synchrotron radiation (see Section 2.1.2), we showed that the modeling of the synchrotron component alone is degenerate. This is not the case when we add the information from the SSC radiative component. We recall that the number of free parameters for a spherical plasmoid in the jet, filled with an electron population parametrized as broken-power-law function, is nine. Assuming that γ min,Max are respectively low and high enough that they cannot be constrained by the data, we are left with seven model parameters. With the information from both SED components, the number of observables is six: the frequencies and luminosities of the synchrotron and SSC peaks, and the two indexes of the synchrotron distribution before and after the peak. We are left with just one missing constraint, which can be the self-consistency of the electron distribution, or can be the observed fastest variability timescale τ var . It can be related to the size of the emitting region R and the Doppler factor δ via a causality argument A one-zone SSC model can describe well the SEDs of HBLs (assuming that the radio emission comes from a larger, less dense region in the jet). An example of an SSC modeling of the HBL Markarian 421 is shown in the top-left panel of Figure 2. In the recent years, a number of algorithms have been published to fully constrain the SSC model parameter space, and have been applied successfully to HBLs [63,[94][95][96][97]. Typical best-fit solutions are: where L e = 2πR 2 cΓ 2 u e is the electron luminosity in the observer's frame, and u e = m e c 2 dγ e γ e N(γ e ) is the electron energy density. The SSC scenario fails to reproduce the SED of other blazar sub-classes: it cannot reproduce the SEDs of LBLs, nor FSRQs; it can reproduce the SEDs of EHBLs at the expenses of a high Doppler factor (δ ≥ 50) and of a high value of γ min ≥ 10 3 [36,98].

Synchrotron-Self-Compton: Time Signatures
Now that we have a model which is able to describe the blazar SEDs (at least for HBLs), we can add the time information to further test it and constrain the parameter values. As for the synchrotron case, we can compute the evolution of the SED with time by solving Equation (7) (including inverse-Compton losses) and calculating both synchrotron and SSC emission as a function of time. One difficulty of SSC time-dependent modeling is that the inverse-Compton cooling term depends on the density of the synchrotron photons, which depends on the integral of the electron distribution itself as a function of time (which is what we want to solve). In presence of significant SSC losses the differential equation becomes thus non-linear and the solution is significantly more complex [99,100]. Blazar light-curves in the γ-ray band have been successfully fit with the SSC model [79,101]. Similarly to synchrotron radiation, time-lags and hysteresis are also expected in SSC light-curves due to the energy dependency of both the acceleration and cooling terms [102]. But contrarily to what is observed in X-rays, there has not been any significant detection of time-lags or of hysteresis in the γ-ray band [103,104]. This negative result is likely due to the intrinsic lower sensitivity of gamma-ray telescopes compared to X-ray ones, and not to a failure of the model.  Another important prediction of the SSC scenario is the correlation between the variability of the two SED components, which are expected to vary in concert. If the variability is simply due to an achromatic change in the particle density in the emitting region, the SSC component is expected to vary as the square of the synchrotron component (due to both the change in the particle distribution, and in the synchrotron photons). A correlation plot between X-rays and γ-rays is thus expected to show a simple quadratic relation. If on the other hand the variability is energy-dependent, the expected correlations are more complex, and depend on the energy bands of the observations, the SED peak frequencies, the energy-dependency of the injection [105]. Correlation studies among different wavebands are indeed a powerful tool to test radiative models when the statistics in an individual band does not allow light-curve fitting [69,[106][107][108][109][110].

External Inverse-Compton
Synchrotron photons are not the only low-energy photons that can be up-scattered to γ-rays. The SMBH environment is very bright, and there is a huge variety of photon fields that can go through inverse-Compton scattering with electrons in the jet. This radiative process is usually called External-inverse-Compton (EIC) and depends strongly on the properties of the external field and on the exact location of the emitting region r (measured from the SMBH). The computation of the EIC component is similar to the SSC one, but requires as first step the transformation of the external photon field distribution into the jet's frame. For blazar modeling, the relevant photon fields, starting from the AGN core outwards, are: • the SMBH accretion disk is a reservoir of thermal photons. EIC emission over disk photons is however hindered due to the Doppler deboosting [111]: the plasma in the jet is travelling away from it with Lorentz factor Γ, and thus the disk photon field energy density, in the reference frame of the plasmoid, is strongly suppressed by a factor Γ −3 . EIC scattering over disk photons can thus dominate the overall γ-ray emission only if the emitting region is located close to the disk itself. In this case the plasma can see part of the disk at angles large enough to reduce the de-boosting. • the BLR produces bright emission lines (Lyα being the dominant one) that can serve as target photon field for EIC emission [112]. In this case the emission is boosted or deboosted in the reference frame of the emitting region as a function of r and r BLR . It is important to underline that the BLR is structured, and that different emission lines are produced at different distances r BLR,i from the SMBH. The size of the BLR can be expressed as a function of its luminosity (or other proxies such as the luminosity of the continuum at 5100 Å, or the disk luminosity [6,113]). The EIC emission is thus a superposition of several components coming from the different lines [114], with their relative strength depending on r.
A second radiative component coming from the BLR is the photon field from the accretion disk which can be Thomson-scattered by electrons in the BLR and thus seen boosted in the reference frame of the blob. This component depends also on the optical depth τ BLR of the BLR. A parameter which is not so well known is the aperture angle α BLR of the BLR: it can have a spherical geometry (α BLR = 0) or be flattened over the disk (α BLR → π/2) [115,116]. The assumed BLR geometry has a direct impact on the photon field seen in the blob reference frame, and in the limit of a flat BLR the computation of the associated EIC component becomes similar to the one from the accretion disk. • the emission from the dusty torus is thermal and can be well described with a single-temperature black-body distribution with temperature T 1000 K . The torus location can also be parametrized as a function of the disk luminosity r torus = 2.5 × 10 18 L disk /10 45 erg s −1 cm [117]. Similarly to the BLR aperture angle, the aperture angle of the torus is also not well known. It is commonly assumed to be of around π/4. As with the BLR photon field, the blob of plasma sees a constant photon density from the torus for r r torus , and then a rapidly decreasing photon field while crossing r torus and leaving the torus behind [118]. An example of an EIC modeling of the FSRQ 3C279 is shown in Figure 3. • A ubiquitous photon field is the one from the Cosmic Microwave Background (CMB): electrons in the jet can also scatter these photons to higher energies. This EIC emission has the important properties of being the only one dependent on the redshift of the source. EIC over the CMB is not so efficient when the emitting region is located close to the SMBH, because there are always stronger photon fields. But it can become an important radiative process to explain X-ray and γ-ray emission from SMBH jets at larger scales ( [119], but see [120]). lSince the CMB intensity scales as (1 + z) 4 , this emission process can play a more important role for high-redshift quasars [121]. • as we discussed before, structured jets can explain the synchrotron emission from radio up to optical-X-rays. The emission from the other regions of the jet can thus serve as target photon field for EIC scattering. In configurations such as the spine-sheath scenario (a fast inner jet surrounded by a slower layer), the emission from the external jet can be boosted in the γ-ray emitting region due to the relative motion of the two components [122][123][124].

Energy Budget of the Emitting Region
An important characteristic of all radiative models is its energy budget: how much energy is required to reproduce the observations? and how is energy distributed among the various components? The total jet power (in erg s −1 ) can be expressed as an explicit function of the various energy densities as L jet = 2πR 2 cΓ 2 (u B + u e + u ph ) (24) where the factor of 2 takes into account that there are always two jets. Equipartition between the various radiative components is always attractive especially in physical systems at equilibrium, because it directly provides a minimum power solution. It is thus interesting to note that SSC solutions for HBLs are not near equipartition, the emitting region being always particle dominated [125,126]. On the other hand, EIC solutions for FSRQs show near-equipartition conditions between particles, magnetic field, and external photons [127].

Electron-Positron Pair Production
The last leptonic process that is relevant to calculate the radiative output of SMBH jets is pair-production via photon-photon interaction [128] This process is the inverse of pair annihilation, in which one electron and one positron annihilate into 0.511 MeV photons. Pair-production can occur only if the energy of the two photons is higher than the rest mass of the leptons that are produced, which is ε 1 ε 2 = hν 1 hν 2 /m 2 e c 4 > 1. The cross-section of the interaction is [129] where s = ε 1 ε 2 (1 − cos ϑ)/2, and ϑ the interaction angle. The cross-section can be reasonably well approximated as This approximation is known as the δ-approximation. The pair-production process is more efficient when the energy of the two photons are inversely proportional. The two photons disappear in the interaction, and while for low-energy photon fields the disappearance of a single photon is negligible, for the high-energy photon field a single photon can contribute significantly to the overall energy flux. For γ-rays, pair-production can result in significant absorption, which can be expressed via the following absorption coefficient (similarly to Equation (5); note that the following equation is correct only if an angle-averaged cross-section is used, otherwise an additional integration over the angle is required): where n ph (ε ) is the number density of the low-energy photon field, which can be computed from the intensity The most relevant γ-γ absorption term for γ-ray blazars is external to the source, and it is represented by the integrated emission from all galaxies (stars and dust) that permeates the Universe, the Extragalactic Background Light (EBL). The EBL acts as an absorber for TeV photons in the Universe, and this absorption effect increases with the distance from Earth and with the energy of the γ-ray photon [130]. The result is a softening of the spectrum of TeV blazars, and an horizon for TeV astronomy. The TeV horizon is defined as the energy-dependent redshift for which the opacity τ EBL = 1. Although the EBL significantly hinders TeV observations of distant blazars (the most distant source observed in the TeV band with Cherenkov telescopes is at z = 0.944, although detected only up to 0.175 TeV [131]), this absorption effect can be used to measure the EBL itself, opening up the path to γ-ray cosmology [132]. The pairs produced in the interaction with the EBL are not lost, and can produce further emission along the line of sight mainly via inverse-Compton scattering over CMB photons, although energy losses from plasma instabilities are also investigated [133]. This emission can be detected both as a new radiative component in the GeV band [134], or as extended gamma-ray pair halos [135], both modulated by the strength of the magnetic field seen by the pairs. None of these features has been detected as per today [136][137][138], and these observations can be used to put constraints on the inter-galactic magnetic field [139,140].
Going back to the inverse-Compton scattering in SMBH jets, all photon fields which are inverse-Compton scattered to γ-ray energies also act as an absorber to the same γ-ray photons via pair-production. Let us start by looking at the SSC case. Under the δ-function approximation for σ γ-γ , and by forcing an equality to the causality relation (Equation (22)), it is possible to write τ γ-γ < 1, where the pair-production opacity depends only on the Doppler factor δ SSC (and physical observables) [141], and it is thus possible to provide a lower limit on δ SSC where ν 0 = 1.6 × 10 40 /ν γ Hz is the frequency of the low-energy photons that act as primary absorber for photons at ν γ , and α is the index of the flux density F ν at ν 0 . For typical blazars parameters, Equation (30) results in δ SSC ≥ 5 − 10. The detection of TeV photons from HBLs can thus be used to put a strong limit on δ, constraining further the parameter space of the SSC model.

On the Location of the γ-ray Emitting Region in FSRQs
For the EIC case, the opacity to γ-γ pair production depends on the location r of the emitting region with respect to the external fields. The opacity is expressed as the integral from r outwards of the absorption coefficient: The presence of the BLR is very relevant for opacity calculations: if r < r BLR , radiation from the BLR will significantly absorb γ-ray photons. Assuming that the most relevant emission line is Lyα (13.6 eV), this will cause an absorption feature at 100 GeV. Detailed calculations including all relevant BLR lines and the thermal continuum from the disk scattered by the BLR [114,[142][143][144][145] show that the very detection of γ-ray photons with energies greater than 100 GeV from FSRQs is enough to put the location of the emitting region beyond the BLR. Observations of FSRQs at E > 100 GeV confirmed this scenario and proved that, at least for the sources that have been detected at high energies, the emitting region is at or beyond r BLR [146][147][148][149][150][151][152].
The dusty torus also leaves an absorption imprint on the γ-ray spectrum. Given its typical temperature of 1000 K, the pair-production absorption on the torus photons is expected to be seen in the TeV band in γ-ray FSRQs. Due to the simultaneous absorption on the EBL, and to the fact that FSRQs have an SED peak in the 100-MeV range with very few TeV photons, such a cut-off has never been observed as per today [153].

Pair Injection
The electron-positron pairs that have such an important role for γ-ray absorption in SMBH jet are injected in the emitting region, and can be treated as a secondary injection of leptons. They will behave exactly as the primary ones, and will radiate synchrotron and inverse-Compton photons, reaching their own equilibrium distribution. The big advantage is that while we do not know the details of the primary accelerator and we simply parametrized it, we have the information on the energy distribution of the pairs injected via pair-production [154]: where n ph (ε 1 ) and n ph (ε 2 ) are the photon densities of the two photon fields. The injection term can then be used to compute the steady-state distribution and the associated synchrotron and inverse-Compton emission following Equation (7). For single-zone blazar leptonic models, emission by this secondary leptonic population is usually negligible, hidden by the primary one.

Proton-Synchrotron Emission
Let us now consider the case of a jet in which there are also protons (and possibly higher-Z nuclei) together with leptons. As with electrons, the proton energy distribution is defined between γ p,min and γ p,Max as a power-law function: Going back to our spherical plasmoid in the jet, the first radiative process that needs to be studied is synchrotron radiation by protons.With respect to electrons, the emission from protons is suppressed and needs to be compensated by increasing either the magnetic field value, the particle density, or both. It is interesting to test the possibility that proton-synchrotron radiation is responsible for the high-energy SED component in blazars. Indeed, proton-synchrotron radiation would have the same index α p = (n p − 1)/2 as SSC radiation, assuming that the injection index of protons is similar to the one of electrons. Following Equation (19) (replacing the electron mass with the proton mass) it is possible to have a proton-synchrotron peak frequency at around 10-100 GeV, consistent with HBL SEDs, assuming γ p,Max 10 9 , B 10-100 G, and δ = 10-50. The proton-synchrotron solution for HBLs lies thus in a completely different part of the parameter space with respect to the SSC solution, with a much stronger (a thousand times higher) magnetic field [155][156][157]. The number of free parameters of the proton-synchrotron scenario is much larger than the SSC one, due to the fact that the proton distribution adds 5 additional free parameters (assuming an additional γ p,break due to cooling, and the two indexes n p,1/2 ), for a total of 14 model parameters. Such a model is clearly degenerate, although the number of free parameters can be considerably reduced by assuming some physically motivated constraints, such as: self-consistency of break Lorentz factors, co-acceleration of electrons and protons (and thus equality of n e,1 and n p,1 ), constraint on γ p,Max via equation of acceleration and cooling, constraint on R via causality arguments. But even imposing all these constraints the model remain degenerate, and and we can only constrain allowed regions in the parameter space [98].

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

Photo-Meson Production
The photo-meson (or photo-pion) process is the production of pions in proton-photon interactions: The photo-meson cross-section is dominated by resonances (∆ and N baryons), direct pion production and multiple pion production. For astrophysical applications, the relevant information is the production of pions, both neutral and charged. In order to produce pions, the photon energy in the proton frame has to be higher than m π c 2 + (m π /2m p )c 2 145 MeV. In the emitting-region frame, (145/γ p ) MeV, meaning that protons with Lorentz factor γ p ≥ 10 7 can photo-meson-produce over ultra-violet photons.
Pions then decay primarily (∼ 100% branching ratio) as π 0 → 2γ π + → µ + + ν µ → e + + ν e +ν µ + ν µ π − → µ − +ν µ → e − +ν e + ν µ +ν µ (35) Photo-meson production has thus one key property: neutrinos are produced together with photons, and they can escape the emitting region without suffering from any additional absorption or energy losses. Detection of neutrinos from an AGN is thus a smoking-gun signal for the presence of relativistic protons in the jet. This implies that AGNs can accelerate protons to high energies, and can be responsible for the flux of particles that we call, once they reach Earth, cosmic rays. Joint neutrino and photon observations of AGNs can be used to constrain three key parameters of the model: how many protons are accelerated, what is their energy distribution, and what is their maximum energy.
Together with the neutrinos, photo-meson production injects photons and electron/positron pairs in the emitting region. The spectrum of the injected pions has a maximum energy equal to the proton maximum energy but reduced by the inelasticity κ 0.2. The maximum energy of photons from the π 0 decay can be as high as 0.1γ p . For protons with γ p = 10 7 , it translates into photons with E 10 PeV. These photons do not reach us: they are absorbed via pair-production both in the jet, and in the path to Earth. The absorption occurring in the source produces an electron/positron pair that starts radiating (and cooling). Given the high energy of the first photon, the synchrotron radiation by these secondary leptons is at energies high enough that the synchrotron photons can suffer a second pair-production, and so on. Photons from the π 0 decay can thus trigger a synchrotron-supported pair-cascade in the emitting region (or, if the soft photon field energy density is larger than the magnetic one, an inverse-Compton supported pair-cascade). Emission from the cascade can be computed calculating the pair distribution at equilibrium and their synchrotron radiation. When τ γ-γ 1 the cascade is saturated, and the resulting synchrotron emission is a power-law with index α = 1, or in terms of energy flux p = 0: the energy transferred from protons to pions and then to high-energy photons, is then redistributed to low energies.
Electrons and positrons produced in the decay of charged pions go through a process similar to the photons from π 0 : they have extremely high-energy, and their synchrotron radiation is also in the PeV band. They thus trigger as well a pair-cascade, which behaves exactly like the one from π 0 .
Another important emission process from photo-meson interactions is the synchrotron radiation by µ ± . Although not stable leptons, muons in the jet can survive long enough to radiate synchrotron photons before decaying to electrons/positrons. To compute the muon energy distribution at equilibrium, the muon decay term must be added to Equation (7) as follows: where τ dec = 2.2 × 10 −6 seconds. The corresponding integral solution for the muon distribution at equilibrium is: Synchrotron radiation by muons emerges as a third radiative component, at higher energies compared to the synchrotron radiation by the parent protons. In some parts of the parameter space it can represent the dominant radiative emission in the TeV band. Another important consequence is that every time the muon synchrotron radiation is important, it means that muons can cool efficiently before decaying into electrons/positrons, and this energy loss has to be included when calculating their injection in the emitting region.
An interesting property of photo-meson interactions is the creation of neutrons, which can escape the emitting region without interacting with magnetic fields, nor producing Bethe-Heitler pairs (see Section 3.2.2). These neutrons can transfer a significant amount of energy at much larger distances downstream the jet [158,159]. In presence of dense photon fields they can trigger additional photo-meson production, or they can decay into protons (with a life-time of γ n × 880 seconds) and radiate again synchrotron photons in the presence of magnetic fields. This "neutral beam" model has drawn attention for its possibility to naturally produce two separate emitting regions in the jet, physically separated but casually connected [160,161].
The main problem we are facing with photo-meson production is that the exact computation of the energy distribution of all secondary particles injected in the emitting region is a complex numerical task. The best approach is to compute them via Monte-Carlo simulations [162]. Several approaches to fit the results of Monte-Carlo simulations and provide a parametrization of the injected secondary particles over a wide parameter space are particularly useful [163,164].

Bethe-Heitler Pair-Production
Bethe-Heitler pair-production is the following process It is a process which is in competition with the photo-meson production, although it happens at lower energies. The threshold for this interaction is lower than the photo-meson one by a factor m e /m π 0.004. In terms of proton Lorentz factors, protons with γ p ≥ 10 5 can pair-produce over 10 eV photons.
For low-energy protons, the Bethe-Heitler pair production is the main proton-photon interaction, but as soon as the proton energy becomes greater than the energy threshold for photo-meson production, the latter takes over. As expected, pairs injected via this process have a lower energy compared to pairs from the photo-meson [165]. The energy distribution of the injected Bethe-Heitler pairs can be computed analytically [163,166,167], or using a Monte-Carlo approach [168,169]. Their fate is the same as all other leptons in the emitting region: they radiate synchrotron and inverse-Compton photons, and if their energy is high enough (and if the low-energy photon field is dense enough) they can trigger a pair-cascade.

Hadronic and Lepto-Hadronic Models
Blazar hadronic models are a large family of models which share the presence of a parent population of protons that are accelerated in the SMBH jet. But as discussed above, from the electromagnetic point-of-view the only true hadronic radiative process is proton-synchrotron emission. Solutions in which the high-energy emission is produced by secondary leptons produced in p-γ interactions are usually referred to as lepto-hadronic models.
Pure hadronic (proton-synchrotron) models can fit blazar SEDs [155][156][157]170], but they face a major problem: the proton density required to fit the data is rather high, and when expressed in terms of L p it gets close or clearly much higher than the Eddington luminosity of the SMBH [117]. Although L Edd should be treated as order-of-magnitude estimate for the available accretion power, if the power in the jet required to fit the data is super-Eddington by several orders of magnitudes it means that either the relation between accretion and jet needs to be revisited, or, more likely, that this particular solution is not viable. This energy-crisis for hadronic models is particularly true for FSRQs [171]. For low-luminosity HBLs and EHBLs it is not the case, and hadronic solutions with L jet < L Edd can be found [98,172]. Even if the high-energy SED component is dominated by proton-synchrotron radiation, emission by secondary leptons from p-γ interactions can emerge (namely for compact and dense solutions) in the X-rays and in the TeV band [173][174][175]. In terms of energy budget of the emitting region, proton-synchrotron solutions are out of equipartition, with u B ≥ u p . An example of a hadronic modeling of the HBL Mrk 421 is shown in the second plot of Figure 2.
In the lepto-hadronic solutions, the proton-synchrotron radiative component is suppressed, and this can be easily achieved by reducing the magnetic field to B ≤ 1 G. We are back in the part of the parameter space where the SSC emission dominates and indeed the simplest lepto-hadronic model is an SSC solution loaded with relativistic protons that produce secondary leptons via p-γ interactions over the electron synchrotron photon field (this solution is called one-zone lepto-hadronic model, to highlight the fact that all radiative mechanisms are coming from a single emitting region in the jet, and there external photon fields are negligible). The emission by secondary leptons emerges again in the X-rays (as Bethe-Heitler component) and in the TeV band (as photo-meson component) [172]. This kind of lepto-hadronic solutions are typically more promising in terms of ν output than pure hadronic solutions. As with hadronic solutions, lepto-hadronic solutions also face energetic issues: especially if the goal of the modeling is to maximize the neutrino output, the required jet power can quickly become very high. A possible way-out to this issue is to use as target photon field not the electron synchrotron radiation, but external photons. By increasing the density of the target photon field, it is possible to achieve the same p-γ output lowering the required proton power. The discovery of the blazar TXS0506 + 056 as first neutrino blazar candidate [176], indicates that a lepto-hadronic solution over an external field may be the favored scenario [177][178][179][180][181][182]. In terms of energy budget of the emitting region, mixed lepto-hadronic are usually particle dominated, although the result depends on the exact shape of the proton distribution (a soft/hard distribution has a significant impact on u p , but a marginal effect on the photon and neutrino emission). An example of a mixed lepto-hadronic modeling of the HBL Mrk 421 is shown in the third plot of Figure 2.
Hadronic radiative models are of particular interest any time leptonic ones face difficulties. One of the most interesting cases is the unusual SED of the nearby radio-galaxy Centaurus A, which shows a unique third radiative component emerging at 100 GeV [183]. This additional component can be naturally explained in a hadronic or lepto-hadronic scenario [184][185][186].
Among blazar hadronic radiative models, it is worth mentioning another potential contribution produced not in the jet, but rather in the path from the source to the observer. If AGNs accelerate cosmic rays, a hadronic beam can be launched emitting photons and neutrinos while travelling towards Earth [187,188]. So far no evidence for this hadronic contribution has been detected. For a recent study of propagation effects comparing both hadronic cascades and electromagnetic cascades (the ones triggered by the absorption over the EBL) see [189].

Hadronic Models: Time Signatures
From the spectral point of view both leptonic and hadronic radiative models provide similarly good fit to current electromagnetic observations, but they can in principle be distinguished from their temporal behavior. Time-dependent hadronic models require the solution of a system of coupled differential equations of the kind: for each species X (protons, photons, neutrinos, leptons), where Q X (t, E) and L X (t, E) represents the injection and loss terms, respectively. From a numerical point of view, hadronic time-dependent codes are significantly more complex than leptonic ones, and have been extensively developed only in the recent years [169,[190][191][192]. As with leptonic models they can be tested on data in two main ways: fitting light-curve profiles, and studying multi-wavelength correlations. For proton-synchrotron models, although the two SED components are produced by two distinct particle populations, it is still possible to reproduce the observed correlations assuming that electrons and protons share the same acceleration mechanism [193]. On the other hand, absence of correlations can also be reproduced and indeed, one of the best applications for blazar hadronic models, are the so-called orphan flare, bright γ-ray flares without multi-wavelength counterparts, which cannot be reproduced in a leptonic one-zone scenario [194][195][196].
Time-dependent models have also the advantage to allow the study of non-linear effects in the development of the hadronic pair-cascades, in which photons from the cascade become targets for pair-production and can result in a unique pulsating signature [197].

Proton-Proton Interactions
The final hadronic process that needs to be taken into account is pion production in proton-proton interactions. As with the photo-meson production, a highly relativistic proton interacts with a target low-energy proton. For astrophysical applications, the only relevant process is the production of neutral pions π 0 and η mesons, which then decay into π ± . The result of the interaction is then where Y and Z are any other particle produced in the interaction. The decay of charged pions is responsible for the associated neutrino emission. As with the photo-meson process, the outcome of proton-proton interactions is better studied via Monte-Carlo simulations, although useful parametrizations have been developed [198,199]. In blazar's jets, proton-proton interactions represent usually a subdominant contribution. In order to have a significant p-p output a very proton-loaded jet is needed [200]. An alternative scenario is represented by interactions of the jet with obstacles, such as BLR clouds, or stars [201][202][203].

Polarization Signatures
An important observable that has the potential to uniquely determine the radiation mechanism, as well as the physical properties of the jet, is polarization. A specific signature of synchrotron radiation is indeed that the emission is linearly polarized. For a power-law electron energy distribution with index n, the maximum polarization fraction is which is equal to 0.69 for n = 2. Detection of polarized radiation in blazars is thus a key prove that what is observed from radio to optical is synchrotron radiation [204]. The polarization angle can be used to map the orientation of the magnetic field in the jet. Recent multi-wavelength campaigns have shown not only that optical polarization is a characteristic of blazars, but that both the polarization fraction and the polarization angle show drastic variability, often correlated with flares at other wavelengths [205]. For an updated review of optical polarization in AGN jets, see [206] in this special issue. Variability in the optical polarization can also be used to put constraints on the geometry of the emitting regions, and in particular on the relative role of the small plasmoid (the single-zone that emits from optical to gamma-rays) and the larger jet [207,208]. Inverse-Compton scattering is also expected to be polarized, depending on the polarization of the target photon field [209]. X-ray and gamma-ray polarimetry have the potential to significantly constrain blazar radiative models: proton-synchrotron radiation is expected to be significantly much more polarized than SSC radiation or EIC radiation (which is no polarized at all if the external field is also not polarized). Leptonic and hadronic radiative processes can thus be distinguished using this unique observable [210,211].

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

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

Abbreviations
The following abbreviations are used in this manuscript: