Mid-infrared variability of γ-ray emitting blazars

ABSTRACT

Using data from the Wide-field Infrared Survey Explorer, we studied the mid-infrared (mid-IR) 3.4 μm (W1-band) and 4.6 μm (W2-band) flux variability of γ-ray emitting blazars. Our sample consists of 460 flat spectrum radio quasars (FSRQs) and 575 BL Lacertae (BL Lac) objects. On intraday time-scales, the median amplitude of variability (σ_m) for FSRQs is 0.04|$^{+0.03}_{-0.02}$| and 0.05|$^{+0.03}_{-0.02}$| mag in W1 and W2 bands. For BL Lacs, we found median σ_m in W1(W2) bands of 0.04|$^{+0.01}_{-0.02}$| (0.04|$^{+0.02}_{-0.02}$|⁠) mag. On long time-scales, for FSRQs we found a median σ_m of 0.44|$^{+0.28}_{-0.27}$| and 0.45|$^{+0.27}_{-0.27}$| mag in W1 and W2 bands, while for BL Lacs, the median values are 0.21|$^{+0.18}_{-0.12}$| and 0.22|$^{+0.18}_{-0.11}$| mag in W1 and W2 bands. From statistical tests, we found FSRQs to show larger σ_m than BL Lacs on both intraday and long time-scales. Among blazars, low synchrotron peaked sources showed larger σ_m compared to intermediate synchrotron peaked and high synchrotron peaked sources. The larger σ_m seen in FSRQs relative to BL Lacs on both intraday and long time-scales could be due to them having the most powerful relativistic jets and/or their mid-IR band coinciding with the peak of the electron energy distribution. BL Lacs have low power jets and the observational window too traces the emission from low-energy electrons, thereby leading to low σ_m. In both FSRQs and BL Lacs predominantly a bluer when brighter behaviour was observed. No correlation is found between σ_m and black hole mass.

galaxies: active, BL Lacertae objects: general, galaxies: jets, infrared: galaxies

1 INTRODUCTION

The extragalactic γ-ray sky is dominated by the blazar category of active galactic nuclei (AGNs) as evident from the Fermi-γ-ray space telescope (Atwood et al. 2009) observations since its launch in 2008 (Abdollahi et al. 2020). AGNs are believed to be powered by accretion of matter on to supermassive black holes located at the centres of galaxies (Lynden-Bell 1969; Rees 1984). Blazars, a peculiar category of AGN, which comprises both flat spectrum radio quasars (FSRQs) and BL Lacs objects (BL Lacs) emit radiation over the entire accessible electromagnetic spectrum extending from low-energy radio to high-energy TeV γ-ray energies predominantly by non-thermal emission processes. This classification of blazars into FSRQs and BL Lacs is based on the rest-frame equivalent width (EW) of their optical emission lines with BL Lacs having EW < 5 Å (Stickel et al. 1991; Stocke et al. 1991). Presence of either weak emission lines or featureless spectrum in BL Lacs is thought to be because of relativistic beaming, owing to their relativistic jets being aligned closely to the observer. However, according to Ghisellini et al. (2011), blazars can be divided into FSRQs and BL Lacs based on the luminosity of their broad emission lines (L_BLR) relative to the Eddington luminosity (L_Edd) with BL Lacs having L_BLR/L_Edd < 5 × 10⁻⁴, where |$L_{\mathrm{ Edd}} = 1.3 \times 10^{38} (\frac{M_{\mathrm{ BH}}}{M_{\odot}})$| erg s⁻¹, M_BH is the mass of the black hole. Apart from the noticeable differences in the optical spectra, both FSRQs and BL Lacs have flat radio spectra at GHz frequencies with the spectral index |$\alpha \, \lt $| 0.5 (S_ν ∝ ν^−α) and show superluminal motion in the radio band (Jorstad et al. 2005). They exhibit rapid flux variations over the electromagnetic spectrum on a range of time-scales from minutes to years (Stalin et al. 2004b; Paliya et al. 2015, 2016; Rani, Stalin & Rakshit 2017). They are highly polarized in optical which is also found to vary with time (Andruchow, Romero & Cellone 2005; Rakshit et al. 2017; Rajput et al. 2019). However, in the large-scale jet structure, FSRQs and BL Lacs do differ with FSRQs being the beamed counterparts of the luminous Fanaroff–Riley type II (FR II) radio galaxies (Fanaroff & Riley 1974) and BL Lacs being the beamed counterparts of the less luminous FR I radio galaxies.

The broad-band spectral energy distribution (SED) of blazars has a typical two hump structure. The low-energy hump peaking between infrared (IR) and X-rays is known to result from synchrotron emission process. The high-energy hump peaks in the MeV to TeV range and its origin is a matter of intense debate, and two competing models are available in the literature to explain the high-energy hump in blazars. In the one zone leptonic emission model, the high-energy hump is explained by inverse Compton (IC) process. The seed photons for the IC scattering can either originate from within the jet, called the synchrotron self-Compton (Konigl 1981; Marscher & Gear 1985; Ghisellini & Maraschi 1989) or external to the jet, called the external Compton process (Begelman et al. 1987; Melia & Konigl 1989; Dermer, Schlickeiser & Mastichiadis 1992). Alternatively, the high-energy hump can also be explained by hadronic process (Böttcher et al. 2013). Based on the peak of the synchrotron emission, blazars are further classified (Abdo et al. 2010) into low synchrotron peaked blazars (LSP; |$\nu ^{S}_{\mathrm{ peak}} \lt 10^{14}$| Hz), intermediate synchrotron peaked blazars (ISP; 10¹⁴ Hz < |$\nu ^{S}_{\mathrm{ peak}} \lt 10^{15}$| Hz), and high synchrotron peaked blazars (HSP; |$\nu ^{S}_{\mathrm{ peak}} \gt 10^{15}$| Hz).

