SMBH binary candidate PKS J2134−0153: possible multi-band periodic variability and inter-band time lags

ABSTRACT

Studying the periodic flux-variation behaviour of blazars is vital for probing supermassive black hole binaries and the kinematics of relativistic jets. In this work, we report the detection of multi-band possible periodic variations of the blazar PKS J2134−0153, including infrared (⁠|$1.6(\pm 0.4)\times 10^3$| day) and optical (⁠|$1.8(\pm 1)\times 10^3$| day) wavelengths. The periods in the infrared and optical bands are statistically consistent with the period in the radio band (⁠|$P_{\mathrm{Radio}} = 1760\pm 33$| days, obtained from our previous work). Moreover, flux variations in different bands are correlated with evident inter-band time delays, and the time lags of infrared and optical emission with respect to radio emission are |$(3.3\pm 2.3)\times 10^{2}$| days and |$(3.0\pm 2.3)\times 10^{2}$| days, respectively. The cross-correlations indicate a common origin of radio, infrared, and optical emission. The relative positions between emission regions of infrared and optical emission and radio emission are estimated according to the time lags, i.e. |$0.37\pm 0.26$| and |$0.33\pm 0.26$| pc. The relative distances seem to be quantitatively consistent with the theoretical prediction.

1 INTRODUCTION

Active galactic nuclei (AGNs) emit strong non-stellar emission in the Universe. Blazars are an important subclass of radio-loud AGNs, the relativistic jets of which point towards our Earth. Blazars emit a broad range emission from radio to |$\gamma$|-ray (e.g. Urry & Padovani 1995; Ulrich, Maraschi & Urry 1997). The central engines of blazars are often difficult to resolve directly because of their small physical scales and the large distance from Earth (i.e. the high redshift); their corresponding angular scales are too small to resolve using current facilities. Multi-band light variations offer an alternative possibility to resolve the central engine in the time domain (e.g. Blandford 1996; Fan et al. 1998). The flux variations provide valuable insights into various aspects of blazar studies, including structures, physical parameters of supermassive black holes (SMBHs, with masses larger than |$10^6\,\mathrm{ M}_{\odot }$|⁠), and emission processes (e.g. Ulrich et al. 1997; Gupta 2018). For instance, the minimum variability time-scale (minutes to hours) can constrain the emission-region sizes of blazars effectively. Combined with theoretical models and measured time lags between two different light curves via the reverberation mapping (RM) technique (e.g. Blandford & McKee 1982), one can constrain the physical parameters of relativistic jets, such as the magnetic field (e.g. Tavecchio, Maraschi & Ghisellini 1998).

Supermassive black hole binaries (SMBHBs), which include two orbiting SMBHs, are fascinating sources. The merging of SMBHBs produces gravitational waves (GWs: e.g. Thorne & Braginskii 1976; Haehnelt 1994; Jaffe & Backer 2003; Dou et al. 2022) detectable by low-frequency GW experiments (e.g. the Pulsar Timing Array (PTA): Hobbs et al. 2010; the Laser Interferometer Space Antenna (LISA): Amaro-Seoane et al. 2023). Spatially correlated signals consistent with the ultra-low-frequency GW background are detected by pulsar timing arrays (e.g. Agazie et al. 2023; EPTA Collaboration et al. 2023; Reardon et al. 2023; Xu et al. 2023). Although searching for SMBHBs is attractive, the direct identification of close SMBHBs (< 1 pc) is highly challenging because of the angular resolution limit of current telescopes (e.g. Dou et al. 2022).

Searching for periodicity in multi-band light curves is an attractive method for identifying SMBHB candidates. Year-like periodic flux variations can probe for the presence of SMBHB systems (e.g. Begelman, Blandford & Rees 1980; Lehto & Valtonen 1996; King et al. 2013; Ackermann et al. 2015; Liu et al. 2015, 2019; Bon et al. 2016; Charisi et al. 2016; Severgnini et al. 2018; Kovačević et al. 2020; Zhang & Wang 2022; Kun et al. 2023). A notable example is OJ 287, which shows a quasi 12-year periodicity in its light curves (e.g. Valtonen et al. 2006, 2008; Gupta et al. 2022; Sinitsyna & Sinitsyna 2022; Valtonen et al. 2024). This periodicity is explained by a SMBHB system where the secondary SMBH crosses the accretion disc of the primary SMBH, leading to the observed quasi-periodic variations (e.g. Dey et al. 2018).

Quasi-periodic behaviour over year-like time-scales can be explained by models with jet precession (e.g. Romero et al. 2000; Stirling et al. 2003; Caproni, Abraham & Monteiro 2013) or helical jets (e.g. Conway & Murphy 1993; Tateyama et al. 1998; Mohan & Mangalam 2015). In the SMBHB scenario, jet precession is driven by gravitational torque from the secondary black hole in misaligned orbits and the intrinsic helical jet structure could be driven by a coiled magnetic field (e.g. Vlahakis & Königl 2004; Miller et al. 2006). As a result, the viewing angle between the observer and the jet changes periodically, which is responsible for the observed periodic variations (e.g. Villata & Raiteri 1999). It should be noted that jet precession may happen in a single SMBH system. According to the disc-driven precession model (e.g. Sarazin, Begelman & Hatchett 1980; Lu 1990; Lu & Zhou 2005), misalignment between the SMBH’s spin axis and the accretion disc and the Lense–Thirring effect make the inner disc and jet precess. The periods produced by this model are often much longer than 10 years (for more details, see Section 5.1). Hence, it is important to find AGNs that show short-period (less than several years) periodic flux variations.

A number of works have been devoted to identifying SMBHB candidates via periodic variations. Graham et al. (2015) reported a possible close supermassive black-hole binary in PG 1302−102 with optical periodicity (the period is |$P = 1884 \pm 88$| days, but also see Liu, Gezari & Miller 2018, questioning this source). Sesana et al. (2018) selected 111 SMBH binary candidates from approximately 250 000 quasar light curves observed by the Catalina Real-time Transient Survey (CRTS) due to their potentially periodic nature. Bhatta (2018) proposed that the blazar J1043+2408 is an SMBH binary candidate because of its periodic flux variation. Li et al. (2018) reported a |$\sim$| 6.1-year quasi-periodicity of the blazar S5 0716+714, which is also an SMBH binary candidate. Serafinelli et al. (2020) found that the Seyfert 1.5 galaxy Mrk 915 emerges as a candidate SMBHB from the Swift Burst Alert Telescope (BAT) 105-month light curves of 553 AGNs. In our previous work, we reported a |$\sim$|4.69-year period in the radio-band light curve of BL Lac PKS J2134−0153 (Ren et al. 2021, hereafter R21; see also O’Neill et al. 2022, hereafter ON22), and we speculate that this is a close SMBH binary candidate with a small orbital separation of |$\sim 0.006$| pc.

