Microvariability in BL Lacertae: “Zooming” into the Innermost Blazar Regions

Abstract

In this work, we present the results of our multi-band microvariability study of the famous blazar BL Lac. We performed microvariablity observations of the source in the optical VRI bands for four nights in 2016. We studied the intranight flux and spectral variability of the source in detail with an objective to characterize microvariability in the blazars, a frequently observed phenomenon in blazars. The results show that the source often displays a fast flux variability with an amplitude as large as ~0.2 magnitude within a few hours, and that the color variability in the similar time scales can be characterized as “bluer-when-brighter” trend. We also observed markedly curved optical spectrum during one of the nights. Furthermore, the correlation between multi-band emission shows that in general the emission in all the bands are highly correlated; and in one of the nights V band emission was found to lead the I band emission by ~13 min. The search for characteristic timescale using z-transformed auto-correlation function and the structure function analyses reveals characteristic timescale of ~50 min in one of the R band observations. We try to explain the observed results in the context of the passage of shock waves through the relativistic outflows in blazars.

1. Introduction

Blazars, a subclass of radio galaxies with their jets pointing close to the earth, are known to exhibit some of the extreme properties such as high luminosity, rapid flux and polarization variability and broadband non-thermal emission. These extreme properties are often attributed to the Doppler boosted emission from the relativistic outflows emanating from the region close to the central engine (e.g., [1]). One of the most important properties that characterize blazars is their variability over a broad range of temporal and spatial frequencies. The variability in general is found to be aperiodic in nature; however, recently, detection of quasi-periodic oscillations in various emission frequencies and timescales have been reported (see [2,3,4]). In particular, low amplitude, rapid variability in the timescales less than a day is widely known as intraday/intranight or microvariability. The terms microvariability and intranight usually refer to the intraday variability as observed from the ground based optical telescopes. Following causality argument, the emission regions where microvariability is produced should be highly compact (sub-)volumes, most likely, close to the central engine. Such compact emission zones cannot be spatially resolved by any current instruments, therefore multifrequency microvariability studies could be one of the most powerful tools that allows us to “zoom” into the innermost regions at the interface of the supermssive black hole and the base of the blazar jets; and thereby constrain the nature of prevalent physical processes, e.g., particle acceleration and energy dissipation mechanism, magnetic field geometry, jet content, etc., that contribute to the origin of microvariability in radio-loud (RL) AGNs. To explain the phenomenon several scenarios have been proposed most of which relate the source of the microvariability to a broad range of possible physical processes on-going both in the accretion disk and the jet. The scenarios include emission regions revolving around the central source, variable obscuration, various magnetohydrodynamic instabilities, propagation of the shock waves down the turbulent jet, projection effects resulting from the orientation of the relativistic jets relative to the line of sight (e.g., [5,6,7,8]). However, in spite of substantial observational efforts accompanied by comprehensive theoretical discussion on the subject, the details of the underlying processes responsible for microvariabity, and in general variability, still remain elusive.

Optical variability in blazars on intranight timescales has been extensively studied over the past several decades (for earlier review, see [7]). In a sample of flat-spectrum radio sources that have a compact Very Long Baseline Interferometry (VLBI) structure, Quirrenbach et al. [9] detected intraday variations with amplitudes ranging up to 25%. Similarly, Heidt & Wagner [10] presented the results of the microvariability study in a large sample of radio-selected BL Lac objects which reported the detection of microvariability in the sources with a large (~0.8) duty cycle (DC) and a typical peak-to-peak variability amplitude as large as ~30%.

In a study of several BL Lac objects consisting of the radio-selected BL Lac objects (RBLs) and X-ray-selected BL Lac objects (XBLs), Bai et al. [11] found that XBLs show more frequent microvariability than RBLs. Kraus et al. [12] observed that intraday variations are common in RL AGNs and that the flux variations are often accompanied by similar variations in the linear polarization. Similarly, during the study of long-term light curves of four blazars, including the source BL Lac. Howard et al. [13] found that the occurrences of microvariability were temporally correlated with long-term optical activity of the source. In other words, microvariability is correlated with the flux gradients rather than specific flux states. Rani et al. [14] reported the detection of microvariability with large DC in a sample of bright low-energy-peaked BL Lac sources (LBLs). Similarly, during their study of multiband optical flux and color variations of blazars on intraday and short term timescales of a few months, Gaur et al. [15] found that the spectra of BL Lac sources often showed bluer-when-brighter trend whereas in flat spectrum radio quasars (FSRQ) the redder-when-brighter trend appears more frequent. Bachev et al. [16] found evidence of quasi-periodic oscillations at very low-amplitude levels during the short-term optical monitoring program of 13 blazars. Furthermore, symmetry analysis of the longterm optical light curves of the blazar BL Lac resulted a positive correlation between flare durations and peak fluxes [17]. As an attempt to characterize microvariability in the famous BL Lac source S5 0716+714 robustly, the source was intensely studied using extended whole earth blazar telescope (WEBT) campaigns (see [8,18,19,20]). Earlier multiband intranight variability in the source S5 0716+714 has been studied by several authors (e.g., [21,22]). More recently, in the hard X-ray regime the properties of intraday variability of a sample of radio-loud AGNs are discussed by Bhatta et al. [23].