Blazars have been studied extensively for flux variability in different wavelengths at different time-scales such as the optical (Stalin et al. 2006), X-ray (Rani et al. 2017), UV (Edelson et al. 1991; Edelson 1992), and radio (Liodakis et al. 2018). However, our knowledge on the IR variability characteristics of blazars is very limited (Kozłowski et al. 2016), though few individual sources have been studied (Carnerero et al. 2015; Zhang et al. 2015; Gabányi, Moór & Frey 2018). Recently, Mao, Zhang & Yi (2018) investigated the long-term mid-IR variability of blazars using about four years of data, however, there is no report yet in literature on their mid-IR variability characteristics on intraday time-scales. Furthermore, there is no comparative study of the mid-IR variability of the different subclasses of blazars such as LSP, ISP, and HSP available in the literature. Studies of IR variability are indeed important to understand the contribution of jet, accretion disc, and torus to the observed IR emission. As the mid-IR variability study of blazars is very limited, it is important to carryout such a study for a clear picture of their mid-IR variability. We have therefore carried out a systematic study on the mid-IR flux variability of a large sample of blazars with the following objectives (i) to characterize the mid-IR variability characteristics of γ-ray emitting blazars in general, on both intraday and long time-scales, (ii) to see for similarities and differences between the mid-IR variability characteristics of FSRQs and BL Lacs, and (iii) to have a comparative analysis of the mid-IR variability characteristics of LSP, ISP, and HSP blazars. Our sample of blazars for this study was taken from the third catalogue of AGN by Ackermann et al. (2015). This is the first systematic study of the mid-IR flux variability characteristics of γ-ray emitting blazars on intraday time-scales. We present the sample and data used in this study in Section 2 and the analysis in Section 3. The results are discussed in Section 4 followed by the summary in the final section.

2 SAMPLE AND DATA

The sample of blazars used in this study was taken from the third catalogue of AGN detected by the Fermi Large Area Telescope (3LAC, Ackermann et al. 2015). Our initial sample consists of 1099 sources of which 467 are FSRQs and 632 are BL Lacs. As the prime motivation of this work is to characterize the mid-IR variability of Fermi blazars, we searched for mid-IR counterparts to our initial sample of blazars in the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010) all sky catalogue. Since its launch and until 2012, WISE mapped the sky in four mid-IR bands namely, W1 (3.4 μm), W2 (4.6 μm), W3 (12 μm), and W4 (22 μm). Its cryogenics failed in 2012 and post-2012, WISE carried out observations in only two bands, W1 and W2. The images from WISE observations have spatial resolution of 6.1, 6.4, 6.5, and 12 arcsec, in W1, W2, W3, and W4 bands, respectively. With WISE making about 15 orbits per day, it is natural to get many photometric points on a single object in a day. The data from WISE were released in two separate catalogues namely, the AllWISE¹ source catalogue (Prior to the cryogenic failure) and NEOWISE (Near-Earth Object Wide-field Infrared Survey Explorer) catalogue² (Post Cryogenic failure). The magnitudes given in WISE are in the Vega system without any corrections for Galactic extinction. Thus, the multi-epoch photometry available in AllWISE and NEOWISE catalogues can be used to investigate the mid-IR flux variability properties of Fermi blazars.

We cross-correlated our initial sample of 1099 blazars selected from Ackermann et al. (2015) with the AllWISE source catalogue with a search radius of 2 arcsec. Our cross-correlation yielded 1035 sources. The distribution of our sample of 1035 sources in the WISE colour–colour diagram is shown in Fig. 1. These 1035 sources form our sample for mid-IR variability study. Of these 1035 sources, 575 are BL Lacs and 460 are FSRQs. FSRQs cover the redshift from 0.189 to 3.10, while BL Lacs cover the redshift between 0.034 and 1.72. The distribution of redshifts taken from Ackermann et al. (2015) for our final sample of FSRQs and BL Lacs are given in Fig. 2. In our sample, about 10 per cent of FSRQs and 50 per cent of BL Lacs do not have redshift measurements. Further dividing our sample, based on the peak of the synchrotron emission in their broad-band SED, we have 565 LSPs, 207 ISPs, and 243 HSP sources. A total of 20 sources in our sample, do not have subclassification in Ackermann et al. (2015). The summary of the different types of sources used in this study is given in Table 1.

Figure 1.

The distribution of the sources used in this study in the WISE colour–colour diagram. FSRQs (red filled circles) and BL Lacs (blue filled circles) are separated in the W1–W2 versus W2–W3 colour–colour diagram. A colour version is available online.

Figure 2.

The redshift distribution of FSRQs and BL Lacs used in this work.

Table 1.

Details of the sources used in this study.

Type	Number	z range	Median z
BL Lac	460	0.189–3.10	1.106
FSRQ	575	0.034–1.72	0.291
HSP	565	0.085–1.25	0.770
ISP	207	0.046–2.19	0.046
LSP	243	0.034–1.60	0.203

For these 1035 sources, we also looked into the availability of data in NEOWISE (Mainzer et al. 2011). Of these 1035 sources, NEOWISE data were available for 914 sources. For sources having data only in AllWISE, the observed duration spans a period of about a year between MJD 55203 and MJD 55593. However, for sources that have observations both in AllWISE and NEOWISE the duration of observations covers a period of about 7 yr between MJD 55203 and MJD 57735. Also, for observations within a day, there can be points as large as the number of orbits made by WISE, thereby enabling us to study variability on both intraday time-scales (of the order of hours) and long time-scales (of the order of years). A sample light curve of the source (a BL Lac J1748.6+7005 at z = 0.77) in W1 and W2 bands spanning about 7 yr of observation is given in Fig. 3. In Fig. 4, we present an expanded one-day light curve of the same source. For most of the sources in our sample, observations are sparse in W3 and W4 bands having ‘null’ entries at many epochs. Therefore, for any further analysis of variability, only two photometric bands were considered namely W1 and W2. Even when using W1 and W2 bands for variability analysis, to ensure use of good photometric measurements, following Rakshit et al. (2019), the following conditions were imposed