In this work, we report on a multi-wavelength (including infrared, optical, and |$\gamma$|-ray) study of PKS J2134−0153. We also find inter-band correlations with evident time delays. Hence, we also constrain the emission-region sizes of PKS J2134−0153. This work is formatted as follows. In Section 2, we describe the multi-band light curves of PKS J2134−0153. In Section 3, we introduce the data analysis methods. In Section 4, we present the results. In Section 5, we present a discussion. Conclusions are made in Section 6.

Throughout this work, a |$\Lambda$|CDM cosmological model with |$\Omega _{\mathrm{m}} = 0.27$|⁠, |$\Omega _{\mathrm{\Lambda }} = 0.73$|⁠, |$H_{\mathrm{0}} = 70\,\,\mathrm{km\,\,s^{-1}\,\,Mpc^{-1}}$| is adopted for calculating the luminosity distance.

2 LIGHT CURVES

BL Lac PKS J2134−0153 (⁠|$z = 1.285$|⁠: see e.g. Rector & Stocke 2001; Truebenbach & Darling 2017) has been observed in various bands, including radio, infrared, optical, and |$\gamma$|-rays, from 2008–2023. The multi-band data we used to investigate the nature of the source are detailed in the following subsections. For gamma-ray data, please refer to Appendix  A.

2.1 Radio observations

PKS J2134−0153 has been observed in the radio band. The radio observations are discussed in detail in our previous work (R21; see also ON22) and will not be re-analysed in this work. In R21, we perform Lomb–Scargle periodogram (LSP: Lomb 1976; Scargle 1982) and weighted wavelet Z-transform (WWZ: Foster 1996) techniques to analyse the 15-GHz radio light curve obtained by the Owens Valley Radio Observatory (OVRO) 40-m telescope. The cadence and duration of the analysed light curve are |$\sim$| 7 days and |$\sim$| 11.4 years in the observed frame, respectively. We find that the radio-flux variations are periodic with an observed-frame period of |$4.69 \pm 0.14$| years, and can be fitted well by a sinusoidal function (R21). ON22 also use LSP and WWZ to analyse the observed-frame 45.1-year radio light curve of the same target. Their data are collected by the Haystack Observatory (15.5 GHz, with a cadence of |$\sim 1$| month in the observed frame), OVRO (15 GHz), and the University of Michigan Radio Astronomy Observatory (14.5 GHz, with a cadence of |$\sim 3$| months in the observed frame). The period (⁠|$2.082\pm 0.003$| years in the rest frame or |$4.757\pm 0.007$| years in the observed frame) obtained by ON22 is close to our result. The radio-band light curve and the best-fitting sinusoidal function obtained by R21 are shown in the top panel of Fig. 1. This periodic feature suggests that PKS J2134−0153 may be a candidate SMBHB. The disc-driven precession model involving a single SMBH can hardly account for the observed period (see Section 5.1).

2.2 Infrared observations

PKS J2134−0153 has been observed by the Wide-field Infrared Survey Explorer (WISE: Wright et al. 2010) and the Near-Earth Object Wide-field Infrared Survey Explorer (NEOWISE: Mainzer et al. 2014) in the infrared bands. WISE was launched on 2009 December 14, and is a 40-cm diameter infrared telescope designed to perform a repetitive all-sky survey at four infrared wavelength bands including |$W1$|⁠, |$W2$|⁠, |$W3$|⁠, and |$W4$| (centred at 3.4, 4.6, 12, and 22 |$\mathrm{\mu{\rm m}}$| in the observed frame). The survey was shut down on 2010 September 29 because the cooling material was depleted. The WISE spacecraft was reoperated in the all-sky survey model in |$W1$| and |$W2$| bands on 2013 December 13, and this is known as the NEOWISE survey. Each sky region is observed by the WISE telescope every 6 months.

The infrared-band observations of PKS J2134−0153 were obtained from the NASA/IPAC Infrared Science Archive¹ with a maximum matching radius of 3 arcsec. For the |$W1$| and |$W2$| bands, 21 and 239 pointings were obtained by WISE (between 2010 May 15 and 2010 November 11) and NEOWISE (from 2014 May 16 to 2022 November 1), respectively. The time span of these observations is 4553 days (⁠|$\sim$| 12.5 yr) in the observed frame. The observations include 20 epochs, and each spans 1–7 days. The time interval between two adjacent epochs is roughly 180 days, except for the gap from 2010 September 29 to 2013 December 13. Following Sheng et al. (2017, 2020) and Li et al. (2023), the bad observations, which were affected by low image quality (‘qi_fact’ |$\lt 1$|⁠), scattered moonlight (‘moon_masked’ |$\gt 0$|⁠), and the satellite’s small separation from the South Atlantic Anomaly (‘saa_sep’ |$\lt 5$|⁠), are removed. Finally, we retain 162 data points (the blue and green dots in the second panel of Fig. 1) in 19 epochs. For each epoch, we rebin the WISE observations because we are interested in the long-term variations. The observations are not considered if the total number of the data points in an epoch is less than three. The binned |$W1$| and |$W2$| light curves are shown as black and purple stars in the second panel of Fig. 1. It can be seen that there are significant flux variations in the two bands, and the variations appear quasi-periodic.

In ON22, the infrared (W1 and W2 bands) light curves of this source were obtained from WISE and NEOWISE (up to MJD 59150, i.e. 730 days shorter than our compiled WISE light curves). They did not perform periodic analyses of the infrared light curves. They calculated the cross-correlation of the infrared and OVRO 15-GHz radio light curves and obtained a |$\sim$| 250–350 day time lag between the two light curves (see fig. 6 of ON22). We measure the infrared–radio time lag independently because our requested WISE light curves are longer than those of ON22. Our infrared–radio time lag is statistically consistent with that of ON22 (see Section 4.2).

2.3 Optical observations

The light curves in optical bands include the V, g, r, i, and c bands. The V-band observations were obtained by the CRTS² (Drake et al. 2009), from 2005 May 3–2013 October 27 (MJD 53494–56593), spanning 3099 days (⁠|$\sim$| 8.5 years) in the observed frame, and have a total of 415 data points (the green points in the bottom panel of Fig. 1). The g-, r-, and i-band observations were obtained from the Zwicky Transient Facility³ (ZTF: Bellm et al. 2019; Masci et al. 2019). The start and end times of the g-, r-, and i-band light curves (see Fig. B1) are from MJD 58263–60244 (1981 days), MJD 58256–60244 (1988 days), and MJD 58280–58771 (491 days), respectively. The corresponding total number of data points is 326, 416, and 42, respectively. The observations in the c band are carried out with the Asteroid Terrestrial-impact Last Alert System⁴ (ATLAS: Heinze et al. 2018; Tonry et al. 2018) from MJD 57246–60292, spanning 3046 days (⁠|$\sim 8.3$| years) in the observed frame, with a total of 611 data points. The original optical light curves are shown in Fig. B1.

