The Unanticipated Phenomenology of the Blazar PKS 2131–021: A Unique Supermassive Black Hole Binary Candidate

The 45.1 yr radio light curve of PKS 2131−021 is shown in Figure 1. It is immediately apparent to the eye that there are three distinct epochs, which are demarcated by the vertical black lines in Figure 1: In epoch 1 a strong 1.5-cycle periodic oscillation is seen in the Haystack data (green squares), of which the last half-cycle is also seen in the UMRAO data (brown triangles). There follows a 20 yr period (epoch 2) in which the oscillation of epoch 1 is not seen. In 2003 we enter the third epoch, in which a strong oscillation is again seen and which has very nearly the same period (within 2%) and also the same phase (within 10% of a cycle) as the periodic oscillation of epoch 1. This analysis by inspection is well corroborated by both the Lomb–Scargle (LS; Lomb 1976; Scargle 1982) and weighted wavelet Z-transform (WWZ; Foster 1996) analyses we carried out (see Section 3).

Zoom In Zoom Out Reset image size

Haystack Observations: The 15.5 GHz (λ1.9 cm) data (O'Dea et al. 1986) are shown by the green square data points in Figure 1. The bandwidth was 1200 MHz, and the primary beamwidth (FWHM) was 22. These data were obtained at roughly 1-month intervals with the 120 ft diameter telescope of the Haystack Radio Observatory located in Westford, MA. Details of the observing procedures (Dent et al. 1974b; Balonek 1982) are briefly summarized here. The observations were taken with the dual-feed beam-switched system, which reduces the effect of atmospheric fluctuations. Pointing corrections were applied, and the source antenna temperature was measured by an "on–off" procedure. Visual inspections of a chart record were used to identify observations affected by interference, which were then deleted. Corrections for atmospheric extinction and elevation-dependent gain were applied. The primary flux density calibrators were the compact H ii region DR 21, 3C 274, and 3C 123. A correction of a few percent was applied to account for partial resolution of DR 21 and 3C 274 by the Haystack observing beam.

UMRAO Observations: The UMRAO light curve is shown by the brown triangular data points in Figure 1. The long-term 14.5 GHz total flux density data for PKS 2131−021 were obtained with the 26 m equatorially mounted, prime-focus University of Michigan paraboloid as part of the UMRAO monitoring program that operated from the mid-1960s until 2012.5 and included observations of total flux density and linear polarization at three centimeter-band wavelengths for hundreds of blazars in a dedicated AGN program. A description of the general observation and calibration procedures is given in Aller et al. (1985). In brief, at 14.5 GHz dual, rotating, linearly polarized feed horns placed symmetrically about the prime focus were employed and alternated beams on the source; these fed a broadband uncooled high electron mobility transistor (HEMT) amplifier with a bandwidth of 1.68 GHz. Each observation shows the average from 1 day of a series of on–on measurements obtained by alternating between rotating feed horns during an interval of 30–40 minutes. The cadence was one of these averages per week. The adopted flux density scale (Baars et al. 1977) uses Cas A as the primary calibrator, accounting for its measured decay rate (Dent et al. 1974a). Observations of secondary calibrators were interleaved approximately every 2 hr throughout each 14.5 GHz observing run to monitor pointing and gain changes. The standard deviation associated with each daily flux density observation is computed from the system noise temperature and the number of individual on–on measurements made on the particular day. The standard error estimates shown include the effects of measurement noise, the errors introduced by uncertainties in the pointing corrections applied to the observations, and the uncertainties in determining the antenna gain as a function of time. The goal of the UMRAO program was to detect and follow blazar flares, so in general the cadence of the observations was set by the variability state (high cadence during flaring). However, as a member of the UMRAO BL Lac sample (Aller et al. 1999), PKS 2131−021 was also routinely observed approximately every 3 months. Additionally, during the last decade of operation, quasi-simultaneous observations (within a week) were matched to MOJAVE epochs for calibration of the Very Long Baseline Array (VLBA) data.

OVRO Observations: The OVRO light curve is shown by the blue dot data points in Figure 1. The OVRO data were part of the OVRO 40 m Telescope Monitoring Program (Richards et al. 2011). The telescope uses off-axis dual-beam optics in which the beamwidth (FWHM) is 157'. The cryogenic receiver uses an HEMT amplifier and is centered at 15 GHz with 2 GHz equivalent noise bandwidth. Gain fluctuations and atmospheric and ground contributions are removed with the double switching technique (Readhead et al. 1989), where one of the beams is always pointed on the source. Until 2014 May the two beams were rapidly alternated using a Dicke switch. In 2014 May a new pseudo-correlation receiver replaced the old receiver. Since then a 180° phase switch has been used to alternate the beams. A temperature-stable noise diode is used for relative calibration to compensate for gain drifts. The primary flux density calibrator is 3C 286 with an assumed value of 3.44 Jy (Baars et al. 1977), and DR 21 is used as a secondary calibrator source. Details of the observation and data reduction (Richards et al. 2011) cover the absolute calibration and the uncertainties, which include both the thermal fluctuations in the receiver and systematic errors that have been added in accordance with a rigorous procedure (Richards et al. 2011).

MRO Observations: The 22 and 37 GHz observations were made with the 13.7 m diameter Aalto University Metsähovi radio telescope, which is a radome-enclosed antenna with Cassegrain optics in Finland ($60^\circ 13^{\prime} 04^{\prime\prime}$ N, $24^\circ 23^{\prime} 35^{\prime\prime}$ E). The receivers have HEMT front ends operating at room temperature. The bandwidth at both frequencies is 1 GHz, and the beamwidths (FWHM) are 40 and 24 at 22 and 37 GHz, respectively. The observations are Dicke-switched ON–ON observations, alternating the source and the sky in each feed horn. A typical integration time to obtain one flux density data point is between 1200 and 1800 s. The detection limit of the system at 22/37 GHz is on the order of 0.2 Jy under optimal conditions, but it is heavily weather dependent. The flux density scale is set by observations of the H ii region DR 21. Sources NGC 7027, 3C 274, and 3C 84 are used as secondary calibrators. A detailed description of the data reduction and analysis is given in Teraesranta et al. (1998). The error estimate in the flux density includes the contribution from the measurement rms and the uncertainty of the absolute calibration.

PKS 2131−021 was observed with the VLBA at 15 GHz on numerous occasions between 1995 and 2012 as part of the 2 cm Survey (Kellermann et al. 1998) and MOJAVE (Lister et al. 2018) programs. The last eight maps made on this program, six of which overlap the OVRO monitoring observations, are shown in Figure 2. These show the basic core−jet structure of PKS 2131−021 on parsec scales. The results have been analyzed for 24 VLBA epochs (Lister et al. 2019; including 6 from other projects obtained from the NRAO data archive) by fitting Gaussian model features to the interferometric visibilities and tracking their evolution in flux density (Figure 3) and separation (Figure 4). The VLBA images have a restoring beam with typical FWHM dimensions ∼1 mas × 0.5 mas, but the positional determination accuracy of the individual model features is approximately 10 times better. The parsec-scale radio morphology of the source consists of a strong core feature with brightness temperature varying between 10^10.6 and 10^11.9 K and a jet structure that extends 2.5 mas (21 pc projected) to the east. The core accounts for roughly half of the total 15 GHz flux density of the source (Figure 3). The jet and core are linearly polarized, with fractional polarization levels increasing from 5% to 15% downstream. These values are typical of other core-dominated blazars in the MOJAVE sample (Lister & Homan 2005). Observations in 2006 (Hovatta et al. 2012) showed low Faraday rotation measure, with values ranging up to 131 rad m⁻², and the electric vectors are generally aligned with the jet direction, indicative of a magnetic field oriented perpendicular to the flow. The one exception was a jet feature (component 3) that had electric vectors oriented in the transverse jet direction.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 3 we show the light curves of the individual milliarcsecond jet components measured by MOJAVE (Lister et al. 2019). These overlap the 2005 peak in the sine-wave signal seen by UMRAO and the 2010 peak seen by both UMRAO and OVRO. Note that the core component and the component nearest to the core (component 5) are the components that produce the periodic signal. The sum of these two components is shown by the solid line in Figure 3. In Figure 4 we show the separations of the individual components from the core. Component 5 is stationary (μ = 3.0 ± 1.9 μas yr⁻¹; Lister et al.2019) at a separation of 0.3 mas (0.2 pc) from the core. We note that, while we might expect all of the periodic variability to arise in the core, in view of the small separation of component 5 from the core, it is not unlikely that the core and component 5 will vary in phase on the timescale of the observed periodicity.

Based on near-simultaneous OVRO and VLBI observations, the VLBI flux density is ≳98% of the single-dish flux density, so that, at most, ∼25 mJy of flux density is missing on small angular scales and is distributed on large scales. PKS 2131−021 has the typical blazar jetted-AGN morphology of a one-sided jet with a flat-spectrum core at one end of a steep-spectrum jet. The components in the jet have been modeled (Lister et al. 2019) with the results shown in Figures 3 and 4. The parsec-scale jet of PKS 2131−021 showed significant activity during the period between mid-1993 and 2000, with several features emerging from the core and moving downstream at apparent superluminal speeds. No new moving features appeared in the jet between 2000 and 2013. As can be seen in Figure 3, at all epochs except the first, the core dominates the flux density. From 2001 March 15 a stationary component (#5) is seen at a distance of 0.3 mas from the core, and together the core plus component 5 dominate the light curve from then on. The sinusoidal flux density fluctuations seen in epoch 3 are clearly due to these two components and dominated by the core. Thus, the periodicity we observe in PKS 2131−021 is coming primarily from unresolved components closer to the central engine and the base of the jet than the structures observed with 15 GHz VLBI.

Kinematic analysis (Lister et al. 2019) identified five distinct features in the jet, each having a different apparent speed. The innermost feature (id = 5) had no detectable motion over a 12 yr period, while the speeds of the other features ranged from 2.6c ± 0.2c (id = 4) to 20c ± 2.6c (id = 2). Most of the moving features display accelerated and/or nonradial motions on the sky that are exaggerated by projection effects. Given the maximum apparent speed, the jet viewing angle, θ, is constrained to lie within 57 of the line of sight, and an analysis based on the median core brightness temperature yields θ = 38 and a Doppler factor of δ = 14 (Homan et al. 2021).

The exact dates of emergence of the individual moving features are difficult to determine since that would require extrapolations of unknown accelerations close to the core. Also, in the case of components 3 and 4, there is blending with other jet emission that made it impossible to measure their positions when they were located near the core. Simple constant-velocity radial fits to components 1 and 2 give ejection dates of 1993.7 ± 0.4 and 1998.6 ± 0.5, respectively. For component 3, the first five epochs (when the feature was moving approximately radially away from the core) imply an ejection date of 1999.2 ± 0.3. In the case of the strongly accelerating component 4, its time of ejection can only be narrowed down to the interval 1997–2000.

The radio structure of PKS 2131−021 on different angular scales is shown in Figure 5. On arcsecond scales PKS 2131−021 has an unusual structure, which is most clearly seen in the 4.86 GHz VLA map of Rector & Stocke (2003)—reproduced here in Figure 5(b), which shows two linear features extending to about 8' from one side of the nucleus, and nothing on the opposite side of the nucleus. An EVN-MERLIN map by Cassaro et al. (2002) shows this bifurcated structure extending inward to within 100 mas of the core. Both gravitational lensing and precession have been suggested as possible causes for the strange structure. Another possibility, which has not previously been suggested to the best of our knowledge, are twin jets, such as are seen in 3C 75 (Owen et al.1985). It is interesting that two of the possible explanations for the large-scale structure of PKS 2131−021 (precession and twin jets) could well require an SMBHB.

Zoom In Zoom Out Reset image size

The MRO 22 GHz and 37 GHz observations are shown in Figure 6(a). We have cross-correlated the MRO 37 GHz observations with the OVRO 15 GHz light curve, with the results shown in Figure 6(b). There is clearly a correlation, as can also be seen by visual examination of the light curves, and it also appears that the 15 GHz variations lag those at 37 GHz, but it is hard to determine the lag—it could be anywhere between ∼0 and ∼200 days.

Zoom In Zoom Out Reset image size

We have extracted the Wide Field Infrared Explorer (WISE) data in the 2.8–3.8 μm and 4.1–5.2 μm bands. These are shown in Figures 6(c) and (e). The cross-correlation between the flux densities of these two bands with the radio light curve can be seen in Figures 6(d) and (f). The cross-correlation function is rather noisy, but it likely shows a real correlation when it is at, or just above, the 2σ level. If this is to be believed, then the radio variations lag the infrared variation by ∼250–350 days.

We observed PKS 2131−021 on 2021 October 6 (HST) with the Low Resolution Imaging Spectrometer (LRIS) mounted on the Keck I telescope. The observations were carried out in dark time, under photometric conditions with a median seeing disk FWHM of 05. A 900 s exposure was obtained at an airmass of 1.08, with the D560 dichroic, the 400/3400 blue-side grism, and the 400/8500 red-side grating (central wavelength 7830 Å). Spectral flux and telluric-absorption calibration was carried out using an observation of the standard star Feige 34. Data reduction followed standard procedures using the lpipe software (Perley 2019). The software was used to perform bias subtraction using the overscan levels, flat-fielding using dome-flat exposures, automatic cosmic-ray rejection, wavelength calibration using internal comparison-lamp exposures and sky emission lines, optimal sky-line subtraction, and optimal extraction of the spectral trace. No significant difference was observed between our spectrum of PKS 2131−021 and a previous spectrum obtained on 1995 November 20 (Rector & Stocke 2001). We confirm the redshift of 1.285, together with the BL Lac nature of the object.

The archival optical light curve of the blazar is presented in Figure 7. It consists of 119 V-band epochs from the Catalina Real-Time Transient Survey (CRTS; Drake et al. 2009) and 203 g-band epochs from the Zwicky Transient Facility (ZTF; Graham et al. 2019; Masci et al. 2019). The CRTS data were collected between 2005 May and 2015 December using the 0.7 m Catalina Schmidt Telescope located north of Tucson, Arizona. The ZTF data were taken between 2018 May and 2021 September with a wide-field camera mounted on the Palomar 48-inch Oschin (Schmidt) telescope.

Zoom In Zoom Out Reset image size

We searched for periodic variations in the optical data. We ran generalized Lomb–Scargle (GLS; Lomb 1976; Scargle 1982) and analysis of variance period-searching algorithms (Schwarzenberg-Czerny 1989) on the CRTS and ZTF data, and we did not find any significant periodicity in the range 1–1850 days. The optical light curve is not correlated with the radio data.

The WWZ is a Z-transform (Ragazzini & Zadeh 1952) that uses the technique of wavelet analysis to detect time-dependent periodicity in data (Foster 1996). We have performed a thorough WWZ analysis ¹⁹ of the light curves of the three epochs and all combinations thereof. These showed very clearly that there are three distinct epochs. This is most easily seen in the WWZ plot covering the whole 45.1 yr light curve shown in Figure 8. The power across this period varies strongly, and the raw WWZ plot (top panel) does not have the dynamic range to enable us to see the key features very clearly—more details are given in Appendix B. For this reason, we calculated time-normalized WWZ estimates by dividing the WWZ power in each 100-day time bin by the maximum power in the time bin in order to bring out the details of any changes in period (bottom panel). This WWZ plot supports the conclusion we arrived at by visual inspection of the light curve that there are three distinct epochs of activity in PKS 2131−021 over the 45.1 yr span of our observations, with a periodic oscillation of approximately the same period dominating in epochs 1 and 3 and a periodic oscillation of about twice that period dominating in epoch 2 (see Table 2). It is interesting to note that the boundary between epochs 1 and 2 is at the point where the WWZ power shifts from one period to another, whereas the boundary between epochs 2 and 3 is later than the point where the WWZ power shifts from one period to another. We ascribe the dominant variability in epoch 2 to red noise. It is possible that the transition from epoch 2 to epoch 3 occurred earlier, but we decided to retain the boundary we picked out by eye since the periodicity of epoch 3 is very clear from this point on and not nearly as clear before that, so that we suspect that all of the interval we demarcate as epoch 2 is dominated by red noise. Furthermore, we tried shifting the boundary between epoch 2 and epoch 3 between MJD 51,200 and MJD 53,800, and we found that it had little effect on the results, as shown in the notes to Table 1. Clearly, the precise choice of the boundary between epoch 2 and epoch 3 is immaterial for the purposes of this study.

Zoom In Zoom Out Reset image size

Table 2. Probabilities and Significance Levels of GLS Tests Computed from Simulations with Matched Red-noise Tail

Test	Test	GLS ${{ \mathcal P }}_{\mathrm{peak}}$	Period Range (ΔP)	Total	Number of Simulations	p-value	Significance
Number		(max = 1)	(days)	Simulations	That Pass Test		(σ)
1.1	epoch 1 ${{ \mathcal P }}_{\mathrm{peak}}$, ${p}_{\mathrm{sim}}\leqslant {p}_{\mathrm{peak}}$	0.71	All	10,000	1446	1.45 × 10−1	1.06
1.2	epoch 2 ${{ \mathcal P }}_{\mathrm{peak}}$, ${p}_{\mathrm{sim}}\leqslant {p}_{\mathrm{peak}}$	0.76	All	10,000	632	6.32 × 10−2	1.53
1.3	epoch 3 ${{ \mathcal P }}_{\mathrm{peak}}$, ${p}_{\mathrm{sim}}\leqslant {p}_{\mathrm{peak}}$	0.81	All	100,000	40	4.0 × 10−4	3.35
2	epoch 1 ${{ \mathcal P }}_{\mathrm{lim}}$, ΔPepoch 3	≥0.50	1661.9–1858.9	10,000	197	1.97 × 10−2	2.06
3	1.3+2	⋯	⋯	⋯	⋯	7.88 × 10−6	4.32
4	3+phase	⋯	⋯	⋯	⋯	1.58 × 10−6	4.66

Note. ${ \mathcal P }$ is the GLS power, and ΔP is the range of periods included in the test. In Tests 1.1, 1.2, and 1.3 we count simulations at all periods with p-values less than the p-value of the peak in that epoch. Since the significance of the peak in epoch 3 was so high, we needed 100,000 simulations to determine the significance. The low individual significance levels for epochs 1 and 2, shown in Tests 1.1 and 1.2, indicate that they are easily produced by red noise. Test 2 is the number of simulations of epoch 1 that were seen in the ±3σ wide period window centered on the epoch 3 period (P = 1760.4 days).

Download table as:  ASCIITypeset image

Since the results of the WWZ analysis were entirely consistent with those from our GLS analysis and least-squares sine-wave analysis, we do not reproduce further details of our WWZ analysis here other than to give the periods associated with the maximum power level in our WWZ analysis for the three epochs (see the notes to Table 1).

We carried out a GLS periodogram analysis (Lomb 1976; Scargle 1982) using the GLS () periodogram ²⁰ (Zechmeister & Kürster 2009) on each epoch of the light curve of PKS 2131−021, as well as all combinations thereof. In order to distinguish it from the period of the periodicity, P, we use the symbol ${ \mathcal P }$ to denote the GLS power. Following Scargle (1982), we evaluate the GLS at frequencies between f_low = 1/T, where T is the total time of the corresponding epoch, and f_high = N/(2T), corresponding to one over twice the average time interval between measurements, with N the number of data points. The selection of a pseudo-Nyquist frequency as the highest frequency does not have any effect on the results, as all relevant peaks are found at much lower frequencies. We sample this frequency range uniformly in steps of Δf = 1/(ζ T). For our analysis we conservatively chose ζ = 10, since with this value we can be confident that we have adequately sampled periods up to the duration of the observations. For comparison, Scargle (1982) used ζ = 5.

We also tested the effect of the frequency grid spacing with coarser and finer grids, ζ = 5 and ζ = 20, and found only slight differences compared to the values for ζ = 10. The variations with ζ were far too small to affect our conclusions in the analysis presented here.

The results of the GLS analyses of epochs 1, 2, and 3, both independently and combined, are shown in Figure 9 and Table 2. There are about 1.5 cycles of the periodic variation in the Haystack+UMRAO data and about 3 in the UMRAO+OVRO data. While this is a small number of cycles, our simulation procedure takes long-term fluctuations into account via the red tail in the spectrum, and so we are confident that our significance levels are robust. Clearly, epoch 3 has the statistically strongest periodicity, and the similarity of the period with that of epoch 1, as well as the phase correspondence, which are what immediately caught our attention when we first saw the Haystack observations, when taken together, are highly significant, as we will show.

Zoom In Zoom Out Reset image size

In the combined GLS analysis in the bottom panel of Figure 9 the largest peak corresponds to the epoch 3 peak, the second-largest peak corresponds to the epoch 1 peak, and the third-largest peak corresponds to the epoch 2 peak. All three of these peaks in the combined analysis of the bottom panel are slightly shifted in period relative to their period values in the individual epoch GLS analyses. Given that the epoch 3 peak is the only peak detected at greater than 3σ significance, and that the spectrum has a powerful red-noise component, the peaks other than the epoch 3 peak should not be regarded as physically significant unless supported by other data, as is true for the epoch 1 peak, but not for any of the other peaks seen in the combined analysis. This, plus the WWZ analysis presented in Figure 8, illustrates the potential danger of carrying out GLS, or WWZ, analyses on astronomical light curves without adequately simulating the variability in the object and, in addition, applying rigorous statistical criteria.

We fitted the light curves in the three epochs with a sine-wave function:

$$\begin{eqnarray}&&{S}_{\mathrm{model},i}=A\sin ({\phi }_{i}-{\phi }_{0})+{S}_{0},\end{eqnarray} \tag{ 1 }$$

where ϕ_i = 2π(t_i − t₀)/P is the phase of the ith data point. Our model has five free parameters: P—period, A—amplitude, ϕ₀—phase of the sine wave, S₀—mean flux, and σ₀—characteristic amplitude of an intrinsic AGN variability. We keep t₀ = 51,000 fixed (which is in the midpoint between the start of epoch 1 and end of epoch 3 observations). We find the best-fitting parameters by maximizing the following likelihood function:

$$\begin{eqnarray}&&\mathrm{ln}{ \mathcal L }=-\displaystyle \frac{1}{2}\sum _{i=1}^{N}\left(\displaystyle \frac{{({S}_{i}-{S}_{\mathrm{model},i})}^{2}}{{\sigma }_{i}^{2}+{\sigma }_{0}^{2}}+\mathrm{ln}({\sigma }_{i}^{2}+{\sigma }_{0}^{2})\right).\end{eqnarray} \tag{ 2 }$$

The uncertainties are estimated using the emcee sampler (Foreman-Mackey et al. 2013) and represent the 68% confidence range of the marginalized posterior distribution. The σ_i are taken from the observations, assuming independent Gaussian errors, and the intrinsic variability is represented by Gaussian white noise of variance ${\sigma }_{0}^{2}$. We have not used the red-noise spectrum because to include it in the fit would mean including the full covariance matrix rather than just the diagonal term. We know that including a red-noise term of the type seen in epoch 2 can change the period at the ∼10% level over periods of a few years (see Section 3.5).

Fit results, separately for epochs 1, 2, and 3, are presented in Table 1. Note that the uncertainties of ϕ₀ are relatively large because there are strong degeneracies between period and phase, especially as the adopted t₀ is relatively far from the actual observations. If we chose t₀ to be closer to epoch 1 (epoch 3) observations, the respective uncertainties of ϕ₀ would be much smaller. The degeneracy between the period and phase reflects the fact that, over time, information about the phase is lost.

We also simultaneously fitted both epoch 1 and epoch 3 data by fitting a single wave, with the same period and phase extending across both epochs, but permitting the amplitude and offset between epoch 1 and epoch 3 to vary. Results of this joint fit are presented in the last column of Table 1.

Prior to MJD 45,500 (epoch 1), the strong sinusoidal variation seen by Haystack+UMRAO has period P₁ =1729.1 ± 32.4 days. In epoch 2 (45,500 < MJD < 52,850), the UMRAO data show sinusoidal variability, but no periodicity at the frequency or in phase with that seen in the first and third epochs. We therefore think that this is a red-noise phenomenon, and we will show that it fails our 3σ criterion by a substantial margin. In epoch 3 (MJD > 52,850) strong sinusoidal variation returns with period P₃ = 1760.4 ± 5.3 days. In addition, the phase of the sinusoidal variation from epoch 3, when extrapolated back to epoch 1, matches to within 10% of P, or 20% of P/2, which is the relevant number here since we associate the extrapolated curve with the nearest peak. Combining the Haystack, UMRAO, and OVRO data from epoch 1 and epoch 3 yields a period P₁₃ = 1737.9 ± 2.6 days, i.e., the fractional uncertainty in period δ P/P ∼ 1.5 × 10⁻³.

To calculate the significance of the power levels observed in epochs 1, 2, and 3, we count the number of simulations in which the p-value, ${p}_{\mathrm{sim}}$, of the strongest peak in the simulation is less than the observed peak p-value, p_peak, for that epoch. The procedure is explained in more detail in Appendix A. The results are shown in Table 2. We see that the peaks in epoch 1 and epoch 2 are significant only at the 1.06σ and 1.53σ levels, respectively (Tests 1.1 and 1.2). However, the peak in epoch 3 is significant at the 3.35σ level (Test 1.3). This makes it clear that PKS 2131−021 has a periodicity in epoch 3 that is not due to red noise, which must therefore be associated with some physical mechanism in the object.

We next turn to the good agreement in the periods, straddling the 20 yr gap in periodic fluctuations, seen in epochs 1 and 3 (see Table 1), which, to our knowledge, is unprecedented in observations of AGNs. The question we address is this: having observed the periodicity in epoch 3, which is what first drew our attention to this source, what is the probability of observing the same periodicity in epoch 1, at an acceptable LS power level and within the 3σ period window of our epoch 3 period?

We first determine what we deem to be an acceptable GLS power level, ${{ \mathcal P }}_{\mathrm{lim}}$. From Table 2, we see that, of the three periodicities we have been discussing, the lowest power level is that of the peak in epoch 1, at ${{ \mathcal P }}_{\mathrm{peak}}=0.71$. From Figure 9 we see that, apart from the largest peak in each epoch, all the other peaks in power are below the ${ \mathcal P }=0.4$ power level. We therefore select a threshold level of power of ${{ \mathcal P }}_{\mathrm{lim}}=0.5$ as a power level suitably above the level of sidelobes in the GLS spectrum and suitably low that we do not miss any strong sinusoidal features in the light curves. We will therefore count all simulations that give an GLS power level ${{ \mathcal P }}_{\mathrm{lim}}\geqslant 0.5$ in epoch 1, in our test of the significance of the epoch 1 periodicity, given the epoch 3 periodicity.

We now turn to the range of periods to be considered. From Table 1 we see that the uncertainty in the difference between P₁ and P₃ is $\sqrt{{32.4}^{2}+{5.3}^{2}}=32.8$ days. Thus, the 3σ uncertainty in the difference is 98.5 days. Since we wish to calculate the probability of observing the periodicity in epoch 1, given that it had been observed in epoch 3, we will count simulations in epoch 1 where the strongest peak lies in the period range P₃ ± 98.5 days, i.e., from 1661.9 to 1858.9 days. So our two criteria are GLS power ${ \mathcal P }\geqslant {{ \mathcal P }}_{\mathrm{lim}}=0.5$ and 1661.9 days < P < 1858.9 days.

The probability of observing in epoch 1 a peak within the ±3σ window of the periodicity observed in epoch 3 and with ${ \mathcal P }\geqslant 0.5$ is shown in Test 2 of Table 2 to be 1.97 × 10⁻², which is significant at the 2.06σ level. We may now combine Test 1.3 and Test 2 to determine the probability of observing both the observed periodicity in epoch 3 and a periodicity in epoch 1 with period within the ±3σ window of the periodicity found in epoch 3. This is determined in Test 3 by multiplying p(Test 1.3) × p(Test 2) = 7.88 × 10⁻⁶, which is significant at the 4.3σ level. When we add the correspondence in phase to within one-fifth of a half-cycle, the probability drops to 0.2 × p(Test 3) = 1.58 × 10⁻⁶, which is significant at the 4.6σ level.

We believe that these significance estimates, based as they are on simulations that have the same PSD and the same pdf as the observed light curves in epochs 1, 2, and 3, and therefore with the same red-noise tail as PKS 2131−021, are robust, and therefore that this analysis demonstrates conclusively that the periodicities observed in epochs 1 and 3 are connected and that there is a physical process that maintains this period over our 45.1 yr observing period even when it is not manifested in the light curve.

From Table 1 we see that the observed difference in period between the epoch 3 data and the epoch 1 + epoch 3 data is significant at the 3.8σ level, so this might indicate that the period is slowly changing. Separate sine-wave fits to the first half and the second half of the OVRO data show a difference of ∼10% in period. PKS 2131−021 is a blazar and, as such, highly variable even in the absence of the periodicity we have been discussing, as is easily seen from the deviations from sinusoidal variations observed in epochs 1 and 3 and throughout epoch 2. Thus, we ascribe the apparent changes in period to random red-noise variability in this blazar, such as that seen in epoch 2. Indeed, because of the strong red-noise component in their variability, such variations in the observed period in radio light curves of an SMBHB in a blazar like PKS 2131−021 are to be expected, but these should average out in the long term.

The reader may have noticed the apparent harmonic relationship between the epoch 2 and epoch 1+3 periodicities. In this paper we do not wish to draw undue attention to this for two reasons:

• 1.  
  The epoch 2 periodicity has significance of only 1.53σ, as shown in Table 2. We do not discuss this periodicity since we think any discussion runs counter to the approach we are adopting in this paper. We encourage workers in this field to carefully assess red noise before entering discussions of apparent periodicities.
• 2.  
  There is no obvious phase relationship in between the real, highly significant, periodic signals we see in epochs 1 and 3 and the low-significance apparent periodicity in epoch 2.

We reemphasize here our deliberate decision to consider only signals that are demonstrably not due to red noise, i.e., only those having significance greater than 3σ, unless supported by other evidence.

3.4.1. The "Look Elsewhere" Effect

The fits of the sine waves shown in Figure 1, while good, are not perfect. There is a clear epoch from 2015 to 2017 when the sine wave is systematically above the data, and similarly in 2020 the data are systematically above the sine wave. If the two halves of the OVRO data are analyzed separately, the periods differ by ∼10%. Does this mean that PKS 2131−021 is simply another QPO? The unique properties we have drawn attention to above suggest that it is not.

We now provide an illustration, based on the variability during epoch 2, that refines the period derived for epoch 3 to P = 1737.6 ± 3.6 days, shown by the gray curve in Figures 10(a) and (b), in good agreement with the period and phase observed over the whole 45.1 yr. This is merely an example of what must be happening, under our hypothesis of an SMBHB. We go through the following argument to illustrate the effect that we think the nonperiodic varying components must be having on the observations. Of course, this does not prove that PKS 2131−021 is an SMBHB. But the exercise is illuminating, and in our view it does therefore strengthen the case for an SMBHB on the grounds of plausibility.

Zoom In Zoom Out Reset image size

We have fitted the period when the sinusoidal variations we are investigating were absent using the epoch 2 data between MJD 45,500 (1983 June 15) and MJD 52,850 (2003 July 30) and found that a 6° polynomial, shown by the brown curve in Figure 10(a), follows the slowly varying component well, whereas lower-degree polynomials do not. We adopted the least-squares sine-wave fit to epoch 1 + epoch 3 and subtracted this, adjusted for the amplitudes of the Haystack and OVRO periodic signals, from the OVRO data to give the residual light curve shown by the black points in Figure 10(a). This shows the flux density variations in PKS 2131−021, under the SMBHB hypothesis and in the absence of the sinusoidal signal. We then adopted the same approach with the OVRO data as with the UMRAO data and fitted a 6° polynomial to the residual light curve to give the black curves shown in Figures 10(b) and (c). These show what the slowly varying components of the source were doing, under the SMBHB hypothesis and apart from the periodic signal. Correcting for the slowly varying components, we derive the gray points shown in Figure 10(c). The sine-wave fit to the corrected OVRO data is shown by the gray curve in Figures 10(a) and (b), which has period P = 1737.6 ± 3.6 days. Thus, as expected, it has a period closer to the period of the epoch 1 + epoch 3 sine-wave fit (P = 1737.9 ± 2.6 days) than the original sine-wave fit to the epoch 3 data (P = 1760.4 ± 5.3 days). In addition, and as expected, the phase during epoch 1, as shown by the gray curve in Figure 10(a), is much closer to that of the joint epoch 1 + epoch 3 fit (black curve) than the original epoch 3 fit (blue dotted curve). This example simply shows how small, slow variations in the nonperiodic signal can have a significant effect on the period and phase of the sinusoidal variations in the periodic signal in PKS 2131−021, a point that, while obvious, is worth investigating to illustrate the levels of correction that are needed to achieve much better coherence. The fractional changes in the flux density required for perfect coherence, which would be grossly overfitting the data, are shown by the black points in Figure 10(c). Note that they are nearly all <10%.

Itemization of the above procedure: ²¹

• 1.  
  We fitted separate least-squares sine waves to the epoch 1 and epoch 3 data.
• 2.  
  We subtracted these sine waves from the Haystack data of epoch 1 and the OVRO data of epoch 3.
• 3.  
  We assumed that the sinusoid represented all of the SMBHB emission, and therefore that the SMBHB signal was now removed from the Haystack data and the OVRO data.
• 4.  
  We assume that epoch 2 is already devoid of SMBHB signal. Thus, in Figure 10 the black data in epoch 1, the brown data in epoch 2, and the black data in epoch 3 represent the data with the SMBHB sinusoid subtracted.
• 5.  
  We fitted a six-order polynomial to the black data in epoch 3 in Figure 10. The six-order polynomial is thus assumed to represent the red-noise variation of the SMBHB.
• 6.  
  We then subtracted the polynomial fit given by the black line from the original data, which gives the gray points.
• 7.  
  Finally, we did a least-squares sine-wave fit to obtain the period for the gray points.

Summary: Our motivation in carrying out the above procedure is as follows. We have, quite deliberately, "made the data fit the SMBHB model" in order to provide an illustration of what we think the nonsinusoidal random flux density variations very likely were during the OVRO observations, and then to compare these to the nonsinusoidal random variation during epoch 2. This does not prove that PKS 2131−021 is an SMBHB, but, given the significance of the similarities in period and phase in epochs 1 and 3, as shown in Table 2, this is almost certainly the case. Thus, understanding and anticipating the effects of the nonsinusoidal variability in a jetted SMBHB is bound to be important as this field opens up.

It is interesting to note that the periodic fluctuations occur at the times of highest flux density in the overall light curve, so that it appears that the periodic signal is in addition to the usual flux density level in this source. Thus, on our SMBHB model we are proposing that PKS 2131−021 has an extremely stable period and that the small changes in period that we see are due to other sources of variation that distort the pure sine-wave variation. We also propose that the emission associated with the periodic variability can turn off and on in a time that is short compared to the period, and that the amplitude of the periodic signal can vary significantly between epochs of periodic emission.

Our observations thus lead to three strong predictions based on the SMBHB hypothesis: (1) we expect the future long-term averaged period to agree with our measured period spanning 45.1 yr of observations (P = 1737.9 ± 2.6 days); (2) in the short term we expect to see small (∼10%) variations in the period, due to the corrupting effects of the other varying emission components in this object; and (3) we expect the SMBHB sinusoidal signal to continue to appear and disappear at random times.

It is clearly important to investigate the stability of the periodic variations seen in PKS 2131−021, since this has potential implications for the SMBHB hypothesis, under which the period should decrease with time. For the purposes of comparison we will denote the periods determined in our least-squares sine-fitting analysis in the different epochs (see Table 1) by P₁ = 1729.1 ± 32.4 days, P₂ = 3779.1 ± 46.0 days, P₃ = 1760.4 ± 5.3 days, and P₁₃ = 1737.9 ± 2.6 days.

The most interesting of these is P₁₃ since it represents a single least-squares sine-wave fit across the 45.1 yr span. On the hypothesis that PKS 2131−021 is an SMBHB, this represents our best estimate of the observed orbital period. On the SMBHB hypothesis we therefore predict that, in future long-term monitoring, this period will be maintained until the effects of gravitational radiation produce a noticeable reduction in the period.

We therefore predict on the SMBHB hypothesis that the next minimum is occurring now (2021 October), and we will be able to confirm this in the next 6 months. The next maximum is predicted in 2024 February, and the following minimum in 2026 July. It should be possible to confirm or disprove the SMBHB hypothesis within the next 5–10 yr, provided that we do not go back into an epoch 2–like phase, in which this periodicity vanishes.

We now consider the stability of the period as can best be determined from the observations. We have already seen that ∼10% variations occur over time spans of a few years, which we ascribe to the corrupting effects of the variability in PKS 2131−021 that are unrelated to the underlying periodicity we seek to study. Taking their associated uncertainties into account, we see that δ P₁₃ = P₃ − P₁ = 31.3 ± 32.8 days. So the 3σ window is δ P₁₃ = −67.2 → + 129.8 days.

On the SMBHB hypothesis we are interested in negative δ P₁₃, and in this case our fractional 3σ limit is δ P₁₃/P₁₃ =−67.2/1737.9 = −3.9 × 10⁻². The midpoint of epoch 1 is MJD = 44,109.8, and that of epoch 3 is MJD = 56,023.6, so the span between the midpoints at which P₁ and P₃ were measured is 11,913.8 days, or δ t = 32.6 yr. Thus, the 3σ upper limit on the fractional period decrease per year is δ P₁₃/P₁₃ δ t = − 1.19× 10⁻³ yr⁻¹. This is the maximum rate of fractional period decrease per year consistent with our observations. From Maggiore (2008), using our rest-frame quantities, we derive a chirp-mass upper limit of 5.4 × 10⁹ M_⊙. This may be compared with the slightly lower upper limit on the chirp mass derived from the GW limit discussed in Section 5.

In this section, before presenting our model for PKS 2131−021 and its implications, we summarize briefly our approach and the logical thread of our statistical tests. The key steps in our approach are as follows:

• 1.  
  The recognition of the role of red noise in the γ-ray, optical, and radio light curves of blazars, as exemplified in the work of Vaughan et al. (2016) and Covino et al. (2019). As a consequence, we require strong statistical evidence that any apparent periodicity is not simply the result of red noise in the light curve, and we apply a 3σ threshold that must be met by any apparent periodicity before it should be discussed. Without such a threshold, workers in this field will expend a great deal of time and energy seeking physical explanations for apparent periodicities that are the result of random fluctuations and have no other physical significance. Consideration of such phenomena will obscure the physics and hinder progress in the field.
• 2.  
  The recognition that there are three distinct epochs in the 45.1 yr radio light curve of PKS 2131−021, as indicated by the periodicities seen in Figure 1. The WWZ analysis of Section 3.1, including Figure 8, and the analysis of different MJD dates for the transition from epoch 2 to epoch 3 given in Section 3.1, and in the notes to Table 1, in addition to our original visual inspection, are our justification for this.
• 3.  
  Starting with epoch 3, the epoch that initially drew our attention to the periodicity in PKS 2131−021, we carried out a GLS analysis making no assumptions about periodicities, and we considered all periodicities via the "look elsewhere" effect. We found that the probability at, or above, the observed power level in the GLS analysis is 4 × 10⁻⁴, i.e., it is significant at the 3.35σ level, and thus it is above the threshold of point #1 above.
• 4.  
  We then investigated the sinusoidal periodicity of epoch 1 as follows: we selected a ±3σ-wide periodicity window, based on the uncertainties in the least-squares periodicities in epoch 1 and 3, and determined the probability of a chance periodicity at or above the GLS power observed in epoch 1 in the window of this width and centered on the periodicity of epoch 3. This probability is 1.97 × 10⁻², i.e., it is significant at the 2.06σ level.
• 5.  
  Since the data from epochs 3 and 1 are independent, the probability of the agreement in periods observed in epoch 1 and epoch 3 is given by the product of the two probabilities and is 7.88 × 10⁻⁶, i.e., it is significant at the 4.32σ level.
• 6.  
  Finally, we took into account the agreement, to within 20% of half a period, of the phase of the epoch 3 periodicity extrapolated to epoch 1. This has a random probability of 0.2. Thus, the probability of the phase agreement in addition to the periodicity agreement between epoch 1 and epoch 3 is 1.58 × 10⁻⁶, i.e., it is significant at the 4.66σ level. This leaves us in no doubt that the periodicity observed in epochs 1 and 3 is a significant physical property of PKS 2131−021 and is not simply a random manifestation of red noise.