The possible long-term periodic variability of the extremely luminous quasar WISE J090924.01+000211.1

Abstract

The extremely luminous infrared galaxy (ELIRG), WISE J090924.01|$+$|000211.1 (hereafter WISE J0909|$+$|0002, |$z=1.87$|⁠) is an extraordinary object with a quasar aspect. This study performed monitoring observations of WISE J0909|$+$|0002 with the 105|$\:$|cm Murikabushi telescope, Okayama and Akeno 50|$\:$|cm telescopes/MITSuME (⁠|$g^{\,\prime}$|⁠, |$R_{\rm c}$|⁠, and |$I_{\rm c}$| bands), and the SaCRA 55|$\:$|cm telescope/MuSaSHI (r, i, and z bands) for three years. We obtain the following results by combining the UV/optical light curves of the CRTS, Pan-STARRS, and ZTF archive data, and our observational data: (1) the light curves of WISE J0909|$+$|0002 present quasi-periodic (sinusoidal) oscillations with the rest-frame period of |$\sim \!660$|–|$689\:$|d; (2) the structure functions of WISE J0909|$+$|0002 do not show a damped random walk (DRW) trend; (3) the mock DRW light curves present periodic-like trend on rare occasions in 10000 simulations; (4) the relativistic boost scenario is favored, since the relation between variability amplitude and power-law slope ratio is consistent with the theoretical prediction of this scenario, and a substantial parameter space exists between the inclination angles and the black hole mass; (5) the circumbinary disk model is insufficient to explain the spectral energy distribution of our target; (6) the significant radio flux density of WISE J0909+0002 is not detected from the VLA FIRST Survey, thus the radio jet precession scenario is ruled out. The Doppler boost scenario is, from our results, likely as a cause of the periodic variability; consequently the quasi-periodic oscillations in WISE J0909|$+$|0002 is possibly interpreted by a supermassive black hole binary. Additional observations to investigate the continuity of the periodic trend would bring new insights into mechanisms of the quasi-periodic oscillations and/or ELIRGs.

1 Introduction

1.1 Periodic variability of active galactic nuclei

The flux variability of active galactic nuclei (AGNs) is a universal phenomenon from X-ray to radio wavelengths. Attempts have been made to describe the AGN variability as the damped random walk (DRW) model (e.g., Kelly et al. 2009; Kozłowski et al. 2010; MacLeod et al. 2010). However, a periodic flux variability with a light curve that was well fitted by a sinusoidal function with a period of |$1884\:$|d was discovered in the quasar PG 1302|$-$|102 (Graham et al. 2015b).

Such quasi-periodic oscillations (QPOs) have been actively studied, for instance, the Pan-STARRS1 Medium Deep Survey data with |$griz$| bands identified 26 candidates of periodic light curves (Liu et al. 2019). A recent study identified 106 (DRW as a null hypothesis) and 86 (single power-law as a null hypothesis) candidates of periodic quasars from the photometric data of the Zwicky Transient Facility (Chen et al. 2024). Overall, the detection rate of quasars with QPOs is |$\sim \!0.01\%$|–|$1.1\%$| (Charisi et al. 2016; Liu et al. 2019; Chen et al. 2020) As examples of individual objects, periodicities of between half a year and several years were discovered in the following AGNs: SDSS J025214.67|$-$|002813.7 (⁠|$z=1.53$|⁠; Liao et al. 2021), Q J0158-4325 (⁠|$z=1.29$|⁠; Millon et al. 2022), PKS 2131|$-$|021 (⁠|$z=1.285$|⁠; O’Neill et al. 2022), SDSS J132144|$+$|033055 (⁠|$z=0.269$|⁠; Zhang 2022), and SDSS J143016.05|$+$|230344.4 (⁠|$z=0.08$|⁠; Dou et al. 2022; Hoshi et al. 2024). The light-curve modeling for SDSS J025214.67|$-$|002813.7 was improved by considering periodic components (Covino et al. 2022).

While the mechanism of the periodic light curves is not well known, a binary supermassive black hole (BSBH) within a milli- to sub-parsec scale is the most likely candidate of that phenomenon (e.g., Komossa 2006; Graham et al. 2015b; D’Orazio et al. 2015; Duffell et al. 2020). That is to say, AGNs with periodic light curves are possibly sites emitting gravitational waves (e.g., Agazie et al. 2023). The BSBH system is thought to yield an orbital motion of primary and secondary mini-accretion disks with a relativistic Doppler boost (D’Orazio et al. 2015). The Doppler boost scenario is partially supported by previous studies (e.g., Charisi et al. 2018; Xin et al. 2020). However, the All-sky Automated Survey for Supernovae (ASAS-SN) light curve of PG 1302|$-$|102 presented a perturbation in its periodicity, which was identified as a false-positive signal and does not support the BHSH theory (Liu et al. 2018; Kovačević et al. 2019).

Other plausible scenarios of the periodic light curves are, for instance, accretion-disc precession (e.g., Eracleous et al. 1995; Storchi-Bergmann et al. 2003) and microlensing by binary stars. The former scenario was denied by the analysis of Zhang (2022). The latter was also ruled out in the lensed quasar Q J0158-4325, since this scenario required a |$\sim \!2000\:$|yr period from dynamical simulation of a binary star with microlensing (Millon et al. 2022). At present there is no clear consensus on the mechanism of QPOs in AGNs. In order to overcome this situation, we focused on extremely luminous infrared galaxies, as they are believed to be in the final stage of galaxy (and may be SMBH) merging (ELIRGs; Tsai et al. 2015; Toba et al. 2018, 2020).

1.2 ELIRGs and our target

Galaxy mergers are an important phase in the evolution of galaxies. It is known that ultra-luminous infrared (IR) galaxies (ULIRGs) and hyper-luminous IR galaxies (HyLIRGs; Rowan-Robinson 2020) indicate the evidence of the final stage of galaxy interactions. Within the last decade, a series of ELIRGs with the infrared luminosity of |$\gt \!10^{14}\, L_\odot$| have been discovered, and these objects also exhibit a sign of galaxy interaction with thick dust. Since such infrared-luminous galaxies are considered to be the peak stage of the co-evolution of galaxies (e.g., Hopkins et al. 2008; Blecha et al. 2018; Yutani et al. 2022), ELIRGs are important research targets. Toba et al. (2021) newly identified WISE J090924.01|$+$|000211.1 (RA = |${09^{\rm h}09^{\rm m}24{^{\rm s}}}$|⁠, Dec = |$00^{\circ }$||${02^{\prime }11^{\prime \prime }}$|⁠, |$z=1.87$|⁠; hereafter WISE J0909|$+$|0002) as an ELIRG and reported its interesting properties: (1) the infrared luminosity is |$L_{\rm IR} = 1.79 \times 10^{14}\, L_\odot$|⁠, (2) the SDSS spectrum presents broad emission lines, and the estimated black hole mass is extremely massive with a mass of |$7.4\times 10^9\, M_{\odot }$|⁠, and (3) the extremely high infrared luminosity with quite low column density may imply that this object is in the midst of a short-lived phase. As shown in figure 1, WISE J0909|$+$|0002 presents an exceedingly high bolometric luminosity compared to most quasars. Namely, the ELIRG WISE J0909|$+$|0002 includes an extremely luminous type 1 quasar. However, the internal structure of WISE J0909|$+$|0002 is still unclear due to its outlandish properties.

Fig. 1.

Redshift dependence of quasar bolometric luminosity from the Sloan Digital Sky Survey (SDSS) Data Release 7 (Shen et al. 2011). The contours indicate that the distribution is concentrated at |$z \sim 1.5$|⁠, |$L_{\rm bol}\sim 46.5$|⁠. The star, circle, and triangle represent the data points for WISE J0909|$+$|0002 from Toba et al. (2021), the periodic quasars J024703.24|$-$|010032.0 from Chen et al. (2020), and PG 1302|$-$|102, respectively.