We adopt the CRTS, ZTF, and ATLAS data to construct a 6798-day (observed frame) long synthetic V-band light curve as follows.

• The optical spectrum (plate-mjd-fiberid=4384-56105-0988, shown in Fig. B2) of PKS J2134−0153 was obtained from Sloan Digital Sky Survey (SDSS) DR16 (Ahumada et al. 2020). The filter response curves of the CRTS V band, ZTF g band, and ATLAS c band are obtained from the Filter Profile Service.⁵

• We convolve the SDSS spectrum with the aforementioned filter response curves to obtain the corresponding magnitudes for the CRTS V band (⁠|$m_{\mathrm{SDSS},\ \mathrm{CRTS}-V}$|⁠), ZTF g band (⁠|$m_{\mathrm{SDSS,\ ZTF}-g}$|⁠), and ATLAS c band (⁠|$m_{\mathrm{SDSS,\ ATLAS}-c}$|⁠).

• We calculate the magnitude differences, |$\Delta m_{\mathrm{SDSS,}\ g-V} =m_{\mathrm{SDSS,\ CRTS-}V} - m_{\mathrm{SDSS,\ ZTF-}g}=-0.21$| mag, and |$\Delta m_{\mathrm{SDSS,\ }c-V} =m_{\mathrm{SDSS,\ CRTS-}V}\ -\ m_{\mathrm{SDSS,\ ATLAS-}c}=-0.04$| mag.

• We add |$\Delta m_{\mathrm{SDSS},\ gV}$| and |$\Delta m_{\mathrm{SDSS},\ cV}$| to the original ZTF g and ATLAS c light curves, respectively.

• The offset-corrected ZTF g and ATLAS c light curves and the CRTS V-band light curve are appended to build a new synthetic V-band light curve.

The new synthetic V-band light curve has 1352 data points from 2005 May 3–2023 December 14 (MJD 53494–60292, see the bottom panel of Fig. 1). Similarly to the infrared band light curves, we rebin the new synthetic V-band light curve using a 180-day bin (the dark red stars in the bottom panel of Fig. 1). It is worth noting that ON22 analysed the CTRS data (approximately from MJD 53491–57357) and part of the ZTF data (approximately from MJD 58239–59458) previously and did not find any periodicity. Here, we build longer light curves with better cadences than ON22 from the CTRS, ZTF, and ATLAS datasets. For a detailed periodicity analysis, see Section 4.1.

3 DATA ANALYSIS

We search for periodic variation in the multi-band light curves of PKS J2134−0153 via the LSP (Lomb 1976; Scargle 1982). We use damped random walk (DRW) simulations to estimate the significance of the period. In addition, we use the cross-correlation function (CCF) to calculate the inter-band time lags.

3.1 Lomb–Scargle periodogram

The LSP is a commonly used method in time-series analysis, and is based on the discrete Fourier transform (DFT) method (e.g. Bhatta 2018; Prince et al. 2023). The LSP is suitable for detecting and characterizing periodic variations in unevenly sampled time series and avoiding the interpolation of data gaps (e.g. Li et al. 2017; Bhatta 2018; VanderPlas 2018). In this analysis, we run the period search on a frequency grid. To construct the frequency grid, one must specify the frequency limits (⁠|$f_{\mathrm{min}}$|⁠, |$f_{\mathrm{max}}$|⁠) and the grid spacing (determined by the total number of required periodogram evaluations, |$N_{\mathrm{eval}}$|⁠). According to VanderPlas (2018), the minimum frequency adopted is |$f_{\mathrm{min}} = 1/t_{\mathrm{span}}$|⁠, where |$t_{\mathrm{span}}$| is the time span of the light curve. The maximum frequency adopted is |$f_{\mathrm{max}} = 1/(2\left\langle \delta t \right\rangle)$|⁠, where |$\left\langle \delta t \right\rangle$| represents the mean sampling time (five days) of the synthetic V-band light curve. The total number of required periodogram evaluations is

and we adopt |$n_{\mathrm{0}}$|=10 (e.g. Urry & Padovani 1995; Debosscher et al. 2007; Richards et al. 2012; R21). Frequencies are linearly and uniformly distributed. The light curves are fitted by the following sinusoidal function:

where t is the sequence of observation times, |$f_{\mathrm{best}}$| is the best frequency (corresponding to the maximum power) fixed by the LSP method, and A, |$\varphi$|⁠, and C represent the amplitude, phase, and offset of the simulated light curve, respectively. The best-fitting model parameters are obtained via the standard minimization |$\chi ^{2}$| method (LMFIT: Newville et al. 2014). Then, for each band, the best-fitting sinusoidal function can be obtained from the best-fitting parameters (i.e. the dashed red curves in Fig. 1).

3.2 Significance estimation

The significance levels of LSP power should be estimated carefully. Red noise can easily mimic periodic signals for unevenly sampled light curves and cause large false peaks in the LSP power. Meanwhile, AGN light curves can be well described by the DRW (a popular red noise model: Kelly, Bechtold & Siemiginowska 2009; Kozłowski et al. 2010; MacLeod et al. 2012; Zu et al. 2013; Secunda, Jiang & Greene 2024), the power spectral density (PSD) of which is

The free parameters are the normalization (C) and the breaking frequency (⁠|$f_{\mathrm{break}}$|⁠). The PSD of the DRW process is dominated by white noise at the low-frequency end and a random walk at the high-frequency end. The DRW model has been applied to fit the light curves of blazars (e.g. Zhang et al. 2022). We use the celerite (Foreman-Mackey et al. 2017) code to fit the DRW model to each light curve and obtain the best-fitting parameters (i.e. |$f_{\mathrm{break}}$| and C) by maximizing the corresponding celerite likelihood (via scipy: Virtanen et al. 2020). Second, we use celerite to generate 20 000 evenly sampled light curves according to the best-fitting DRW parameters, and the time duration of each one is the same as the observed light curve. For each of the 20 000 generated light curves, we then use the linear interpolation method to construct a light curve that shares the same sampling as the observed one. That is, like real observations, each simulated light curve is unevenly sampled. Third, the simulated 20 000 unevenly sampled light curves are analysed by the LSP method and 20 000 simulated LSP power spectra are obtained. Again, just like the PSDs of real light curves, the PSD of each simulated light curve is affected by the same window function and uneven sampling. Finally, we calculate the 84.13th, 97.72th, and 99.87th percentiles of the 20 000 simulated LSP power for each frequency, which correspond to the |$1\sigma$|⁠, |$2\sigma$|⁠, and |$3\sigma$| significance levels. The periodic analysis of the radio observations, which suggests that the period is |$1760\pm 33$| days (R21; ON22), provides important prior knowledge of the expected period in the optical and IR bands. Hence, we only focus on a period range of 1661–1859 days (i.e. the |$3\sigma$| lower and upper limits of the radio period) and do not need to consider the ‘Look Elsewhere’ effect (e.g. Gross & Vitells 2010). The corresponding period is insignificant if the observed LSP power peak is lower than the simulated LSP power spectra. Our code to perform the significance analysis can be downloaded freely from the supplementary material for this article.