BL Lacertae (RA = 22^h02^m 43^s.3, Dec = + 42${}^{\circ }$16’40” and z = 0.0686), the eponymous blazar originally mistaken for a star, is one of the widely studied blazars across a wide range of electromagnetic frequencies using most of the currently available instruments. Because of its highly pronounced optical variability, it has been one of the most favorable targets for a large number of multifrequency studies (e.g., [24,25,26,27,28]). On the intranight timescales, the source was found to show rapid variability [29]; and on similar timescales, rapid flux and color variations with a trend to become bluer when brighter have been reported by several authors (e.g., [30,31,32,33,34,35,36]). The optical power spectral density (PSD) in the intranight timescale can be characterized as red noise behavior with a slope index of 2 [37]. Cross-correlation between longterm optical light curve and radio hardness ratio revealed a radio time delay of more than 3 months (see [38]). The VLBA images showed superluminal features with the apparent speeds ranging from 3.9–13.5c [39]. In addition, the Very Energetic Radiation Imaging Telescope Array System (VERITAS) observations of the source revealed one of the most rapid TeV $\gamma$-ray flares with decay timescale of ~13 min accompanied by the emergence of a superluminal component from the VLBA core as well as changes in the optical polarization angle [40].

In this paper, we present the results of our microvariability study of the blazar BL Lac in VRI photometric bands during four nights in 2016. The work represents a part of long-term and broader effort to characterize the nature of microvariablity in a sample of radio-loud AGNs (for detail see [41]). We organize our presentations in the following way: In Section 2, the optical observations and the data processing are discussed. In Section 3, we elaborate on our flux and spectral variability analysis that includes cross-correlation study on the light curves along with auto-correlation and structure function analyses. Finally, we present our discussions and conclusions in Section 4.

2. Observations and Data Processing

Most of the microvariablity observations utilized in the work were acquired with the Southeastern Association for Research in Astronomy (SARA) suite of telescopes which includes the SARA 0.9-m at Kitt Peak and the SARA Jacobus Kapteyn Telescope (JKT) at La Palma Canary islands. To extract the source brightness in magnitudes along with related uncertainties, the images were reduced and processed using the software MIRA1. The standard procedures for aperture photometry were followed after the images were corrected for bias, dark, and flat-fielding. Apertures of about 2–4 arcseconds were chosen so as to minimize the flux scatter in the comparison stars in the same CCD frame. The comparison stars were selected from the Heidelberg AGN finder chart2. Star B in the field was excluded as we suspected it to be slightly variable. To determine if the microvariability is real or instrumental in nature, we also routinely monitor the flux of the comparison stars; and to minimize the sky effects microvariability observations are performed on the nights with stable photometric conditions such as seeing. The test for microvariability was performed by comparing the flux changes in the source with the flux changes in the comparison stars on the same CCD frame and following the statistical procedure discussed in [42].

We have been monitoring the source since past several years. However, in this paper, we particularly discuss its variability properties in the year 2016. We observed the source continuously for a few hours in the nights dated 7 August, 12 August, 30 September and 5 October 2016 (hereafter, Night 1, Night 2, Night 3 and Night 4, respectively) for micro-variability study. The source brightness in magnitudes was converted into the flux in mJy units by using the zero points for the Johnson-Cousins UBVRI photometric system as given in Table A2 of [43], and to calculate the optical spectra the fluxes were interstellar-extinction corrected using the extinction magnitudes for various filters listed in the NED.3

3. Analysis and Results

