Possible Quasi-periodic Modulation in the z = 1.1 Gamma-Ray Blazar PKS 0426–380

Blazars are one type of active galactic nuclei (AGNs) that aim their jets almost directly at Earth. Blazar emission is generally dominated by nonthermal radiation over all frequencies ranging from radio to TeV γ-rays. The typical multiwavelength spectral energy distribution of a blazar is characterized by two distinct bumps. It is generally accepted that the first bump peaking in infrared to X-ray frequencies is synchrotron emission from relativistic electrons in the jet. The second bump peaking in the MeV to GeV band could be produced via inverse Compton (IC) scattering of synchrotron photons (i.e., synchrotron self-Compton [SSC]; e.g., Maraschi et al. 1992; Tavecchio et al. 1998; Finke et al. 2008; Yan et al. 2014) and external photons (i.e., external Compton [EC]; e.g., Dermer & Schlickeiser 1993; Sikora et al. 1994; Kang et al. 2014) by the same population of relativistic electrons that produce the synchrotron emission.

According to the emission lines' features, blazars are classified as BL Lacertae objects (BL Lacs; having weak or no emission lines) and flat-spectrum radio quasars (FSRQs; having strong emission lines). FSRQs are usually the low synchrotron-peaked blazars (i.e., the synchrotron peak frequency ${\nu }_{{\rm{s}}}\lt {10}^{14}$ Hz). There exist low, intermediate, and high synchrotron-peaked BL Lacs (LBLs, IBLs, and HBLs, respectively, defined by whether ${\nu }_{{\rm{s}}}\lt {10}^{14}$ Hz, ${10}^{14}\,\lt {\nu }_{{\rm{s}}}\ (\mathrm{Hz})\lt {10}^{15}$, or ${\nu }_{{\rm{s}}}\gt {10}^{15}$ Hz; Abdo et al. 2010). γ-rays from HBLs can be produced by SSC, while γ-rays from FSRQs are usually attributed to EC. It is very likely that the jet properties of LBLs are similar to that of FSRQs.

By modeling blazar spectra, emission mechanisms and physical properties of the relativistic jets can be determined (e.g., Zhang et al. 2012; Ghisellini et al. 2014). Besides, the information of violent variation, including variability timescale and the profile of light curves, also puts constraints on the jet properties of blazars. Although it seems that blazar variability is usually aperiodic, many efforts have still been made to search for the periodic variabilities in blazars (e.g., Kidger et al. 1992; Bai et al. 1998, 1999; Fan & Lin 2000; Xie et al. 2008; Li et al. 2009; Urry 2011; King et al. 2013; Gupta 2014). Naturally, these studies mainly focus on radio and optical–UV bands because of the abundant data.

Quasi-periodic oscillations (QPOs) have been detected in the X-ray radiation of black hole (BH) X-ray binaries (Remillard et al. 2002). The X-ray QPOs were found to have two types, low-frequency QPOs and high-frequency QPOs (Remillard & McClintock 2006), which are thought to originate in the inner accretion disk of a BH. Remillard & McClintock (2006) suggested an inverse linear relation between QPO frequency and BH mass. Recently Pan et al. (2016) and Zhang et al. (2017b) analyzed the X-ray data of narrow-line Seyfert 1 galaxy 1H 0707-495 and found a possible QPO. The QPO in 1H 0707-495 follows the QPO frequency–BH mass relation suggested by Remillard & McClintock (2006). It is interesting that this relation spans from stellar-mass to supermassive BHs. It may indicate that accretion onto a stellar-mass BH is comparable to accretion onto a supermassive BH.

In this work we focus on QPOs in blazars. It should be pointed out that blazar emission is dominated by nonthermal radiation from jets, and the X-ray QPOs mentioned above were found in thermal radiation from accretion disks. Some possible QPO signals were found in radio, optical, and X-ray bands. The period has two types, the short term of several tens of days and the long term of several years (e.g., Rieger 2007). QPOs in blazars may give us more insight into the physics of blazars.

Thanks to the Large Area Telescope (LAT) on board Fermi (Abdo et al. 2009), running in an all-sky coverage monitoring every ∼3 hr (Abdo et al. 2009; Atwood et al. 2009), we have a long-term view of the variability of a large sample of γ-ray blazars. Fermi-LAT has been collecting data for over 8 yr, allowing us to search for quasi-periodic variability with a timescale of a few years in gamma-ray flux (e.g., Sandrinelli et al. 2014, 2016; Ackermann et al. 2015; Zhang et al. 2017a).