3.3 Cross-correlation function

We calculate the possible time lags between the multi-band light curves of PKS J2134−0153 via the public python version of the interpolated cross-correlation function (ICCF: Peterson et al. 1998, 2004), pyccf (Sun, Grier & Peterson 2018), which is one of the commonly used methods in time-series analysis of AGNs. The code uses the linear interpolation method to deal with unevenly sampled light curves and calculate the cross-correlation coefficient as a function of the time lag for two light curves. The ICCF is evaluated for a time lag (⁠|$\tau$|⁠) range from |$-4000$| days to 4000 days. The searching step |$\Delta \tau$| should be smaller than the median sampling time |$\left\langle \delta t \right\rangle$| of the light curves. We adopt |$\Delta \tau =60$| days in this work. As will be demonstrated in Section 4.2, the ICCF has two strong peaks. For each peak, we calculate its corresponding centroid using only the ICCF for the time lags around the peak. In this work, we adopt the centroid of the cross-correlation function (⁠|$\tau _{\mathrm{cent}}$|⁠) using only time lags with |$r \gt 0.8 r_{\mathrm{max}}$|⁠, where |$r_{\mathrm{max}}$| represents the peak value of the CCF. The 1|$\sigma$| confidence ranges can be calculated via the popular random subset selection and flux redistribution method. pyccf adopts the Student’s t statistic to calculate the p-value of each peak; the p-value is the probability of observing similar or stronger cross-correlations under the no-correlation hypothesis. If the p-value is larger than 0.01, the peak is statistically insignificant.

4 RESULTS

4.1 Lomb–Scargle results

We perform LSP analysis of the light curves in the infrared and V bands of PKS J2134−0153. We compute LSPs from minimum (⁠|$f_{\mathrm{min}}$|⁠) to maximum frequencies (⁠|$f_{\mathrm{max}}$|⁠). The values of |$f_{\mathrm{min}}$| and |$f_{\mathrm{max}}$| for each light curve are obtained according to the description in Section 3.1. The total number of periodogram frequencies |$N_{\mathrm{eval}}$| is set according to equation (1), while |$f_{\mathrm{min}}$|⁠, |$f_{\mathrm{max}}$|⁠, and |$N_{\mathrm{eval}}$| for each light curve are shown in Table 1.

The LSP power spectra in the synthetic unbinned V band are shown in the left panel of Fig. 2. Some ‘fake’ PSD peaks are likely caused by the irregular sampling (the pink curves in Fig. 2). We test the significance of peaks using the method described in Section 3.2. There are several peaks that cannot be attributed entirely to irregular sampling effects. All of these peaks have significance levels |$\lesssim 2.5\sigma$|⁠. The significance level is based on the fact that the possibility of the DRW process reproducing the observed peak, or a more statistically significant peak, is |$1{{\ \rm per\ cent}}$|⁠. Peaks with periods of |$\sim 100$| days are not observed in the radio data and are neglected in subsequent analysis. The highest peak corresponds to the period |$P_{\mathrm{V}} = 1.8(\pm 1)\times 10^3$| days and seems unable to explain by irregular sampling. This period is very close to the period detected in the radio light curve (Section 2.1). Hence, we will focus on this tentative detection in the V band.

It is worth noting that ON22 have performed a similar periodic analysis. Unlike our work, ON22 only considered CRTS and ZTF observations, with 322 data points and a baseline of |$\sim$| 6000 days. Hence, they did not find any significant peaks in their PSD. We think that ATLAS data can fill the gap between CRTS and ZTF observations, hence improving the PSD estimation. Indeed, if we neglect the ATLAS observations and the new ZTF data released in 2022 and 2023 in our synthetic V-band light curve, we also cannot detect any statistically significant periods. Our synthetic 180-day bin V-band light curve is fitted with a sinusoidal function (i.e. equation 2, where |$f_{\mathrm{best}}$| is fixed to |$1/P_{\mathrm{V}}$|⁠) and the best-fitting parameters are |$A_{\mathrm{V}} = 0.36$|⁠, |$\varphi _{\mathrm{V}} = 1.3$|⁠, and |$C_{\mathrm{V}} = 17.97$|⁠, respectively. The best-fitting curve is shown in the bottom panel of Fig. 1. The |$\chi ^{2}$| of the sinusoidal fit is 51. We also fit the data with a constant function and the |$\chi ^{2}$| is 123. In addition, we calculated the Akaike information criterion (AIC) for the sinusoidal and constant fits to be 19 and 46, respectively. Statistically speaking, if the AIC difference between two models is larger than 10, the model with the larger AIC value has relatively little support (e.g. Burnham, Anderson & Huyvaert 2011). Hence, the sinusoidal function is preferred over the constant fit. It should be noted that the AIC difference does not appear to be large enough for the sinusoidal fit to dominate completely.

The IR variations appear to be periodic. The LSP power spectrum for the unbinned |$W1$| light curve is shown in the right panel of Fig. 2. The LSP power spectrum suggests that there is also a |$P_{\mathrm{IR}} = 1.6(\pm 0.4)\times 10^3$| day period in the IR light curve. However, the confidence of this peak is about 2|$\sigma$| due to the small number of IR data points. We performed the same sinusoidal fit (⁠|$f_{\mathrm{best}}$| is fixed to |$1/P_{\mathrm{IR}}$|⁠) to the binned infrared (⁠|$W1$|⁠) light curve (the middle panel of Fig. 1). The best-fitting parameters are |$A_{\mathrm{IR}} = -0.53$|⁠, |$\varphi _{\mathrm{IR}} = 1.27$|⁠, and |$C_{\mathrm{IR}} = 12.43$|⁠. The |$\chi ^{2}$| of the sinusoidal fit is 2400 and the |$\chi ^{2}$| of the constant function fit is 5067 (the same fitting analysis was performed for the |$W2$| band; see the middle panel of Fig. 1). Similarly, we calculate the AIC values of 97 and 108 for the sinusoidal and constant fits of the W1-band light curve, respectively. Thus, the light curve in the IR band can be fitted better by the sinusoidal function (e.g. Burnham et al. 2011). The value of |$1.6(\pm 0.4)\times 10^3$| days is statistically consistent with the period in the radio (⁠|$P_{\mathrm{Radio}} = 1760\pm 33$| days: R21; ON22) and optical bands. Hence, it is very likely that the same fluctuating synchrotron radiation produces radio, infrared, and optical variations.

