Multi-wavelength variability study of the classical BL Lac object PKS 0735+178 on timescales ranging from decades to minutes

TL;DR

This study analyzes the multi-wavelength variability of the BL Lac object PKS 0735+178 over 23 years, revealing a consistent red noise behavior in optical variability and flatter gamma-ray spectra, advancing understanding of blazar stochastic processes.

Contribution

First to construct an optical power spectrum of a blazar over timescales from minutes to decades using combined multi-instrument data.

Findings

Optical variability follows a single power-law red noise across all timescales.

Gamma-ray power spectral slope (~1) is flatter than optical and radio (~2).

No low-frequency flattening or high-frequency cut-offs detected in the power spectra.

Abstract

We present the results of our power spectral analysis for the BL Lac object PKS 0735+178 utilizing the Fermi-LAT survey at high-energy $\gamma$-rays, several ground-based optical telescopes, and single-dish radio telescopes operating at GHz frequencies. The novelty of our approach is that, by combining long-term and densely sampled intra-night light curves in the optical regime, we were able to construct for the first time the optical power spectrum of the blazar for a time domain extending from 23 years down to minutes. Our analysis reveals that: (i) the optical variability is consistent with a pure red noise, for which the power spectral density can well be approximated by a single power-law throughout the entire time domain probed; (ii) the slope of power spectral density at high-energy $\gamma$-rays ($\sim 1$), is significantly flatter than that found at radio and optical frequencies ($\sim 2$) within the corresponding time variability range; (iii) for the derived power spectra we did not detect any low-frequency flattening, nor do we see any evidence for cut-offs at the highest frequencies down to the noise floor levels due to measurement uncertainties. We interpret our findings in terms of a model where the blazar variability is generated by the underlying single stochastic process (at radio and optical frequencies), or a linear superposition of such processes (in the $\gamma$-ray regime). Along with the detailed PSD analysis, we also present the results of our extended (1998-2015) intra-night optical monitoring program and newly acquired optical photo-polarimetric data for the source.

Tables

Table 1. Table 1 Observational summary of the INV data and data analysis

Date of obs.	Tel.	Dur.	$N_p$	$\Delta m_{t-s1}$	$\Delta m_{t-s2}$	$\Delta m_{s1-s2}$	${\sigma }_{s1-s2}$	SDs1-s2	$\psi$	$F-test$:value (status)†			INV†	Ref.
		(hr)		(mag)	(mag)	(mag)	(1$e$-2 mag)	(1$e$-2 mag)	(1$e$-2 mag)	(t-s1)	(t-s2)	(s1-s2)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)	(15)
1998 Dec 26	ST	7.8	49	$-$0.531	$-$0.407	0.124	0.4	0.6	2.6	1.38(N)	1.05(N)	0.75(N)	N	(a)
1999 Dec 30	ST	7.4	64	$-$0.988	$-$1.055	$-$0.067	0.4	0.5	2.1	0.49(N)	0.68(N)	0.62(N)	N	(a)
2000 Dec 24	ST	6.0	42	$-$1.676	$-$1.737	$-$0.061	0.4	0.5	2.2	1.26(N)	0.70(N)	0.70(N)	N	(a)
2001 Dec 24	ST	7.3	38	$-$1.115	$-$1.306	$-$0.191	0.3	0.4	1.8	1.46(N)	0.98(N)	0.52(N)	N	(a)
2003 Dec 20	HCT	6.0	38	$-$0.774	0.045	0.819	0.2	0.3	1.5	1.76(PV)	2.13(PV)	0.80(N)	PV	(b)
2004 Dec 10	ST	6.2	30	0.286	$-$0.227	$-$0.513	0.2	0.3	2.1	4.48(V)	3.54(V)	1.17(N)	V	(b)
2004 Dec 23	ST	5.9	13	0.196	$-$0.310	$-$0.506	0.2	0.3	1.5	5.59(V)	3.35(PV)	1.15(N)	V	(b)
2005 Jan 2	ST	4.9	22	0.108	$-$0.414	$-$0.522	0.2	0.3	1.2	0.85(N)	1.65(N)	0.81(N)	N	(b)
2005 Jan 5	ST	5.2	26	0.269	0.111	$-$0.158	0.1	0.2	1.4	4.45(V)	3.09(V)	1.08(N)	V	(b)
2005 Jan 9	ST	7.1	30	0.325	0.172	$-$0.153	0.1	0.2	1.5	3.50(V)	3.66(V)	0.90(N)	V	(b)
2005 Nov 9	ST	4.3	19	0.449	0.362	$-$0.086	0.1	0.3	1.2	3.00(PV)	2.32(PV)	1.74(N)	PV	(b)
2006 Nov 16	ST	5.0	21	0.580	1.550	0.970	0.2	0.4	1.4	1.72(N)	0.82(N)	2.49(N)	N	(b)
2006 Nov 29	ST	6.5	28	0.947	0.430	$-$0.517	0.2	0.3	1.9	1.31(N)	1.27(N)	1.00(N)	N	(b)
2006 Dec 17	ST	6.5	28	1.013	0.505	$-$0.507	0.2	0.3	2.2	1.87(N)	1.93(PV)	1.45(N)	N	(b)
2007 Jan 11	Ch	3.9	90	1.053	1.924	0.870	0.5	0.6	4.8	0.59(N)	0.58(N)	0.64(N)	N	(c)
2007 Dec 15	ST	7.1	29	0.364	0.201	$-$0.163	0.1	0.2	2.6	4.62(V)	4.99(V)	1.35(N)	V	(b)
2007 Dec 16	ST	7.1	29	0.281	$-$0.228	$-$0.509	0.2	0.2	1.7	2.16(N)	1.31(N)	0.46(N)	N	(b)
2008 Nov 22	ST	6.0	29	1.037	0.908	$-$0.129	0.2	0.2	1.4	0.64(N)	0.74(N)	0.53(N)	N	(b)
2009 Jan 20	ST	4.0	66	0.459	2.043	1.584	0.2	0.7	5.5	6.01(V)	7.47(V)	5.26(V)	N	(d)
2009 Dec 8	ST	6.9	31	0.455	0.326	$-$0.129	0.4	0.5	3.2	1.65(N)	1.23(N)	0.91(N)	N	(e)
2011 Jan 5	ST	6.8	32	1.284	0.953	$-$0.331	0.3	0.4	3.8	1.13(N)	1.16(N)	0.51(N)	N	(e)
2011 Nov 29	ST	6.1	29	1.113	0.613	$-$0.500	0.2	0.3	1.9	1.18(N)	0.65(N)	0.51(N)	N	(e)
2012 Dec 20	ST	6.9	115										N	(f)
2013 Jan 7	ST	5.6	22	0.535	0.383	$-$0.152	0.3	0.4	4.9	6.43(V)	7.00(V)	0.91(N)	V	(e)
2015 Feb 15	AOJU	4.1	30	0.630	1.143	0.513	0.8	1.1	9.1	1.99(PV)	1.86(PV)	0.81(N)	PV	(e)

Table 2. Table 2 The PSD analysis for the light curves of PKS 0735+178

Light curve	Range	$N_{obs}$	$\Delta T_{\mathrm{obs}}$	$\Delta T_{\mathrm{intep}}$	$T_{\mathrm{obs}}$	$\log \left(P_{\mathrm{stat}}\right)$	$\log \left({\nu }_k\right)$ range	$\beta \pm err$
			(day)	(day)		((rms/mean)2 day)	(day${}_{}{}^{-1})$	bin
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)
Fermi-LAT (LT)	0.1-200 (GeV)	168	15	0.5	6.9 (yr)	$+$0.50	$-$3.4 to $-$1.5	1.1$\pm$0.2
OVRO (LT)	15 (GHz)	263	3	0.5	5.1 (yr)	$-$2.24	$-$3.3 to $-$1.0	1.8$\pm$0.1
UMRAO (LT)	14.5 (GHz)	359	15	0.5	33 (yr)	$-$1.00	$-$4.1 to $-$1.6	2.1$\pm$0.1
UMRAO (LT)	8.0 (GHz)	299	15	0.5	31 (yr)	$-$1.06	$-$4.1 to $-$1.6	2.2$\pm$0.1
UMRAO (LT)	4.8 (GHz)	220	15	0.5	30 (yr)	$-$0.99	$-$4.0 to $-$1.6	2.6$\pm$0.1
Optical (LT)	R-band	677	1	0.5	22.6 (yr)	$-$2.27	$-$3.8 to $-$1.1	2.0$\pm$0.1
2004 Dec 10 (IN)	R-band	30	0.008	0.0004	6.2 (hr)	$-$7.23	$+$0.6 to $+$1.8	2.8$\pm$0.7
2004 Dec 23 (IN)	R-band	13	0.017	0.0004	5.9 (hr)	$-$6.83	$+$0.6 to $+$1.3	4.1$\pm$0.3
2005 Jan 5 (IN)	R-band	26	0.008	0.0004	5.2 (hr)	$-$7.23	$+$0.8 to $+$1.8	2.2$\pm$0.3
2005 Jan 9 (IN)	R-band	30	0.008	0.0004	7.1 (hr)	$-$7.21	$+$0.5 to $+$1.7	1.5$\pm$0.2
2007 Dec 15 (IN)	R-band	29	0.008	0.0004	7.0 (hr)	$-$7.10	$+$0.5 to $+$1.7	2.3$\pm$0.5
2013 Jan 7 (IN)	R-band	22	0.008	0.0004	5.6 (hr)	$-$6.73	$+$0.5 to $+$1.6	3.3$\pm$0.4