We independently found the periodic variability trend of WISE J0909|$+$|0002 from archival data and our monitoring observations. This study provides QPOs of this object and the related analysis of its light curves. We describe the details of photometric data of WISE J0909|$+$|0002 from archives and observation with the observation collaboration of the Optical and Infrared Synergetic Telescopes for Education and Research (OISTER) in section 2. Section 3 presents the flux variability trend WISE J0909|$+$|0002, and we then discuss the inferred physical properties of this object in section 4. The conclusion of our study is summarized in section 5. A |${\rm \Lambda }$|CDM cosmology with |$H_{0}= 70\:$|km|$\:$|s|$^{-1}\:$|Mpc|$^{-1}$|⁠, |$\Omega _{m} = 0.3$|⁠, and |$\Omega _{\Lambda } = 0.7$| is adopted in this study.

2 Archive data and our observation

2.1 Archive data

In order to measure the flux variability of WISE J0909|$+$|0002, we use the archival photometric data as follows: the Catalina Real-Time Transient Survey (CRTS; Drake et al. 2009; Mahabal et al. 2011), the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS; Chambers et al. 2016), and the Zwicky Transient Facility (ZTF). In addition, we conducted a follow-up observation of WISE J0909|$+$|0002 with the OISTER collaboration. Table 1 summarizes details of those archival and observational data. Details of those photometric data are explained in the following sub-subsections.

Table 1.

List of archival and our observational data of WISE J0909|$+$|0002.

			m†	|$\sigma _{m}$|‡	Observation epoch
	Filter	N*	(mag)	(⁠|$10^{-2}$| mag)	(MJD)
Archival data
CRTS	V	330	16.46	7.0	53494.1–56397.1
Pan-STARRS	g	16	16.83	0.5	55239.4–56715.4
	r	7	16.72	0.5	56315.4–56662.6
	i	21	16.32	0.4	55577.4–57021.5
	z	15	16.39	0.7	55284.3–56991.6
	y	15	16.37	1.0	55523.5–56991.6
ZTF	g	5	16.72	1.0	58430.5–59271.3
	r	42	16.59	1.0	58439.3–59638.2
Our observations
Murikabushi/MITSuME	|$g^{\prime }$|	35	16.75	1.0	59257.6–59619.6
	|$R_{\rm c}$|	35	16.68	0.9	59257.6–59619.6
	|$I_{\rm c}$|	34	16.38	1.0	59257.6–59619.6
Akeno and Okayama/MITSuME	|$g^{\prime }$|	33	16.70	3.0	60079.5–60401.5
	|$R_{\rm c}$|	42	16.60	3.0	60079.5–60401.5
	|$I_{\rm c}$|	33	16.45	4.0	60079.5–60401.5
SaCRA/MuSaSHI	r	4	16.60	3.0	60269.7–60419.5
	i	4	16.42	2.0	60269.7–60419.5
	z	4	16.36	2.0	60269.7–60419.5

*

Number of data points.

†

Median magnitudes over the observation period.

‡

Median photometric errors over the observation period.

2.1.1 CRTS data

The CRTS operates three telescopes: the Catalina Sky Survey (CSS) 0.7|$\:$|m Schmidt, the Mount Lemmon Survey (MLS) 1.5|$\:$|m Cass in Arizona, and the Siding Springs Survey (SSS) 0.5|$\:$|m Schmidt in Australia. Those telescopes sweep sky fields of up to |$\sim \!2500\:$|deg|$^2$| with four observations per visit on a clear night. The archival data taken by CRTS are converted to the V-band magnitude based on 50–100 G-type stars using 2MASS data (Drake et al. 2013). The CSS light curve of WISE J0909|$+$|0002 covers the period from MJD 53494 to 56397 with 330 data points. During this period, the light curve showed a QPO trend (see subsection 3.1).

2.1.2 Pan-STARRS data

Pan-STARRS carries out three-day cadence |$grizy$| band observations with the 1.8|$\:$|m telescope. The |$grizy$| bandwidths are listed in Tonry et al. (2012). We converted “psfFlux” (in units of Jy) into an AB magnitude, of which the zero-point flux is |$3631\:$|Jy. The epoch of the photometric data for WISE J0909|$+$|0002 ranges from MJD 55239.4 to 57021.5.

2.1.3 ZTF data

The ZTF survey scans the entire northern sky with two-day cadence observations, using the 48-inch Samuel Oschin Schmidt telescope with the ZTF g, r, and i filters (e.g., Bellm et al. 2019). We obtained five (for ZTF g) and 42 (for ZTF r) photometric data without bad-data flags (i.e., catflags |$=0$|⁠) from the ZTF Data Release 20. The epoch of those data covers the observation period from MJD 58430.5 to 59271 (or MJD 58439.4 to 59638.2) for the ZTF g (or r) band. The ZTF magnitudes were converted to those of Pan-STARRS by Medford, Lu, and Schafly (2020).

2.2 Monitoring observations

From 2021 February to 2022 March we performed |$g^{\,\prime}$| (⁠|$4770\,$|Å), |$R_{\rm c}$| (⁠|$6492\,$|Å), and |$I_{\rm c}$| (⁠|$8020\,$|Å) bands follow-up observations (⁠|$\sim \!1\:$|d to one week cadence) of WISE J0909|$+$|0002 with the 105|$\:$|cm Murikabushi telescope/MITSuME. From 2023 May to 2024 April we observed the same object with the Okayama and Akeno 50|$\:$|cm telescopes/MITSuME with the |$g^{\prime }$|⁠, |$R_{\rm c}$|⁠, and |$I_{\rm c}$| bands (Kotani et al. 2005; Yatsu et al. 2007; Shimokawabe et al. 2008) and the SaCRA 55|$\:$|cm telescope/MuSaSHI (Oasa et al. 2020) with r, i, and z bands as the OISTER collaboration to confirm the continuous periodicity flux variability of our target. The magnitudes of WISE J0909|$+$|0002 obtained from those telescopes and instruments were converted to the AB magnitude based on the bandpass transformations of Tonry et al. (2012) and the relation between the Vega and AB magnitude (Blanton & Roweis 2007). We used photometric data of the signal-to-noise ratio of |$\gtrsim \!10$| taken with those telescopes and removed any data with a wrong weather condition or astronomical seeing.

3 Data analysis and results

3.1 Light curves of WISE J0909|$+$|0002

First, we applied sinusoidal models to the light curves of WISE J0909|$+$|0002 as follows:

$$\begin{eqnarray} F(t) = A\sin \left(\frac{2\pi t}{P_{\rm Fit}} + \phi \right) + b, \end{eqnarray}$$

(1)

where A, |$P_{\rm Fit}$|⁠, |$\phi$|⁠, and b are the best-fitting amplitude, period, phase, and average magnitude, respectively. The Python package, |${\tt curve\_fit}$|⁠, was used for calculations of these parameters in each band. The y- and |$R_{\rm c}$|-band parameters were omitted from the best-fitting results, since the amplitude of these two bands were consistent with those errors (i.e., |$A \lt 2\sigma$| uncertainty). Secondly, we combined all bands’ light curves into one, since our light curves include a sparse period from MJD 57021 to 58430. In order to connect those light curves, the best-fitting amplitudes and median values of each band were adjusted with respect to the CRTS V-band data; hereafter we define this curve as the combined light curve. Finally, we calculated the signal-to-noise ratio (SNR), |$\xi =A^2/(2\sigma _{\rm res}^2)$|⁠, of our light curves, where |$\sigma _{\rm res}^2$| is the variance of the residues between the sinusoidal curve model and each data point (Horne & Baliunas 1986). If |$\xi \gt \!0.5$|⁠, the periodical signal is significantly larger than the residual noise (see Chen et al. 2020).

