Blazar variability power spectra from radio up to TeV photon energies: Mrk 421 and PKS 2155-304

We present the results of the power spectral density (PSD) analysis for the blazars Mrk\,421 and PKS\,2155$-$304, using good-quality, densely sampled light curves at multiple frequencies, covering 17 decades of the electromagnetic spectrum, and variability timescales from weeks up to a decade. The data were collected from publicly available archives of observatories at radio from OVRO, optical and infrared (B, V, R, I, J, H, and K-bands), X-rays from the {\it Swift} and the {\it Rossi} X-ray Timing Explorer, high and very high energy $\gamma-$rays from the {\it Fermi} and Very Energetic Radiation Imaging Telescope Array System as well as the High Energy Stereoscopic System. Our results are: (1) the power-law form of the variability power spectra at radio, infra-red and optical frequencies have slopes $\sim$1.8, indicative of random-walk type noise processes; (2) the power-law form of the variability power spectra at higher frequencies, from X-rays to very high energy \,$\gamma$-rays, however, have slopes $\sim$1.2, suggesting a flicker noise type process; (3) there is significantly more variability power at X-rays, high and very high energy $\gamma$-rays on timescales $\lesssim$ 100 days, as compared to lower energies. Our results do not easily fit into a simple model, in which a single compact emission zone is dominating the radiative output of the blazars across all the timescales probed in our analysis. Instead, we argue that the frequency-dependent shape of the variability power spectra points out a more complex picture, with highly inhomogeneous outflow producing non-thermal emission over an extended, stratified volume.


