Open Access
Issue
A&A
Volume 634, February 2020
Article Number A120
Number of page(s) 11
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/201935117
Published online 20 February 2020

© F. Ait Benkhali et al. 2020

Licence Creative Commons
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Open Access funding provided by Max Planck Society.

1. Introduction

Blazars belong to the most luminous and variable extragalactic sources in the Universe. They represent a special sub-class of active galactic nuclei (AGNs) characterized by a relativistic plasma jet oriented very close to the line of sight at angles θ  ≤  1/Γ, where Γ is the bulk Lorentz factor (e.g., Urry & Padovani 1995; Ackermann et al. 2015a). Blazar sources are known to be variable across all wavelengths from radio frequencies to teraelectronvolt γ-ray energies, and on a wide range of timescales from sub-minute to several years (e.g., Böttcher & Chiang 2002; Aharonian et al. 2007). Their multi-wavelength spectral energy distributions (SED) often exhibit two bumps. The first, low-energy bump (peaking at infrared to X-ray frequencies) is usually interpreted as synchrotron emission from highly relativistic electrons within the jet. The second bump on the other hand, peaking in the HE range, has frequently been related to inverse Compton emission (IC) up-scattering of various soft photon fields (mostly synchrotron or external thermal radiation), although hadronic interactions may contribute as well.

Periodic variability in the light curves of blazars has been investigated extensively in the radio and optical bands (e.g., Fan et al. 2002; Kadler et al. 2006; King et al. 2013; Bhatta et al. 2016). In this context, the two prominent sources OJ 287 and 3C 279 are worth mentioning, with long-term optical periods of ∼12 and 29.6 years being reported, respectively (Valtonen et al. 2006; Li et al. 2009). In the high-energy gamma-ray regime only five blazars (PG 1553+113, PKS 2155−304, PKS 0301−243, PKS 0426−380 and PKS 0537−441) have been reported so far to exhibit long-term quasi-periodic variability (these blazars are therefore referred to as QPVs) in their γ-ray fluxes with significance higher than 3σ (Ackermann et al. 2015b; Sandrinelli et al. 2014, 2016; Zhang et al. 2017a,b; Prokhorov & Moraghan 2017). Apart from PKS 0537−441, these QPVs have been seen over the whole length of the light curve. In the majority of these cases, attempts have been made to relate the detected periodic behavior at different wavelengths and on various timescales to the motion of two supermassive black holes in a binary system, general helical jet structures, shocks, or instabilities of the disk or jet-plasma flow (e.g., Rieger 2004; Wiita 2011; Sobacchi et al. 2017; Caproni et al. 2017; Tavani et al. 2018; Yan et al. 2018).

In the present paper we use Fermi-LAT γ-ray data between 2008 August 4 and 2018 March 15 to re-evaluate the five blazars mentioned above and to investigate an additional source, PKS 0449−439. Results of a periodicity search based on the generalized Lomb-Scargle periodogram (GLSP) and the Weighted Wavelet Z-transform (WWZ) are presented. The signifiance of the inferred periods is then estimated on the basis of Bootstrap resampling and light-curve simulations using the Timmer & Koenig (1995) and Emmanoulopoulos algorithms (Timmer & Koenig 1995; Emmanoulopoulos et al. 2013)

The paper is organized as follows. Section 2 provides an exemplary illustration of the performed Fermi-LAT data analysis procedures to generate γ-ray light curves based on a binned likelihood method and also presents the aperture photometry technique. In Sect. 3 we describe the GLSP and WWZ techniques used here to search for periodicity. Significance estimation and variability characterization is presented in Sect. 4. Finally, Sect. 5 provides a summary and discussion of the results.

2. Observations and analysis

2.1. Fermi-LAT likelihood analysis

Fermi-LAT on board the Fermi satellite is an electron-positron pair-conversion γ-ray detector sensitive to photon in the energy range between ∼20 MeV and 500 GeV. The LAT has a wide field of view of ∼2.4 sr and observes the entire sky every two orbits. It has a point spread function (PSF) < 0.8° and the largest effective area of approximately 8000 cm2 above 1 GeV. Fermi-LAT data were downloaded from the Fermi-LAT Data Server1 and events were selected between 2008 August 4 and 2018 March 15 (MET from 239557418.0 to 542821500.0), covering about ∼9.5 years with energy range from 100 MeV to 500 GeV in a circular region of interest (ROI) of 15° centered on the position of each source (Atwood et al. 2013).

A binned maximum likelihood analysis (Mattox et al. 1996) is performed using the Pass 8 (P8R2) algorithms and employing the standard Fermi Science Tools v10r0p5 software package. We follow the standard procedure provided by the Fermi Science Support Center (FSSC) to reduce the data. In addition, photons coming from zenith angles larger than 90° were all rejected to reduce the background from gamma rays produced in the atmosphere of the Earth (albedo) with P8R2_SOURCE_V6 instrument response functions (IRFs). We performed standard quality cuts in accordance with the Pass8 data analysis criteria. The background emission was modeled using the Galactic and isotropic diffuse emission gll_iem_v06_iso and P8R2_SOURCE_V6_v06 files two. All sources from the third LAT source catalog (3FGL; Acero et al. 2015) within the ROI are included in the model to ensure a satisfactory background modeling.

A binned-likelihood analysis with a bin size of 0.1° was performed. The spectral indices were allowed to vary for sources located within a radius of 5° around the position of the investigated source. The Test Statistic value (TS) defined as TS = −2ln(L0/L1) was used to determine the source detection significance, with the threshold set to TS = 25 (∼5σ). The significance of a source detection is given by (Abdo et al. 2010). In the case of PKS 0447−439 for example, the likelihood analysis reveals a point source with a high statistical significance TS ≃ 29305 corresponding to 171σ. Both energy and temporal bins with TS < 9 (or ∼3σ) are set as upper limits throughout the paper.