Figures 2 and 3 show each band light curve and combined light curve, respectively. The z- and |$I_{\rm c}$|-band light curves exhibited very different trends from the other curves. We note that the difference in the phase |$\phi$| in each band is not a real value but just fitting results, since there is the disparity in observation epochs and cadences between each band. Table 2 lists the parameters of all the bands’ light curves and the combined light curve. The SNRs indicate significant periodic signals in all cases (i.e., |$\xi \gt \!0.5$|⁠).

Fig. 2.

Light curves of WISE J0909|$+$|0002 with the CRTS V, Pan-STARRS (⁠|$grizy$|⁠), ZTF (g and r), Murikabushi telescope, Akeno and Okayama telescope/MITSuME (g, |$R_{\rm c}$|⁠, and |$I_{\rm c}$|⁠), and SaCRA/MuSaSHI (r, i, and z). The solid sinusoidal curves present the best-fitting results for these light curves. The shadowed region indicates |$1\sigma$| errors in the variability amplitude of the light curves. The y- and |$R_{\rm c}$|-band fitting results are not displayed, since their parameters are statistically insufficient.

Fig. 3.

Combined light curve of WISE J0909|$+$|0002. The symbols in this figure are the same as those in figure 2.

Table 2.

Sinusoidal-model parameters for the light curves.

	A*	|$P_{\rm Fit}$|†		b§
Band	(mag)	(d)	|$\phi$|‡	(mag)	|$\xi$|‖
CRTS V	0.064 ± 0.005	654.09 ± 15.24	−119.97 ± 4.26	16.449 ± 0.003	0.59
g	0.085 ± 0.011	661.80 ± 12.80	−117.65 ± 3.77	16.727 ± 0.006	1.62
r	0.077 ± 0.010	648.09 ± 16.31	−121.91 ± 4.99	16.620 ± 0.006	2.76
i	0.087 ± 0.005	686.28 ± 7.10	−111.10 ± 1.90	16.400 ± 0.004	18.66
|$I_{\rm c}$|	0.067 ± 0.046	378.56 ± 40.58	−165.89 ± 41.44	16.464 ± 0.034	0.84
z	0.084 ± 0.036	1232.69 ± 74.26	−31.82 ± 6.05	16.352 ± 0.024	2.28
Combined light curve	0.061 ± 0.004	666.03 ± 3.24	−116.71 ± 0.93	16.453 ± 0.002	0.68

*

Variability amplitude.

†

Rest-frame periodicity.

‡

Phase.

§

Average magnitude.

‖

SNR, |$A^2/(2\sigma _{\rm res}^2)$|⁠.

3.2 Lomb–Scargle periodgram

To evaluate the rest-frame periods of the combined and individual bands’ light curves, we calculated the power spectrum of the archival and our observational light curves using the Lomb–Scargle method (Lomb 1976; Scargle 1982; VanderPlas 2018) with its corresponding package, |${\tt LombScargle}$| from |${\tt astropy}$|⁠. This method has the advantage of being able to process signals sampled at indefinite intervals and with missing samples. Whether the obtained power peaks are pure noise or not was determined by estimating the false alarm probability (FAP) with the bootstrap method. In this study, the periodicity is significant if the following conditions are satisfied: (i) the FAP at a power peak is significantly small (⁠|$\lt \!1.0\%$|⁠), (ii) the monitoring length of each band data is longer than their periods, and (iii) the time interval without observations is shorter than one period of the light curve. It was found that the CRTS V-, g-, and r-band light curves (hereafter, three-band light curves) satisfied criteria (i), (ii), and (iii).

Figure 4 presents the Lomb–Scargle power of the combined light curve and three-band light curves. Despite the different observation epochs of the three-band light curves, their rest-frame periods (⁠|$\sim \!660$|–|$689\:$|d) were consistent with each other. The FAPs at the peak of the combined and three-band light curves (⁠|$\lt \!0.001$| for all cases) indicate that those light curves are significantly periodic. Therefore the light curves of WISE J0909|$+$|0002 obtained from MJD 53494.1 (CRTS V) to MJD 60419.5 (SaCRA/MuSaSHI r, i, and z) had likely continued the periodic variability for at least |$\sim \!6.6\:$|yr in the quasar rest frame, corresponding to |$\sim \!3.6$| cycles.

Fig. 4.

Lomb–Scargle periodogram of (a) the combined light curve, (b) CRTS V-, (c) g-, and (d) r-band light curves including the PanSTARRS, ZTF, and our observation data. The vertical lines present the rest-frame day of their power peaks.

3.3 Auto-correlation analysis

We performed an auto-correlation function (ACF) analysis, which evaluates the correlation between a signal and its delayed replica, so as to verify the periodicity of the combined light curve. Since the combined light curve is unevenly sampled time-series data, we carried out a z-transformed discrete correlation function (ZDCF; Alexander 1997) analysis by using |${\tt pyzdcf}$|⁠;¹ this method is applicable to uneven sampling data. The ACF shall not be applied to three-band light curves that are sparser than the combined light curve.

If a signal is periodic, its ACF is known to decay with periodicity (Jung 1993). We applied an exponentially decaying cosine model to the ACF:

$$\begin{eqnarray} {\rm ACF}(\tau )=A\exp (-\lambda \tau )\cos (\omega \tau ), \end{eqnarray}$$

(2)

where A, |$\lambda$|⁠, |$\tau$|⁠, and |$\omega$| are the amplitude, decay rate, time lag, and angular frequency, respectively (see also Chen et al. 2020). Figure 5 shows the ACF with the ZDCF method. In the combined light curve, the error rate between the ACF period (⁠|$2\pi /\omega =659.6\:$|d) and the period estimated from the sinusoidal model in table 2 (or the Lomb–Scargle periodogram: |$689.4\:$|d) is |$1.0\%$| (or |$4.3\%$|⁠); the ACF analysis exhibited a periodicity which is comparable to the previous sections.

Fig. 5.

ACF of the combined light curve (solid line) and the best-fitting decayed cosine model (dashed line) with the ZDCF method.

3.4 The DRW parameters and structure functions

While the periodicities of the three-band light curves determined by the Lomb–Scargle method are consistent with each other, the possibility with which they indicate false-positive QPOs must be excluded. The DRW process, which is well used to model AGN light curves, is known as an aperiodic stochastic process (e.g., Kelly et al. 2009; Kozłowski et al. 2010; MacLeod et al. 2010). If the light curves of our target (figures 2 and 3) are false-positive QPOs, they would behave similarly to the DRW or other processes. Here we examine whether the three-band light curves present the DRW trend. It should be noted that the DRW is not necessarily a fundamental process underlying AGN variability (Kozłowski 2016).

The DRW model is well described by the first-order continuous autoregressive [CAR(1)] process. In this process, the components of the covariance matrix S is expressed as:

$$\begin{eqnarray} S_{ij}=\frac{\tau \sigma ^2}{2}\exp \left({-\frac{|t_i-t_j|}{\tau }}\right), \end{eqnarray}$$

(3)

where |$t_i-t_j$|⁠, |$\tau$|⁠, and |$\sigma$| are the time separation, damping timescale, and variability amplitude, respectively. In order to evaluate the parameters |$\sigma$| and |$\tau$|⁠, we utilized the JAVELIN (Just Another Vehicle for Estimating Lags In Nuclei; Zu et al. 2011) code,² which employs the Markov Chain Monte Carlo (MCMC; Foreman-Mackey et al. 2013) method for determining posterior distributions of time lags from the reverberation mapping and physical parameters, assuming the DRW process.