To study the multiband optical variability properties of a sample of radio-loud sources, we have been frequently monitoring them for the past several years4. In general, the observations can be of the two types: single observation of the source for the purpose of long-term monitoring and microvariability observations, i.e., continuously tracking the source throughout the night. The left panel of Figure 1 presents the multi-band light curves spanning ~five months showing the mean magnitude of the nightly photometric observations of BL Lac in 2016. The figure shows that the source, displaying dramatic variability, brightens by ~1 mag (~2.5 times) in V and R band within a timescale of a month. The right panel of the figure shows the spectra in the optical bands (VRI) which can be approximated as power-law shapes with a mean spectral slope of 1.1. However, it can be seen that some of the spectra are markedly different from the power-law shapes, e.g., on the 164th day, revealing significant spectral breaks. We examined all the observations carefully to make sure that such spectral breaks were not resulted due to extrinsic factors including observing conditions.

In addition to the long-term flux points, we performed microvariability observed in four nights as mentioned earlier. The details of the multiband micorvariability analysis of the blazar BL Lac and the results following are discussed below.

3.1. Multi-Band Micro-Variability

The intranight variability in the photometric bands VRI exhibited by the source during the four nights are presented in Figure 2. The observation date, photometric filter, the mean magnitude and the observation duration are listed in the 1st, 2nd, 3rd and 4th column of

Table 1. Summary of the microvariability program in the year 2016.

Date	Band	Mean Mag.	Duration (h)	VA (mag)	F${}_{\mathit{var}}$ (%)
2016-08-07	R	13.48 ± 0.002	4.3	0.08	1.47 ± 0.03
2016-08-12	V	13.986 ± 0.003	3.8	0.13	3.70 ± 0.05
2016-08-12	R	13.359 ± 0.002	3.8	0.11	3.37 ± 0.05
2016-09-30	V	14.099 ± 0.018	6.2	0.08	1.61 ± 0.03
2016-09-30	I	12.670 ± 0.003	6.2	0.07	1.35 ± 0.03
2016-10-05	V	13.5640 ± 0.003	2.9	0.19	4.75 ± 0.07
2016-10-05	R	12.970 ± 0.002	2.9	0.17	4.30 ± 0.07
2016-10-05	I	12.193 ± 0.002	2.9	0.14	3.79 ± 0.07

, respectively. To quantify the observed variability, we estimated variability amplitude (VA) indicating peak-to-peak oscillation, and fractional variability (FV) representing mean variability. The amplitude of the peak-to-peak variations was estimated by using the relation given in [10],

$$\mathrm{VA}=\sqrt{{(A_{max}-A_{min})}^2-2{\sigma }^2},$$

(1)

where $A_{max}$, $A_{min}$, and $\sigma$ are the maximum, minimum, and the mean of the magnitude errors in the light curves, respectively. Similarly, for a mean flux of $〈F〉$ with $S^2$ variance and $〈{\sigma }_{err}^2〉$ mean squared uncertainties, the fractional variance, as in [44], is given as

$$F_{var}=\sqrt{\frac{S^2-〈{\sigma }_{err}^2〉}{{〈F〉}^2}}.$$

(2)

The error in the fractional variability can be expressed as

$${\sigma }_{F_{var}}=\sqrt{F_{var}^2+\sqrt{\frac{2}{N}\frac{{〈{\sigma }_{err}^2〉}^2}{{〈F〉}^4}+\frac{4}{N}\frac{〈{\sigma }_{err}^2〉}{{〈F〉}^2}{F}_{var}^2}}-F_{var}$$

(3)

(see [45]). The multiband VA and FV for all the observations are listed in the 5th and 6th column of Table 1, respectively. These quantities show that the source displays multiband fast variability by changing the flux by as large as ~0.2 R mag (or ~20%) within a timescale of ~3 h. In addition, we note that for the given night, both VA and $F_{\mathrm{var}}$ are consistently larger for the higher energy optical band implying larger variability for the higher frequency emission.

Color Variability: Bluer-When-Brighter

The relative source flux changes in different bands give rise to the intranight color evolution. Here color refers to the quasi-simultaneous (simultaneous within 1 min) difference in the magnitudes in V and I (or R) filters. Figure 3 shows that, along with flux variability, the source also exhibits considerable color variability (~0.05–0.20 mag) within a few hours. Figure 4 shows colors plotted against the corresponding V ( R for Night 3) band magnitudes along with the linear fits. To quantify the relation between the flux in magnitudes and color variations, we calculated the linear Pearson’s correlation coefficient between them. The coefficients for Night 2, Night 3 and Night 4 are 0.35, 0.50 and 0.70, respectively, at the confidence level of more than 99.99%. This indicates that the relation can be taken as bluer-when-bright trend which appears less, moderately and strongly pronounced during Night 2, Night 3 and Night 4, respectively. Similar observations were made by [46] in their study of the source on intranight timescales.

3.2. Multi-Band Cross-Correlation

