EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 672 (April 2023) A&A, 672 (2023) A86 Full HTML
Open Access
Issue
A&A
Volume 672, April 2023
Article Number A86
Number of page(s) 31
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/202244754
Published online 04 April 2023

© The Authors 2023

Licence Creative CommonsOpen 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.

This article is published in open access under the Subscribe to Open model.

Open access funding provided by Max Planck Society.

1. Introduction

Super-massive black holes (SMBHs) dwell at the centers of all galaxies. When accreting material, they are observed from Earth as active galactic nuclei (AGNs). In a minority of AGNs, the accretion of matter onto the black hole is associated with the launching of a pair of relativistic jets of plasma along the polar axis (see Blandford et al. 2019, for a recent review). Studying the inflows and outflows of matter on and coming from black holes provides indirect access to information on the physics of these compact objects.

The low-mass analogs of SMBHs are stellar-mass black holes produced in star collapses. Such low-mass black holes, when situated in a binary system, have also been observed accreting surrounding matter and launching plasma jets. These systems are observationally classified as X-ray binaries and microquasars. A peculiar characteristic of X-ray binaries is the presence of quasi-periodic oscillations (QPOs) in their X-ray light curves, with periods on the order of 0.1−10 Hz (see, e.g., the review by Ingram & Motta 2019, and references therein). Such a period is consistent with the innermost stable orbit of the black hole and suggests that QPOs can be used to access the region closest to the horizon and to study the physics of black hole accretion. Given the analogy between stellar-mass black holes and SMBHs, we would expect to see QPOs in the latter, although with much longer periods due to the larger size of the Schwarzschild radius that is scaled linearly with the mass of the compact object. Several observational claims have been made for QPOs in the X-ray light curves of Seyfert galaxies (radio-quiet AGNs), with periods on the order of several hours (Gierliński et al. 2008; Alston et al. 2015). The much longer periods involved in AGNs make it difficult to unequivocally measure QPOs in their light curves; thus, an uninterrupted sampling over several periods is needed in order to obtain a statistically significant measurement.

Over the last decade, the Fermi-LAT instrument (Atwood et al. 2009) has revolutionized our view of the γ-ray sky in the 100 MeV−100 GeV band thanks to its unprecedented sensitivity and all-sky monitoring observing strategy. Fourteen years after its launch, we have now access to continuous light curve data on hundreds of AGNs. In the GeV band, the extra-galactic sky is dominated by a peculiar AGN class: blazars, which are AGNs whose relativistic jet is closely aligned with our line of sight. The relativistic Doppler effect boosts the emission and makes these type of AGNs among the brightest sources in the Universe. From an observational point of view, blazars are divided into two sub-classes: BL Lacertae objects (BL Lacs for short), with weak or absent broad emission lines in their optical/UV spectrum, and flat spectrum radio quasars (FSRQs), with strong broad emission lines (Urry & Padovani 1995).

The search for periodicities in the light curves of Fermi-LAT blazars have been an active research topic since the beginning of the mission but the first claim had to wait six years while the data were gathered: Ackermann et al. (2015) presented the first evidence (at about 3σ) for a periodic two-year modulation in the light curve of the blazar PG 1553+113, coincident in period with a QPO seen at longer wavelengths. Interestingly, this first hint for a QPO in a blazar γ-ray light curve is at frequencies much lower than the ones detected in X-rays in Seyfert galaxies, suggesting that its origin is intrinsically different. Since 2015, several research groups have investigated QPOs in Fermi-LAT light curves resulting in several positive claims. All γ-ray QPO candidates in blazars shown a period of the order of years. A notable exception is represented by the highly significant (5.2σ) QPO seen in PKS 2247−131 which has a period of about a month and, most importantly, seems to be a transient phenomenon (Zhou et al. 2018).

The details of the radiative mechanism(s) at the origin of the γ-ray emission in blazars are a subject of investigation, but there is consensus on the fact that this high-energy emission is produced in the relativistic jet and far from the SMBH (see Cerruti 2020, for a recent review on blazar emission models). The physical mechanisms for the production of QPOs in γ-ray blazars will thus be intrinsically different from the ones at work in radio-quiet AGNs and, rather, linked to the physics of jets. Several models have been developed to explain these quasi-periodicities. They can be associated to the movement of plasmoids in the jet along helical paths (Rieger 2004) or related to the precession of the jet itself. This precession could be driven by the gravitational perturbation of another SMBH (Begelman et al. 1980; Sobacchi et al. 2017), meaning that QPOs will provide key constraints on SMBH binaries in the Universe. SMBH binaries are expected to form during galaxy mergers and, thus, it is clear how important their detection is for our understanding of Galaxy and SMBH evolution through the history of the Universe. Alternatively, the periodicities could be related to a regular change in the accretion of matter onto the SMBH that is then translated into a QPO in the γ-ray emission in the jet. In this case as well, SMBH binaries can be the source of this periodicity via a perturbation of the accretion flow (Valtonen et al. 2009). The SMBH binary can also imprint a QPO in the light curve via gravitational stresses by one of the black holes on the jet of the other one (Tavani et al. 2018). The most promising SMBH binary candidate identified thanks to a twelve-year QPO in its optical light curve is the blazar OJ 287 (Sillanpaa et al. 1988). The source is a known γ-ray emitter, but its famous 12-yr QPO has not been seen in its Fermi-LAT light curve due to the long period, although there have been claims of a γ-ray QPO, with a period of about 300 days (Kushwaha et al. 2020).

The goal of this work is to systematically study the light curves of bright Fermi-LAT AGNs in order to identify QPO candidates. In order to get access to transient QPOs, and QPOs with varying periods that could be hidden in a time-integrated power spectral density, we make use of the continuous wavelet transform. The paper is organized as follows. In Sect. 2, we discuss the selection of the targets. In Sect. 3, we provide the details of the Fermi-LAT data analysis. In Sect. 4, we describe the wavelet analysis, whose results are presented in Sect. 5. The discussion and the conclusions are in Sects. 6 and 7.

2. Sources

The only selection criterion applied is for the source to be bright enough to have continuous Fermi-LAT light curves with time bins of seven days or one month. The main issue is to avoid the effect of large gaps in the continuous wavelet transform (see Sect. 4). This aspect, together with the need to limit the total computing time, drives the choice to limit the analysis to the 34 brightest sources in the 4LAC catalog (data release 2, Ajello et al. 2020, sorting them by integral energy flux above 100 MeV and removing sources with a test significance lower than 100). The least bright source has an integral flux above 100 MeV equal to 8.95 × 10−11 erg cm−2 s−1. The only exception in our source selection is represented by PKS 2247−131, which has been manually added to our list due to the high-significance QPO detected in this source by Zhou et al. (2018). The source was not detected by Fermi-LAT before entering in an active state in 2016 and thus its long-term average flux is biased towards lower values. If we consider only its emission after 2016, this source is definitely among the brightest γ-ray AGNs and clearly fulfills our requirement of having an uninterrupted light curve. This manual addition brings the total of our sample to 35 AGNs. The sources we included in our analysis are listed in Table 1 by their right ascension (RA). The sample consists of 19 FSRQs, 15 BL Lacs, and 1 radio galaxy. For each source, together with the source name and 4FGL name, we provide the source class, coordinates in J2000 and the redshift (the three last quantities as provided in the 4LAC catalog).

Table 1.

Selected Fermi-LAT sources with the catalog name, blazar type, coordinate in J2000, and redshift extracted from the Fermi-LAT 4FGL catalog.

3. Fermi-LAT data analysis

The Large Area Telescope (LAT) on board the Fermi Gamma-ray Space Telescope is a pair conversion telescope sensitive to γ-rays in the energy band from 20 MeV to 500 GeV. Since its launch on June 11, 2008, Fermi-LAT scans the sky every ∼3 h, regularly monitoring the γ-ray emission from different sources. More details on the Fermi-LAT instrument are described in Atwood et al. (2009) and references therein.

In the current study, γ-ray data from the observation of the considered sources from August 8, 2008 to April 4, 2021 (from MJD 54686 to MJD 59308) were downloaded and analyzed using the standard analysis procedure provided by the Fermi-LAT Collaboration. The data for each of the considered sources are analyzed in an identical manner. The Pass 8 data in the energy range from 100 MeV to 500 GeV are analyzed using Fermi ScienceTools (1.2.1) and the P8R3_SOURCE_V3 instrument response functions. We selected only events within a maximum zenith angle of 90° to reduce contamination by photons from Earth’s atmosphere and we used the filter expression (DATAQUAL > 0) and (LAT CONFIG == 1) to select good time intervals. The photons from a circular region with a radius of 12° around each source under consideration were selected. Then, the events were binned within a 16.9° ×16.9° square region of interest (ROI) into 0.1° ×0.1° pixels and into 37 equal logarithmically spaced energy bins. The model for which the likelihood is computed includes a combination of point-like and diffuse sources and standard templates describing the diffuse emission from the Galaxy (gll_iem_v07) and the isotropic γ-ray background (iso_P8R3_SOURCE_V3_v1). The model file for each source is created using the Fermi-LAT fourth source catalog (4FGL; Abdollahi et al. 2020), where all sources falling within ROI+5° were all included with the same spectral models as in the catalog. The normalization parameters of the diffuse Galactic and the isotropic component, and both the normalization and spectral parameters of the point sources within the ROI were set as free parameters, while those of the sources outside ROI were fixed to the catalog values. The binned likelihood analysis implemented in the gtlike tool is used to optimize the parameters to best match the observations for the whole time period.