In addition to the above analysis, we estimated the structure functions (e.g., di Clemente et al. 1996; Vanden Berk et al. 2004; de Vries et al. 2005; Wilhite et al. 2008; MacLeod et al. 2010; Sun et al. 2018; De Cicco et al. 2022) for the three-band light curves expressed as:

$$\begin{eqnarray} SF(\Delta t)=\sqrt{\frac{1}{N(\Delta t)}\sum _{i\lt j}(\Delta m_{ij})^2}, \end{eqnarray}$$

(4)

where |$N(\Delta t)$| is the number of data points in a time-separation bin |$\Delta t\, [=|t_i-t_j|/(1+z)]$|⁠, and |$\Delta m_{ij}$| is the magnitude difference between two epochs. We adopted the structure function definition given by de Vries et al. (2005), which can be evaluated even if the photometric uncertainties are larger than the mean variability amplitude. If the behavior of |$SF$| follows (or does not follow) the DRW process, it would present asymptotic (or jiggly) distribution with increasing time. Alternatively, the light curves of WISE J0909|$+$|0002 are difficult to express through the DRW process, if the asymptotic values of structure functions |$SF_{\infty }=\sigma \sqrt{\tau }$| significantly deviate from that of normal quasars (see MacLeod et al. 2010).

Figure 6 displays the two-dimensional posterior distributions of |$\sigma$| and |$\tau$| and the structure functions of the combined light curve and three-band light curves. Using those parameters, we estimated the asymptotic value, |$SF_{\infty }$|⁠, and compared our results with the values in figure 3 of MacLeod et al. (2010). Table 3 summarizes |$\sigma$|⁠, |$\tau$|⁠, and |$SF_{\infty }$|⁠. The asymptotic values |$\log SF_{\infty }$| of the three-band light curves significantly deviate from those listed in MacLeod et al. (2010): |$\log SF_{\infty }\sim \!-0.70$|⁠. While |$\log SF_{\infty }$| of the combined light curve (⁠|$\log SF_{\infty }\sim \!-0.57$|⁠) was consistent with that of radio-loud quasars in MacLeod et al. (2010), the |$\tau$|–|$\sigma$| distribution differed from their results (figure 6a). MacLeod et al. (2010) also calculated |$SF_{\infty }$| for quasars with significant periodicity (see figure 18 of that paper); the |$SF_{\infty }$| values in this study were comparable to that study. Moreover, the structure functions exhibited the jiggly distribution against the trend of the DRW process (figure 6e). At least in the above analysis, the combined light curve and three-band light curves unlikely showed false-positive QPOs.

Fig. 6.

Two-dimension posterior distributions of variability amplitude |$\sigma$| and damping timescale |$\tau$| for (a) combined light curve, (b) CRTS V-, (c) g-, and (d) r-band light curves. The brighter (or dimmer) contours present denser (or thinner) regions of these distributions. The circle and triangle indicate no-radio and radio-loudest quasar subsamples, respectively (table 2 in MacLeod et al. 2010). Panel (e) shows the structure functions of the combined light curve and the three-band light curves (inverse triangle: combined light curve, filled circles: CRTS V band, triangles: g band, and squares: r band).

Table 3.

Variability amplitude, damping timescale, and |$SF_{\infty }$| of WISE J0909|$+$|0002.

	|$\sigma$|	|$\tau$|	|$\log SF_{\infty }$|*
	(mag|$\:$|d|$^{-1/2}$|⁠)	(d)	(mag)
Combined light curve	0.050|$^{+0.005}_{-0.005}$|	28.74|$^{+13.24}_{-11.94}$|	|$-0.571^{+0.123}_{-0.163}$|
CRTS V	0.043|$^{+0.013}_{-0.008}$|	420.97|$^{+424.01}_{-193.54}$|	|$-0.054^{+0.158}_{-0.223}$|
g	0.068|$^{+0.019}_{-0.014}$|	354.96|$^{+331.50}_{-155.72}$|	0.107|$^{+0.250}_{-0.226}$|
r	0.067|$^{+0.026}_{-0.014}$|	435.22|$^{+546.31}_{-234.49}$|	0.145|$^{+0.319}_{-0.270}$|
No radio†	0.0163 ± 0.0002	201.83 ± 2.32	−0.634 ± 0.003
Radio-loud†	0.019 ± 0.002	199.52 ± 22.97	−0.57 ± 0.02

*

The asymptotic value of a structure function equal to |$\sigma \sqrt{\tau }$|⁠.

†

Non-radio and radio-loudest quasar subsamples from table 2 in MacLeod et al. (2010).

3.5 Simulations of light curves with single power-law and damped random walk

Pure red noise such as the single power-law (SPL) and DRW processes can produce quasi-periodic light curves over a few cycles (Vaughan et al. 2016). We generated mock light curves with the same baseline as the combined light curve of the SPL and DRW simulations, using |${\tt astroML}$|⁠,³ so as to compare the sinusoidal model with the pure red noise. As shown in table 4, we adjusted the power-law index, |$\alpha$|⁠, of the SPL model in the range from |$-0.5$| to |$-2.0$|⁠. The DRW parameters, |$\tau$| and |$SF_{\infty }$|⁠, were set referring to the results in table 3. Since simulated-DRW curves with |$\tau \lesssim 100\:$|d exhibited extremely large amplitude (e.g., |$10\:$|mag), we restricted the |$\tau$| range to be from 200 to |$600\:$|d. 10000 light curves were generated for each SPL or DRW parameter combination in table 4.

Table 4.

Parameters of the SPL and DRW processes and the model comparison with the sinusoidal model and red noise.

Figure			|$\tau$|	|$SF_{\infty }$|
number	|$\alpha$|*	|$\Delta \mbox{BIC}_{\rm Data-SPL}$|†	(d)	(mag)	|$\Delta \mbox{BIC}_{\rm Data-DRW}$|‡
7a	−0.5	−138.15	200	0.1	−293.33
7b	−0.7	−274.63	200	0.2	−1234.78
7c	−1.0	−581.73	400	0.1	101.72
7d	−1.2	−1476.40	400	0.2	−2080.08
7e	−1.5	−1910.87	600	0.1	313.03
7f	−2.0	−2615.59	600	0.2	−1747.82

*

Power-law index.

†

The maximum value of BIC|$_{\rm sin,data}- \mbox{BIC}_{\rm sin,SPL}$| in the 10000 simulations.

‡

The maximum value of BIC|$_{\rm sin,data}- \mbox{BIC}_{\rm sin,DRW}$| in the 10000 simulations.

In order to compare the suitability to the sinusoidal model (black solid line in figure 3) of the combined light curve and the SPL and DRW model curves, we evaluated the Bayesian information criterion (BIC)⁴ under the Gaussian error model:

$$\begin{eqnarray} {\rm BIC} = N\ln \left[\frac{\sum _{i=1}^N(y_i - f_{\rm i, sinusoidal})^2}{N} \right] +k\ln (N), \end{eqnarray}$$

(5)

where N is the number of data, |$y_i$| is the ith observed or simulated data, |$f_{\rm i, sinusoidal}$| is the ith model flux, and k is the number of free parameters. We estimated the BIC differences |$\Delta \mbox{BIC}_{\rm Data-SPL}\ (= \mbox{BIC}_{\rm sin,data}- \mbox{BIC}_{\rm sin,SPL}$|⁠) and |$\Delta \mbox{BIC}_{\rm Data-DRW}\ (= \mbox{BIC}_{\rm sin,data}- \mbox{BIC}_{\rm sin,DRW}$|⁠), where BIC|$_{\rm sin,data}$|⁠, BIC|$_{\rm sin,SPL}$|⁠, and BIC|$_{\rm sin,DRW}$| are the BIC between the best-fitting sinusoidal model and the combined light curve, mock SPL, and DRW curves, respectively. If |$\Delta \mbox{BIC}$| is less than |$-10$| (e.g., Chen et al. 2020; Liao et al. 2021), the combined light curve is more suitable with the sinusoidal model than mock-light curves by the SPL or DRW simulations.