• The χ² per degree of freedom of the single-exposure profile fit in both W1 and W2 bands should be less than 5.

• The number of components used in the fit of the point spread function (PSF) of a source should be less than 3.

• The best quality single-exposure image frames are not affected by known artefacts and are not actively deblended.

• The number of data points in a day must be at least 5 in both W1 and W2 bands. This condition was utilized only for the analysis of variability on intraday time-scales.

Figure 3.

Sample light curves spanning about seven years for the BL Lac (J1748.6+7005) in W1 (red) and W2 (blue) bands. A colour version is available online.

Figure 4.

Expanded one-day light curves in W1 (red) and W2 (blue) bands for the source J1748.6+7005 shown in Fig. 3. A colour version is available online.

3 ANALYSIS OF VARIABILITY

The conditions imposed in Section 2 were used on the selection of data for variability analysis. For studying long-term variability, we used all such selected photometric points. However, for the analysis of variability on intraday time-scales we separated the photometric points into groups (hereafter called days). A one-day light curve includes all photometric points that have time gaps less than 1.2 d between any two consecutive photometric points. Also, to get rid of cosmic rays affecting our photometric data we employed a 3σ clipping to remove outliers in the one-day and long-term light curves following Rakshit et al. (2019). Thus, good quality photometric measurements were used in our analysis of mid-IR variability on intraday and long time-scales.

3.1 Variability amplitude

To characterize the variability shown by a source both on intraday and long time-scales, we calculated the amplitude of variability (σ_m) as described in Sesar et al. (2007). For each light curve, σ_m was calculated by removing the measurement uncertainty from the variance of the light curve. According to Sesar et al. (2007) σ_m is given by

$$\begin{eqnarray*} \sigma _{m} = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\sqrt{\Sigma ^2 - \epsilon ^2} & \text{if $\Sigma \gt \epsilon $,} \\ 0 & \text{otherwise} \end{array}\right. \end{eqnarray*}$$

(1)

where Σ is the variance of the light curve defined as

$$\begin{eqnarray*} \Sigma = \sqrt{\frac{1}{n-1} \sum _{n=1}^{N} (m_i - \lt m\gt)^2} \end{eqnarray*}$$

(2)

Here, m_i is the magnitude of the ith point and <m > is the weighted average. And, ϵ is defined as

$$\begin{eqnarray*} \epsilon ^{2} = \frac{1}{n} \sum _{n=1}^{N} \epsilon _i^{2}+\epsilon _s^{2} \end{eqnarray*}$$

(3)

Here, ϵ_i is the measurement uncertainty of the ith point and ϵ_s is the corresponding systematic uncertainty. For W1 and W2 bands, the systematic uncertainties are 0.024 and 0.028 mag, respectively (Jarrett et al. 2011). The systematic errors were added in quadrature to the measurement uncertainties to get the total error on each photometric measurement.

3.1.1 Intraday variability amplitude

On intraday time-scales, the variability amplitude σ_m in mag was calculated using equation (1). We found median σ_m values of 0.04|$^{+0.03}_{-0.02}$| and 0.05|$^{+0.03}_{-0.02}$| mag for FSRQs in W1 and W2 bands. Similarly, for BL Lacs, we found median σ_m values of 0.04|$^{+0.01}_{-0.02}$| and 0.04|$^{+0.02}_{-0.02}$| mag in W1 and W2 bands, respectively. The upper and lower uncertainties in the median σ_m values were determined such that 15.87 per cent of σ_m values have σ_m > σ_m(median) + σ_m(upper error) and 15.87 per cent of σ_m values have σ_m < σ_m (median) − σ_m (lower error). This corresponds to 1σ error for a Gaussian distribution. The histogram and cumulative distribution of σ_m for FSRQs and BL Lacs in W1 and W2 bands are shown in Fig. 5. The median values of σ_m in W1 and W2 bands seem indistinguishable within errors for both FSRQs and BL Lacs. However, the two sample Kolmogorov–Smirnov (KS) test indicates that there is difference in the variability between W1 and W2 bands on intraday time-scales. For FSRQs, the KS test gave a D-value of 0.113, with a null hypothesis (there is no difference in variability between W1 and W2 bands) probability (p) of 2.14 × 10⁻¹⁹, while for BL Lacs, from KS test we found a D-value of 0.213 with a p of 8.79 × 10⁻⁹⁸. Thus, on intraday time-scales, there is difference in the mean variability amplitude between W1 and W2 bands in both FSRQs and BL Lacs. From our analysis, we found FSRQs to show similar median amplitude of variability to BL Lacs in both W1 and W2 mid-IR bands. However, a two sample KS test carried out for the distribution of σ_m in the W1 band in FSRQs and BL Lacs showed that the two distributions are indeed different with a D-statistic of 0.211 and a null-hypothesis (the distribution of σ_m in W1 band for FSRQs and BL Lacs are drawn from the same population) p of 3.60 × 10⁻⁷⁹. Similarly in W2 band too, from KS test we found that the distribution of σ_m are different between FSRQs and BL Lacs with a D-statistic value of 0.121 and a p-value of 6.86 × 10⁻²⁶. A summary of the results on intraday variability analysis is included in Table 2.

Figure 5.

Histogram and cumulative distribution of σ_m on intraday time-scales for FSRQs (dashed line) and BL Lacs (solid line) in W1 (left-hand panel) and W2 (right-hand panel) bands.

Table 2.

Results of the analysis of variability.

Type	Number	σm ± σ (intraday time-scale)		σm ± σ (long time-scale)		Duty cycle
		W1 (mag)	W2 (mag)	W1 (mag)	W2 (mag)	W1 (per cent)	W2 (per cent)
