Gamma-Ray Astrophysics in the Time Domain

The last few years have seen gamma-ray astronomy maturing and advancing in the field of time-domain astronomy, utilizing source variability on timescales over many orders of magnitudes, from a decade down to a few minutes and shorter, depending on the source. This review focuses on some of the key science issues and conceptual developments concerning the timing characteristics of active galactic nuclei (AGN) at gamma-ray energies. It highlights the relevance of adequate statistical tools and illustrates that the developments in the gamma-ray domain bear the potential to fundamentally deepen our understanding of the nature of the emitting source and the link between accretion dynamics, black hole physics, and jet ejection.


Introduction
The last decade has seen tremendous experimental progress in gamma-ray astronomy much beyond simple source detection (e.g., see [1][2][3] for a review). In many cases, detailed spectral and timing characterization have become possible allowing one to probe deeply into the nature and physics of these sources. In particular, gamma-ray astronomy has by now matured and progressed further in the field of time-domain astronomy, utilizing source variability on timescales over many orders of magnitudes, from a decade down to a few minutes and shorter. Instruments such as the Fermi Large Area Telescope (Fermi-LAT), the High-Altitude Water Cherenkov Gamma-Ray Observatory (HAWC) or the First G-APD Cherenkov Telescope (FACT), for example, have opened up the possibility for unbiased long-term (timescales up to several years) studies of bright astrophysical objects in the high energy (HE; >100 MeV) and the very high energy (VHE; >100 GeV) domain, respectively (e.g., [4][5][6]), while modern Imaging Atmospheric Cherenkov Telescopes (IACTs) have demonstrated their excellent capabilities to characterize VHE flaring states down to below sub-hour timescales.
This paper focuses on some of the key issues and conceptual developments concerning the timing characteristics of AGN (including rapid variability, log-normal flux distributions, power-law noise, and quasi-periodic oscillations) at gamma-ray energies.