As a result, |$\Delta \mbox{BIC}$| was below |$-10$| in most cases. It was difficult to reproduce the light curve close to the best-fitting sinusoidal model with the SPL simulations. However, in the DRW simulations with |$(\tau, SF_{\infty })=(400, 0.1)$| and |$(600, 0.1)$|⁠, |$\Delta \mbox{BIC}_{\rm Data-DRW}$| exceeded |$-10$| on rare occasions with a probability of |$0.02\%$| [for |$(400, 0.1)$|] and |$0.08\%$| [for |$(600, 0.1)$|]. Figure 7 displays the best-fitting sinusoidal curve and the mock SPL and DRW curves, which showed the maximum |$\Delta \mbox{BIC}$| in the 10000 simulations. For instance, the DRW curve in figure 7e presented a sinusoidal-like feature. Therefore, it is difficult to completely deny that the light curve of WISE J0909|$+$|0002 is the pure red noise in our simulations. More long-term monitoring observations would bring us more robust validations for QPOs in WISE J0909|$+$|0002.

Fig. 7.

Mock light curves with the SPL (dashed lines) and DRW models (dash–dotted lines), showing the maximum |$\Delta \mbox{BIC}$| in the 10000 simulations. The solid line is the mean-subtracted (best-fitting) sinusoidal model for the combined light curve. See table 4 for the SPL and DRW parameters corresponding to each figure.

3.6 Power spectral density analysis

We confirm the periodicity of the combined light curve by the power spectral density (PSD) approach. The PSD is useful tool to determine whether signals are noisy or periodic. We calculated the PSD, |$P(f)$|⁠, for the combined light curve according to the following equation:

$$\begin{eqnarray} P(f_i) = \frac{2\Delta \tau }{N} |F(f_i)|^2, \end{eqnarray}$$

(6)

and

$$\begin{eqnarray} |F(f_i)|^2=\biggr | \sum ^N_{j=1}x_j\cos (2\pi f_i t_j) \biggr |^2 + \biggr | \sum ^N_{j=1} x_j\sin (2\pi f_i t_j) \biggr |^2, \end{eqnarray}$$

(7)

where N, |$\Delta t=(t_N-t_1)/N$|⁠, |$x_j$|⁠, and |$f_i=i/(N\Delta t)$| are the number of data, the rest-frame time separation, the mean-subtracted flux at time |$t_j$|⁠, and the ith frequency with |$i=1$|⁠, 2, ..., |$N/2$|⁠, respectively (e.g., Uttley et al. 2002; Li et al. 2019). For modeling the PSD, we considered the (i) SPL (aperiodic)

$$\begin{eqnarray} P_{\rm SPL}(f) = Af^{-\alpha }, \end{eqnarray}$$

(8)

(ii) DRW (aperiodic)

$$\begin{eqnarray} P_{\rm DRW}(f) = \frac{A}{1+(f/f_{\rm b})^2}, \end{eqnarray}$$

(9)

and (iii) periodic with aperiodic models

$$\begin{eqnarray} P_{\rm P+AP}(f) = \frac{A_{\rm p}}{\sqrt{2\pi }\omega _{\rm p}}\exp \left[{-\frac{\left(f-f_{\rm p}\right)^2}{2\omega _{\rm p}^2}}\right]+P_{\rm AP}, \end{eqnarray}$$

(10)

where A, |$\alpha$|⁠, |$f_{\rm b}$|⁠, |$A_{\rm p}$|⁠, |$f_{\rm p}$|⁠, and |$\omega _{\rm p}$| are free parameters, and |$P_{\rm AP}$| is |$P_{\rm SPL}(f)$| or |$P_{\rm DRW}(f)$|⁠. For model selection, we calculated the BIC difference, |$\Delta \mbox{BIC}$|⁠, with respect to the SPL model.

Figure 8 presents the PSD of the combined light curve with the SPL, the DRW, and both SPL and DRW with the periodic model. Although the PSD distribution presented a peak at |$802.2\:$|d, the sparse sampling of the combined light curve resulted in a deviation from the period obtained in subsections 3.1–3.3 (i.e., 660 to |$689\:$|d). Table 5 lists the free parameters of the PSD models. In table 5, the BIC difference, |$\Delta \mbox{BIC}$|⁠, between the SPL and the other models indicates that the SPL |$+$| periodic model is the most preferable of the above models.

Fig. 8.

PSDs of the combined light curve (solid lines). Dashed lines present the best-fitting lines for (a) the SPL, (b) DRW, (c) SPL |$+$| periodic, and (d) DRW |$+$| periodic model. Shadowed areas indicate the |$1\sigma$| error regions for each model. Vertical shadowed regions show the frequency corresponding to the 660 to |$689\:$|d period obtained in subsections 3.1–3.3.

Table 5.

Model parameters of the PSD.

			|$\log f_{\rm b}$|		|$\log f_{\rm p}$|	|$\log \omega _{\rm p}$|
Model	|$\log A$|	|$\alpha$|	(d|$^{-1}$|⁠)	|$\log A_{\rm p}$|	(d|$^{-1}$|⁠)	(d|$^{-1}$|⁠)	|$\Delta \mbox{BIC}$|*
SPL	2.44 ± 0.17	−0.67 ± 0.06	—	—	—	—	0
DRW	4.76 ± 0.05	—	−2.66 ± 0.07	—	—	—	−46.04
SPL|$+$|Periodic	2.71 ± 0.07	−0.40 ± 0.04	—	1.93 ± 0.01	−2.930 ± 0.004	−3.74 ± 0.02	−945.62
DRW|$+$|Periodic	3.56 ± 0.04	—	−1.30 ± 0.06	1.95 ± 0.01	−2.928 ± 0.003	−3.716 ± 0.003	−917.11

*

The BIC difference relative to the SPL model.

4 Discussion

We found that the light curves of WISE J0909|$+$|0002 were likely to show multi-year QPOs in the previous section. The period of the combined light curves (⁠|$\sim \!660$|–|$689\:$|d in the rest frame) is somewhat longer compared with periodic quasars listed in table 2 of Charisi et al. (2016). After that, we discuss the validity of the following scenarios as mechanisms of QPOs: the Doppler boost, circumbinary disk, and radio-jet precession.

4.1 Doppler boost model

4.1.1 Variability amplitude ratio

Here, we examine the relativistic Doppler boost model qualitatively by using multicolor light curves (e.g., D’Orazio et al. 2015; Charisi et al. 2018) taken with the CRTS, Pan-STARRS, and the OISTER project. In the Doppler boost model, photon frequencies are modulated by the factor |$D= [\gamma (1-\beta _{\parallel })]^{-1}$|⁠, where |$\gamma =({1-\beta ^2})^{-1/2}$|⁠, |$\beta =v/c$|⁠, and |$\beta _{\parallel }$| are the Lorentz factor, the three-dimensional velocity v in the units of the speed of light c, and velocity along our line of sight, respectively. The velocity component, |$\beta _{\parallel }$|⁠, is equal to |$\beta \cos {\varphi }\sin {i}$|⁠, where |$\varphi$| is the orbital phase and i is the orbital inclination. The spacial photon density is the Lorentz invariant and is proportional to |$F_{\nu }/\nu ^3$|⁠, where |$F_{\nu }$| is the apparent flux as a function of the frequency |$\nu$|⁠. Namely, the flux |$F_{\nu }$| is expressed as |$F_{\nu }=D^{3-\alpha _{\nu }}F^0_{\nu }$|⁠. Assuming that intrinsic flux |$F^0_{\nu }$| obeys the power law (⁠|$F^0_{\nu }\propto \nu ^{\alpha _\nu }$|⁠), the variability amplitude to first order in β is described as follows:

