A sub-pc BBH system in SDSS J1609+1756 through optical QPOs in ZTF light curves

SDSS J1609+1756  is collected as target of the manuscript, due to two main reasons. First, as a candidate of off-nucleus AGN in Ward et al. (2021), SDSS J1609+1756  has peculiar broad Balmer lines with shifted red shoulders which can be well explained by broad line emissions from central two independent BLRs related to a central BBH system, quite similar as what we have recently discussed in Zhang (2021d). Second, after checking long-term variabilities of SDSS J1609+1756, apparent optical QPOs can be detected to support an expected BBH system.

The 4.45 years-long ZTF g/r/i-band light curves of SDSS J1609+1756  are collected and shown in top left panel of Fig 1, with MJD-58000 from 203 (Mar. 2018) to 1828 (Jul. 2022). Meanwhile, long-term light curves of SDSS J1609+1756  are also checked in CSS, PanSTARRS and ASAS-SN, and shown in Fig. 2 without apparent variabilities. There are 390 reliable data points included in the CSS light curve, however, as shown in left panel of Fig. 2, more than 92% of the 390 data points are lying between the range of mean magnitude plus/minus corresponding 1RMS scatters, indicating not apparent variabilities in the CSS light curve, probably due to lower quality of light curves. There are 1 reliable data point and 20 reliable data points included in the ASAS-SN V/g-band light curves shown in middle panel of Fig. 2. However, due to he quite large time steps with mean value about 200 days, it is not appropriate to search QPOs in ASAS-SN light curves. Similar as ASAS-SN g-band light curve, there are only 10 data points, 14 data points, 21 data points, 11 data points and 13 data points included in the PanSTARRS g/r/i/z/y-band light curves shown in right panel of Fig. 2, respectively, with large mean time steps about 266 days, 217 days, 184 days, 20 days and 182 days. Therefore, the PanSTARRS light curves are not considered in the manuscript.

Finally, the long-term ZTF light curves are mainly considered and there are no further discussions on variabilities from the other Sky Survey projects.

Based on the high-quality long-term optical light curves from ZTF, the following four methods are applied to detect optical QPOs in SDSS J1609+1756, similar as what have been done in Penil1 et al. (2020).

Due to short time durations of ZTF light curves, it is necessary and interesting to check whether the determined optical QPOs was mis-detected QPOs tightly related to central intrinsic AGN activities of SDSS J1609+1756, although four different mathematical methods applied in the section above can provide robust evidence to support the detected optical QPOs in SDSS J1609+1756. Similar as what we have recently done in Zhang (2022a, c) to check probability of mis-detected QPOs in AGN activities in SDSS J0752 and in SDSS J1321, the following procedure is applied.

As well discussed in Kelly, Bechtold & Siemiginowska (2009); Kozlowski et al. (2010); MacLeod et al. (2010); Zu et al. (2013, 2016); Zhang & Feng (2017); Takata, Mukuta & Mizumoto (2018); Moreno et al. (2019); Sheng, Ross & Nicholl (2022), the known Continuous AutoRegressive (CAR) process and/or the improved Damped Random Walk (DRW) process can be applied to describe the fundamental AGN activities (Rees, 1984; Ulrich et al., 1997; Madejski & Sikora, 2016; Baldassare et al., 2020; Burke et al., 2021). Here, based on the DRW process, the public code JAVELIN (Just Another Vehicle for Estimating Lags In Nuclei) (Kozlowski et al., 2010; Zu et al., 2013) is firstly applied to describe the ZTF r-band light curve, with two process parameters of intrinsic characteristic variability amplitude and timescale of $\sigma$ and $\tau$. The best descriptions are shown in left panel of Fig. 8. And corresponding MCMC determined two dimensional posterior distributions of $\sigma$ and $\tau$ are shown in right panel of Fig. 8, with the determined $\ln(\tau/days)\sim 4.95\pm 0.35$ ($\tau\sim 140_{-40}^{+60}$ days) and $\ln(\sigma/(mag/dyas^{1/2}))\sim-1.45\pm 0.25$. Here, because of the totally same periodicities through the ZTF r-band light curve by different methods in the section above, the r-band rather than the g-band light curve is considered in the section. Certainly, Fig. 8 also shows corresponding results through the g-band light curve.

Then, probability of mis-detected QPOs from DRW process described intrinsic AGN variabilities can be estimated as follows, through applications of the CAR process discussed in Kelly, Bechtold & Siemiginowska (2009):