The model file obtained from the full time analysis is used to compute the light curves. We produced the light curves binned into weekly and monthly intervals by applying unbinned likelihood analysis in the 0.1−300 GeV energy range with the appropriate quality cuts mentioned above.

4. Wavelet analysis

4.1. Continuous wavelet transform

In this work we make use of the continuous wavelet transform (CWT) technique, which is the convolution of a time series with a dilated and translated wavelet function, to analyze the time frequency properties of the light curves. Acting as a band-pass filter, the CWT maps the power of any particular periodic behavior at different times in the time-frequency space. This is notably useful, since the CWT technique not only gives access to the frequencies of potential QPOs, but also when the periodicities appear and end and how they evolve in time.

We used the Python implementation PyCWT provided by Torrence & Compo (1998) and as mother wavelet the Morlet wavelet. This wavelet consists of a plane wave modulated by a Gaussian and is shown in Appendix A. The wavelet power spectrum is defined as the square of the amplitude of the wavelet coefficient. The global wavelet spectrum can then be computed as the time-average of the wavelet spectrum. The global wavelet spectrum provides an unbiased and consistent estimation of the true power spectrum of a time series (Percival 1995). Because of the finite length of the time series, border effects occur at the edges of the wavelet power spectrum. The cone of influence (COI) is defined as the region of the wavelet power spectrum in which edge effects become important. It is calculated as suggested by Torrence & Compo (1998).

In Appendix A, we show the CWT map and global wavelet spectrum of a periodic signal with two periods at 0.1 and 0.02 s (Fig. A.2) and a Dirac delta signal (Fig. A.3). The CWT has an intrinsic resolution in reconstructing a periodic signal that can be visualized as broad bands in the CWT map. We estimate this uncertainty as half of the full-width at half-maximum (FWHM) of the global wavelet spectrum. For the Dirac delta function, a vertical feature appears in the wavelet power spectrum which reveals the response of this technique when a flare-like signal is present in the light curve. This response of the CWT to flares is of particular relevance for our study. A high power in the CWT by itself does not mean that a periodicity is present in the light curve: all CWT maps have to be inspected to visually confirm that the power is distributed horizontally and is not related to a flaring behavior.

4.2. Significance estimation

To determine the significance and confidence levels of the analysis, we simulated artificial light curves following the work by Emmanoulopoulos et al. (2013), using the Python version provided by Connolly (2015). This algorithm generates artificial light curves having the same power spectral density (PSD) and probability distribution function (PDF) as the original light curve. It represents an improvement of the procedure of Timmer & Koenig (1995), which produces normally distributed time series from a given PSD. The algorithm of Emmanoulopoulos et al. (2013) is able to precisely reproduce light curves which match both the PSD and the PDF of a given observed light curve – or a theoretical model, where the PSD estimate is performed using a maximum likelihood methodology and assuming a smoothly bending power-law model plus a constant.

For each Fermi-LAT light curve, we produced 10 000 artificial light curves. In Appendix B, we show a comparison of the PDF and PSD of one of the simulated light curve to the original one. The global wavelet spectrum is computed for each simulated light curve, such that a histogram of the power spectrum can be produced at every period or scale. We then fit the histogram with a χ2 distribution with k degrees of freedom:

χ ( x , k ) 2 = 1 2 k / 2 Γ ( k / 2 ) x k / 2 − 1 exp ( − x / 2 ) . $$ \begin{aligned} \chi (x,k)^2 = \frac{1}{2^{k/2} \Gamma (k/2)} x^{k/2-1} \exp (-x/2). \end{aligned} $$(1)

The confidence levels are obtained by using the percentiles of the power for each scale, which define the global significance of the results. An example of the χ2 fitting to the histogram is included in Appendix C, as well as the resulting confidence levels for the AGN S5 0716+714.

When searching for significant results in physics, it is necessary to take into account the “look-elsewhere effect”, which can be quantified in terms of a specific number of trials. We estimated the post-trial confidence levels of the global and local wavelet spectrum following the work presented by Auchère et al. (2016). Taking into account the fact that the bins in the wavelet spectrum are not statistically independent, these authors showed that the post-trial probability, PG, can be computed from the pre-trial one PL as:

P G = 1 − ( 1 − P L a ) n , $$ \begin{aligned} P_G = 1-(1-P_L^{a})^n, \end{aligned} $$(2)

where a and n are empirically derived coefficients, which are related to the number of bins in the spectrum and the resolution of the CWT that we are using (δj), and that are specific for each mother wavelet. The resolution δj in our analysis is 1/12. Different parameterization for the global and local wavelet spectrum are needed, since in the local one, trials are made in the time-period bins of the power spectrum map, whereas after the time-average, only the period bins account for the global wavelet spectrum. Auchère et al. (2016) parameterized a = 0.810 (Nout δj)0.011 and n = 0.491 (Nout δj)0.926 for the local wavelet spectrum, and a = 0.805 + 0.45 × 2−Sout δj and n = 1.136 (Sout δj)1.2 for the global one, where Nout and Sout are the time-period bins and the period bins of the wavelet spectrum outside the COI, respectively. These coefficients are valid for time series affected by power-law noise, but adequate only for a wavelet analysis using the Morlet wavelet (for more details see Auchère et al. 2016).

Additionally, we are considering 35 sources in this work and two time-binnings for each of them. Therefore, for each time-binning, we computed the number of bins, Nout and Sout, and we multiplied them by the number of sources in order to account the total number of trials carried out. One special case is the AGN PKS 2247−131, whose light curve is cut at MJD 57427, since there was no detection by Fermi-LAT before that date. Hence, the number of time-period bins, Nout, for this source is lower.

The total number of bins Nout, T and Sout, T are thus computed as:

N out , T = ( 34 × N out , 1 month + N out , 1 month PKS 2247 − 131 ) + ( 34 × N out , 7 days + N out , 7 days PKS 2247 − 131 ) = 1698893 , S out , T = ( 34 × S out , 1 month + S out , 1 month PKS 2247 − 131 ) + ( 34 × S out , 7 days + S out , 7 days PKS 2247 − 131 ) = 4833 . $$ \begin{aligned} \begin{aligned} N_{\rm out,T} =&(34 \times N_{\rm out,1\,month} + N^\mathrm{PKS\,2247{-}131}_{\rm out,1\,month})\\&+ (34 \times N_{\rm out,7\,days} + N^\mathrm{PKS\,2247{-}131}_{\rm out,7\,days}) = 1698893, \\ S_{\rm out,T} =&(34 \times S_{\rm out,1\,month} + S^\mathrm{PKS\,2247{-}131}_{\rm out,1\,month})\\&+ (34 \times S_{\rm out,7\,days} + S^\mathrm{PKS\,2247{-}131}_{\rm out,7\,days}) = 4833. \end{aligned} \end{aligned} $$(3)

5. Results

In this section, we present the results for those light curves showing possible QPOs and we discuss in more detail the five sources exhibiting the most significant features: 4C +01.02, PKS 0537−441, S5 1044+71, B2 1520+31, and PKS 2247−131. Furthermore, results that are consistent with previous QPO searches are also found in several sources, however at a lower level of significance: PKS 0426−380, S5 0716+714, Mrk 421, PKS 1424−418, Mrk 501, and PKS 2155−304. These results are also commented on briefly in this section, while their CWT maps are provided together with all the remaining sources in Appendix D. Table 2 lists the candidate QPOs found by CWT.

Table 2.

QPO candidates identified by the CWT of the Fermi-LAT light curves in time bins of 30 days and 7 days.