The appearance of new sources within the ROI could in principle strongly influence the γ-ray spectrum and light curve of the investigated source. In order to investigate this we also searched for possible new γ-ray sources within a FOV of 15°. Our analysis of 9.5 years of Fermi-LAT data indicates additional transient sources beyond the 3FGL catalog; these were appropriately taken into account in the full analysis, i.e., for each additional source we sequentially add a new point source with a standard spectral definition (Power law) and maximize the likelihood as a function of its flux.

We produced the γ-ray SED and light curves of each source through the binned maximum likelihood-fitting technique, with gtlike to determine the flux and TS value for each time bin. The effect of energy dispersion below 300 MeV is also accounted for in the analysis; to do so, we enabled the energy dispersion correction.

As an example, the SED of PKS 0447−439 is shown in Fig. 1 for the full time and energy range using ten energy bins. We performed a maximum likelihood analyses, exploring different spectral forms for the whole energy range, namely single power-law (PL), broken power-law (BPL), log-parabola (LP), and power-law super exponential cutoff (PLSEC). In the case of PKS 0447−439, a likelihood ratio test comparison yields TS = −2log(LPL/LBPL)≃53.8, thus preferring BPL over PL at a significance level of 7.3σ. The best fit indices are Γ1 = 1.84 ± 0.01 and Γ2 = 2.82 ± 0.17 with a break at an energy Eb = (42.8 ± 1.8) GeV. The obtained results of the maximum likelihood analysis are summarized for the different models in the Table 1.

thumbnail Fig. 1.

SED of PKS 0447−439 in the energy band of 100 MeV to 500 GeV as extracted from the complete data set (2008−2018) along with a Broken Power Law fit (solid red line). The regions show the 1σ confidence intervals resulting from the fit with BPL. The resulting break energy is Eb = 42.84 ± 1.78.

Open with DEXTER

Table 1.

Parameters obtained for the fit of the Fermi-LAT energy spectrum of PKS 0447−439 using power law (PL), log parabola (LP), broken power law (BPL) and PLSuperExpCutOff (PLSEC) spectral models.

2.2. Fermi-LAT γ-ray light curves

We generate γ-ray light curves for all our six sources using the individual best-fit results (i.e., a BPL model in the case of PKS 0447−439); see Table 2 for details. For comparison, two methods are employed: the maximum likelihood optimization and the aperture photometry. In the analysis with the first method the γ-ray light curves are produced by using separate, equal time bins of one month (30 days). We apply a maximum likelihood fitting technique by running the ScienceTools gtlike to extract the flux and TS value for each time bin. In this computationally more intensive procedure, the γ-events are selected in a circular ROI of 15° radius centered on the position of the source. The resulting light curve for PKS 0447−439 is shown in Fig. 2 and has an average γ-ray flux of ⟨ϕ⟩=(7.62 ± 0.13)×10−8 cm2 s−1. The source exhibits flux variability during the whole observational period, with flux levels of peak-to-peak oscillations changing by a factor four. The calculated fractional variability index is Fvar = (47.91 ± 1.16)% (Vaughan et al. 2003; Edelson et al. 2002).

thumbnail Fig. 2.

γ-ray light curve of PKS 0447−439 from August 2008 to December 2017 in the energy range between 100 MeV and 500 GeV in intervals of 30 days. The dotted line gives a fit with a constant flux. There is clear evidence for variability. Vertical error bars indicate the one-sigma error bars. The gray shaded regions show the 1σ, 3σ, and 5σ confidence intervals resulting from the fit with a constant function.

Open with DEXTER

Table 2.

Spectral parameters obtained for the best fit of Fermi-LAT data for each source from the energy range between 100 MeV and 500 GeV using four different spectral models: PL, BPL, LP and PLSEC.

The aperture photometry method on the other hand provides a model-independent measure of the γ-ray flux and is less computationally demanding. It also enables the use of short time bins whereas the maximum likelihood technique requires that time bins contain sufficient photons for the analysis. We select a very small aperture radius of 1° to exclude most background γ-events and to focus on the events that are most likely associated with the source itself. We produce γ-ray light curves in the range between 100 MeV and 500 GeV using a weekly binning. The exposure of each time bin is determined with the ScienceTools gtexposure. Figure 3 shows the resultant γ-ray light curve for PKS 0447−439 obtained with this method.

thumbnail Fig. 3.

γ-ray light curve of PKS 0447−439 from 2008 to 2018 in the energy range (0.1 − 500) GeV obtained using the aperture photometry technique for a radius of 1° around the source. A binning of one week is employed. Error bars correspond to 1σ error bars. The gray shaded regions show the 1σ, 3σ, and 5σ confidence intervals resulting from the fit with a constant function.

Open with DEXTER

In order to verify that the aperture photometry (AP) approach gives results which are consistent with the results obtained using the binned likelihood (BL) analysis, we performed diagnostic tests on each of the investigated blazars. The BL method is the preferred procedure for most types of Fermi-LAT analysis. Taking for granted that BL flux values are indeed more accurately measured than the AP flux, we examine the correlation between the AP and BL calculated flux to assess the performance of the AP method. This correlation is plotted for PKS 0447−439 in Fig. 4 using monthly binned light curves and the calculated value of the Pearson correlation coefficient is ρ ∼ 0.936.

thumbnail Fig. 4.

Aperture photometry flux in comparison with the binned likelihood flux calculated for PKS 0447−439 (logarithmic scale). The red solid line represents a fit of the data with a linear function and the gray shaded regions show the resulting 3σ and 5σ confidence intervals. Error bars indicate the 1σ error interval. No background subtraction is used in the aperture photometry method.

Open with DEXTER

2.3. Gaussian versus log-normal flux distributions