$$\begin{eqnarray} \frac{\Delta F_{\nu }}{F_{\nu }}=(3-\alpha _{\nu })\beta \cos {\varphi }\sin {i}. \end{eqnarray}$$

(11)

For instance, the variability amplitude in the CRTS V band is |$0.064\:$|mag, which corresponds to |$\Delta F_{\nu }/F_{\nu }=0.064$|⁠. We adopt the phase |$\cos {\varphi }=1$|⁠, since our observation timescale is far shorter than the coalescence timescale of BSBHs. From equation (11), the ratio of variability amplitude between two bands is written as follows:

$$\begin{eqnarray} \frac{A_s}{A_l}=\frac{3-\alpha _{\nu , s}}{3-\alpha _{\nu , l}}, \end{eqnarray}$$

(12)

where A is the variability amplitude of a sinusoidal curve model. The symbol s (or l) is a value for the shorter (or longer) wavelength.

Using the SDSS DR17 (Abdurro’uf et al. 2022) spectrum data, we estimated the power-law index in each band. The SDSS performed two-epoch spectroscopic observations on MJD 51929 and MJD 55532 for WISE J0909|$+$|0002. We estimated the power-law slopes of those SDSS spectra to verify the Doppler boost scenario by assuming constant variability amplitude in each band throughout the observation epochs. Table 6 summarizes the power-law slopes in each band. Figure 9 displays the SDSS spectrum of WISE J0909|$+$|0002 and the ratios of variability amplitude and power-law slopes in each band combination.

Fig. 9.

Panels (a) and (c): SDSS spectrum of WISE J0909|$+$|0002 (black solid line). The dotted line presents the model spectrum without emission lines. The power-law curves of the g-, CRTS V-, and r-band ranges are shown as solid lines with their filter names. Characteristic emission lines are labeled. Panels (b) and (d): Relation between the variability amplitude ratio and the power-law index ratio from equation (12). The theoretical line based on the Doppler boost scenario is shown with the black solid line.

Table 6.

Spectral power-law slopes of WISE J0909|$+$|0002.

	g band	CRTS V	r band
MJD 55532
|$\alpha _{\rm \nu }$|*	|$-0.558$| ± 0.027	|$-0.502$| ± 0.011	|$-0.966$| ± 0.017
MJD 51929
|$\alpha _{\rm \nu }$|*	|$-0.764$| ± 0.016	|$-0.583$| ± 0.013	|$-0.111$| ± 0.021

*

Power-law index of flux as a function of frequency (⁠|$=-\beta _{\rm \lambda }-2$|⁠), where |$\beta _{\lambda }$| is the power of the flux as a function of wavelength.

Three of six cases in figures 9b (⁠|$g/r$| and |$V/r$|⁠) and 9d (⁠|$g/r$|⁠) were consistent with the theoretical prediction of the Doppler boost model within a |$1\sigma$| accuracy. Thus the scenario is somewhat positive from the aspect of the variability amplitude ratios.

According to the results of Charisi et al. (2018), only |$\sim \!20\%$| (or |$\sim \!37\%$|⁠) of quasars with periodic light curves showed the Doppler boost signature in the NUV (or FUV) region. Integrating the results of Chen et al. (2020) and Liao et al. (2021), one out of five periodic quasars showed evidence supporting the Doppler boost scenario by the variability amplitude ratios. Thus, we found that WISE J0909|$+$|0002 is likely a rare QPO sample, showing the Doppler boost trend. However, whether the Doppler boost theory is supportive depends on the observation wavelength as reported in the Swift observations for PG 1302|$-$|102 (Xin et al. 2020): additional X-ray and/or MIR (e.g., Jun et al. 2015) observations would be needed for detailed verifications of this scenario.

4.1.2 Parameter space of inclination angle and black hole mass

We also examine the Doppler boost model by exploring parameter spaces of the black hole mass |$M_{\rm BH}$|⁠, the binary black hole mass ratio |$q =M_2/M_1 (\le 1)$|⁠, and the orbital inclination i (cf. D’Orazio et al. 2015; Chen et al. 2020; Liao et al. 2021; Song et al. 2021). In a binary black hole, the velocity of a secondary disk is written using q:

$$\begin{eqnarray} v_2=\left(\frac{1}{1+q}\right)\left(\frac{2\pi GM_{\rm BH}}{P}\right)^{1/3}, \end{eqnarray}$$

(13)

where G is the gravitational constant and P is the orbital period of a BSBH system. Using equations (11) and (13) with |$\cos {\varphi }=1$| and a primary black hole velocity |$v_1=-qv_2$|⁠, the inclination angle i can be expressed as

$$\begin{eqnarray} i&=& \arcsin {\left[(1+q)\left(\frac{2\pi GM_{\rm BH}}{P}\right)^{-1/3}\frac{c}{3-\alpha _\nu }\frac{\Delta F_{\nu }}{F_{\nu }}\frac{1}{f} \right] }, \end{eqnarray}$$

(14)

and

$$\begin{eqnarray} f&=& f_2-q(1-f_2), \end{eqnarray}$$

(15)

where |$f_2$| is the ratio of luminosity contribution from the secondary black hole. We selected the following mass ratios: |$q=0$| (extreme case), 0.11, 0.43 (the characteristic values with two power peaks of mass accretion rates in a hydrodynamical simulations; Farris et al. 2014), and 1.0 (extreme case). Figure 10 presents the |$M_{\rm BH}$|⁠, i, and q parameter space of the Doppler boost theory for each band (see also Chen et al. 2020; Liao et al. 2021). The inclination |$i=90^\circ$| (or |$i=0$|⁠) indicates an edge-on (or a face-on) view to accretion disks. Under the conditions of the black hole mass of |$7.4\times 10^9\, M_{\odot }$| and the mass ratio range of |$0.11\le q \le 0.43$| at |$f_2=0.8$|⁠, the inclination angle of WISE J0909|$+$|0002 allows |$i \sim \!10^\circ$| (i.e., near the face-on angle).

Fig. 10.

Parameter space of the black hole mass, orbital inclination, and black hole mass ratio for the (a) CRTS V, (b) g, and (c) r bands. The top, second, third, and bottom curves correspond to mass ratios of |$q=0, 0.11, 0.43$|⁠, and 1.0, respectively. The dashed (or solid) curves show the case where the secondary-disk luminosity contribution, |$f_2$|⁠, is 1.0 (or 0.8). Shadowed regions indicate the ranges of the inclination angle such that the ranges of the mass ratio become |$0.11\le q \le 0.43$| at |$f_2=0.8$|⁠. The vertical solid line and shadows represent the black hole mass of WISE J0909|$+$|0002 (⁠|$= 7.4\times 10^{9}\, M_{\odot }$|⁠) and its |$1\sigma$| error range, respectively.