We identified a total of 36 QPO candidates in 24 out of the 35 selected sources. We considered candidates that have a post-trial significance larger than 3σ in at least one of the two binnings and show at least three complete cycles once fitted. Long-term QPO candidates with periods longer than and around one year are observed at various significance levels in five sources: S5 0716+714, S5 1044+71, Mrk 421, Mrk 501, PKS 2155−304, and CTA 102. On the other hand, month-long QPO candidates with periods smaller than ∼300 d can be identified in the following 20 sources: 4C +01.02, 4C +28.07, NGC 1275, PKS 0402−362, PKS 0426−380, PKS 0447−439, PKS 0454−234, PKS 0537−441, 1H 1013+498, S5 1044+71, 4C +21.35, 3C 273, 3C 279, PKS 1424−418, PKS 1510−089, B2 1520+31, CTA 102, PKS 2247−131, 3C 454.3, and PMN J2345−1555. We reported no QPO detection in the γ-ray emission of sources 3C 66A, PKS 0235+164, 4C +55.17, B3 1343+451, PKS 1424+240, PKS 1502+106, PG 1553+113, 4C +38.41, 1ES 1959+650, and BL Lac and PKS 2326−502.

The wavelet power spectra of the five more significant sources are shown in Figs. 1–5. Additionally, to help in visualizing the QPO candidates identified by the CWT, we fit periodic functions with periods equal to the ones computed via the CWT to the Fermi-LAT light curves. A detailed description of these results is given below.

thumbnail Fig. 1.

CWT map for the monthly binned (upper) and the weekly binned light curve (bottom) for 4C +01.02 (left). In each subplot, the panels represent (a) Fermi-LAT light curve, (b) wavelet power spectrum and (c) global wavelet power spectrum. The solid coloured contours in (b) and the dashed coloured lines in (c) are the confidence levels (1 to 5σ in black, blue, green, violet, and red). Monthly binned and weekly binned light curves with the fitted periodic signal in red dashed lines in the upper panels and the relative error in the bottom panels (right).

thumbnail Fig. 2.

CWT map for the monthly binned (upper) and the weekly binned light curve (bottom) for PKS 0537−441 (left). In each subplot, the panels represent (a) Fermi-LAT light curve, (b) wavelet power spectrum and (c) global wavelet power spectrum. The solid coloured contours in (b) and the dashed coloured lines in (c) are the confidence levels (1 to 5σ in black, blue, green, violet, and red). Monthly binned and weekly binned light curves with the fitted periodic signal in red dashed lines in the upper panels and the relative error in the bottom panels (right).

thumbnail Fig. 3.

CWT map for the monthly binned (upper) and the weekly binned light curve (bottom) for S5 1044+71 (left). In each subplot, the panels represent (a) Fermi-LAT light curve, (b) wavelet power spectrum and (c) global wavelet power spectrum. The solid coloured contours in (b) and the dashed coloured lines in (c) are the confidence levels (1 to 5σ in black, blue, green, violet, and red). Monthly binned and weekly binned light curves with the fitted periodic signal in red dashed lines in the upper panels and the relative error in the bottom panels (right).

thumbnail Fig. 4.

CWT map for the monthly binned (upper) and the weekly binned light curve (middle) for B2 1520+31 (left). In each subplot, the panels represent (a) Fermi-LAT light curve, (b) wavelet power spectrum and (c) global wavelet power spectrum. The solid coloured contours in (b) and the dashed coloured lines in (c) are the confidence levels (1 to 5σ in black, blue, green, violet, and red). Monthly binned and weekly binned light curves with the fitted periodic signal in red dashed lines in the upper panels and the relative error in the bottom panels (right). The third CWT is computed with the weekly binned light curve cut from MJD 54697 to MJD 56100. And on the right, the respective fitted light curve shows the ∼39 d period.

thumbnail Fig. 5.

CWT map for the monthly binned (upper) and the weekly binned light curve (middle) for PKS 2247−131 (left). In each subplot, the panels represent (a) Fermi-LAT light curve, (b) wavelet power spectrum and (c) global wavelet power spectrum. The solid coloured contours in (b) and the dashed coloured lines in (c) are the confidence levels (1 to 5σ in black, blue, green, violet, and red). Monthly binned and weekly binned light curves with the fitted periodic signal in red dashed lines in the upper panels and the relative error in the bottom panels (right). In the bottom fitted light curve we show a possible reappearance of the candidate QPO with ∼34 d period at MJD 58050, adding three more cycles to the series.

5.1. 4C +01.02

The wavelet analysis shows two significant QPO candidates in the light curves of 4C +01.02, a distant FSRQ at z = 2.099 (see Fig. 1a,c). No previous report of QPO analysis has been found for this source in the literature. Both QPOs have month-long periods, with the longer one centered at 268 ± 55 d (268 ± 54 d), and the shorter one at 123 ± 26 d (122 ± 26 d), in the monthly (weekly) binned light curves. The significance of the longer QPO is over 5σ in both time binnings, whereas the shorter QPO has a significance ∼4.7σ for the monthly binned light curve, and larger than 5σ for the weekly binned light curve. The fitted light curves in Fig. 1b,d present four complete cycles for the first QPO, starting from around MJD 56900 to MJD 58000 (August 2014 to September 2017), and five complete cycles for the second QPO, between MJD 57300 and MJD 57900 approximately (October 205 and May 2017). Tinier structures also appear at lower periods, although with more vertical and ambiguous shapes, which are inevitably produced by small flares (see for example the signals found between MJD 57500 and MJD 58000 at periods < 50 d). We did not include these among our results due to their shortness and unclear shape, but this is a very subjective criteria and these features could be investigated further with finer time bins.

5.2. PKS 0537–441

We identified a significant QPO candidate in this BL Lac object at z = 0.892 (see Fig. 2a,c). The wavelet spectrum shows a clear horizontal feature centered at period 285  ±  67 d (286  ±  73 d) in the monthly (weekly) binned light curve, at a significance above 5σ in the CWT map of both time binnings. This result is consistent with Sandrinelli et al. (2016), who claimed a candidate QPO at ∼280 d in the first few years of Fermi-LAT data, using the Lomb Scargle periodogram (LSP) technique. We can see four complete cycles of oscillation during a period of high γ-ray flux, starting from the beginning of Fermi-LAT observations till around MJD 55700 (May 2011). The cycles are also shown in the fitted light curves in Fig. 2b,d.

5.3. S5 1044+71

We identified, for the first time, two highly significant QPO candidates in S5 1044+71, a moderately distant (z = 1.15) FSRQ. The two CWTs are shown in Fig. 3a,c. A very evident long-term oscillation of period 1133 ± 229 d (1127 ± 226 d) at a confidence level of ∼4.9σ (∼4.6σ) for the monthly (weekly) binned light curve appears in the wavelet map, emerging around MJD 56000 (March 2012). Secondly, a short month-long QPO emerged during the last flaring state of the source, between MJD 58650 and MJD 59000 (June 2019 and September 2020), with a period estimated to be 116 ± 33 d (117 ± 38 d) for the monthly (weekly) binned light curve at a significance of more than 5σ. Three complete cycles of the year-long QPO can be observed also in the fitted light curve in Fig. 3b, which corresponds to the large flaring states of the FSRQ. The last maximum occurred at around MJD 58800 (November 2019), so we expect the next maximum to occur at approximately MJD 59900 (November 2022). The absence of a clear fourth peak at around MJD 55500 might indicate additional modulation1. Regarding the short period QPO, four cycles are indicated in the fitted light curve in Fig. 3d, which are the three narrower flares in the last flaring state and a fourth more suppressed one after MJD 59000. This kind of structured flares are not unique in this last flaring state, but could also appear in the two previous maxima. However, these signals exhibit more complex and vertical shapes, which leads to ambiguity. If real, these could also suggest a trend toward longer periods. We can also notice in the weekly CWT maps hints (2σ at most) of several horizontal features at longer periods, which nonetheless remain consistent with noise once corrected for trials.

5.4. B2 1520+31

Three possible month-long QPOs were identified by the CWT of the light curves of this distant (z = 1.489) FSRQ, as shown in Fig. 4a,c,e. The CWT maps of this source are a good example of the inherent difficulty of this analysis technique when dealing with blazar light curves: rapid flares result in vertical structures in the map that overlap with the horizontal bands we are interested in and require visual inspection of all peaks that appear in the global wavelet spectrum. The first QPO candidate is found in both light curves at around 176  ±  48 d (179 ± 42 d) for the monthly (weekly) binned light curve with a significance above 5σ in both cases. This QPO candidate is of particular interest because its period seems to increase with time. The second and third candidate QPOs can be better identified in the CWT map of the weekly binned light curve, with periods of 71 ± 15 d and 39 ± 11 d, both exceeding a 5σ significance. A zoomed-in CWT is shown in Fig. 4e, which allows us to better visualize the signal of the two very-short-period QPO candidates. The 71 ± 15 d period is compatible with the one reported by Gupta et al. (2019) with a period of ∼71 d, by analysing the first four years of Fermi-LAT data of B2 1520+31, using the LSP and weighted wavelet z-transform (WWZ) techniques. This source has also been studied by Tarnopolski et al. (2020): although some signal can be seen in their analysis at around ∼70 d, it remains below the 3σ interval they computed in their study.