Rapid Variability at VHE Energies
Radio-loud AGN can be highly variable gamma-ray emitters. This not only applies to the strongly-Doppler-boosted blazar sources, but also to misaligned jet sources such as radio galaxies. VHE doubling timescales down to a few minutes (∆t obs ∼ 2-3 min) are known for the VHE blazars Mkn 501 [7] and PKS 2155-304 [8] and at HE for the flat spectrum radio quasar (FSRQ) blazar 3C279 [9]. Interestingly, rapid VHE activity has been seen in radio galaxies as well, i.e., from day-scale for M87 [10,11], see Figure 1, to intra-day (∆t obs ∼ 10 h) for NGC 1275 [12], down to minute-scale for IC 310 [13]; see [14]. Rapid variability is usually taken to indicate extreme jet conditions. Causality arguments, for example, imply that the emission would need to arise in a very compact region of comoving size ∆r < ∼ δc∆t obs = 10 15 (δ/10)(∆t obs /1 h) cm (where δ is the Doppler factor), possibly close to the black hole at distances d < ∼ 2cγ 2 b ∆t obs 2 × 10 16 (γ b /10) 2 (∆t obs /1 h) cm (where γ b ≥ 1 is the jet bulk Lorentz factor) only, i.e., typically less than a few hundred Schwarzschild radii away. Despite its compactness, the emitting region often needs to be highly luminous to account for the observed gamma-ray output, and this can impose a severe constraint on possible models (see below).
In some cases such as for PKS 2155-304, the VHE light curve during the high state is sufficiently resolved to suggest the presence of multiple (possibly interacting) emitting zones. The discriminating potential will strongly increase with the upcoming Cherenkov Telescope Array (CTA). Figure 2 shows a simulation based on the strong VHE flare in 2006, demonstrating the capability to probe sub-minute timescales. Significant day-scale activity is evident. Since the jet in M87 is believed to be misaligned (suggestive of modest Doppler beaming only) and its black hole light-crossing scale large (r s /c > ∼ 0.5 days), this VHE activity is extremely fast. The curve represents a fit by an exponential function, indicating doubling times of ∼ 1.7 days and ∼ 0.6 days during rise and decay, respectively. From [11]. In the extragalactic jetted-source context, rapid variability comes along with a general challenge. In the conventional picture of an emission region embedded in a uni-directional jet flow, the observed variability imposes a constraint on the size of the central engine that is independent of Doppler boosting (e.g., [16,17]). To see this, note that the observed timescale ∆t obs is related to the comoving (flow) timescale ∆t by ∆t obs = (1 + z)∆t /δ, where δ = 1/(γ b [1 − β cos i]) is the Doppler factor, and that by causality ∆t > ∼ ∆r /c. Length scales in the rest frame of the central engine and the comoving (flow) frame, on the other hand, are related by ∆r = ∆r /γ b (length contraction). Noting that the characteristic length scale in the rest frame of the central engine driving the jet is ∆r ∼ r s ≡ 2GM BH /c 2 , one finds for the generic variability timescale: using that δ = 2γ b in the head-on (i = 0) approximation. The ∼ 3-min VHE variability observed in PKS 2155-304 (z = 0.116) would then imply a black hole mass of M BH < ∼ 4 × 10 7 M only, much smaller than what is inferred from the host-galaxy luminosity relation (e.g., [18]). Faster variability could only be achieved if the jet collided with a small obstacle (size r c r s ). Yet, in this case, the jet power that could be tapped for producing radiation is much smaller, roughly by a factor (r c /r s ) 2 , so that typically, very high jet powers would be needed to account for the observed VHE emission.
In principle, however, the dynamics and the geometry of the emission region(s) could be much more complex. Salvati et al. [19], for example, have argued that if the emission is due to a series of conical shocks, fast blazar variability does not necessarily require extreme jet parameters. In addition the observational findings have triggered new conceptual developments in which rapid variability is related to, e.g., moving subregions or turbulence within the jet, black hole magnetospheric processes, or the propagation of non-linear electromagnetic waves (e.g., [20][21][22][23][24][25]). Major reference scenarios currently include black hole horizon gaps, jets-in-jet, or jet interaction models (e.g., [26][27][28]); see Figure 3 for an illustration. Each of these has its own challenge from, e.g., limited power extraction in the case of magnetospheric gaps and unusual magnetizations in the case of reconnection-induced mini-jets to the requirement of maximal jet powers in the case of jet interaction models (see [14] for a recent discussion). At the current stage, the most promising way to gain further insights into the physical origin of variability seems to be to go beyond classical minimum variability considerations and homogeneous one-zone approaches and to take into account the full timing characteristics of the emission, as outlined below. Methodologically, this draws on the general concept that the variability in AGN can be fairly well understood as a stochastic process. An observed time series (i.e., a light curve) is then seen as a realization of the underlying stochastic process that is sampled by the observation. 1 1 There are two common approaches to deal with irregular sampling in stochastic light curves, i.e., by means of extensive Monte Carlo simulations of artificial light curves [29,30] or by means of likelihood-based approaches using special parameterized stochastic models in the time domain such as, e.g., the first-order continuous-time autoregressive (Ornstein-Uhlenbeck) process (e.g., [31,32]). The latter offer an efficient means to extract information from large time-domain datasets, in particular the power spectral density (PSD). While autoregressive models are becoming increasingly popular (cf. [33,34]), general caveats concern the limitations of linear (e.g., ARMA, CARMA) models to extract physical information from (very probably) non-linear systems (e.g., [35]); see also Section 2.2. Particle acceleration in unscreened ("gap") electric fields close to the black hole horizon triggers an electron-positron pair cascade, which is accompanied by inverse Compton up-scattering of the ambient disk photon to the VHE regime. VHE variability then becomes a possible signature of jet formation [26]. Middle: Jets-in-jet model [27], where a variety of "mini-jet" features ("plasmoids") is produced by magnetic reconnection events within the main jet flow. This could lead to an additional velocity component relative to the main flow (Γ r ) and allow for a favorable orientation (i.e., increased Doppler boosting and time scale reduction) with respect to the observer. Right: Illustration of a hadronic model, where interactions of the jet with a small, massive obstacle (star or cloud) facilitate shock-acceleration and introduce a sufficient target density to allow for efficient pp-collisions [28].