FSRQ	460	0.04|$^{+0.03}_{-0.02}$|	0.05|$^{+0.03}_{-0.02}$|	0.44|$^{+0.28}_{-0.27}$|	0.45|$^{+0.27}_{-0.27}$|	78.82	53.06
BL Lac	575	0.04|$^{+0.01}_{-0.02}$|	0.04|$^{+0.02}_{-0.02}$|	0.21|$^{+0.18}_{-0.12}$|	0.22|$^{+0.18}_{-0.11}$|	89.59	61.47
LSP	565	0.05|$^{+0.05}_{-0.03}$|	0.05|$^{+0.05}_{-0.03}$|	0.40|$^{+0.30}_{-0.22}$|	0.41|$^{+0.29}_{-0.22}$|	79.51	54.12
ISP	207	0.04|$^{+0.04}_{-0.02}$|	0.04|$^{+0.05}_{-0.02}$|	0.22|$^{+0.19}_{-0.11}$|	0.24|$^{+0.21}_{-0.12}$|	92.79	82.29
HSP	243	0.04|$^{+0.03}_{-0.02}$|	0.04|$^{+0.04}_{-0.02}$|	0.14|$^{+0.11}_{-0.08}$|	0.16|$^{+0.10}_{-0.09}$|	82.89	43.71

3.1.2 Long-term variability amplitude

The distribution of long-term variability amplitudes in W1 and W2 bands for FSRQs and BL Lacs are shown in Fig. 6. For FSRQs in W1 band, we found σ_m to range between 0.020 and 1.550 mag with a median value of 0.44|$^{+0.28}_{-0.27}$| mag. Similarly for W2 band we found σ_m to range between 0.016 and 1.186 mag with a median value of 0.45|$^{+0.27}_{-0.27}$| mag. From two sample KS test applied to the distribution of σ_m between W1 and W2 bands, we found a D-statistic of 0.044 with a p-value of 0.90. Thus in FSRQs, there is no difference in variability between W1 and W2 bands. In the case of BL Lacs, in W1, we found σ_m to range between 0.008 and 1.350 mag with a median of 0.21|$^{+0.18}_{-0.12}$| mag. Similarly in W2 band, we found σ_m to lie in the range between 0.010 and 1.249 mag with a median of 0.22|$^{+0.18}_{-0.11}$| mag. A two sample KS test to the distribution of σ_m between W1 and W2 bands in BL Lacs gave a D-statistic of 0.070 with a p-value of 0.210. Thus in long term, the amplitude of flux variations between W1 and W2 bands are found to be similar in both FSRQs and BL Lacs. A two sample KS test applied to the distribution of σ_m in W1(W2) bands between FSRQs and BL Lacs yielded a D-statistics of 0.42 (0.43) and a p-value of 1.10 × 10⁻³³ (2.23 × 10⁻²⁸). Thus, on long time-scales too, FSRQs showed larger amplitude variations than BL Lacs in both W1 and W2 bands. A summary of the results of the long-term variability analysis is included in Table 2.

Figure 6.

Histogram and cumulative distribution of σ_m on long time-scales for FSRQs (dashed line) and BL Lacs (solid line) in both W1 (left-hand panel) and W2 (right-hand panel) bands.

3.2 Flux variability on subsamples of blazars

We also divided our sample of blazars into different subclasses based on the peak frequency of the synchrotron component in their broad-band SED such as LSP, ISP, and HSP and analysed the amplitude of variability in them both on intraday and long time-scales. The results are included in Table 2 for both intraday and long time-scales, respectively. The histogram and cumulative distribution of σ_m for different subclasses of blazars are given in Fig. 7 for intraday time-scales and Fig. 8 for long time-scales. On long time-scales LSP blazars showed the largest σ_m followed by ISP and HSP sources in both W1 and W2 bands. This is evident from the cumulative distribution function (CDF) plots in the bottom panels of Fig. 8, where the CDFs of LSPs are systematically at higher σ_m values than those of ISPs and HSPs, while the CDFs of ISPs are larger than HSPs. This is expected from the results in the previous section as majority of LSP sources are FSRQs, while most of HSP sources are BL Lacs. On intraday time-scales, while LSP sources are more variable compared to ISP and HSP sources, the amplitude of variability is indistinguishable between ISP and HSP sources. Here too, in the CDFs shown in the bottom panels of Fig. 7, LSPs have higher σ_m values than ISPs and HSPs, while the CDFs of ISPs and HSPs are indistinguishable.

Figure 7.

Distribution of σ_m on intraday time-scales for LSP, ISP, and HSP objects.

Figure 8.

Histogram and cumulative distribution of σ_m on long time-scales for the different subclasses of blazars in both W1 and W2 bands.

3.3 Ensemble structure function

In addition to the variability of blazars characterized by σ_m described in the previous sections, blazar variability can also be described by structure function (SF), which shows the dependency of variability as a function of time lag (Simonetti, Cordes & Heeschen 1985). SF can be calculated for individual AGN having light curve with multiple epochs of observations, taking the magnitude difference for each pair of time lags in a light curve. It can also be calculated for a group of blazars, known as ensemble SF, allowing us to obtain the mean variability behaviours of the population. Since WISE data have sparse sampling, we studied the mean variability of different classes of blazars using the ensemble SF following (di Clemente et al. 1996)

$$\begin{eqnarray*} \mathrm{SF} = \sqrt{\frac{\pi }{2} \langle|\Delta m|\rangle ^2 - \langle\sigma ^2_n\rangle}, \end{eqnarray*}$$

(4)

where |Δm| = m_i − m_j, is the magnitude difference of any two epochs (i, j) separated by the time difference or lag Δτ = t_i − t_j. |$\sigma ^2_n=\sigma ^2_i+\sigma ^2_j$| is the sum of squares of the uncertainties in magnitudes m_i and m_j. To calculate SF, the time bins were selected to have equal interval in logarithmic observer frame time lag with at least 2000 data points in each bin. Note that a few of our sources do not have redshift measurements in the literature, thus, SF was calculated in the observer frame.