So far, three blazars (PG 1553+113, PKS 2155–304, and PKS 0537–441) have been reported having quasi-periodic variability in gamma-ray fluxes with the significance of $\geqslant 3\ \sigma$ (Sandrinelli et al. 2014, 2016; Ackermann et al. 2015; Zhang et al. 2017a). In particular, Zhang et al. (2017a) reported a γ-ray quasi-periodic variability in PKS 2155–304 with a significance of 4.9σ. The three blazars are HBLs. HBL usually has a "clean" environment around the jet, and the γ-rays are produced by SSC. It is interesting to search for quasi-periodic variabilities in other subclasses of blazars, e.g., FSRQs.

In this work we analyze the gamma-ray data from the distant FSRQ PKS 0426–380 in the interval between 2008 August and 2016 December. We find a significant quasi-periodicity in the γ-ray flux. The significance of the 3.35 ± 0.68 yr period is $\simeq 3.6\sigma$. This paper is organized as follows: we describe Fermi-LAT data analysis procedures in Section 2, in Section 3.1 we show results, and summaries and discussions are given in Section 4.

The blazar PKS 0426–380 is classified as an FSRQ (Healey et al. 2007; Ghisellini et al. 2011; Sbarrato et al. 2012). PKS 0426–380 has broad lines visible in the low emission state (Sbarufatti et al. 2005). PKS 0426–380 has a large γ-ray luminosity (i.e., ${L}_{\gamma }\gt {10}^{48}\ \mathrm{erg}\ {{\rm{s}}}^{-1}$) and a large ratio between the broad-line region (BLR) luminosity and the Eddington luminosity (i.e., ${L}_{\mathrm{BLR}}/{L}_{\mathrm{Edd}}\sim {10}^{-3}$; Ghisellini et al. 2011). PKS 0426–380 is a distant blazar with the redshift z = 1.11 (Heidt et al. 2004; Sbarufatti et al. 2005).

The public γ-ray data of PKS 0426–380 are obtained from the observations of the LAT. Fermi-LAT is an electron–positron pair production detection sensitive to photon energies from ∼20 MeV to 300 GeV.⁶ The LAT has a largest effective area of $\sim 8000\,{\mathrm{cm}}^{2}$ at 1 GeV, a field of view of ∼2.4 sr, and a point-spread function of $\lt 0\buildrel{\circ}\over{.} 8$ above 1 GeV.

The data are selected between 2008 August 4 and 2016 December 15 (MJD 54,682.66–57,737.66), with an energy range from 100 MeV to 500 GeV. The LAT data analysis employs the Fermi Science Tools version v10r0p5 package. We follow the standard procedure provided in the Fermi Science Support Center (FSSC),⁷ to reduce the data. We select the events in a circular region of interest (ROI) of 15° radius centered on the position of PKS 0426–380. To minimize the contamination of γ-rays from the Earth limb, we exclude the events with zenith angle $\gt \,90^\circ$. The diffuse γ-rays emission of the Galactic and extragalactic are modeled using two files: gll_iem_v06.fit and iso_P8R2_SOURCE_V6_v06.txt. We run the Science Tools gtmktime to select the good time intervals. The SOURCE class of photon-like events of the new Pass 8 data (Atwood et al. 2013) is used with the instrument response functions (IRFs) P8R2_SOURCE_V6.

A binned maximum likelihood algorithm implemented in gtlike is applied between the model file and the events of observation. The data from a $20^\circ \times 20^\circ$ square ROI are divided into a spatial pixel size of $0\buildrel{\circ}\over{.} 1\times 0\buildrel{\circ}\over{.} 1$. In energy dimensionality, they are binned into 30 logarithmically equal bins. We use the script make3FGLxml.py to generate the model file.⁸ The file is composed of the information (the parameters of the spectrum and the spatial positions) of all known 3FGL sources (Acero et al. 2015). We employ gtlike to derive the best-fitting results, including flux, photon index, and test statistic (TS) value. The TS is defined as $\mathrm{TS}=-2\mathrm{ln}({{ \mathcal L }}_{0}/{{ \mathcal L }}_{1})$, where ${{ \mathcal L }}_{0}$ and ${{ \mathcal L }}_{1}$ are, respectively, the maximum likelihood value of the model with and without the source, and it describes the significance of the source. We save the best-fitting results into a new model file. The following analysis of light curves is based on this new model file. The results in this paper are only given with statistical errors.

3.1.  $\gamma$-ray Light Curves

Fermi-LAT collaboration provides two techniques to generate light curves, i.e., the maximum likelihood optimization and the exposure-weighted aperture photometry (Corbet et al. 2007; Kerr 2011).