INTRODUCTION
Blazars form an extreme class of active galactic nuclei (AGN), for which the total radiative output is dominated by a highly magnetized and non-stationary nuclear relativistic jet launched from the center of a massive elliptical galaxy (e.g., Urry & Padovani 1995). The blazar family includes flat-spectrum radio-loud quasars (FSRQs) and BL Lacertae (BL Lac) objects, depending on the prominence of emission lines in the optical spectrum. Their spectral energy distribution (SED) usually exhibits two broad peaks, resulting from non-thermal processes. Within the framework of the 'leptonic' scenario, electron-positron (e±) pairs accelerated to GeV/TeV energies are believed to emit radiation from radio-to-optical/UV frequencies (sometimes up to X-rays in the case of BL Lac objects) via the synchrotron process, while the X-ray-to-TeVγ−ray emission arises from inverse-E-mail: arti@oa.uj.edu.pl (AG) Comptonization by the particle population producing the synchrotron radiation of the seed photons which are either produced locally (Synchrotron-Self Compton; SSC), or externally (External Compton; EC) to the jet (e.g., Madejski & Sikora 2016, and references therein). Alternatively, in the 'hadronic' scenario for the blazar emission, the higherfrequency emission peak is believed to originate due to protons accelerated to ≥ PeV-EeV energies and producing γrays via either direct synchrotron process, or meson decay and synchrotron emission by the secondaries produced in proton-photon interactions (e.g., Böttcher et al. 2013, and references therein). The location of the synchrotron emission peak in the SED of BL Lac objects divides them further into high-frequency-peaked BL Lacs (ν peak > 10 15 Hz; HBL), intermediate-frequency-peaked BL Lacs (ν peak 10 14−15 Hz; IBL), and low-frequency-peaked BL Lacs (ν peak < 10 14 Hz; LBL), with FSRQs sharing the space with the LBLs (e.g, Fossati et al. 1998;Ghisellini et al. 2017).
As a class, blazars are known to display strong vari-ability ranging from radio to VHE γ−ray energies and on multiple timescales (e.g., Aller et al. 1999;Wagner & Witzel 1995;Falomo et al. 2014;Abdo et al. 2010) with a factor of few intensity changes on timescales as short as minutes which are found be especially prominent at high energies (Cui 2004;Aharonian et al. 2007;Albert et al. 2007;Ackermann et al. 2016;Zhu et al. 2018). Highly efficient particle acceleration, either due to the formation of shocks and turbulence in the jet flow (Agudo et al. 2018;Marscher 2014;Hughes et al. 2011), or annihilation of the magnetic field lines at the magnetic reconnection sites (e.g., Giannios 2013;Sironi et al. 2015;Petropoulou et al. 2016) are considered to be prime candidates for the energy dissipation processes. The origin of rapid flux changes and particularly their relation to those observed on longer timescales of days-to-years is still widely debated. In this regard, the large amplitude variability on the shortest timescales poses additional challenge because jet bulk Lorentz factors needed to overcome the photon opacity barrier (Γ > 50; e.g., Begelman et al. 2008) are rather too extreme to be reconciled with the currently favored models of AGN jet acceleration (McKinney & Blandford 2009).
The power spectral density (PSD) of a blazar light curve typically has a power-law shape, P(ν k ) ∝ ν −β k , where ν k (≡timescale −1 ) is the temporal frequency and β is the slope. This indicates that the blazar variability is a stochastic processes. Moreover, the typical slopes, β ∼ 1-3, on timescales from years (sometimes even decades) down to days (and in some cases, minutes), indicate that this stochastic process is of correlated colored-noise type (β = 1; flicker/pink and β = 2; random-walk/red), unlike the uncorrelated whitenoise type process which give rise to flatter (β ∼0) PSD slopes (Press 1978). Whether the stochastic process remains stationery on the timescales probed is crucial for understanding the physical processes operating on widely different spatial scales (e.g., Kushwaha et al. 2016;Giebels & Degrange 2009). Studies of different blazars performed at different wavelengths have indeed revealed breaks in the slopes of PSDs, which could be attributed to the size of the emission zone or, alternatively, to cooling timescales of the particles (e.g., Abramowski et al. 2010;Sobolewska et al. 2014;Isobe et al. 2015;Kastendieck et al. 2011;Abdalla et al. 2017;Ryan et al. 2019;Żywucka et al. 2020). The emerging picture of the broad-band blazar variability is that at higher energies (i.e., X−ray, and GeV and TeV γ−rays) is characterized by flicker/pink-noise processes (i.e., Abdalla et al. 2017;Abdo et al. 2010;Isobe et al. 2015), whereas at lower energies (GHz band radio and optical frequencies) damped/red-noise type processes dominate (i.e., Ciprini et al. 2007;Kastendieck et al. 2011;Park & Trippe 2014;Max-Moerbeck et al. 2014b;Nilsson et al. 2018).
In Goyal et al. (2017Goyal et al. ( , 2018, we presented the first attempts to uniformly characterize the statistical properties of the variability at radio-to-γ−rays, on timescales ranging from decades to weeks, by analyzing the GHz band radio, optical, X-ray and γ−ray light curves for the well-known blazars PKS 0735+178 and OJ 287. Our results conform to the findings mentioned above: the PSD shape at higher energies (β ∼1) differs markedly from that at lower energies (β ∼ 2). We interpreted our findings in terms of a single stochastic process, or a linear superposition of such, giving rise to the observed broad-band variability. We emphasized, that different statistical characters of the low and high energy variability, is difficult to explain within the framework of the widely discussed 'one-zone' models, which invoke the same particle (electron) population within a single well-defined and homogeneous emission volume to account for both synchrotron and IC emission components, and in the framework of which one should expect the same slopes of the variability PSDs across all the wavelengths (Finke & Becker 2014. Even though our results find support from the lack of clear correlations between the variations at different frequencies as inferred from cross-correlation studies of blazar light curves (e.g., Max-Moerbeck et al. 2014a;Lindfors et al. 2016), we note that many studies do find correlated variability among different frequencies, along the expectations of the widely used SSC scenario invoked for blazar sources (e.g., Hovatta et al. 2015;Fuhrmann et al. 2016;Ramakrishnan et al. 2016) Due to extreme photon deficiencies at TeV energies, augmented by the absorption process on the extragalactic background light for cosmologically distant sources, and sensitivity constraints of the currently operating Cherenkov telescopes with limited duty cycles, so far 73 blazars have been detected in the TeV range 1 ; note that such problems are much less severe in the case of the all-sky monitoring by the Fermi-LAT at GeV energies (Ackermann et al. 2015). Blazars are often targeted at TeV energies during short flares from lower energies; therefore, the observations are usually biased towards the flaring states (see, Ahnen et al. 2017;Şentürk et al. 2013, and references therein). The blazars Mrk 421 and PKS 2155−304 are exceptions to this rule: thanks to their relatively high flux at TeV energies, these blazars could be successfully monitored on a nightly basis by the Very Energetic Radiation Imaging Telescope Array System (VERITAS) Acciari et al. 2014) and the High Energy Stereoscopic System (H.E.S.S.) observatories -2011Abdalla et al. 2017). For this reason, the TeV measurements in these cases can be taken as representative for energy dissipation processes operating in blazar jets at the highest energies in general, and not only during the highest flux states. Therefore, we focus in the present study on determination of the statistical properties of the stochastic processes giving rise to broadband variability up to TeV γ−ray energies, using good-quality, roughly decadelong, multiwavelength light curves.
For this, we have assembled from the publicly available archives of both ground-based and space-borne observatories, light curves at very high energy (VHE) γ−rays from the VERITAS and the H.E.S.S.(>200 GeV; Acciari et al. 2014;Hinton 2004), high energy (HE) γ−ray from the Fermi-Large Area Telescope (LAT (0.1-300 GeV); Thompson 2004), X−rays from the Rossi X-ray Timing Explorer (RXTE) Proportional Counter Array (PCA (3-20 keV); Bradt et al. 1993) and from the Swift-X−Ray Telescope (XRT (0.3-10 keV); Gehrels et al. 2004), infra-red (IR; J, H, and K-band), optical (B, V, R, and I-band) and 15 GHz radio frequency from the Owens Valley Radio Observatory (OVRO). This has enabled us to present, for the first time, the PSD analysis of the multiwavelength light curves which defines/constrains the nature of stochastic variability of a blazar all the way up to VHE γ−ray energies and on timescales ranging from weeks to a decade.
Mrk 421 is an archetypal blazar which due to its proximity (J2000.0 R.A. = 11 h 04 m 27. s 314, Dec. =+38 • 12 31. 80; redshift = 0.03002; de Vaucouleurs et al. 1991) and the high flux, was the first blazar to be detected at TeV energies by the 10-m Whipple telescope (Punch et al. 1992). 2155−304, (J2000.0 R.A. = 21 h 58 m 52. s 065, Dec. =-30 • 13 32. 11; z = 0.116; Bechtold et al. 2002), is one of the first few blazars to be detected at VHE γ−rays (Chadwick et al. 1999). Mrk 421 has been a subject of intense monitoring ever since its discovery in the VHE band, in particular, the relation between the variations at X-ray and VHE γ−rays have been probed widely (B lażejowski et al. 2005;Fossati et al. 2008;Hovatta et al. 2015;Ahnen et al. 2016;González et al. 2019). For the blazar Mrk 421, no quasi-periodic oscillations (QPOs) in the decade-long optical and HE γ−ray light curves were reported in the study of Sandrinelli et al. (2017) ; 2004-2012), partial HE γ−ray (Fermi-LAT; and X−ray (Swift-XRT;2006 light curves have been presented in Abdalla et al. (2017) and Kapanadze et al. (2014). PKS 2155−304 is also known for a log-normal distribution of TeV flux densities, and also for a change in the PSD slope, such that β ∼ 1 and ∼ 2 on timescales longer and shorter than 1 day, respectively). This indicates a dominance of multiplicative over additive energy dissipation processes manifesting in the VHE regime (Aharonian et al. 2007;Abdalla et al. 2017;Chevalier et al. 2019). Interestingly, correlated variability between VHE γ−ray and X−ray (Costamante 2008), as well as VHE γ−rays and optical frequencies have also been claimed for this blazar (Aharonian et al. 2009). Recently, hints of QPOs, with a characteristic timescale of ∼1-2, yrs have been noted in the HE γ-ray and optical light curves (see, Sandrinelli et al. 2016;Zhang et al. 2017;Chevalier et al. 2019). Possible quasi-periodicities on similar timescales have also been reported in the decade-long optical/IR light curves (e.g., Sandrinelli et al. 2014b;Zhang et al. 2014). In addition, QPOs on comparatively shorter timescales have been claimed for this blazar, in its X−ray light curve (∼4.6 hrs ; Lachowicz et al. 2009) and in the polarized optical emission (∼13 and 30 min ; Pekeur et al. 2016).
In Section 2 we describe in more detail all the data assembled here, and the data reduction procedures. The data analysis and the results obtained are outlined in Sections 3, and 4, respectively, while a discussion and our main conclusions are presented in Section 5.  imaging Cherenkov telescope and the VERITAS array at energies >400 GeV for the period 1995-2009 (Acciari et al. 2014). PKS 2155−304 has been a target of regular monitoring at the H.E.S.S. observatory due to its favorable location in the southern hemisphere and relatively high flux which, together with the telescope's large collecting area enables relatively high accuracy of measurements on a nightly basis. Fig. 2 Abdalla et al. 2017). We note that the presented light curve does not include the portion of the data when the blazar underwent an extreme flare when its flux level increased by a factor of ∼100 (Aharonian et al. 2007), compared to average flux levels of this blazar reported here (Abdalla et al. 2017). The details regarding the observations and data analysis can be found in Acciari et al. (2014) and Abdalla et al. (2017).

HE γ−rays: Fermi-LAT
Fermi-LAT data for the field containing the targets were analysed for the 0.1 -300 GeV band, with an integration time of seven days from August 2008 until March 2018 (Atwood et al. 2009). The details of the LAT data analysis can be found in Goyal et al. (2017Goyal et al. ( , 2018 and below we briefly recall the main features. We have generated the source light curve by performing the unbinned likelihood analysis using Fermi ScienceTools 10r0p5 with latest instrument response function p8r2 source v6 source event selection and zenith angle < 90 •3 . All photons were counted within 20 • region centered at the blazar position to account for the broad point spread function (PSF) of the Fermi-LAT at the desired energies (0.1-300 GeV). The analysis method 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 • radius for each time bin using tasks 'gtltcube' and 'gtexpmap'. Next, we computed the diffuse source response (task 'gtdifrsp') and finally modeled the data with the maximum-likelihood method (task 'gtlike'). In this last step, we included all point sources togather with the targets, Mrk 421 and PKS 2155−304 (189 and 167 other point sources, respectively) inside the ROI from the 3FGL catalog, as well as the standard templates for diffuse emission from our Galaxy (gll iem v06.fits) and the isotropic γ-ray background (iso p8r2 source v6 v06.txt) (Acero et al. 2016). To account for the contamination in the source flux due to variable point-sources within the broad PSF of the LAT (∼5 • at 100 MeV), we checked for variable point sources within 5 • of the target from the 3FGL source-list in the Fermi All-sky Variability Analysis (FAVA; Abdollahi et al. 2017). In the modeling, we fixed the photon indices and fluxes of all the point sources within the ROI other than the target and the FAVA sources at their 3FGL values. Lastly, the γ-ray spectrum was modeled with a power-law and logparabolic functions for Mrk 421 and PKS 2155−304, respectively, as recommended by the 3FGL catalog by keeping the normalization and the spectral parameters free. Finally, a criterion of test statistic TS ≥ 10 was set to consider a measurement to be a successful detection ≥ 3σ (Abdo et al. 2009). The generated light curves consisting of successful detections are shown in Fig. 1(b) and 2(b) for Mrk 421 and PKS 2155−304, respectively.

X−rays: RXTE-PCA and Swift-XRT
We downloaded all of 1,207 individual spectra for Mrk 421 from the RXTE quick look analysis from the heasarc website 4 . The spectra were binned for 20 points in energy bins and fitted with a simple power-law in xspec with a χ 2 statistics. In the spectral analysis, the Galactic absorption corresponding to neutral hydrogen column density (N H, Gal = 1.91 × 10 20 cm −2 ) were fixed in the direction of the source (Kalberla et al. 2005). The unabsorbed energy fluxes were obtained within 3-20 keV energy range. The measured fluxes were then averaged using 1 d bins. Fig. 1(c) shows the resulting light curve. For the blazar PKS 2155−304, the 457 individual spectra are sporadically observed between 1995-2012, therefore, we do not include them in the present analysis.
For the targets, we have analysed the archival data from the Swift-XRT observatory which consists of 1,070 (Mrk 421) and 214 (PKS 2155−304) pointed observations made between September 2005 and January 2019. The details of the XRT data analysis can be found in Goyal et al. (2018) and below we briefly recall the main features. For the data reduction, we used the latest version of the calibration database (CALDB) and version 6.19 of the heasoft package 5 . For each data set, we used the level 2 cleaned event files of the 'photon counting' (PC) and 'window timing' (WT) data acquisition mode generated using the standard xrtpipeline tool. For the PC mode data, the source and background spectra were generated using a circular aperture with appropriate region sizes and grade filtering using the xselect tool. The source spectra were extracted using an aperture radius of 47 , centered at the source position, while four source-free regions of 118 radius each were used to estimate the background spectrum. For the WT mode data, the source region of 47 circular aperture and an annular region with inner and outer radii of 187 and 281 , respectively, for background estimation were selected to estimate the source and the background spectra. Appropriate corrections were made in the 'BACKSCAL' keyword in the spectra to account for the different sizes of source and background regions, following the XRT tutorial 6 . The ancillary response matrix was generated using the task xrtmkarf for the exposure map generated by xrtexpomap. All the source spectra were then binned for 20 points and corrected for the background using the task grppha. All the observations were checked for the recommended pile-up limit for the PC (∼0.5 counts s −1 ) and WT (∼100 counts s −1 ) modes. In almost all the PC mode data, the count rate exceeded the pile-up limit for both blazars and we made appropriate corrections to the source count by modeling the XRT PSF with a Kings function following the standard procedure 7 . On a few occasions, the source count exceeded the pile-up limit for the WT mode data for Mrk 421 and appropriate corrections, i.e., rejecting the central 2-pixels while extracting the source spectrum, were applied following Romano et al. (2006). In none of the observations, however, did the source count rate exceed the recommended pile-up limit for PKS 2155−304.
For each exposure, we used routines from the X−ray data analysis software ftools and xspec to calculate and to subtract the X−ray background model from the data. Spectral analysis was performed between 0.3 and 10 keV energy range by fitting a simple power-law moderated by the Galactic absorption with the corresponding N H, Gal fixed to 1.91 ×10 20 and 1.48 × 10 20 cm −2 in the directions of the blazars Mrk 421 and PKS 2155−304, respectively (Kalberla et al. 2005). We used the 1 d binning interval to average the unabsorbed 0.3-10 keV fluxes, to construct the source light curve for Mrk 421 and PKS 2155−304, respectively (panel d of Fig The B, V, R, J, and K band SMARTS light curves for the source covers the period from ∼2008 until 2015. The V, R, I, J, H and K band light curves provided by S14 span 2005 through 2012. We refer the reader to Bonning et al. (2012) and Sandrinelli et al. (2014a) for details regarding the monitoring programmes and the data analyses. The data at V, R, J, and K bands have a great deal of overlap between the two programmes and they complement each other well at V, R, J frequencies, with the pairs of magnitudes within 0.1 mag on the same nights, consistent with the calibration uncertainties provided by the two programmes. Therefore, we have combined the two datasets at V, R, and J frequencies to produce the daily binned light curves. The K-band light curve from the SMARTS programme has a magnitude offset ≥0.1 mag from that of provided by S14 on the same nights. On visual inspection, the K-band SMARTS light curve appears more variable than the S14 light curve on similar timescales. Therefore, in this study, we present the K-band PSD provided by S14 data alone. Next, all the B, V, R, I, J, H, and K magnitudes were converted to fluxes using m zero ×10 −0.4×m , where m zero (=4063 Jy, 3636 Jy, 3064 Jy, 2635 Jy, 1590 Jy, 1020 Jy, and 600 Jy for the B, V, R, I, J, H, and K-bands, respectively) refers to zero point magnitude flux of the photometric system (Glass 1999 Uttley et al. (2002). We describe the main features below. The periodogram of an evenly sampled light curve f (t i ) observed at discrete times t i , consisting of N data points and the total monitoring duration T, is defined as the squared modulus of its DFT, where the mean µ is subtracted from the flux values, f (t i ), in order to remove the zero frequency power. The DFT of a real function is a complex symmetric function with values at both negative and positive Fourier frequencies, ranging from −N/2 + 1, −N/2 + 2, .., −1, 0, +1, +2, ..., +N/2 − 1, +N/2, in the case of an even number of data points (Press et al. 1992). The DFT is computed for evenly spaced frequencies ranging from total duration of the light curve down to Nyquist sampling frequency (ν Nyq ). Specifically, the frequencies corresponding to ν k = k/T with k = 1, ..., N/2, ν Nyq = N/2T, and T = N(t k − t 1 )/(N − 1) are considered. The power P(ν k ) at the given frequency ν k is then obtained as the rms-squarednormalized periodogram The normalized periodogram as defined in Eq. 2 corresponds to total excess variance when integrated over positive frequencies. Moreover, a constant amount of power (corresponding to white noise type process; β ∼ 0) is expected to contribute in the estimated variability power due to measurement uncertainities alone. Such, 'noise floor levels', are estimated as (e.g., Isobe et al. 2015;Vaughan et al. 2003) where, σ 2 st at = j=N j=1 ∆ f (t j ) 2 /N is the mean variance of the measurement uncertainties on the flux values ∆ 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. Because the observed light curves are not evenly sampled, we scale the noise floor level using typical (mean) sampling intervals (see also, Vaughan et al. 2003). The above described method requires, however, a regular sampling of the light curve at evenly spaced time intervals, otherwise the spectral window function corresponding to the sampling times gives non-zero response in the Fourier-domain(e.g., Deeming 1975). This will result in false powers in the DFT of the light curve (see Appendix A for discussion). Therefore, in order to perform the DFT of the observed light curves, we obtained the regular sampling only by linearly interpolating between the two consecutive observed data points on timescales typically 5-7 times smaller than the original observed sampling interval (see Table 1). The distortions in the PSDs due to the discrete sampling and the finite duration of the light curve are described in detail in the Appendix of Goyal et al. (2017) and the references therein. In our analysis, the resulting PSDs are generated down to the typical (mean) Nyquist sampling frequency of the observed light curves using the 'Hanning' window function (see also Max-Moerbeck et al. 2014b).
The periodogram obtained using equation (2), is known as the 'raw' periodogram. By definition, it consists of independently distributed χ 2 variables with two DOF as it is the sum of squares of Gaussian distributed variables (Timmer & Koenig 1995). This means that each PSD estimate has a standard deviation around the true value that equals the value itself, providing a noisy estimate of the spectral power (Papadakis & Lawrence 1993;Vaughan et al. 2003). Therefore, in order to obtain a reliable estimate of spectral power, a number of PSD estimates should be averaged where the mean spectral power estimate can resemble a Gaussian distribution provided enough number of data points are available for averaging. According to Papadakis & Lawrence (1993), averaging about 20 data points in logarithmic space constrains the PSD estimate within 1σ confidence interval. However, following this method would essentially remove a large number of data points from lower frequencies, instead, we bin the 'raw' periodograms by averaging with a constant factor of 1.6 in frequency range (i.e., with a fixed step in logarithmic space) and evaluating the mean power at the representative frequency taken as the geometric mean of each bin, following Isobe et al. (2015) and Goyal et al. (2017). This value is chosen to have at least two periodograms in each frequency bin except for the first bin.
The computed periodograms P(ν k ), for a noise-like process, are scattered around true value following a χ 2 distribution with 2 degrees of freedom (Papadakis & Lawrence 1993;Timmer & Koenig 1995;Vaughan et al. 2003). Therefore, the true power-spectrum is given as The transformation to log-log space, offsets the the observed periodograms as: This offset is constant and is equal −0.25068 which is the expectation value of χ 2 distribution with 2 DOF in loglog space (Papadakis & Lawrence 1993). Finally, this value is added in computing the binned periodogram estimates.
Since the aim of the present study is to derive shapes of PSDs, the best-fit PSD parameters are determined using Monte Carlo (MC) simulations, following the 'power spectral response' (PSRESP) method, introduced by Uttley et al. (2002) and followed by many others, including Chatterjee et al. (2008) proach, a large number of light curves are simulated with a known underlying power-spectral shape. Each light curve is then rebinned to the observed sampling pattern and interpolated to have even sampling for the DFT application.
The DFT of such light curve gives the distorted PSD due to various effects of rebinning, red-noise leak and aliasing.
Averaging large number of such PSDs gives the mean of the distorted model (input) power spectrum. The standard deviation around the mean gives errors on the modeled power spectrum. The goodness of fit of the model is estimated by computing two functions, similar to χ 2 , defined as: [log 10 P sim (ν k ) − log 10 P obs (ν k )] 2 ∆log 10 P sim (ν k ) 2 and, [log 10 P sim (ν k ) − log 10 P sim,i (ν k )] 2 ∆log 10 P sim (ν k ) 2 where, log P obs and log P sim,i are the observed and the sim-ulated Log-binned periodograms, respectively. log P sim and ∆log P sim are the mean and the standard deviation obtained by averaging 1,000 PSDs; k represents the number of frequencies in the log-binned power spectrum (ranging from ν min to ν max ), while i runs over the number of simulated light curves for a given β.
Where the χ 2 obs determines the minimum χ 2 for the model compared to the data and the χ 2 dist values determine the goodness of the fit corresponding to the χ 2 obs . Note that χ 2 obs and χ 2 dist are not the same as that of standard χ 2 distribution because log 10 P obs (ν k )'s are not normally distributed variables since the number of power spectrum estimates averaged in each frequency bin is small (<20 for the first few frequency bins). Therefore, reliable goodness of fit is computed using the distribution of χ 2 dist values. For this, the χ 2 dist is sorted in ascending order. The probability, or p β , that a given model can be rejected is then given by the percentile of χ 2 dist distribution above which χ 2 dist is found to be greater than χ 2 obs for a given β (also known as the success fraction; Chatterjee et al. 2008). A large value of p β represents a good-fit in the sense that a large fraction of random realizations of the model (input) power spectrum are able to recover the shape and slope of the intrinsic (of which the observed PSD makes one realization) PSD. Therefore, this analysis essentially uses the MC approach toward a frequentist estimation of the quality of the model compared to the data. This is a well-known approach to describe the goodness of fit in the absence of well-understood fit statistics (see, for details, Press et al. 1992).
In the present study, the light curves are simulated following the method of Emmanoulopoulos et al. (2013). This method has the advantage of preserving the probability density function of the flux distribution as well as the underlying power spectral shape, and not just the shape as given by Timmer & Koenig (1995). In the simulations, we have assumed single power-law PSDs with a given β (to reproduce the PSD shape) and supplied mean and standard deviation of the logarithmically transformed flux values (to reproduce the flux distribution) which is, in general, a valid assumption for blazar light curves (Abdalla et al. 2017;Chevalier et al. 2019;Liodakis et al. 2017;Kushwaha et al. 2017). In our simulations we have computed the mean and the standard deviation σ directly from the data as opposed to obtaining them from fitting a Gaussian function to the distribution. The computed mean and σ from the data were found to be within 95% confidence intervals of the mean and σ obtained from the fitting. This is because the analysed light curves are relatively well-sampled and the fluxes are obtained with good precision (few percent measurement accuracies); therefore, any distortions due to limited sensitivity and finite sampling of the light curve are not significant. Each simulated light curve was then rebinned to have the same sampling pattern as that of the observed light curve. Finally, the measurement errors in the simulated flux values were incorporated by adding a Gaussian random variable with mean 0 and standard deviation equal to the mean error of the measurement uncertainties on the observed flux values (Meyer et al. 2019). In such a manner, 1,000 light curves are simulated in the β range 0.2 to 3.0, with a step of 0.1 for each observed light curve. The periodograms are derived for each simulated light curve in the same manner as that of the observed light curve. The probability distribution curves as a function of β for the analysed PSDs are given in the Appendix B. The best-fit PSD slope for the observed PSD is given by the one with the highest p β value. The errors on the best-fit PSD slope is obtained by fitting a Gaussian function to the p β curve and gives the full-width at half maximum (FWHM; Chatterjee et al. 2008;Bhatta et al. 2016a). This gives a roughly 98 percent confidence interval on the best-fit PSD slopes.
All observed PSDs, along with the best-fit PSDs and the maximum p β , are summarized in Table 1. The PSDs are displayed in Figs. 3 and 4 for the actual duration of the corresponding light curves, down to the observed (mean) sampling intervals. In our analysis, we have not subtracted the constant noise level (shown by the dashed horizontal lines in the figures), as some of the data points are below this level. In all the cases, the probability is higher than 10 percent (except for the OVRO light curve of the blazar Mrk 421 where it is 8.4 percent), meaning that the rejection confidence (1-p β value) is lower than 90 percent for the modeled best-fit PSDs. This means that the single powerlaw PSD shape provides a good fit to the PSDs studied here. Furthermore, we also compute the square fractional variability by multiplying ν k with the corresponding P(ν k ) for the best-fit spectral shapes; such estimates are equivalent to fractional variability, F var = r ms mean , where rms is the standard deviation of the light curve (Vaughan et al. 2003). In such a representation, one can explicitly compare variability amplitudes at different wavelengths on different variability timescales.

RESULTS
Here we present the results of our PSD analysis of multiwavelength light curves of the TeV blazars Mrk 421 and PKS 2155−304 on timescales ranging from about a decade down to weeks or days. For the analysis we have used goodquality, roughly decade-long light curves sampled around once a week at 12 different frequencies, covering ∼17 decades of the electromagnetic spectrum. Our main findings are: (i) For Mrk 421, the VHE, HE γ−ray and X-ray PSDs on the timescales from years to months, can all be well represented by a single power-law function with the slopes β 1.1 ± 0.5, 1.1 ± 0.4, 1.1 ± 1.6 and 1.3 ± 0.7, respectively (panels a, b, c and d of Fig. 3 and Table 1). Similarly, for PKS 2155−304 and on the corresponding timescales, the PSDs are represented by a single power-law function with the slopes β 0.6 ± 1.4, 1.3 ± 0.6, and 1.3 ± 2.1, respectively (panels a, b and c of Fig. 4 and Table 1). The variability thus is indicative of essentially "pink/flicker noise" type at these higher frequencies of the electromagnetic spectrum. We note that the PSD slopes derived here for the Xray light curves are in good agreement with β ∼1.2-1.5 reported in Isobe et al. (2015); Chatterjee et al. (2018) for the blazar Mrk 421. Similarly, the PSD slopes derived at VHE and HE γ−ray light curves by means of the DFT method, are in good agreement with β ∼1.2 and ∼1.1 reported in (3) the total duration of the observed light curve; (4) the number of data points in the observed light curve; (5) the minimum sampling interval for the observed light curve; (5) the maximum sampling interval for the observed light curve; (7) the mean sampling interval for the observed light curve (light curve duration/number of data points); (8) the sampling interval for the interpolated light curve; (9) the noise level in PSD due to the measurement uncertainty; (10) the temporal frequency range covered by the binned logarithmic power spectra; (11) the best-fit power-law slope of the PSD along with the corresponding errors representing 98 per cent confidence intervals (see Sec. 3); (12) corresponding p β .  (Table 1). Abdalla et al. (2017), who used other methods for the PSD characterization.
(ii) The gathered optical-IR monitoring data for PKS 2155−304 reveal large changes in intensities on vastly different timescales. The roughly daily sampled (barring the interval when the target was close to the Sun), light curves at B, V, R, I, J, H, and K bands closely resemble one another, with a factor of a few changes in intensity on weeks/months timescales, and a few percent on the timescales of days. The overall optical and IR PSDs, derived on the timescales from years to weeks, can be well represented within errors by a single power-law function of  (Table 1).
slope β 1.8, meaning a "pure red noise" type of the source variability close to peak of the synchrotron component of the electromagnetic spectrum (panels d, e, f, g, h, i, and j of Fig. 4, and Table 1).
(iii) The PSD slope derived using the 15 GHz radio light curve for the blazar Mrk 421 on timescales ranging from years to weeks is β ∼1.6±0.3, indicative of a statistical char-acter of red noise type of source variability at these frequencies (panel e of Fig. 3 and Table 1).
(iv) An explicit comparison of squared fractional variability (i.e., ν k P(ν k ) versus ν k ) between higher (VHE and HE γ−ray) energy bands and lower (radio-to-optical) of the electromagnetic spectrum, reveals more variability power at the higher frequency bands on timescales ≤100 days (panels a and b of Fig. 5). Simiar behaviour is shown by variability at X-rays energies for Mrk 421 as compared to radio energies (panel a of Fig. 5). The X-ray PSD for PKS 2155−304 could not be constrained on such short variability timescales due to the very sparse sampling of the X-ray light curve.
(v) The PSDs generated using our long-duration optical, IR, and Fermi-LAT light curves do not reveal the QPOs reported in Sandrinelli et al. (2016) and Zhang et al. (2017) for PKS 2155−304 as the observed power spectra are distributed within 1σ confidence regions of best-fit PSDs. This could be either due to the marginal significance of the detected features ( 3σ), or it could be indicative of a transitory nature of the QPOs (see, Vaughan et al. 2016;Bhatta et al. 2016b, for a detailed discussion). We note, however, that we have not followed Lomb-Scargle periodogram or Weighted Wavelet Z-transform methods which are used to detect periodicities in the original publications which could potentially cause the difference.

DISCUSSION AND CONCLUSIONS
Based on our PSD analyses of multi-band (GHz band radio, R-band, and Fermi−LAT) light curves of the two LBLs, namely PKS 0735+178 and OJ 287 (Goyal et al. 2017(Goyal et al. , 2018, we have argued that the observed broad-band emission of the blazars is generated in the extended volumes of structured, highly turbulent jets (as opposed to the scenarios postulating a well-defined single "emission zone"). That is because the Fourier decomposition of the analysed light curves did not reveal any single well-defined variability timescale, but instead a wide range of such, manifesting as the colored-noise type of the derived power spectra. We emphasize, that a single/broken power-law shape of the blazar PSDs is a general finding of the timing analysis performed by various authors (e.g., Kataoka et al. 2001;Abramowski et al. 2010;Abdo et al. 2010;Kastendieck et al. 2011;Sobolewska et al. 2014;Park & Trippe 2014;Isobe et al. 2015;Kushwaha et al. 2016;Abdalla et al. 2017;Chatterjee et al. 2018), and that the presence of QPOs in blazar lightcurves is still an open issue due to typically low significance of the claimed features (see in this context Covino et al. 2019, and references therein; also, see the discussion on the 12-year-long periodicity in the optical light curve of OJ 287 in Appendix C of Goyal et al. 2018). Moreover, in order to account for the specific colours of the noise observed, in Goyal et al. (2017Goyal et al. ( , 2018 we have proposed that the variability in the synchrotron-dominated spectral region is driven by a single stochastic process with relatively long relaxation timescale τ 1 100 days, whereas in the IC-dominated spectral domain, by a linear superposition of stochastic processes with vastly different relaxation timescales ranging from τ 1 down to τ 2 1 day. In particular, following the discussion in Kelly et al. (2009Kelly et al. ( , 2011, we speculated that stochastic (Gaussian) fluctuations in the local jet conditions (e.g., velocity fluctuations in the jet plasma, and/or magnetic field variations) lead to energy dissipation over a wide range of spatial scales. The resulting acceleration of jet particles, and their following radiative response, is however delayed with respect to the input perturbations, and this forms the damped/red-noise (i.e., β ∼ 2) segment of the PSD at the synchrotron-dominated frequencies on the timescales shorter than τ 1 , which may be identified with some global MHD timescale characterizing the extended jet volume. Over the IC-dominated spectral segment, on the other hand, due to density inhomogeneities in the local photon populations available for the upscattering, the PSD is modified further to form the flickering/pink noise (β ∼ 1) down to the relaxation timescale τ 2 , which should be of the order of the maximum light travel timescale across the jet where the observed emission is being produced.
Our new results for the HBLs Mrk 421 and PKS 2155−304 are based on good-quality, densely sampled, decade-long multiwavelength light curves (except of X-rays in case of PKS 2155−304) on timescales ranging from weeks to a decade. We find that the PSDs above the noise floor level (which becomes dominant at high temporal frequencies), are well represented by single power-laws with slopes β ∼ 2 at synchrotron-dominated frequencies (radio and B, V, R, I, J, H, and K-bands), and β ∼ 1 at the IC-dominated frequencies (i.e., VHE and HE γ− rays). The power-law slope is also β ∼ 1 at X-rays, even though in these sources, within the single-zone SSC scenario, the keV photons are expected to be generated in the synchrotron process by the high energy tail of the particle distribution while the HE and VHE energy emission results from the up-scattering of seed photons in the EC part of the spectrum (see Section 1; Abramowski et al. 2012;Petropoulou 2014). Note, on the other hand, that due to the particularly sparse sampling of the XRT monitoring data, the X-ray power spectrum of PKS 2155−304 could be constrained only down to the variability timescales of ∼ 100 days, while in the multi-wavelength comparison Figure 5(b), we see that up to this timescale the variability amplitudes at X-rays are comparable to those seen at optical, and also γ-ray, photon energies; that is to say, with the available X-ray monitoring data for the blazar, we do not probe the variability timescales at which the difference between γ-ray (inverse-Compton) and optical (synchrotron) power spectra are significant.
The other important result comes from a systematic comparison of the squared fractional variability at the selected synchrotron and IC light curves. On timescales of months ( 100 days, down to the sampling intervals of the γ-ray light curves, i.e. ∼weeks), the observed variability amplitudes in the γ-ray range are significantly larger -up to one order of magnitude -than those observed at radio, IR or optical frequencies. For the blazar PKS 2155−304, quite surprisingly, within this exact time domain, the variability amplitudes seem larger at HE γ-rays than in the VHE γray range (1σ error bars), although it should be stressed out here once again that the sampling of the light curves within these three bands is very different (half-a-year observing window for the H.E.S.S. instruments versus a continuous monitoring with the Fermi-LAT). This behaviour is not shown by Mrk 421 where the variability amplitudes at X-rays (RXTE and Swift), HE and VHE γ-ray range seem to follow closely each other on timescales ≤100 days. Whether this tentative finding contrasts with the variability behavior of other blazars is to be investigated by means of a systematic analysis of multiwavelength light curves that covers many cycles of "quiescence" and "flaring" states for a larger sample of sources. At this point we only mention that in the framework of the leptonic scenarios for BL Lacs, the GeV emission results from IC scattering of seed photons that proceeds in mainly in the Thompson regime, while in the case of the TeV emission the Klein-Nishina effects become typically more relevant, depending on the exact shape of the source broad-band spectrum, and this may be one of the factors shaping the amplitudes of the observed flux changes (see in this context, e.g., Aleksić et al. 2015). All in all, the analysis results for the HBLs Mrk 421 and PKS 2155−304 presented in this paper, which are in accord with the results for the LBLs PKS 0735+178 and OJ 287 presented in Goyal et al. (2017Goyal et al. ( , 2018, strengthen our main conclusion, namely that commonly employed onezone scenarios for the blazar emission cannot account for all the complexity of the blazar variability observed on various timescales. Indeed, in the framework of the simplest onezone SSC model, different statistical patterns (in particular, red vs. pink noise-type PSDs) are not expected at different frequencies within the synchrotron-dominated spectral region (IR/optical vs. X-rays), neither should they appear between the synchrotron and IC-dominated spectral regions (IR/optical vs. γ−rays; see Finke & Becker 2014). This, in our opinion, points out to a more complex picture of blazars, with a highly inhomogeneous outflow producing non-thermal emission over an extended, stratified volume, as in this way stochastic, uncorrelated fluctuations in the local jet conditions, leading to the locally enhanced energy dissipation, could result in a correlated (colored noise-type) variability of the overall emission continuum (see in this context the discussion in Kelly et al. 2009, 2011, andreferences therein). In other words, while the one-zone emission scenarios could still be applied in a meaningful way when dealing with welldefined, and possibly shorter flaring events, they should be considered with caution when modelling the blazar emission continua averaged over longer monitoring epochs, or when comparing the blazar flux changes on vastly different timescales (say, years vs. hours and minutes).
Finally, we note that Mrk 421 and PKS 2155−304 are the only blazars for which broad-band variability power spectra could be generated right up to TeV γ−ray energies, without the VHE measurements being biased towards high flux states. The present work is also relevant to future studies of blazar variability up to VHE using the upcoming Cherenkov Telescope Array (CTA), which is expected to raise the number of currently known TeV-emitting blazars by an order of magnitude (Cherenkov Telescope Array Consortium et al. 2017). This paper has made use of up-to-date SMARTS optical/near-infrared light curves that are available at www.astro.yale.edu/smarts/glast/home.php.
The Fermi LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariatà l'Energie Atomique and the Centre National de la Recherche Scientifique / Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden.
Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France. This work performed in part under DOE Contract DE-AC02-76SF00515.
The support of the Namibian authorities and of the University of Namibia in facilitating the construction and operation of H.E.S.S. is gratefully acknowledged, as is the support by the German Ministry for Education and Research (BMBF), the Max Planck Society, the French Ministry for Research, the CNRS-IN2P3 and the Astroparticle Interdisciplinary Programme of the CNRS, the U.K. Particle Physics and Astronomy Research Council (PPARC), the IPNP of the Charles University, the South African Department of Science and Technology and National Research Foundation, and by the University of Namibia. We appreciate the excellent work of the technical support staff in Berlin, Durham, Hamburg, Heidelberg, Palaiseau, Paris, Saclay, and in Namibia in the construction and operation of the equipment.
This research has made use of data from the OVRO 40-m monitoring program (Richards, J. L. et al. 2011, ApJS, 194, 29) which is supported in part by NASA grants NNX08AW31G, NNX11A043G, and NNX14AQ89G and NSF grants AST-0808050 and AST-1109911 The spectral window function is a diagnostic of deleterious effects of Fourier transform when subjected to unevenly sampled time series. It is defined to be the Fourier transform of the sampling times, which in case of evenly sampled time series, is 1 at zeroth Fourier frequency and 0 otherwise and moreover is unknown for unevenly sampled data series (Deeming 1975). Figs A1 and A2 show the distributions of sampling interval of the analysed light curves and the spectral window functions of the respective data series (panels a-e and f-j, respectively for Mrk 421 and panels a-j and k-t, respectively, for PKS 2155−304). As shown, the distribution of sampling intervals is neither sharply peaked and concentrated in one bin (as is the case for evenly sampled data series) nor uniform over the intervals (as for randomly sampled data series). For the Fermi-LAT data, where the fluxes were measured using one-week integration times, the distribution is at intervals corresponding to integer multiples of seven days, as expected. The spectral window function for the data series is computed down to typical (mean) Nyquist sampling frequencies (Table 1). In each case, the spectral window function shows considerable power (other than 0) at all the Fourier frequencies precisely due to the uneven sampling of the light curves which could be mistaken for real power in the periodograms. Moreover, a mild peak at frequencies corresponding to ∼200-300 days is noted in the spectral window functions of the VERITAS, H.E.S.S., optical and IR light curves. This peak is due to systematic (seasonal) gaps attained in the light curves for the source from the ground-based observatories when the target is too close to the Sun. A more prominent peak at ∼28 days are noted in the spectral window functions of the VERITAS and the H.E.S.S. light curves because of systematics gaps due to full Moon periods during which the above facilities do not operate. This peak is more prominent than the ∼300 day peak because the data series contains more periods than the yearly gaps in the decade-long light curves. The spectral window function for the RXTE-PCA and the Swift-XRT light curve also show a mild peak around ∼300 days, perhaps because of scheduling of these particular observations. The spectral window function for the Fermi-LAT light curve shows a fairly uniform response, indicative of no systematic gaps. Therefore, a careful examination of spectral window function is essential as it can reveal periodicities that occur precisely due to regular gaps in the sampling of the time series which could be mistaken as a QSOs. We note that aliasing folds back roughly equal amount of power to the frequency range probed which is dominated at the Nyquist frequency (Uttley et al. 2002), thus could not create peaks in the spectral window function.

APPENDIX B: PROBABILITY DISTRIBUTION CURVES FOR THE ANALYSED PSDS
The probability that a single power-law fit with a given simulation (input) β describe the PSDs within the β range 0.2-3.0 is shown in Figs. B1 and B2 for the blazars Mrk 421 and PKS 2155−304, respectively. We note that p β is higher than 10 percent in all the cases except for OVRO PSD for the blazar Mrk 421 for which it is 8.4 percent. In general, the high p β values for the periodograms indicate that the single power-law fits describe the PSDs well. For clarity, we have truncated the histograms above the sampling intervals with only two or fewer data points. We note that the sampling interval extends up to 210 days, 21 days, 620 days, 288 days, and 324 days for the VERTIAS, Fermi-LAT, RXTE-PCA, Swift-XRT and the OVRO light curves, respectively (Table 1). (f-j): Corresponding spectral window functions.   Figure B1. Probablity, p β , that a single power-law slope is acceptable to describe the PSDs for the blazar Mrk 421. The filled circles show the fraction for which χ 2 dist,i > χ 2 obs was found out of 1,000 simulations for a fixed β. The solid lines show the Gaussian fits to the distributions.  Figure B2. Probablity, p β , that a single power-law slope is acceptable to describe the PSDs for the blazar PKS 2155−304. The filled circles show the fraction for which χ 2 dist,i > χ 2 obs was found out of 1,000 simulations for a fixed β. The solid lines show the Gaussian fits to the distributions.