The almost continuous detection of the investigated sources over time bins of seven days provides a good opportunity to investigate whether the probability density function (PDF) for the long-term γ-ray emission reveals any preference for a Gaussian (normal) or a log-normal flux distribution. It is expected that such a feature offers important clues for understanding the central engine of a blazer (e.g., Uttley et al. 2005; Shah et al. 2018; Rieger 2019). A log-normal distribution of γ-ray fluxes, for example, would suggest that the mechanism driving the variability is multiplicative in nature rather than additive. Such multiplicative behavior could possibly be related to accretion disk fluctuations (Lyubarskii 1997; Arévalo & Uttley 2006) or the particle acceleration process itself (Sinha et al. 2018). As a possible caveat we note that the PDF of stochastic processes with power spectral density (PSD) steeper than index −1 can show deviations from normality owing to divergence of power at low frequencies (e.g., Vaughan et al. 2003).

We explore the PDF and quantify this in terms of fluxes using flux-histograms. For each source the weekly γ-ray fluxes were distributed in a histogram of fluxes. We fit all histograms in log-scale, with Gaussian G(ϕ) and log-normal L(ϕ) distribution functions given by

(1)

and

(2)

respectively, where σ and μ are the standard deviation and the mean of the distribution, respectively. For illustration, the obtained flux histograms for PKS 0447−439, PG 1553+113, and PKS 2155−304 are shown in Fig. 5. We compute the Anderson Darling test (AD) statistics (Anderson & Darling 1954) for each of the light curves as coming from a Gaussian or a log-normal distribution. The AD allows one to determine whether one can accept or reject the null hypothesis that a sample is drawn from a population that follows a particular distribution; in our case a normal distribution. The results obtained are shown in Table 3. The obtained (AD test) values provide strong evidence in favor of rejecting the null hypothesis of an underlying normal distribution as the AD test value is greater than the relevant critical value (e.g., Stephens 1974).

thumbnail Fig. 5.

Normalized histograms of γ-ray photon fluxes fitted with a log-normal (solid red line) and Gaussian (blue solid line), respectively. Apart from PG 1553+113, all sources show a clear preference for a log-normal distribution.

Open with DEXTER

Table 3.

Best fit parameters for the log-normal and normal (Gaussian) distribution, respectively; the AD test statistics and are the Anderson Darling test and the reduced χ2, respectively, [*] in units of photon cm−2 s−1.

Except for PG 1553+113, all sources exhibit a clear deviation from a Gaussian distribution. The flux distributions are instead compatible with being log-normal, which suggests that the underlying variability is of a nonlinear, multiplicative origin. Similar results are obtained from the D’Agostino’s K2 and the Shapiro–Wilk test; see Table 4.

Table 4.

K2 Test of D’Agostino tests the null hypothesis that a sample comes from a normal distribution.

3. Searching for periodicity

In order to search for periodicity in blazars we apply two widely used methods to light curves: the Lomb-Scargle periodogram (LSP) and the WWZ.

3.1. Lomb-Scargle periodogram

The LSP is a commonly used algorithm for detecting and characterizing periodicity in unevenly sampled light curves (Lomb 1976; Scargle 1982). The standard normalized LSP is equivalent to fitting sine waves of the form y(t) = A cos(ωt)+B sin(ωt), and is defined for an uneven, simple time series (ti, yi ) as

(3)

The standard LSP however has some limitations. On the one hand, it does not take measurement errors into account; on the other hand, in the calculation, the mean of the flux is subtracted, which assumes that the mean of the data and the mean of the fitted sine function are the same. To overcome those limitations we instead use the GLSP (Zechmeister & Kürster 2009). As an example, the calculated GLSP for PKS 0447−439 using monthly γ-ray light curves (generated based on the binned likelihood) is shown in Fig. 8. A strong peak at around a period of 945 ± 40 days is apparent in this figure. We have estimated the uncertainty of the calculated period based on the half-widths at half-maxima (HWHM) of Gaussian fits to the profile at the position of the highest peak. An evaluation of confidence levels determined from Monte Carlo simulations of colored noise is explained in detail in Sect. 4.

For comparison, the GLSP method has also been applied to the light curves based on aperture photometry. The possible interference of nearby sources (including background contribution) is evaluated for each source, as shown in Figs. 6 and 7 for the case of PKS 0447−439. The results agree well with those obtained based on the binned likelihood method, as evident from Table 5.

thumbnail Fig. 6.

Smoothed count map (logarithmic scale) of the PKS 0447−439 region in energy range from 100 MeV up to 500 GeV as seen by Fermi-LAT based on data between August 2008 and December 2017. The color bar has units of counts per pixel and the pixel dimensions are 0.1 × 0.1°. The locations of the background region and surrounding sources (labeled P1 = new source beyond the 3FGL catalog, P2 = 3FGL J0438.8−4519, P3 = 3FGL J0437.2−4713 and P4 = 3FGL J0455.7−4617) are shown as circles. The magenta circle shows the location of PKS 0447−439 and the yellow circles represent the background region.

Open with DEXTER

thumbnail Fig. 7.

Generalized LSPs obtained based on light curves generated using aperture photometry for PKS 0447−439 and surrounding sources (P1,P2,P3,P4) shown as circles in Fig. 6.

Open with DEXTER

Table 5.

Periodicity results derived from the GLSP method for the monthly and weekly γ-ray light curves, respectively, and including the intrinsic (redshift-corrected) period. [*] in units of days.

The observed period Pobs. for a blazar at redshift z is related to its intrinsic period by Pint. = Pobs./(1 + z). Table 5 shows the GLSP results for all sources, except PKS 0537−441, based on both monthly (generated based on the binned likelihood) and weekly (aperture photometry) binned light curves. Year-type HE periodicity is found for all five sources, with the results based on aperture photometry or binned likelihood method being in good agreement. For the blazar PKS 0537−441 no significant period was found from a search of the entire light curve. Furthermore, the source is estimated to have a red-noise power spectrum (Covino et al. 2019), which suggests a greater possibility of detecting a fake peak. This source is therefore not included in the following. We confirm however that during the initial approximately three-year high state, a peak can be seen on a periodicity timescale of ∼280 d (see e.g., Sandrinelli et al. 2016). However, a longer observation window showing a larger number of cycles is needed, given the timescale of periodicity (Vaughan et al. 2016), for this result to be confirmed.