Equations

$$F_{\nu }^{\alpha }=\frac{V_{\mathrm{obs}}}{V_{\exp }}=\frac{V_{t-s}}{⟨{\eta }^2{\sigma }_{t-s}^2⟩},$$

$$\psi =\sqrt{(D_{\max }-D_{\min }{)}^2-2{\sigma }^2},$$

$$DC=100\%\frac{{\sum }_{j=1}^n{N}_j(1/\Delta t_j)}{{\sum }_{j=1}^n(1/\Delta t_j)},$$

$$P({\nu }_k)=\frac{2T}{{\mu }^2{N}^2}|F({\nu }_k){|}^2,$$

$$P_{stat}=\frac{2T}{{\mu }^{2}N}{\sigma }_{\mathrm{stat}}^2.$$

$$\log [P({\nu }_k)]=\log [(P_{true}({\nu }_k)]+\log [\frac{{\chi }^2}{2}].$$

$$P(\nu_{k})=\frac{2\,T}{\mu^{2}N^{2}}\,\left\{\Bigg{[}\sum_{i=1}^{N}f(t_{i})\,\cos(2\pi\nu_{k}t_{i})\Bigg{]}^{2}+\Bigg{[}\sum_{i=1}^{N}f(t_{i})\,\sin(2\pi\nu_{k}t_{i})\Bigg{]}^{2}\right\}\,,$$

$$P_{f}(2\pi\nu_{k})=\frac{1}{2}\left\{\frac{\Bigg{\{}\sum_{i=1}^{N}f(t_{i})\cos[2\pi\nu_{k}(t_{i}-\omega)]\Bigg{\}}^{2}}{\sum_{i=1}^{N}f(t_{i})\cos^{2}[2\pi\nu_{k}(t_{i}-\omega)]}+\frac{\Bigg{\{}\sum_{i=1}^{N}f(t_{i})\sin[2\pi\nu_{k}(t_{i}-\omega)]\Bigg{\}}^{2}}{\sum_{i=1}^{N}f(t_{i})\sin^{2}[2\pi\nu_{k}(t_{i}-\omega)]}\right\}\,,$$

$$\omega=\arctan\Bigg{\{}\frac{1}{2(2\pi\nu_{k})}\,\frac{\sum_{i=1}^{N}\sin[2(2\pi\nu_{k})t_{i}]}{\sum_{i=1}^{N}\cos[2(2\pi\nu_{k})t_{i}]}\Bigg{\}}\,,$$

$$r(\tau )=\frac{E\{[f(t)-\mu ]\times [f(t+\tau )-\mu ]\}}{{\sigma }_f^2},$$

$$P({\nu }_k)=2\sum _{i=1}^{N}r(\tau )\times \cos (2\pi {\nu }_k\tau ).$$

Full text

22affiliationtext: Astronomical Institute of St. Petersburg State University, Petrodvorets 198504, Russia33affiliationtext: Centre for Excellence in Basic Sciences (CEBS), University of Mumbai campus (Kalina), Mumbai 400098, India44affiliationtext: Department of Physics, The College of New Jersey, 2000 Pennington Rd., Ewing, NJ 08628-0718, USA55affiliationtext: Aryabhatta Research Institute of Observational Sciences (ARIES), Manora Peak, Nainital 263002, India66affiliationtext: Instituto de Astrofísica de Andalucía (CSIC), Apartado 3004, E–18080 Granada, Spain

1. Introduction

The blazar class of active galactic nuclei (AGN) includes BL Lacertae objects (BL Lacs) and high polarization quasars (HPQs, which are a subset of flat spectrum radio quasars, FSRQs, characterized by high fractional optical polarization degree ${\rm PD_{opt}}>3\%$; see Urry & Padovani, 1995). For these, the total radiative energy output is dominated by the broad-band, non-thermal emission produced in relativistic jets. These jets are launched by supermassive black holes surrounded by magnetized accretion disks, that are located at the centers of (usually) massive elliptical galaxies, and when pointing close to our line of sight the emission is substantially Doppler boosted (Begelman et al., 1984; de Young, 2002; Meier, 2012). The radio-to-optical/X-ray segment of the blazar emission continuum is due to the synchrotron radiation of ultrarelativistic leptons (hereafter “electrons”, for simplicity), while the high-frequency X-ray-to-$\gamma$-ray segment is most widely believed to be due to the inverse-Comptonization of various circumnuclear photon fields (produced both internally and externally to the outflow) by the jet electrons.

Blazars display a strong variability from radio to $\gamma$-rays on timescales ranging from decades down to hours, or even minutes. The observed flux changes are often classified broadly into the three major types, namely a ‘long-term variability’ (LTV, with the corresponding timescales of decades-to-months), a ‘short-term variability’ (STV; weeks-to-days), and an intra-night/day variability (INV/IDV; timescales less than a day; see, e.g., Wagner & Witzel, 1995; Ulrich et al., 1997; Falomo et al., 2014). Rapid – hour- and minute-long – blazar flares are especially pronounced at higher energies, particularly in the $\gamma$-ray regime, with intensity changes of up to even a few magnitudes (Aharonian et al., 2007; Aleksić et al., 2011; Ackermann et al., 2016; Foschini et al., 2011; Saito et al., 2013; Rani et al., 2013). The origin of such dramatic behavior, and in particular its relation to the similarly rapid but smaller-amplitude ‘microvariability’ observed at lower wavelengths, (including INV in the optical range), is still being widely debated.

Several competing scenarios have been proposed to explain the rapid variability of blazar sources. Some include extrinsic causes, such as gravitational lensing (where the observed flux changes are attributed to a lensing of light rays by a foreground compact object; e.g., Schneider & Weiss, 1987; Gopal-Krishna & Subramanian, 1991; Webb et al., 2000) or, in the radio domain, interstellar scintillation (see Melrose, 1994, for a review). Others involve an intrinsic origin, including some that are purely geometrical in nature, such as ‘the lighthouse effect’, where the precession of a jet results in the differential forward beaming of the emission (Camenzind & Krockenberger, 1992; Gopal-Krishna & Wiita, 1992). A set of promising intrinsic origin models, which may be favored because of polarimetric measurements revealing the accompanying changes in the fractional polarization and polarization angle, involves various plasma instabilities leading to the formation of shocks and turbulence in the outflow, which then heat and accelerate the jet particles (e.g., Narayan & Piran, 2012; Subramanian et al., 2012; Marscher, 2014; Saito et al., 2015, and references therein). Other intrinsic models invoke the annihilation of magnetic field lines of opposite polarity, transferring the energy from the field to particles at the reconnection sites and during the subsequent evolution of the reconnected magnetic field (e.g., Giannios et al., 2009; Sironi et al., 2015; Guo et al., 2016).

At optical frequencies, the LTV of non-blazar AGN, including radio-quiet quasars (RQQs), radio-intermediate quasars (RIQs), lobe-dominated radio-loud quasars (LDQs), and low-polarization flat spectrum radio quasars (LPQs), is characterized by much smoother, yet comparable in amplitude, intensity changes to those observed in blazar sources (e.g., Stalin et al., 2005; Mushotzky et al., 2011; Edelson et al., 2013). This finding is surprising, keeping in mind that the bulk of the optical emission of these non-blazar type active galaxies is understood to originate in accretion disks, rather than relativistic jets (e.g., Wiita, 2006). On the other hand, a clear dichotomy exists on intra-night timescales, where the blazar class shows a considerably higher optical INV in both amplitude, $\psi$ (see equation 2 below), and duty cycle, DC (equation 3), than non-blazar objects (e.g., Gopal-Krishna et al., 2003; Stalin et al., 2005; Gupta & Joshi, 2005; Villforth et al., 2009; Goyal et al., 2010). More recently, based on a systematic study using 262 intra-night light curves, each monitored for a duration of $\geq$ 4 hr, it was shown that the optical INV duty cycle is DC $\sim 40\%$ for the blazar class, while it is $\sim 5\%$, $\sim 11\%$, $\sim 3\%$, and $\sim 10\%$ for the RQQs, RIQs, LDQs, and LPQs, respectively, at least whenever clear INV amplitudes of $\psi\geq 3\%$ are considered (Goyal et al., 2013a).

Even though the physical processes giving rise to the flaring emission of blazars remain debatable, considerable progress has been made in characterizing the statistical properties of blazar variability at different wavelengths, and in different time domains. In particular, it has been demonstrated repeatedly that the power spectral density (PSD) of blazar light curves is, in general, of a power-law form (Simonetti et al., 1985; Kataoka et al., 2001; Papadakis et al., 2003; Brinkmann et al., 2003; Aharonian et al., 2007; Ciprini et al., 2007; Chatterjee et al., 2008, 2012; Abdo et al., 2010a; Carini et al., 2011; Kastendieck et al., 2011; Edelson et al., 2013; Nakagawa & Mori, 2013; Sobolewska et al., 2014; Revalski et al., 2014; Park & Trippe, 2014; Max-Moerbeck et al., 2014; Aleksić et al., 2015; Isobe et al., 2015; Kasliwal et al., 2015). A physical process with such a variability power spectrum, denoted hereafter as $P(\nu_{k})=A\,\nu_{k}^{-\beta}$, where $\nu_{k}$ is the temporal frequency (corresponding to the timescale $1/\nu_{k}$), $A$ is the normalization constant, and $\beta$ is the spectral slope, is called white noise when $\beta=0$, flicker (pink) noise when $\beta=1$, and Brownian (red) noise when $\beta=2$ (Press, 1978). The PSD integrated over some variability frequency range is then a measure of the variance of the underlying signal in the time series within the corresponding range of variability timescales. Breaks in the slope or in the normalization of a PSD may appear, signaling characteristic/critical variability timescales in the system. In the case of blazars, various segments of radio, optical, X-ray, and $\gamma$-ray power spectra within the variability time domains from years to days (and in some instances, even sub-hour timescales), are characterized by spectral slopes $1\leq\beta<3$, meaning that the variability amplitude increases with increasing variability timescale. Rarely, however, have blazar PSDs been analyzed in a systematic way at different wavelengths across the electromagnetic spectrum and over a truly broad range of temporal frequencies. It is important to note that colored noise-type power spectra are expected to flatten on longer variability timescales (to preserve the total finite variance), and to cut-off at frequencies corresponding to the shortest variability timescale in a system. The detection of such cutoffs in blazar periodograms would be of a primary importance for constraining the physics of blazar jets; however, such detections may be hampered by the finite duration of available monitoring blazar data, on the one hand, and statistical fluctuations resulting from the measurement errors, on the other hand.

In this work, we present our analysis and interpretation of the multi-wavelength (radio, optical, and high-energy $\gamma$-ray), and particularly long-span (decades to minutes) variability power spectrum of the BL Lac object PKS 0735$+$178. This source is singled out from the blazar class by its persistently weak intra-night optical variability: using 17 nights data spanning over 11 years of optical monitoring of this blazar, Goyal et al. (2009) estimated the INV DC as $\sim 0\%$ for $\psi>3\%$. Note that while such a low duty cycle is not unusual for non-blazar AGN (see above), it is surprising for highly polarized blazars such as PKS 0735$+$178. In this study, we present the results of our extended intra-night monitoring programme of the source, now consisting a total of 25 nights and spanning over 18 yr of observations (1998–2015).

PKS 0735$+$178 (J2000.0 R.A. $\rm=07^{h}38^{m}07\fs 39$, Dec. $\rm=+17\arcdeg 42\arcmin 19\farcs 99$; redshift $z=0.45\pm 0.06$; Nilsson et al., 2012) is an otherwise typical example of a low-frequency-peaked BL Lac object detected in the GeV photon energy range (Abdo et al., 2010b). It is highly polarized in the optical band (${\rm PD_{opt}}>3\%$; Gabuzda et al., 1996; Wills et al., 2011), and exhibits a flat-spectrum radio core with a superluminal pc-scale radio jet, both characteristic of blazars in general (Wills et al., 1992; Lister et al., 2016). Despite its pronounced optical and radio variability on yearly timescales (Webb et al., 1988; Ciprini et al., 2007), the source is relatively quiet in X-rays (Madejski & Schwartz, 1988), and the existing rather sparse X-ray monitoring data preclude a meaningful power spectral analysis. Within the high-energy $\gamma$-ray regime covered by Fermi-LAT, the blazar is detected at a high significance level only with weekly or longer temporal binning (Nolan et al., 2003; Abdo et al., 2010a). Upper limits for the PKS 0735$+$178 emission in the very high-energy $\gamma$-ray domain (photon energies $>100$ GeV) have been recently provided by the VERITAS Collaboration (Archambault et al., 2016).

As revealed by very long baseline interferometry (VLBI) imaging, the radio jet of the blazar underwent several dramatic structural changes between 1981 and 2001 on the scale of 10 pc, altering between a ‘staircase jet’ with a highly bent trajectory (circa 1981 and 1995), and a straight jet with a linear trajectory (1985 and 2001); the last morphological transition coincided with the optical flaring but was accompanied by only a mild increase in the radio intensity (Britzen et al., 2010). By inspecting milli-arcsec scale resolution images from the MOJAVE database111http://www.physics.purdue.edu/astro/MOJAVE/sourcepages/0735+178.shtml and the Boston University database222https://www.bu.edu/blazars/VLBA_GLAST/0735.html, we confirmed that the jet trajectory has remained linear from 2001 until 2015. Faraday rotation gradients and circular polarization have been detected in the radio core of PKS 0735$+$178 (Gabuzda et al., 2008). The jet magnetic field structure revealed by the radio polarization maps is rather complex (though again, not unusual for a low-frequency-peaked BL Lac), consisting of magnetic field lines mostly perpendicular to the jet axis near the core region, but parallel to the jet further along the outflow; this change suggests either a helical configuration with changing pitch angle or strong interactions with the ambient medium leading to velocity shears and jet bends (Agudo et al., 2006).

In Section 2 we describe our data and its reduction. Our analysis and results are given in Section 3 and Section 4 provides a discussion of the results and our main conclusions

2. Data acquisition

2.1. Optical: intra-night

The vast majority of the intra-night observations were carried out using the 104-cm Sampurnanand telescope (ST) located at the Aryabhatta Research Institute of observational sciencES (ARIES), Naini Tal, India. The ST has Ritchey-Chrétien (RC) optics with a f$/$13 beam (Sagar, 1999). The detector was a cryogenically cooled $2048\times 2048$ chip mounted at the Cassegrain focus. This chip has a readout noise of 5.3 e*-/pixel and a gain of 10 e-$/$Analog to Digital Unit (ADU) in slow readout mode. Each pixel has a dimension of 24 $\mu$m2* which corresponds to 0.37 arcsec2 on the sky, covering a total field of $13^{\prime}\times 13^{\prime}$. Our observations were carried out in $2\times 2$ binned mode to improve the signal-to-noise ratio.

We similarly used the 201-cm Himalayan Chandra Telescope (HCT) at the Indian Astronomical Observatory (IAO), located in Hanle, India. This telescope is also of the RC design but has a f$/$9 beam at the Cassegrain focus333http://www.iiap.res.in/$\sim$iao. The detector was a cryogenically cooled $2048\times 4096$ chip, of which the central $2048\times 2048$ pixels were used. The pixel size is 15 $\mu$m2, so that the image scale of 0.29 arcsec$/$pixel covers an area of $10^{\prime}\times 10{{}^{\prime}}$ on the sky. The readout noise of CCD is 4.87 e*-/pixel and the gain is 1.22 e-*$/$ADU. The CCD was used in the unbinned mode.

Lastly, we also employed the 50-cm Cassegrain telescope in the Astronomical Observatory of the Jagiellonian University (AOJU), located in Kraków, Poland. This telescope is also of the RC design with f$/$6.7 beam at the Cassegrain focus. The detector was a thermoelectric cooled $1024\times 1024$ chip, corresponding to an image scale of 0.70 arcsec$/$pixel, covering a total of $12^{\prime}\times 12{{}^{\prime}}$ on the sky. The CCD was used in $2\times 2$ binned mode.

The seeing mostly ranged between $\sim 1^{\prime\prime}.5$ to $\sim 3^{\prime\prime}$, as determined using three sufficiently bright stars on each CCD frame. All the observations were made using the R filter, as the CCD responses are maximized in this band. The field positioning was adjusted in order to have within each CCD frame at least two, and usually three, comparison stars. For all the telescopes, bias frames were taken intermittently, and twilight sky flats were also obtained.

The pre-processing of the images (bias subtraction, flat-fielding and cosmic-ray removal) was done by applying the standard procedures in the Image Reduction and Analysis Facility (IRAF)444http://iraf.noao.edu/ software package. The instrumental magnitudes of the target AGN and the stars (all point-like) in the image frames were determined by the aperture photometry using APPHOT. The magnitude of the target AGN was measured relative to a few apparently steady comparison stars present on the same CCD frame; relative star-star magnitudes were also recorded. In this way, the Differential Light Curves (DLCs) for the AGN and the comparison stars were derived. Out of the resulting three star-star DLCs, we selected the steadiest star-star DLC (based on the lowest variance) for testing the INV of the blazar monitored on a given night. These chosen stars are hereafter named ‘s1’ and ‘s2’, and the corresponding target-star and star-star DLCs are denoted as ‘t-s1’, ‘t-s2’, ‘s1-s2’, respectively. Basic information about the comparison stars (apparent magnitudes, optical colors, positions on the sky) is given in Goyal et al. (2009).

The comparison stars we used are typically within about one magnitude of the target AGN; note that avoiding large differences in brightness is of crucial importance for minimizing the possibility of a spurious INV detection (e.g., Cellone et al., 2007). In this context, spurious variability on account of different second-order extinction coefficients for the AGN and their comparison stars can also be problematic if the target AGN and the comparison stars have very different optical colors. However, as shown by Carini et al. (1992) and Stalin et al. (2004), for color differences of up to 1-2 mag, the differential extinction of photons traveling through varying airmasses do not influence significantly the derived INV parameters of BL Lac objects, given their typical flux measurement uncertainties ($\simeq 0.1-0.2\%$, including the observations analyzed here). For each night, an optimum aperture radius for the photometry was chosen by identifying the minimum dispersion in the star-star DLC, starting from the median seeing (i.e., full width at half maxima) value on that night to four times that value. For a given night, we selected the aperture size which yielded the minimum dispersion in the steadiest star-star DLC. This also set the threshold for the INV detection on that night (see Goyal et al., 2013b). Typically, the selected aperture radius was $\sim 4^{\prime\prime}$ and the effective seeing was $\sim 2^{\prime\prime}$.

Our entire intra-night data are summarized in Table 1, and the newly acquired DLCs are shown in Figure 1.

2.2. Optical: long/short-term

The nightly averaged R-band photo-polarimetric data for the period 2005 October 22 to 2015 September 27 were obtained using (i) the 70-cm AZT-8 telescope of the Crimean Astrophysical Observatory (Nauchnij, Ukraine), (ii) the 40-cm LX-200 telescope of the St. Petersburg State University (St. Petersburg, Russia), (iii) the 1.8-m Perkins telescope of the Lowell Observatory (Flagstaff AZ, USA), (iv) the 1.54-m Kuiper and 2.3-m Bok telescopes of the Steward Observatory (Mt. Bigelow AZ and Kitt Peak AZ, USA), and (v) the 2.2-m telescope of the Calar Alto Observatory (Calar Alto, Spain) within the MAPCAT program555http://www.iaa.es/$\sim$iagudo/research/MAPCAT/MAPCAT.html (see Larionov et al. 2008 for the description of the program and analysis). The photometric data were supplemented with the nightly averaged R-band optical data from 1993 until 2005, obtained from Ciprini et al. (2007). These nightly averaged data have quoted photometric uncertainty of the order of $\sim 5-10\%$, arising mainly from large calibration errors in the estimated magnitudes of the stars in the field (Ciprini et al., 2007). It is a standard procedure to use only one or two stars to scale the magnitude of a blazar to the standard system. In case of intra-night data, however, as a very high photometric accuracy is required to detect variability down to $1\%$ amplitudes, we standardized our comparison stars (Table 2 of Goyal et al., 2009) using all the standard stars in the field of the target, as listed in Ciprini et al. (2007); this resulted in the accuracy of $\leq 0.2-0.5\%$. Finally, for a given R-band magnitude $M_{\rm R}$ the R-band flux (in Jy) was derived as $3064\times 10^{-0.4\times M_{\rm R}}$, where 3064 Jy is the zero point magnitude flux of the photometric system (Glass, 1999); the errors in R-band fluxes were derived using standard error propagation (Bevington & Robinson, 2003).

2.3. High energy $\gamma$-rays: long-term

We have analyzed the Fermi-LAT data for the field containing PKS 0735+178 from 2008 August through 2015 September, and produced a source light curve between 0.1 and 200 GeV with an integration time of 15 days. We have performed the unbinned likelihood analysis using Fermi ScienceTools-v10r0p5 with p8r2_source_v6 source event selection and instrument response function, for the $20^{\circ}$ region centered at the blazar, following the Fermi tutorial666http://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/. The procedure starts with the selection of good data and time intervals (using the tasks ‘gtselect’ and ‘gtmktime’ with selection cuts evclass=128 evtype=3), followed by the creation of an exposure map in the region of interest (ROI) with $30^{\circ}$ radius for each time bin (tasks ‘gtltcube’, ‘gtexpmap’ while counting photons within zenith angle $<90^{\circ}$). We then computed the diffuse source response (task ‘gtdifrsp’), and finally modeled the data through a maximum-likelihood method (task ‘gtlike’). In this last step, we used a model that includes PKS 0735+178 and 170 other sources inside the ROI (according to the third Fermi Large Area Telescope source catalog, 3FGL; Ackermann et al., 2015). The model also takes into account the diffuse emission from our Galaxy and the extragalactic $\gamma$-ray background777gll_iem_v06.fits and iso_p8r2_source_v6_v06.txt(Acero et al., 2016). In the modelling, we followed the usual method and fixed the spectral indices and fluxes of all the point sources within the ROI, other than the target, at their 3FGL values. The $\gamma$-ray spectrum of PKS 0735+178 was modeled with a simple power law. We considered a measurement to be a successful detection for a test statistic TS $\geq$ 10, which corresponds to a signal-to-noise ratio $\geq 3\sigma$ (Abdo et al., 2009).

2.4. Radio: long-term

The radio data were obtained from the University of Michigan Radio Astronomy Observatory (UMRAO) 26m dish at 4.8, 8.0, and 14.5 GHz, and the 40-m Telescope at the Owens Valley Radio Observatory (OVRO) at 15 GHz. The UMRAO fluxes at 4.8, 8.0, and 14.5 GHz were typically sampled twice per month (Aller et al., 1985), from 1980 February 16 – 2010 April 16, from 1977 March 11 – 2008 April 5, and from 1977 August 19 – 2011 January 24, respectively, while the OVRO light curve at 15 GHz was sampled twice a week (Richards et al., 2011), during the period from 2009 march 17 to 2014 May 8. The discussion on the corresponding observing strategy and calibration procedures can be found in Aller et al. (1985) for the UMRAO data, and in Richards et al. (2011) for the OVRO data.

3. Data analysis and Results

3.1. Optical microvariability

In our analysis, we used the $F-$test for assigning the INV detection significance. The $F-$statistic compares the observed variance $V_{\rm obs}$ to the expected variance $V_{\rm exp}$. The null hypothesis of no variability is rejected when the ratio

[TABLE]

exceeds a critical value for a chosen significance level $\alpha$, for a given number of degrees of freedom (DOF) $\nu$; here $V_{t-s}$ is the variance of the ‘target-star’ DLC, $\langle\sigma_{t-s}^{2}\rangle$ is the mean of the squares of the (formal) rms errors of the individual data points in the ‘target-star’ DLC. Since this method requires flux or magnitude estimates along with their error estimates, it is important to determine the photometric errors accurately. As emphasized in several independent studies, the photometric errors returned by $APPHOT$ are significatnly underestimated (Stalin et al. 2004 and references therein). Goyal et al. (2013b) reported the latest attempt to determine this under-estimation factor, $\eta$, using an unprecedented data set consisting of 262 steady star-star DLCs. They find $\eta=1.54\pm 0.05$ which confirms the previous estimates by the same group which were based on much smaller data sets. Goyal et al. (2013b) also showed that the determination of $\eta$ is quite insensitive to the magnitude difference between the pair of objects used for deriving a DLC, as long as it did not exceed 1.5 mag. Thus, $\eta=$1.54 has been used in the present analysis to scale up the IRAF photometric magnitude errors (see also Gopal-Krishna et al. 1995; Villforth et al. 2010a). Note that the standard expression for $F$ is given by $F_{\nu_{1},\nu_{2}}^{\alpha}=\sigma_{1}^{2}/\sigma_{2}^{2}$, where $\sigma_{1}$ and $\sigma_{2}$ are the two variances with the corresponding DOF, $\nu_{1}$ and $\nu_{2}$. In our analysis, we have simplified this definition since $\nu_{1}$ = $\nu_{2}=\nu$ is also the number of DOF for the ‘star-star’ DLC.

The significance level set for a given test determines the expected number of false positives, which is an indicator of the robustness of the test. We have chosen two significance levels, $\alpha=$ 0.01 and 0.05, corresponding to $p-$values of $\gtrsim$ 0.99 and $\gtrsim$ 0.95, respectively. Recall that the smaller the value of $\alpha$ is, the less likely it is for the variability to occur by chance. Thus, in order to claim a genuine INV detection, i.e., to assign a ‘variable’ designation (V), we stipulate that the computed statistic value is above the critical value corresponding to $p>0.99$ (i.e., $\alpha=$ 0.01) for a given degree of freedom ($\nu=N_{p}-1$, where $N_{p}$ stands for the number of data points in a given DLC). We assign a ‘probable variable’ designation (PV) when the computed test statistic value is found to be between the critical values at $\alpha=$ 0.01 and 0.05; otherwise, a ‘non-variable’ (N) designation is assigned to a DLC. All the three DLCs, i.e., $t-s1,t-s2$ and $s1-s2$, are subjected to the $F-$test analysis. In a few cases, the INV status was different for the two blazar-star DLCs, and this indicated a small amplitude variation of one or the other comparison stars. Since such small amplitude variations in star-star DLCs are difficult to ascertain, we only ascribed a “V” status if both blazar-star DLCs gave a “V” status; otherwise, we quote a“PV” or “N” status if the star-star DLC itself turned to be variable. The analysis results are summarized in Table 1.

Following Romero et al. (1999) the peak-to-peak INV amplitude was calculated as

[TABLE]

with $D_{\rm min/max}$ denoting the minimum/maximum in the differential light curve of the source, and $\sigma^{2}=\eta^{2}\,\langle\sigma^{2}_{i}\rangle$ where $\eta=1.54$ and $\sigma_{i}$ is the nominal error associated with each data point.

The INV DC was computed according to

[TABLE]

where $\Delta t_{j}=\Delta t_{j,\,{\rm obs}}\,(1+z)^{-1}$ is the duration of the monitoring session of a source on the $j^{th}$ night, corrected for the cosmological redshift $z$, and $N_{j}$ is set equal to 1 if INV was detected, and otherwise to 0 (Romero et al., 1999; Stalin et al., 2004). This estimate is essentially the ratio of the number of nights a source is found to be variable to the total number of nights it was monitored. Note that, since for a given source the monitoring times on different nights are not necessarily equal, the evaluation of the DC has been appropriately weighted by the actual monitoring duration $\Delta t_{j}$. The computed INV DC for the entire data set, consisting of 25 nights spanning over 18 years (1998–2015) is 22% (37% if three PV cases are included). The INV DC for the stronger nightly variations, where $\psi>3\%$, is 4% (10% if one PV case is included).

In order to validate our analysis procedure, we have performed a sanity check by computing the number of ‘Type 1 errors’, or the false positives, for our data set. A false positive arises due the rejection of a true null hypothesis by a test, when applied to a non-varying DLC (i.e., the inability to discern a non-variable object as non-variable). Assuming a priori that the star-star DLCs are steady, the outcome of the statistical test applied to these should be consistent with the expected number of false positives for the assumed value of $\alpha$. The number of false positives depends only on the number of the star-star DLCs examined and the value of $\alpha$ chosen for the test. Thus, if the number of false positives is found to be significantly different from the expected number, either the test is not robust or the the measurement errors are not properly accounted for (see, e.g., de Diego, 2010; Goyal et al., 2013a, for the robustness of statistical tests in microvariability studies).

We note that for our data set consisting of 25 steady star-star DLCs, the means of the expected numbers of false positives are $\simeq 0.3$ and $\simeq 1.3$ for $\alpha=0.01$ and 0.05, respectively. Since the distribution of false positives is expectedly binomial, for $\alpha=0.01$ the number of false positives should in fact be scattered between 0 and 2, and for most of the cases around $\simeq 0.3\pm 0.5$. Similarly, with $\alpha=0.05$, the number of false positives should lie between 0 and 5, and should largely cluster at $\simeq 1\pm 1$. Meanwhile, the observed numbers of false positives reported by the application of the $F-$test (see column 13 of Table 1) is 1 for $\alpha=0.01$ and 1 for $\alpha=0.05$. The good agreement between the expected and the observed numbers of false positives provides validation for our analysis procedure.

3.2. Multi-wavelength long-term variability

Figure 2 presents the long-term, multi-wavelength, and photo-polarimetric light curve of PKS 0735$+$178. Note that a brief discussion of the long-term optical and radio variability of the source (until 2008) was given in Goyal et al. (2009). The figure includes the 0.1–200 GeV integrated $\gamma$-ray flux $S_{\gamma}$ (panel 2a), the optical R-band photometric flux $S_{\rm opt}$, polarization degree ${\rm PD_{opt}}$, and polarization angle $\chi_{\rm opt}$ (panels 2b, 2c, and 2d, respectively), as well as the radio $S_{\rm rad}$, ${\rm PD_{rad}}$, $\chi_{\rm rad}$, and the spectral index $\alpha_{\rm rad}$ (panels 2e–2h). The radio spectral index, defined here as $S_{\nu}\propto{\nu^{\alpha}}$, was calculated by a linear regression analysis of the flux values at the three GHz frequencies in log-log space (with measurements at various radio frequencies performed within 14 days considered as simultaneous). The optical polarization measurements (prior to 2005) were obtained from Smith et al. (1987), Visvanathan & Wills (1998), Villforth et al. (2009), and Wills et al. (2011). Figure 3 shows the expanded long-term light curve consisting of our newly acquired data for the period 2005–2015.

3.2.1 Power spectral analysis

The PSDs of the observed light curves of PKS 0735+178 were generated using the standard methods outlined in Uttley et al. (2002) and Max-Moerbeck et al. (2014), which can be summarized as follows. The PSD of an evenly sampled light curve $f(t_{i})$ with mean $\mu$ and a total duration of $T$, consisting of $N$ data points, is defined as the rms-normalized periodogram,

[TABLE]

where $|F(\nu_{k})|^{2}$ is the squared modulus of the discrete Fourier transform (DFT) calculated after subtracting the mean from flux measurements. In order to achieve a regular sampling of the light curve at evenly spaced time intervals, when our data were necessarily obtained at irregular intervals, we linearly interpolate between the two consecutive observed data points on the timescales typically 15-20 times smaller than the original observed sampling interval. Reasons for subtracting the mean and performing data interpolation in the case of colored noise-type sources are discussed in Appendix A.

By definition, the integration of the periodogram over positive frequencies yields the total excess variance. Meanwhile, the noise floor levels corresponding to the variability power due solely to statistical fluctuations are estimated following Isobe et al. (2015) as

[TABLE]

In the above, $\sigma_{stat}^{2}=\sum_{j=1}^{j=N}\Delta f(t_{j})^{2}/N$ is the mean variance of the measurement uncertainties on the flux values $\Delta f\!(t_{j})$ in the observed light curve at times $t_{j}$, with $N$ denoting the number of data points in the original light curve. In the case of Gaussian errors, the noise floor level is scaled using typical sampling intervals of 1 day for the long term optical light curve, 3 days for the OVRO light curve, and 15 days for the UMRAO light curves (see also Vaughan et al., 2003, Appendix A). Note that in the case of the lower frequency UMRAO data, in particular the 4.8 GHz ones, the calculated noise floor level is relatively high when compared with the corresponding variability power (see Fig. 5, bottom panel), and is not obviously reflected as a high-frequency plateau in the raw periodogram; this is solely due to the fact that in this case the errors on flux measurements were effectively over-estimated during the first quarter of the survey program.

PSDs generated using this method are distorted due to the effects of the discrete sampling and the finite duration of the light curve, known as aliasing and red noise leakage, respectively. The net effect of aliasing is to fold back the power from high frequencies to lower frequencies; at the same time, the red noise leak causes extra power at higher frequencies through the side-lobes of a sampling window function (Press et al., 1992). The aliasing is reduced for steeper power-law power spectra (Papadakis & Lawrence, 1993; Uttley et al., 2002), while the red-noise leak can be minimized using a proper sampling window function (see Max-Moerbeck et al., 2014, for a detailed discussion). In our analysis of the PSD, we used the ‘Hanning’ window function as it has the lowest side lobes and the fastest fall-off rate as a function of frequency when compared to other frequently employed window functions (see Appendix A for a detailed discussion).

The periodogram obtained using equation (4), with the discrete set of frequencies $\nu_{k}=k/T$ with $k=1,...,N/2$, is known as the ‘raw’ periodogram. It consists of independently distributed $\chi^{2}$ variables with two DOF. This means that the dispersion in each estimate around the true value equals the true value itself, providing a noisy estimate of the spectral power (Papadakis & Lawrence, 1993; Vaughan et al., 2003). To circumvent this problem, we have also constructed ‘binned’ logarithmic periodograms, following the procedure outlined in Papadakis & Lawrence (1993), Vaughan (2005), and Isobe et al. (2015). We binned our periodograms with a constant factor of 1.6 in frequency (and by a factor of 1.5, in case of periodograms derived for the long-term optical light curve) and evaluated the mean power at the representative frequency taken as the geometric mean of each bin. The derived periodogram displays constant variance around the true underlying power spectrum, as if we are observing a noise process following the $\chi^{2}$ distribution with two DOF. For the binned logarithmic periodogram the variance decreases by a factor of $0.310/M$, where $M$ is the number of data points in each bin. Furthermore, the true power spectrum is related to the observed power spectrum as $P(\nu_{k})=(\chi^{2}/2)\,P_{true}({\nu_{k}})$. In order to derive the PSD slope we fit the power law function using a least-squares fit method in log-log space. Since the scatter is multiplicative in the linear space, it is additive in log-log space and identical at each frequency, so

[TABLE]

Hence, the expectation value of the periodogram in log-log space is not the same as the expectation value of the power spectrum: there is a bias between the two values. This bias is a constant because of the shape of the $\chi^{2}$ distribution in log-log space. The expectation value of $\log[\chi^{2}/2]=-0.25068$ (Papadakis & Lawrence, 1993). Therefore, this value is added in the estimate of the binned periodogram.

All the generated PSDs are summarized in Table 2, and presented in Figures 4 (Fermi-LAT), 5 (UMRAO and OVRO), and 6 (optical R-band) for the actual durations of the corresponding light curves, down to the observed sampling intervals. We have not subtracted the constant noise floor level (shown by the dashed horizontal lines in the figures), as some of the data points are below this level. For the analysis, the logarithmic binned PSDs were fitted with a single power law model $P(\nu_{k})\propto{\nu_{k}^{-\beta}}$ using linear regression with weighted error in log-log space; the results of the fitting, along with the errors calculated as the rms residuals between the model and the data, as well as the corresponding variability frequency ranges of the analyzed light curves, are summarized in Table 2.

Figure 6 presents a composite PSD for the optical data set analyzed in this paper, including the long-term monitoring data (1993–2015; 23 years down to nightly sampling), as well as the intra-night data with confirmed INV detections (spanning 4–8 h down to $\sim 10-15$ min sampling; see Table 2). The fitted power-law slope is $2.0\pm 0.1$ for the long-term segment of the optical PSD, and ranges from $\sim 1.5$ up to $\sim 4.0$ for the individual intra-night data sets. The composite optical PSD (long-term + intra-night) covers an unprecedented frequency range of nearly 6 dex for a blazar source. This range, from, $\sim 10^{-9}-10^{-3}$ Hz, is primarily due to the high photometric accuracy achieved in our intra-night monitoring program. In all the cases, however, the normalization of intra-night PSDs turns out to be consistent with a simple extrapolation of the red-noise ($\beta\sim 2$) optical PSD from lower temporal frequencies.

Figure 7 presents the composite multi-wavelength PSDs of PKS 0735+178 corresponding to the variability timescales $>10$ days. This explicitly shows the similarities between the radio and optical bands and the clear difference of the $\gamma$-ray band from the others for which we have enough data to compute sensible PSDs.

4. Discussion and Conclusions

The main findings from our analysis of the multi-wavelength, multi-epoch, photo-polarimetric data for PKS 0735+178, can be summarized as follows:

The INV duty cycle of PKS 0735$+$178 in the extended dataset consisting of 25 nights (22 from our monitoring program and three from the literature; Figure 1), over the time span of 18 years, is DC $\sim 22\%$ for all the variability amplitudes $\psi$, but DC $\sim 4\%$ when only the nights showing $\psi>3\%$ are considered. This value of INV DC is $\sim$ 5 –10 times smaller than that typically observed in blazar sources, indicating that PKS 0735$+$178 remains in a state of optical quiescence on intra-night timescales, in conformity with our previous study (Goyal et al., 2009). 2. 2.

The newly obtained long-term Fermi-LAT and optical data reveal weakly-correlated, a-factor-of-a-few, flux changes in the source on month-like timescales (Figure 3); meanwhile, the archival radio and optical data show uncorrelated, larger-amplitude (factor of several) variability over the timescale of years (Figure 2). 3. 3.

The gathered photo-polarimetric data reveal large variations in the optical PA, traced during the period 2005–2015 on timescale of months to years; those changes seem erratic, with no clear repetitive pattern and no correlation with either the total optical flux, the optical PD, or even with the overall level of the optical polarization degree ($\lessgtr 3\%$; Figure 3). Interestingly, the optical PD of this blazar has frequently been found to be below 3%. 4. 4.

The PSDs of PKS 0735$+$178 at radio (GHz) and optical frequencies, on the timescales from years to weeks, resemble each other closely: both can be well represented by a single power-law function with the slope $\beta\simeq 2$, meaning a “pure red noise” type of the source variability at the lower frequencies of the electromagnetic spectrum. On the other hand, the high-energy $\gamma$-ray PSD, within the corresponding temporal frequency range, is best fitted with a much flatter slope of $\beta\simeq 1$, corresponding to a “pink noise” behavior (Figure 7); no low-frequency flattenings in the PSDs have been detected. 5. 5.

The optical PSDs on intra-night timescales (for the nights during which INV has been detected), are characterized by a range of slopes mostly between 1.5 and 3.3 with the exception of one instance where it was $\sim$4 (however, this was from a shorter observation with relatively long integration times so that it was derived using only 3 points in the binned periodogram). This implies a non-stationarity of the variability process in the source, at least on hourly timescales; in all the cases, however, the normalized intra-night PSDs, is consistent with a simple extrapolation from the red-noise ($\beta\sim 2$) optical PSD observed on longer timescales (Figure 6, Table 2). No high-frequencycut-offs have been detected in the intra-night power spectra down to the noise floor levels.

4.1. Polarization and microvariability

The very low INV DC of PKS 0735+178 reported previously in the literature (Goyal et al., 2009) was puzzling, in view of the persistently high optical polarization degree claimed for the source (Gabuzda et al., 1996; Wills et al., 2011). Our extended monitoring has confirmed the low INV DC but, at the same time, revealed that the ${\rm PD_{opt}}$ of the blazar frequently drops below $3\%$ – for about one-third of the monitoring time, based on the available optical polarimetry.

Traditionally, the degree of optical polarization has been used to differentiate between blazars and LPQs (Angel & Stockman, 1980; Moore & Stockman, 1984). In blazar sources, states with very high $\rm PD_{opt}$ reaching even $\sim 40-50\%$, are not uncommon (e.g., Fan, 2002; Larionov et al., 2013), indicating highly ordered magnetic fields within the jet regions where the observed optical synchrotron emission is being produced; in the case of LPQs, on the other hand, the optical polarization is presumed to be reduced predominantly through the more significant contribution of (unpolarized) thermal emission from an accretion disk to the radiative output of a source at optical wavelengths (e.g., Scarpa & Falomo, 1997). Yet, PKS 0735+178 is a typical example of a BL Lac object lacking any strong emission lines, while those would be expected to be visible in the presence of pronounced disk continuum emission (Nilsson et al., 2012).

Although some blazars are known to alter between the “blazar-like” and the “LPQ-like” optical polarization levels (e.g., Fugmann, 1988), the frequency of such transitions is poorly known, presumably due to the lack of regular polarimetric observations. Only a few of the brightest sources have been subjected to long-term optical photo-polarimetric monitoring. For example, the BL Lac object OJ 287 during the period 2005–2009, displayed $\rm PD_{opt}<3\%$ for only $\lesssim 0.5\%$ of the total observed time span (Villforth et al., 2010b). Hence, one may conclude that PKS 0735$+$178 is a peculiar blazar in showing both a relatively low microvariability duty cycle, and frequent states of low optical polarization.

Based on a comparative study of quasars of the HPQ and LPQ types, Goyal et al. (2012) concluded that the pronounced and frequent INV is correlated with the overall high degree of the optical polarization, rather than with the degree of relativistic beaming, which seems to play a secondary role. Specifically, blazars and quasars with similar degrees of relativistic beaming (i.e., with a similar radio core prominence; e.g., Orr & Browne, 1982), exhibit considerably and systematically different INV DCs. This is in agreement with our new results summarized above for PKS 0735$+$178, suggesting a basic link between the magnetic field ordering in the jet, as reflected in the optical polarization degree, and the production of rapid flux changes.

We finally note in this context that the VLBI radio jets of HPQs show $\chi_{\rm rad}$ of the core/inner jet region parallel to the jet direction, and typically well aligned with $\chi_{\rm opt}$, while the VLBI radio jets of LPQs are characterized by $\chi_{\rm rad}$ misaligned with respect to the jet axis, and uncorrelated with $\chi_{\rm opt}$ (Lister & Smith, 2000). For PKS 0735$+$178, we observe erratic changes of $\chi_{\rm opt}$ between 0 and 180 deg, uncorrelated with $\chi_{rad}$ and with no special relation to the VLBI jet direction in the source (the VLBI jet position angle is $\sim 70$ deg; Britzen et al., 2010). Such a behavior is hard to reconcile with a “grand design” helical magnetic field structure inferred for the source by Gabuzda et al. (2006) and Agudo et al. (2006).

4.2. Power spectra

The LT radio PSDs of blazar sources analyzed in the literature within the temporal frequency range $\sim 10^{-9}-10^{-7}$ Hz, are best represented by single power-laws with slopes $\beta\sim 1.5-2.5$ (Park & Trippe, 2014; Max-Moerbeck et al., 2014). Specifically, the PSD of the BL Lac object PKS 2155$-$304, while consistent with a power-law function $\beta\sim 1.8$, shows a flattening at frequencies below $\sim 10^{-3}$ d*-1*, suggesting a transition from a red noise to a white noise type of variability on the timescale of a few years (Kastendieck et al., 2011). For PKS 0735$+$178, Ciprini et al. (2007) analyzed the long-term optical monitoring data (1970–2004) using the structure function method and found the associated PSD slope $\beta$ to be between 1.5 and 2.0 on the timescales ranging from 33 years down to weeks; at the same time, they found a variety of PSD slopes, $\beta\sim 1.5-2.3$, for densely-sampled light curves from separate observing seasons (covering the timescales of a few months to a few days). Also, the optical PSDs of blazars detected with the Kepler telescope are characterized by $\beta\sim 1.5-2.0$ within the temporal frequency range $\sim 10^{-7}-10^{-4}$ Hz (Revalski et al., 2014; Edelson et al., 2013).

The X-ray PSDs of BL Lacs, on the other hand, are typically consistent with a broken power-law model with the slopes $\beta\sim 2-3$ above the break frequency $\sim 1$ d*-1*, and $\beta\sim 1-2$ below the break (Simonetti et al., 1985; Kataoka et al., 2001; Zhang et al., 2002; Brinkmann et al., 2003; Shimizu & Mushotzky, 2013; Isobe et al., 2015). Unfortunately, PKS 0735$+$178 is quite weak in X-rays, so very few observations of it have been made and it is not possible to compute a sensible X-ray PSD for it. The Fermi-LAT PSD slopes for the bright BL Lacs and FSRQs, have been estimated by Abdo et al. (2010a) as $\beta\sim 1.7\pm 0.3$ and $\sim 1.4\pm 0.3$, respectively, for temporal frequencies from $\sim 10^{-8}$ to $10^{-6}$ Hz (see also in this context Nakagawa & Mori, 2013; Ackermann et al., 2010). Sobolewska et al. (2014), who modeled in detail the 4 yr-long Fermi-LAT light curves of the brightest blazar sources down to week-long sampling intervals, claimed the corresponding PSD slopes to be typically $\beta\simeq 1$, in a very good agreement with our result for PKS 0735$+$178.

Rather limited work has been carried out to compare the PSD properties at different wavelengths of the electromagnetic spectrum for a given blazar source. Most notably, Chatterjee et al. (2008), by comparing the LT PSD of the luminous FSRQ 3C 279 at X-ray (3–20 keV), optical (R-band), and radio (GHz) frequencies, covering the total time span of $\sim 11$ years (temporal frequencies from $10^{-8.5}$ Hz down to $10^{-5}$ Hz), demonstrated that each PSD can be well fitted by a single power-law with $\beta\sim 2.3$ (X-rays), $\sim 1.7$ (optical), and $\sim 2.3$ (radio). A more extensive study by Chatterjee et al. (2012), consisting of six blazars observed with Fermi-LAT and optical/near infra-red telescopes, albeit only for a total duration $\sim 1$ year (temporal frequency range $\sim 10^{-7.5}-10^{-5.5}$ Hz), revealed that, on average, the blazar PSDs at different frequencies are roughly consistent with $\beta\sim 1.6$.

Only a few blazar studies have addressed the issue of the PSD characteristics on intra-night timescales. Papadakis et al. (2003) estimated the optical PSD slopes for BL Lacertae to be $\beta\simeq 1.87\pm 0.16$ on five observing nights within the frequency range $\sim 10^{-4}-10^{-3}$ Hz. In the case of a highly variable BL Lac object 0716+714, Azarnia et al. (2005) noted rather flat optical PSD slopes with $\beta\sim 0.9-1.3$ on 10 monitoring nights for the frequency range $\sim 10^{-5}-10^{-3}$ Hz; meanwhile, Carini et al. (2011) estimated $\beta\sim 2-3$ on five nights and within the same temporal frequency range (see also Mocanu & Sándor, 2012). The wide range of PSD slopes cited above on intra-night timescales for 0716$+$716, implies either some significant inconsistency between the different analysis methods applied, or the non-stationarity of blazar variability on hourly timescales. The latter case echoes our results for PKS 0735$+$178.

Finally, in the very high energy $\gamma$-ray domain, Aharonian et al. (2007) reported the PSD slope of $\beta\sim 2$ for $\nu_{k}\sim 10^{-4}-10^{-2}$ Hz during the famous TeV outburst of PKS 2155$-$304.

The PSD slopes derived from our detailed analysis of the LT high energy $\gamma$-ray (0.1–200 GeV), optical (R-band), and radio (GHz frequencies) data on PKS 0735$+$178, on the timescales ranging from years to weeks/days, suggest that the statistical character of the $\gamma$-ray flux changes is different from that of the radio and optical flux changes: there is increasingly more variability power in $\gamma$-ray fluctuations on comparable temporal frequencies when going to shorter and shorter variability timescales (see Fig. 7). Our finding is in agreement with the Fermi-LAT blazar data analysis presented by Sobolewska et al. (2014). This result is somewhat surprising since, at least within the framework of the standard one-zone leptonic models for the broad-band blazar emission, which is supported to some extent by the multi-wavelength correlations sometimes observed (e.g., Hovatta et al., 2014; Arshakian et al., 2012), the power spectrum of the synchrotron emission component (radio and optical photon energies) should be the same as, or eventually flatter than, the power spectrum of the inverse-Compton component (Finke & Becker, 2014).

As discussed by Kelly et al. (2009), a broken power-law form of the PSD with flat low-frequency segment $P(\nu_{k})\propto const$ and the high-frequency slope $\beta=2$, lacking any peaks indicative of a (quasi-)periodicity, can be understood in terms of a source variability being driven by an underlying stochastic process. In particular, Kelly et al. (2009) proposed that such a variability can be modelled as a first-order autoregressive process (Gaussian Ornstein-Uhlenbeck process; OU for short), in which the source emissivity responds to some input noise (Gaussian white noise, by assumption), with a given relaxation timescale $\tau_{1}$. For such a situation, the source variability at temporal frequencies above the break $\nu_{k_{1}}\equiv(2\pi\tau_{1})^{-1}$ is of the red noise type, and below the break, of the white noise type.

Kelly et al. (2011) discussed a more complex case of a linear superposition of OU processes with two very different relaxation timescales, $\tau_{1}>\tau_{2}$, resulting in the intermediate pink noise ($\beta=1$) segment in the source PSD, in between of the white noise ($\beta=0$) below $\nu_{k_{1}}$ and the red noise ($\beta=2$) above $\nu_{k_{2}}$.

The shallower slopes of the high-energy $\gamma$-ray PSD (as compared to the optical and radio PSDs) could be therefore understood by hypothesizing, following Kelly et al. (2011), that the $\gamma$-ray variability in PKS 0735$+$178 is shaped by a linear superposition of two types of stochastic processes: (i) the same process driving synchrotron (optical and radio) variability, with the relaxation time $\tau_{1}$ larger than decades (since we do not recover the transition from the red noise to the white noise in the optical PSD up to timescales of $\sim 20$ years); and (ii) the additional process relevant only in the inverse-Compton ($\gamma$-ray) domain, characterized by a relaxation timescale $\tau_{2}$ shorter than, or at most of the order of, days (the shortest variability timescales covered by our analysis of the Fermi-LAT PSD).

Note that in this scenario there need not be some instability operating within the jet on the timescale $\tau_{1}$ which is driving the synchrotron variability – the “drivers” can instead be random (stochastic) fluctuations in local jet conditions. Such fluctuations dissipate their energy and accelerate plasma particles, thus creating fluctuations in the distribution of ultra-relativistic jet electrons. The acceleration process itself and the radiative response of the accelerated electrons are, however, delayed with respect to the input perturbations by $\tau_{1}$, so that in the time domain below $\tau_{1}$ the flux changes are smoothed out and damped (forming the red noise segment of the PSD). The jet adjusts, forgetting about the input perturbations, only on timescales longer than $\tau_{1}$, for which the flux changes become uncorrelated (forming the white noise segment of the PSD).

At this point we can only speculate about the nature of the input perturbations inferred above. They could be identified with, for example, local fluctuations in the plasma bulk velocity, leading to the formation of a complex system of internal shocks and plasma turbulence, possibly passing through and interacting with a stationary reconfinement shock further down the jet (Malzac, 2014; Marscher, 2014; Calafut & Wiita, 2015; Pollack et al., 2016). Alternatively, input disturbances could be related to a turbulent magnetic field inherited from strongly magnetized accretion disk and carried out by the jet (see in this context Igumenshchev et al. 2003; McKinney et al. 2012; also Lohfink et al. 2013). In both cases, a relatively long relaxation timescale $\tau_{1}$ would correspond to some global (magneto-)hydrodynamical timescale characterizing an extended blazar emission region.

Meanwhile, in the $\gamma$-ray regime, the additional process should be related to inhomogeneities in the local (radiating fluid rest frame) energy densities of soft photons, inverse-Comptonized to higher energies by the jet electrons in the neighboring cells of the outflow. A natural “relaxation” timescale in this case should be simply a light crossing timescale for the jet region where the bulk of the observed emission is being produced. For $\tau_{2}\lesssim$ day, one then has the corresponding jet radii $R_{2}\sim c\tau_{2}\,\langle\delta_{\rm j}\rangle\lesssim 0.03$ pc, and the distances from the core $r_{2}\sim R_{2}\,\langle\Gamma_{\rm j}\rangle\leq 1$ pc, where, for order-of-magnitude estimates only, we assumed a conical jet with a half-opening angle $\theta_{j}\simeq 1/\langle\Gamma_{\rm j}\rangle$, and the volume-averaged jet bulk Lorentz and Doppler factors, $\langle\Gamma_{\rm j}\rangle\simeq\langle\delta_{\rm j}\rangle\simeq 30$.

In the interpretation mentioned above, which we will discuss in detail in future work, we must revisit the overly simplistic assumption of a single homogeneous emission zone in blazar jets. Instead, we propose that all the observed blazar emission, variable on timescales of years, months, days, and hours — i.e., both the long-term large-amplitude variability and the microvariability — is generated by an underlying single stochastic process (radio and optical frequencies), or a linear superposition of such processes ($\gamma$-ray regime), within a highly non-uniform portion of the outflow extending from the jet base up to the $\lesssim$ pc-scale distances.

Authors wish to thank the referee for making several constructive comments on the manuscript. AG thanks Bindu Rani, Alok C. Gupta, Rumen Bachev, Margo Aller, Hugh Aller, Paul Smith and Svetlana Jorstad for kindly providing data in electronic form. AG, ŁS and MO acknowledge support from the Polish National Science Centre (NCN) through the grant 2012/04/A/ST9/00083. AG also acknowledges partial support from 2013/09/B/ST9/00026 and MS acknowledges the support of 2012/07/B/ST9/04404. VL acknowledges the support of Russian RFBR grant 15-02-00949 and St. Petersburg University research grant 6.38.335.2015. IA acknowledges support by a Ramón y Cajal grant of the Ministerio de Economía y Competitividad (MINECO) of Spain. Acquisition and reduction of the MAPCAT data was supported in part by MINECO through grants AYA2010-14844, AYA2013-40825-P, and AYA2016-80889-P, and by the Regional Government of Andalucía through grant P09-FQM-4784. The MAPCAT observations were carried out at the German-Spanish Calar Alto Observatory, which is jointly operated by the Max-Plank-Institut für Astronomie and the Instituto de Astrofísica de Andalucía-CSIC. PJW is grateful for hospitality at KIPAC, Stanford University, during a sabbatical visit. This research has made use of data from the University of Michigan Radio Astronomy Observatory which has been supported by the University of Michigan and by a series of grants from the National Science Foundation, most recently AST-0607523, and NASA Fermi grants NNX09AU16G, NNX10AP16G, and NNX11AO13G. The OVRO 40m Telescope Fermi Blazar Monitoring Program is supported by NASA under awards NNX08AW31G and NNX11A043G, and by the NSF under awards AST-0808050 and AST-1109911. Data from the Steward Observatory spectropolarimetric monitoring project were used. This program is supported by Fermi Guest Investigator grants NNX08AW56G, NNX09AU10G, NNX12AO93G, and NNX15AU81G. The paper uses optical photometric and polarimetric data from the BU blazar monitoring programme. The research at BU was supported in part by NASA grants NNX14AQ58G and NNX15AR34G. The Fermi-LAT Collaboration acknowledges support from a number of agencies and institutes for both development and the operation of the LAT as well as scientific data analysis. These include NASA and DOE in the United States, CEA/Irfu and IN2P3/CNRS in France, ASI and INFN in Italy, MEXT, KEK, and JAXA in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council, and the National Space Board in Sweden. Additional support from INAF in Italy and CNES in France for science analysis during the operations phase is also gratefully acknowledged.

Appendix A Interpolation of unevenly sampled light curves

Here we discuss in more detail the problem of interpolation of an unevenly sampled light curve of a red noise-type astrophysical source, presenting in particular the results of our simulations using various techniques for estimating source power spectra, and various window functions. We conclude that, even though the distortion of the signal by the interpolation procedure is inevitable, it may be nonetheless properly recognized, while on the other hand no choice of a particular window function in unevenly sampled data can substitute for linear interpolation. That is to say, power spectral density derived for unevenly sampled colored-noise data without linear interpolation, does not resemble the real periodogram regardless on a choice of a window function (see in this context Deeming, 1975; Uttley et al., 2002; Max-Moerbeck et al., 2014).

The model power spectrum derived using Fourier decomposition of a discrete signal sampled $N$ times at fixed $\Delta t$ time intervals, and observed for a finite amount of time $T=N\times\Delta t$, is a convolution of a true power spectrum and a spectral window. In the finite time series, for which the end points of a light curve are not matched (i.e., when the number of periods in the acquisition is not an integer in the frequency domain), variability power “leaks” in the periodogram from low to high frequencies, and this effect is known as a “red-noise leakage”. This leakage can be minimized by multiplying the time series with a proper window function, which reduces the amplitude of the discontinuities at the boundaries of a given time series (Press et al., 1992). On the other hand, due to the discrete sampling of a light curve, for each variability frequency $\nu$ an additional power is added in the derived periodogram at variability frequencies $2\nu_{\rm Nyq}\pm\nu$, $4\nu_{\rm Nyq}\pm\nu$, …, where $\nu_{\rm Nyq}=(2\,\Delta t)^{-1}$ is the Nyquist frequency. The net effect of such “aliasing” is known a priori for evenly sampled data, while it has to be investigated a posteriori (in the Fourier domain of the spectral window) in the case of uneven sampling (Deeming, 1975).

In our test simulations an artificial light curve is generated using the method of Timmer & Koenig (1995), assuming a pure red-noise ($\beta=2$) power spectrum; the data points are evenly sampled with a sampling period of one day, and the total duration of the simulated time series is 1,000 days. The power spectrum is obtained using the discrete Fourier transform (DFT) with three various window functions applied — rectangular, Bartlett (triangular), and Hanning window functions — and later also using the Lomb-Scargle periodogram (LSP) method and the Fourier transform (FT) of the autocorrelation function (ACF).

For completeness, we note that the power spectral density (PSD) of an evenly sampled light curve $f(t_{i})$ is defined as the rms-normalized squared modulus of the DFT,

[TABLE]

computed for evenly spaced frequencies $\nu_{k}=k/T$ with $k=1,...,N/2$, where $T=N(t_{k}-t_{1})/(N-1)$ and $\mu$ is the mean flux. The LS periodogram (Scargle, 1982), on the other hand, computes the power spectrum as

[TABLE]

where $\omega$ is defined by

[TABLE]

so that the data are equivalent to sum of sines and cosines. Finally, the ACF of the analyzed time series is given by

[TABLE]

where $E\{F\}$ is the expectation value of the function $F$, while $\sigma_{f}$ denotes the standard deviation of the flux values around the mean $\mu$, and $\tau$ is the ‘lag’ with which the data points are separated (Edelson & Krolik, 1988). The Fourier transform of this ACF gives the spectral power (Deeming, 1975),

[TABLE]

At first, we consider a light curve with a global monotonic trend, in order to demonstrate the effect of a power leakage from the side-lobes of the window function. The original simulated light curve, i.e. the light curve before subtracting the mean, is shown in Figure  8(a); in the panel (b) of the figure, we present the PSDs derived using the LSP method, and the DFT with rectangular, Bartlett, and Hanning window functions for this original light curve. In Figure  8(c), the analogous PSDs derived for the mean-subtracted light curve are given, all in arbitrary units. As shown, the general shape of the PSD is recovered well in the derived peridograms, with slopes $\beta=1.83-1.94$ estimated by fitting straight lines in the log-log space. Here the difference with the true PSD, $\beta=2$, is exactly due to the red noise leakage effect. In addition, one may see that in the periodograms derived for the original light curve (i.e., the light curve without mean subtraction), ripples arise in the low-frequency segments of the power spectra derived using the DFT methods, depending on the particular shapes of the window functions applied. In principle, this kind of a distortion of the power spectrum should be present only around the 0’th frequency in the power spectrum (as is the case with rectangular window). We found out, however, that with the usage of Bartlett or Hanning windows, lower frequencies (other than the 0’th frequency) are affected as well. These ripples/distortions are minimized in the periodogram derived for the mean-subtracted light curve using the DFT with the Hanning window: as discussed in detail in Max-Moerbeck et al. (2014), the Hanning window has the lowest sidelobe level and higher fall-off rate as compared to rectangular and Bartlett windows.

We repeated the analysis for a light curve without any global monotonic trend. The results of the simulations are presented in Figure  8(d-f). As shown, in this case the red noise leakage effect is less pronounced (although still present), but the low-frequency ripples still arise in the power spectra calculated using the DFT-Bartlett and DFT-Hanning methods for the original light curve. Again, such distortions are minimized in the case of the DFT-Hanning method applied to the mean-subtracted light curve.

Next, we demonstrate that the LS periodogram and Fourier transform of the ACF do not reproduce correctly steep red noise-type power spectra in the case of unevenly sampled data. Generally speaking, in the case of the LS method that is because of a small number of degrees of freedom (DOF) available to characterize the variability at the highest frequencies. On the other hand, in the periodograms derived using Fourier transform of the ACF, one cannot deconvolve the effect of the spectral window, which becomes more and more prominent for less and less evenly spaced data. In order to demonstrate this, we simulated the light curve with $\beta=2$, containing 1,000 data points with a sampling period of one day; we then subtract the mean and keep 30% of the data selected at random times to mimic an unevenly sampled dataset; the simulated light curves are presented in the Figure 9(a). Figure 9(b) shows the corresponding power spectra derived using the LS method and the Fourier transform of the ACF; the Figure 9(c) gives the response of the spectral window for the Fourier transform of the ACF. As shown, the derived PSDs (slopes $\beta=$ 0.13 and $-$0.37 for the LS periodogram and the Fourier transform of the ACF, respectively), are very different from the true PSD ($\beta=2$). The flattening of the derived power spectrum at high variability frequencies in the case of the LS periodogram is due to the aforementioned DOF problem, which introduces artificial power in the high frequency range of the spectrum. In the case of the Fourier Transform of the ACF, which is free of such a “missing data points” problem, the flattening is instead due to a substantial amount of variability power provided by the spectral window. Note in this context that the spectral window does always introduce some extra power in the peridograms derived based on the Fourier transform, but in case of evenly sampled data there is no such extra power at frequencies other than the zeroth Fourier frequency. For unevenly sampled data this extra power becomes, however, significant when compared to the true variability power, particularly in the highest frequency range; as a result the derived power spectra are artificially flattened. Fourier transforming of the ACF instead of the initial light curve does not provide any solution in this respect, since the window function cannot be deconvolved from the derived power spectrum.

Finally, we argue that the only method enabling one to produce a robust derivation of the power spectrum for unevenly sampled data, must rely on a linear interpolation of the light curve. In order to do so, we start from the unevenly sampled time series shown in the left panel of Figure  9, which is then interpolated on a regular grid with a sampling interval of 0.5 day. The Figure 10(a) shows this original unevenly sampled time series along with the linearly interpolated time series; the corresponding power spectra derived using the LS method and the DFT method with three different window functions are presented in the Figure 10(b). The PSD slopes are estimated via linear fitting in the log-log space down to the Nyquist frequency of the original time series ($\sim 6$ days). As shown, all the methods considered reproduce now the true PSD reasonably well. However, keeping in mind the fact that for the time series with a global monotonic trend the DFT-Hanning method provides the best results (in particular against the red noise leakage effect), we conclude that this method is the most robust one.

\twocolumngrid