with $\epsilon(t)$ as a white noise process with zero mean and variance equal to 1. Here, the value 19.21, the mean value of the ZTF r-band light curve of SDSS J1609+1756, is set to be the mean value of $LC_{t}$. Then, a series of 100000 simulating light curves [$t$, $LC_{t}$] are created, with $\tau$ randomly selected from 100 days to 200 days (range of the determined $\tau$ of SDSS J1609+1756, after considering uncertainties). Then, variance $\tau\sigma_{*}^{2}/2$ of CAR process created light curve is set to be 0.168 which is the variance of ZTF r-band light curve of SDSS J1609+1756. And, time information $t$ are the same as the observational time information of ZTF r-band light curve of SDSS J1609+1756. And similar uncertainties $LC_{t}\times\frac{LC_{err}}{LC_{r}}$ are simply added to the simulating light curves $LC_{t}$, with $LC_{r}$ and $LC_{err}$ as the ZTF r-band light curve and corresponding uncertainties of SDSS J1609+1756.

Among the 100000 simulated light curves, there are 263 light curves collected with reliable mathematical determined periodicities according to the following four simple criteria. First, the simulated light curve can be well described by equation (1) with corresponding $\chi^{2}/dof$ smaller than 8 ($\chi^{2}/dof\sim 5.8$ in Fig. 1). Second, the simulated light curve has apparent GLS periodogram determined peak with corresponding periodicity smaller than 500 days with confidence level higher than $5\sigma$. Third, the simulated light curve has apparent ACF and WWZ method determined peak with corresponding periodicity smaller than 500 days with confidence level higher than $5\sigma$. Fourth, the determined periodicity by equation (1) is consistent with the GLS periodogram determined periodicity (which are similar as the ACF and WWZ method determined periodicities) within 10 times of the determined uncertainties by applications of equation (1). Left panel of Fig. 9 shows properties of the determined $\chi^{2}/dof$ related to the best fitting results determined by equation (1) and the GLS periodogram determined periodicities (which are similar as the ACF and WWZ method determined periodicities) of the 263 simulated light curves. And right panel of Fig. 9 shows one of the 263 light curves with best fitting results by equation (1).

Therefore, even without any further considerations, it can be confirmed that the probability is only 0.26% (263/100000) (confidence level higher than $3\sigma$) to support mis-detected QPOs in CAR process simulated light curves related to AGN activities. Furthermore, accepted uncertainty $\Delta_{T}=5$ days of optical periodicity in ZTF r-band light curve, if to limit periodicities within range from $344\pm 5\times\Delta_{T}$ in the simulated light curves, there are only 112 simulated light curves collected, leading to the probability only about 0.11% (112/100000) (confidence level higher than $3.2\sigma$) to support mis-detected QPOs in AGN activities. In other words, the confidence level higher than $3\sigma$ to confirm the optical QPOs not from intrinsic AGN activities in SDSS J1609+1756, although short time durations in ZTF light curves, after well considering effects of AGN activities described by CAR process.

Fig. 10 shows the SDSS spectrum with PLATE-MJD-FIBERID=2200-53875-0526 of SDSS J1609+1756  collected from SDSS DR16 (Ahumada et al., 2021). The apparent red-shifted shoulders, marked in Fig. 10, can be found in broad Balmer emission lines, as shown in Ward et al. (2021). And the shoulders should be more apparently determined after measurements of emission lines in SDSS J1609+1756  by the following emission line fitting procedure.

Multiple Gaussian functions can be applied to well simultaneously measure the emission lines around H$\beta$ within rest wavelength from 4600Å  to 5100Å  and around H$\alpha$ within rest wavelength from 6150Å  to 6850Å  in SDSS J1609+1756, similar as what we have recently done in Zhang (2021a, b, c); Zhang & Zhao (2022b). Three broad and one narrow Gaussian functions are applied to describe broad and narrow components in H$\alpha$ (H$\beta$). And two narrow and two broad Gaussian functions are applied to describe [O iii]$\lambda 4959,5007$Å  doublet, for the core components and for the components related to shifted wings. And one Gaussian function is applied to describe He ii$\lambda 4687$Å  emission line. And six Gaussian functions are applied to describe the [O i]$\lambda 6300,6363$Å, [N ii]$\lambda 6549,6583$Å  and [S ii]$\lambda 6716,6731$Å  doublets. As the following shown best fitting results, it is not necessary to consider components related shifted wings in the [O i], [N ii] and [S ii] doublets. And, a power law function is applied to describe the AGN continuum emissions underneath the emission lines around H$\beta$ (H$\alpha$).

