Fourier Analysis of Blazar Variability

# 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
###### Abstract

Blazars display strong variability on multiple timescales and in multiple radiation bands. Their variability is often characterized by power spectral densities (PSDs) and time lags plotted as functions of the Fourier frequency. We develop a new theoretical model based on the analysis of the electron transport (continuity) equation, carried out in the Fourier domain. The continuity equation includes electron cooling and escape, and a derivation of the emission properties includes light travel time effects associated with a radiating blob in a relativistic jet. The model successfully reproduces the general shapes of the observed PSDs and predicts specific PSD and time lag behaviors associated with 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 to flat-spectrum radio quasars (FSRQs), where there are hints that some of the predicted features have already been observed. We also find that FSRQs should have steeper PSD power-law indices than BL Lac objects at Fourier frequencies  Hz, in qualitative agreement with previously reported observations by the Fermi Large Area Telescope.

BL Lacertae objects: general — quasars: general — radiation mechanisms: nonthermal — galaxies: active — galaxies: jets
slugcomment: ApJ, accepted; July 12, 2019

## 1. Introduction

Blazars, active galactic nuclei (AGN) with jets moving at relativistic speeds aligned with our line of sight, are the most plentiful of the identified point sources observed by the Fermi Large Area Telescope (LAT; Abdo et al., 2010a; Nolan et al., 2012). Their spectral energy distributions (SEDs) are dominated by two components. At lower frequencies there is a component produced via synchrotron emission, peaking in the radio through X-rays. There is also a high energy component in -rays, most likely due to Compton-scattered emission. The seed photons for Compton scattering could be the synchrotron photons themselves (synchrotron self-Compton or SSC; Bloom & Marscher, 1996), or they could be produced outside the jet (external Compton or EC) by the accretion disk (Dermer & Schlickeiser, 1993, 2002), the broad line region (Sikora et al., 1994), or a dust torus (Kataoka et al., 1999; Błażejowski et al., 2000). The Doppler-boosting due to the combination of relativistic speed and a small jet inclination angle amplifies the observed 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 curves of blazars on long time scales. For some sources, -ray data have been supplemented by high cadence observations in the radio through very-high energy (VHE) -rays creating unprecedented light curves with few gaps in wavelength or time (e.g., Abdo et al., 2011b, c). Blazar variability is often characterized by power spectral densities (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 representations of the Fourier transform without phase information. Although the LAT can provide long baseline, high cadence light curves, it has difficulty probing short timescales for all but the brightest flares. However, shorter-timescale variability may be observed with optical, 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 energy channels are computed from light curves (e.g., Zhang, 2002). The PSDs of light curves at essentially all wavelengths resemble power-laws in frequency, , with typically , usually steeper than PSDs found from Seyfert galaxies (Kataoka et al., 2001). Despite the popularity of PSDs for characterizing variability, their theoretical motivation has not been thoroughly explored (although see Mastichiadis et al., 2013).

In this paper, our goal is to bridge the gap between theory and observations by exploiting a powerful new mathematical approach for the modeling and interpretation observed PSDs and time lags. Theoretically, the variability of blazars is often described by a continuity equation (e.g., Chiaberge & Ghisellini, 1999; Li & Kusunose, 2000; Böttcher & Chiang, 2002; Chen et al., 2011, 2012). This equation describes the evolution of electrons in a compact region of the jet, which is homogeneous by assumption. The electrons are injected as a function of time and energy and the electron distribution evolves due to energy loss and escape. The expected electromagnetic emission observers might detect can be compared with observations (e.g., Böttcher & Reimer, 2004; Joshi & Böttcher, 2007). This can allow one to explore individual flares; however, in this paper, we take a different approach, by studying the electron continuity equation in the Fourier domain. This allows the exploration of individual flares, as well as the study of aggregated long-timescale variability of sources using Fourier transform-related quantities such as PSDs and phase and time lags. We focus here on long-timescale variability, including multiple epochs of flaring and quiescence. A similar study applying the same Fourier transform concept to the modeling of time lags in accreting Galactic black hole candidates has been recently carried out by Kroon & Becker (2014).

We begin in Section 2, by defining the Fourier transform and its inverse, and various other functions used throughout the rest of the paper. In Section 3, we explore analytic solutions to the continuity equation in the Fourier domain, and present solutions in terms of PSDs. In Section 4, we explore the solutions in terms of time lags between different electron Lorentz factor “channels” as a function of Fourier frequency. We explore the expected synchrotron, SSC, and EC PSDs and Fourier frequency-dependent time lags in Section 5, including light travel time effects due to the emitting region’s finite size. We discuss applications to some PSDs and time lags in the literature in Section 6, and conclude with a discussion of our simplifying assumptions and observational prospects in Section 7. Several of the detailed derivations are relegated to the appendices.

## 2. Definitions

Several definitions of the Fourier transform and associated quantities are used throughout the literature. Here we make the definitions used in this paper explicit. For a real function , we define the Fourier transform by

 ~x(f)=∫∞−∞dt x(t) e2πift =∫∞−∞dt x(t) eiωt , (1)

where . We will indicate the Fourier transform by a tilde, the Fourier frequency by , and the angular frequency by . We define the inverse Fourier transform by

 x(t)=∫∞−∞df ~x(f) e−2πift=12π∫∞−∞dω ~x(ω) e−iωt . (2)

We define the PSD

 S(f)=|~x(f)|2=~x(f)~x∗(f) (3)

where the asterisk indicates the complex conjugate. We will make use of the related representation of the Dirac -function

 δ(t−t0)=∫∞−∞df e2πif(t−t0) . (4)

We will use two versions of the Heaviside function, a “step” function defined by

 H(x)={1x>00otherwise , (5)

and the two-sided Heaviside “top-hat” function,

 H(x;a,b)={1a

In Section 3.5 we will use the lower incomplete Gamma function given by

 γg(a,x)=∫x0dy ya−1e−y . (7)

## 3. Continuity Equation in Fourier Space

### 3.1. General Solution

The evolution of electrons in a nonthermal plasma “blob” can be described by a continuity equation given by (e.g., Chiaberge & Ghisellini, 1999; Li & Kusunose, 2000; Böttcher & Chiang, 2002; Chen et al., 2011, 2012)

 ∂Ne∂t+∂∂γ[˙γ(γ,t)Ne(γ;t)]+Ne(γ;t)tesc(γ,t)=Q(γ,t) , (8)

where gives the number of electrons with Lorentz factor between and at time . Here is the rate at which electrons lose or gain energy, is the escape timescale, and is the rate at which electrons are injected in the jet. In the simple model presented here, we will assume that size of the plasma blob does not change with time, so that adiabatic losses can be neglected. We will also assume that the electron distribution in the blob is homogeneous and isotropic, and that variations occur throughout the blob simultaneously. This is a common and useful assumption, although it is somewhat unphysical, as in order for the blob to be causally connected, variations cannot propagate through the blob faster than the speed of light .

Assuming that and are independent of , we can take the Fourier transform of both sides of the continuity equation leading to

 −2πif~Ne(γ,f)+∂∂γ[˙γ(γ)~Ne(γ,f)]+~Ne(γ,f)tesc(γ) =~Q(γ,f) (9)

where is the Fourier transformed source term. This is a linear ordinary first-order differential equation with a relatively simple solution. It is shown in Appendix A that if

 ~Ne(γ,f)=1|˙γ(γ)|∫∞γdγ′ ~Q(γ′,f) ×exp[−∫γ′γdγ′′|˙γ(γ′′)|(1tesc(γ′′)−iω)] . (10)

If is independent of and cooling is from synchro-Compton processes, so that , then one can perform the integral in the exponent, and

 γ2~Ne(γ,f)=1ν exp[−1νγ(1tesc−iω)] ×∫∞γdγ′ ~Q(γ′,f) exp[1νγ′(1tesc−iω)] . (11)

### 3.2. Green’s Function Solution

Consider an instantaneous injection of monoenergetic electrons with Lorentz factor at . Then in Equation (3.1) and one gets

 γ2~Ne(γ,f) =Q0ν exp[−1ν(1γ−1γ0)(1tesc−iω)] ×H(γ0−γ) . (12)

The PSD is

 S(γ,f) =|γ2~Ne(γ,f)|2 =[Q0ν]2exp[−2νtesc(1γ−1γ0)]H(γ0−γ) . (13)

The result is white noise for all electron Lorentz factors () and Fourier frequencies ().

### 3.3. Colored Noise

Since the PSDs of blazars resemble colored noise, and electrons are generally thought to be injected as power laws in , one might expect that

 ~Q(γ,f) =Q0(f/f0)−a/2γ−q ×H(γ;γ1,γ2)H(f;f1,f2) (14)

where is some constant frequency and . That is, in the jet, shocks will occur randomly which accelerate and inject particles as a power-law distribution in between and with index . We will deal only with frequencies in the range . These limits are needed for the PSD to be normalized to a finite value. Frequencies greater than the inverse of the blob’s light crossing timescale are particularly unphysical, although we allow this for two reasons. First, it allows us to compare with other theoretical studies that allow variations faster than the light crossing timescale (e.g., Chiaberge & Ghisellini, 1999; Zacharias & Schlickeiser, 2013). Second, our blob is already unphysical, since we allow variations throughout the blob simultaneously in the blob’s comoving frame. The normalization constant is related to the time-averaged power injected in electrons over a time interval by

 Q0=2πΔt⟨Linj⟩mec2G√I2r+I2i−2IrI0+I20 . (15)

A derivation of this equation and definitions of the quantities , , , and can be found in Appendix B. With , given by Equation (3.3), Equation (3.1) can be rewritten as

 γ2~Ne(γ,f) =Q0(f/f0)−a/2exp[−1νγ(1tesc−iω)]νq−2 ×(1tesc−iω)1−q∫umaxumindu uq−2 eu (16)

where

 umin=1νγ2(1tesc−iω) (17)

and

 umax=1νmax(γ,γ1)(1tesc−iω) . (18)

### 3.4. Electron Injection Index q=2

It is instructive to look at the case where . In this case, the remaining integral in Equation (3.3) can easily be performed analytically. Then

 γ2~Ne(γ,f) =Q0(f/f0)−a/21/tesc−iωexp[−1νγ(1tesc−iω)] ×[eumax−eumin] , (19)

and the PSD is

 S(γ,f) =|γ2~Ne(γ,f)|2 =exp[−2νγtesc]Q20(f/f0)−a(1t2esc+ω2)−1 ×{exp[2νγ2tesc]+exp[2νmax(γ,γ1)tesc] −2exp[1νtesc(1max(γ,γ1)+1γ2)] ×cos[ων(1max(γ,γ1)−1γ2)]} . (20)

We identify asymptotes for the PSD for , Equation (3.4). For these asymptotes we assume .

1. If , and then

 S(γ,f)≈Q20(f/f0)−aν2γ2 . (21)
2. If then

 S(γ,f)≈Q20(f/f0)−a−2f20π2sin2(πfνγ) . (22)
3. If then

 S(γ,f) ≈Q20(f/f0)−a1/t2esc+(2πf)2exp[−2tesc(1νγ−tcool)] (23) ×{1−exp[−tcooltesc]cos[2πftcool]}

where we define .

The electron PSD resulting from Equation (3.4) is plotted in Figure 1 for parameters described in the caption, which are fairly standard ones for FSRQs. We use here, which represents an instantaneous injection of power-law particles at , to more easily display the observable features, a number of which are present. For the , 170, and curves, where ), 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πtesc)−1 .

This is in agreement with asymptote 3 above. In general, by examining asymptote 3, it is clear that for low a break in the PSD will be found at a frequency of . Since the PSD measures periodic variability, this indicates that for low , periodic variability on timescales less than the escape timescale is less preferred. This is because electrons will always be escaping at a single timescale, which does not vary with time. One can also see in the curve in Figure 1 at high structure related to the cosine seen in asymptote 3, with local minima at integer multiples of .

In the and curves, at high () 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−1cool ,

in agreement with asymptotes 1 and 2. Thus, variability on timescales less than the cooling timescales will be less periodic, since cooling on those smaller timescales will always be present. It is also clear that at high minima from the term in asymptote 2 occur at integer multiples of . Since at high values of , where the cooling is strongest, cooling on timescales will be immediate, and so periodic variability on timescales that are integer multiples of will also be strongly avoided.

### 3.5. Electron Injection Index q≠2

In this case, there is no simple analytic solution to equation (3.3), although it can be written with the incomplete Gamma function as

 γ2~Ne(γ,f) =Q0(f/f0)−a/2exp[−1νγ(1tesc−iω)] ×νq−2(iω−1tesc)1−q ×[γg(q−1,−umax)−γg(q−1,−umin)] . (24)

We compute the function numerically, and results can be seen in Figure 2 for . This confirms the features seen in the case are also seen for other values of , although the minima at high are not as pronounced.

## 4. Time Lags

The phase lag between two electron Lorentz factor “channels”, and as a function of Fourier frequency () can be calculated from

 Δϕ(γa,γb,f)=arctan{YI(γa,γb,f)YR(γa,γb,f)]} (25)