We phase-fold the V-band and IR light curves according to their LSP spectral peaks. First, the phase is calculated as |$((\mathrm{MJD-MJD_{ref}}){{\ \rm per\ cent}}P)/P$|⁠, where |${{\ \rm per\ cent}}$| and P represent the modulo operator and period, respectively. Second, the light curves are rearranged according to the phases. Third, rearranged light curves are rebinned with a phase bin size of 0.05 for the optical data. For the WISE IR bands, we rebin the phase light curves by the WISE observed epoch. The results are shown in Fig. 3. It appears that the light curves show sinusoidal variations on top of some short-term random variations. In addition, the phase-folded V- and IR-band light curves are fitted by the sinusoidal and constant functions. We calculate ΔAIC values of 187 (for the V band) and 60 (for the IR band). Hence, the constant model is less supported by the data than the sinusoidal model (e.g. Burnham et al. 2011). Compared with the radio variations, sinusoidal variations are less prominent in the phase-folded optical light curve, possibly because the accretion disc also produces a significant fraction of optical emission and this disc component varies statistically, as commonly seen in radio-quiet AGNs (Ulrich et al. 1997).

Overall, possible periodic signals are present in the PKS J2134−0153 optical and infrared light curves with significance levels of |$\sim 2\!-\!2.5\sigma$|⁠. Future observations with longer durations can test and verify our results.

In addition, we obtain the |$\gamma$|-ray flux variations of PKS J2134−0153 from the Pass 8 database of the Fermi data server⁶ (see Appendix  A). We perform LSP analysis of the |$\gamma$|-ray light curve and cannot detect any statistically significant periods.

4.2 Time lag

The time lags among the multi-bands of PKS J2134−0153 are calculated using the ICCF (see e.g. Peterson et al. 1998) method described in Section 3.3. The parameters adopted for pyccf (Sun et al. 2018) are listed in Table 2.

ON22 calculate the time lag between the infrared and radio light curves through the cross-correlation function. They find that the infrared variation leads the radio variation by |$\sim$||$250\!-\!350$| days, and the confidence level of this value is more than |$2\sigma$|⁠. In this work, we re-estimate the IR–radio time lag because our IR light curves are longer than that of ON22 by about 600 days.

Due to the observational strategy of WISE, the gap between two adjacent epochs is about 180 days (see Section 2.2 and Fig. 1). Similarly, the 180-day binned light curves in the V band (the dark red stars in the bottom panel of Fig. 1) and radio band are adopted. The ICCF method requires the interpolation of unevenly sampled light curves, and the interpolation for binned data is more statistically robust than for unbinned data. We confirm that the time-lag measurements remain unchanged if the unbinned light curves are adopted.

There are cross-correlations between the radio and optical/WISE light curves. The time lags between the 180-day binned synthetic V-band and binned radio light curves are calculated from |$-4000$| days to 4000 days with a step of |$\Delta \tau$| = 60 days. High cross-correlation coefficients might be unreliable because of statistical fluctuations. Following Grier et al. (2018), we apply a prior weight to the cross-correlation coefficients, defined by |$\omega (\mathrm{\tau })=\left[ N(\mathrm{\tau })/N(\mathrm{0}) \right] ^{2}$|⁠, where |$N(\mathrm{\tau })$| represents the number of overlapping points at |$\tau$| (see column 7 of Table 2) and |$N(\mathrm{0})$| is the number of overlapping points at |$\tau = 0$|⁠. Then, the weighted cross-correlation coefficients (⁠|${\rho }^{\prime }$|⁠) can be calculated, i.e. |${\rho }^{\prime } = \rho \omega$| (see column 8 of Table 2). The weighted ICCF is shown in the top panel of Fig. 4. There are two distinct peaks. For each one, we use the ICCF at lag ranges of [0, 1000] or [|$-2000$|⁠, |$-500$|] to calculate the time lags and their uncertainties; the corresponding time lags are listed in Table 2, i.e. |$(3.0\pm 2.3)\times 10^2$| and |$(-1.3\pm 0.4)\times 10^3$| days. The cross-correlation between the binned optical–radio light curves is statistically significant because the p-values for the two ICCF peaks are |$6\times 10^{-6}$| and 0.03, respectively. To assess the significance, we also calculate the weighted ICCF between the observed radio and DRW simulated V-band light curves (see Section 3.2); the results are shown in the bottom panel of Fig. 4. We simulate 20 000 optical band light curves by the DRW model, and calculate the mock ICCFs between these simulated binned optical light curves and the observed radio light curve, respectively. The grey dotted line represents the |$99{{\ \rm per\ cent}}$| upper limit of the 20 000 mock ICCFs. We find that the possibility of obtaining the observed ICCF peak or an even higher peak is again |$1{{\ \rm per\ cent}}$|⁠. We also calculate the p-value (hereafter |$p_{\mathrm{mock}}$|⁠) for the highest peak in each mock ICCF. The distribution of the 20 000 cases of |$p_{\mathrm{mock}}$| is shown as a blue histogram in Fig. 5. The solid purple curve represents the cumulative distribution of |$p_{\mathrm{mock}}$|⁠. The red vertical line corresponds to the median of these |$p_{\mathrm{mock}}$|⁠. The black vertical line corresponds to the p-value (hereafter |$p_{\mathrm{opt\ vs\ Radio}}$|⁠) of the CCF peak between the real optical and radio-band light curves (the blue curve in the top panel of Fig. 4). The probability of having |$p_{\mathrm{mock}}\lt p_{\mathrm{opt\ vs\ Radio}}$| is only |$0.38{{\ \rm per\ cent}}$|⁠. Hence, we argue that the optical–radio cross-correlation is likely to be statistically significant.

We also obtain significant cross-correlations between the radio and WISE light curves. The time lag, maximum correlation coefficient, and corresponding p-value for the WISE|$W1$| band are listed in Table 2. The positive peak we obtained is statistically consistent with that of ON22. Just like the radio–optical ICCF, the radio–IR ICCF also shows two significant peaks and the time-lag difference for the two peaks is |$\simeq 1700$| days, which is close to the variational period (see Section 4.1). Moreover, the radio–optical and radio-IR time lags are statistically similar. Indeed, the weighted ICCF between the binned optical–IR light curves has a high weighted cross-correlation coefficient (⁠|$\rho ^{\prime } = 0.81$|⁠), and the corresponding time lag is |$0\pm 140$| days (Fig. 6 and Table 2). The cross-correlation between the binned optical–IR light curves is statistically significant because the p-value for this peak is |$1\times 10^{-5}$|⁠. The significant correlation with the zero-day lag suggests a common origin of the IR and optical emission.

