Radio Galaxies at VHE energies

Radio Galaxies have by now emerged as a new $\gamma$-ray emitting source class on the extragalactic sky. Given their remarkable observed characteristics, such as unusual gamma-ray spectra or ultrafast VHE variability, they represent unique examples to probe into the nature and physics of AGN in general. This review provides a compact summary of their observed characteristics at very high $\gamma$-ray energies (VHE; $>100$ GeV) along with a discussion of their possible physics implications. A particular focus is given to a concise overview of fundamental concepts concerning the origin of variable VHE emission, including recent developments in black hole gap physics.


Introduction
The current decade has seen a tremendous progress in the extragalactic Gamma-Ray Astronomy. Numerous new sources have been discovered by the current generation of instruments, sometimes with highly unexpected and extreme characteristics. More than 2900 of the identified or associated high energy (HE, > 100 MeV) sources in the Fermi-LAT 8-year Point Source List (FL8Y) 1 , are active galactic nuclei (AGN) of the blazar class. In the very high energy (VHE, > 100 GeV) domain the detection of about 70 AGN is currently summarised in the TeVcat catalog 2 . Again, most of these sources are of the blazar type, i.e. radio-loud AGN such as BL Lac objects where the jet is thought to be inclined at very small viewing angles i to the line of sight. This results in substantial Doppler-boosting of their intrinsic jet emission, S(ν) = D a S (ν ) where D = 1/[Γ b (1 − β b cos i)] denotes the Doppler factor and Γ b = (1 − v 2 b /c 2 ) −1/2 the jet bulk Lorentz factor and typically a ≥ 2, privileging their detection on the sky. Nevertheless, non-blazar AGN such as Radio Galaxies (RGs), while less occurrent, have in the meantime solidly emerged as a new gamma-ray emitting source class as well. With their jets misaligned and associated Doppler boosting effects modest, they enable unique insights into often hidden regions and processes. This review aims at a compact summary of their properties and highlights their role in facilitating theoretical progress in AGN physics. The unification model of radio-loud AGNs postulates that RGs are viewed at a substantial inclination i to the jet axis so that the broad-line optical emitting regions become obscured by a dusty component ("torus" or warped disk) in Narrow Line RGs (NLRGs) such as in Cen A or M87 [22,29]. Depending on their radio structure, RGs have early on been divided into Fanaroff-Riley I and II sources (FR I, FR II) [49], the former one (FR I) encompassing lower radio luminosity, edge-darkened sources and the latter one (FR II) higher luminous, edge-brightened sources where the radio lobes are dominated by bright 1 http://fermi.gsfc.nasa.gov/ssc/data/access/lat/fl8y/ 2 http://tevcat.uchicago.edu arXiv:1810.05409v2 [astro-ph.HE] 16 Nov 2018 hot spots. Various considerations suggest that the high-power FR II sources might be accreting in a "standard" (geometrically thin, optically thick) mode, while most FR I sources are probably supported by a radiatively inefficient accretion flow (RIAF) [56,147]. The general relationship between the blazar and RG class is complex. Urry & Padovani (1995) have described BL Lacs as beamed FR I RGs [145], though evidence exists that the parent population of BL Lac objects contains both FR I and FR II sources [23,82]. A more detailed view might be to posit that X-ray loud BL Lacs (mostly HBLs, peaking in UV/X-rays) are preferentially associated with FR I, while radio-loud BL Lacs (mostly LBL, peaking in the infrared) could show a mixture of FR I and FR II morphologies [115].