WISE J0909|$+$|0002 is very unusual since it is not only an ELIRG with likely periodic variability, but also contains a broad absorption line (BAL; FWHM |$\ge 2000\:$|km|$\:$|s|$^{-1}$|⁠; Weymann et al. 1991) quasar (Moravec et al. 2017), which is thought to eject outflow gas from a viewing angle closer to the edge-on of the accretion disk (e.g., Proga & Kallman 2004; Nomura et al. 2013). Assuming that the BAL quasar fraction (⁠|$\sim \!10\%$|–|$40\%$|⁠; Weymann et al. 1991; Reichard et al. 2003; Gibson et al. 2009; Allen et al. 2011) reflects the outflowing angle from the edge-on view (i.e., |$i \gtrsim \!90^{\circ }-40^{\circ }=50^{\circ }$|⁠), our result (⁠|$i \sim \!10^\circ$|⁠) has a range of considerably smaller angles than the above physical picture of BAL quasars. Basing on our analysis for the relativistic boost hypothesis, the inclination range of WISE J0909|$+$|0002 permits the claim that states outflows bringing BALs can be ejected from intermediate or close to face-on angles (Matthews et al. 2017). Namely, the physical picture of inclination angles based on the detection rate of BALs is inconsistent if |$f_2\gt \!0.8$|⁠. WISE J0909|$+$|0002 probably ejected outflows near the face-on angle in the process of evolution to be a BAL quasar (e.g., Farrah et al. 2007; Lipari et al. 2009; Bischetti et al. 2023) and/or an ELIRG under a high |$f_2$| ratio.

4.2 Circumbinary disk model

As explained in the previous studies, the circumbinary accretion disk (CBD) model is one of the plausible candidates of periodic variability in AGNs. This model consists of the radiation from the CBD, primary, and secondary minidisks and predicts a UV–optical–IR spectral cut-off, which indicates gapped or truncated binary disk structures (Roedig et al. 2014; Farris et al. 2015). In Guo et al. (2020), six out of the 138 candidates of periodic quasars identified by Graham et al. (2015a) and Charisi et al. (2016) exhibited the spectral cut-offs predicted by the CBD model in the UV–optical–IR SEDs. The cut-off temperature is described as

$$\begin{eqnarray} T_{\rm cut} \sim 2.0\times 10^4 \left[\dot{m}\left(\frac{\eta }{0.1}\right)^{-1}\left(\frac{M_{\rm BH}}{10^{8}\, M_{\odot }}\right)^{-1}\left(\frac{a}{10^2\, R_{\rm g}}\right)^{-3}\right]^{1/4}, \nonumber\\ \end{eqnarray}$$

(16)

where |$\dot{m}$|⁠, |$\eta$|⁠, a, and |$R_{\rm g}$|  |$(= GM_{\rm BH}/c^2)$| are the accretion rate in Eddington units, radiative efficiency, the semi-major axis of a binary black hole, and the gravitational radius, respectively. The minimum value of a notch appears around |$T_{\rm notch}$|⁠, where |$T_{\rm notch}\sim \!2^{3/4}T_{\rm cut}$|⁠. Table 7 lists the physical parameters for considering the CBD model with respect to the spectral energy distributions (SEDs) of our target.

Table 7.

Physical properties of WISE J0909|$+$|0002.

|$L_{\rm bol}$|*	|$M_{\rm BH}$|†			a‖	|$T_{\rm cut}$|#	|$\lambda _{\rm cut}$|**	|$T_{\rm notch}$|††	|$\lambda _{\rm notch}$|‡‡
(⁠|$10^{47}\:$|erg|$\:$|s|$^{-1}$|⁠)	(⁠|$10^9\, M_{\odot }$|⁠)	|$\dot{m}$|‡	|$\eta$|§	(ld)	(⁠|$10^3\:$|K)	(⁠|$10^3\,$|Å)	(⁠|$10^4\:$|K)	(⁠|$10^3\,$|Å)
4.3 ± 0.6	7.4 ± 0.3	0.4 ± 0.1	0.1	29.0 ± 0.4	7.2 ± 0.5	4.0 ± 0.3	1.2 ± 0.1	2.4 ± 0.1

*

Bolometric luminosity (Toba et al. 2021).

†

Blackhole mass (Toba et al. 2021).

‡

Eddington ratio (Toba et al. 2021).

§

Radiative efficiency.

‖

Semi-major axis of a binary black hole estimated from Song et al. (2021).

#

Cut-off temperature.

**

Cut-off wavelength based on Wien’s law.

††

Notch temperature.

‡‡

Notch wavelength based on Wien’s law.

We obtained the SED (in units of mJy) of WISE J0909|$+$|0002 (Toba et al. 2021) and the composite spectra of 259 type 1 quasars (Richards et al. 2006). The galactic extinction is corrected for the SED by Toba et al. (2021). Here, we employed the K-corrected SED to convert the flux into luminosity. Figure 11 displays the normalized SED of WISE J0909|$+$|0002. According to equations (6), (7), and (8) of Roedig, Krolik, and Miller (2014), we applied the CBD model to the infrared-to-FUV flux. From figure 11a, the SED of WISE J0909|$+$|0002 exhibited a different trend from the CBD model for both |$q=0.11$| and 0.43 (Farris et al. 2014) at |$f_2=0.8$|⁠: the cut-off and notch predicted in the CBD model are not shown in our target. This is similar to the result that the periodic quasars PSO J334.2028|$+$|01.4075 and SDSS J025214.67|$-$|002813.7 did not clearly exhibit the CBD feature (Foord et al. 2017, 2022).

Fig. 11.

Rest-frame SED of WISE J0909|$+$|0002 (squares) obtained from Toba et al. (2021). Panel (a): SED normalized by the rest-frame |$3550\,$|Å flux. The solid line with emission lines indicates the rest-frame UV flux obtained from the SDSS DR18 catalog. Solid and dash-dotted curves indicate the CBD model flux with mass ratios, |$q=0.11$| and 0.43, respectively. The flux ratio of the secondary accretion disk, |$f_2$|⁠, is fixed at 0.8. The left (or right) vertical line shows the positions of the cut-off (or notch) wavelength that are expected from the CBD model. Panel (b): Flux of WISE J0909|$+$|0002, and SDSS J025214.67|$-$|002813.7 (diamonds, Liao et al. 2021; Foord et al. 2022) normalized by each |$\sim \!10^{14.93}\:$|Hz flux. The solid, dashed, and dash-dotted lines denote the composite spectra of all quasars (“R06-All”), infrared luminous quasars (“R06-IL”), and infrared dim quasars (“R06-ID”), respectively (Richards et al. 2006).

In addition to the above argument, the NUV and FUV SEDs of WISE J0909|$+$|0002 (⁠|$\sim \!10^{15.5}\:$|Hz flux in figure 11b) have a steep deficit, which plausibly originates from a BSBH accretion mode (Yan et al. 2015) or a reddened flux (Leighly et al. 2014) as discussed in Mrk 231. Liao et al. (2021) found that the quasar SDSS J025214.67|$-$|002813.7 presented |$\sim \!1.2\:$|dex dropping off near |$\sim \!1400\,$|Å (figure 11b). Then Foord et al. (2022) concluded that the |$\sim \!1400\,$| Å flux deficit of SDSS J025214.67|$-$|002813.7 is better interpreted by a reddened radiation from a single accretion disk [|$A_V= 0.17$| and |$R(V)=2.54$|] rather than a CBD model. The estimated dust extinction of WISE J0909|$+$|0002 is |$E(B-V)=0.13$|⁠, which can be translated into |$A_V = 0.4$| by assuming |$R_V = 3.1$| (Toba et al. 2021). This value is consistent with that of SDSS J025214.67|$-$|002813.7. However, |$A_V= 0.4$| is not particularly large and is below the value required for the dust extinction of |$A_V \sim \!7\:$|mag (Veilleux et al. 2013). Furthermore, the rest-frame wavelength range from the NUV to FUV bands overlaps with the Lyman limit system (Moravec et al. 2017). Hence, the observed deficit of NUV and FUV flux densities is likely due to intrinsic and/or IGM absorptions by neutral hydrogen.

4.3 Precession models

