The X-ray outburst of PG 1553+113: A precession effect of two jets in the supermassive black hole binary system

The Large Area Telescope (Fermi-LAT), the primary instrument on the Fermi Gamma-ray Space Telescope (Fermi) mission, is an imaging, wide field-of-view, high-energy $\gamma$-ray telescope, covering the energy range from below 20 MeV to more than 300 GeV (Atwood et al., 2009). We collected $\gamma$-ray data of PG 1553+113 from the Fermi-LAT data website¹¹1https://fermi.gsfc.nasa.gov/cgi-bin/ssc/LAT/LATDataQuery.cgi obtained during 8 years from Jan. 1, 2012 to Jan. 1, 2020.

The data were analyzed with standard unbinned likelihood tutorials using the Fermi Science Tools Fermitools version 1.0.10. The photons event data were extracted by gtselect while selecting the source class events centered on the target from the region of interest (ROI) with the radius of $10^{\circ}$. The good time intervals were selected with the gtmktime procedure. A counts map of ROI was created by gtbin and the exposure map was generated by gtltcube and gtexpmap.

We have adopted the current Galactic diffuse emission model (gll_iem_v07.fits) and the model for the extragalactic isotropic diffuse emission iso_P8R3_SOURCE_V2_v1.txt. The source model XML file was created with Python code make4FGLxml.py, and the diffuse source responses were created by gtdiffrsp. Finally, the gtlike procedure was run to obtain flux. The threshold of the test statistics (TS) value is set to be 10 in this work.

XMM-Newton is an X-ray observatory, and its main instrument is the European Photon Imaging Camera (EPIC), consisting of two MOS detectors and a pn camera which operate in the 0.2–12.0 keV energy range. Totally, 24 sets of observations of PG 1553+113 obtained from 2001 to 2020 were employed and analysed with the software Science Analysis System (SAS version 18.0). We followed the recommended data analysis threads by XMM-Newton Science Operation Centre, cifbuild and odfingest procedures for the preparation, and then xmmextractor for the extractions of the light curves and spectra. The evselect was run when we extracted the light curves on 0.3–2.0 keV and 2.0–10.0 keV.

The 40 m telescope at the Owens Valley Radio Observatory (OVRO), is monitoring more than 1800 blazars (Richards et al., 2011). We obtain the 15 GHz data covering from 2012 January to 2020 January from the public data archived website⁵⁵5http://www.astro.caltech.edu/ovroblazars/index.php?page=home.

In this model, a SMBHB system is considered and each black hole carries its own jet and the movement of the black holes on the orbit causes the quasi-periodic variation on the light curves (Villata et al., 1998; Tavani et al., 2018).

Motivated by the work of Villata et al. (1998), we obtain a function of the angle between the line of sight and the jet direction

where the subscripts 1 and 2 correspond to the first and the second jet, respectively, and $\theta_{obs,1}$, $\theta_{obs,2}$ are the angles between the line of sight and the jet direction, $\phi_{1}$ and $\phi_{2}$ are the angles between the line of sight and the orbital angular momentum/spin of the individual black hole (assuming that the spin parallel to the orbital angular momentum), $\alpha_{1}$ and $\alpha_{2}$ are the angles between the orbital angular momentum and the jet direction, $\omega$ is the orbital angular velocity and $\psi$ is the phase position.

The flux of the jet can be expressed as $F(t)=\delta^{3}(t)F_{0}$, where $F_{0}$ is the flux in the rest frame, and $\delta(t)$ is the Doppler factor which can be written as $\delta(t)=\frac{1}{\gamma(1-\beta\cos\theta_{obs})}$, where $\beta$ is the bulk velocity and $\gamma=\left(1-\beta^{2}\right)^{-\frac{1}{2}}$ is the Lorentz factor. In addition, before fitting, we transform the flux into luminosity through $L=4\pi D^{2}F$, where $D$ is the luminosity distance. In our model, the observed luminosity can be obtained as

where $L_{1}$ and $L_{2}$ are the luminosity of jet 1 and jet 2, the subscript 1, 2 of $\delta$ refers to the Doppler factor of the two jets, and $L_{b}$ is the contribution from other sources such as the disc, corona, etc. For simplification, the precession of the disc, gravitational redshift, and intrinsic variety of accretion flow is not considered here, therefore, for convenience, we assume that $L_{b}$ is not a time-dependent physical quantity in our model.