Radio Galaxies as VHE emitters -Experimental Status
In the HE range Fermi-LAT has detected about 20 RGs [e.g., 129]. Six of them are also known as VHE emitters (see Fig. 1 This emergence of RGs as a new VHE emitting source class is particularly interesting. Given the substantial misalignment of their jets (i > 10 • ), RGs are commonly thought to be characterized by rather modest Doppler boosting only (bulk Doppler factor D ≤ a few). If, following simple unification considerations, the nuclear emission of FR I type RGs is interpreted as "misaligned BL Lac type" (i.e., of a jet-related, homogeneous synchrotron self-Compton [SSC] origin, yet with small Doppler factor) [45], only a few sources should become detectable at GeV energies (which seems indeed to be the case), and almost none at TeV energies. The discovery of RGs as a new VHE emitting source class thus points to a more complex situation, and promises new insights into some of the fundamental (and often hidden) non-thermal processes in γ-ray emitting AGN. The following aims to provide a short summary of the experimental source characteristics:

PKS 0625-354
The detection of VHE emission from PKS 0625-354 (z = 0.055) above 200 GeV (at a level of ∼ 6σ in 5.5 h of data) has been recently reported by H.E.S.S. [70]. No significant variability is found in the data. The VHE spectrum extends up to ∼ 2 TeV and is compatible with a rather steep power law of photon index Γ γ ∼ −2.8 ± 0.5. The VHE power is moderate with an apparent isotropic luminosity of the order of L VHE ∼ 5 × 10 42 erg/s. Both leptonic and hadronic SED interpretations seem possible [70]. PKS 0625-354 is thought to harbour a black hole of mass M BH ∼ 10 9 M [99] that is probably accreting in an inefficient mode. The source is known as a low excitation line radio-loud AGN, but being a transitional FR I / BL Lac object its proper classification has been debated. Recent findings are favouring its classification as a BL Lac object with non-modest Doppler boosting [101,114,153]. It inclusion in the list of (misaligned) "radio galaxies" may thus have to be re-considered, limiting possible inferences as to the physical origin of its non-thermal emission based on current data [e.g., 51,70].

3C 264
The most recent addition to the RG list has been the FR I source 3C 264 (d ∼ 95 Mpc) seen by VERITAS (with a significance level of ∼ 5.4σ in 12 h of data) [103]. Given an estimated black hole mass M BH ∼ 5 × 10 8 M [127], the VHE luminosity appears to be moderate (∼ 1% of the Crab Nebula) with an isotropic equivalent L(> 300 GeV) ∼ 10 42 erg/sec. 3C 264 has been included in the 3FHL catalog, that lists Fermi-LAT sources which are significantly detected above 10 GeV [14]. The reported VHE flux level seems roughly compatible with a simple power law extrapolation based on the 3FHL results (FHL photon index of −1.65 ± 0.33). There are indications, though, that the VHE spectrum is relatively hard (when compared to other VHE sources) with a photon index close to Γ γ ∼ −2.3. The source shows a low, weakly variable VHE flux along with some month-scale variations. While 3C 264 is known for rapidly evolving knot structures in its jet up to some hundred of parsecs [95], no major knot activity has been observed around the time of the VERITAS observations. Given the previously noted unclear classification of PKS 0625-354, 3C 264 may be the most distant RG detected at VHE so far.

NGC 1275
NGC 1275 (3C 84), the central Perseus cluster RG at a distance of ∼ 70 Mpc, has been detected at VHE energies above 100 GeV by MAGIC, initially (based on data between 2009-2011) at moderate flux levels (∼ 3% of the Crab Nebula) and with a very steep VHE spectrum (photon index of Γ γ ∼ −4.1 if characterized by a power law) extending up to ∼ 650 GeV [18,19]. When HE (Fermi-LAT) and VHE data are combined, the average ("quiescent") γ-ray spectrum appears compatible with either a log-parabola or a power-law with a sub-exponential cut-off, suggestive of a common physical origin and of a peak or cut-off around several GeV. More recently, MAGIC has reported the detection of strong VHE activity with flux levels increased by up to a factor of 50 around December 31, 2016 and January 1, 2017 (reaching ∼ 1.5 of the Crab Nebula or an isotropic equivalent L V HE ∼ 10 45 erg/sec)[96]. Significant day-scale variability has been observed, with the flux doubling timescales as short as ∆t obs ∼ 10 h. The VHE SED measured up to > 1 TeV shows a curved shape (cf. Fig. 2), compatible with an exponential cut-off around a few hundred GeV. The possibility of a joined HE-VHE fit along with day-scale variability, suggests that the HE-VHE emission originates in a (possibly, single) compact zone. The physical nature of this emission is not yet clear, though magnetospheric processes have been favoured over mini-jets-and jet-cloud-interaction scenarios [96]. The central engine in NGC 1275 hosts a black hole of mass M BH ∼ (3 − 4) × 10 8 M [152,154] and exhibits a pc-scale radio jet orientated at i ∼ 30 − 60 • [53,148]. Its inferred jet power is of the order of L j ∼ (0.5 − 1) × 10 44 erg/sec [52,154]. The high ratio L V HE /L j ∼ 10 thus raises questions for a magnetospheric origin of the gamma-ray flare emission, cf. [79] (see also below), unless strong short-term magnetic flux increase would occur. On the other hand, a homogeneous SSC interpretation, assuming the sub-pc scale jet to be weakly misaligned (i < ∼ 20 • ), would be in tension with the inferred jet inclination on pc-scales. This could perhaps be alleviated if the emitting component would e.g. follow a non-straight trajectory that relaxes with distances, and/or if the jet has some internal structure (e.g., spine-shear) allowing for multiple contributions and a more complex inverse Compton interplay [143]. Opacity constraints may pose a severe problem, though (see below). At the moment detailed modelling seems required before firm conclusions can be drawn.

Centaurus A
As the nearest AGN (d 3.7 Mpc) Centaurus A (Cen A) belongs to the best studied extragalactic sources. Its central engine hosts a black hole of mass (0.5 − 1) × 10 8 M [e.g., 107] and emits (assuming a quasar-type SED) a bolometric luminosity of L bol ∼ 10 43 erg/sec [151]. This is much less than the expected Eddington luminosity L Edd and suggests that accretion in its inner part might occur in a radiatively inefficient mode [50,97]. At radio frequencies Cen A has revealed a peculiar morphology including a compact radio core, a sub-pc scale jet and counter-jet, a one-sided kpc-scale jet and inner lobes, up to giant outer lobes with a length of hundreds of kiloparsec. VLBI observations indicate that Cen A is a "non-blazar" source with its inner jet misaligned by i ∼ (12 − 45) • based on TANAMI jet-counter jet flux ratio measurements, and characterized by moderate bulk flow speeds in the radio band of u j < 0.5 c only [e.g., 102]. At VHE energies Cen A has been the second RG detected by H.E.S.S. [10]. A recent, updated analysis based on more than 200h of data shows that the VHE emission extends from 250 GeV up to ∼ 6 TeV and is compatible with a single, rather hard power-law of photon index Γ γ −2.5 ± 0.1 [71]. The source is relatively weak with an equivalent apparent isotropic luminosity of L(> 250 GeV) (1 − 2) × 10 39 erg/s. No significant VHE variability has been found, neither on monthly or yearly timescales, so that an extended origin or contribution (i.e., within the angular resolution ∼ 0.1 • of H.E.S.S., corresponding to ∼ 5 kpc) of the VHE emission cannot per se be discarded. At HE energies, both the core region (i.e., within ∼ 0.1 • ) and the giant lobes of Cen A have been detected by Fermi-LAT [2,3,140,156]. Results concerning the later are indicating that HE lobe emission substantially extends beyond the radio maps. The HE emission of the lobes (which is most likely due to leptonic IC-CMB and IC-EBL, possibly with some additional hadronic pp) is of a particular interest as it provides model-independent information about the spatial distribution of the non-thermal electrons. Fermi-LAT has by now reported extended HE emission from only two RGs, Cen A and Fornax A (d ∼ 20 Mpc) [7]. The core region of Cen A, on the other hand, was initially detected up to 10 GeV (at a level of 4σ) based on ten months of data, with the HE spectrum at that time seemingly compatible with a single power law with photon index Γ p = −2.67 ± 0.1. While this HE power index is very close to the VHE one, a simple extrapolation of the HE power-law was soon found to under-predict the fluxes measured at TeV energies. The comparison was based on non-simultaneous HE and VHE data, but the absence of variability in both energy bands suggested that the discrepancy might be real. Refined analyses based on larger data sets have in the meantime found intriguing evidence for an unusual spectral hardening of the core spectrum by ∆Γ ∼ 0.5 around a few GeV [39,128]. The most recent analysis, involving contemporary VHE and HE data, finds (at a level of 4σ) that the HE spectral index changes around E b 2.8 GeV from Γ γ −2.7 (below E b ) to about −2.3 (above E b ), respectively [71], see Fig. 3. For AGN spectral steepening at gamma-ray energies is a familiar feature that can be related to classical Figure 3. The gamma-ray core spectrum of Cen A above 100 MeV based on 8 yr of Fermi-LAT and more than 200h of H.E.S.S. data. The spectrum shows an unusual spectral hardening at E b 2.8 GeV, with photon index changing by ∼ 0.4 ± 0.1 (assuming a broken power law), see Ref. [71] for details. This spectral feature is most naturally attributed to a second emission component that emerges towards highest energies and that allows to smoothly connect the HE emission (above E b ) with the VHE one. constraints on the acceleration and radiation efficiencies. The observed spectral hardening in Cen A is unusual in this regard; in a "misaligned BL Lac approach" it is best understood as related to the presence of an additional emission component beyond the conventional single-zone SSC-contribution that often satisfactorily describes the SED in blazars. Apart from circumstantial evidence for the blazar Mkn 501 [106,133], Cen A is the first source where spectral results provide clear evidence for the appearance of a physically distinct component above a few GeV. Unfortunately, Cen A is a rather weak γ-ray emitting source, which significantly limits the possibilities to further probe its variability characteristics, particularly above the break. This makes it difficult to observationally disentangle the true nature of the second component with current data. In principle a variety of different (not mutually exclusive) interpretations as to its astrophysical origin are conceivable. Related proposals in the literature operate on different scales (from a few r g to several kpc) and include: (i) (rotational) magnetospheric models that are based on leptonic inverse Compton (IC) processes in an under-luminous accretion environment [121,122], (ii) inner (parsec-scale and below) jet scenarios that invoke differential IC scattering in a stratified jet [57], multiple SSC-emitting components moving at different angles to the line of sight [84] or photo-meson (pγ) interactions of UHE protons in a strong photon field [76,111,130] along with lepto-hadronic combinations [42,117]; alternatively, the hardening could be related to γ-ray induced pair-cascades in a strong disk photon field [137], a dusty torus-like region [126] or the overall host photon field [139]. Moreover, the limited angular resolution of Fermi-LAT and H.E.S.S. (∼ 5 kpc), and the fact that no significant statistical evidence for variability has been found so far, also allows for (iii) scenarios where the emission arises on larger scales; extended scenarios in this context include the interaction of energetic protons with ambient matter (pp) in its kpc-scale region [128], the overall γ-ray contribution of a supposed population of millisecond pulsars [39], or the IC contribution by its kpc-scale jet via up-scattering off various photon fields (e.g., host galaxy starlight or CMB) [68,138], up to more extraordinary explanations invoking the self-annihilation of dark matter particles of mass ∼ 3 TeV within a central dark matter spike [39]. While, given current knowledge, not all of these models are equally likely, and all of them come with some challenges [see e.g., 125], further observational input (such as evidence for VHE variability or extension, the latter now possibly been seen [131]) is needed to better constrain them and help disclosing the real nature of this new component.

IC 310
The Perseus Cluster RG IC 310, located at a distance of d ∼ 80 Mpc (z=0.019), has received particular attention in recent times. The source, originally detected by MAGIC during a campaign in 2009-2010 [17], has shown extreme VHE variability during a strong flare in November 2012, revealing VHE flux variations on timescales as short as ∆t 5 min [20], see Fig. 4. The 2012 VHE flare spectrum appears compatible with a single, hard power law of photon index Γ γ > ∼ − 2 (and possibly as low as ∼ −1.5) over a range from 70 GeV to 8.3 TeV, with no indications of any internal absorption [see also 12]. The source can reach high VHE flux levels, corresponding to an isotropic-equivalent luminosity of L V HE 2 × 10 44 erg s −1 . IC 310 is commonly believed [e.g., 20] to harbour a black hole of mass M BH 3 × 10 8 M [but see also Ref. 32, for a ten times smaller estimate] and has for some time been classified as a head-tail RG. The apparent lack of jet bending along with more recent indications for a one-sided pc-scale radio jet inclined at i < ∼ 38 • suggests, however, that IC 310 is a transitional source at the borderline dividing low-luminosity RGs and BL Lac objects [77]. The extreme VHE variability along with the high VHE power ( > ∼ L Edd /200) and the hard γ-ray spectrum are surprising findings for a misaligned source. Based on a variety of considerations, including the orientation of its jet (probably i ∼ [10 − 20] • ) as well as kinetic jet power and timing constraints, Aleksić et al. [20] have disfavoured several alternative models for rapid VHE variability such as magnetic reconnection [e.g., 62] or jet-cloud/star interaction [e.g., 27]. This inference is, however, less robust as has been shown later on [cf. 11, for details]. Nevertheless, the fact that the VHE flux varies on timescales ∆t much shorter than the light travel time across the black hole horizon, r g (3 × 10 8 M )/c = 25 min, has been interpreted as evidence for the occurrence of gap-type particle acceleration on sub-horizon scales, i.e. in unscreened electric field regions ("gaps") of height h 0.2r g [e.g., 20,73]. Questions concerning such an interpretation are related to the fact that the characteristic VHE power of a (steady) gap scales with the jet power, L V HE ∼ L j (h/r g ) a , a = 2 − 4 [79], the latter of which is known to be rather modest on average for IC 310, i.e. L j ∼ 10 43 erg s −1 [cf. also , 134]. Unless strong (short-term) magnetic flux increases would occur, the expected gap output would under-predict the VHE fluxes measured during the flaring state. IC 310 has subsequently (after November 2012) shown a rather low TeV emission state with a steeper spectrum (Γ ∼ −2.4) measured up to ∼ 3 TeV and with little evidence for variability. The multi-wavelength SED during this state appears to be satisfactorily reproducible with a one-zone SSC model using parameters that are comparable to those found for other misaligned, γ-ray emitting AGN [12].

M87
The Virgo Cluster galaxy M87 (NGC 4486) has been the first extragalactic source detected at VHE energies [8]. Classified as low-excitation, weak-power FR I source, M87 hosts one of the most massive black holes of M BH (2 − 6) × 10 9 M [e.g., 149], and is thought to be accreting in a radiatively inefficient (RIAF) mode [118]. Given its proximity at a distance of d 16.4 Mpc [34] and its large mass-scale r g , M87 has become a prominent target to probe jet formation scenarios with high-resolution radio observations down to scales of tens of gravitational radii [e.g., 15,48,67,81,98]. Its sub-parsec scale radio jet appears misaligned by an angle i ∼ (15 − 25) • and shows a rather complex structure, seemingly compatible with a slower, mildly relativistic (β ∼ 0.5c) layer and a faster moving, relativistic spine (Γ b ∼ 2.5)[see e.g., 98]. Indications of a parabolic jet shape suggest that the jet initially experiences some external confinement as by a disk wind [59]. In general the inferred jet seeds and inclinations are consistent with rather modest Doppler factors D < ∼ a few (for review, see e.g. [124]). At VHE energies, M87 is well known for its rapid day-scale variability (flux doubling time scales ∆t obs ∼ 1 d) during active source states, and a rather hard, featureless photon spectrum compatible with a single power law (of index Γ γ = −2.2 ± 0.2 in high, and somewhat steeper Γ γ ∼ −2.6 in low states) extending from ∼ 300 GeV up to ∼ 10 TeV [4,5,9,16,21]. Both, the observed rapid VHE variability and the hard VHE spectrum are remarkable features for a misaligned AGN, and reminiscent of those seen in IC 310. Based on the first 10 months of data, Fermi-LAT has reported HE gamma-ray emission from M87 up to 30 GeV [1] with a photon spectrum then seemingly compatible with a single power-law of index Γ γ = −2.26 ± 0.13 and comparable to the one(s) in the VHE high states. Nevertheless, a simple extrapolation of this HE power-law to the VHE regime turned out to be insufficient to account for the flux levels measured during the TeV high states (up to equivalent levels of L(> 350 GeV) ∼ 5 × 10 41 erg/s, e.g. [21]), suggesting that the high states might be accompanied by the emergence of an additional component [124]. No evidence for significant flux variations (down to timescales of 10 days) has been found during these early HE observations, though on experimental grounds the occurrence of shorter-timescale variations cannot per se be excluded. Similar spectral results have been reported in the 3FGL catalog (4 yr of data), with the HE spectrum below 10 GeV compatible with a single power-law of Γ γ = −2.04 ± 0.07 [6], but with indications for a possible change above 10 GeV. The most recent analysis based on ∼ 8 yr of Fermi-LAT data reports evidence for month-type HE variability and indications for excess emission over the standard power-law model above ∼ 10 GeV, similar to earlier findings in Cen A [13], see also Fig. 5. When viewed in an HE-VHE context, these findings are most naturally explained by an additional emission component that dominates the highest-energy part of the spectrum and allows for a smooth HE-VHE spectral connection. As the HE spectrum extends to about 100 GeV without indications for a cut-off and the VHE thresholds reach down to about 200 GeV, variability seen with high statistics at VHE can be used to constrain the nature of this additional component. This contrasts with Cen A where no significant VHE variability has been found yet. For M87 current findings do support proposals in which the emission arises on innermost jet scales and below. Figure 5. Gamma-ray SED for M87 based on ∼ 8 yr of Fermi-LAT data including the different observed VHE states. The average ("regular") spectrum shows a break in the SED around ∼ 10 GeV, suggestive of an additional HE component. The break appears masked in the "high state" by flaring above ∼ 10 GeV. The current situation now in principle allows for a smooth connection of HE and VHE states. This suggests that the nature of this additional component is constrained by the observed VHE variability. The light grey curves for the two components are intended to guide the eyes only. Following Ref. [13].
Light travel time arguments in fact point to a compact VHE emission region (R < c∆t obs D) in M87 of a size comparable to the Schwarzschild radius r s = (0.6 − 1.8) × 10 15 cm of its black hole. Similar as for Cen A, a variety of models have been introduced to account for this, cf. Fig. 6 for an exemplary illustration (see Refs. [27,54,61,84,87,116,117,120,142]). The interested reader is referred to Refs. [123,124] for a more detailed description and discussion of them. M87 has been repeatedly active over the past ten years, with VHE high states being detected in 2005, 2008 and 2010, and an elevated one (flux levels 2-3 times higher than average) in 2012. Interestingly, during all high states, day-scale VHE variability has been found. The 2012 monitoring data by VERITAS do not reveal a bright flare, but the light curve indicates VHE variability on timescales of (at least) weeks, suggesting that the often called "quiescent" state also shows some longterm evolution [31], cf. also [13]. No major VHE flare has been seen since then, though hints for day-scale variability in 2013 have been reported [25]. As the angular resolution of current VHE instruments is limited (to scales of ∼ 25 kpc for M87), coordinated VLBI radio observations, capable of probing down to scales of tens of gravitational radii, have been performed during the 2008, 2010 and 2012 high VHE states. These results indicate that the TeV emission is accompanied by (delayed) radio core flux enhancements, supporting proposals that the VHE emission originate at the jet base very near to the black hole [5,15,65,66]. The radio-VHE correlation along with the required compactness of the VHE zone have served as a strong motivation to explore plasma injection via gap-type magnetospheric processes close to the black hole [38,87,113]. Concerning its rapid VHE variability M87 shows some similarities with IC 310, though its associated VHE luminosity output is by more than a factor of 100 smaller than the one for IC 310. EC starlight photons (e.g. Stawarz+2006) Figure 6. Possible scenarios for the origin of the variable VHE emission in M87 with exemplary references. Day-scale variability favours models on scales of the inner jet and below. "IC" refers to inverse Compton radiation.

Models for the HE-VHE γ-ray emission
In the following we shortly describe and comment on some recent theoretical trends and developments concerning the origin of the VHE emission beyond conventional jet (one-zone) synchrotron self-Compton (SSC) models. We distinguish between scenarios aimed at addressing the highly variable VHE part (as in M87) and those focusing on an apparently steady VHE part (as for Cen A).

Variable VHE and black hole gap models
The activation of Blandford-Znajek (BZ) outflows requires continuous injection of plasma in the magnetospheric region enclosed between the inner and outer light cylinders, the origin of which is yet an open issue. To fully screen out the magnetosphere, the plasma injection rate must be sufficiently high to maintain the density everywhere in the magnetosphere above the Goldreich-Julian (GJ) value. Whether the dense accretion flow surrounding the black hole can provide sufficient charges for complete screening is unclear at present; direct feeding seems unlikely, as charged particles would have to cross magnetic field lines on a timescale much shorter than the accretion time in order to reach the polar outflow. Plasma injection by virtue of macroscopic instabilities, that might lead to re-arrangement of the magnetic field configuration, is a possibility, however, the growth time of such instabilities may be much longer than the gap evacuation time, ∼ r g /c, and it is quite likely that even if occasional injection of plasma into the casual section of the magnetospheric does occur, the plasma density may not be sustained above the required level at all times. An alternative charge supply mechanism is pair creation on magnetic field lines via annihilation of MeV photons that emanate from a RIAF during low accretion states, or from a putative accretion disk corona during intermediate states. To estimate the γ-ray luminosity required for complete screening we note that the density of pairs thereby created is roughly n ± σ γγ n 2 γ r g /3 [87], where σ γγ is the pair production cross section, n γ 10 22 m −1R−2 γ l γ cm −3 is the density of MeV photons, l γ = L γ /L Edd and R =Rr g are, respectively, the Eddington ratio and radius of the radiation source, and m = M/M is the black hole mass in solar mass units. This should be compared with the GJ density, n GJ = ΩB/(2πec) = 2 × 10 11 B 8 (Ω/ω H )m −1 cm −3 , expressed here in terms of the angular velocity of magnetic surfaces Ω, the angular velocity of the black hole ω H c/2r g , the strength of the magnetic field near the horizon B = 10 8 B 8 Gauss, and the magnitude of the electron charge e > 0. The requirement n ± > n GJ implies For sources accreting in the RIAF regime, the gamma-ray luminosity can be related to the accretion rate, albeit with a large uncertainty, using an ADAF model. The strength of the magnetic field near the horizon also scales with the accretion rate, roughly as B 10 9 (ṁ/m) 1/2 G, whereṁ =Ṁ/Ṁ Edd , andṀ Edd = 10L Edd /c 2 . These two relations can be combined to yield the critical accretion rate below which the magnetosphere is expected to be starved, For M87, where m = 6 × 10 9 (up to a factor of two), this implies starvation atṁ < ∼ 10 −4 , which seems to be above its inferred accretion rate. In sources accreting well above the critical ADAF rate, the accretion flow is anticipated to be cold, and the emission spectrum is unlikely to extend to energies above the electron mass. Gamma-rays may, nonetheless, originate from a tenuous corona, if present. No reliable constraints on the spectrum and luminosity of this coronal component have been imposed thus far. In principle, it could be that in sources that accrete at relatively high rates the magnetic field is much higher than in RIAF sources, while the gamma-ray luminosity (of the corona) is smaller. If indeed true, it could mean that gap emission in such objects may be more intense than in RIAF sources. Equation (2) has been employed to show [73,79,87] that under conditions likely to prevail in many stellar and supermassive black hole systems, the annihilation rate of disk photons is insufficient to maintain the charge density in the magnetosphere at the GJ value, giving rise to formation of spark gaps [33,35,72,85]. It has been further pointed out [73,74,79,87,91,105,122] that the gap activity may be imprinted in the high-energy emission observed in these sources, whereby the variable TeV emission detected in M87 [5,9] and IC310 [20] was speculated to constitute examples of the signature of magnetospheric plasma production on horizon scales [73,79,85,87,105,113]. In what follows we provide a concise overview of recent black hole gap models.

Stationary gap models
Stationary models tacitly assume that a gap forms in a localized region of the magnetosphere outside which the ideal MHD condition (that is, E · B = 0) prevails. To gain insight it is instructive to derive the gap electrostatic equation in flat spacetime first. The generalization to Kerr spacetime then readily follows. In general, the angular velocity of magnetic field lines, Ω = Ωẑ, is conserved along magnetic surfaces only in regions where the ideal MHD condition is satisfied. However, under the assumption that the gap constitutes a small disturbance in the global magnetosphere, the variation of Ω across the gap (due to the finite potential drop) can be ignored. It is then convenient to transform to a rotating coordinate system, t = t, r = r, θ = θ and ϕ = ϕ − Ω t, here in spherical coordinates, in which the electric and magnetic fields are given by B = B, E = E + v × B in terms of their components in the non-rotating frame (using geometric units, c = 1), where v = Ω × r is the tangential velocity of the magnetic flux tube. It can be readily shown that in the rotating frame Gauss's law takes the form where ρ e is the charge density and ρ GJ = −∇ · (v × B)/4π denotes the Goldreich-Julian (GJ) density. In dipolar and split monopole geometries this simplifies to the well known result ρ GJ = −Ω · B/2π.
Equation (3) indicates that a charge density ρ e = ρ GJ is needed to screen out a parallel E || . The generalization to Kerr geometry is straightforward. In Boyer-Lindquist coordinates, (t, r, θ, ϕ), the electric field component in the rotating coordinate frame is F µt = F µt + ΩF µϕ , where F µν denotes the electromagnetic tensor, and we henceforth adopt the metric signature (−, +, +, +). The GR equations can be recast in a form similar to that in flat spacetime upon using quantities measured in the local ZAMO (zero angular momentum observer) frame. For illustration, we adopt a split monopole geometry, defined by the potential A ϕ = Ψ H (1 − cos θ), with Ψ H being the magnetic flux on the horizon. The non-corotating electric field, measured by a ZAMO, is then given by E r = √ AF rt /Σ in terms of the metric functions (in geometric units) Σ = r 2 + a 2 cos 2 θ and A = (r 2 + a 2 ) 2 − a 2 ∆ sin 2 θ, with ∆ = r 2 + a 2 − 2Mr, M and a being the black hole mass and specific angular momentum (spin parameter), respectively. The general relativistic Gauss's law reads: where is the general relativistic (GR) version of the GJ density, α = (Σ∆/A) 1/2 is the lapse function and ω = 2aMr/A is the ZAMO angular velocity. In flat spacetime √ A = Σ = r 2 , and Eq. (4) reduces to Eq. (3) in spherical coordinates. The charge density is related to the proper densities, n ± , and the time component of the 4-velocities, u 0 ± , of the e ± pairs through ρ e = e(n + u 0 (4) is subject to the condition E = 0 at the inner and outer gap boundaries, which follows from the assertion that the flow outside the gap is in a force-free state. This implies that E r does not change sign, and that |E r | has a maximum inside the gap (see inset in Fig. 7). An immediate consequence is that ρ GJ must change sign inside the gap [73,89]. Thus, the gap should form around the null surface on which ρ GJ vanishes. Such a null surface is a distinctive feature of a Kerr black hole that results from the frame dragging effect (see Fig. 7) 3 . Now, from the above it is seen that at the null surface the gap electric field scales as |E r | ∝ h 2 with the gap width h, and the potential as V ∝ h 3 . The maximum electric current flowing through the gap is limited by j r = |ρ GJ (r 1 )|c ∼ |dρ/dr| 0 h/2, where r 1 is the outer gap boundary and the subscript 0 designates values at the null surface. Thus, the maximum power that can be tapped, L gap 2πr 2 H jV, and, hence, the maximum gamma-ray luminosity that can be emitted by the gap, scale as h 4 . A more precise expression for L gap is derived in Ref. [73], see also [79]. Since the gap closure condition restricts the multiplicity inside the gap to unity, the gap width h increases with decreasing pair production opacity. The salient lesson is that in steady gap models the output power of the gap increases steeply with decreasing disk luminosity. Detectability of gap emission favours low luminosity sources. We shall see shortly that this is not true in case of intermittent gaps. A quantitative treatment of gap emission requires inclusion of plasma dynamics, radiation back-reaction, Compton scattering and pair production. Details are given in Refs. [73][74][75]89]. It is found that the characteristic spectrum produced in the gap consists of two components; curvature emission that peaks at sub-TeV energies, and inverse Compton emission that peaks at 10 to 100 TeV, depending on the black hole mass. Refs. [73][74][75]91] predict that gap emission from stellar and supermassive black holes should be detected by upcoming experiments.

Time-dependent models
The main difficulty with steady gap solutions for the null surface is that they exist only under highly restricted conditions, that may not apply to most objects [89]. The reason is that the gap closure condition (i.e., unit multiplicity) imposes a relation between the global magnetospheric current 3 A null surface exists in a pulsar outer gap for other reasons (through the gap width) and the pair production opacity, which in reality are two independent and unrelated quantities. An exemplary illustration for the existence regime of steady gap solutions is displayed in Fig. 8, for a supermassive black hole of mass m = 10 9 and a power-law disk emission spectrum. Similar results are obtained for stellar mass black holes. This example indicates that in practice, steady gaps can form only in sources with extremely low Eddington ratios. Additional difficulty stems from an inconsistency between the flow direction at the outer gap boundary and the ideal MHD flow below the stagnation surface [58,89]. Finally, the stability of steady gaps is questionable. The main conclusion is that sparking of starved magnetospheric regions is likely to be inherently intermittent. Attempts to construct 1D time-dependent models have been reported recently [40,90]. The analysis described in Ref. [90] was performed using a newly developed, fully GR (in Kerr geometry) particle-in-cell code that implements Monte-Carlo methods to compute the interaction of pairs and gamma-rays with the soft photons emitted by the accretion flow. As in the steady models mentioned above, the gap is assumed to be a small disturbance that does not affect the global magnetospheric structure. The evolution of the electric field is governed by the equation that generalizes the flat spacetime model derived in Ref. [86]. Here j r = e(n + u r + − n − u r − ) is the radial electric current density, expressed in terms of the proper electron (n − ) and positron (n + ) densities and their 4-velocities u − and u + , and J 0 represents the global magnetospheric current, which is an input parameter of the model in addition to Ω. The initial condition, E r (t = 0, r), is obtained upon solving Eq. (4) at the beginning of each run. For details the reader is referred to Ref. [90]. The analysis indicates that when the Thomson length for collision with disk photons becomes smaller than the gap width, i.e. τ 0 > 1, screening of the gap occurs, following a prompt discharge phase that exhausts the initial energy stored in the gap, through low amplitude, rapid plasma oscillations that produce self-sustained pair cascades, with quasi-stationary pair and gamma-ray spectra. An example is exhibited in Figs. 9 Figure 8. Maximum Eddington ratio below which local steady gap solutions exist, versus normalized magnetospheric current, 2π J 0 /(Σ null ΩB H cos θ), where Σ null is the value of Σ on the null surface, for a supermassive black hole of mass M BH = 10 9 M , magnetic field B = 10 4 G, a radiation source of size R = 30r g and a power law spectrum, I ν ∝ ν −2 , with a low energy cutoff at ν min = 10 12 Hz, and inclination angle θ = 30 • . (See Ref. [89] for further details). Figure 9. Evolution of the gap electric field in a magnetosphere of a supermassive black hole with the same parameters as in Fig. 8: The left panel shows the electric flux Φ E = √ AE r at the initial time t = 0. The right panel shows the electric flux at time t = 20r g /c, roughly the onset of the quasi-steady oscillations. Note the change of scale on the vertical axis between the two panels. The amplitude of E r during the quasi-steady oscillations is about 10 −3 of its initial magnitude at the null surface. The fiducial opacity in this case is τ 0 = 10. and 10. The gamma-ray spectrum emitted from the gap peaks in the TeV band (Fig. 10), with a total luminosity that constitutes a fraction of about 10 −5 of the corresponding Blandford-Znajek power. This seems different than the spectrum predicted by steady gap models. As those simulations are demanding, only a small range of parameters was explored in [90]. Future studies should investigate how the emission properties depend on the black hole mass, the target radiation spectrum and the magnetic field strength. A caveat of the model described above is that it does not account for the effect of the ideal MHD fields on the dynamics of the plasma. An attempt to include such additional forces has been made in Ref. [40], where it has been argued that the gap dynamics is cyclic in nature rather than quasi-steady. However, the analysis in this work is based on a mirror gap model in flat spacetime, that inherently assumes the presence of a current sheet (or current side flow) at the stagnation surface (by reversal of the toroidal magnetic field there). Whether this model is a reasonable representation of a starved Kerr black hole magnetosphere remains to be demonstrated. A cyclic behaviour is naively anticipated in a global model that links the gap and the magnetosphere together [89]. How this would affect the spectrum and light curves of the HE gap emission is yet an open issue that can only be studied using 2D simulations. ct/r g =20.06 Figure 10. Evolution of the electron (blue), positron (red) and photon (yellow) spectra for the case shown in Fig. 9. The photon energy is given in m e c 2 units. The peak of the quasi-steady gamma-ray spectrum in this example (rightmost panel) is at ∼ 1 TeV.

Variable VHE and inner jet models
Rapid variability on timescales of days and below implies a compact emitting region, and is often taken to indicate that this emitting region is located at sub-parsec jet scales (< 10 4 r g ) and below. Promising jet-related proposals in this context (beyond black hole gaps) include the following (cf. Fig. 11):

Spine-shear scenarios
The complex interplay of black hole-and disk-driven outflows along with environmental interaction (e.g., entrainment) is likely to give rise to a non-uniform flow topology where different jet layers possess different bulk flow speeds. In its simplest realization a fast (Γ 1 1) spine is taken to be surrounded by a slower-moving (Γ 2 < Γ 1 ) sheath or layer [e.g., 57]. This is essentially a two-zone model, yet with parameter-space constrained by the requirement of an efficient radiative interplay. Depending on the viewing angle we may see different parts of the jet (spine or sheath), resulting in different emission characteristics for misaligned AGN when compared with those of blazars. In particular, given their larger inclination angles, the VHE emission in RGs would usually be dominated by the stronger (!) Doppler-boosted emission from the layer. VHE variability then essentially imposes constraints on the size of the layer. Spine-sheath models can in principle accommodate fast variability and facilitate an increased gamma-ray luminosity as each component sees the (external) radiation of the other amplified by the relative motion between them (Γ rel = Γ 1 Γ 2 [1 − β 1 β 2 ]), thereby enhancing the IC contribution of both its spine and its layer [57,142,143]. In order for this interplay to work efficiently however, the emitting zones need to be (quasi) co-spatial. This in turn leads to transparency issues as long as one aims at a simultaneous reproduction of the infrared-optical part of the SED, since the large gamma-ray opacity associated with the intense (weakly Doppler-boosted) radiation field of the spine usually results in considerable γγ-absorption and a steep slope at TeV energies. Hard TeV spectra thus cannot simply be accounted for with such models [e.g., 96]. Nevertheless, an internal velocity stratification (shear) bears the potential to solve several issues regarding the unification of BL Lacs and RGs [44,94,132]. The noted emission models represent a simplified first approach, that focus on the radiative interplay only and do not yet take any acceleration effects [e.g., 119] nor time-dependencies into account. Decoupling of the infrared and VHE emission parts would alleviate some tensions, though it would leave a larger part of the parameter space unconstrained. Similar to other two-zone approaches, the non-thermal emission can in many cases be reproduced under equipartition conditions which may be counted in its favour [144]. Time-dependent extensions are certainly needed to further assess its potential in the context of RG modelling.

Reconnection -mini-jets and plasmoids
Relativistic jets are generally expected to be initially magnetically dominated (magnetization σ m 1), and this has motivated studies where the field energy is released via collisionless magnetic reconnection allowing for non-thermal particle acceleration and subsequent high energy gamma-ray (IC) production [e.g. 41,[60][61][62]78,104,135,150]. An open question at present concerns whether formation of current sheets on sufficiently small scales can result from current driven instabilities induced during the propagation of the jet [e.g., 63,100,108], or may inherently form during injection of the jet, e.g., due to advection of asymmetric magnetic field or magnetic loops by the disk [e.g., 88,109]. The observed collimation profile of M87 seems to indicate that in this source the jet is kink stable [59]. Upcoming EHT observations on horizon scale may shed more light on this issue. Reconnection-related models of the type proposed in [61] assume that relativistic (Petschek-type) reconnection occurs in the jets of RGs (taking them to be electron-proton dominated with σ m ∼ 100), leading to efficient electron acceleration and the generation of variable VHE gamma-rays via inverse Compton (SSC and/or EC) scattering [e.g., 46]. Moreover, efficient reconnection allows for an additional relativistic velocity component of the ejected plasma (Γ r √ σ m ) relative to the mean bulk flow of the jet, and directed at some angle to it. Given multiple, localised reconnection sites ("mini-jets") within the jet, strong differential Doppler boosting effects (D 1) could become possible even for RG sources whose main jet direction are substantially misaligned. This could then provide a simple explanation for ultra-fast variability in misaligned AGN. Caveats concern whether such high magnetisations should indeed be expected for electron-proton (disk-driven) jets, and whether the impact of a magnetic guide field, that would lead to weak dissipation only [93], can be neglected on the anticipated scale of VHE gamma-ray production. Nevertheless, reconnection is particularly interesting as it could facilitate non-thermal particle energization in those magnetized part of the jets where diffusive shock acceleration is inefficient [e.g., 83], and ensure a rough equipartition between magnetic fields and radiating particles in the reconnection downstream (emitting) region [135]. In the absence of a guide field, relativistic reconnection results in a flat (hard) power-law tail n(γ) ∼ γ −α with slope α < 2 for σ m > 10, most probably approaching α ∼ 1 in the extreme relativistic case [64,150,158]. For large Lundquist numbers the reconnection layer becomes unstable to resistive (tearing) instabilities, causing its fragmentation into local magnetic islands or plasmoids. More recent studies have thus concentrated on characterizing the (fractal-like nature of) plasmoid formation in relativistic pair plasmas [112,136].

Jet-star/cloud interactions
The ubiquity of stars and gas in the central region of the AGN host galaxy can lead to frequent jet-star/cloud interactions, facilitate jet entrainment and mass loading by stellar winds [e.g., 110], and give rise to knotty X-ray structures as e.g. seen in the large-scale jet of Cen A [155]. Over the last few years a variety of studies have explored its impact on the gamma-ray emission characteristics in AGN, both for its steady and transient/flaring states [e.g. 24,26,30,36,37,47,80,146,157]. In the RG context, one interesting application concerns the possible contribution of inelastic proton-proton (pp) collisions in the generation of variable VHE emission. While AGN jets are commonly considered not to carry enough target matter (n p ) to allow for efficient pp-collisions (given its long timescale t pp 10 15 /n p sec), interactions of a red giant star or a massive, dense gas cloud (size r o ) with the base of the jet (radius r j ) could occasionally introduce a high amount of matter, trigger shock acceleration and potentially drive rapid VHE activity [26,27]. Model calculations in the case of M87 [28] suggest that such a scenario could account for the observed VHE characteristics (including day-scale variability) if the jet would be powerful enough and a sufficiently large fraction (∝ r 2 o /r 2 j ) of it could be channeled into VHE γ-ray production. Given the observed large opening angle and transversal dimension r of the milli-arcsec radio jet in M87, the latter is not obvious. In fact, simple models of this type often need jet powers in excess of current estimates. In principle, however, these power constraints could be somewhat relaxed if the jets possess a spine-shear-type configuration with most of the energy flux concentrated into a narrow core. This would then bear some similarities with the set-up discussed in Sec. 3.2.1. In addition, the effective size r 0 (if related to shocks at quite some distance from the star) could well be larger than the initial size of the obstacle [e.g., 30,37]. While the availability of suitable stellar orbits on sub-parsec scales limits the possible recurrence (frequency) of short, star-driven VHE flaring events, multiple collisions along the jet seem unavoidable. Depending on conditions in the host galaxy, this could also result in a detectable, steady γ-ray contribution at VHE energies; a recent calculation of the cumulative emission from multiple jet -stellar wind interactions in M87, however, suggests this to be too low to account for its overall VHE flux levels [146].

Steady VHE and extended jet models
The detection of extended X-ray emission from the large-scale jets in AGN by Chandra 4 has raised the possibility that these jets are also steady sources of VHE γ-rays. Electron synchrotron radiation is by now the favoured interpretation for the X-ray emission [e.g., 55,69], indicating the presence of highly energetic electrons with Lorentz factors up to γ ∼ 10 8 (100 µG/B) 1/2 D −1/2 , where D is the Doppler factor and B the large-scale magnetic field strength. While diffusive shock acceleration in the jet could in principle facilitate such energies, localized (shock-type) acceleration at knots is usually not sufficient given the fact that there is little evidence for e.g. the inter-knot regions to have significantly steeper spectra than the adjacent knots as one would expect in the case of synchrotron cooling. This may point to the operation of a continued or distributed acceleration mechanisms such a stochastic or shear particle acceleration [e.g., 92]. Figure 12 shows a SED result of a recent spectral analysis in M87 [141]. The difference of the radio and X-ray spectral indices supports a synchrotron as opposed to an IC origin. The multi-TeV electrons in the jets will then also up-scatter the host starlight, dust, EBL or CMB photons to γ-ray energies, resulting in a weak and steady VHE contribution that may become detectable in nearby RG. IC up-scattering of CMB (u CMB = 4.2 × 10 −13 erg cm −3 ; ν p 1.6 × 10 11 Hz), for example, leads to a non-reducible multi-TeV contribution at flux levels of f IC ∼ (u IC /u B ) f x ∼ 10 −3 f x . For M87 however, the total X-ray flux of the knots is of the order f x ∼ 10 −12 erg cm −2 s −1 only, implying an IC-CMB contribution at flux levels well below current detections. Up-scattering of dust or starlight photons on the other hand, is expected to yield a much higher VHE contribution that could be probed with CTA [68]. Given the (current) absence of VHE variability in Cen A, an extended (leptonic) origin of its VHE emission certainly cannot be excluded. In addition, the cumulative IC emission from multiple jet-star collisions in the kpc-scale jet [30] might further contribute. Evidence for a possible VHE extension has in fact been recently reported for Cen A [131]. We note that an experimental verification of extended VHE emission would support the notion that the large-scale jets in AGN could make a relevant contribution to the TeV background when compared to the highly-boosted VHE emission from blazar cores [55].

Conclusions
The experimental progress over the last decade has led to the discovery of radio galaxies (RGs) at γ-ray energies, revealing exceptional features such as spectral hardening or ultra-fast variability. With their jets misaligned and related Doppler boosting effects only modest, RGs are offering unique insights into physical mechanisms and environments (e.g., the plasma physics of jets or the black hole vicinity) that are otherwise difficult to access. The unexpected spectral hardening at gamma-ray energies seen in Cen A, for example, points to the emergence of a new physical component beyond the conventional SSC-type one, with current interpretations ranging from the smallest (sub-pc) to the largest (kpc) jet scales, up to extended dark matter scenarios. On the other hand, rapid VHE variability on timescales shorter or comparable to the light travel time across the horizon of the black hole, as e.g. seen in IC 310 and M87, provides evidence for a highly compact emitting zone, probably at the very origin of the jet itself (black hole gaps), or internally-(reconnection) or externally (star collision)-induced. These and related experimental findings at gamma-ray energies have been vital in triggering important conceptual progress in black hole -jet physics and generated a variety of promising research avenues. Along with further theoretical efforts, dedicated observational studies are currently needed to e.g. clarify source classification and to better resolve their timing characteristics. The upcoming CTA array [43] will have the potential to probe deeper into the spectral and variability characteristics of RGs, and thereby allow to fundamentally advance our understanding of the AGN phenomena in general.