Fourier Analysis of Blazar Variability
aa r X i v : . [ a s t r o - ph . GA ] J un ApJ, accepted; June 20, 2018
Preprint typeset using L A TEX style emulateapj v. 12/16/11
FOURIER ANALYSIS OF BLAZAR VARIABILITY
Justin D. Finke and Peter A. Becker U.S. Naval Research Laboratory, Code 7653, 4555 Overlook Ave. SW, Washington, DC, 20375-5352 School of Physics, Astronomy, and Computational Sciences, MS 5C3, George Mason University, 4400 University Drive, Fairfax, VA
ApJ, accepted; June 20, 2018
ABSTRACTBlazars display strong variability on multiple timescales and in multiple radiation bands. Theirvariability is often characterized by power spectral densities (PSDs) and time lags plotted as functionsof the Fourier frequency. We develop a new theoretical model based on the analysis of the electrontransport (continuity) equation, carried out in the Fourier domain. The continuity equation includeselectron cooling and escape, and a derivation of the emission properties includes light travel timeeffects associated with a radiating blob in a relativistic jet. The model successfully reproduces thegeneral shapes of the observed PSDs and predicts specific PSD and time lag behaviors associatedwith variability in the synchrotron, synchrotron self-Compton (SSC), and external Compton (EC)emission components, from sub-mm to γ -rays. We discuss applications to BL Lacertae objects and toflat-spectrum radio quasars (FSRQs), where there are hints that some of the predicted features havealready been observed. We also find that FSRQs should have steeper PSD power-law indices thanBL Lac objects at Fourier frequencies . − Hz, in qualitative agreement with previously reportedobservations by the
Fermi
Large Area Telescope.
Subject headings:
BL Lacertae objects: general — quasars: general — radiation mechanisms: non-thermal — galaxies: active — galaxies: jets INTRODUCTION
Blazars, active galactic nuclei (AGN) with jets movingat relativistic speeds aligned with our line of sight, arethe most plentiful of the identified point sources observedby the
Fermi
Large Area Telescope (LAT; Abdo et al.2010a; Nolan et al. 2012). Their spectral energy dis-tributions (SEDs) are dominated by two components.At lower frequencies there is a component produced viasynchrotron emission, peaking in the radio through X-rays. There is also a high energy component in γ -rays,most likely due to Compton-scattered emission. Theseed photons for Compton scattering could be the syn-chrotron photons themselves (synchrotron self-Comptonor SSC; Bloom & Marscher 1996), or they could be pro-duced outside the jet (external Compton or EC) bythe accretion disk (Dermer & Schlickeiser 1993, 2002),the broad line region (Sikora et al. 1994), or a dusttorus (Kataoka et al. 1999; B la˙zejowski et al. 2000). TheDoppler-boosting due to the combination of relativisticspeed and a small jet inclination angle amplifies the ob-served flux, shifting the emission to higher frequencies,and decreasing the variability timescale.The LAT monitors the entire sky in high-energy γ -rays every 3 hours, providing well-sampled light curvesof blazars on long time scales. For some sources, γ -raydata have been supplemented by high cadence observa-tions in the radio through very-high energy (VHE) γ -rays creating unprecedented light curves with few gapsin wavelength or time (e.g., Abdo et al. 2011b,c). Blazarvariability is often characterized by power spectral den-sities (PSDs; e.g., Abdo et al. 2010b; Chatterjee et al.2012; Hayashida et al. 2012; Nakagawa & Mori 2013;Sobolewska et al. 2014), which are essentially represen-tations of the Fourier transform without phase infor- justin.fi[email protected] mation. Although the LAT can provide long base-line, high cadence light curves, it has difficulty probingshort timescales for all but the brightest flares. How-ever, shorter-timescale variability may be observed withoptical, X-ray or VHE instruments (e.g., Zhang et al.1999; Kataoka et al. 2001; Zhang et al. 2002; Cui 2004;Aharonian et al. 2007; Rani et al. 2010). Less often,Fourier-frequency dependent time lags between two en-ergy channels are computed from light curves (e.g.,Zhang 2002). The PSDs of light curves at essentially allwavelengths resemble power-laws in frequency, S ( f ) ∝ f − b , with typically b ∼ −
3, usually steeper than PSDsfound from Seyfert galaxies (Kataoka et al. 2001). De-spite the popularity of PSDs for characterizing variabil-ity, their theoretical motivation has not been thoroughlyexplored (although see Mastichiadis et al. 2013).In this paper, our goal is to bridge the gap between the-ory and observations by exploiting a powerful new math-ematical approach for the modeling and interpretationobserved PSDs and time lags. Theoretically, the variabil-ity of blazars is often described by a continuity equation(e.g., Chiaberge & Ghisellini 1999; Li & Kusunose 2000;B¨ottcher & Chiang 2002; Chen et al. 2011, 2012). Thisequation describes the evolution of electrons in a compactregion of the jet, which is homogeneous by assumption.The electrons are injected as a function of time and en-ergy and the electron distribution evolves due to energyloss and escape. The expected electromagnetic emissionobservers might detect can be compared with observa-tions (e.g., B¨ottcher & Reimer 2004; Joshi & B¨ottcher2007). This can allow one to explore individual flares;however, in this paper, we take a different approach, bystudying the electron continuity equation in the
Fourierdomain . This allows the exploration of individual flares,as well as the study of aggregated long-timescale variabil-ity of sources using Fourier transform-related quantities Finke & Beckersuch as PSDs and phase and time lags. We focus hereon long-timescale variability, including multiple epochsof flaring and quiescence. A similar study applying thesame Fourier transform concept to the modeling of timelags in accreting Galactic black hole candidates has beenrecently carried out by Kroon & Becker (2014).We begin in Section 2, by defining the Fourier trans-form and its inverse, and various other functions usedthroughout the rest of the paper. In Section 3, we ex-plore analytic solutions to the continuity equation inthe Fourier domain, and present solutions in terms ofPSDs. In Section 4, we explore the solutions in terms oftime lags between different electron Lorentz factor “chan-nels” as a function of Fourier frequency. We explore theexpected synchrotron, SSC, and EC PSDs and Fourierfrequency-dependent time lags in Section 5, includinglight travel time effects due to the emitting region’s finitesize. We discuss applications to some PSDs and time lagsin the literature in Section 6, and conclude with a dis-cussion of our simplifying assumptions and observationalprospects in Section 7. Several of the detailed derivationsare relegated to the appendices. DEFINITIONS
Several definitions of the Fourier transform and associ-ated quantities are used throughout the literature. Herewe make the definitions used in this paper explicit. Fora real function x ( t ), we define the Fourier transform by˜ x ( f ) = Z ∞−∞ dt x ( t ) e πift = Z ∞−∞ dt x ( t ) e iωt , (1)where i = −
1. We will indicate the Fourier transformby a tilde, the Fourier frequency by f , and the angularfrequency by ω = 2 πf . We define the inverse Fouriertransform by x ( t ) = Z ∞−∞ df ˜ x ( f ) e − πift = 12 π Z ∞−∞ dω ˜ x ( ω ) e − iωt . (2)We define the PSD S ( f ) = | ˜ x ( f ) | = ˜ x ( f )˜ x ∗ ( f ) (3)where the asterisk indicates the complex conjugate. Wewill make use of the related representation of the Dirac δ -function δ ( t − t ) = Z ∞−∞ df e πif ( t − t ) . (4)We will use two versions of the Heaviside function, a“step” function defined by H ( x ) = (cid:26) x >
00 otherwise , (5)and the two-sided Heaviside “top-hat” function, H ( x ; a, b ) = (cid:26) a < x < b . (6)In Section 3.5 we will use the lower incomplete Gammafunction given by γ g ( a, x ) = Z x dy y a − e − y . (7) CONTINUITY EQUATION IN FOURIER SPACE
General Solution
The evolution of electrons in a nonthermal plasma“blob” can be described by a continuity equation givenby (e.g., Chiaberge & Ghisellini 1999; Li & Kusunose2000; B¨ottcher & Chiang 2002; Chen et al. 2011, 2012) ∂N e ∂t + ∂∂γ [ ˙ γ ( γ, t ) N e ( γ ; t )] + N e ( γ ; t ) t esc ( γ, t ) = Q ( γ, t ) , (8)where N e ( γ ; t ) dγ gives the number of electrons withLorentz factor between γ and γ + dγ at time t . Here˙ γ ( γ, t ) is the rate at which electrons lose or gain energy, t esc ( γ, t ) is the escape timescale, and Q ( γ, t ) is the rateat which electrons are injected in the jet. In the simplemodel presented here, we will assume that size of theplasma blob does not change with time, so that adia-batic losses can be neglected. We will also assume thatthe electron distribution in the blob is homogeneous andisotropic, and that variations occur throughout the blobsimultaneously. This is a common and useful assump-tion, although it is somewhat unphysical, as in orderfor the blob to be causally connected, variations cannotpropagate through the blob faster than the speed of light c . Assuming that ˙ γ and t esc are independent of t , we cantake the Fourier transform of both sides of the continuityequation leading to − πif ˜ N e ( γ, f ) + ∂∂γ [ ˙ γ ( γ ) ˜ N e ( γ, f )] + ˜ N e ( γ, f ) t esc ( γ )= ˜ Q ( γ, f ) (9)where ˜ Q ( γ, f ) is the Fourier transformed source term.This is a linear ordinary first-order differential equationwith a relatively simple solution. It is shown in AppendixA that if ˙ γ ≤ N e ( γ, f ) = 1 | ˙ γ ( γ ) | Z ∞ γ dγ ′ ˜ Q ( γ ′ , f ) × exp " − Z γ ′ γ dγ ′′ | ˙ γ ( γ ′′ ) | (cid:18) t esc ( γ ′′ ) − iω (cid:19) . (10)If t esc is independent of γ and cooling is from synchro-Compton processes, so that ˙ γ = − νγ , then one canperform the integral in the exponent, and γ ˜ N e ( γ, f ) = 1 ν exp (cid:20) − νγ (cid:18) t esc − iω (cid:19)(cid:21) × Z ∞ γ dγ ′ ˜ Q ( γ ′ , f ) exp (cid:20) νγ ′ (cid:18) t esc − iω (cid:19)(cid:21) . (11) Green’s Function Solution
Consider an instantaneous injection of monoenergeticelectrons with Lorentz factor γ at t = 0. Then ˜ Q ( γ, f ) = Q δ ( γ − γ ) in Equation (11) and one gets γ ˜ N e ( γ, f ) = Q ν exp (cid:20) − ν (cid:18) γ − γ (cid:19) (cid:18) t esc − iω (cid:19)(cid:21) × H ( γ − γ ) . (12)ourier Analysis of Blazar Variability 3The PSD is S ( γ, f ) = | γ ˜ N e ( γ, f ) | = (cid:20) Q ν (cid:21) exp (cid:20) − νt esc (cid:18) γ − γ (cid:19)(cid:21) H ( γ − γ ) . (13)The result is white noise for all electron Lorentz factors( γ ) and Fourier frequencies ( f ). Colored Noise
Since the PSDs of blazars resemble colored noise, andelectrons are generally thought to be injected as powerlaws in γ , one might expect that˜ Q ( γ, f ) = Q ( f /f ) − a/ γ − q × H ( γ ; γ , γ ) H ( f ; f , f ) (14)where f is some constant frequency and a ≥
0. Thatis, in the jet, shocks will occur randomly which ac-celerate and inject particles as a power-law distribu-tion in γ between γ and γ with index q . We willdeal only with frequencies in the range f ≤ f ≤ f .These limits are needed for the PSD to be normal-ized to a finite value. Frequencies greater than the in-verse of the blob’s light crossing timescale are partic-ularly unphysical, although we allow this for two rea-sons. First, it allows us to compare with other theoret-ical studies that allow variations faster than the lightcrossing timescale (e.g., Chiaberge & Ghisellini 1999;Zacharias & Schlickeiser 2013). Second, our blob is al-ready unphysical, since we allow variations through-out the blob simultaneously in the blob’s comovingframe. The normalization constant is related to the time-averaged power injected in electrons h L inj i over a timeinterval ∆ t by Q = 2 π ∆ t h L inj i m e c G p I r + I i − I r I + I . (15)A derivation of this equation and definitions of the quan-tities G , I r , I i , and I can be found in Appendix B. With˜ Q ( γ, f ), given by Equation (14), Equation (11) can berewritten as γ ˜ N e ( γ, f ) = Q ( f /f ) − a/ exp (cid:20) − νγ (cid:18) t esc − iω (cid:19)(cid:21) ν q − × (cid:18) t esc − iω (cid:19) − q Z u max u min du u q − e u (16)where u min = 1 νγ (cid:18) t esc − iω (cid:19) (17)and u max = 1 ν max( γ, γ ) (cid:18) t esc − iω (cid:19) . (18) Electron Injection Index q = 2It is instructive to look at the case where q = 2. In thiscase, the remaining integral in Equation (16) can easily -7 -6 -5 -4 -3 Fourier Frequency [Hz]10 -9 -8 -7 -6 -5 -4 -3 -2 -1 PS D [ H z - ] γ = 78 γ = 170 γ = 10 γ = 10 Fig. 1.—
The electron PSD from Equation (20) resulting from aninstantaneous flash ( a = 0) of electrons injected with a power-lawenergy index q = 2. Here we set t esc = 10 s, ν = 3 . × − s − , h L inj i = 10 erg s − , ∆ t = 1 year, γ = 10 , γ = 10 . Dashedlines indicate f = t − for each curve, and the dotted line indicates f = (2 πt esc ) − . be performed analytically. Then γ ˜ N e ( γ, f ) = Q ( f /f ) − a/ /t esc − iω exp (cid:20) − νγ (cid:18) t esc − iω (cid:19)(cid:21) × [ e u max − e u min ] , (19)and the PSD is S ( γ, f ) = | γ ˜ N e ( γ, f ) | = exp (cid:20) − νγt esc (cid:21) Q ( f /f ) − a (cid:18) t + ω (cid:19) − × ( exp (cid:20) νγ t esc (cid:21) + exp (cid:20) ν max( γ, γ ) t esc (cid:21) − (cid:20) νt esc (cid:18) γ, γ ) + 1 γ (cid:19)(cid:21) × cos (cid:20) ων (cid:18) γ, γ ) − γ (cid:19)(cid:21)) . (20)We identify asymptotes for the PSD for q = 2, Equa-tion (20). For these asymptotes we assume γ ≪ γ .1. If 1 / ( νt esc ) ≪ γ , and 2 πf /ν ≪ γ then S ( γ, f ) ≈ Q ( f /f ) − a ν γ . (21)2. If 1 / ( νt esc ) ≪ γ ≪ πf /ν then S ( γ, f ) ≈ Q ( f /f ) − a − f π sin (cid:18) πfνγ (cid:19) . (22)3. If γ ≪ / ( νt esc ) then S ( γ, f ) ≈ Q ( f /f ) − a /t + (2 πf ) exp (cid:20) − t esc (cid:18) νγ − t cool (cid:19)(cid:21) (23) × (cid:26) − exp (cid:20) − t cool t esc (cid:21) cos[2 πf t cool ] (cid:27) Finke & Beckerwhere we define t − = ν max( γ, γ ).The electron PSD resulting from Equation (20) is plot-ted in Figure 1 for parameters described in the cap-tion, which are fairly standard ones for FSRQs. We use a = 0 here, which represents an instantaneous injectionof power-law particles at t = 0, to more easily display theobservable features, a number of which are present. Forthe γ = 78, 170, and 10 curves, where γ ≪ ( νt esc ) − ),a break in the power law from S ( γ, f ) ∝ f − a to S ( γ, f ) ∝ f − ( a +2) is apparent, and the break frequency is at f ≈ (2 πt esc ) − . This is in agreement with asymptote 3 above. In gen-eral, by examining asymptote 3, it is clear that for low γ a break in the PSD will be found at a frequency of f = (2 πt esc ) − . Since the PSD measures periodic vari-ability, this indicates that for low γ , periodic variabilityon timescales less than the escape timescale is less pre-ferred. This is because electrons will always be escapingat a single timescale, which does not vary with time. Onecan also see in the γ = 170 curve in Figure 1 at high f structure related to the cosine seen in asymptote 3, withlocal minima at integer multiples of t cool .In the γ = 10 and γ = 10 curves, at high γ ( γ ≫ ( νt esc ) − ) the PSD will transition from S ( γ, f ) ∝ f − a to S ( γ, f ) ∝ f − ( a +2) , but in this case the transition is at f = t − , in agreement with asymptotes 1 and 2. Thus, variabilityon timescales less than the cooling timescales will be lessperiodic, since cooling on those smaller timescales willalways be present. It is also clear that at high f min-ima from the sin term in asymptote 2 occur at integermultiples of t − . Since at high values of γ , where thecooling is strongest, cooling on timescales t cool will beimmediate, and so periodic variability on timescales thatare integer multiples of t cool will also be strongly avoided. Electron Injection Index q = 2In this case, there is no simple analytic solution toequation (16), although it can be written with the in-complete Gamma function as γ ˜ N e ( γ, f ) = Q ( f /f ) − a/ exp (cid:20) − νγ (cid:18) t esc − iω (cid:19)(cid:21) × ν q − (cid:18) iω − t esc (cid:19) − q × [ γ g ( q − , − u max ) − γ g ( q − , − u min )] . (24)We compute the function numerically, and results can beseen in Figure 2 for q = 2 .
5. This confirms the features -7 -6 -5 -4 -3 Fourier Frequency [Hz]10 -9 -8 -7 -6 -5 -4 -3 -2 -1 PS D [ H z - ] γ = 78 γ = 170 γ = 10 γ = 10 Fig. 2.—
The same as Fig. 1, except that we set the injectionpower-law index q = 2 . seen in the q = 2 case are also seen for other values of q ,although the minima at high γ are not as pronounced. TIME LAGS
The phase lag between two electron Lorentz factor“channels”, γ a and γ b as a function of Fourier frequency( f ) can be calculated from∆ φ ( γ a , γ b , f ) = arctan (cid:26) Y I ( γ a , γ b , f ) Y R ( γ a , γ b , f )] (cid:27) (25)where [ γ a ˜ N e ( γ a , f )][ γ b ˜ N e ( γ b , f )] ∗ = Y R ( γ a , γ b , f ) + iY I ( γ a , γ b , f ) . (26)This implies that Y R ( γ a , γ b , f ) = Re [ γ a ˜ N e ( γ a , f )] Re [ γ b ˜ N e ( γ b , f )]+ Im [ γ a ˜ N e ( γ a , f )] Im [ γ b ˜ N e ( γ b , f )] (27)and Y I ( γ a , γ b , f ) = Re [ γ b ˜ N e ( γ b , f )] Im [ γ a ˜ N e ( γ a , f )] − Re [ γ a ˜ N e ( γ a , f )] Im [ γ b ˜ N e ( γ b , f )] . (28)The time lag can be calculated from the phase lag,∆ T ( γ a , γ b , f ) = ∆ φ ( γ a , γ b , f )2 πf . (29)For our solution, Equation (19) the time lag is∆ T ( γ a , γ b , f ) = 12 πf arctan (cid:26) Z I ( γ a , γ b , f ) Z R ( γ a , γ b , f ) (cid:27) (30)ourier Analysis of Blazar Variability 5 -7 -6 -5 -4 -3 Fourier Frequency [Hz]-10 ti m e l a g [ s ec ] γ a = 100, γ b = 120 γ a = 170, γ b = 10 γ a = 170, γ b = 78 γ a = 10 , γ b = 10 Fig. 3.—
The electron time lags with various values of γ a and γ b . Parameters are the same as in Figure 1. where Z I ( γ a , γ b , f ) = exp (cid:20) − νt esc (cid:18) γ a + 1 γ b (cid:19)(cid:21) × ( exp (cid:20) νt esc γ (cid:21) sin (cid:20) ων (cid:18) γ a − γ b (cid:19)(cid:21) + exp (cid:20) νt esc (cid:18) γ , γ a ) + 1max( γ , γ b ) (cid:19)(cid:21) × sin " ων γ , γ b ) − γ b + 1 γ a − γ , γ b ) ! − exp (cid:20) νt esc (cid:18) γ , γ a ) + 1 γ (cid:19)(cid:21) × sin (cid:20) ων (cid:18) γ a − γ , γ a ) + 1 γ − γ b (cid:19)(cid:21) − exp (cid:20) νt esc (cid:18) γ + 1max( γ , γ b ) (cid:19)(cid:21) × sin (cid:20) ων (cid:18) γ a − γ + 1max( γ , γ b ) − γ b (cid:19)(cid:21) ) (31)and Z R ( γ a , γ b , f ) is the same as Z I ( γ a , γ b , f ) exceptwith cos in place of sin. Several examples of electrontime lags can be seen in Fig. 3.For f ≪ νγ a / (2 π ) and f ≪ νγ b / (2 π ) the time de-lay will be approximately independent of frequency, withvalue ∆ T ( γ a , γ b , f ) ≈ ν A I ( γ a , γ b ) A R ( γ a , γ b ) (32) where A I ( γ a , γ b , f ) = (cid:18) γ a − γ b (cid:19) exp (cid:20) − νt esc (cid:18) γ a + 1 γ b (cid:19)(cid:21) + 1 γ b exp (cid:20) − νt esc γ b (cid:21) − γ a exp (cid:20) − νt esc γ a (cid:21) (33)and A R ( γ a , γ b , f ) = 1 + exp (cid:20) − νt esc (cid:18) γ a + 1 γ b (cid:19)(cid:21) − exp (cid:20) − νt esc γ b (cid:21) − exp (cid:20) − νt esc γ a (cid:21) , (34)where we also assumed γ < γ a ≪ γ and γ < γ b ≪ γ . At these low values of f , the lags are positive for γ a > γ b indicating the smaller γ lags behind the larger γ . This is due to the fact that electrons with smaller γ will take longer to cool than those with larger γ . If also( νt esc ) − ≪ γ a and ( νt esc ) − ≪ γ b then∆ T ( γ a , γ b , f ) ≈ ν (cid:18) γ a − γ b (cid:19) . (35)An example of this can be seen in Fig. 3, with the γ a =10 , γ b = 10 curve. If γ a ≪ ( νt esc ) − ≪ γ b then∆ T ( γ a , γ b , f ) ≈ t esc ( − νt esc γ b − νt esc γ a exp (cid:20) − νt esc γ a (cid:21)) . (36)In Fig. 3, an example can be seen with the γ a = 170, γ b = 10 curve.If γ a ≪ ( νt esc ) − and γ b ≪ ( νt esc ) − then∆ T ( γ a , γ b , f ) →
0. In Fig. 3, the curve that most closelyapproximates this is the γ a = 100, γ b = 120 curve.At high f , The behavior is quite complex. By inspect-ing Equation (31), we can see that the important fre-quencies are those that make the sine or cosine terms goto 0. They will be the integer or half-integer multiplesof f = t − ,a = νγ a (37) f = t − ,b = νγ b (38) f = ( t cool ,a − t cool ,b ) − . (39)There are many local minima and maxima at the integeror half-integer multiples of these values. EMISSION AND LIGHT TRAVEL TIME EFFECTS
In the previous sections we have explored the PSD andtime lags for the electron distribution. However, what isobserved is the emission of these electrons, through syn-chrotron or Compton-scattering. We will assume thatthe emitting region is spherical with co-moving radius R ′ and homogeneous, containing a tangled magnetic field ofstrength B . The blob is moving with a relativistic speed βc ( c being the speed of light), giving it a bulk Lorentzfactor Γ = (1 − β ) − / . The blob is moving with an Finke & Beckerangle θ to the line of sight giving it a Doppler factor δ D = [Γ(1 − β cos( θ ))] − . Although we make the simplify-ing assumption the blob is homogeneous, with variationsin the electron distribution taking place throughout theblob simultaneously, since it has a finite size, photonswill reach the observer earlier from the closer part thanthe farther part, and thus one must integrate over thetime in the past, t ′ . This is similar to the “time slices”of Chiaberge & Ghisellini (1999). We note again thatFourier frequencies higher than the inverse of the lightcrossing timescale are unphysical in the simple model wepresent here. Synchrotron and External Compton
Taking the light travel time into account, in the δ -approximation the observed νF ν flux at observed energy ǫ (in units of the electron rest energy) from synchrotronor Compton scattering of an external isotropic monochro-matic radiation field as a function of the observer’s time t is F ǫ ( t ) = K (1 + z ) t lc δ D Z R ′ /c dt ′ N e (cid:18) γ ′ ; tδ D z − t ′ (cid:19) (40)where t lc = 2 R ′ (1 + z ) cδ D (41)is the light crossing time in the observer’s frame. Aderivation of the light travel time effect can be foundin Appendix C. For synchrotron emission, K = K sy = δ D πd L cσ T u B γ ′ sy , (42)and γ ′ = γ ′ sy = s ǫ (1 + z ) δ D ǫ B . (43)The Thomson cross-section is σ T = 6 . × − cm ,the Poynting flux energy density is u B = B / (8 π ), ǫ B = B/B c where B c = 4 . × G, the redshift of thesource is z , and the luminosity distance to the source is d L . For external Compton (EC) scattering, K = K EC = δ D πd L cσ T u γ ′ EC , (44) γ ′ = γ ′ EC = 1 δ D s ǫ (1 + z )2 ǫ (45)(Dermer & Menon 2009). Here the external radiationfield energy density and photon energy (in units of theelectron rest energy) are u and ǫ , respectively. This ap-proximation is valid in the Thomson regime, i.e., when γ ′ . (Γ ǫ ) − . Primed quantities refer to the frame co-moving with the emitting region. The cooling rate pa-rameter is ν = 43 m e c cσ T ( u B + Γ u ) (46)where we ignore the effects of SSC cooling. It is shown in Appendix D that the Fourier transformof Equation (40) is˜ F ǫ ( f ) = K (1 + z )2 πif t lc δ D ˜ N e (cid:18) γ ′ , (1 + z ) fδ D (cid:19)(cid:26) exp (cid:20) πif (1 + z ) R ′ cδ D (cid:21) − (cid:27) . (47)Equation (47) implies the PSD of the synchrotron orEC flux is S ( ǫ, f ) = | ˜ F ǫ ( f ) | = K (1 + z ) ( πf t lc δ D ) (cid:12)(cid:12)(cid:12)(cid:12) ˜ N e (cid:18) γ ′ , (1 + z ) fδ D (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) sin ( πf t lc ) . (48)If the emitting region is very compact, i.e., if R ′ ≪ cδ D (2 f (1 + z )) − , then S ( ǫ, f ) ≈ K (1 + z ) ( πf t lc δ D ) (cid:12)(cid:12)(cid:12)(cid:12) ˜ N e (cid:18) γ ′ , (1 + z ) fδ D (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) ( πf t lc ) S ( ǫ, f ) ≈ K (1 + z ) δ (cid:12)(cid:12)(cid:12)(cid:12) ˜ N e (cid:18) γ ′ , (1 + z ) fδ D (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) . (49)So in the case of a very compact emitting region, thelight travel time effects play no part in the PSD, as onewould expect.The observed flux PSDs from synchrotron and EC onewould expect are shown in Figure 4, calculated fromEquations (48) and (20). The parameter values are thesame as in Figure 1, with additional parameters givenin the caption. The seed photon source is assumed tobe Ly α photons, presumably from a broad-line region.Again, parameters are chosen to be consistent with thoseone would expect from an FSRQ. Synchrotron PSDs areshown for 12 Hz, 12 µ m (the WISE
W3 filter’s centralwavelength), and 0.648 µ m (the central wavelength ofthe Johnson R band). EC PSDs are shown for 0.1 and1.0 GeV, which are within the Fermi -LAT energy range.Frequencies lower than 10 Hz can be affected by syn-chrotron self-absorption, which is not considered in thispaper. X-rays are not shown as they are likely dominatedby SSC emission, which is considered in Section 5.2 be-low. Also note that in real FSRQs the 12 µ m PSD couldsuffer from contamination from dust torus emission (e.g.,Malmrose et al. 2011) and the R band could suffer fromcontamination from accretion disk emission, and we donot take either of these possibilities into account.Figure 4 shows many of the features seen in Figure1. For photons generated from electrons with low γ ′ thePSDs show a break from S ( ǫ, f ) ∝ f − a to S ( ǫ, f ) ∝ f − ( a +2) at approximately f = (2 πt esc ) − , as seen in the 10 Hz and 0.1 GeV PSDs. For syn-chrotron, if u B ≪ Γ u , as is usually the case for FSRQs,ourier Analysis of Blazar Variability 7this regime occurs when ν ≪ ν cr , sy = 10 Hz (cid:18) δ D Γ (cid:19) (cid:18) Γ30 (cid:19) − × (cid:18) u − erg cm − (cid:19) − (cid:18) t esc s (cid:19) − (cid:18) B (cid:19)
11 + z . (50)For EC, this regime occurs when m e c ǫ ≪ E cr , EC = 2 GeV (cid:18) δ D Γ (cid:19) (cid:18) Γ30 (cid:19) − × (cid:18) u − erg cm − (cid:19) − (cid:18) t esc (cid:19) − × (cid:18) ǫ × − (cid:19)
11 + z . (51)If u B ≫ Γ u then ν cr , sy = 5 × Hz (cid:18) δ D (cid:19) (cid:18) t esc s (cid:19) − × (cid:18) B (cid:19) −
11 + z (52)and E cr , EC = 1 TeV (cid:18) δ D (cid:19) (cid:18) t esc (cid:19) − (cid:18) B (cid:19) − × (cid:18) ǫ × − (cid:19)
11 + z . (53)For photons generated from electrons with high γ ′ ,( ν ≫ ν cr , sy or m e c ǫ ≫ E cr , EC for synchrotron or EC,respectively) minima are seen at integer multiples of f = t − , as seen in the 12 µ m, R band and 1.0 GeVPSDs, as well as a break of 2 at approximately f = t − .There is an additional feature in Figure 4 not seen in Fig-ure 1 related to the light travel timescale. In agreementwith Equation (48), sin minima can be seen at integermultiples of f = 1 /t lc which appear in the PSDs at allphoton energies. Additionally, there is a break of 2 at ap-proximately this frequency, in agreement with Equations(48) and (49).In Figure 5 we plot observed PSDs for the parameters a = 1. This demonstrates that our model can reproduceany color noise in a synchrotron or EC PSD for an appro-priate choice of a , and that the features described aboveare preserved for different values of a .For synchrotron or EC emission, in our simple model,the time delays will play no role in the Fourier frequency-dependent time lags (Section 4). It is easy enough to seeby inspecting Equation (47) that terms associated withthe light travel time will cancel when calculating timelags. However, one must be careful to shift the time lagand frequency into the observed frame by multiplying by(1 + z ) /δ D and δ D / (1 + z ), respectively. An example ofobserved time lags can be seen in Fig. 6 for synchrotronemission and EC emission. Synchrotron Self-Compton -5 -4 -3 Fourier Frequency [Hz]10 -22 -20 -18 -16 -14 -12 -10 -8 -6 f l ux PS D [ ( e r g c m - s - H z - ) ] Hz 12 µ m 0.648 µ m0.1 GeV 1.0 GeV Fig. 4.—
The flux PSD computed from Equations (48) and (20)using the same parameters as in Figure 1. Additional parametersare δ D = Γ = 30, B = 1 G, u = 10 − erg cm − , ǫ = 2 × − , R ′ = 10 cm, and z = 1. At this redshift with a cosmology( h, Ω m , Ω Λ ) = (0 . , . , . d L = 2 × cm. The observedphoton frequency, wavelength, or energy is shown, along with t − for each curve (dashed lines), (2 πt esc ) − (dotted line), and t − lc (dashed-dotted line), all computed in the observer’s frame. -5 -4 -3 Fourier Frequency [Hz]10 -34 -32 -30 -28 -26 -24 -22 -20 -18 f l ux PS D [ ( e r g c m - s - H z - ) ] Hz 12 µ m 0.648 µ m0.1 GeV 1.0 GeV Fig. 5.—
The same as Fig. 4 except that we set a = 1. -6 -5 -4 -3 -2 Fourier Frequency [Hz]-10000100020003000400050006000 ti m e l a g [ s ec ] µ m, 0.365 µ m12 µ m, 0.648 µ m Fig. 6.—
The observed time lags as a function of observed Fourierfrequency. Parameters are the same as in Figure 1, with parametersthe same as Figures 1 and 4. The 0.648 µ m (central wavelengthof R band) and 0.365 µ m (central wavelength of U band) are fromsynchrotron, while the 0.1, 1.0, and 10 GeV emission is from EC. Finke & BeckerThe synchrotron-producing electrons will also Comp-ton scatter the synchrotron radiation they produce, lead-ing to SSC emission. Again, assuming the blob is homo-geneous, and taking into account light travel time effects,we have F SSCǫ ( t ) = K SSC (1 + z ) t lc δ D Z R ′ /c dt ′ Z min[ ǫ ′ ,ǫ ′− ]0 dǫ ′ i ǫ ′ i × N e s ǫ ′ ǫ ′ i ; tδ D z − t ′ ! N e s ǫ ′ i ǫ B ; tδ D z − t ′ (54)where K SSC = δ D cσ R ′ u B πd L V ′ (cid:18) ǫǫ B (cid:19) / , (55) V ′ = 43 πR ′ (56)is the blob volume in the comoving frame, and ǫ ′ = (1 + z ) ǫδ D (57)(Dermer & Menon 2009). Following the same procedurefor synchrotron and EC, for the Fourier transform we get˜ F SSCǫ ( t ) = K SSC (1 + z )2 πit lc f δ D (cid:26) exp (cid:20) πif (1 + z ) R ′ cδ D (cid:21) − (cid:27) × Z ∞−∞ df ′ Z min[ ǫ ′ ,ǫ ′− ]0 dǫ ′ i ǫ ′ i × ˜ N e s ǫ ′ ǫ ′ i ; (1 + z ) fδ D − f ′ ! ˜ N e s ǫ ′ i ǫ B ; f ′ . (58)Note, however, that we ignore the effects of SSC cool-ing. In this case, ˙ γ ( γ ) would be dependent on the elec-tron distribution, γ N ( γ ; t ) (or γ ˜ N ( γ, f )) leading to anon-linear differential equation. This is treated in de-tail for the continuity equation by Schlickeiser (2009);Schlickeiser et al. (2010); and Zacharias & Schlickeiser(2010, 2012a,b). The SSC PSD is S SSC ( ǫ, f ) = | ˜ F SSCǫ ( t ) | = (cid:18) K SSC (1 + z )2 πt lc f δ D (cid:19) sin ( πf t lc ) | I ( ǫ, f ) | (59)where I ( ǫ, f ) = Z ∞−∞ df ′ Z min[ ǫ ′ ,ǫ ′− ]0 dǫ ′ i ǫ ′ i × ˜ N e s ǫ ′ ǫ ′ i ; (1 + z ) fδ D − f ′ ! ˜ N e s ǫ ′ i ǫ B ; f ′ . (60) At low frequencies, f ≪ ( πt lc ) − and f ′ ≪ ( πt ′ esc ) − , S SSC ( ǫ, f ) ≈ K SSC (1 + z ) δ ( Q t esc ) (cid:18) δ D ǫ B ǫ (1 + z ) (cid:19) × [ C ( a, ǫ min , ǫ max ) C ( a, f min , f max )] (61)where C ( a, ǫ min , ǫ max ) = Z ǫ max ǫ min dǫ ′ i ǫ ′ i × ( exp " − νt esc r ǫ ′ i ǫ ′ + r ǫ b ǫ ′ i ! − exp " − νt esc r ǫ ′ i ǫ ′ − exp (cid:20) − νt esc r ǫ b ǫ i (cid:21) + 1 ) , (62) C ( a, f min , f max ) = f a × (cid:26) ( f − a max − f − a min )(1 − a ) − a = 1ln( f max /f min ) a = 1 , (63) f min = max (cid:20) f , (1 + z ) fδ D − f (cid:21) , (64) f max = min (cid:20) f , (1 + z ) fδ D − f (cid:21) , (65) ǫ min = max[ γ ǫ B , ǫ ′ /γ ] , (66)and ǫ max = min[ ǫ ′ , ǫ ′− , γ ǫ B , ǫ ′ /γ ] . (67)If f > (1 + z ) f /δ D ≫ f , then at low frequencies theSSC PSD will go as S SSC ( ǫ, f ) ∝ f − (2 a − . Synchrotron and SSC PSDs are shown in Figure 7. Theparameters were chosen to be those one would expectfrom a high-peaked BL Lac object. At low frequency,they agree with the asymptote above. Here, the PSDsare flatter (the PSD power-law index is smaller) thanfor synchrotron or EC. At high frequency the behav-ior is complex, but the features from the light crossingtimescale are apparent. There are no features associatedwith the cooling timescale of the electrons which produceSSC photons, since the SSC emission for a particular ob-served energy will be produced by a broad range of elec-trons. This is in contrast to synchrotron or EC wherethe photons at a particular energy can be approximatedas being produced by a single electron Lorentz factor. APPLICATIONS
The PSDs of blazars are almost always power-laws, S ( ǫ, f ) ∝ f − b (i.e., colored noise), although there issometimes evidence that they deviate from this. As seenin Section 5, our model can reproduce this, since in ourourier Analysis of Blazar Variability 9 -5 -4 -3 Fourier Frequency [Hz]10 -16 -14 -12 -10 -8 -6 -4 -2 f l ux PS D [ ( e r g c m - s - H z - ) ] Hz 12 µ m 1 keV0.658 µ m1.0 TeV10 GeV Fig. 7.—
The synchrotron and SSC flux PSD computed fromEquation (59). Parameters are the same as in Figure 4, except u = 0 and z = 0 .
1, giving d L = 1 . × cm with a cosmologywhere ( h, Ω m , Ω Λ ) = (0 . , . , . t cool for the electrons that produce thosephotons is shown as the dashed lines. The dotted curve indicatesthe frequency (2 πt esc ) − and the dashed-dotted line indicates thefrequency t − lc , all computed in the observer’s frame. model at low frequency for synchrotron or EC the PSDgoes as S sy/EC ( ǫ, f ) ∝ f − a and for SSC S SSC ( ǫ, f ) ∝ f − (2 a − , where a is the parameter from Equation (14).For an appropriate choice of a , it can reproduce anypower-law PSD. At higher frequencies, our model pre-dicts features in PSDs that deviate from a strict singlepower-law. In this section, we explore some of the appli-cations of our model to observed PSDs from the litera-ture. The VHE Gamma-ray PSD of PKS 2155 − In the PSD measured by HESS from PKS 2155 − S ( ǫ, f ) ∝ f − , out to f & − Hz. There does appear to be aminimum feature at f ≈ . × − Hz. It is not clearif there is a break in the PSD at higher frequencies thatthis. If this feature is associated with the light crossingtimescale, then R ′ ≈ ct lc δ D z ) = 6 . × (cid:18) δ D (cid:19) cm . (68)Such a high Doppler factor is needed for these flares toavoid γγ attenuation (Begelman et al. 2008; Finke et al.2008). Since the γ -rays from this source are likely as-sociated with the SSC mechanism for this source, thisassociation in unambiguous. There are no PSD featuresassociated with a cooling timescale for SSC, as shown inSection 5.2. The X-ray Timing Properties of Mrk 421
Zhang (2002) constructed PSDs from
BeppoSAX dataon Mrk 421. They construct PSDs with data from twoenergy intervals: 0.1–2 keV, and 2–10 keV (see theirFigure 4). Both PSDs show a minimum feature at f ≈ . × − Hz. Although the feature is tentative,since it is only dependent on one point for each PSD,and there are no error bars, the fact that the feature isat the same frequency in both energy bands is a hint that the feature is associated with the light crossing timescale.In this case, t lc = 2 . × s = 6 . R ′ ≈ ct lc δ D z ) = 1 . × (cid:18) δ D (cid:19) cm . (69)Zhang (2002) also provides time lags as a function ofFourier frequency between the two energy bands, 0.1–2 keV and 2–10 keV (see their Figure 5), giving us anopportunity to compare them with the results of Sec-tion 4. Zhang (2002) finds time lags at f . − Hzthat are approximately constant at ∆ T ≈ s. Theapproximate independence of the lag with frequency im-plies that the lag is in the regime where ( νt esc ) − ≪ γ a and ( νt esc ) − ≪ γ b , and the lag ∆ T ′ can be approxi-mated by Equation (35). Note that this equation givesthe lag in the jet comoving frame; in the observer’s frame,∆ T = ∆ T ′ (1 + z ) /δ D . The X-rays for Mrk 421 are likelyproduced by synchrotron emission and the external en-ergy density ( u ) is likely to be negligible for a BL Lacobject, so one can combine Equation (35) with Equations(43) and (46) to get B = B c ( δ D ǫ a ǫ b ) / ( z ) m e c ( ǫ / b − ǫ / a )8 cσ T u Bc ∆ T ) / (70)where u Bc = B c / (8 π ). With z = 0 .
03 for Mrk 421 and m e c ǫ a = 0 . m e c ǫ b = 2 keV, B = 0 . (cid:18) δ D (cid:19) − / (cid:18) ∆ T s (cid:19) − / G . (71) The Gamma-Ray PSDs of FSRQs and BL LacObjects
Nakagawa & Mori (2013) used more than four years of
Fermi -LAT data to compute the PSD of 15 blazars. EachPSD is fit with either a single or broken power-law model.Their values of b from their fits, where S ( ǫ, f ) ∝ f − b fromthe single power-law fit or lower index from the brokenpower-law fit if that fit is statistically significant are givenin Table 1. We neglect the FSRQ S4 1030+61, which haspoor statistics. There is a clear separation between b forFSRQs and BL Lac objects. All BL Lac objects have b ≤ .
6, while all FSRQs except for PKS 1222+216 have b > .
7. PKS 1222+216 seems to be an outlier in termsof its PSD power-law index, although its b is still largerthan for any of the BL Lac objects. We also computethe value of a from our model (recall its definition inEquation [14]) needed to reproduce the values of b . ForFSRQs, presumably emitting by EC, this is just a = b EC , while for BL Lac objects, presumably emitting by SSC,this is a = b SSC + 22(recall Section 5.2). The mean values of a for FSRQs andBL Lac objects are within one standard deviation (S.D.)of each other. The values of b from Nakagawa & Mori(2013) are in agreement with our theory with values of a that cluster around a ∼ TABLE 1
Fermi -LAT PSD power-law indices( b ) from Nakagawa & Mori (2013)and the values of a from ourmodel needed to reproduce them. Object b a
FSRQs4C +28.07 0 . ± .
23 0.93PKS 0426 − a . ± .
47 1.16PKS 0454 −
234 0 . ± .
27 0.78PKS 0537 − . ± .
64 0.86PKS 1222+216 0 . ± .
21 0.653C 273 1 . ± .
27 1.303C 279 1 . ± .
35 1.23PKS 1510 −
089 1 . ± .
30 1.103C 454.3 1 . ± .
24 1.00PKS 2326 −
502 1 . ± .
44 1.26mean 1.01 1 . . ± .
44 1.22Mrk 421 0 . ± .
21 1.19PKS 2155 −
304 0 . ± .
33 1.29BL Lac 0 . ± .
47 1.21mean 0.49 1 . PKS 0426 −
380 and PKS 0537 − Based on the traditional classification, PKS 0537 − −
380 are classified as BL Lac ob-jects. However, based on the new classification byGhisellini et al. (2011) they are considered FSRQs (seealso Sbarrato et al. 2012). We use the more recent clas-sification. For a discussion see D’Ammando et al. (2013)for PKS 0537 −
441 and Tanaka et al. (2013) for PKS0426 − γ -ray emission over a long timescale. The PSDscould be computed from these light curves. One expectsthat if the γ -rays are emitted by EC, they would have thesame b as synchrotron ( b sy = b EC = a ). But if the γ -raysare produced by SSC emission, one would expect b to beless steep than for the γ -rays compared to the optical,with the relation between the SSC and synchrotron PSDindices given by b SSC = 2 b sy − a −
2. Of coursethis could be complicated by emission from a thermalaccretion disk unrelated to the jet emission.
The Gamma-Ray PSD of 3C 454.3
The PSD of 3C 454.3 from Nakagawa & Mori (2013)shows a break at frequency f brk = 1 . × − Hz, cor-responding to a timescale of 6 . × s = 190 hours= 7.9 days. Their broken power-law fit to the PSD,where S ( ǫ, f ) ∝ f − b at f < f brk and S ( ǫ, f ) ∝ f − b at f > f brk shows b ≈ b ≈
3. This is a break ofabout 2, which is what is expected from our theory. Thetimescale could correspond to the light crossing, cooling,or escape timescales. If one interprets it as the light- crossing timescale, t lc , then R ′ ≈ ct lc δ D z ) = 1 . × (cid:18) δ D (cid:19) cm . (72)Interpreting it as the cooling timescale, t cool , and assum-ing δ D = Γ, the external radiation field is u ≈ m e c cσ T Γ t ′ cool γ ′ = 9 . × − (cid:18) Γ30 (cid:19) − (cid:18) E
100 MeV (cid:19) − / × (cid:18) ǫ × − (cid:19) / erg cm − = 6 . × − (cid:18) Γ30 (cid:19) − (cid:18) E
100 MeV (cid:19) − / × (cid:18) ǫ × − (cid:19) / erg cm − (73)where E is the observed photon energy. The first line as-sumes the seed photon source is a dust torus with temper-ature 1000 K; the second assumes it is Ly α , presumablyfrom the broad line region. Both numbers give ratherlow values for u .However, interpreting the break as either the cool-ing timescale or the light crossing timescale is prob-lematic, since variations on timescales much shorterthan this have been observed from 3C 454.3, includingdecreases on much faster timescales (Ackermann et al.2010; Abdo et al. 2011a). But one could also associatethe break with the escape timescale for electrons in theblob, as shown in Section 5.1. This break will occur at f = (2 πt esc ) − , so if t lc = t esc , the break will still be at afrequency 2 π lower than the one related to the light cross-ing timescale. Furthermore, the escape timescale could inprinciple be longer than the light crossing timescale, sincemagnetic fields in the blob would curve the electron’spath and decrease the time it takes to escape. We notethat in Figure 4, the break in the 0.1 GeV PSD is indeedassociated with the escape timescale, f = (2 πt esc ) − ,showing that this is at least plausible. If the break in thePSD of 3C 454.3 (Nakagawa & Mori 2013) is due to elec-tron escape, then the escape timescale in the observer’sframe will be t esc = 7 . / (2 π ) = 30 hours, and in thecomoving frame, t ′ esc = 20 days (cid:18) δ D (cid:19) . (74)How could one distinguish between these interpreta-tions? One possibility would be to observe the PSDsat more than one waveband. If the break is due tothe light-crossing timescale, the break frequency shouldbe present independent of the waveband. The es-cape timescale break could also be independent of fre-quency if the escape timescale is energy independent,as it is in our model. The cooling timescale shouldbe energy-dependent, and thus the break frequency willbe different in different wavebands. For 3C 454.3, thelight-crossing timescale interpretation is disfavored sincesmaller timescale fluctuations are present in its lightcurve (e.g., Ackermann et al. 2010; Abdo et al. 2011a).ourier Analysis of Blazar Variability 11 Optical PSDs of Blazars
Chatterjee et al. (2012) compute R band PSDs for6 blazars based on about 200-250 days of continuousdata. Their PSDs have power-law indices are signifi-cantly steeper than those from the same objects’ γ -rayPSDs from Nakagawa & Mori (2013). The exception isPKS 1510 − b = 0 . +0 . − . , significantly flatter than the γ -ray PSD.Our theory predicts that synchrotron and EC emissionshould have the same PSD slopes if produced by the sameelectron energies, and all but one of their sources are FS-RQs, which are expected to emit γ -rays by EC. One pos-sible reason for the discrepancy could be the contamina-tion in the optical by the accretion disk. Another possi-bility is that the time intervals used by Chatterjee et al.(2012) are significantly shorter than the ones used byNakagawa & Mori (2013). As Chatterjee et al. (2012)point out, the large number of bright flares in theirtime interval for PKS 1510 −
089 could be the cause ofits especially flat R band PSD power-law index. Itcould also be that the different analysis methods usedby Chatterjee et al. (2012) and Nakagawa & Mori (2013)could lead to different results. Finally, it could be thatone of the assumptions of our theory is just not correct.The Kepler mission, with its excellent relative pho-tometry and short timescale sampling is well-suited formeasuring high-frequency PSDs. Wehrle et al. (2013) re-ported the
Kepler
PSDs of several radio-loud AGN, andfound no departure from a single power-law up to ∼ − Hz, above which white noise dominates. Edelson et al.(2013) explored the
Kepler
PSD of the BL Lac objectW2R1926+42 and found a “bending” power-law pro-vided a good fit to its PSD, with “bend frequency” corre-sponding to ≈ . Hz according to Edelson et al.(2013), making it an intermediate synchrotron peaked(ISP) object by the classification of Abdo et al. (2010c).However, its optical SED appears to be dominated byaccretion disk emission, implying its synchrotron peak isprobably at . . Hz, which would make it an LSP.The
Kepler light curve could have a contribution fromboth the thermal accretion disk emission and the non-thermal jet emission, making interpretation of its PSDdifficult. If the optical band is dominated by synchrotronemission, its status as an LSP implies that the elec-trons that produce its optical emission are in the regime γ ′ ≫ ( νt ′ esc ) − , meaning the “bend frequency” could beassociated with the light-crossing timescale or the cool-ing timescale. If it is associated with the light-crossingtimescale, the size of the emitting region is R ′ ≈ . × (cid:18) δ D (cid:19) cm . (75)If it is associated with the cooling timescale, the cool-ing is dominated by EC, and δ D = Γ, then u ≈ . × − (cid:18) Γ30 (cid:19) − / (cid:18) B (cid:19) / × (cid:18) λ obs (cid:19) − / erg cm − (76)where λ obs is the observed wavelength. If the cooling is dominated by synchrotron, then the cooling timescalecan be used to estimate the magnetic field, B ≈ . (cid:18) δ D (cid:19) − / (cid:18) λ obs (cid:19) − / G . (77) The X-ray PSD of 3C 273
An X-ray PSD of 3C 273 based on data combined from
RXTE , EXOSAT and other instruments was reported byMcHardy (2008). Similar to the γ -ray PSD of 3C 454.3,the PSD shows two power-laws with b ≈ . b ≈ . f brk ≈ . × − Hz. The break is close to 2,and the break frequency corresponds to a timescale of1 . × s ≈
280 hours ≈
12 days. The interpretationfor this break is more difficult than the γ -ray PSD for3C 454.3, since it is not clear whether the X-rays areproduced by SSC, EC, or even by a hot corona at thebase of the jet. If the X-ray emission is dominated byEC, the most likely interpretation of the break is withthe escape timescale, since the electrons generating theEC emission will almost certainly have Lorentz factors γ ≪ ( νt esc ) − . In this case the escape timescale in theobserver’s frame is t esc = 12 days / (2 π ) = 46 hours. Quasi-Periodic Oscillations
A number of Quasi-periodic oscillations (QPOs) havebeen reported in the X-ray and optical PSDs of blazars.These QPOs could be associated with the maxima athalf-integer values of either the light crossing timescale( t lc ) or the cooling timescale ( t cool ). See for exampleFigure 4, where maxima in the 0.648 µ m or 1.0 GeVPSDs could be confused with QPOs in noisy PSDs. Ifthis interpretation is correct, one could distinguish be-tween these possibilities ( t lc or t cool ) by observing thePSDs at more than one wavelength. If the QPO appearsin at the same frequency independent of wavelength, itwould argue for a t lc interpretation. If QPOs are foundat different frequencies at different wavelengths, it arguesfor the t cool interpretation. Assessing the significance ofQPOs in red noise PSDs can be subtle (e.g., Vaughan2005; Vaughan & Uttley 2006).As an example, we look at the claimed QPO in re-ported in XMM-Newton observations of PKS 2155 − ∼ . × − Hz and ∼ . × − Hz. Lachowicz et al. (2009) also examined thePSDs in the energy bands 0.3–2 keV and 2–10 keV,and found the maxima were significant in the softband but not in the hard band, although they werestill found in the hard band. Gaur et al. (2010) andGonz´alez-Mart´ın & Vaughan (2012) find this QPO inonly one of many
XMM-Newton observations of thissource. If the maxima are significant and found in bothenergy bands, this argues for a light crossing time inter-pretation, with timescale t lc = 3 . × s and R ′ ≈ ct lc δ D z ) = 3 . × (cid:18) δ D (cid:19) cm . (78)This is larger than the size from the HESS observa-tion (Section 6.1). The XMM-Newton observations weretaken on 2006 May 1, and the HESS observations on 2006July 28, which could account for the discrepancy. The2 Finke & Beckeremitting region size could have been different at differ-ent times. A number of claimed significant detections ofQPOs from the literature are listed in Table 2. DISCUSSION
We have presented a new theoretical formalism formodeling the variability of blazars, based on an analyt-ical solution to the underlying electron continuity equa-tion governing the distribution of radiating electrons ina homogeneous blob moving out in the jet. The anal-ysis was carried out in the Fourier domain, and theresults are therefore directly comparable with obser-vational Fourier data products such as the PSDs andtime/phase lags. This formalism assumes emission froma jet closely aligned with the line of sight, so that theemission produced in the comoving frame is Dopplershifted into the observer’s frame. Internal shocks in thejet randomly accelerate electrons to high energies, whichare then injected into an emitting region at random in-tervals with that variability characterized by a power-law in Fourier space. The observable radiation producedby the electrons is affected by cooling, electron escape,and the light-travel time across the blob. The modelmakes specific predictions regarding the the PSDs andFourier-dependent time lag components resulting fromsynchrotron, EC, and SSC emission, and it successfullyreproduces the characteristics of the colored noise seenin nearly all blazars.The study presented here is a first attempt at exam-ining blazar variability with this formalism. As such, itmakes a number of simplifying assumptions.1. It assumes the only thing which varies with time ina blazar is the rate at which electrons are injectedinto the emitting region. All other parameters—the magnetic field strength ( B ), the size of theemitting region ( R ′ ), the electron injection power-law index ( q ), the jet’s angle to the line of sight( θ ), and so on—are assumed not to vary. Althoughthis is likely an over-simplification, we note thatthe PSDs of PKS 0537 −
441 can be explained byonly varying the electron distribution of the source(D’Ammando et al. 2013), so that in some casesthis may be justified.2. In computing emission, we have used simplified δ -function expressions for the synchrotron andCompton emission, and assumed Thomson scatter-ing. Using more accurate expressions, in particular,the full Compton cross-section for the energy lossescould lead to interesting effects (Dermer & Atoyan2002; Moderski et al. 2005; Sikora et al. 2009;Dotson et al. 2012).3. The calculations neglect SSC cooling, which isquite difficult to model analytically (Schlickeiser2009; Zacharias & Schlickeiser 2010, 2012a,b).This would likely not be important for FSRQs,where the EC component probably dominates thecooling, but could be important for BL Lac objectsthat do not have a strong external radiation field.4. We have neglected the details of the accelerationmechanism. Such a mechanism may produce inter-esting features in PSDs that could be observable. 5. Although we take light travel-time effects into ac-count, we assume all parts of the blob vary simul-taneously, which is obviously not the case.In the future, we will perform more detailed analyseswhich explore these more complicated cases.What are the best wavebands to observe blazars andcompare with the theory outlined in this work? Ob-servations at (electromagnetic) frequencies lower than ∼ Hz would not be useful, since here the emis-sion is likely dominated by the superposition of manyself-absorbed jet components (Konigl 1981). For low-synchrotron peaked (LSP) blazars (including almost allFSRQs; Finke 2013), The optical and GeV γ -rays wouldhave features at high Fourier frequencies due the rapidenergy losses of electrons that produce this radiation,and thus this emission could be quite interesting to ob-serve. Observing such short timescales with the Fermi -LAT could be difficult, since for this instrument one mustusually integrate over fairly long timescales ( & a fewhours for all but the brightest sources) to get a signif-icant detection. Bright flares and adaptive light curvebinning may be helpful in producing accurate LAT PSDsat high Fourier frequencies (Lott et al. 2012). Significantdetections can be made in the optical with shorter in-tegration times ( ∼ a few minutes or even fractions of aminute), leading to better PSDs at high frequencies (e.g.,Rani et al. 2010) although in FSRQs this could be con-taminated by thermal emission from an accretion disk.Observing LSPs at wavebands (e.g., infrared) whereemission is dominated by less energetic electrons couldprobe the escape timescale. PSDs produced from simul-taneous light curves at multiple frequencies would be ex-tremely helpful for verifying the predictions of this paper.At low Fourier frequencies the PSDs should have essen-tially identical power-law shapes at all wavebands, andat higher Fourier frequencies features associated with thelight-crossing timescale could be identified, and shouldbe the same in all wavebands. This could serve as astrong test for the theory presented in this paper. Wealso predict that all breaks in observed synchrotron andEC PSD power laws should be by 2, i.e., from ∝ f − a to ∝ f − ( a +2) , and more gradual breaks could be observedin SSC PSDs. Kataoka et al. (2001) observed smallerbreaks in the X-ray PSDs of Mrk 421, Mrk 501, andPKS 2155 −
304 based on
ASCA and
Rossi X-ray TimingExplorer observations. However, note that they did notobtain acceptable fits to their PSDs with broken power-laws with all parameters left free to vary. Also, it maybe that including SSC cooling could modify the the PSDso that a breaks other than 2 are possible, although thatis beyond the scope of this paper.For high synchrotron-peaked blazars (HSPs), almostall of which are BL Lac objects, the γ -ray emission isexpected to be from SSC. We do not predict any fea-tures in SSC PSDs from cooling or escape, althoughfeatures from the light crossing time are still expected.These predictions could be tested with Fermi -LAT andVHE γ -ray instruments such as MAGIC, HESS, VER-ITAS, and the upcoming CTA. For synchrotron emis-sion, however, features from cooling and escape shouldbe present. An X-ray telescope could be used to probeemission from the highest energy electrons potentiallyseeing cooling features. The proposed Large Observa- ourier Analysis of Blazar Variability 13
TABLE 2Claimed QPOs from blazars reported in the literature.
Authors Object Bandpass QPO frequency [Hz]Espaillat et al. (2008) a
3C 273 0.75–10 keV 3 . × − Lachowicz et al. (2009) PKS 2155 −
304 0.2–10 keV 6 × − Gupta et al. (2009) S5 0716+714 V and R band variousRani et al. (2010) S5 0716+714 R band 1 . × − This claimed QPO is disputed by Mohan et al. (2011) andGonz´alez-Mart´ın & Vaughan (2012). tory for X-ray Timing (LOFT) spacecraft could be veryuseful for this (Donnarumma et al. 2013).
LOFT willprovide excellent timing coverage ( ∼ ms) with a rela-tively large effective area. PSDs produced from opticallight curves of HSPs could be used to probe the escapetimescale. As with LSP blazars, PSDs produced by si-multaneous light curves from multiple wavebands would be extremely helpful for the study of HSP blazars.We are grateful to the anonymous referee for insight-ful suggestions that helped improve the discussion andpresentation, and to C. Dermer for useful discussions.J.D.F. was supported by the Office of Naval Research. APPENDIX A. SOLUTION TO THE FOURIER TRANSFORM OF THE CONTINUITY EQUATION
Here we solve the differential equation − πif ˜ N e ( γ, f ) + ∂∂γ [ ˙ γ ( γ ) ˜ N e ( γ, f )] + ˜ N e ( γ, f ) t esc ( γ ) = ˜ Q ( γ ; f ) (A1)for ˙ γ ≤
0. This equation can be rearranged, ∂∂γ [ ˙ γ ( γ ) ˜ N e ( γ, f )] + (cid:20) t esc − iω (cid:21) ˜ N e ( γ, f ) = ˜ Q ( γ ; f ) (A2)recalling that ω = 2 πf . One can multiply both sides byexp (cid:20)Z ∞ γ dγ ′ | ˙ γ ( γ ′ ) | (cid:18) t esc − iω (cid:19)(cid:21) (A3)and then rearrange the left side so that ddγ (cid:26) ˙ γ ( γ ) ˜ N e ( γ, f ) exp (cid:20)Z ∞ γ dγ ′ | ˙ γ ( γ ′ ) | (cid:18) t esc − iω (cid:19)(cid:21)(cid:27) = ˜ Q ( γ ; f ) exp (cid:20)Z ∞ γ dγ ′ | ˙ γ ( γ ′ ) | (cid:18) t esc − iω (cid:19)(cid:21) . (A4)Integrating both sides gives | ˙ γ ( γ ) | ˜ N e ( γ, f ) exp (cid:20)Z ∞ γ dγ ′ | ˙ γ ( γ ′ ) | (cid:18) t esc − iω (cid:19)(cid:21) = Z ∞ γ dγ ′ ˜ Q ( γ ′ ; f ) exp (cid:20)Z ∞ γ ′ dγ ′′ | ˙ γ ( γ ′′ ) | (cid:18) t esc − iω (cid:19)(cid:21) . (A5)Solving this for ˜ N ( γ, f ) results in˜ N e ( γ, f ) = 1 | ˙ γ ( γ ) | Z ∞ γ dγ ′ ˜ Q ( γ ′ ; f ) exp " − Z γ ′ γ dγ ′′ | ˙ γ ( γ ′′ ) | (cid:18) t esc ( γ ′′ ) − iω (cid:19) . (A6) B. NORMALIZATION OF ELECTRON INJECTION FUNCTION
The time average of the total power injected in electrons is h L inj i = m e c ∆ t Z ∞ dγ γ Z ∆ t dt Q ( γ, t ) (B1)where ∆ t is the length of the time interval over which the electrons are injected. Using Equation (2), h L inj i = m e c ∆ t Z ∞ dγ γ Z ∞−∞ df ˜ Q ( γ, f ) Z ∆ t dt exp[ − πif t ] . (B2)4 Finke & BeckerSubstituting Equation (14), ˜ Q ( γ ; f ) = Q ( f /f ) − a/ γ − q H ( γ ; γ , γ ) H ( f ; f , f ) , (B3)and integrating over γ , one gets h L inj i = m e c ∆ t Q G ( q, γ , γ ) Z f f df ( f /f ) − a/ Z ∆ t dt exp[ − πif t ] (B4)where G ( q, γ , γ ) = (cid:26) ( γ − q − γ − q ) / ( q −
2) for q = 2ln( γ /γ ) for q = 2 . (B5)Performing the integral over time gives h L inj i = m e c − πi ∆ t Q G ( q, γ , γ ) Z f f dff ( f /f ) − a/ [exp( − πif ∆ t ) − . (B6)This can be rewritten as h L inj i = m e c − πi ∆ t Q G ( q, γ , γ ) [ I r ( a, f , f ) − iI i ( a, f , f ) − I ( a, f , f )] (B7)where I r ( a, f , f ) = Z f f dff ( f /f ) − a/ cos(2 πf ∆ t ) , (B8) I i ( a, f , f ) = Z f f dff ( f /f ) − a/ sin(2 πf ∆ t ) , (B9)and I ( a, f , f ) = Z f f dff ( f /f ) − a/ = (cid:26) /a [( f /f ) − a/ − ( f /f ) − a/ ] a = 0ln( f /f ) a = 0 . (B10)Multiplying Equation (B7) by its complex conjugate gives |h L inj i| = (cid:20) m e c Q G π ∆ t (cid:21) [ I r + I i − I r I + I ] , (B11)or, solving for Q , Q = 2 π ∆ t h L inj i m e c G p I r + I i − I r I + I . (B12) C. LIGHT TRAVEL TIME
Let us take a cylindrical blob, as shown in cross section in Figure 8. The blob has length 2 R , and everywhere withinthe blob radiation is emitted simultaneously as a function of time t as g ( t ). The entire blob has length R and is dividedinto N individual segments, each with length ∆ x = 2 R/N . The radiation emitted by each individual segment is acorresponding fraction of the whole, g ( t ) /N = g ( t )∆ x/ R . The radiation an observer co-moving with the blob sees atany given time t obs will be a sum over the individual segments at that time, h ( t obs ) = g ( t obs − x /c ) ∆ x/ (2 R ) + g ( t obs − x /c ) ∆ x/ (2 R ) + (C1) . . . + g ( t obs − R/c ) ∆ x/ (2 R )or h ( t obs ) = 12 R N X j =1 g ( t obs − x j /c ) ∆ x . (C2)As N → ∞ , h ( t obs ) → R Z R dx g ( t obs − x/c ) . (C3)ourier Analysis of Blazar Variability 15 xto observer ∆ x 2R Fig. 8.—
Sketch of geometry of emitting blob for the purpose of computing light travel time effects. This sketch is in the frame co-movingwith the blob. The blob of length 2 R is divided into N pieces each with length ∆ x . Since t = x/c , h ( t obs ) = c R Z R/c dt g ( t obs − t ) . (C4)Let us now move to a frame where the blob is moving relative to the observer so that she or he sees a time t obs = t ′ obs (1 + z ) /δ D , where now all lengths and times in the co-moving frame will be primed. Then h ( t obs ) = c R ′ Z R ′ /c dt ′ g ( t ′ obs − t ′ ) (C5)= c R ′ Z R ′ /c dt ′ g (cid:18) t obs δ D z − t ′ (cid:19) . This is similar to the “time slices” of Chiaberge & Ghisellini (1999). We use a cylindrical normal blob geometryfor simplicity here, although a spherical one would be more consistent with the SSC calculation. Since the actualgeometry of an emitting blob is not known, differences should not be too great. A similar derivation for a sphericalgeometry is given by Zacharias & Schlickeiser (2013). D. FOURIER TRANSFORM INCLUDING LIGHT TRAVEL TIME
In this appendix we derive Equation (47) from Equation (40). From the definition of the inverse Fourier transform,Equation (2), N e ( γ ; t ) = Z ∞−∞ df ˜ N e ( γ, f ) exp( − πif t ) . (D1)Putting this in Equation (40) and rearranging gives F syǫ ( t ) = K (1 + z ) t lc δ D Z ∞−∞ df ′ ˜ N e ( γ ′ , f ′ ) exp (cid:20) − πif ′ δ D t z (cid:21) Z R ′ /c dt ′ exp(2 πif ′ t ′ ) , (D2)recalling that t lc = 2 R ′ (1 + z ) / ( cδ D ). Performing the integral over t ′ gives F syǫ ( t ) = K (1 + z )2 πit lc δ D Z ∞−∞ df ′ f ′ ˜ N e ( γ ′ , f ′ ) exp (cid:20) − πif ′ δ D t z (cid:21) (cid:26) exp (cid:20) πif ′ R ′ c (cid:21) − (cid:27) . (D3)6 Finke & BeckerThe Fourier transform of the synchrotron flux light curve is defined as (see Equation [1])˜ F syǫ ( f ) = Z ∞−∞ dt F ǫ ( t ) exp(2 πif t ) . (D4)Substituting Equation (D3) for F ǫ ( t ) in this equation gives˜ F syǫ ( f ) = K (1 + z )2 πit lc δ D Z ∞−∞ df ′ f ′ ˜ N e ( γ ′ , f ′ ) (cid:26) exp (cid:20) πif ′ R ′ c (cid:21) − (cid:27) (D5) × Z ∞−∞ dt exp (cid:20) πit (cid:18) f − f ′ δ D z (cid:19)(cid:21) . The integral over t has the form of a Dirac δ -function (Equation [4]), so˜ F syǫ ( f ) = K (1 + z )2 πit lc δ D Z ∞−∞ df ′ f ′ ˜ N e ( γ ′ , f ′ ) (cid:26) exp (cid:20) πif ′ R ′ c (cid:21) − (cid:27) δ (cid:18) f − f ′ δ D z (cid:19) . (D6)Using the well-known property for δ functions, δ (cid:18) f − f ′ δ D z (cid:19) = 1 + zδ D δ (cid:18) f ′ − (1 + z ) fδ D (cid:19) , (D7)one can perform the integral over f ′ in Equation (D6) to get˜ F syǫ ( f ) = K (1 + z )2 πif t lc δ D ˜ N e (cid:18) γ ′ , (1 + z ) fδ D (cid:19) (cid:26) exp (cid:20) πif (1 + z ) R ′ cδ D (cid:21) − (cid:27) . (D8)This is Equation (47). REFERENCESAbdo, A. A., et al. 2010a, ApJS, 188, 405—. 2010b, ApJ, 722, 520—. 2010c, ApJ, 716, 30—. 2011a, ApJ, 733, L26—. 2011b, ApJ, 736, 131—. 2011c, ApJ, 727, 129Ackermann, M., et al. 2010, ApJ, 721, 1383Aharonian, F., et al. 2007, ApJ, 664, L71Begelman, M. C., Fabian, A. C., & Rees, M. J. 2008, MNRAS,384, L19B la˙zejowski, M., Sikora, M., Moderski, R., & Madejski, G. M.2000, ApJ, 545, 107Bloom, S. D., & Marscher, A. P. 1996, ApJ, 461, 657B¨ottcher, M., & Chiang, J. 2002, ApJ, 581, 127B¨ottcher, M., & Reimer, A. 2004, ApJ, 609, 576Chatterjee, R., et al. 2012, ApJ, 749, 191Chen, X., Fossati, G., B¨ottcher, M., & Liang, E. 2012, MNRAS,424, 789Chen, X., Fossati, G., Liang, E. P., & B¨ottcher, M. 2011,MNRAS, 416, 2368Chiaberge, M., & Ghisellini, G. 1999, MNRAS, 306, 551Cui, W. 2004, ApJ, 605, 662D’Ammando, F., et al. 2013, MNRAS, 431, 2481Dermer, C. D., & Atoyan, A. M. 2002, ApJ, 568, L81Dermer, C. D., & Menon, G. 2009, High Energy Radiation fromBlack Holes: Gamma Rays, Cosmic Rays, and NeutrinosDermer, C. D., & Schlickeiser, R. 1993, ApJ, 416, 458—. 2002, ApJ, 575, 667Donnarumma, I., Tramacere, A., Turriziani, S., Costamante, L.,Campana, R., De Rosa, A., & Bozzo, E. 2013, ArXiv:1310.6965Dotson, A., Georganopoulos, M., Kazanas, D., & Perlman, E. S.2012, ApJ, 758, L15Edelson, R., Mushotzky, R., Vaughan, S., Scargle, J., Gandhi, P.,Malkan, M., & Baumgartner, W. 2013, ApJ, 766, 16Espaillat, C., Bregman, J., Hughes, P., & Lloyd-Davies, E. 2008,ApJ, 679, 182Finke, J. D. 2013, ApJ, 763, 134Finke, J. D., Dermer, C. D., & B¨ottcher, M. 2008, ApJ, 686, 181Gaur, H., Gupta, A. C., Lachowicz, P., & Wiita, P. J. 2010, ApJ,718, 279Ghisellini, G., Tavecchio, F., Foschini, L., & Ghirlanda, G. 2011,MNRAS, 414, 2674 Gonz´alez-Mart´ın, O., & Vaughan, S. 2012, A&A, 544, A80Gupta, A. C., Srivastava, A. K., & Wiita, P. J. 2009, ApJ, 690,216Hayashida, M., et al. 2012, ApJ, 754, 114Joshi, M., & B¨ottcher, M. 2007, ApJ, 662, 884Kataoka, J., et al. 1999, ApJ, 514, 138—. 2001, ApJ, 560, 659Konigl, A. 1981, ApJ, 243, 700Kroon, J. J., & Becker, P. A. 2014, ApJ, 785, L34Lachowicz, P., Gupta, A. C., Gaur, H., & Wiita, P. J. 2009,A&A, 506, L17Li, H., & Kusunose, M. 2000, ApJ, 536, 729Lott, B., Escande, L., Larsson, S., & Ballet, J. 2012, A&A, 544,A6Malmrose, M. P., Marscher, A. P., Jorstad, S. G., Nikutta, R., &Elitzur, M. 2011, ApJ, 732, 116Mastichiadis, A., Petropoulou, M., & Dimitrakoudis, S. 2013,MNRAS, 434, 2684McHardy, I. 2008, in Blazar Variability across theElectromagnetic SpectrumModerski, R., Sikora, M., Coppi, P. S., & Aharonian, F. 2005,MNRAS, 363, 954Mohan, P., Mangalam, A., Chand, H., & Gupta, A. C. 2011,Journal of Astrophysics and Astronomy, 32, 117Nakagawa, K., & Mori, M. 2013, ApJ, 773, 177Nolan, P. L., et al. 2012, ApJS, 199, 31Rani, B., Gupta, A. C., Joshi, U. C., Ganesh, S., & Wiita, P. J.2010, ApJ, 719, L153Sbarrato, T., Ghisellini, G., Maraschi, L., & Colpi, M. 2012,MNRAS, 421, 1764Schlickeiser, R. 2009, MNRAS, 398, 1483Schlickeiser, R., B¨ottcher, M., & Menzler, U. 2010, A&A, 519,A9+Sikora, M., Begelman, M. C., & Rees, M. J. 1994, ApJ, 421, 153Sikora, M., Stawarz, L., Moderski, R., Nalewajko, K., &Madejski, G. M. 2009, ApJ, 704, 38Sobolewska, M. A., Siemiginowska, A., Kelly, B. C., & Nalewajko,K. 2014, ApJ, in press, arXiv:1403.5276Tanaka, Y. T., et al. 2013, ApJ, 777, L18Vaughan, S. 2005, A&A, 431, 391Vaughan, S., & Uttley, P. 2006, Advances in Space Research, 38,1405 ourier Analysis of Blazar Variability 17ourier Analysis of Blazar Variability 17