An investigation into the relation among emission in different energy bands can offer an important insight regarding the structure of the blazar emission regions, the dominant radiative processes involved and distribution of the emitting particles. The discrete correlation function (DCF) described in [47] is one of the most extensively used methods to investigate the cross-correlation between two time series with uneven spacing (see [19]). The unbinned DCF can be estimated as,

$$UDC{F}_{ij}=\frac{\left(x_i-\overline{x}\right)\left(y_j-\overline{y}\right)}{\sqrt{\left({\sigma }_x^2-e_x^2\right)\left({\sigma }_y^2-e_y^2\right)}}$$

(4)

where $\overline{x}$ and $\overline{y}$, ${\sigma }^2$, and $e^2$ correspond to the mean, variance and the uncertainties in the magnitudes of the two light curves, respectively. These discrete pairs can be binned of bin width comparable to the sampling widths of the light curves. Then the average DCF including the M pairs within the given bin width is written as,

$$DCF\left(\tau \right)=\frac{1}{M}UDC{F}_{ij}$$

(5)

The sampling distributions of DCF, however, are highly skewed making the estimation of DCF uncertainties by the sample variances unreliable. The z-transformed DCFs (ZDCFs) follow closely normal distribution with known mean and standard deviation; this makes the estimation of ZDCF uncertainties more robust (For details see [48]).

We calculated the ZDCF between V and R band observations during Night 2 and between V and I band observation during Night 3 and Night 4, and estimated the likelihood of the ZDCF peaks and the associated uncertainties using the publicly available software5. For Night 2 and Night 4, the emission in the two bands were found to be highly correlated at the zero lag. However, for Night 3, a soft lag (I band emission lagging behind V band emission) of ${13}_{-5}^{+6}$ min (with correlation coefficient = 0.85 and likelihood = 0.74) was witnessed as shown in Figure 5.

Search for Characteristic Timescales

To search for the possible characteristic timescales present in the multiband light curves of the blazar BL Lac, we employed two widely used methods of time series: discrete auto-correlation function (ACF) and structure function (SF). Similar studies of blazars involving Fermi/LAT light curves were presented by [49]; and [38] used auto-correlation function (ACF) and SF to find the possible periodicities in the longterm radio light curve of the source. An expression for the discrete auto-correlation function (ACF) can be obtained by setting $x=y$ in the expressions given by Equations (4) and (5). The ZACFs, similar to ZDCF as discussed earlier, for all eight intranight light curves are presented in the eight panels of Figure 6. The figure shows presence of characteristic timescales in the R band light curve on Night 1 (Panel a) as represented by the well-resolved secondary peaks of at at $0.8\pm 0.02$ h with ZACF coefficient $0.79\pm 0.03$ and likelihood of 0.63.

To further investigate the presence of the characteristic timescales, we also analyzed the multiband observations employing first order structure function [50]. SFs are also widely used in the astronomical time series analysis; and for a signal $x\left(t\right)$ and a signal offset $\tau$, the first order structure function is given by

$$SF\left(\tau \right)=〈{[x\left(t\right)-x\left(t+\tau \right)]}^2〉$$

(6)