Fitted light curves are presented in the right hand side of Fig. 4, where B shows a periodic function of ∼176 d with six cycles between MJD 54900 and MJD 56000 (March 2009 and March 2012); in Fig. 4d we can see 14 oscillations of ∼71 d spanning from MJD 55000 to MJD 56000 (June 2009 to March 2012); and Fig. 4E shows a fit with 17 cycles of period ∼39 d between MJD 55300 and 55970 (April 2010 and February 2012). A forth QPO could be also considered at a period halfway between ∼71 d and ∼176 d, which is visible in both the CWT map and the global wavelet spectrum of the weekly binned light curve. Furthermore, the ∼39 d period seems to also be present before the above mentioned time span, between around MJD 54800 and MJD 55000.

5.5. PKS 2247–131

The wavelet analysis of PKS 2247−131, a BL Lac object at z = 0.22, shows two QPO candidates (see Fig. 5a,c). The first of the QPO candidates is found to be at around 217 ± 38 d (214 ± 43 d), at above a 5σ confidence level in the CWT of monthly (weekly) binned light curves. This QPO candidate seems to span at least from MJD 57600 to MJD 58500 approximately (July 2016 to January 2019) in the CWT map. Just by looking at the light curve, we can notice the tentative periodic oscillations of this source, and this can be confirmed by the fitted function shown in Fig. 5b, presenting six peaks that include the one at around MJD 58750.

More interestingly, we identified a QPO candidate in the time interval around MJD 57600 to MJD 57900, approximately (July 2016 to May 2017), with a much shorter period. We can appreciate this further in Fig. 5c, which shows the QPO centered around 34 ± 13 d, with a significance larger than 5σ and only noticeable in the full-range weekly binned light curve. Tentative oscillations can also be seen, with at least seven full cycles, shown in Fig. 5d with the fitted periodic function. In addition, with a more detailed evaluation one can notice that the ∼34 d period reappears after MJD 58000 (September 2017). This is also shown in the fitted light curve in Fig. 5d, with three more cycles spanning MJD 58050 to MJD 58170, approximately.

In November 2018, two years after the Fermi-LAT announcement of detection of this BL Lac type blazar, Zhou et al. (2018) presented the first claim of a month-scale QPO in the γ-emission of PKS 2247−131. They found a relatively short, month-scale oscillation at period 34.5 ± 1.5 d, which can be indicative of the presence of an SMBH binary in the center of this blazar. This discovery conducted us to explore the CWT of PKS 2247−131 more deeply. Small features like this one also appear in several other sources which (given in Appendix D), but we do not systematically inspect all of them. The short QPO candidate in our analysis is compatible with the aforementioned one proposed by Zhou et al. (2018). On the other hand, no previous claims for the longer QPO candidates at periods of around 220 d have been reported in the literature.

5.6. Other sources

In the remaining 27 sources showing QPO candidates, we found several candidates compatible with previously reported claims in the literature. The CWT maps and fitted light curves of these AGNs can be found in Appendix D.

S5 0716+714. We found a QPO candidate in both the one-month and seven-day binned light curves of this BL Lac object, centered at 325 ± 75 d (324 ± 77 d), at ∼2.4σ (∼3.2σ) in the wavelet spectrum of the monthly (weekly) binned light curve. This feature only appears in the time interval between MJD 56000 and MJD 57500, approximately (March 2012 to April 2016). Our result is compatible with the analyses from several authors: Prokhorov & Moraghan (2017) reported a QPO evidence at ∼346 d, Li et al. (2018) detected a periodic feature at ∼344 d, Bhatta & Dhital (2020) found a QPO candidate at around 340 d and Peñil et al. (2020) identified a period at ∼0.9 yr. However, Covino et al. (2019) and Żywucka et al. (2022) claim the absence of any periodic emission in this object. Bhatta & Dhital (2020) and Peñil et al. (2020) also reported a possible QPO at longer period around 1000 d. In our CWT analysis, a high power spectrum feature does appear and it is centered at around 1000 d, spanning almost the full time series. However, this result is only significant before applying the trial correction and is compatible with noise once the “look-elsewhere effect” is taken into account.

Mrk 421. Bhatta & Dhital (2020) reported an oscillating γ-ray emission at ∼280 d, compatible with our QPO candidate at 300 ± 64 d (300 ± 65 d), with a significance exceeding 5σ in both time binnings. This QPO candidate is also not a persistent one, having been seen between MJD 55800 and MJD 57000 approximately (August 2011 to December 2014). In our analysis, both maps show a complex structure, with a possible period change at the end.

Mrk 501. Bhatta (2019) found a candidate QPO at period ∼230 d, compatible with the one we identify at 315 ± 98 d (326 ± 76 d) with a significance ∼2.9σ in the monthly binned light curve, but a high significance above 5σ in the weekly binned one. This QPO candidate is also temporary, emerging from MJD 55800 to MJD 57000 approximately, and reappearing between ∼MJD 57800 to ∼MJD 58700.

We want to remark the following three sources which also present hints of QPOs that are compatible with some claims in the literature, although they are only significant before trial corrections.

PKS 1424−418. Bhatta & Dhital (2020) and Yang et al. (2021a) claimed a flux oscillation at period ∼353 d and ∼355 d, respectively, in the light curve of this FSRQ. The wavelet analysis, however, shows that this QPO candidate is compatible with noise, with a post-trial significance at ∼1σ for both monthly and weekly binned light curves. On the other hand, we detected a much shorter QPO candidate centered at ∼94 ± 25 d (90 ± 22 d) in the CWT of the monthly (weekly) binned light curve, at above a 5σ significance. This month-long QPO candidate emerged from MJD 56100 to MJD 56500, approximately (June 2012 to August 2013), showing five complete oscillations in the fitted light curve.

PG 1553+113. The first detection of a periodic behavior in the γ-emission of this high-synchrotron-peak BL Lac object was claimed by Ackermann et al. (2015), with a period of ∼798 d. Some more recent studies found a periodic feature at similar values: Sandrinelli et al. (2018) reported two peaks at same frequency ∼780−810 d, Covino et al. (2020) detected a period of ∼790 d, Yang et al. (2021b) claimed a possible evidence at period ∼800 d, and Peñil et al. (2020) found high significance level at a period of ∼803 d. Our results show that this QPO candidate is significant in pre-trial significance (at a confidence level of ∼3σ (∼2σ) for the monthly (weekly) binned light curve), but it becomes less than 1σ after trial correction (consistent with the results by Ait Benkhali et al. 2020).

PKS 2155−304. We identify a low-significance QPO candidate for this BL Lac object centered at around 334 ± 107 d (341 ± 106 d), with a significance of ∼2.2σ (∼3.5σ) in the monthly (weekly) binned light curve. The QPO candidate clearly shows a time-dependent behavior, with a decreasing frequency over time. This source has been extensively studied in the past and several QPO claims have been made: Zhang et al. (2017a) detected a possible QPO at period ∼640 d; Bhatta & Dhital (2020) found two quasi-periodic features, one at 610 ± 51 d and the other at ∼260 d; Covino et al. (2020) found a QPO at ∼610 d, albeit not significant enough; Peñil et al. (2020) claimed a QPO at ∼1.7 yr; Tarnopolski et al. (2020) found a QPO at 610 ± 42 d; and Żywucka et al. (2022) detected a QPO at ∼612 d. With respect to the results of our analysis, the QPO found by CWT technique is close to the short QPO candidate presented by Bhatta & Dhital (2020), whereas the long-term feature at ∼600 d, also noticeable in the CWT maps, is not significant after the trial correction. The case is the same for the long ∼1000 d QPO candidate in S5 0716+714.

6. Discussion

We identified a total of 36 QPO candidates in 24 sources, with the longest one found in the light curve of S5 1044+71 at ∼1130 d and the shortest in the light curve of PKS 2247−131 at ∼34 d. Many of the candidates have a period between one month and one year. This is a new result since not many month-year-long QPOs have been reported in the literature. It might be reveal the fact that middle-term QPOs are actually frequent in nature. On the other hand, only four candidate QPOs with period of around one month are detected in our results. This can be explained by the fact that these kinds of short periods, if not lasting continuously in time, would appear in the CWT map as small structures that are difficult to recognize and disentangle from vertical structures produced by flares. A good example of this is PKS 2247−131, whose analysis was inspired by the work published by Zhou et al. (2018). More examples of short-leaving features with periods of the order of one month can be seen in the CWT maps of B2 1520+31, 1H 1013+498, and 3C 279 (see Sect. 5.4 and Appendix D). We recall that our threshold for selecting candidates is a post-trial significance of > 3σ in at least one of the time-bins and exhibiting more than three cycles. Thus, we cannot exclude the possibility that some of the candidates are spurious and come from the analysis itself.