where

 [γ2a~Ne(γa,f)][γ2b~Ne(γb,f)]∗= YR(γa,γb,f)+iYI(γa,γb,f) . (26)

This implies that

 YR(γa,γb,f)=Re[γ2a~Ne(γa,f)]Re[γ2b~Ne(γb,f)] +Im[γ2a~Ne(γa,f)]Im[γ2b~Ne(γb,f)] (27)

and

 YI(γa,γb,f)=Re[γ2b~Ne(γb,f)]Im[γ2a~Ne(γa,f)] −Re[γ2a~Ne(γa,f)]Im[γ2b~Ne(γ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 (3.4) the time lag is

 ΔT(γa,γb,f)=12πfarctan{ZI(γa,γb,f)ZR(γa,γb,f)} (30)

where

 ZI(γa,γb,f) =exp[−1νtesc(1γa+1γb)] ×{exp[2νtescγ2]sin[ων(1γa−1γb)] +exp[1νtesc(1max(γ1,γa)+1max(γ1,γb))] ×sin[ων(1max(γ1,γb)−1γb −exp[1νtesc(1max(γ1,γa)+1γ2)] ×sin[ων(1γa−1max(γ1,γa)+1γ2−1γb)] −exp[1νtesc(1γ2+1max(γ1,γb))] ×sin[ων(1γa−1γ2+1max(γ1,γb)−1γb)] } (31)

and is the same as except with in place of . Several examples of electron time lags can be seen in Fig. 3.

For and the time delay will be approximately independent of frequency, with value

 ΔT(γa,γb,f)≈1νAI(γa,γb)AR(γa,γb) (32)

where

 AI(γa,γb,f) =(1γa−1γb)exp[−1νtesc(1γa+1γb)] +1γbexp[−1νtescγb] −1γaexp[−1νtescγa] (33)

and

 AR(γa,γb,f) =1+exp[−1νtesc(1γa+1γb)] −exp[−1νtescγb]−exp[−1νtescγa] , (34)

where we also assumed and . At these low values of , the lags are positive for 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 and then

 ΔT(γa,γb,f)≈12ν(1γa−1γb) . (35)

An example of this can be seen in Fig. 3, with the , curve. If then

 ΔT(γa,γb,f)≈tesc{1−1νtescγb −1νtescγaexp[−1νtescγa]} . (36)

In Fig. 3, an example can be seen with the , curve.

If and then . In Fig. 3, the curve that most closely approximates this is the , curve.

At high , The behavior is quite complex. By inspecting Equation (4), we can see that the important frequencies are those that make the sine or cosine terms go to 0. They will be the integer or half-integer multiples of

 f=t−1cool,a=νγa (37) f=t−1cool,b=νγb (38) f=(tcool,a−tcool,b)−1 . (39)

There are many local minima and maxima at the integer or half-integer multiples of these values.

## 5. Emission and Light Travel Time Effects

In the previous sections we have explored the PSD and time lags for the electron distribution. However, what is observed is the emission of these electrons, through synchrotron or Compton-scattering. We will assume that the emitting region is spherical with co-moving radius and homogeneous, containing a tangled magnetic field of strength . The blob is moving with a relativistic speed ( being the speed of light), giving it a bulk Lorentz factor . The blob is moving with an angle to the line of sight giving it a Doppler factor . Although we make the simplifying assumption the blob is homogeneous, with variations in the electron distribution taking place throughout the blob simultaneously, since it has a finite size, photons will reach the observer earlier from the closer part than the farther part, and thus one must integrate over the time in the past, . This is similar to the “time slices” of Chiaberge & Ghisellini (1999). We note again that Fourier frequencies higher than the inverse of the light crossing timescale are unphysical in the simple model we present here.

### 5.1. Synchrotron and External Compton

Taking the light travel time into account, in the -approximation the observed flux at observed energy (in units of the electron rest energy) from synchrotron or Compton scattering of an external isotropic monochromatic radiation field as a function of the observer’s time is

 Fϵ(t)=K(1+z)tlcδD∫2R′/c0dt′ Ne(γ′;tδD1+z−t′) (40)

where

 tlc=2R′(1+z)cδD (41)

is the light crossing time in the observer’s frame. A derivation of the light travel time effect can be found in Appendix C. For synchrotron emission,

 K=Ksy=δ4D6πd2LcσTuBγ′3sy , (42)

and

 γ′=γ′sy=√ϵ(1+z)δDϵB . (43)

The Thomson cross-section is , the Poynting flux energy density is , where , the redshift of the source is , and the luminosity distance to the source is . For external Compton (EC) scattering,

 K=KEC=δ6D6πd2LcσTu0γ′3EC , (44)
 γ′=γ′EC=1δD√ϵ(1+z)2ϵ0 (45)

(Dermer & Menon, 2009). Here the external radiation field energy density and photon energy (in units of the electron rest energy) are and , respectively. This approximation is valid in the Thomson regime, i.e., when . Primed quantities refer to the frame co-moving with the emitting region. The cooling rate parameter is

 ν=43mec2cσT(uB+Γ2u0) (46)

where we ignore the effects of SSC cooling.

It is shown in Appendix D that the Fourier transform of Equation (40) is

 ~Fϵ(f)=K(1+z)2πiftlcδD~Ne(γ′,(1+z)fδD) {exp[4πif(1+z)R′cδD]−1} . (47)

Equation (5.1) implies the PSD of the synchrotron or EC flux is

 S(ϵ,f) =|~Fϵ(f)|2 =K2(1+z)2(πftlcδD)2∣∣ ∣∣~Ne(γ′,(1+z)fδD)∣∣ ∣∣2sin2(πftlc) . (48)

If the emitting region is very compact, i.e., if , then

 S(ϵ,f)≈K2(1+z)2(πftlcδD)2∣∣ ∣∣~Ne(γ′,(1+z)fδD)∣∣ ∣∣2(πftlc)2 S(ϵ,f)≈K2(1+z)2δ2D∣∣ ∣∣~Ne(γ′,(1+z)fδD)∣∣ ∣∣2 . (49)

So in the case of a very compact emitting region, the light travel time effects play no part in the PSD, as one would expect.

The observed flux PSDs from synchrotron and EC one would expect are shown in Figure 4, calculated from Equations (5.1) and (3.4). The parameter values are the same as in Figure 1, with additional parameters given in the caption. The seed photon source is assumed to be Ly photons, presumably from a broad-line region. Again, parameters are chosen to be consistent with those one would expect from an FSRQ. Synchrotron PSDs are shown for 12 Hz, 12 m (the WISE W3 filter’s central wavelength), and 0.648 m (the central wavelength of the Johnson band). EC PSDs are shown for 0.1 and 1.0 GeV, which are within the Fermi-LAT energy range. Frequencies lower than Hz can be affected by synchrotron self-absorption, which is not considered in this paper. X-rays are not shown as they are likely dominated by SSC emission, which is considered in Section 5.2 below. Also note that in real FSRQs the 12 m PSD could suffer from contamination from dust torus emission (e.g., Malmrose et al., 2011) and the band could suffer from contamination from accretion disk emission, and we do not take either of these possibilities into account.

Figure 4 shows many of the features seen in Figure 1. For photons generated from electrons with low the PSDs show a break from

 S(ϵ,f)∝f−a

to

 S(ϵ,f)∝f−(a+2)

at approximately

 f=(2πtesc)−1 ,

as seen in the Hz and 0.1 GeV PSDs. For synchrotron, if , as is usually the case for FSRQs, this regime occurs when

 ν ≪νcr,sy=1013 Hz (δDΓ)(Γ30)−3 ×(u010−3 erg cm−3)−2(tesc105 s)−2(B1 G)11+z . (50)

For EC, this regime occurs when

 mec2ϵ ≪Ecr,EC=2 GeV(δDΓ)2(Γ30)−2 ×(u010−3 erg cm−3)−2(tesc105)−2 ×(ϵ02×10−5)11+z . (51)

If then

 νcr,sy =5×1015 Hz(δD30)2(tesc105 s)−2 ×(B1 G)−311+z (52)

and

 Ecr,EC =1 TeV (δD30)2(tesc105)−2(B1 G)−4 ×(ϵ02×10−5)11+z . (53)

For photons generated from electrons with high , ( or for synchrotron or EC, respectively) minima are seen at integer multiples of , as seen in the 12 m, band and 1.0 GeV PSDs, as well as a break of 2 at approximately . There is an additional feature in Figure 4 not seen in Figure 1 related to the light travel timescale. In agreement with Equation (5.1), minima can be seen at integer multiples of which appear in the PSDs at all photon energies. Additionally, there is a break of 2 at approximately this frequency, in agreement with Equations (5.1) and (5.1).

In Figure 5 we plot observed PSDs for the parameters . This demonstrates that our model can reproduce any color noise in a synchrotron or EC PSD for an appropriate choice of , and that the features described above are preserved for different values of .

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 see by inspecting Equation (5.1) that terms associated with the light travel time will cancel when calculating time lags. However, one must be careful to shift the time lag and frequency into the observed frame by multiplying by and , respectively. An example of observed time lags can be seen in Fig. 6 for synchrotron emission and EC emission.

### 5.2. Synchrotron Self-Compton

The synchrotron-producing electrons will also Compton scatter the synchrotron radiation they produce, leading to SSC emission. Again, assuming the blob is homogeneous, and taking into account light travel time effects, we have

 FSSCϵ(t) =KSSC(1+z)tlcδD∫2R′/c0dt′ ∫min[ϵ′,ϵ′−1]0dϵ′iϵ′i ×Ne⎛⎝√ϵ′ϵ′i;tδD1+z−t′⎞⎠ Ne⎛⎝√ϵ′iϵB;tδD1+z−t′⎞⎠ (54)

where

 KSSC=δ4Dcσ2TR′uB12πd2LV′(ϵϵB)3/2 , (55)
 V′=43πR′3 (56)

is the blob volume in the comoving frame, and

 ϵ′=(1+z)ϵδD (57)

(Dermer & Menon, 2009). Following the same procedure for synchrotron and EC, for the Fourier transform we get

 ~FSSCϵ(t) ×∫∞−∞df′ ∫min[ϵ′,ϵ′−1]0dϵ′