PDF Shape and Log-Normality
Important insights can be obtained by characterizing the underlying process driving the variability by its observed probability density function (PDF). The PDF can in principle be estimated by fitting a model function to the histogrammed data. As a simple comparison, one can study the distribution of VHE fluxes (e.g., the number of runs as a function of the measured photon flux or the logarithm thereof) and evaluate the appropriateness of a fit by a Gaussian distribution as representative of a normal random process. Figures 4 and 5 provide examples in the case of PKS 21455-304 and Mkn 501, where a clear preference for the logarithm of the flux to be Gaussian (normal) distributed, i.e., for log-normality, has been found [36,37]. 2 Similar results are obtained by studying the relationship between the mean flux and the absolute rms (root mean square) variability (e.g., [38,39]), 3 where a linear relationship is known to provide evidence for a log-normal distribution of the fluxes [40,41]. 2 A random variable X is log-normally distributed if log(x) obeys a normal (Gaussian) distribution. The PDF of such a random variable is of the form f (x) = 1 where µ is the mean and σ the standard deviation. 3 Defined as the square-root of the light curve variance and related to the square-root of the integral of the PSD, see below, over the observable frequency range, i.e., PSD(ν) · dν 1/2 .  is preferred over a Gaussian. The difference is not as significant, though (see also [42]), and one could also speculate whether the sum of two Gaussians might provide a better characterization of the log(flux)-distribution. From [37].
Where statistics permits, current evidence suggests that the gamma-ray fluxes in blazars are preferentially log-normally distributed (cf. also [43] in the case of Mkn 421 and [44] for the case of bright Fermi blazars). Lognormal flux variability is not unusual, but in fact is well-known for accreting galactic sources such as X-ray binaries [40], where it has been linked to the underlying accretion process. The current findings are of interest since in a lognormal process, the fluctuations of the flux are on average proportional, or at least correlated, to the flux itself. This rules out additive processes for the origin of the observed variability and instead favors multiplicative ones. 4 If this is true, then additive models (e.g., shot-noise or a simple superposition of many "mini-jets") are no longer adequate to describe the observed variability behavior, and multiplicative, cascade-like scenarios need to be invoked. 4 For a stationary stochastic process X that results from a multiplication of N random subprocesses x i , X = ∏ x i , the logarithm of X is equivalent to the sum of the logarithm of the individual x i , i.e., log X = log x 1 + log x 2 + ... + log x N . By the central limit theorem, this sum must approach a normal (Gaussian) distribution for N → ∞. Astrophysically, N does not have to be large to achieve a good log-normal distribution [45].
In principle, several scenarios for the origin of log-normal flux distributions could be envisaged: (i) Similar as for X-binaries, the inferred log-normality at gamma-ray energies in blazars could mark the influence of the accretion disk on the jet (e.g., [46,47]). Independent density fluctuations in the accretion disk on local viscous timescale, i.e., t v (r) ∼ (1/α) (r/h) 2 (r/r g ) 3/2 r g /c in the case of a standard disk (where r g := GM BH /c 2 , h is the disk height, and α the viscosity parameter), provide one possible realization. Provided damping is negligible, these fluctuations can propagate inward and couple together to produce a multiplicative (power-law noise) behavior in the innermost disk part [48][49][50]. If this behavior is efficiently transmitted to the jet (particle injection rate), the gamma-ray emission could be modulated accordingly [18]; see Figure 6. This requires, amongst others, that the timescale for particle acceleration and radiative losses within the jet be correspondingly small [18]. Log-normality at different wavelengths and over long timescales might then well be feasible (i.e., up to timescales t v (r d ) where r d is the characteristic disk radius). Short timescales (on horizon light-crossing times), however, could only be achieved in the case of a thick (inner) disk (where h ∼ r); see also below. Figure 6. A possible accretion-disk origin for log-normality. Density fluctuations in the disk, arising on local viscous timescales, can propagate inwards and couple together so as to produce a multiplicative behavior in the accretion rate. If this is transmitted to the jet and picked up by the acceleration mechanism, the resultant VHE emission might share in the underlying log-normal characteristics.
(ii) Cascade-related emission processes (e.g., proton-induced synchrotron cascades, cf. [51], or magnetospheric inverse-Compton pair production cascades, cf. [26]) might possibly lead to some log-normal flux distributions. Relevant constraints in these cases arise, however, from the energy bands in which log-normality has been seen (e.g., in the optical, X-ray, and gamma-ray range for PKS 2155-304 [38]) and the timescales over which log-normality has been found (i.e., from sub-hour to yearly timescales at gamma-ray energies). The latter are expected to be limited by the gap travel time for magnetospheric processes and the dynamical or escape time for hadronic cascades.
(iii) Alternatively, log-normality could be associated with the acceleration process itself, e.g., with random fluctuations in the particle acceleration rate [52]. In the case of diffusive shock acceleration, for example, the accelerated particle distribution contains powers in t acc , i.e., n(γ) ∝ t acc γ −1−t acc /t esc (1 − γ/γ max ) t acc /t esc −1 (1/γ 0 − 1/γ max ) −t acc /t esc , where t acc ∼ κ/u 2 s denotes the acceleration and t esc the escape timescale, respectively [53], γ 0 is the injection Lorentz factor, u s is the shock speed, and γ max is the maximum Lorentz factor determined by radiative losses. If diffusion κ is characterized by Gaussian perturbations such that t acc = t acc,0 + ∆t acc , then the fractional variability ∆n(γ)/n(γ) becomes a linear combination of Gaussian and log(γ)-terms [52]. Depending on the energy scale, the accelerated particle number density could thus resemble a lognormal distribution. Figure 7 provides an illustration of this effect. Log-normality in the radiating particle distribution could then well result in a log-normal flux distribution in the case of synchrotron or external inverse Compton. Log-normality in this case, however, would be energy-dependent, with its significance becoming weaker towards lower energies and disappearing for energies close to a threshold (γ → γ 0 ). PDF studies at different frequencies (cf. [54] for PKS 2155-304) could be of help to distinguish between the noted scenarios. In principle, the above does not preclude the possibility that log-normality in different states (long term, low level, and short term, flare) arises due to different processes. Figure 7. Log-normality in the accelerated particle number density arising from random fluctuations in the particle acceleration rate at a shock front. The figure shows the histogram of the particle number density at γ = 10 3 based on simulations in which Gaussian perturbations are introduced into t acc . The dashed green line represents the best-fitting Gaussian and the solid blue line the best-fitting log-normal PDF. Particle injection with γ 0 = 10 has been assumed. From [52].