In the analysis with the maximum likelihood optimization, the light curves are first built by using a 1-month bin. We employ an unbinned maximum likelihood fitting technique and run the Science Tools gtlike to obtain the flux and TS value for each time bin. In this step, the events are selected in a circular ROI of 15° radius centered on the coordinates of the target. The model is the same as the new model file mentioned in Section 2, except the parameters of the spectrum shape are frozen. The light curve is shown in Figure 1. To test the impacts of different time bins on the results, we also build the light curve with a 10-day bin over 0.1 GeV (Figure 2).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the analysis with the exposure-weighted aperture photometry, we generate light curves with a modified version of aperture photometry. We exclude the events during the period when the target is within 5° of the Sun. The probability of each photon from a specific source is calculated by the Science Tools gtpsrcrob with the best-fitting model and IRFs P8R2_SOURCE_V6. Then we sum the probability of the events within a 3° radius circle region centered on the position of PKS 0426–380 over 100 MeV. The 2.5 day bin is used. Then the exposure of each time bin is determined with the Science Tools gtexposure. The light curve weighted by the exposure of each time bin is shown in Figure 3.

Zoom In Zoom Out Reset image size

3.2. Searching for $\gamma$-ray Quasi-periodic Modulation

The Lomb–Scargle periodogram (LSP; Lomb 1976; Scarle 1982) is a widely used tool for searching for period of variability. However, this method does not take into account the fluctuations of periodic signal with time. Therefore, we also use the weighted wavelet Z-transform (WWZ) technique with a Morlet mother function (Foster 1996) for testing the result of the LSP. In Figures 1–3 we also show the corresponding LSP and WWZ powers. All the powers show an obvious peak at ∼1222 days (i.e., 3.35 ± 0.68 yr). The uncertainty of the signal is evaluated based on the half-width at half-maximum of the Gaussian fitting at the position of the power peak. We first apply a simple method for testing the significance of peaks in the periodogram of red-noise data (Vaughan 2005) to all the highest power peaks in Figures 1–3, and we derive a similar significance of $\sim 3.5\ \sigma$. For obtaining the precise significance of the signal, we should simulate light curves based on the obtained best-fitting result of power spectral density and the probability density function of observed variation (e.g., Zhang et al. 2017a). The details of the simulation and significance estimation are provided by Emmanoulopoulos et al. (2013), Ackermann et al. (2015), and Bhatta et al. (2016). Following the procedure, we simulate 3 × 10⁴ light curves with the DELCgen program, and then we evaluate the significance of the signal. The significance of the signal is $\simeq 3.6\,\sigma$.

Photons above 100 GeV from PKS 0426–380 are found in our analysis (also see Tanaka et al. 2013). In order to avoid the effect of the absorption by extragalactic background light, we build the light curves in the energy range of 1–50 GeV (Figure 4). Again we find a strong signal near the period of ∼1222 days in the corresponding LSP and WWZ powers.

Zoom In Zoom Out Reset image size

Using the phase-resolved likelihood analysis method, we fold the gamma-ray photons from the $20^\circ \times 20^\circ$ square region centered on the position of the target. We calculate a phase for every event in the region according to the period and the arrival time (with phase zero corresponding to MJD 54,682.66). The pulse interval is divided into 15 segments. For each segment, we use a binned likelihood technique to obtain the flux with the same model file as the generating light curve. The folded light curve is shown in Figure 5. We fit the folded light curve with a constant flux and derive the reduced ${\chi }_{\min }^{2}=3156.5/14$. The null hypothesis that the phase-resolved light curve is steady can be rejected at a far-exceeding 5σ confidence level. The folded light curve in Figure 5 also confirms the signal mentioned above and indicates that the gamma-ray flux varies with phase.

Zoom In Zoom Out Reset image size

3.3. Searching for Optical–UV Quasi-periodic Modulation

We have tried to collect long-term multiwavelength archive data of PKS 0426–380. We obtain the optical–UV data covering from 2005 August to 2013 April provided by the ASI Science Data Center (see Figure 6). The optical–UV flux is analyzed using the same techniques as for the γ-ray flux but does not show any significant periodic signal.

Zoom In Zoom Out Reset image size

We have produced the γ-ray light curves of distant FSRQ PKS 0426–380 with the latest LAT data. We have found a significant quasi-periodic variability in the γ-ray flux, with a period of 3.35 ± 0.68 yr and a significance of $3.6\sigma$. More multiwavelength data are needed to confirm the periodicity in PKS 0426–380. If the γ-ray quasi-periodic variability is true, the next γ-ray flux peak will occur in 2019.

The archive optical–UV data covering from 2005 August to 2013 April have been also analyzed, but no quasi-periodic signal is found. This may be due to the incompleteness of the data. On the other band, optical–UV emission from FSRQs may be dominated by accretion disk photons. This complicates the confirmations of γ-ray quasi-periodic variability in FSRQs.