In the fitting, we only fixed the orbital angular velocities ($\omega=\frac{2\pi}{P}$) of both black holes. To be consistent with the case of $\gamma$-ray, the period was frozen as $P=2.2~{}\text{yr}$. Additionally, $L_{b}$ was frozen as $8.0\times 10^{45}\text{erg}~{}\text{s}^{-1}$ which is the minimum value in the light curve. Meanwhile, with assumption $\phi_{1}=\phi_{2}$, other parameters in the model were freed. However, Villata et al. (1998) assumed the two jets share the same values of the parameters when they tried to explain the quasi-periodic double peak structure in the optical light curves of OJ 287. In our work, we do not set the parameters of the two jets (except for the orbital angular velocities and $\phi_{1}=\phi_{2}$) to be equal. The Markov Chain Monte Carlo (MCMC) emcee package (Foreman-Mackey et al., 2013) was used in the fitting and 500 walkers and 10000 steps were run⁶⁶6The chain was accessed by using the EnsembleSampler.get_chain in python package emcee. Within this method, we obtained the parameter values as a function of the step numbers in the chain. When the parameters do not change much, the chain can be considered stable. . We adopt the uniform prior as $0.0<\phi_{1}<10^{\circ}$, $0.0<\alpha_{1},\alpha_{2}<10^{\circ}$, $0.0<\psi<2\pi$, $0.0<\gamma_{1},\gamma_{2}<50.0$, $0.0<L_{1},L_{2}<5.0\times 10^{44}\text{erg}~{}\text{s}^{-1}$. The fitted result is shown in Figure 2. More details for the fitted parameters in our model are shown in Table 1.

The fitting is unsatisfactory, with a reduced chi-square of 348.41. However, we did not expect a perfect result based on the X-ray light curve, mainly because of the simple model and the lack of continuous observations for this object. Nevertheless, we had tried to free the period parameter in the fitting, but in the results, the fitted period of 2.8 years is inconsistent with the case of $\gamma$-ray (Ackermann et al., 2015). Tavani et al. (2018) showed the consistency between the R-band and X-ray light curves in the long-term observations. Additionally, Sandrinelli et al. (2018) and Tavani et al. (2018) reported a similar behavior in the long-term $\gamma$-ray and the R-band light curve. Therefore, we can assume that the $\gamma$-ray and X-ray light curves are correlated to some extent. Hence, in fitting, it is reasonable to fix the periodicity to 2.2 years from continuously monitoring $\gamma$-ray data (Caproni et al., 2017; Sandrinelli et al., 2018; Cavaliere et al., 2019).

As mentioned in the last subsection, we only consider the geometric effect in the system in this work. The profile of the X-ray light curve of PG 1553+113 is fairly complicated, while the fitted model in this work is rather simple and it can not describe all the data points in the X-ray. We have tried to improve the model to fit the X-ray light curve and some geometric effects were considered in the test. Although the statistical fitting result is slightly improved, the data sampling still plays an important role. The sparse sample can be fitted well by such complex models may not indicate the true nature of the system. To further investigate the detail of the variability, more high cadence observations are needed. The origin of the X-ray is very complex and there exist some other mechanisms that would modulate the variability. Based on the assumption that the quasi-periodic variation in the X-ray light curve is origin from the emission of the jets, we only explore the effect of the two jets’ precession. However, according to Sobacchi et al. (2017), the Lense-Thirring precession would modulate the oscillation of angle with the time scale from days to years and this may affect on the X-ray variability. Additionally, we did not consider the intrinsic variation of the jets, so our model can not mirror the activity in the jets. For these reasons, the predicted light curve may not fit the observed data very well.

Although the data points in the light curve can not be described well by the two-jet precession scenario, the main and weak flares can be explained well. In order to predict variability accurately, more aspects need to be taken into account. However, in this work, we mainly focus on the effect caused by the precession of two jets in the SMBHB system.

The mass of each black hole was not directly estimated in the previous investigations. We can only constrain the orders of magnitude of the two black holes masses through the fitting results. Following the relation between the power of jets and the luminosity of the discs, the mass of each black hole can be estimated. According to Ghisellini et al. (2014),

where $P_{rad}$ is the radiative power of the jet, and $L_{disc}$ is the disc luminosity. Assuming that X-ray originated from synchrotron emission, the radiative power of the jet can be written as

where $L_{bol,jet}$ is the bolometric luminosity of the jet which can be estimated following the method of Runnoe et al. (2012).