Quasi-periodic modulations observed in the high-energy γ-ray fluxes of AGNs should be related to the relativistic jets launched by these objects or to the process feeding the jet itself. Two main subgroups of QPOs have often been discussed in the literature in the context of the oscillation period: short intra-day and year-long QPOs. The former are though to be associated to pulsational accretion flow instabilities (Honma et al. 1992), despite the fact that longer periods are observed in magnetohydrodynamic simulations for slowly spinning SMBHs (McKinney et al. 2012). The origin of long-period QPOs, on the other hand, could be related to jet precession and geometry, as well as (possibly) the presence of a binary SMBH system (see Ackermann et al. 2015, and references therein). In simple SMBH systems, jet precession, rotation, and helical structure, in the presence of a sufficiently strong magnetic field, yield to the observable periodicity from the change of the line of sight. Moreover, these variabilities could also appear in binary SMBH systems due to periodic perturbations of the secondary compact object to the accretion disk and jet. The orbital period of the binary systems is in the range of several years, which is compatible with the year-long period QPOs. The most significant, multi-year, QPO candidate in our analysis, seen at more than 4σ in the blazar S5 1044+71, has a period of about 1100 d (3.5 yr); thus, this γ-ray source a high-significance SMBH binary candidate in the Universe.

No QPOs at an intermediate time scale were reported until the publication of Zhou et al. (2018), who claimed the detection of ∼34.5 d QPO in the light curve of PKS 2247−131. This monthly modulation suggested a helical jet structure due to the short period. However, these authors also noted that the helical structure could be driven by the orbital period of a secondary SMBH. By considering the time compression due to the high Doppler factor of the emitting region, the observed period will be shortened with respect to the physical period in the host galaxy reference frame and even then, it may still be on the order of typical orbital periods for close binary SMBHs. In our analysis of PKS 2247−131, an additional QPO (∼200 d) seems to be existing in the form of large γ-ray flares during the last few years, and could be used to provide further constraints on the model.

Below, we highlight three important findings from our analysis that could be used as key observables to understand the origin of QPOs in blazar γ-ray light curves.

We did not detect persistent QPOs, lasting for the whole observing period. None of the global wavelet power spectra show significant features after trial correction. It is possible for S5 1044+71 to be considered a persistent QPO candidate if we make the hypothesis that the first maximum is suppressed due to additional modulation. Despite this potential exception, all other QPO candidates in our analysis are transient ones.

We identified some cases where the QPO candidate shows period shifts. This can be seen in the CWT maps of B2 1520+31 and, at a lower level of significance, of PKS 2155−304, where one of the QPO candidates is decreasing in frequency. In hunting for a possible explanation, we might consider changes in the inclination angle of jet precession. However, the reasonable timescales involved should not produce a transition as fast as the observed one. A much more reasonable possibility is that we are observing the helical geometry of a jet with an intrinsic opening angle, which would naturally lead to a slowdown of the QPO.

Lastly, we see several occurrences of multiple QPO candidates occurring simultaneously and with harmonic periods. One of the examples is 4C +01.02, where two QPOs with different frequencies overlap in a specified time interval. The longer period is approximately two times the shorter one, suggesting the presence of resonances in the emission. Simultaneous QPO candidates at harmonic ratios (within errors) can be seen in PKS 0402−362, 1H 1013+498, B2 1520+31, and CTA 102. Harmonics are common in QPO analysis of X-ray binaries (Ingram & Motta 2019), so it is not surprising to see them also in AGNs. They indicate that the various QPOs share the same origin and are not due to independent physical processes, each one imprinting its particular QPO on the light curve. Theoretical models aimed at explaining QPOs in γ-ray light curves of blazars can thus be further constrained by studying harmonics and harmonic ratios.

Furthermore, no QPO candidates were identified in the γ-ray emission of 3C 66A, PKS 0235+164, 4C +55.17, B3 1343+451, PKS 1424+240, PKS 1502+106, PG 1553+113, 4C +38.41, 1ES 1959+650, BL Lac, and PKS 2326−502. Three of them, namely, PG 1553+113, BL Lac, and 3C 454.3, were previously reported showing periodic oscillations.

PG 1553+113. In the previous section, we discuss how the wavelet analysis does show a increase in power spectrum at a period similar to the one reported by Ackermann et al. (2015), Sandrinelli et al. (2018), Covino et al. (2020), Yang et al. (2021b), and Peñil et al. (2020; ∼800 d), but with a very low post-trial significance, that does not even reach 1σ. In order to check the consistency of our result with the literature, we made a test with a reduced light curve (removing the last cycle). The result shows that although the significance rises a bit, it still remains below 3σ post trial. The major effect must come from the trial correction, which was not done when analyzing a single source. And furthermore, when increasing the duration of the light curve, the number of trials also increases. Before trial correction, we find a significance that is not much different, as compared to other works in the literature.

BL Lac. Sandrinelli et al. (2017, 2018) identified a candidate QPO with period of ∼680 d in BL Lac’s light curve. Nevertheless, our results show that this possible QPO is compatible with red noise, as was also found by Covino et al. (2019) and Peñil et al. (2020).

3C 454.3. Sarkar et al. (2021) claimed a > 4σ QPO at period ∼47 d analyzing the one-day binned light curve of 3C 454.3, lasting from MJD 56800 to MJD 57250. We can observe in our CWT maps of this source that there is a strong vertical signal before MJD 56000 and no further significant feature appears afterwards. Thus, our CWT maps of 3C 454.3 show no evidence of such a QPO, even for the seven-day binned light curve.

Finally, some other γ-ray AGN have been claimed to show QPOs with periods of the order of years, but were not included in our sample only because they did not pass our original cut on brightness. For completeness, we mention here the cases of PKS 0301−243 (Zhang et al. 2017b), with a period of 2.1  ±  0.3 yr, PKS 0521−36 (Zhang et al. 2021), with a period of about 1.1 yr, PKS 0601−70 (Zhang et al. 2020), with a period of 1.22  ±  0.06 yr, and OJ 287, with a period of about 314 d (Kushwaha et al. 2020).

Furthermore, we recall that our work is limited by two main choices: the target list and time binning. Extending the target list towards less bright objects is certainly possible, although it comes with a major price: at some point, the Fermi-LAT light curves of these fainter AGNs will start to show non-detections in individual bins, and the presence of these zeros is problematic for the CWT technique that is sensitive to discontinuities. A simple solution will then be to also investigate longer time bins. Finer time bins could also be investigated for very bright γ-ray flares, giving access to QPO searches on time-scales of days. Both options come with a price, namely, an increase in the number of trials. As a final caveat, our results and, in particular, the determination of the respective significances depend on the Monte Carlo simulations of artificial light curves. In our study, we work under the assumption that the PSDs of the original light curves can be reconstructed by fitting them with a smoothly bended powerlaw function, but this choice is far from unique. For a comprehensive study of Fermi-LAT PSD using different methods, we refer to the recent work by Tarnopolski et al. (2020).

7. Conclusions

The search for QPOs in the light curves of AGNs is a major research topic in astrophysics, providing additional constraints on the physics of the SMBHs that power these systems. Long-term QPOs, with periods on the order of months and beyond, are particularly difficult to identify due to the need of highly-sampled and unbiased light curves over long periods of time. The Fermi-LAT γ-ray telescope, thanks to its monitoring capabilities, is an ideal instrument for performing such a study.

We analyzed a 13-year-long (from August 2008 to April 2021) Fermi-LAT γ-ray light curve in two different time binnings (7 and 30 days) of 35 bright γ-ray AGNs. By using the CWT technique, we systematically searched for QPO candidates in this data set. In order to compute the confidence levels of the QPO candidates, 10 000 simulated artificial light curves were generated for each light curve and the histograms of global power spectrum at each period scale were fitted to a χ2 function. We corrected for the trial effect in our analysis by estimating the trial number due to the number of sources and the number of time-period bins (in the wavelet power spectrum) or the number of period bins (in the global wavelet spectrum), outside the COI, and following the parameterization given by Auchère et al. (2016).

In this way, 36 QPO candidates in 24 sources are identified (at various significance levels), with periods ranging from one month to several years. Our most significant, multi-year QPO candidate is in the blazar S5 1044+71, with a period of about 1100 d.

We confirmed some previously claimed γ-ray candidate QPOs in the sources PKS 0537−441, B2 1520+31, PKS 2247−131, S5 0716+714, Mrk 421, Mrk 501, and PKS 2155−304, while new possible QPO candidates were also identified in 4C +01.02, 4C +28.07, NGC 1275, PKS 0402−362, PKS 0426−380, PKS 0447−439, PKS 0454−234, 1H 1013+498, S5 1044+71, 4C +21.35, 3C 273, 3C 279, PKS 1424−418, PKS 1510−089, B2 1520+31, CTA 102, PKS 2247−131, 3C 454.3, and PMN J2345−1555.