3.2. Weighted wavelet Z-transform

The GLSP provides an excellent tool for the periodicity analysis of light curves with unevenly spaced sampling. Nevertheless it does not account for the possibility that in some astrophysical systems quasi-periodic oscillations may develop that vary significantly in frequency and amplitude over a specified period of time. In such cases, the WWZ method turns out to be a more convenient technique for detecting and quantifying such variations (Foster 1996; Han & van der Baan 2013), and has been applied in the timing analysis of AGN light curves at different wavelengths (e.g., Bhatta et al. 2016; Mohan & Mangalam 2015). The method is based on a similar concept like LSP, where sinusoidal functions are used to fit the data; however, the waves can now be localized in both time and frequency domain to account for the possible transient nature of the quasi-periodic oscillations (QPOs; Torrence & Compo 1998; Bravo et al. 2014).

The WWZ is based on a weighed projection onto three trial functions, y(t) = ∑iyiϕi(t),

(4)

where the “best-fit” coefficients yi are the coefficients for which the model function y(t) best fits the data. For each projection, the statistic weight is given by

(5)

where c is a constant that determines how rapidly the Morlet wavelet decays, and is usually chosen to be close to 0.0125. The WWZ power can then finally be defined as

(6)

with Vx and Vy as weighted variations of the data and the model function, respectively, and Neff representing the effective number of data points (for more details, see Foster 1996).

The WWZ powers for PKS 0447−439 are shown in Fig. 9 as a function of both observing time (horizontal axis) and period (vertical axis). The peaks in the power characterize the strength and duration of a possible quasi-periodic modulation in the data. The WWZ indicates a characteristic period at (937.3 ± 152.7) days (around 2.6 years) with peak power of about 18.5 and no significant changes over time. The constancy of the period indicates that the quasi-periodicity is probably driven by a physical process that is stable over the duration of the observation. The results of the WWZ method are shown in Table 6 and the results agree well with those from the GLSP.

thumbnail Fig. 8.

GLPS obtained using the monthly light curves of PKS 0447−439 is shown in cyan (circles); the red solid line, green line, and orange line represent the 3σ, 2.5σ, and 2σ confidence levels, respectively, calculated based on simulations of 50 000 light curves using the method of Timmer & Koenig (1995).

Open with DEXTER

thumbnail Fig. 9.

Two-dimensional WWZ of the γ-ray light curve of PKS 0447−439 based on Fermi-LAT data from August 2008 to December 2017. The color scale represents the Z-statistics of the WWZ power of a certain period at a given time. The panel shows the signature of a quasi-periodic oscillation at a (observed) period of ∼2.6 years without any significant changes over time.

Open with DEXTER

Table 6.

Summary results of the strongest periods for the WWZ method based on the monthly γ-ray light curves, including observed and intrinsic period.

4. Significance and uncertainty estimation

The effects related to irregular sampling of a light curve and the noisy nature of the periodogram can in some situations lead to the generation of false (artificial) periods in the GLSP that could be mistaken as a real periodic signal of the source. For this reason, it is important to take such effects into account when searching for periodicities and to analyze their impact on the GLSP carefully.

The variability of AGN light curves often exhibits a colored-noise-like behavior with PSD characterized by a simple power law of the form PSD(ν) ∼ νβ where ν is the (temporal) frequency and β the power-law index.

4.1. Power spectrum response method

In order to determine the appropriate form of the underlying colored noise needed as an input for simulating the light curves, we first applied the power spectrum response (PSRESP; Uttley et al. 2002) method, which is a widely used technique for the characterization of AGN power spectra (e.g., Chatterjee et al. 2008; Bhatta et al. 2016). The method attempts to fit the binned periodogram with different realisations of a given PSD model in order to determine which model maximizes the probability that the observed PSD can be reproduced. Our implementation of the PSRESP method is described in full detail in Chatterjee et al. (2008). Selected details are as follows.

We consider a simple power-law model for the underlying power spectrum of the form

(7)

where A is the amplitude of the model at the reference frequency ν0, β corresponds to the power-spectral slope, and Cnoise is a constant which is fixed at the Poisson noise level for the light curve. We simulate N = 1000 light curves starting from the underlying model and using the Timmer & Koenig algorithm (see below), and re-sample these with the same sampling interval of the observed light curves. We then calculate the periodogram of the observed light curve (P(ν)obs.) and that of each of the simulated light curves (P(ν)sim, i where i = 1, …N). The statistic is calculated from the underlying model average and observed PSD of each light curve, with

(8)

(9)

where νmin and νmax are the minimum and maximum frequencies measured by P(ν)obs., respectively. We calculate the for each simulated PSD P(ν)sim, i over all frequencies between νmin and νmax with respect to the data. We then compare with the distribution and count the number m of for which is smaller than . Finally, we calculate the success fraction m/M (goodness of fit), which gives the probabilities of a model being accepted, for a range of β from 0.0 to 2.0 with a step size of 0.05. The obtained results according to the PSRESP method are summarized in Table 7 and shown in Fig. 10.

thumbnail Fig. 10.

Success fraction vs. slope (β) for all three PSDs (Fermi-LAT). The success fractions indicate the goodness of fit obtained from the PSRESP method (see text). (a) Weekly, (b) monthly light curves.

Open with DEXTER

Table 7.

Slopes β obtained from fitting the γ-ray power spectra with a simple power-law model, and the corresponding success fraction calculated based on the PSRESP method (Uttley et al. 2002) (see Sect. 4.1) (the success fractions indicate the goodness of fit obtained from the PSRESP method).

4.2. MC-simulations of colored-noise light curves

The Timmer and Koenig (TK) algorithm (Timmer & Koenig 1995) is a commonly used method to produce artificial light curves. This technique allows for the generation of nondeterministic (stochastic Gaussian-distributed) time series from a given underlying PSD model by randomizing both the phase and the amplitude of the Fourier components. Limitations could arise however for light curves that exhibit strong deviations from Gaussian distributions (e.g., a burst-like behavior).