5 DISCUSSION

5.1 SMBH binary system

We find possible periodic variations in the infrared and optical bands; their observed-frame periods are |$P_{\mathrm{IR}} = 1.6(\pm 0.4)\times 10^3$| and |$P_{\mathrm{V}} = 1.8(\pm 1)\times 10^3$| days, respectively. The two periods are statistically consistent with the period in the radio band (⁠|$P_{\mathrm{Radio}} = 1760\pm 33$| days: R21; ON22). Moreover, our ICCF analyses suggest that the cross-correlation between the radio and infrared/optical variations is significant. The optical period and the cross-correlation between the optical and IR emission are not reported by R21 or ON22.

The infrared emission has two possible origins: (1) the electron synchrotron radiation in the jet (e.g. Ghisellini & Tavecchio 2009; Böttcher et al. 2013) and (2) the reprocessed emission of the dusty torus. As proposed by Massaro et al. (2011), the origin of the IR emission can be understood by the |$W2-W3$| (3.22 mag) and |$W1-W2$| (1.08 mag) colours, where the mean magnitudes of the |$W1$|⁠, |$W2$|⁠, and |$W3$| bands are 12.37, 11.29, and 8.07 mag, respectively. Following Massaro et al. (2011, their fig. 3), the |$W2-W3$| versus |$W1-W2$| of our target is also consistent with a non-thermal jet origin.

We can estimate the rest-frame infrared and optical luminosities of this source by

where |$D_{\mathrm{L}}$| is the luminosity distance, |$\nu$| is the frequency, and |$f_{\nu }$| is the flux density. For the infrared bands, |$f_{\nu } = f_{\nu 0} 10^{(-m/2.5)}$| and the zero-point flux |$f_{\nu 0}$| for |$W1$| and |$W2$| (Jarrett et al. 2011) are 309.54  and 171.787 Jy, respectively. The average magnitude of the |$W1$| band is |$m_{\mathrm{W1}} = 12.37\ \mathrm{mag}$|⁠. The rest-frame wavelength for |$W1$| is |$3.4\, \mu{\rm m}/(1+z) \simeq 1.5\, \mu{\rm m}$| and the corresponding infrared luminosity is |$L_{\mathrm{1.5\, \mu{\rm m}}} = 3.1 \times 10^{46}\ \mathrm{erg\ s^{-1}}$|⁠. For the optical band, the source flux density can calculated by |$m_{\nu } = -2.5\log _{10}(f_{\nu })-48.6$|⁠; the average magnitudes of the g and r bands are |$m_{g} = 18.07\ \mathrm{mag}$|⁠, |$m_{r} = 17.53\ \mathrm{mag}$|⁠, respectively; the rest-frame wavelengths for g and r are |$4770\, \mathrm{\mathring{A}}/(1+z)=2087\, \mathrm{\mathring{A}}$| and |$6231\, \mathrm{\mathring{A}}/(1+z)=2727\, \mathrm{\mathring{A}}$|⁠, respectively; the corresponding rest-frame luminosities are |$L_{\mathrm{2087\, \mathring{A}}} = 1.4 \times 10^{46}\ \mathrm{erg\ s^{-1}}$| and |$L_{\mathrm{2727\, \mathring{A}}} = 1.7 \times 10^{46}\ \mathrm{erg\ s^{-1}}$|⁠, respectively. The luminosity at |$3000\ \mathrm{\mathring{A}}$| rest-frame (⁠|$L_{3000\, \mathrm{\mathring{A}}}= 1.8\times 10^{46}\ \mathrm{erg\ s^{-1}}$|⁠) is calculated by linearly extrapolating |$L_{\mathrm{2087\, \mathring{A}}}$| and |$L_{\mathrm{2727\, \mathring{A}}}$|⁠. The bolometric luminosity of an AGN disc can be estimated from |$L_{\mathrm{bol}} = 5L_{3000\, \mathrm{\mathring{A}}} = 9.0\times 10^{46}\ \mathrm{erg\ s^{-1}}$| (e.g. Richards et al. 2006), which is slightly larger than |$L_{\mathrm{1.5\, \mu{\rm m}}}$|⁠. It is worth noting that the optical luminosity may be underestimated if dust extinction is important.

Given that the infrared and optical periods are statistically consistent with the radio period, we propose that the synchrotron radiation of the same electron population produces periodic features in the optical, infrared, and radio bands. The variation period has several possible origins: for instance, the disc-driven precession model (e.g. Sarazin et al. 1980; Lu 1990) due to the Lense–Thirring effect and misalignment between the SMBH’s spin axis and the accretion disc. As a result, the jet precesses around the central supermassive black hole. The precession period can be estimated by (Lu & Zhou 2005)

where |$M_{\mathrm{abs}}$| is the absolute magnitude of the source in the B band, |$\alpha$| is the dimensionless viscosity parameter, a is the dimensionless specific angular momentum, M is the black hole mass, and |$\eta \sim 0.1$| is the radiative efficiency. |$M_{\mathrm{abs}}$| can be calculated using |$M_{\mathrm{abs}}=M_{\mathrm{B}}-2.5\log \left({L_{\mathrm{B}}/{\rm L}_{\odot } } \right)$|⁠, where |$L_{\mathrm{B}}$| is the luminosity of the B band (the rest-frame wavelength is |$\sim 4400\, \mathrm{\mathring{A}}$|⁠), |$M_{\mathrm{B}} = 5.44$| is the B-band absolute magnitude of the Sun (Willmer 2018), and |${\rm L}_{\odot }$| is the B-band luminosity of the Sun. If we adopt |$L_{\mathrm{B}} = 2.5\times 10^{46}\ \mathrm{erg\ s^{-1}}$|⁠, which is estimated by linearly extrapolating |$L_{\mathrm{2087\, \mathring{A}}}$| and |$L_{\mathrm{2727\, \mathring{A}}}$|⁠, and we can calculate |$M_{\mathrm{abs}}= -27.31$| mag. For a typical mass of |$M_{\mathrm{BH}}\simeq 10^{8}\ {\rm M}_{\odot }$|⁠, |$\alpha =0.4$| (King, Pringle & Livio 2007), and |$a=0.9$| (i.e. a highly spinning SMBH), the corresponding precession period is |$10^{5.2}$| years. Hence, it is hard to reconcile this precession period with the observed ones in the optical and IR, unless the spin is close to zero, |$L_{\mathrm{B}}\simeq 1.1\times 10^4L_{\mathrm{1.5\, \mu{\rm m}}}$|⁠, or the dimensionless viscosity parameter |$\alpha$| is⁷|$\simeq 10^{-4}$|⁠.