The possible physical origins of these quasi-periodical emissions are the precession of the AGN jet, its helical structure, and changes in the accretion flow. These scenarios can be due to the presence of a second SMBH in the system, opening up the window for studying SMBH binaries. Since the orbital periods of the SMBH binaries are of the order of several years, sources showing long-term QPOs are naturally suspected as binary candidates. In particular, we put forward the blazar S5 1044+71 as the most promising SMBH candidate in our sample, due to its high significance as a QPO candidate with a period of about 1100 d. Shorter, month-long QPOs could also be related to close SMBH binaries, once the period is corrected for the Doppler factor of the jet. In our analysis, we identified a peculiar behavior in some QPO candidates that might hint towards a QPO origin related to the jet geometry, that is, a varying (slowing-down) QPO frequency, seen in B2 1520+31. We also put forward the possibility that some simultaneous QPO candidates, such as the ones seen in 4C +01.02, are indeed harmonics, with a single physical mechanism at their origin.

With this, we can conclude that the CWT technique is a very powerful tool sensitive to any periodic oscillations in the light curves, considering an appropriate time binning. It has a major advantage over other statistical tools in being sensitive to transient QPO and period-shifting QPOs. However, it also reacts strongly to flares, resulting in vertical features in the map, and is influenced notably by border effects at large periods. This leads to the requirement of a visual inspection of all CWT maps to avoid misleading results. A major point that ought to be highlighted is the trial-correction of significances. Due to the large number of scales probed, the number of sources, and the number of time binnings investigated, the number of trials is pretty large and a careful correction has to be implemented to avoid false positives. While this correction is clearly mandatory for systematic studies as ours, it is also necessary for studies that analyze a single source without justifying how a particular target has been selected among the Fermi-LAT catalog.

A natural perspective for a future study would be a multi-wavelength analysis of quasi-periodic emissions from the selected interesting γ-ray sources. First, the identification of QPOs at multiple wave bands with the same periodicity will automatically boost the significance of the detection or, if the QPO is only seen in γ-rays, it can help identify false positives. Second, and more importantly, multi-wavelength observations will constrain the theoretical models that aim to explain these QPOs. Long-term, unbiased, multi-wavelength monitoring campaigns over several years, as complex and expensive as they might seem, are the key to identifying AGN QPOs and understanding their origin.

1

During the writing of this manuscript, an analysis of S5 1044+71 was presented in a pre-print by Wang et al. (2022), claiming a three-year modulation (∼3.06 ± 0.43 yr) at a significance level of ∼3.6σ. This result is in total agreement with ours.

Acknowledgments

The authors wish to thank Jonathan Biteau and Pablo Peñil for fruitful inputs and discussions, and Josep Maria Paredes and Bruno Khélifi for comments on the early draft, as well as the anonymous referee for their constructive comments that significantly improved the manuscript. This work has been possible thanks to the open source code for the simulation of artificial light curves provided by Emmanoulopoulos et al. (2013), its python version provided by Connolly (2015) and the open source Python library PyCWT elaborated by Torrence & Compo (1998). M.C. has received financial support through the Postdoctoral Junior Leader Fellowship Programme from la Caixa Banking Foundation, grant No. LCF/BQ/LI18/11630012. N.S. acknowledges the support by the Science Committee of RA, in the frames of the research projects No. 20TTCG-1C015 and 21T-1C260.

References

  1. Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33 [Google Scholar]
  2. Ackermann, M., Ajello, M., Albert, A., et al. 2015, ApJ, 813, L41 [NASA ADS] [CrossRef] [Google Scholar]
  3. Ait Benkhali, F., Hofmann, W., Rieger, F. M., & Chakraborty, N. 2020, A&A, 634, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Ajello, M., Angioni, R., Axelsson, M., et al. 2020, ApJ, 892, 105 [NASA ADS] [CrossRef] [Google Scholar]
  5. Alston, W. N., Parker, M. L., Markevičiūtė, J., et al. 2015, MNRAS, 449, 467 [CrossRef] [Google Scholar]
  6. Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071 [CrossRef] [Google Scholar]
  7. Auchère, F., Froment, C., Bocchialini, K., Buchlin, E., & Solomon, J. 2016, ApJ, 825, 110 [Google Scholar]
  8. Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Nature, 287, 307 [Google Scholar]
  9. Bhatta, G. 2019, MNRAS, 487, 3990 [NASA ADS] [CrossRef] [Google Scholar]
  10. Bhatta, G., & Dhital, N. 2020, ApJ, 891, 120 [NASA ADS] [CrossRef] [Google Scholar]
  11. Blandford, R., Meier, D., & Readhead, A. 2019, ARA&A, 57, 467 [NASA ADS] [CrossRef] [Google Scholar]
  12. Cerruti, M. 2020, Galaxies, 8, 72 [NASA ADS] [CrossRef] [Google Scholar]
  13. Connolly, S. D. 2015, arXiv e-prints [arXiv:1503.06676] [Google Scholar]
  14. Covino, S., Sandrinelli, A., & Treves, A. 2019, MNRAS, 482, 1270 [NASA ADS] [CrossRef] [Google Scholar]
  15. Covino, S., Landoni, M., Sandrinelli, A., & Treves, A. 2020, ApJ, 895, 122 [NASA ADS] [CrossRef] [Google Scholar]
  16. Emmanoulopoulos, D., McHardy, I. M., & Papadakis, I. E. 2013, MNRAS, 433, 907 [NASA ADS] [CrossRef] [Google Scholar]
  17. Gierliński, M., Middleton, M., Ward, M., & Done, C. 2008, Nature, 455, 369 [CrossRef] [Google Scholar]
  18. Gupta, A. C., Tripathi, A., Wiita, P. J., et al. 2019, MNRAS, 484, 5785 [NASA ADS] [Google Scholar]
  19. Honma, F., Matsumoto, R., & Kato, S. 1992, PASJ, 44, 529 [NASA ADS] [Google Scholar]
  20. Ingram, A. R., & Motta, S. E. 2019, New Astron. Rev., 85, 101524 [Google Scholar]
  21. Kushwaha, P., Sarkar, A., Gupta, A. C., Tripathi, A., & Wiita, P. J. 2020, MNRAS, 499, 653 [NASA ADS] [CrossRef] [Google Scholar]
  22. Li, H. Z., Jiang, Y. G., Yi, T. F., et al. 2018, Ap&SS, 363, 45 [NASA ADS] [CrossRef] [Google Scholar]
  23. McKinney, J. C., Tchekhovskoy, A., & Blandford, R. D. 2012, MNRAS, 423, 3083 [Google Scholar]
  24. Peñil, P., Domínguez, A., Buson, S., et al. 2020, ApJ, 896, 134 [CrossRef] [Google Scholar]
  25. Percival, D. P. 1995, Biometrika, 82, 619 [CrossRef] [Google Scholar]
  26. Prokhorov, D. A., & Moraghan, A. 2017, MNRAS, 471, 3036 [NASA ADS] [CrossRef] [Google Scholar]
  27. Rieger, F. M. 2004, ApJ, 615, L5 [NASA ADS] [CrossRef] [Google Scholar]
  28. Sandrinelli, A., Covino, S., & Treves, A. 2016, ApJ, 820, 20 [NASA ADS] [CrossRef] [Google Scholar]
  29. Sandrinelli, A., Covino, S., Treves, A., et al. 2017, A&A, 600, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Sandrinelli, A., Covino, S., Treves, A., et al. 2018, A&A, 615, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Sarkar, A., Gupta, A. C., Chitnis, V. R., & Wiita, P. J. 2021, MNRAS, 501, 50 [Google Scholar]
  32. Sillanpaa, A., Haarala, S., Valtonen, M. J., Sundelius, B., & Byrd, G. G. 1988, ApJ, 325, 628 [NASA ADS] [CrossRef] [Google Scholar]
  33. Sobacchi, E., Sormani, M. C., & Stamerra, A. 2017, MNRAS, 465, 161 [NASA ADS] [CrossRef] [Google Scholar]
  34. Tarnopolski, M., Żywucka, N., Marchenko, V., & Pascual-Granado, J. 2020, ApJS, 250, 1 [NASA ADS] [CrossRef] [Google Scholar]
  35. Tavani, M., Cavaliere, A., Munar-Adrover, P., & Argan, A. 2018, ApJ, 854, 11 [NASA ADS] [CrossRef] [Google Scholar]
  36. Timmer, J., & Koenig, M. 1995, A&A, 300, 707 [NASA ADS] [Google Scholar]
  37. Torrence, C., & Compo, G. P. 1998, Bull. Am. Meteorol. Soc., 79, 61 [Google Scholar]
  38. Urry, C. M., & Padovani, P. 1995, PASP, 107, 803 [NASA ADS] [CrossRef] [Google Scholar]
  39. Valtonen, M. J., Nilsson, K., Villforth, C., et al. 2009, ApJ, 698, 781 [NASA ADS] [CrossRef] [Google Scholar]
  40. Wang, G. G., Cai, J. T., & Fan, J. H. 2022, ApJ, 929, 130 [NASA ADS] [CrossRef] [Google Scholar]
  41. Yang, J., Cao, G., Zhou, B., & Qin, L. 2021a, PASP, 133, 024101 [CrossRef] [Google Scholar]
  42. Yang, S., Yan, D., Zhang, P., Dai, B., & Zhang, L. 2021b, ApJ, 907, 105 [NASA ADS] [CrossRef] [Google Scholar]
  43. Zhang, P.-F., Yan, D.-H., Liao, N.-H., & Wang, J.-C. 2017a, ApJ, 835, 260 [NASA ADS] [CrossRef] [Google Scholar]
  44. Zhang, P.-F., Yan, D.-H., Zhou, J.-N., et al. 2017b, ApJ, 845, 82 [NASA ADS] [CrossRef] [Google Scholar]
  45. Zhang, P.-F., Yan, D.-H., Zhou, J.-N., Wang, J.-C., & Zhang, L. 2020, ApJ, 891, 163 [NASA ADS] [CrossRef] [Google Scholar]
  46. Zhang, H., Yan, D., Zhang, P., Yang, S., & Zhang, L. 2021, ApJ, 919, 58 [NASA ADS] [CrossRef] [Google Scholar]
  47. Zhou, J., Wang, Z., Chen, L., et al. 2018, Nat. Commun., 9, 4599 [NASA ADS] [CrossRef] [Google Scholar]
  48. Żywucka, N., Tarnopolski, M., Marchenko, V., & Pascual-Granado, J. 2022, High Energy Astrophysics in Southern Africa 2021, online at https://pos.sissa.it/cgi-bin/reader/conf.cgi?confid=401, 8 [Google Scholar]