Given the preference for log-normality, we also simulate 5 × 104 light curves using the method proposed by Emmanoulopoulos et al. (2013, E13) to obtain best-fitting results from the PSRESP method. The latter E13 method is able to account for a general PDF, that is, it can match both the PSD and the PDF of an observed light curve, thus relaxing restrictions of the TK method. This does in fact better comply with our previous findings of nonGaussianity in Sect. 2.3.

The results are shown in Fig. 8 using the (main value of the) best fit slope β. In many cases the detected periods are close to or above 3σ. Given available data there is rather little difference between the results based on the Timmer & Koenig (1995) and Emmanoulopoulos et al. (2013) methods. The outcome is however obviously dependent on the slope β. Within the PSRESP inferred range quite different results can be obtained as shown in Table 9, suggesting that a narrowing-down of the PSD slope will be most relevant for assessing the real significance of the detected periods.

Table 8.

Significance values of the observed periods for the local and the global method using the Timmer & Koenig (1995) and Emmanoulopoulos et al. (2013) algorithm, based on the monthly and weekly γ-ray light curves.

Table 9.

Significance values of the observed periods using either the local or the global method with respect to the monthly and weekly γ-ray light curves of PG 1553+113 and PKS 0447−489, respectively.

4.3. Statistical confidence of the GLSP-detected period

4.3.1. Local significance

In the presence of an a priori expectation of periodicity either from theory or from entirely independent observations, one can compute the statistical significance of the detected period at that position. This is the so-called “local significance” as referenced in several papers (e.g., Bell et al. 2011). Therefore, using the local method, the frequency channels corresponding to the GLSP-detected period were sought and the power spectra in which the peak power at this period exceeded the observed value were recorded. Finally, the fraction of occurrences with power greater than the detected period is the probability of a false positive signal resulting from random red noise in the observed light curve.

4.3.2. Global significance

In the majority of the cases of reported periodicities, particularly in gamma-rays, we do not have strong a priori indications of the expected periodicity. In such a circumstance, it is statistically more rigorous to evaluate the “global” rather than the local significance (e.g., Bell et al. 2011). This constitutes computing the fraction of occurrences of larger peaks at any period within a reasonable range of timescales (dependent on the cadence properties of the light curve), relative to the detected one. This ensures that we factor in false positives from this larger range of timescales rather than a specific period. This so-called “look-elsewhere” effect (or also “multiple testing problem” in statistics) can be quantified in terms of a trial factor (e.g., Lyons 2008; Gross & Vitells 2010). This is defined as the ratio of the probability of detecting a peak (or period/excess) at some fixed frequency to the probability of detecting it anywhere in the (tested) range. The look-elsewhere effect has been explored and factored in for detection of resonant peaks in particle physics and indeed other fields including astronomy. We therefore re-evaluated the probability that the power of any observed peak is equal to or greater than a selected value somewhere in the periodogram and calculated the global significance (see e.g., Vaughan 2010). The results are shown in Table 8 and reveal that in the absence of other physical reasons for restricting the period range, the QPV evidence is strongly reduced, with none of the sources reaching 3σ significance. These findings are in line with similar indications in Covino et al. (2019).

5. Discussion and conclusions

There are an increasing number of reports suggesting the presence of periodicities of years in duration in the γ-ray light curves of Fermi-LAT-detected blazars. Given the still limited duration of the light curves (≲10 yr), care needs to be taken however to properly assess the significance of these periods against a colored-noise background. Vaughan et al. (2016) have shown, for example, that clear phantom periodicities can be found even in pure noise data, with typical periods corresponding to (1.5 − 2.5) cycles over the available data.

Our systematic investigation performed in this study reveals that four out of the six investigated blazars show long-term QPO indications (local significance ≳3σ) in their Fermi-LAT light curves with an intrinsic period around 1.6 yr. As we have shown, there is some uncertainty as to the significance of these periodicities. Depending on the inferred best-fit PSD slope, the (local) significance can encompass a range from ∼2.6σ to > 5σ (method 1); see Table 8. The QPO significance is however strongly reduced in terms of global significance, suggesting that claims of long-term periods should be treated with caution.

While all sources except PG 1553+113 show clear indications for log-normality in their distribution of fluxes, incorporation of the appropriate simulation method (Timmer & Koenig 1995 or Emmanoulopoulos et al. 2013) does not have a strong impact on the significance evaluation. Improving the PSD characterization, that is, by narrowing down the range for the PSD slope, will instead be more relevant to better assess the significance.

Though the inferred periods are only tentative, one could speculate over the possible physical mechanisms capable of causing these longer periodicities. It is interesting, for example, that the intrinsic periods appear quite similar, clustering around P ≈ 1.6 yr for sources at different redshifts (with the caveat of some possible redshift uncertainty for PKS 0447−439).

Perhaps one of the most natural scenarios for the origin of long-term periodic variability is related to the orbital motion in a supermassive binary black hole system (SBBHS). In principle, a SBBHS phase is likely to occur in radio-loud AGNs (hosted by elliptical galaxies) at some stage. However, the periods inferred here appear too short to plausibly relate them to orbital SBBH motion. Gravitation radiation would lead to coalescence on a characteristic timescale of only tgrav  ≃  7  ×  103 (1 + q)1/3/[q (MBH/5  ×  108M)5/3] (P/2 yr)8/3 yr, meaning that for typical systems with q := m/M ≥ 0.05 the presumed source state would be highly unlikely. In fact, one alternatively expects orbital periods of the order of ∼10 yr for SBBHSs in blazars (Rieger 2007). This would then also be compatible with pulsar timing constraints on the inferred gravitational background (Holgado et al. 2018). As orbital motion usually allows for the shortest periodic driving (e.g., Rieger 2004), a direct SBBH origin of the observed periods appears less likely. This does not argue against the presence of SBBHSs in blazars in general, but simply suggests that caution should be taken when directly relating periodicities of the order of P ∼ 1 yr to such systems.