The second class of models involves an SMBH binary. In such a system, precession of the jet can be induced because of the misalignment of the orbital plane of the secondary SMBH with the primary SMBH’s accretion disc (e.g. Abraham 2000; Romero et al. 2000; Caproni & Abraham 2004; De Paolis, Ingrosso & Nucita 2004; Liu & Chen 2007; Caproni et al. 2017; Serafinelli et al. 2020; Casey-Clyde et al. 2022; Witt et al. 2022) or accretion disc warping (e.g. Britzen et al. 2018). Hence, the angle between the jet and our sightline changes periodically on the orbital time-scale. The same mechanism is one of the physical explanations for the geometric helical jet that explains the periodic variations in a number of blazars (e.g. Conway & Murphy 1993; Bahcall et al. 1995; Tateyama et al. 1998; Zhou et al. 2018; Jorstad et al. 2022; Zhang & Wang 2022). Based on the long-term (⁠|$\sim 45$| yr in observed frame) radio variations, ON22 propose that PKS J2134−0153 hosts a binary SMBH and the orbital motion of the primary SMBH (along with the jet) should be responsible for the radio light curve, with the separation of the two SMBHs only about |$10^{16}$| cm (also see R21). In this work, we find that the infrared and optical periods are statistically consistent with the radio period. In addition, there is a significant correlation between the V and infrared emission with significant time lags. These new results suggest that the driving mechanism for periodicity for the optical and IR emission regions is similar to that in the radio emission region, albeit there is about |$\sim 1$| pc separation between the optical/IR and radio emission regions. That is, the jet precession acts as a rigid body. In summary, our periodicity and cross-correlation analyses reveal more details about jet precession in the SMBH binary candidate PKS J2134−0153 (e.g. R21; ON22).

5.2 The position of emission regions

As shown in Section 4.1, the periodic variations in the infrared and optical bands (⁠|$P_{\mathrm{IR}} = 1.6(\pm 0.4)\times 10^3$| days, |$P_{\mathrm{V}} = 1.8(\pm 1)\times 10^3$| days) are statistically consistent with the radio period. The significance level of |$P_{\mathrm{V}}$| is more than 2.5|$\sigma$|⁠. ON22 reported a strong correlation (with correlation coefficient |$\sim$|1) between the infrared- and radio-band light curves. We find significant cross-correlations between infrared, optical, and radio variations, since the p-values are much smaller than 0.01 (see Section 4.2). In summary, the infrared, radio, and a significant fraction of the optical variable fluxes are likely to be produced by the same precession or helical jet.

The time lags in optical, infrared, and radio variations suggest that the emission regions of optical/infrared are different from radio emission regions. The optical and infrared emission regions may be overlapped significantly. Max-Moerbeck et al. (2014) proposed a model to interpret the time lags. In their model, a moving disturbance moves across the jet core and induces multi-band variations in the corresponding emission regions; the infrared emission regions are closer to the centre SMBH than the radio emission regions, and the V-band emission regions should be closer to the centre SMBH than the infrared band. Hence, negative time lags should be excluded. In the subsequent analysis, we adopt |$(3.3\pm 2.3) \times 10^2$| days as the time lag between infrared and radio emission and |$(3.0\pm 2.3) \times 10^2$| days as the time lag between optical and radio emission. The relative distance between different emission regions can be calculated by the time lags between the corresponding light curves, i.e.

where |$\Gamma$| is the bulk Lorentz factor, |$\delta$| is the Doppler factor, |$\beta c$| is the bulk jet speed, c represents the speed of light, |$T_{\mathrm{lag}}$| is the time lag, and z is the redshift (e.g. Pushkarev, Kovalev & Lister 2010; Max-Moerbeck et al. 2014). Following Pushkarev et al. (2010), we can use the |$\delta$|–|$\beta _{\mathrm{app}}$| relation of Cohen et al. (2007) and equation (6) to obtain

where |$\beta _{\mathrm{app}}$| and |$\theta$| are the apparent bulk jet speed and the viewing angle, respectively. We obtained the apparent bulk jet speed |$\beta _{\mathrm{app}}=0.19\pm 0.12$| from the Monitoring Of Jets in Active galactic nuclei with Very Long Baseline Interferometry (MOJAVE) Very Long Baseline Interferometry (VLBI) experiments (e.g. Lister & Homan 2005; Lister et al. 2019). The viewing angle is assumed to be |$3{_{.}^{\circ}}6$|⁠, a typical value for blazars (Hovatta et al. 2009; Pushkarev et al. 2009). The time lags between different bands are listed in Table 2. Hence, the relative distances between infrared and radio regions and between optical and radio regions are |$0.37\pm 0.26$| and |$0.33\pm 0.26$| pc, respectively.

The radio core size can be determined by the angular diameter, |$\theta _{\mathrm{core}}$|⁠,

where |$d_{\mathrm{A}}=1742.7\,\,\mathrm{Mpc}$| is the angular diameter distance for a |$\Lambda$| cold dark matter cosmology (e.g. Komatsu et al. 2011) and redshift |$z=1.285$|⁠. |$\theta _{\mathrm{core}}$| represents the angular diameter of the radio core, |$\alpha _{\mathrm{int}}$| is the intrinsic opening angle. For PKS J2134−0153, |$\theta _{\mathrm{core}}$| is measured by MOJAVE (Lister & Homan 2005; Lister et al. 2019), yielding the value |$\theta _{\mathrm{core}} = 0.3\,\,\mathrm{mas}$| (Kovalev et al. 2005). We adopt |$\alpha _{\mathrm{int}} = 2{_{.}^{\circ}}4\pm 2{_{.}^{\circ}}0$|⁠, the mean value for BL Lacs (Pushkarev et al. 2009). According to equation (8), we find |$d_{\mathrm{core}}=60\pm 11\,\,\mathrm{pc}$|⁠. The relative distances between infrared and radio and between optical and radio are much smaller than the core size.

These size measurements can be used to test jet models. For the synchrotron self-absorption model (for details, see Lobanov 1998; Kovalev, Lobanov & Pushkarev 2008a; O’Sullivan & Gabuzda 2009; Tramacere et al. 2022), the positions of the different emission regions scale as

where r and |$\nu$| are the position of the emission region and the frequency, respectively. According to equation (9), the distances between the different emission regions can be calculated as