In Fig. 9, we plot SF against observer frame time lag for BL Lacs (red) and FSRQs (blue). The error bars in SF were calculated via error propagation following Vanden Berk et al. (2004). The figure shows FSRQs are significantly more variable than BL Lacs consistent with the result obtained by σ_m calculation. SF increases gradually from a time lag of 1 to ∼200 d and become flatter at higher time lag. Such a trend has been noted previously by various authors (e.g. Vanden Berk et al. 2004; Welsh, Wheatley & Neil 2011; Kozłowski 2016) from studies of long-term quasar variability in the optical band. To characterize the SF, we fitted a simple power-law model in the form of

$$\begin{eqnarray*} \mathrm{SF} = \left(\frac{\Delta \tau }{\tau _0}\right)^{\gamma }, \end{eqnarray*}$$

(5)

where τ₀ and the power-law slope (γ) are free parameters. The fitting results are given in Table 3. The inferred value of γ is 0.29 ± 0.03 for BL Lacs and 0.25 ± 0.02 for FSRQs. We also fitted a three parameters exponential model following Kozłowski (2016) in the form of

$$\begin{eqnarray*} \mathrm{SF} = \mathrm{SF_{\infty }} \sqrt{ 1- e^{- \left(\frac{|\Delta \tau |}{\tau _c} \right)^{\beta }}}, \end{eqnarray*}$$

(6)

where SF_∞ is the amplitude of variance at long time-scale, τ_c is decorrelation time-scale, and β is the power-law index. The equation becomes a damped random walk (DRW; Kelly, Bechtold & Siemiginowska 2009) for β = 1. We note that the decorrelation time-scales could be affected by the limited length of the light curves as shown from simulations by Emmanoulopoulos, McHardy & Uttley (2010).

Figure 9.

SF against observer frame time lag for BL Lacs (red dots) and FSRQs (blue dots). Best fits of the SF using equation (2, dashed line) and equation (3, solid line) are also shown. The 3.6 μm SF for quasars from Kozłowski et al. (2016) is shown with a black dashed line. A colour version is available online.

Table 3.

Results of model fits to the ensemble SF for various classes of blazars.

	Power-law model		Exponential model
Object class	γ	τ0 (d)	SF∞ (mag)	β	τc (d)
BL Lac	0.29 ± 0.03	1.19 ± 0.52 × 104	0.38 ± 0.02	0.89 ± 0.09	128 ± 46
FSRQ	0.25 ± 0.02	2.25 ± 0.63 × 103	0.66 ± 0.01	0.82 ± 0.03	90 ± 12
HSP	0.30 ± 0.05	3.81 ± 3.26 × 104	0.30 ± 0.07	0.75 ± 0.14	327 ± 410
ISP	0.32 ± 0.02	1.23 ± 0.28 × 104	0.40 ± 0.04	0.79 ± 0.05	375 ± 167
LSP	0.27 ± 0.04	3.01 ± 1.29 × 103	0.57 ± 0.03	0.96 ± 0.13	89 ± 39

Equation (6) provides a better fit to the data than simple power law defined in equation (5) especially at lower lag. The fitted parameters are given in Table 3. We found SF_∞ = 0.38 ± 0.02 for BL Lacs and 0.66 ± 0.01 for FSRQs suggesting a higher variability in FSRQs than BL Lacs. The value of β is about ∼0.85 slightly deviating from DRW. We also plot in the same figure, the 3.6 μm SF of confirmed quasars from Kozłowski et al. (2016) in black dashed line, which lies well below that of the BL Lacs and FSRQs studied here. This could be due to the contribution of relativistic jets to the mid-IR variability of the blazars studied here in comparison to the sample of quasars studied by Kozłowski et al. (2016). In Fig. 10, we plot the SFs of HSP, ISP, and LSP sources. Here too, we found a better fit with the exponential function given in equation (6). We found SF_∞ = 0.30 ± 0.07, 0.40 ± 0.04, and 0.57 ± 0.03 for HSPs, ISPs, and LSPs.

Figure 10.

SFs for HSP (red), ISP (black), and LSP (blue). Lines have the same meaning as in Fig. 9. A colour version is available online.

The SF plots in Fig. 10 show that variability is significantly stronger in LSPs than in ISPs and variability is the lowest in HSPs. Among ISPs and HSPs, variability is stronger in ISPs. This is in agreement with the results obtained from the analysis of the amplitude of flux variations.

3.4 Duty cycle of variability

The data set analysed here for intraday variability is suited to estimate the duty cycle (DC) of mid-IR variability of different classes of blazars. Following Romero, Cellone & Combi (1999), the DC of a particular class of object is given as

$$\begin{eqnarray*} \mathrm{ DC} = 100 \frac{\sum _{j=1}^{n} A_j \times (1/\Delta t_j)}{\sum _{j=1}^{n} (1/\Delta t_j)} {{\ \rm per\ cent}} \end{eqnarray*}$$

(7)

Here, Δt_j = Δt_{j, obs}/(1 + z) is the duration of observation corrected for the cosmological redshift of the source of a jth one-day light curve (as described in Section 3) for a selected class. A_j = 1 or 0, if the source was classified as a variable or non-variable respectively during the observed duration Δt_j. A source is considered variable in a one-day light curve (A_j = 1) if σ_m calculated using equation (1) is greater than zero. For BL Lacs in the W1 band, we found a high DC of 89.59 per cent, while for FSRQs we found a lower DC of variability of 78.82 per cent. In the W2 band too, BL Lacs have a high DC of 61.47 per cent followed by FSRQs of 53.06 per cent. Separating the blazars into different classes based on their SED, we found the ISP sources to have the highest DC of variability and HSPs have the lowest DC of variability in both W1 and W2 bands. The results of the DC of variability in various subclasses of blazars is given in Table 2. In the optical band on intranight time-scales, Stalin et al. (2004a) found the BL Lacs to have a high DC of variability related to FSRQs which was interpreted as due to the close alignment of BL Lac jets with the line of sight related to FSRQs (Stalin et al. 2004a). The results found here on DC of the mid-IR variability on intraday time-scales closely matches with that found by Stalin et al. (2004a) on intranight time-scales in the optical band in spite of differences in the time resolution between the optical and the mid-IR light curves.