Alternatively, QPOs might be related to quasi-periodic changes in accretion flow conditions that are effectively transmitted to the jet, modulating its nonthermal emission properties. Time-dependent modulations of the transition radius rt between an outer cooling-dominated (standard) disk and an inner radiatively inefficient flow (ADAF), for example, could lead to periodic mass flux variations (e.g., Gracia et al. 2003). If one requires the advective timescale tad(rt)∼rg (0.5αc)−1(rt/rg)3/2 to be (at most) comparable to P, with α = 0.25 being the viscosity coefficient and rg = GMBH/c2 the gravitational radius, this would place the transition radius at a characteristic scale of rt ≲ 200 rg(P/1.6 yr)2/3 (5 × 108M/MBH)2/3 for a reference mass of MBH = 5 × 108M. This seems compatible with estimates for the transition radius in BL Lacs (e.g., Cao 2003; Xie et al. 2008). This would then suggest a similar black hole mass-scale for the systems investigated here.

On the other hand, year-type QPVs could perhaps also trace plasma motion in the jet close to its outer jet radius r0. In the lighthouse model (Camenzind & Krockenberger 1992) for example, the disk-related jet is initially rotating, leading to a helical trajectory for a component injected on scales r0 ∼ 10 rL beyond the light cylinder rL ∼ 10 rg of the innermost part of the disk magnetosphere. Angular momentum conservation would imply a characteristic intrinsic period for such a component of P = 2πrL(r0/rL)2/c ≲ 2 (r0/20rL)2(MBH/5 × 108M) yr. While collimation might occur earlier (e.g., Fendt 1997), thereby reducing the jet radius, slightly changing footpoint radii could possibly compensate for this. The lighthouse model was originally designed to account for observed QPOs with Pobs of less than or equal to a few weeks by taking into account travel time effects with respect to an outwardly moving, (single) flaring component. It seems likely however that the fundamental period P might be visible even if the flux were quickly suppressed. It will be interesting to probe this with a larger source sample, as intrinsic periods in this case are expected to be less than 2 yr for a typical black hole mass range.

The fact that the intrinsic periods seem to be clustering around P ∼ 1.6 yr remains particularly interesting and suggestive of a common physical mechanism. The inferred periods are however only tentative. Adding one or two more cycles of data (i.e., ∼2 − 3 yr) is expected to significantly improve the situation and to help to clarify their putative presence and physical implication.


Acknowledgments

Useful comments by an anonymous referee are gratefully acknowledged. We would also like to thank J. Graham, A. Sinha, and S. Wagner (LSW Heidelberg) for helpful comments and suggestions. FMR acknowledges support by a DFG Heisenberg Fellowship RI 1187/6-1.