Appendix A: Wavelet analysis method

Figure A.1 shows the Morlet wavelet used as the mother wavelet in our calculations, with the Fourier transform of the wavelet on the right. Figure A.2 shows the CWT map and global wavelet spectrum of a periodic function with two periods, 0.1 s and 0.02 s. A random normally distributed error of ∼10% to the amplitude of the function is added to simulate the error in the light curve data. The maximum values of the global wavelet have a relative error of ∼1% to ∼2%. On the other hand, the broadness of the humps are much larger than the relative error on the maximum, being ∼20% to ∼40% of the real value. Thus, the main contribution to the period uncertainty is the intrinsic broadness of the features.

thumbnail Fig. A.1.

The mother wavelet used in our analysis. a) Real part of the Morlet wavelet in solid line and the imaginary part in dashed line. b) Fourier transform of the Morlet wavelet.

The dashed line in Fig. A.2 (c) is the 95% significance level provided by PyCWT, which in was not used for the light curves of Fermi-LAT sources in this work. Instead, we performed a Monte Carlo approach to take into account the statistical and variability properties of each light curve, as explained in Appendix B and Appendix C. The cross-hatched area in the CWT map shown in Fig. A.2 (b) corresponds to the cone of influence (COI). The importance of the COI is evident in this figure, where we can visualize the spurious signal at the upper edge and the first periodicity is clearly distorted when it is within the COI.

thumbnail Fig. A.2.

CWT map for a periodic function. Each panel represents: a) periodic function with two periods p1 = 0.02 s and p2 = 0.10 s, b) corresponding wavelet power spectrum map obtained using PyCWTand c) global power spectrum shown as the solid line, with a significance level at 95%, provided by PyCWT as a dashed line. The dashed areas in b) indicate the COI.

A second example is shown in Fig. A.3, where in this case, we can see how the CWT behaves when the input is a Dirac delta function. The vertical shape of the wavelet power spectrum in Fig. A.3 (b), with increasing time interval as the scale increases, reveals the response of CWT to flare-shaped signals.

thumbnail Fig. A.3.

CWT map for a Dirac Delta function. Each panel represents: a) Dirac delta function, b) corresponding wavelet power spectrum map obtained using PyCWT and c) global power spectrum in solid line and the significance level at 95% provided by PyCWT in dashed line. The dashed area in b) indicate the COI.

Appendix B: Statistical and variability properties of artificial light curves

Artificial light curves computed using the Python code provided by Connolly (2015), following the work of Emmanoulopoulos et al. (2013) share both the PDF and the PSD of the original light curve, as can be seen in Fig. B.1.

thumbnail Fig. B.1.

Comparison of the PDF and the PSD of the original light curve of the source S5 0716+714 and the light curve simulated by Connolly (2015), in time bins of 30 (upper) and 7 days (bottom). The solid line in the PSD is a simple power-law fit for a visual check. The actual fitted PSD is a smoothly bending power-law model plus a constant (see Eq. 2 in Emmanoulopoulos et al. 2013).

Appendix C: Significance levels estimation

We fit the global power spectrum histograms, computed using the artificial light curves, for each period with a χ2 function as shown in example histograms in Fig. C.1. We also display the five σ confidence levels, equivalent to those of a Gaussian distribution, namely, 68.3%, 95.4%, 99.7%, 99.994%, and 99.99994%, for both the pre-trial (solid lines) and post-trial (dashed lines) cases.

thumbnail Fig. C.1.

Global power spectrum histogram fitted with a χ2 function. The upper subfigure is the histogram for the monthly binned light curve of S5 0716+714, at the scale corresponding to period ∼197 d, and the bottom subfigure is the histogram for the weekly binned light curve, at the scale corresponding to period ∼260 d. Five pre-trial confidence levels are shown in both graphs corresponding to 1 to 5σ significance in solid lines, whereas the post-trial confidence levels are shown in dashed lines.

Appendix D: CWT for all sources

thumbnail Fig. D.1.

CWT map for monthly binned light curve (left) and weekly binned light curve (right) of 3C 66A and 4C +28.07 and fitted light curves for 4C +01.02.

thumbnail Fig. D.2.

CWT map for monthly binned light curve (left) and weekly binned light curve (right) of PKS 0235+164 and NGC 1275 and fitted light curves for NGC 1275.

thumbnail Fig. D.3.

CWT map for monthly binned light curve (left) and weekly binned light curve (right) of PKS 0402-362 and PKS 0426-380 and the respective fitted light curves.

thumbnail Fig. D.4.

CWT map for monthly binned light curve (left) and weekly binned light curve (right) of PKS 0447-439 and PKS 0454-234 and the respective fitted light curves.

thumbnail Fig. D.5.

CWT map for monthly binned light curve (left) and weekly binned light curve (right) of S5 0716+714 and 4C +55.17 and the fitted light curves for S5 0716+714.

thumbnail Fig. D.6.

CWT map for monthly binned light curve (left) and weekly binned light curve (right) of 1H 1013+498 and the fitted light curves.

thumbnail Fig. D.7.

CWT map for monthly binned light curve (left) and weekly binned light curve (right) of Mrk 421 and 4C +21.35, and the respective fitted light curves.

thumbnail Fig. D.8.

CWT map for monthly binned light curve (left) and weekly binned light curve (right) of 3C 273 and PKS 1424+240, and the fitted light curves for 3C 273.

thumbnail Fig. D.9.

CWT map for monthly binned light curve (left) and weekly binned light curve (right) of B3 1343+451 and 3C 279, and the fitted light curves for 3C 279. The fitted light curve with period ∼101 d (third subfigure on the right hand side column) shows six cycles between around MJD 57700 and MJD 58400. However, we can fit six more cycles if we consider the previous rise in power spectrum during MJD 56500 and MJD 57300, approximately.

thumbnail Fig. D.10.

CWT map for monthly binned light curve (left) and weekly binned light curve (right) of PKS 1424-418 and PKS 1502+106 and the fitted light curves for PKS 1424-418.

thumbnail Fig. D.11.

CWT map for monthly binned light curve (left) and weekly binned light curve (right) of PKS 1510-089 and PG 1553+113 and the fitted light curve for PKS 1510-089.

thumbnail Fig. D.12.

CWT map for monthly binned light curve (left) and weekly binned light curve (right) of 4C +38.41 and Mrk 501 and the fitted light curves for Mrk 501.

thumbnail Fig. D.13.

CWT map for monthly binned light curve (left) and weekly binned light curve (right) of 1E 1959+650 and PKS 2155-304 and the fitted light curves for PKS 2155-304.

thumbnail Fig. D.14.

CWT map for monthly binned light curve (left) and weekly binned light curve (right) of BL Lac and CTA 102 and the fitted light curves for CTA 102.

thumbnail Fig. D.15.