3.5 Colour variability

Blazars are known to show spectral variations in the optical band. It has been thought that FSRQs show a redder when brighter behaviour (RWB; Gu et al. 2006; Bonning et al. 2012). Alternatively, BL Lacs are found to show a bluer when brighter behaviour (BWB; Massaro et al. 1998; Villata et al. 2002; Vagnetti, Trevese & Nesci 2003; Gaur et al. 2012). Departures from this conventional observations have also been noted recently. FSRQs are found to show BWB behaviour (Gu & Ai 2011) and in the FSRQ 3C 345 both RWB and BWB trends were noticed (Wu et al. 2011). There are also reports in which the spectrum of a blazar was found not to change with increasing/decreasing brightness (Stalin et al. 2006).

The available literature is more focused on the colour variability of blazars in the optical and near-IR bands (Bonning et al. 2012; Gaur et al. 2019; Sarkar et al. 2019), but not in the mid-IR bands. Most of these studies were based on nearly quasi-simultaneous observations, without properly taking into account the errors in both colours and magnitudes which could lead to incorrect characterization of spectral variability (Sukanya et al. 2016). In this work, we report on the mid-IR spectral variations in a large sample of blazars.

3.5.1 Intraday colour changes

The advantage of colour variability studies in mid-IR using WISE is that the observations in the different bands are simultaneous. To characterize the colour variability in W1 and W2 bands, we constructed colour–magnitude diagrams, wherein (W1–W2) colour is plotted along the Y-axis and the W1 brightness is plotted along the X-axis. We carried out a linear least-squares fit to the colour–magnitude diagram by taking into account the errors in both the colour and magnitude. The slope of the fit is taken to quantify the spectral change. We used the Spearman rank correlation analysis to probe the correlation between W1–W2 colour against the W1 brightness. The source becomes increasing bluer with increasing W1–W2 colour along the Y-axis and W1 increasing (decreasing in brightness) towards the right. This is the BWB trend as seen in Fig. 11 (right-hand panel). Alternatively, the situation wherein the W1–W2 colour gets smaller with increasing W1 band (decreasing in brightness) is the RWB trend. An example of such a trend is shown in the left-hand panel of Fig. 11. In this work, we adopted the following criteria to characterize the colour variations in blazars. We considered a source to show a BWB trend if the Spearman rank correlation coefficient is larger than 0.3 and probability of no correlation (p) is less than 0.05. Similarly we considered a source to show a RWB trend if the Spearman rank correlation coefficient is less than −0.3 and p is less than 0.05. In the middle top panel of Fig. 11 is shown the distribution of the slopes obtained from linear least-squares fit to the colour–magnitude diagram on intraday time-scales. Clearly the blazars in our sample showed all types of spectral behaviours namely (i) constancy of spectral shape with brightness, (ii) RWB behaviour, and (iii) BWB behaviour. However, the distribution is shifted from zero to positive values (there are more positive than negative values) thereby indicating that on intraday time-scales most of the blazars showed a BWB trend. On intraday time-scales, in BL Lacs, for 1409 light curves, we found a BWB trend, while only on four light curves, we noticed a RWB trend. For FSRQs, 960 light curves showed a BWB trend, while only on five light curves RWB trend was noticed. Thus, on intraday time-scales, we found both FSRQs and BL Lacs predominantly showed a BWB trend, while only on very few instances RWB trend was noticed.

Figure 11.

Distribution of the slopes of the linear least-squares fits done to the colour–magnitude diagram of FSRQs (top middle panel) and BL Lacs (bottom middle panel) on intraday time-scales. The left- and right-hand panels on the top show a typical example of RWB and BWB trends, respectively in FSRQs, while the bottom left and bottom right panels show a sample RWB and BWB trend in BL Lacs, respectively. A colour version is available online.

3.5.2 Long-term colour changes

Using data that spans from about a year to as long as 7 yr, we also studied the colour changes in the long time-scales. Here too, colour–magnitude diagrams were constructed, and linear least-squares fit were carried out to the colour–magnitude diagrams by taking into account the errors in both the colours and magnitudes. The results of the spectral fits are shown in Fig. 12. Similar to the spectral trends noticed on intraday light curves, here too, we found all types of spectral behaviours namely, (i) the spectral shape has not changed with brightness, (ii) the spectrum has become BWB, and (iii) the spectrum has become RWB. An example of a BWB character is shown in the bottom right panel of Fig. 12, while an example of an RWB trend is shown on the bottom left hand panel of Fig. 12. The distribution of the slopes of the spectral fits are shown in the middle panels of Fig. 12. Similar to the spectral variations observed from intraday light curves, here too, the distribution of spectral slopes have more positive than negative values, pointing to the BWB trend in most of the blazars. On long time-scales, we found 174 FSRQs to show a BWB trend, while RWB trend was found in 24 sources. Similarly, for BL Lacs, 160 objects showed a BWB trend, and a small number of 47 objects showed an RWB trend. Thus, on long time-scales majority of FSRQs and BL Lacs showed a BWB spectral behaviour.

Figure 12.

The top middle and the bottom middle panels show the distribution of slopes obtained from the linear least-squares fit to the colour–magnitude diagram on long time-scales for FSRQs and BL Lacs. The top left and top right panels show sample RWB and BWB trend, respectively in FSRQs, while the same for BL Lacs is shown in the bottom left and bottom right panels, respectively. A colour version is available online.

3.6 Correlation of long time-scale σ_m with M_BH