Gamma-ray quasi-periodic modulations with a significance of $\geqslant 3\sigma$ have been reported in HBLs (Sandrinelli et al. 2014, 2016; Ackermann et al. 2015; Zhang et al. 2017a). PKS 0426–380 is the first FSRQ having a significant γ-ray quasi-periodic variability. Although PKS 0426–380 has a longer observed period (${T}_{\mathrm{obs}}\sim 3$ yr) than the three HBLs (${T}_{\mathrm{obs}}\sim 2$ yr), the intrinsic periods (${T}_{\mathrm{int}}$) for PKS 0426–380 and the three HBLs are the same because of ${T}_{\mathrm{obs}}={T}_{\mathrm{int}}(1+z)$. It is well known that FSRQs have a larger BH mass and a higher accretion rate compared to HBLs (e.g., Ghisellini et al. 2011, 2014). It is interesting that FSRQs and HBLs have the same period. It may indicate that the γ-ray quasi-periodic modulations in blazars are not solely linked to the central BH and the accretion disk.

A possible interpretation for the γ-ray QPOs is a varying Doppler boosting for a periodically changing viewing angle, which is related to a helical jet (e.g., Rieger 2004; Komossa & Zensus 2016). The helical structure in the radio jet of AGNs has been detected (e.g., Rieger 2007). One way to explain the helical jet is to involve the presence of a binary supermassive BH (SMBH; e.g., Komossa & Zensus 2016; Sobacchi et al. 2017). In the scenario of binary SMBHs, one can estimate the parameters of the binary system, e.g., the separation between the two SMBHs (e.g., Sobacchi et al. 2017),

$$\begin{eqnarray*}&&D\sim {10}^{16}\displaystyle \frac{1+q}{q}\displaystyle \frac{{T}_{\mathrm{int}}}{2\ \mathrm{yr}}\ \mathrm{cm},\end{eqnarray*}$$

where q is the mass ratio of the secondary SMBH to the primary SMBH. Assuming $q\sim 1$ and using T_int ∼ 2 yr, we obtain D ∼ 0.003 pc. The two SMBHs will merge owing to emission of gravitational waves. The timescale for SMBHs merging is (e.g., Sobacchi et al. 2017)

$$\begin{eqnarray*}&&{t}_{{\rm{m}}}\sim 4\times {10}^{4}q{\left(\displaystyle \frac{q}{1+q}\right)}^{3}\displaystyle \frac{{T}_{\mathrm{int}}}{2\ \mathrm{yr}}\ \mathrm{yr}.\end{eqnarray*}$$

Using $q\sim 1$ and ${T}_{\mathrm{int}}\sim 2\mathrm{yr}$, we obtain ${t}_{{\rm{m}}}\sim 4000\ \mathrm{yr}$.

It should be kept in mind that the binary SMBH scenario is not the only explanation for the γ-ray QPOs. Alternative models, for instance, disk oscillations or Kelvin–Helmholtz instabilities, have been proposed (see Komossa & Zensus 2016, for a review). Nevertheless, we suggest that cumulative multifrequency blazar QPOs can be used to probe the underlying physical mechanism, and they may also shed light on the questions of galaxy mergers. Moreover, multifrequency QPOs may put new constraints on blazar emission models (Zhang et al. 2017a).

We thank the referees for constructive suggestions that significantly improved the manuscript. Part of this work is based on archival data, software, or online services provided by the ASI Science Data Center (ASDC). We acknowledge the financial support from the 973 Program of China under grant 2013CB837000, the National Natural Science Foundation of China (NSFC-11573060 and NSFC-11573026), and the Key Laboratory of Astroparticle Physics of Yunnan Province (no. 2016DG006). The work of D.-h.Y. is supported by the CAS "Light of West China" Program.

In addition to the methods based on Fourier analysis (e.g., LSP and WWZ), an autocorrelation analysis also has been used in light-curve analysis to search for quasi-periodic behavior or variations with a characteristic timescale (e.g., Li et al. 2016). A regularly periodic light curve will have strong correlations at time lags that are multiples of the period. For an evenly spaced time series, the autocorrelation function (ACF) is the fraction of the total variance due to correlated values at time lags. We apply an autocorrelation analysis based on the algorithm of Box & Jenkins (1976) to the γ-ray data. We calculate the ACF with the parametric autoregressive integrated moving average model. In Figure 7, we show the ACF of the 10-day-bin likelihood gamma-ray light curve. It can be seen that two positive correlation peaks appear respectively at the lag times of ∼1200 days and ∼2300 days. The result of autocorrelation analysis is consistent with the result of Fourier analysis.

Zoom In Zoom Out Reset image size