Substituting the fitting results from Table 1 into Equations (5) and (4), and assuming the accretion rate as 0.1 $\dot{M}_{\text{edd}}$ (following Ackermann et al., 2015), where $\dot{M}_{\text{edd}}$ is the Eddington accretion rate, we obtain the masses of the black holes launching jet 1 and jet 2 are $M_{1}\simeq 1.40\times 10^{8}M_{\sun}$ and $M_{2}\simeq 3.47\times 10^{8}M_{\sun}$, respectively. The inferred mass ratio is $q=\frac{M_{1}}{M_{2}}\simeq 0.41$ that is in agreement with Tavani et al. (2018). The total mass of the SMBHB system is consistent with the results of Ackermann et al. (2015) and Cavaliere et al. (2019).

The X-ray energy spectra were analyzed with the redshift power-law model (zpowerlw) through xspec version 12.10.1f. The mean value of the power-law index is $2.38\pm 0.02$ with the maximum and minimum value of $2.98^{+0.47}_{-0.42}$ and $1.81^{+0.27}_{-0.24}$ respectively, and the results can be seen on the middle panel of Figure 3. Centering at the peak time, we define the 80% of the interval between the two predicted time of minimum luminosity as the flare epoch. We can also define three states to describe the X-ray activity: quiescent, weak flares, and main flares, which corresponding to the different epochs in the results of the two-jet precession model.

To study the evolution of the X-ray flux of PG 1553+113, we analyzed the hardness ratio (HR) expressed as

where CR(2.0–10.0 keV) and CR(0.3–2.0 keV) are counts rate of the hard and soft X-ray, respectively. The evolution of HR is shown on the bottom panel of Figure 3. From 2012.0 to 2020.0, the average value of HR is $-0.689\pm 0.004$ with the maximum and minimum value are $-0.468\pm 0.062$ and $-0.838\pm 0.058$ respectively, which exhibits a very soft level in the stage.

We analyse X-ray data acquired by the XMM-Newton Observatory. The main flare spectrum was obtained in September 6, 2001 (ID 0094380801), the weak flare was observed in July 24, 2013 (ID 0727780101), and the spectrum for a quiescent state with observation ID 0810830201 was obtained in August 27, 2019. The spectra were analysed with xspec, and final results are presented in Figure 4. The X-ray energy spectra exhibit different behavior at these epochs. Raiteri et al. (2015) compared the spectral energy distribution (SEDs) for three states, and from their results, the quiescent state SEDs showed the significant different behavior from the flare states in $\gamma$-ray, which may be the contribution of the jet’s activity.

The power-law index and hardness ratio stay at a relatively stable level most of the time. The best linear fit of the correlation between flux and the power-law index are shown in Figure 5, and the results reveal a slightly “harder-when-brighter” feature during the main and weak flares while for the quiescent state, the slightly “softer-when-brighter” tendency can be seen. We explore the correlation between the index and flux with Spearman test. The analysis is processed by the python package scipy.stats.spearmanr. During main flares, the correlation coefficients ($r_{s}$) $<-0.2$ for all cases with p-value $<0.05$. For the case of weak flares, $r_{s}<-0.1$ and the p-value $<0.7$ (the result from XMM-Newton shows the p-value of 0.675). The result of quiescent state reveals $r_{s}<0.2$ with p-value $<0.6$. It should be noted that in the quiescent state, we obtain only one data point from XMM-Newton, hence the correlation test for this instrument can not be applied. More details for the linear fitting and Spearman test are shown on Table 2.

When we fit all the data together, the slop is very different from the separated fitting results. We note that this behavior is related to the distribution of data points. 8 XMM data points for weak flares locate on the low flux and the power-law index $>2.5$. And for the case of XRT, few points locate on such interval. Therefore, the XMM data points decrease the slope when we fit the data from two instruments together.

The similar evolution of the spectral index for both main and weak flares suggest the same radiation mechanism, which may be related to the jets. The different behavior of the quiescent states and flares may be the contribution of the precession of jets. When the jet pointing to the observer, the “harder-when-brighter” tendency in the flares states can be seen, and in the quiescent state, with the fraction of the contribution from the disc increase, the spectra tend to be “softer-when-brighter” (Wang & Jiang, 2020).

In the previous subsections, we argue that the X-ray variation is dominated by the precession of jets, however, the other potential contributions may be considered. We discuss the lighthouse effect, general relativity effect, disc instability, and Lense-Thirring precession below. But the timescale of the variation from these effects is inconsistent with the 2.2-year quasi-periodicity, however, we can not exclude out the modulation of these effects on the long and short timescale variability of PG 1553+113.