General interpretations of the SF features are discussed in [51]. The SF can be treated as equivalent to the power spectral density (PDS) of a light curve calculated in the time domain instead of frequency domain. One of the main advantages of the SF over PSD is that SF method is less sensitive to irregular sampling in the time series and therefore it is relatively free of artifacts similar to windowing and aliasing in frequency domain (e.g., [52,53]). This makes SF one of the favored tools in time series analysis widely applied in the search for periodicity in the AGN light curves (see [54] and references therein). Figure 7 shows the structure function plots for all the microvariability observations of the blazar BL Lac. Although characteristic variability timescales are not obvious in the most of the plots, the SF for R band observations on Night 1 (Figure 7a) shows clear break around ~0.8 h, confirming the presence of such a timescale detected by ACF. The statistical reliability of structure functions have been questioned (for a critical review of SF see [55], however in our case the fact that it matches the result from ZACF gives us confidence that the result is not merely statistical.

4. Discussion and Conclusions

It is seen that the blazar displays multi-band flux variability in the hour-like timescales. Using the causality argument, the minimum variability timescale ${\tau }_{var}$, can be used to estimate the upper limit of the minimum sizes of the emission region (R) as

$$R=\frac{\delta }{\left(1+z\right)}c{t}_{var}$$

(7)

where $\delta$, Doppler factor, is given $\delta ={\left(\mathsf{\Gamma }\left(1-\beta cos\theta \right)\right)}^{-1}$ and , for the velocity $\beta =v/c$, the bulk Lorentz factor can be written as $\mathsf{\Gamma }=1/\sqrt{1-{\beta }^2}$. In blazar jets, the emission regions move with high speeds along the path making an angle with the line of sight. Using a typical variability timescale $t_{var}=1$ h and source red-shift $z=0.0685$ along with a moderate value of Doppler factor $\delta =10$, the Equation (7) can be used to constrain the upper limit of the emission regions size to ~$1\times {10}^{15}$ cm, a size comparable to the solar system. Such small scale fluctuations leading to the rapid intranight multiband variability are most likely intrinsic to the jet processes and often hard to be explained as mere result of extrinsic processes such as geometrical and projection effects often involving variations in Doppler factor of the emission regions. In general, blazar variability should originate within a small, possibly independent sub-volumes of blazar jets, which could be identified with isolated turbulent cells, magnetic reconnection sites, small-scale shocks induced by turbulence within the main jet body (see in this context [8,18,20,56,57]). In particular, microvariability displayed by the source might also arise due to magnetohydrodynamical (MHD) instabilities inducing localized stochastic particle injections on similar time scales within the turbulent jets of the blazar. In addition, there is observational evidence as well as evidence from MHD simulations for the existence of standing shocks in the relativistic jets of radio galaxies, quasars, and blazars (see [58] and references therein). In one of the likely scenarios, passage of jet (turbulent) plasma material through the standing shock(s) and/or the after-effects may provoke a generation of microvariable emission in the jet of BL Lac.

Then, the combined effects of the cooling (synchrotron and adiabatic expansion) on the accelerated particles and the light crossing timescales contribute in shaping the observed fast variability. The frequency dependent variability amplitude i.e., larger variability amplitude at higher frequency,—also reported in previous works on intranight variability (e.g., [31,32,34])—could also be the result of synchrotron emission; assuming constant magnetic field and isotropic distribution of pitch angles, the synchrotron emission (by particles with Lorentz factor $\gamma$) together with $-d\gamma /dt\propto {\gamma }^2$ and $t_{syn}\propto 1/\gamma$ could result dissipation larger amount of energy in shorter timescales. For such calculations, BL Lac being a LSP source, it is also considered that most of the synchrotron radiation is emitted around optical frequency (as characteristic frequency).

In addition to the rapid flux variability, the observed significant color (spectral) evolution on intranight timescale further help us to delve into the emission regions. In particular, bluer-when-brighter (alternatively, harder-when-brighter) trend has been frequently observed in the source during intranight observations (see [34,37]). Such an achromatic behavior, both bluer-when-brighter and redder-when-brighter, have been frequently observed in the blazars (for further discussion see [2,19]). Similarly, the optical spectra are found to be of the power-law type consistent with the scenario in which the synchrotron emission contribute by power-law particle injections. Another interesting observation was that the optical spectra revealed a power-spectral break. Such breaks could arise due to the break in the injected particle distribution. Similar spectral curvatures in the source were also detected during intraday flux changes in Swift X-ray observations [26] and in hard X-ray (3–79 keV) NuSTAR observation (see [23]). Such spectral curvature (or the spectral breaks) often indicate the transition region between the low- and high-energy components of the spectral energy distribution (SED) (see [59]).

The study of the correlation between multiband emission on intranight timescales shows that, in most cases, both V and R and V and I band emission appear well correlated with each other with no evident lag. This suggests that microvariability primarily originates from a single emission zone with the similar population of the particles. However, during one of the nights, we found that the I band variations lag behind V band variations, also noted in the source by previous authors (e.g., [37]) on similar timescale. Such lags between multifrequency emission could be resulted due to frequency dependent opacity in which the higher emission escaped the emission region first followed by the lower energy emission. Alternatively, the observed lag could be resulted due to the separation of the two (high and low energy) emission regions (see [60]). Furthermore, both discrete auto-correlation function and structure function analyses consistently suggest presence of a characteristic timescales of ~0.86 h in the R band observation of Night 1. However we do not find similar detection in the other observations. This implies the transient nature of processes resulting in the characteristic timescale. Such short timescales characteristic timescales might have resulted due to localized magnetohydrodynamical instabilities either in the accretion disk (e.g., [61]) or in the relativistic jets (e.g., [62]). In the disc, such a small timescale processes might occur in the innermost regions whereas in the jets the rest frame timescale might be longer by the Doppler factor of the emission region and thereby leaving door open for a wider range of likely scenarios depending upon the combined effect of the changes in the orientation angle and the bulk speed of the emission region.