where |$r_{\mathrm{Radio} }$| is the absolute position of the radio-band emission region. Hence, their ratio is |$f_\mathrm{th} = \Delta r_{ {V-\mathrm{Radio}}}^{\mathrm{th}} / \Delta r_{ \mathrm{IR-Radio}}^{\mathrm{th}}=(1-\nu _{\mathrm{Radio}}/\nu _{\mathrm{V}}) / (1 - \nu _{\mathrm{Radio}}/\nu _{\mathrm{IR}}) \simeq 1$|⁠. This ratio is consistent with the observed one (⁠|$0.9\pm 0.7$|⁠), and also consistent with the zero-day time lag between the optical and infrared band light curves.

6 CONCLUSIONS

In this work, we have searched systematically for possible periodic variations in the optical and infrared light curves of PKS J2134−0153. Our results are as follows.

(1) Using the Lomb–Scargle periodogram, we have found a periodic signal with a possible period of |$1.6(\pm 0.4)\times 10^3$| days in the infrared band light curve and a |$1.8(\pm 1)\times 10^3$| day possible period in the V band of PKS J2134−0153 (see Section 4.1 or Fig. 2). The two periods (both in the observed frame) are consistent with the period in the radio band (⁠|$P_{\mathrm{Radio}} = 1760 \pm 33$| days: R21; ON22). A similar period was recently reported by Kiehlmann et al. (2024). The SMBH binary system likely generates periodic variations in these three bands.

(2) The distances between various emission regions are calculated according to the inter-band time lags. The relative distance between the infrared and radio emission regions is |$0.37\pm 0.26$| pc. The relative distance between the V and radio emission regions is |$0.33\pm 0.26$| pc (see Section 5.2 and Table 2). The time lag between the infrared and V-band light curves is |$0\pm 140$| days (see Section 4.2, Table 2, and Fig. 6). This suggests that the emission regions of the infrared and optical bands may be very close to each other. The relative distances are quantitatively consistent with the theoretical prediction (Lobanov 1998; Kovalev et al. 2008a ; O’Sullivan & Gabuzda 2009; Tramacere et al. 2022).

The long multi-band monitoring of PKS J2134−0153 is needed to confirm our conclusions in future work.

ACKNOWLEDGEMENTS

GWR and MYS acknowledge support from the National Natural Science Foundation of China (NSFC-11973002; NSFC-12322303), the Natural Science Foundation of Fujian Province of China (No. 2022J06002), and China Manned Space Project grants (No. CMS-CSST-2021-A06; No. CMS-CSST-2021-B11). ND is sincerely grateful for the financial support of the Xingdian Talents Support Program, Yunnan Province (NO. XDYC-QNRC-2022-0613). ZXZ acknowledges support from the National Science Foundation of China (NSFC-12033006; NSFC-12103041). We thank the referee for constructive comments that greatly improved this work. We thank Haikun Li for maintaining the computing resources.

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.

Based on observations obtained with the Samuel Oschin Telescope 48-inch and the 60-inch Telescope at the Palomar Observatory as part of the Zwicky Transient Facility (ZTF) project. ZTF is supported by the National Science Foundation under Grants No. AST-1440341 and AST-2034437 and a collaboration including current partners Caltech, IPAC, the Weizmann Institute for Science, the Oskar Klein Center at Stockholm University, the University of Maryland, Deutsches Elektronen-Synchrotron and Humboldt University, the TANGO Consortium of Taiwan, the University of Wisconsin at Milwaukee, Trinity College Dublin, Lawrence Livermore National Laboratories, IN2P3, University of Warwick, Ruhr University Bochum, Northwestern University, and former partners the University of Washington, Los Alamos National Laboratories, and Lawrence Berkeley National Laboratories. Operations are conducted by COO, IPAC, and UW.

This work has made use of data from the Asteroid Terrestrial-impact Last Alert System (ATLAS) project. The Asteroid Terrestrial-impact Last Alert System (ATLAS) project is primarily funded to search for near-Earth asteroids through NASA grants NN12AR55G, 80NSSC18K0284, and 80NSSC18K1575; byproducts of the NEO search include images and catalogues from the survey area. This work was partially funded by Kepler/K2 grant J1944/80NSSC19K0112 and HST GO-15889, and STFC grants ST/T000198/1 and ST/S006109/1. The ATLAS science products have been made possible through the contributions of the University of Hawaii Institute for Astronomy, the Queen’s University Belfast, the Space Telescope Science Institute, the South African Astronomical Observatory, and the Millennium Institute of Astrophysics (MAS), Chile.

This publication makes use of data products from the Near-Earth Object Wide-field Infrared Survey Explorer (NEOWISE), which is a joint project of the Jet Propulsion Laboratory/California Institute of Technology and the University of Arizona. NEOWISE is funded by the National Aeronautics and Space Administration.

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, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, 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.

DATA AVAILABILITY

The CSS light curve can be downloaded from https://crts.caltech.edu/; the ZTF and WISE light curves can be requested via https://irsa.ipac.caltech.edu; the ATLAS observations can be obtained from https://atlas.fallingstar.com/. The code and data are available in the online supplementary material.

Footnotes

REFERENCES

APPENDIX A: THE |$\gamma$|-RAY LIGHT CURVE OF PKS J2134−0153

The Fermi|$\gamma$|-ray Space Telescope was launched in 2008. One of its instruments is the Large Area Telescope (LAT: Atwood et al. 2009), which carried out long-term monitoring of the |$\gamma$|-ray sky over a wide range of |$\gamma$|-ray energies. The Fermi LAT Light Curve Repository (LCR: Valverde et al. 2022)⁸ contains a public database of variable Fermi-LAT source light curves on various time-scales. The 0.1–100 GeV |$\gamma$|-ray flux variations of PKS J2134−0153 are obtained from LCR. The light curves are rebinned with bin sizes of 3 days, 7 days (weekly), and 30 days (monthly), as shown in Fig. A1. The start and end times of the three light curves are from MJD 54691–60038 (5347 days).

APPENDIX B: THE CONSTRUCTION OF THE SYNTHETIC V-BAND LIGHT CURVE

The CRTS, ATLAS, and ZTF light curves (Fig. B1) are obtained by different telescopes with different filters. We aim to merge them to construct a |$\sim 6000$|-day (observed frame) long optical light curve. Hence, we must correct for the filter differences. Our approach is to convolve the SDSS spectrum of PKS J2134−0153 (Fig. B2) with the filter response curves and obtain the offsets between CRTS-V, ZTF-g, and ATLAS-c (as explained in Section 2.3). We use these three filters because their effective wavelengths are close to each other. Our approach is based on the implicit assumption that the optical colour of PKS J2134−0153 varies weakly. Indeed, the ZTF |$g-r$| colour of PKS J2134−0153 only changes by |$\sim 0.05$| mag, much smaller than the flux-variation amplitudes.