CWT map for monthly binned light curve (left) and weekly binned light curve (right) of 3C 454.3 and PKS 2326-502 and the fitted light curves for 3C 454.3.

thumbnail Fig. D.16.

CWT map for monthly binned light curve (left) and weekly binned light curve (right) of PMN J2345-1555 and the fitted light curves.

All Tables

Table 1.

Selected Fermi-LAT sources with the catalog name, blazar type, coordinate in J2000, and redshift extracted from the Fermi-LAT 4FGL catalog.

In the text
Table 2.

QPO candidates identified by the CWT of the Fermi-LAT light curves in time bins of 30 days and 7 days.

In the text

All Figures

thumbnail Fig. 1.

CWT map for the monthly binned (upper) and the weekly binned light curve (bottom) for 4C +01.02 (left). In each subplot, the panels represent (a) Fermi-LAT light curve, (b) wavelet power spectrum and (c) global wavelet power spectrum. The solid coloured contours in (b) and the dashed coloured lines in (c) are the confidence levels (1 to 5σ in black, blue, green, violet, and red). Monthly binned and weekly binned light curves with the fitted periodic signal in red dashed lines in the upper panels and the relative error in the bottom panels (right).
In the text
thumbnail Fig. 2.

CWT map for the monthly binned (upper) and the weekly binned light curve (bottom) for PKS 0537−441 (left). In each subplot, the panels represent (a) Fermi-LAT light curve, (b) wavelet power spectrum and (c) global wavelet power spectrum. The solid coloured contours in (b) and the dashed coloured lines in (c) are the confidence levels (1 to 5σ in black, blue, green, violet, and red). Monthly binned and weekly binned light curves with the fitted periodic signal in red dashed lines in the upper panels and the relative error in the bottom panels (right).
In the text
thumbnail Fig. 3.

CWT map for the monthly binned (upper) and the weekly binned light curve (bottom) for S5 1044+71 (left). In each subplot, the panels represent (a) Fermi-LAT light curve, (b) wavelet power spectrum and (c) global wavelet power spectrum. The solid coloured contours in (b) and the dashed coloured lines in (c) are the confidence levels (1 to 5σ in black, blue, green, violet, and red). Monthly binned and weekly binned light curves with the fitted periodic signal in red dashed lines in the upper panels and the relative error in the bottom panels (right).
In the text
thumbnail Fig. 4.

CWT map for the monthly binned (upper) and the weekly binned light curve (middle) for B2 1520+31 (left). In each subplot, the panels represent (a) Fermi-LAT light curve, (b) wavelet power spectrum and (c) global wavelet power spectrum. The solid coloured contours in (b) and the dashed coloured lines in (c) are the confidence levels (1 to 5σ in black, blue, green, violet, and red). Monthly binned and weekly binned light curves with the fitted periodic signal in red dashed lines in the upper panels and the relative error in the bottom panels (right). The third CWT is computed with the weekly binned light curve cut from MJD 54697 to MJD 56100. And on the right, the respective fitted light curve shows the ∼39 d period.
In the text
thumbnail Fig. 5.

CWT map for the monthly binned (upper) and the weekly binned light curve (middle) for PKS 2247−131 (left). In each subplot, the panels represent (a) Fermi-LAT light curve, (b) wavelet power spectrum and (c) global wavelet power spectrum. The solid coloured contours in (b) and the dashed coloured lines in (c) are the confidence levels (1 to 5σ in black, blue, green, violet, and red). Monthly binned and weekly binned light curves with the fitted periodic signal in red dashed lines in the upper panels and the relative error in the bottom panels (right). In the bottom fitted light curve we show a possible reappearance of the candidate QPO with ∼34 d period at MJD 58050, adding three more cycles to the series.
In the text
thumbnail Fig. A.1.

The mother wavelet used in our analysis. a) Real part of the Morlet wavelet in solid line and the imaginary part in dashed line. b) Fourier transform of the Morlet wavelet.
In the text
thumbnail Fig. A.2.

CWT map for a periodic function. Each panel represents: a) periodic function with two periods p1 = 0.02 s and p2 = 0.10 s, b) corresponding wavelet power spectrum map obtained using PyCWTand c) global power spectrum shown as the solid line, with a significance level at 95%, provided by PyCWT as a dashed line. The dashed areas in b) indicate the COI.
In the text
thumbnail Fig. A.3.

CWT map for a Dirac Delta function. Each panel represents: a) Dirac delta function, b) corresponding wavelet power spectrum map obtained using PyCWT and c) global power spectrum in solid line and the significance level at 95% provided by PyCWT in dashed line. The dashed area in b) indicate the COI.
In the text
thumbnail Fig. B.1.

Comparison of the PDF and the PSD of the original light curve of the source S5 0716+714 and the light curve simulated by Connolly (2015), in time bins of 30 (upper) and 7 days (bottom). The solid line in the PSD is a simple power-law fit for a visual check. The actual fitted PSD is a smoothly bending power-law model plus a constant (see Eq. 2 in Emmanoulopoulos et al. 2013).
In the text
thumbnail Fig. C.1.

Global power spectrum histogram fitted with a χ2 function. The upper subfigure is the histogram for the monthly binned light curve of S5 0716+714, at the scale corresponding to period ∼197 d, and the bottom subfigure is the histogram for the weekly binned light curve, at the scale corresponding to period ∼260 d. Five pre-trial confidence levels are shown in both graphs corresponding to 1 to 5σ significance in solid lines, whereas the post-trial confidence levels are shown in dashed lines.
In the text
thumbnail Fig. D.1.

CWT map for monthly binned light curve (left) and weekly binned light curve (right) of 3C 66A and 4C +28.07 and fitted light curves for 4C +01.02.
In the text
thumbnail Fig. D.2.

CWT map for monthly binned light curve (left) and weekly binned light curve (right) of PKS 0235+164 and NGC 1275 and fitted light curves for NGC 1275.
In the text
thumbnail Fig. D.3.

CWT map for monthly binned light curve (left) and weekly binned light curve (right) of PKS 0402-362 and PKS 0426-380 and the respective fitted light curves.
In the text
thumbnail Fig. D.4.

CWT map for monthly binned light curve (left) and weekly binned light curve (right) of PKS 0447-439 and PKS 0454-234 and the respective fitted light curves.
In the text
thumbnail Fig. D.5.

CWT map for monthly binned light curve (left) and weekly binned light curve (right) of S5 0716+714 and 4C +55.17 and the fitted light curves for S5 0716+714.
In the text
thumbnail Fig. D.6.

CWT map for monthly binned light curve (left) and weekly binned light curve (right) of 1H 1013+498 and the fitted light curves.
In the text
thumbnail Fig. D.7.

CWT map for monthly binned light curve (left) and weekly binned light curve (right) of Mrk 421 and 4C +21.35, and the respective fitted light curves.
In the text
thumbnail Fig. D.8.

CWT map for monthly binned light curve (left) and weekly binned light curve (right) of 3C 273 and PKS 1424+240, and the fitted light curves for 3C 273.
In the text
thumbnail Fig. D.9.

CWT map for monthly binned light curve (left) and weekly binned light curve (right) of B3 1343+451 and 3C 279, and the fitted light curves for 3C 279. The fitted light curve with period ∼101 d (third subfigure on the right hand side column) shows six cycles between around MJD 57700 and MJD 58400. However, we can fit six more cycles if we consider the previous rise in power spectrum during MJD 56500 and MJD 57300, approximately.
In the text
thumbnail Fig. D.10.

CWT map for monthly binned light curve (left) and weekly binned light curve (right) of PKS 1424-418 and PKS 1502+106 and the fitted light curves for PKS 1424-418.
In the text
thumbnail Fig. D.11.

CWT map for monthly binned light curve (left) and weekly binned light curve (right) of PKS 1510-089 and PG 1553+113 and the fitted light curve for PKS 1510-089.
In the text
thumbnail Fig. D.12.

CWT map for monthly binned light curve (left) and weekly binned light curve (right) of 4C +38.41 and Mrk 501 and the fitted light curves for Mrk 501.
In the text
thumbnail Fig. D.13.

CWT map for monthly binned light curve (left) and weekly binned light curve (right) of 1E 1959+650 and PKS 2155-304 and the fitted light curves for PKS 2155-304.
In the text
thumbnail Fig. D.14.

CWT map for monthly binned light curve (left) and weekly binned light curve (right) of BL Lac and CTA 102 and the fitted light curves for CTA 102.
In the text
thumbnail Fig. D.15.

CWT map for monthly binned light curve (left) and weekly binned light curve (right) of 3C 454.3 and PKS 2326-502 and the fitted light curves for 3C 454.3.
In the text
thumbnail Fig. D.16.

CWT map for monthly binned light curve (left) and weekly binned light curve (right) of PMN J2345-1555 and the fitted light curves.
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.