References

  1. Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJS, 187, 460 [NASA ADS] [CrossRef] [Google Scholar]
  2. Acero, F., Ackermann, M., Ajello, M., et al. 2015, ApJS, 218, 23 [NASA ADS] [CrossRef] [Google Scholar]
  3. Ackermann, M., Ajello, M., Atwood, W. B., et al. 2015a, ApJ, 810, 14 [NASA ADS] [CrossRef] [Google Scholar]
  4. Ackermann, M., Ajello, M., Albert, A., et al. 2015b, ApJ, 813, L41 [NASA ADS] [CrossRef] [Google Scholar]
  5. Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2007, ApJ, 664, L71 [NASA ADS] [CrossRef] [Google Scholar]
  6. Anderson, T. W., & Darling, D. A. 1954, J. Am. Stat. Assoc., 49, 765 [CrossRef] [Google Scholar]
  7. Arévalo, P., & Uttley, P. 2006, MNRAS, 367, 801 [NASA ADS] [CrossRef] [Google Scholar]
  8. Atwood, W., Albert, A., Baldini, L., et al. 2013, 2012 Fermi Symposium Proceedings – eConf C121028 [Google Scholar]
  9. Bell, M. E., Tzioumis, T., Uttley, P., et al. 2011, MNRAS, 411, 402 [NASA ADS] [CrossRef] [Google Scholar]
  10. Bhatta, G., Zola, S., Stawarz, Ł., et al. 2016, ApJ, 832, 47 [NASA ADS] [CrossRef] [Google Scholar]
  11. Böttcher, M., & Chiang, J. 2002, ApJ, 581, 127 [NASA ADS] [CrossRef] [Google Scholar]
  12. Bravo, J. P., Roque, S., Estrela, R., Leão, I. C., & De Medeiros, J. R. 2014, A&A, 568, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Camenzind, M., & Krockenberger, M. 1992, A&A, 255, 59 [NASA ADS] [Google Scholar]
  14. Cao, X. 2003, ApJ, 599, 147 [NASA ADS] [CrossRef] [Google Scholar]
  15. Caproni, A., Abraham, Z., Motter, J. C., & Monteiro, H. 2017, ApJ, 851, L39 [NASA ADS] [CrossRef] [Google Scholar]
  16. Chatterjee, R., Jorstad, S. G., Marscher, A. P., et al. 2008, ApJ, 689, 79 [NASA ADS] [CrossRef] [Google Scholar]
  17. Covino, S., Sandrinelli, A., & Treves, A. 2019, MNRAS, 482, 1270 [NASA ADS] [CrossRef] [Google Scholar]
  18. D’Agostino, R. B. 1970, Biometrika, 57, 679 [Google Scholar]
  19. D’Agostino, R., & Pearson, E. S. 1973, Biometrika, 60, 613 [Google Scholar]
  20. Edelson, R., Turner, T. J., Pounds, K., et al. 2002, ApJ, 568, 610 [NASA ADS] [CrossRef] [Google Scholar]
  21. Emmanoulopoulos, D., McHardy, I. M., & Papadakis, I. E. 2013, MNRAS, 433, 907 [NASA ADS] [CrossRef] [Google Scholar]
  22. Fan, J. H., Lin, R. G., Xie, G. Z., et al. 2002, A&A, 381, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Fendt, C. 1997, A&A, 319, 1025 [NASA ADS] [Google Scholar]
  24. Foster, G. 1996, AJ, 112, 1709 [NASA ADS] [CrossRef] [Google Scholar]
  25. Gracia, J., Peitz, J., Keller, C., & Camenzind, M. 2003, MNRAS, 344, 468 [NASA ADS] [CrossRef] [Google Scholar]
  26. Gross, E., & Vitells, O. 2010, Eur. Phys. J. C, 70, 525 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Han, J., & van der Baan, M. 2013, Geophysics, 78, O9 [CrossRef] [Google Scholar]
  28. Holgado, A. M., Sesana, A., Sandrinelli, A., et al. 2018, MNRAS, 481, L74 [NASA ADS] [CrossRef] [Google Scholar]
  29. Kadler, M., Hughes, P. A., Ros, E., Aller, M. F., & Aller, H. D. 2006, A&A, 456, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. King, O. G., Hovatta, T., Max-Moerbeck, W., et al. 2013, MNRAS, 436, L114 [NASA ADS] [CrossRef] [Google Scholar]
  31. Li, H. Z., Xie, G. Z., Chen, L. E., et al. 2009, PASP, 121, 1172 [NASA ADS] [CrossRef] [Google Scholar]
  32. Lomb, N. R. 1976, Ap&SS, 39, 447 [NASA ADS] [CrossRef] [Google Scholar]
  33. Lyons, L. 2008, Ann. Appl. Stat., 2, 887 [CrossRef] [MathSciNet] [Google Scholar]
  34. Lyubarskii, Y. E. 1997, MNRAS, 292, 679 [NASA ADS] [CrossRef] [Google Scholar]
  35. Mattox, J. R., Bertsch, D. L., Chiang, J., et al. 1996, ApJ, 461, 396 [NASA ADS] [CrossRef] [Google Scholar]
  36. Mohan, P., & Mangalam, A. 2015, ApJ, 805, 91 [NASA ADS] [CrossRef] [Google Scholar]
  37. Muriel, H., Donzelli, C., Rovero, A. C., & Pichel, A. 2015, A&A, 574, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Prokhorov, D. A., & Moraghan, A. 2017, MNRAS, 471, 3036 [NASA ADS] [CrossRef] [Google Scholar]
  39. Rieger, F. M. 2004, ApJ, 615, L5 [NASA ADS] [CrossRef] [Google Scholar]
  40. Rieger, F. M. 2007, Ap&SS, 309, 271 [NASA ADS] [CrossRef] [Google Scholar]
  41. Rieger, F. 2019, Galaxies, 7, 28 [NASA ADS] [CrossRef] [Google Scholar]
  42. Sandrinelli, A., Covino, S., & Treves, A. 2014, ApJ, 793, L1 [NASA ADS] [CrossRef] [Google Scholar]
  43. Sandrinelli, A., Covino, S., & Treves, A. 2016, ApJ, 820, 20 [NASA ADS] [CrossRef] [Google Scholar]
  44. Scargle, J. D. 1982, ApJ, 263, 835 [NASA ADS] [CrossRef] [Google Scholar]
  45. Shah, Z., Mankuzhiyil, N., Sinha, A., et al. 2018, Res. Astron. Astrophys., 18, 141 [NASA ADS] [CrossRef] [Google Scholar]
  46. Shapiro, S. S., & Wilk, M. B. 1965, Biometrika, 52, 591 [Google Scholar]
  47. Sinha, A., Khatoon, R., Misra, R., et al. 2018, MNRAS, 480, L116 [NASA ADS] [CrossRef] [Google Scholar]
  48. Sobacchi, E., Sormani, M. C., & Stamerra, A. 2017, MNRAS, 465, 161 [NASA ADS] [CrossRef] [Google Scholar]
  49. Stephens, M. A. 1974, J. Am. Stat. Assoc., 69, 730 [CrossRef] [Google Scholar]
  50. Tavani, M., Cavaliere, A., Munar-Adrover, P., & Argan, A. 2018, ApJ, 854, 11 [NASA ADS] [CrossRef] [Google Scholar]
  51. Timmer, J., & Koenig, M. 1995, A&A, 300, 707 [NASA ADS] [Google Scholar]
  52. Torrence, C., & Compo, G. P. 1998, Bull. Am. Meteorol. Soc., 79, 61 [NASA ADS] [CrossRef] [Google Scholar]
  53. Urry, C. M., & Padovani, P. 1995, PASP, 107, 803 [NASA ADS] [CrossRef] [Google Scholar]
  54. Uttley, P., McHardy, I. M., & Papadakis, I. E. 2002, MNRAS, 332, 231 [NASA ADS] [CrossRef] [Google Scholar]
  55. Uttley, P., McHardy, I. M., & Vaughan, S. 2005, MNRAS, 359, 345 [NASA ADS] [CrossRef] [Google Scholar]
  56. Valtonen, M. J., Lehto, H. J., Sillanpää, A., et al. 2006, ApJ, 646, 36 [NASA ADS] [CrossRef] [Google Scholar]
  57. Vaughan, S. 2010, MNRAS, 402, 307 [NASA ADS] [CrossRef] [Google Scholar]
  58. Vaughan, S., Edelson, R., Warwick, R. S., & Uttley, P. 2003, MNRAS, 345, 1271 [NASA ADS] [CrossRef] [Google Scholar]
  59. Vaughan, S., Uttley, P., Markowitz, A. G., et al. 2016, MNRAS, 461, 3145 [NASA ADS] [CrossRef] [Google Scholar]
  60. Wiita, P. J. 2011, J. Astrophys. Astron., 32, 147 [NASA ADS] [CrossRef] [Google Scholar]
  61. Xie, Z. H., Hao, J. M., Du, L. M., Zhang, X., & Jia, Z. L. 2008, PASP, 120, 477 [NASA ADS] [CrossRef] [Google Scholar]
  62. Yan, D., Zhou, J., Zhang, P., Zhu, Q., & Wang, J. 2018, ApJ, 867, 53 [NASA ADS] [CrossRef] [Google Scholar]
  63. Zechmeister, M., & Kürster, M. 2009, A&A, 496, 577 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  64. Zhang, P.-F., Yan, D.-H., Zhou, J.-N., et al. 2017a, ApJ, 845, 82 [NASA ADS] [CrossRef] [Google Scholar]
  65. Zhang, P., Yan, D., Liao, N., et al. 2017b, ApJ, 842, 10 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1.