In order to find the correlation between the black hole mass (M_BH) and long time-scale σ_m, we collected M_BH values from the catalogue of virial black hole mass estimates obtained by Kozłowski et al. (2016). We could get M_BH values for a total of 67 objects. For those objects, we checked for the correlation between M_BH and σ_m in both W1 and W2 bands, by applying a linear least-squares fit to the logarithmic variability amplitude σ_m as a function of the logarithmic black hole mass. This is shown in Fig. 13.

Figure 13.

Correlation of σ_m with M_BH. The solid and dashed lines are unweighted linear least-squares fit to the data on FSRQs and BL Lacs, respectively.

For FSRQs in the W1 band, the Spearman rank correlation coefficient is −0.15 with a p-value of 0.11 and for the W2 band, the correlation coefficient was found to be −0.20 with a p-value of 0.06. For BL Lacs we found the correlation coefficient of −0.23 with a p of 0.41 in W1 band, while for W2 band we found a correlation coefficient and p of −0.17 and 0.59, respectively. Thus, in our data, we did not find statistically significant correlation between M_BH and mid-IR amplitude of variability. This is the first investigation of the correlation of σ_m with M_BH for blazars in the mid-IR bands, however, such correlation analysis are already available in other wavelengths in other types of AGN. For example in the optical band for quasars, there is no conclusive evidence on the dependence of variability amplitude with M_BH as there are reports of positive (Wold, Brotherton & Shang 2007; Lu et al. 2019), negative (Kelly et al. 2009), and no correlation (Simm et al. 2016) on time-scales similar to the one probed here in the mid-IR band. Rakshit et al. (2019) studied mid-IR variability of narrow-line Seyfert 1 galaxies and found no correlation between variability amplitude and M_BH.

4 DISCUSSION

4.1 Flux variability

The observations from WISE, that makes about 15 orbits per day is well suited to study the mid-IR flux variations on intraday as well as long time-scales. We utilized this data set to probe the mid-IR variability characteristics of a sample of blazars selected from 3FGL. Our analysis is a first characterization of the mid-IR variability of γ-ray emitting blazars. We found that on intraday time-scales FSRQs showed larger amplitude variability than BL Lacs. In the optical band on time-scales much shorter than a day, BL Lac objects were found to show large amplitude and high DC of variability compared to FSRQs. These observations in the optical band are explained in the context of BL Lacs having closely aligned jets relative to FSRQs (Stalin et al. 2004a). The differences between the optical and mid-IR intraday variability characteristics of blazars could be ascribed to the following: (i) the optical analysis of Stalin et al. (2004a) was based on a limited number of 15 blazars, with each blazar having three or four intraday light curves. The mid-IR variability obtained here is based on the analysis of a larger number of blazars with each blazar having many (>3) intraday light curves and (ii) the time resolution of the optical light curves analysed by Stalin et al. (2004a) is of the order of minutes, while, the time resolution of the mid-IR intraday light curves analysed here is of the order of hours. Mid-IR light curves with time resolution of the order of minutes are needed to make a direct comparison to the results reported by Stalin et al. (2004a) in the optical band.

Thermal IR radiation in AGN is generally believed to be from the accretion disc and the dusty torus outside of the dust sublimation radius having typical scales of around few tenths of parsec (Koshida et al. 2014; Mandal et al. 2018). Dust reverberation mapping observations show that the inner radius (R) of the torus ranges from 0.01 to 0.1 pc (Suganuma et al. 2006). Using this measured size of the torus, one can estimate the variability time-scale (Δt_int) as

$$\begin{eqnarray*} R \lt c \times \Delta t_{\mathrm{ int}} \end{eqnarray*}$$

(8)

Here, Δt_int = Δt_obs/(1 + z) and c is the speed of light. The FSRQs studied here have redshifts between 0.19 and 3.10 and for those sources, using equation (8) and considering an R of 0.1 pc, the observed time-scale of variability ranges from about 140 d (for z = 0.19) to about 480 d (for z = 3.10). In our analysis of intraday variability, the maximum time difference between the data points is 1.2 d. We have found variations using such observations that span about 1.2 d and therefore our observed time-scales of variability are much smaller than 1.2 d. This implies that the observed mid-IR flux variations on intraday time-scales are not from the torus. Using the maximum time difference between the data points (i.e. 1.2 d) as the observed variability time-scale, from equation (8), the mid-IR emission could have come from spatial scales ranging from 2 × 10⁻⁴ (for a source at z = 0.19) to 8 × 10⁻⁴ pc (for a source at z = 3.10). This suggests that the observed mid-IR emission comes from a very compact region corresponding to around few hundreds of Schwarzschild radii considering FSRQs to have black hole masses around 10⁸ M_⊙. The size of the emission region as deduced from from the observed intraday variability time-scale is much smaller than that of the torus and the accretion disc. This suggests that the observed radiation in W1 and W2 bands is dominated by the non-thermal emission from the jets of these sources. This is not surprising, as these sources are γ-ray emitting sources and detected by Fermi. In this scenario, the observed mid-IR variability is well explained by jet-based models (Marscher & Gear 1985; Marscher 2014).

Statistical tests indicate FSRQs to be more variable than BL Lacs in both the W1 and W2 band on intraday and long time-scales. Analysis of the ensemble variability using SF on long time-scales shows differences in variability amplitude between FSRQs and BL Lacs, which confirms the results obtained through amplitude of variability method. Dividing the sample of sources based on the position of the synchrotron peak in their broad-band SED, SF analysis and the analysis of variability amplitude too indicates that LSPs show mid-IR variability with the largest amplitude, followed by ISPs and HSPs. The increased variability of FSRQs relative to BL Lacs and LSPs relative to ISPs and HSPs could be due to a combination of the following two reasons. First, the observed W1 and W2 bands trace the peak of the synchrotron component in the case of FSRQs, while in the case of BL Lacs, it traces the rising part of the synchrotron component and thus the low-energy electron population. Secondly, FSRQs are known to have powerful jets compared to BL Lacs (Gardner & Done 2018). As most of the FSRQs are LSP sources and BL Lacs are HSP sources, LSPs show large amplitude variability compared to ISPs and HSPs. Our results on the mid-IR amplitude of variability on long time-scales is in agreement with what is found in the optical band. From an analysis of optical data in long term, Bauer et al. (2009) and Hovatta et al. (2014) noticed FSRQs to be more variable than BL Lacs. This too points to enhanced contribution of jet emission in the optical band of these sources.