When the model functions are applied, only three restrictions are accepted. First, emission flux of each Gaussian emission component is not smaller than zero. Second, Gaussian components of each doublet have the same redshift. Third, corresponding Gaussian components in broad Balmer line have the same redshift. Then, through the Levenberg-Marquardt least-squares minimization technique (the known MPFIT package), best fitting results to emission lines can be determined and shown in Fig. 11 with $\chi^{2}/dof\sim 1.09$. The measured parameters of each Gaussian component are listed in Table 3. Here, although three Gaussian functions are applied to describe each broad Balmer line, only two Gaussian components have reliable measured parameters which are at least 3times larger than their corresponding uncertainties, indicating two broad Gaussian components are enough to describe each broad Balmer line. Moreover, it is clear that there are red shifted shoulders in each broad Balmer emission line, especially based on the red-shifted broad Gaussian component in each broad Balmer line.

The main objective to discuss emission line properties of broad Balmer lines is that the shoulders can provide information on peak separation of two broad components coming from two BLRs related to an expected central BBH system, which will provide further clues on maximum space separation of the two black holes in the expected BBH system in SDSS J1609+1756. The red-shifted velocity of the shoulder in broad Balmer lines are about $V_{r}=2200\pm 300{\rm km/s}$, based on the red-shift broad Gaussian component in each broad Balmer line. Meanwhile, peak separation $V_{p}$ in units of km/s related to a BBH system with space separation $S$ and with BH masses of $M_{BH1}$ and $M_{BH2}$ of the two BHs can be simply described as

with $i$ as inclination angle of the orbital plane to line-of-sight and $\phi$ as orientation angle of orbital phase. Therefore, considering the observed $V_{r}~{}\sim~{}2200{\rm km/s}~{}<~{}V_{p}=V_{r}+V_{b}$ ($V_{b}$ as blue-shifted velocity of undetected blue-shifted shoulders in broad Balmer lines in SDSS J1609+1756), the maximum space separation $S$ can be simply determined as

. Once there are clear information of central total BH mass which will be discussed in the next section, the maximum space separation can be well estimated.

For a broad line AGN, the virialization assumption accepted to broad line emission clouds is the most convenient method to estimate central virial BH mass as discussed in Peterson et al. (2004); Greene & Ho (2005); Vestergaard & Peterson (2006); Shen et al. (2011). However, considering mixed broad Balmer line emissions, it is difficult to determine the two clear components from central two independent BLRs related to central BBH system in SDSS J1609+1756, after checking the following point.

If the determined two broad Gaussian components in each broad Balmer line shown in Fig. 11 were truly related to a central BBH system, the virial BH mass as discussed in Greene & Ho (2005) of each BH in central region can be estimated as

with $f_{\alpha,r}$ and $\sigma_{\alpha,r}$ ($f_{\alpha,b}$ and $\sigma_{\alpha,b}$) as line flux and line width (second moment) of red-shifted (blue-shifted) broad component in broad H$\alpha$, after considering the more recent empirical R-L relation to estimate BLRs sizes through line luminosity in Bentz et al. (2013). Then, under assumption of a BBH system, shift velocity ratio $R_{v}$ of red-shifted broad component to blue-shifted broad component in broad H$\alpha$ can be estimated as

accepted the measured line widths and line fluxes and uncertainties listed in Table 3.

However, according to measured central wavelengths of the two broad components in broad H$\alpha$ listed in Table 3, the observed shift velocity ratio $R_{v,obs}$ is

with 6564.61 in units of Å  as the theoretical value of central wavelength of broad H$\alpha$ in rest frame, which is quite different from the $R_{v}$. Therefore, although there are two broad Gaussian components determined in broad H$\alpha$, the two components are not appropriate to be applied to determine virial BH masses of central two black holes in the expected BBH system. Similar results can also be found through the two components in broad H$\beta$.

Therefore, rather than virial BH mass estimated through broad line luminosity and broad line width, BH mass estimated by continuum luminosity as shown in Peterson et al. (2004) is determined as

with $L_{1}$ and $L_{2}$ in units of ${\rm erg/s}$ as continuum luminosity at 5100Å  coming from central two BH accreting systems, then upper limit of central total BH mass $M_{BH}=M_{BH1}+M_{BH2}$ can be estimated as

. Then, based on the fitting results shown in right panels of Fig. 11, total continuum luminosity at 5100Å  in rest frame is $L_{t}~{}\sim~{}(1.47\pm 0.02)\times 10^{44}{\rm erg/s}$. Simply accepted $L_{t}~{}=~{}L_{1}~{}+~{}L_{2}$, upper limit of central total BH mass is about $(1.03\pm 0.22)\times 10^{8}{\rm M_{\odot}}$.

Accepted the estimated upper limit of total BH mass, the upper limit of space separation of the expected central BBH system in SDSS J1609+1756  is about

. Meanwhile, considering optical periodicity $\sim$340 days, as discussed in Eracleous et al. (2012), the space separation of the central BBH system can be estimated as