Parameters obtained for the fit of the Fermi-LAT energy spectrum of PKS 0447−439 using power law (PL), log parabola (LP), broken power law (BPL) and PLSuperExpCutOff (PLSEC) spectral models.

Table 2.

Spectral parameters obtained for the best fit of Fermi-LAT data for each source from the energy range between 100 MeV and 500 GeV using four different spectral models: PL, BPL, LP and PLSEC.

Table 3.

Best fit parameters for the log-normal and normal (Gaussian) distribution, respectively; the AD test statistics and are the Anderson Darling test and the reduced χ2, respectively, [*] in units of photon cm−2 s−1.

Table 4.

K2 Test of D’Agostino tests the null hypothesis that a sample comes from a normal distribution.

Table 5.

Periodicity results derived from the GLSP method for the monthly and weekly γ-ray light curves, respectively, and including the intrinsic (redshift-corrected) period. [*] in units of days.

Table 6.

Summary results of the strongest periods for the WWZ method based on the monthly γ-ray light curves, including observed and intrinsic period.

Table 7.

Slopes β obtained from fitting the γ-ray power spectra with a simple power-law model, and the corresponding success fraction calculated based on the PSRESP method (Uttley et al. 2002) (see Sect. 4.1) (the success fractions indicate the goodness of fit obtained from the PSRESP method).

Table 8.

Significance values of the observed periods for the local and the global method using the Timmer & Koenig (1995) and Emmanoulopoulos et al. (2013) algorithm, based on the monthly and weekly γ-ray light curves.

Table 9.

Significance values of the observed periods using either the local or the global method with respect to the monthly and weekly γ-ray light curves of PG 1553+113 and PKS 0447−489, respectively.

All Figures

thumbnail Fig. 1.

SED of PKS 0447−439 in the energy band of 100 MeV to 500 GeV as extracted from the complete data set (2008−2018) along with a Broken Power Law fit (solid red line). The regions show the 1σ confidence intervals resulting from the fit with BPL. The resulting break energy is Eb = 42.84 ± 1.78.

Open with DEXTER
In the text
thumbnail Fig. 2.

γ-ray light curve of PKS 0447−439 from August 2008 to December 2017 in the energy range between 100 MeV and 500 GeV in intervals of 30 days. The dotted line gives a fit with a constant flux. There is clear evidence for variability. Vertical error bars indicate the one-sigma error bars. The gray shaded regions show the 1σ, 3σ, and 5σ confidence intervals resulting from the fit with a constant function.

Open with DEXTER
In the text
thumbnail Fig. 3.

γ-ray light curve of PKS 0447−439 from 2008 to 2018 in the energy range (0.1 − 500) GeV obtained using the aperture photometry technique for a radius of 1° around the source. A binning of one week is employed. Error bars correspond to 1σ error bars. The gray shaded regions show the 1σ, 3σ, and 5σ confidence intervals resulting from the fit with a constant function.

Open with DEXTER
In the text
thumbnail Fig. 4.

Aperture photometry flux in comparison with the binned likelihood flux calculated for PKS 0447−439 (logarithmic scale). The red solid line represents a fit of the data with a linear function and the gray shaded regions show the resulting 3σ and 5σ confidence intervals. Error bars indicate the 1σ error interval. No background subtraction is used in the aperture photometry method.

Open with DEXTER
In the text
thumbnail Fig. 5.

Normalized histograms of γ-ray photon fluxes fitted with a log-normal (solid red line) and Gaussian (blue solid line), respectively. Apart from PG 1553+113, all sources show a clear preference for a log-normal distribution.

Open with DEXTER
In the text
thumbnail Fig. 6.

Smoothed count map (logarithmic scale) of the PKS 0447−439 region in energy range from 100 MeV up to 500 GeV as seen by Fermi-LAT based on data between August 2008 and December 2017. The color bar has units of counts per pixel and the pixel dimensions are 0.1 × 0.1°. The locations of the background region and surrounding sources (labeled P1 = new source beyond the 3FGL catalog, P2 = 3FGL J0438.8−4519, P3 = 3FGL J0437.2−4713 and P4 = 3FGL J0455.7−4617) are shown as circles. The magenta circle shows the location of PKS 0447−439 and the yellow circles represent the background region.

Open with DEXTER
In the text
thumbnail Fig. 7.

Generalized LSPs obtained based on light curves generated using aperture photometry for PKS 0447−439 and surrounding sources (P1,P2,P3,P4) shown as circles in Fig. 6.

Open with DEXTER
In the text
thumbnail Fig. 8.

GLPS obtained using the monthly light curves of PKS 0447−439 is shown in cyan (circles); the red solid line, green line, and orange line represent the 3σ, 2.5σ, and 2σ confidence levels, respectively, calculated based on simulations of 50 000 light curves using the method of Timmer & Koenig (1995).

Open with DEXTER
In the text
thumbnail Fig. 9.

Two-dimensional WWZ of the γ-ray light curve of PKS 0447−439 based on Fermi-LAT data from August 2008 to December 2017. The color scale represents the Z-statistics of the WWZ power of a certain period at a given time. The panel shows the signature of a quasi-periodic oscillation at a (observed) period of ∼2.6 years without any significant changes over time.

Open with DEXTER
In the text
thumbnail Fig. 10.

Success fraction vs. slope (β) for all three PSDs (Fermi-LAT). The success fractions indicate the goodness of fit obtained from the PSRESP method (see text). (a) Weekly, (b) monthly light curves.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.