4.2 Colour variability

Brightness variations in blazars are often accompanied by spectral variations that manifest in colour–magnitude correlation plots or spectral index–magnitude correlation plots. Analysis of spectral variations too are important as it can provide additional clues to the origin of flux variations. Blazars have been studied for spectral variations in the optical and near-IR bands (Bonning et al. 2012; Gaur et al. 2019; Sarkar et al. 2019), however, such studies in the mid-IR band are very limited. For example, from an analysis of CTA 102 mid-IR (W1 and W2) light curves during its optical outburst in 2016 (Kaur & Baliyan 2018), Jiang (2018) noticed a BWB trend within a day. The data set analysed here is unique (larger number of sources and many epochs) and it can provide insights on a statistical basis into the mid-IR colour variations in blazars.

The observed emission in W1 and W2 bands is a combination of thermal emission from the accretion disc and torus and non-thermal synchrotron emission from the relativistic jet (Kozłowski et al. 2016). For most of the γ-ray bright blazars, the bright state in the γ-ray band is accompanied by correlated increase in brightness in the optical as well as the IR bands (Bonning et al. 2012). However, there are exceptions to this general observations known in blazars. There are cases in blazars where the optical, near-IR, and γ-ray emissions are not correlated (Chatterjee et al. 2013; Liodakis et al. 2019; Rajput et al. 2019). From a systematic study of long-term optical flux variations in γ-ray-detected blazars along with a control sample of gamma-ray undetected blazars, it has been found that γ-ray-detected blazars are more variable than γ-ray undetected blazars (Hovatta et al. 2014). Thus, it is clear that the observed flux variations (across the electromagnetic spectrum), both on intraday and long time-scales in γ-ray bright blazars is due to the relativistic jets in them. The observed BWB trend seen in this study can be explained by (i) localized temperature changes in the accretion disc with changes in the accretion rate (Ruan et al. 2014), (ii) increased amplitude of variability at shorter wavelengths (Stalin et al. 2009) which in the one zone synchrotron emission model could happen due to the injection of fresh electrons that have an energy distribution harder than that of the earlier cooled electrons (Kirk, Rieger & Mastichiadis 1998; Mastichiadis & Kirk 2002) and (iii) changes in the Doppler factor in a convex spectrum (Villata et al. 2004) or variation in Doppler factor due to changes in the viewing angle of a curved and inhomogeneous jet (Papadakis, Villata & Raiteri 2007). Considering that the observed flux variations are intrinsic to the source, we rule out geometric effects on the cause of colour variations and instead focus on the intrinsic causes for the observed colour variations. Among blazars, FSRQs have a stronger accretion disc and more powerful jets relative to BL Lacs that have weak accretion disc and less powerful jets (Gardner & Done 2018). In both FSRQs and BL Lacs, predominantly BWB trend is observed and therefore, spectral variations with a BWB trend that could result from temperature changes in the accretion disc with changes in the accretion rate is disfavoured. Flux variability analysis discussed in Section 4.1 unambiguously point to the jet-based origin of the observed IR variability and therefore, the observed colour variations are related to complex processes intrinsic to the jets of both FSRQs and BL Lacs.

5 SUMMARY

Using a large sample of 1035 blazars taken from the third catalogue of AGN detected by Fermi, and cross-matched with the WISE catalogue, we studied the mid-IR variability properties of FSRQs and BL Lacs on both intraday and long time-scales. While blazars have been studied for mid-IR variability on long time-scales (Mao et al. 2018), the study presented here is the first on the mid-IR intraday variability characteristics of different categories of blazars using an extensive data set taken from WISE observations. We quantified variability by calculating the amplitude of variability, σ_m. In addition to flux variability, we also studied the IR colour variations in our sample of sources. The major findings of this present study are summarized below

• All sources in our sample, except three showed flux variations both on intraday and long time-scales.

• On intraday time-scales, we found FSRQs to show larger amplitude flux variations in the mid-IR W1 and W2 bands relative to BL Lac objects. When the sample is divided into different subclasses based on the position of the synchrotron peak in their broad-band SEDs, LSPs showed the largest amplitude of variability, while HSPs and ISPs showed similar variability amplitudes.

• On long time-scales, FSRQs showed large amplitude flux variations compared to BL Lacs in both W1 and W2 bands. However, there is no difference in variability amplitudes between W1 and W2 bands in FSRQs and BL Lacs. Among the various subclasses of blazars, in W1 and W2 bands, LSPs showed the largest amplitude of flux variability and HSPs showed the lowest amplitude of flux variations, while ISP sources sources showed flux variations with the amplitude of variability intermediate between LSPs and HSPs.

• Most of FSRQs and BL Lacs showed a BWB trend, while only a small fraction of them showed an RWB behaviour.

• No correlation was found between the mid-IR amplitude of variability and black hole mass in both FSRQs and BL Lacs.

• From the analysis of intraday light curves, we found BL Lacs to show an increased DC of variability than FSRQs in both W1 and W2 bands.

ACKNOWLEDGEMENTS

We thank the anonymous referee for his/her critical review of our manuscript that helped to improve the presentation significantly. This publication makes use of data products from the WISE, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration.

Footnotes

1

http://wise2.ipac.caltech.edu/docs/release/allwise/

2

http://wise2.ipac.caltech.edu/docs/release/neowise/