with accepted total BH mass $(1.03\pm 0.22)\times 10^{8}{\rm M_{\odot}}$. The space separation $S_{BBH}$ determined by periodicity is well below the upper limit of space separation determined by peak separation $V_{p}$, not leading to clues against assumptions of the central BBH system in SDSS J1609+1756.

Before end of the subsection, one point should be noted. As simply discussed above, through optical periodicity $\sim$340 days, estimated space separation of the central BBH system is about 2.6 light-days in SDSS J1609+1756. Meanwhile, continuum luminosity about $10^{44}{\rm erg/s}$ in SDSS J1609+1756  can lead BLRs sizes to be about 36 light-days through the R-L relation in Bentz et al. (2013), quite larger than $S_{BBH}\sim 2.6$light-days. Therefore, there should be few effects of dynamics of central BBH system on emission clouds of probably mixed BLRs in SDSS J1609+1756, or only apparent effects on emission clouds in inner regions of central mixed BLRs of SDSS J1609+1756. The results can provide further clues to support that it is not appropriate to estimate central BH mass by properties of broad emission lines in SDSS J1609+1756  as well discussed above, and also to support that the shifted velocity ratio discussed above should be not totally confirmed to be related to dynamics of central BBH system. Multi-epoch monitoring of variabilities of broad emission lines should provide further and accurate properties of dynamic structures of central expected BBH system in SDSS J1609+1756  in the near future.

Meanwhile, besides the expected central BBH system, precessions of emission regions with probable hot spots for the optical continuum emissions can also be applied to describe the detected optical QPOs in SDSS J1609+1756. As discussed in Eracleous et al. (1995) and in Storchi-Bergmann et al. (2003), the expected disk precession period can be estimated as

, with $R_{3}$ as distance of optical emission regions to central BH in units of 1000 Schwarzschild radii ($R_{g}$) and $M_{8}$ as the BH mass in units of $10^{8}{\rm M_{\odot}}$. Considering optical periodicity about $T_{\rm pre}\sim 340$ days and BH mass about $(1.03\pm 0.22)\times 10^{8}{\rm M_{\odot}}$ above estimated through the continuum luminosity, the expected $R_{3}$ could be around 0.06 in SDSS J1609+1756.

However, based on the discussed distance of NUV emission regions to central BHs in Morgan et al. (2010) through the microlensing variability properties of eleven gravitationally lensed quasars, the NUV 2500Å  continuum emission regions in SDSS J1609+1756  have distance from central BH as

leading size of NUV emission regions to be about $60R_{g}$. The estimated NUV emission regions have similar distances as the optical continuum emission regions in SDSS J1609+1756  under the disk precession assumption, strongly indicating that the disk precessions of emission regions are not preferred to be applied to explain the detected optical QPOs in SDSS J1609+1756.

Moreover, long-term QPOs can be detected in blazars due to jet precessions as discussed in Sandrinelli et al. (2018); Bhatta (2019); Otero-Santos et al. (2020). However, SDSS J1609+1756  is covered in Faint Images of the Radio Sky at Twenty-cm (Becker, White & Helfand, 1995; Helfand et al., 2015), but no apparent radio emissions. Therefore, jet precessions can be well ruled out to explain the optical QPOs in SDSS J1609+1756.

Furthermore, it is interesting to discuss whether the known relativistic Lense-Thirring precession (Cui, Zhang & Chen, 1998; Wagoner, 2012) can be applied to explain the detected optical QPOs in SDSS J1609+1756. As well discussed in Cui, Zhang & Chen (1998), observed periodicity related to the Lense-Thirring precession is about

with $z$ as redshift, $M_{BH}$ as BH mass in units of $\rm M_{\odot}$, $R_{e}$ in units of $R_{g}$ as distance of emission regions in central accretion disk to central BH and $a_{*}$ (between $\pm$1) as dimensionless BH spin parameter. If accepted central BH mass as $(1.03\pm 0.22)\times 10^{8}{\rm M_{\odot}}$, minimum value of $R_{e}$ as $60R_{g}$ (the estimated value for NUV emission regions above) for optical emission regions and maximum value $|a_{*}|=1$, the minimum $P_{LT,obs}$ can be estimated to be about 5200days, about 13 times larger than the detected optical periodicity about 340days in SDSS J1609+1756. Therefore, the detected optical QPOs in SDSS J1609+1756  are not related to relativistic Lense-Thirring precessions.

Before ending of the manuscript, one point is noted. The SDSS J1609+1756  is collected from the candidates of off-nucleus AGN reported in the literature. However, it is not unclear whether are there tight connections between optical QPOs and off-nucleus AGN, studying on a sample of off-nucleus AGN with apparent optical QPOs could provide further clues on probable intrinsic connections in the near future.