While the Doppler boost scenario (figures 9 and 10) was favored from our results, the CBD model was unlikely (figure 11a), since its model curves was not suitable for the SED of WISE J0909|$+$|0002. It is necessary to verify other plausible scenarios for supporting QPOs such as precessions of accretion disks and/or radio jets (e.g., McKinney et al. 2012; King et al. 2013; Tremaine & Davis 2014). For the accretion-disk precession scenario, Zhang (2022) revealed that the sizes of optical (⁠|$37 R_{\rm G}$|⁠) and NUV emission regions (⁠|$35R_{\rm G}$|⁠) in SDSS J132144|$+$|033055 are very similar to each other: the disk precession scenario is unlikely to be the cause of the periodic flux variability. The radio jet precession hypothesis was also ruled out by Chen et al. (2021), since deep |$6\:$|GHz radio imaging with NSF’s Karl G. Jansky Very Large Array (VLA) of three periodic quasars found that their radio flux was too weak (i.e., radio-quiet quasars) to affect UV and optical regions.

The radio flux of WISE J0909|$+$|0002 was not detected by VLA FIRST (Helfand et al. 2015).⁵ In order to estimate the upper limit of the radio loudness (⁠|$R \equiv f_{\rm 6}/f_{\rm 2500}$|⁠, where |$f_{\rm 6}$| and |$f_{\rm 2500}$| are the rest-frame flux densities at |$6\:$|cm and |$2500\,$|Å, respectively) with a |$1\:$|mJy detection threshold, we converted the radio flux limit at |$1.4\:$|GHz to that at |$5\:$|GHz, using a power-law in a radio-wavelength region, |$F_\nu \propto \nu ^{-0.7}$| (e.g., Condon 1992; Toba et al. 2019). Here, the SDSS spectrum taken in MJD 55532 was adapted for calculating |$f_{\rm 2500}$|⁠. Consequently, the radio loudness is estimated to have an upper limit of |$\le 0.4$|⁠, which is unlikely to affect the UV/optical flux (i.e., radio-quiet quasar). Thus the radio jet precession scenario is rejected for WISE J0909|$+$|0002. In Benke et al. (2023), the jet precession (or clear evidence of BSBH nature) was disfavored (or not found) from the periodic quasar PSO J334.2028|$+$|1.4075. They suggested that a misalignment between the inner jet and the outer lobes possibly originated from a warped accretion disk. Further observations to determine if broad line regions oscillate and/or accretion disks precess would bring new insights into QPO mechanisms or the evolution process of ELIRGs.

5 Conclusion

WISE J0909|$+$|0002 was identified as an ELIRG containing a type-1 quasar (Toba et al. 2021). We investigate the periodic flux variability in the UV range (quasar-rest frame) using archival and observational data obtained from the 105|$\:$|cm Murikabushi telescope, the Okayama and Akeno 50|$\:$|cm telescope/MITSuME and the SaCRA 55|$\:$|cm telescope/MuSaSHI with the OISTER collaboration. Our conclusions are summarized as follows:

1. The CRTS V-, g-, and r-band light curves (figure 2) and the combined light curve (figure 3) show significant and similar periodicity in each band (⁠|$\sim \!660$|–|$689\:$|day) using the SNR, the Lomb–Scargle method (figure 4), and the ACF (figure 5);

2. The light-curve periodicity of WISE J0909|$+$|0002 plausibly continues at least |$\sim \!3.6\:$|yr in the rest frame (figure 3);

3. The |$\tau$|–|$\sigma$| distributions estimated with the JAVELIN code and the structure functions do not present the DRW trend (figure 6), in other words it is unlikely that the periodicity of WISE J0909|$+$|0002 is a false-positive case at least in these analyses.

4. In 10000 simulations, the DRW light curves show periodic-like curves close to the best-fitting sinusoidal model on rare occasions (figure 7), thus the possibility that the light curves of WISE J0909|$+$|0002 are reproduced by the pure red noise cannot completely be rejected.

5. The PSD of the combined light curve supports the SPL |$+$| periodic model (figure 8).

6. The relativistic boost hypothesis is positive from the aspects of the ratio of variability amplitude and power-law slopes (figure 9) and the existence of a significant parameter space between the inclination angles and black hole mass (figure 10).

7. The SED of our target is insufficient to interpret using the CBD model. We cannot find cut-offs and/or notches in the SED predicted by this model (figure 11a). The drop-off of the FUV flux is probably explained by the intrinsic and/or IGM absorption by neutral hydrogen (figure 11b);

8. The radio jet precession scenario is ruled out, since the significant radio flux density of WISE J0909|$+$|0002 is not detected with the VLA FIRST Survey.

Our results generally agree with the relativistic boost scenario. To verify other hypotheses of periodic variability such as the warped-disk precession (e.g., Tremaine & Davis 2014), follow-up (multi-wavelength) observations are needed. In addition, if more longer and more accurate light curves can be obtained, we will be able to strengthen the validity of the periodicity and the relativistic boost for WISE J0909|$+$|0002.

Acknowledgments

We thank the anonymous referee for their valuable suggestions. This research is (partially) supported by the Optical and Infrared Synergetic Telescopes for Education and Research (OISTER) program funded by the MEXT of Japan. This work was (partially) carried out by the joint research program of the Institute for Cosmic Ray Research (ICRR), The University of Tokyo. The MITSuME system was supported by a Grant-in-Aid for Scientific Research on Priority Areas (19047003). This work was partially supported by JSPS KAKENHI Grant Numbers; 23K22537 (YT), 21H01126 and 23K20865 (TM), and 20K14521 (KI). TH thank to the staffs in Ishigakijima Astronomical Observatory and Yaeyama Star Club for the support in preparing the environment to describe this paper.

The CSS survey is funded by the National Aeronautics and Space Administration under Grant No. NNG05GF22G issued through the Science Mission Directorate Near-Earth Objects Observations Program. The CRTS survey is supported by the U.S. National Science Foundation under grants AST-0909182 and AST-1313422. MITSuME system was supported by a Grant-in-Aid for Scientific Research on Priority Areas (19047003).

Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions.

SDSS-IV acknowledges support and resources from the Center for High Performance Computing at the University of Utah. The SDSS website is |$\langle$|www.sdss.org|$\rangle$|⁠.

SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics | Harvard & Smithsonian, the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU)/University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

This study has made use of the NASA/IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

The Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max-Planck Institute for Astronomy, Heidelberg and the Max-Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, the Queen’s University Belfast, the Harvard-Smithsonian Center for Astro- physics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation grant No. AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation.

Supported by the National Science Foundation under Grants Nos. AST-1440341 and AST-2034437 and a collaboration including current partners Caltech, IPAC, the Oskar Klein Center at Stockholm University, the University of Maryland, University of California, Berkeley, the University of Wisconsin at Milwaukee, University of Warwick, Ruhr University, Cornell University, Northwestern University and Drexel University. Operations are conducted by COO, IPAC, and UW.

Footnotes

1

|$\langle$|https://github.com/LSST-sersag/pyzdcf|$\rangle$|⁠.

2

|$\langle$|https://github.com/nye17/javelin|$\rangle$|⁠.

3

|$\langle$|https://www.astroml.org/index.html|$\rangle$|⁠.

4

The BIC is defined as |$-2\ln (\mathcal {L})+k\ln (N)$|⁠, where |$\mathcal {L}$|⁠, k, and N are the likelihood function, the number of free parameters, and the number of data, respectively. On the basis of the Gaussian error model, the first term is described as |$N\ln ({\rm RSS}/N)$|⁠, where |${\rm RSS}$| is the residual sum of squares.

5

|$\langle$|https://third.ucllnl.org/cgi-bin/firstcutout|$\rangle$|⁠.