PSD and Power-Law Noise
In time series analysis, the power spectral density (PSD) represents an important tool for characterizing variability. The PSD provides a measure for the contribution of different timescales to the variability power, i.e., it quantifies the amount of variability power as a function of (temporal) frequency (ν ∼ 1/t). Methodologically, this implies a move from a signal description in the time domain to one in the frequency domain by means of Fourier transformation. In practice, we can only get an estimate of the PSD, and the simplest method is to estimate the PSD by the periodogram (i.e., the squared modulus of the discrete Fourier Transform, properly normalized). Accordingly, a periodicity at a given timescale would lead to a peak in the PSD at the corresponding frequency, while a continuum power spectrum would correspond to non-periodic signals. The latter are typically parameterized as power-law noise with a PSD of the form: where β ≥ 0 is the power index (e.g., β = 0 for white noise, β = 1 for pink or flicker noise, and β = 2 for red or Brownian noise). A broken power-law PSD would then be indicative of a characteristic timescale (t b ) at the corresponding break frequency (ν b ∼ 1/t b ). PSD analysis has by now been performed for a variety of gamma-ray-emitting AGN, particularly in the high-energy Fermi-LAT domain [31,[55][56][57][58]. In the VHE domain, a particularly interesting result relates to PKS 2155-304, where an extended PSD analysis revealed a flicker noise behavior β ∼ 1 for its long-term low-level VHE activity over frequencies corresponding to timescales ≥ 1 d (and with a power index similar to the one seen at HE on timescales ≥ 10 d), cf.  In the simplest interpretation, the PSD slope is considered to be stationary and to follow a power-law with a transition from β ∼ 1 to β ∼ 2 at around ν b ∼ 1/t b , where t b ∼ 1 d; see Figure 9. This would then be comparable to the X-ray PSDs of Seyfert AGNs, which reveal a similar steepening by ∆β 1 at their break frequencies [59]. Interestingly, at X-ray energies, a break time of t b ∼ 1 d has been suggested earlier for PKS 2155-304 [60]. In principle, accretion-disk fluctuations (see above) could also account for a power-law noise behavior [48]. In the case of the radio-quiet Seyfert AGNs, McHardy et al. [59] have reported a simple quantitative relationship between the break time t b -then associated with the characteristic timescale at the inner edge of the accretion disk-the observed (bolometric) luminosity L B , and the black hole mass M BH of the source. Although PKS 2155-304 is not a radio-quiet object, one could speculate that a similar relation applies if the timing properties originate in the accretion flow [39]. For a standard disk (with the accretion rate as a proxy for L B ), this scaling relation becomes (t b /1 d) 0.7 (M BH /10 8 M ) 1.12 /ṁ 0.98 E , where the accretion rateṁ E has been expressed in units of the Eddington rate. Hence, if this relation also applies to the supposed VHE break time of PKS 2155-304 (t b ∼ 1 d), accretion rates close to Eddington would be implied even for the quiescent VHE state. It seems thus more likely that the break time in PKS 2155-304 is related to a change in accretion flow conditions such as a transition from an inner advection-dominated to an outer standard disk configuration [39].
Obviously, Fourier analysis of non-thermal emission models could allow one to gain further insights into the origin of the gamma-ray variability in blazars. This particularly includes recent explorations as to the possible modifications of an underlying (injected) PSD shape by radiation [61,62]. This can be done by starting from the appropriate time-dependent particle transport equation (including the relevant acceleration and radiative loss terms) for N e (γ, t), the number of electrons between γ and γ + dγ at time t. One can then look for solutionsÑ e (γ, f ) of the Fourier-transformed equation assuming injection of power-law noise withQ(γ, f ) ∝ f −β . This can be used to, e.g., evaluate the impact on the expected synchrotron, external Compton (EC), or synchrotron self-Compton (SSC) fluxF (withF as the Fourier transform of the respective flux). Estimating the PSD from the periodogram, one, e.g., finds that at sufficiently low frequencies |F(SSC)( f )| 2 ∝ f −(4β−2) for SSC and |F(EC)( f )| 2 ∝ f −2β in the case of EC [61]. Since the gamma-ray emission in flat spectrum radio quasars (FSRQs) and BL Lacs is thought to be dominated by different radiation processes (EC versus SSC), one might expect differences in the PSD slopes for these sources if the injection (β) would be similar. In particular, if the above relation is applied to the BL Lac object PKS 2155-304, believed to be dominated by SSC, this would suggest injection with β ∼ 1 (for the flare) and β ∼ 0.75 (for the long-term, low-level activity), respectively. The former could possibly be related to the flicker-noise behavior induced by accretion disk fluctuations, while the latter would require some radial dependencies [48]. At optical (R-band) frequencies, the emission of PKS 2155-304 has been reported to obey a PSD with slope ∼ (1.8-2.4) on timescales larger than several days [55] (but see also [54]). In the above approach, the PSD for synchrotron emission roughly follows f −2β and steepens towards higher frequencies ( f > ∼ 1/t esc ) if escape is included [61]. Given current uncertainties in PSD slope determination, it seems thus possible to relate the different slopes at optical and HE/VHE to synchrotron and SSC radiation processes, respectively (cf. also [54]). More advanced modeling in this regard, however, seems desirable.
From a formal point of view (and also relevant to the below), it seems important to keep in mind that parameter estimation based on light curve simulations need to be consistent in as much as they account for both the details of the PDF (e.g., log-normality) and the PSD. 5

Quasi-Periodic Variability
With the availability of continuously-sampled gamma-ray light curves, periodicity analysis has gained new momentum as a possible tool to probe deeper into the astrophysical nature of the sources. In particular, the (unbiased) Fermi-LAT light curves of blazars have been used to search for quasi-periodic oscillations (QPOs) on yearly timescales. A prominent example is provided by the BL Lac objects PG 1553+113 (at a redshift of z ∼ 0.5), where a periodicity of P = 2.2 y in HE gamma-rays (and similarly, also in the optical) over ∼ 9 y has been reported [63][64][65]; see Figure 10 for an illustration. Figure 10. The ∼ 9 y-long light curve of PG 1553+113 in HE gamma-rays (20-d bins; green points) and in the optical (blue points) band. The black line shows the periodicity model, and the vertical line corresponds to some spectrum calculations. From [64]. 5 The widely employed Timmer and König algorithm [29] to simulate light curves assumes a normal (Gaussian) stochastic process, and thus potentially introduces errors in parameter estimation for log-normal processes [30]. It remains to be studied how this affects the results given the quality of available gamma-ray data.
Multi-year quasi-periodic variability in AGN light curves has often been interpreted in the context of supermassive binary black holes (SBBHs), such as the optical 12-y periodicity in the BL Lac object OJ 287 (z = 0.3) (e.g., [72]) or the apparent optical P 5.2-y periodicity in the quasar PG 1302-102 (z = 0.28) [73]. In general, SBBHs are a natural stage of hierarchical galaxy formation in which elliptical galaxies (e.g., the host galaxies of radio-loud AGN) are formed by galaxy mergers [74]. Depending on their evolutionary path, close SBBHs (with separation less than a few parsecs) could potentially still reside in the center of (some) radio-loud AGN (e.g., [75][76][77]). This is in fact supported by recent morphological evidence for SBBH-driven geodetic precession (P ∼ 10 6 -10 7 y) in the radio maps of powerful AGN [78]. Close, accreting SBBH systems are likely to be accompanied by a circumbinary disk surrounding the binary, with circumbinary gas streams feeding mini-disks around each black hole and affecting their evolution (e.g., [79,80]).
The HE findings noted above have led to the emergence of new SBBH interpretations. From a statistical point of view, however, some caution needs to be exercised concerning the significance of the inferred periods, as in most cases (perhaps apart from Mkn 501), only a few cycles are present over the available data. As shown by Vaughan et al. [81], clear phantom periodicities over ∼(2-3) cycles can be well found in pure noise data. Moreover, there are indications that the QPO results may be more dependent on the analysis method employed than initially thought [82]. In terms of a global significance (assuming the absence of other physical reasoning) (cf., e.g., [83]), year-type QPOs are often no longer significant [84]. Increasing the number of possible cycles by continuous (well-sampled) monitoring will be important to better assess the current evidence and qualify its impact. In addition, multi-wavelength comparison (as partly done for some of the sources) could improve the robustness of possible detections, though (physically) periods may not necessarily have to be coinciding for very different wavebands. The HE Fermi-LAT γ-ray observations could in principle be complemented by observations in the VHE domain, though gaps and uneven sampling may represent a more serious issue. Combining different VHE instruments such as FACT and HAWC (e.g., see [85] for the case of Mkn 501) could however improve the situation.
While speculative, one could still try to evaluate possible physical mechanisms for periodicity in AGN. While SBBHs are expected, they may often not provide the best explanation for the noted QPOs: QPOs with physical periods P r = P/(1 + z) ∼ 1 y would imply very close binary systems and thus a significant amount of gravitational wave emission. This would result in a very short gravitational lifetime of the system, T ∼ (10 3 -10 4 ) y, making it unlikely that we should be able to detect many of such systems. Similarly, Pulsar Timing Array observations impose upper limits on the (nano-hertz) gravitational stochastic background. If the binary hypothesis is tested for the HE blazar population, assuming year-like QPOs and using the respective luminosity functions, only a very minor fraction (0.01-0.1%) of BL Lacs and FSRQs could harbor SBBHs [86]. In fact, pervious results are rather suggestive of orbital periods of the order of ∼ 10 y [76]. It seems thus more promising to relate year-type QPOs to other origins such as the helical motion of a component inside a rotating jet, where differential Doppler boosting accounts for a periodic lighthouse effects [87]; or some quasi-periodic modulation in the accretion flow (e.g., induced by modulations of the transition from an ADAF to a standard disk) feeding the jet (e.g., [88]). Constraints on the jet radius suggest possible periods P < 2 y for the former scenario (cf., [89]), while the latter would suggest transition radii r tr ∼ 100r g . At the current stage, it seems important to evaluate the adequacy of such explanations in more detail.
We note that a very interesting, month-type (34.5 d) HE QPO over six cycles has been recently reported for the BL Lac object PKS 2247-131 (z = 0.22) [90]. The QPO behavior follows an outburst of the source in October 2016 and is thought to be associated with a helical jet structure induced by the orbital motion in an SBBH. The observed month-type HE QPO could in principle results from a shortening of the real physical period due to relativistic travel-time effects (e.g., [89]), in which case it would be suggestive of a (physical) orbital period of P 7.1 y [90]. Interestingly, a similar month-type (∼ 23 d) behavior (over 4-5 cycles) at VHE energies has been reported for Mkn 501 during a large flare in 1997 [91,92] and has also been interpreted in an SBBH framework [93].

A Possible Example: PKS 2155-304
Given the above-noted VHE timing characteristics of PKS 2155-304 (log-normality and power-law noise behavior) and the properties of disk fluctuations, one could hypothesize that the VHE variability is driven by density fluctuations in the accretion-disk. If these are efficiently transmitted to the jet, it could lead to a respective, power-law noise spectrum in injection for Fermi-type particle acceleration. One would only be able, however, to observe rapid VHE variations (occurring on timescales as short as ∆t obs 200 s) with the corresponding disk characteristics if these signatures do not get obscured by other processes occurring on longer timescale within the source. As the flux changes seen by an observer will appear to be convolved and thus dominated by the longest timescale, this would require that the (observed) timescales for photons traveling across the radial width of the source and for the relevant radiative loss processes still remain smaller than ∆t obs . As shown in [18], this seems feasible in the case of PKS 2155-304. To account for the observed minimum variability, however, the black hole mass would need to be limited to m BH ≤ 4 × 10 7 M (assuming a thick inner disk); see Equation (1). In particular, the scenario would require a black hole mass that is seemingly smaller than the total central black hole mass inferred from the M BH − L bulge relation. A straightforward way to account for this is to consider a binary black hole system where the jet dominating VHE emission is emitted from the secondary (less massive) black hole [18]; see Figure 11. As noted above, viewing elliptical galaxies as a merger product of spiral galaxies, an SBBH stage is expected to occur at some time. Once the secondary becomes embedded in the outer disk around the primary (more massive) black hole, it will start clearing up an annular gap. Mass supply from the circumbinary disk to the central binary, however, can still continue through tidal, time-dependent (periodically-modulated) gas streams, which penetrate the gap and preferentially feed the secondary (disk) (e.g., [79,94,95]). Numerical simulations in particular reveal that a reversal of mass accretion rates can occur (i.e.,ṁ BH,2 >ṁ BH,1 despite m BH,2 < m BH,1 ), acting towards equal-mass binaries and suggesting that the secondary and its jet can become more luminous than the primary. In such a scenario, a high VHE luminosity, as inferred for PKS 2155-304, could be accounted for, despite the smaller (secondary) black hole mass. Figure 11. A possible binary black hole model for PKS 2155-304, with the binary (black holes plus mini-disks) being surrounded by a circumbinary disk. The tidal gas stream from this disk preferentially feeds the less massive, secondary black hole, whose jet emission dominates the observed VHE radiation spectrum. Both the short VHE variability and the high (apparent) luminosity can thus be accommodated, cf. [18].
When put in context, these considerations suggest that care needs to be exercised in the interpretation of VHE variability constraints to avoid a spurious identification of the emitting (possibly secondary) black hole mass with the total (possibly primary + secondary) central black hole mass of the system. In the case of PKS 2155-305, there are additional pieces of circumstantial evidence, such as a (comparable) small black hole mass inferred from its X-ray variability (PSD) properties (e.g., [96]), indications for optical long-term (∼ 7 y) periodicity [97], and jet bending on pc-scales (e.g., [98]), which may in fact be more easily integrated within a binary framework.

Conclusions
Time variability provides key information about the physics of astrophysical objects beyond that given solely in their energy spectra. As indicated in this paper, this information is vital to properly understand the phenomena and processes involved. Given the progress in gamma-ray astronomy, statistical tools can now be applied to study the timing properties of AGN and to characterize their variability in terms of, e.g., flux distributions (PDF) and (temporal) frequency dependence (PSD), along with characteristic timescales and possible QPOs. This positive development is set to strengthen with upcoming instrumental developments and will make it possible to deepen our understanding of the link between accretion dynamics, black hole physics, and jet ejection.