A Supermassive Binary Black Hole Candidate in Mrk 501

Abstract

Using multifrequency observations, from radio to $\gamma$-rays of the blazar Mrk 501, we constructed their corresponding light curves and built periodograms using RobPer and Lomb–Scargle algorithms. Long-term variability was also studied using the power density spectrum and the detrended function analysis. Using the software VARTOOLS Version 1.40, we also computed the analysis of variance, box-least squares and discrete fourier transform. The result of these techniques showed an achromatic periodicity ≲$229\mathrm{d}$. This, combined with the result of pink-color noise in the spectra, led us to propose that the periodicity was produced via a secondary eclipsing supermassive binary black hole orbiting the primary one locked inside the central engine of Mrk 501. We built a relativistic eclipsing model of this phenomenon using Jacobi elliptical functions, finding a periodic relativistic eclipse occurring every ∼$224\mathrm{d}$ in all the studied wavebands. This implies that the frequency of the emitted gravitational waves falls slightly above $0.1$ mHz, well within the operational range of the upcoming LISA space-based interferometer, and as such, these gravitational waves must be considered as a prime science target for future LISA observations.

1. Introduction

Among active galactic nuclei (AGN), blazars are objects that emit variable non-thermal radiation throughout the electromagnetic (EM) spectrum [1] with their jets pointing at an angle of no more than ∼$30$ from the observer’s line of sight [2]. These extragalactic sources present total luminosities in the range of ${10}^{41}$–${10}^{47}\mathrm{ergs}/\mathrm{s}$ (see, e.g., [3]).

The light curve variabilities in blazars are commonly classified as follows: (a) intraday variability (IDV), corresponding to periods of over a day or less [4]—they are also called intra-night variability or micro-variability [5]; (b) short-term variability (STV), which corresponds to variability over days to several weeks; and (c) long-term variability (LTV) that takes place on timescales of months to years [6].

These variabilities have many explanations that are related to whether the source is jet-dominated or not. In the first case, Marscher and Travis [7] suggest that the variability in compact jets is justified because of the non-thermal emission in blazars. Furthermore, variability studies in the radio and optical wave bands by Camenzind and Krockenberger [8] conclude that shock waves in relativistic collimated flows are responsible for the observations. This idea was further studied by Mohan and Mangalam [9], proposing a general relativistic model of jet variability in AGN, incorporating orbiting blobs in a helical motion along a magnetic surface near the black hole. In this direction, de la Cruz-Hernández and Mendoza [10] interpreted that shock wave emissions inside the jet are caused by a periodic injection velocity of the flow at the base of the jet.

For non-jet-dominated sources, McHardy et al. [11] suggest that the differences in lag between different bands indicate that the variability is produced by the disk. Also, Edelson et al. [12]’s account of rapid luminosity changes that indicates emission regions confined to the inner disk or corona. In this sense, Uttley et al. [13] suggest that the variability is due to processes and accretion rate variations. Also, Stella and Vietri [14] explain that the variability can be attributed to the relativistic dragging of inertial frames around a rapidly rotating disk.

A famous example of a non-jet-dominated source is represented by the quasar OJ 287, which was extensively studied by Lehto and Valtonen [15], who found sharp flares within major outbursts of the optical light curve. The authors proposed a model in which a smaller black hole crosses the accretion disk of a larger black hole during its binary orbit. An extensive analysis of its optical light curve was used to infer this supermassive black hole binary system [16]. The periodicity of this blazar was discovered by analyzing its historical optical light curve, which contains data from more than 100 years, showing repeated bursts at intervals of about $11.65\mathrm{y}$. The best-known model of this periodicity was constructed by Lehto and Valtonen [15], and it consists of a primary black hole—the central engine, with a mass of ∼$17\times {10}^9{\mathrm{M}}_{\odot }$, surrounded by an accretion disk and a secondary black hole with a mass of ∼${10}^8{\mathrm{M}}_{\odot }$, orbiting the primary and intersecting the accretion disk on each orbit, causing tidally induced mass fluxes from the accretion disk to the primary black hole.

Of particular relevance to the studies carried out in the present article is the case for postulating the existence of a secondary supermassive black hole orbiting a primary one. The first proposal of this kind was made by Begelman et al. [17] in order to account for periodic or quasi-periodic oscillations.

An extensive survey to find periodic light curves in optical light was carried out by Graham et al. [18]. Of the $247,000$ studied light curves, it was found that the one corresponding to PG 1302-102 shows a periodicity of $1884\pm 88\mathrm{d}$. The authors assumed that this was due to the existence of a secondary black hole “eclipsing” the primary one, and they concluded that the system is separated by less than a parsec. More recently, Tavani et al. [19] found 2.2 yr QPOs in the $\gamma$-rays band of the blazar PG 1553+113 and once again proposed it as a supermassive black hole binary system.

Relevant periodicities of other AGN appear in the literature. Quite important to mention is the work by  Li et al. [20], who report long-term variability of ∼14 yr in the optical continuum of the nucleus of NGC 5548. For this same object, Bon et al. [21] found a periodicity ∼43 yr. Also, Li et al. [22] studied a possible ∼20 yr periodicity in long-term optical photometric and spectral variations of the nearby radio-quiet Active Galactic Nucleus Ark 120, and Chen et al. [23] published a sample of quasar candidates with periodic variations from the Zwicky Transient Facility.

A very interesting and quite well studied blazar is Markarian 501 (Mrk 501). It is a BL Lac object with several periodicities reported in the literature: (1) A periodicity of $23\mathrm{d}$ was reported during a flare detected in X-rays and $\gamma$-rays and modeled as a supermassive binary black hole (see, e.g., [24] and references therein). (2) In the same frequencies, a periodicity of $72\mathrm{d}$ was found by Rödig et al. [25]. (3) A periodicity of $630\mathrm{d}$ was discovered by Wang et al. [26] in X-rays. (4) Finally, Bhatta [27] found a $332\mathrm{d}$ periodicity in the Fermi-LAT $\gamma$-rays light curve. All of these reported periodicities are not achromatic, and they represent different databases in time. Although interesting, they are not relevant to the study presented in this article, which used long-term databases in four different frequencies.

Mrk 501 has a redshift of $z=0.034$ (∼$456\mathrm{Mly}$∼$140\mathrm{Mpc}$) with R.A. $={16}^{\mathrm{h}}{53}^{\mathrm{m}}52.2^{\mathrm{s}}$, Dec. = $+39°{45}^{\prime }{37}^{″}$. It was discovered by Quinn et al. [28] using the Whipple Imaging Atmospheric Cherenkov Telescope (IACT). It has been monitored since 1996 in various frequencies: radio [29], optical [30], and $\gamma$-rays [31,32].

In this article, we report a mean periodicity of $224.07\pm 0.22\mathrm{d}$ in multi-frequency observations of Mrk 501. The dataset in radio was obtained from the Owens Valley Radio Observatory (OVRO) (https://sites.astro.caltech.edu/ovroblazars/, accessed on 27 June 2020), the optical dataset from the American Association of Variable Star Observers (AAVSO) (https://www.aavso.org/, accessed on 12 September 2021), the X-rays dataset is from the Neil Gehrels Swift Observatory (Swift) (https://swift.gsfc.nasa.gov/, accessed on 8 October 2021) and for $\gamma$-rays, the data were taken from the Fermi Gamma-Rays Space Telescope (FGST, also FGRST) (https://fermi.gsfc.nasa.gov/, accessed on 13 March 2020). These datasets and their corresponding processing (reduction) is explained in Section 2. In Section 3, we describe different methods to find this multifrequency periodicity, and in Section 3.4, we model this periodic behavior, assuming a supermassive binary black hole using Jacobi elliptical functions, which offer good representations of eclipses that produce occultations (as they occur for binary stars or exoplanets eclipsing their central star) and magnifications (such as the ones that are produced by massive relativistic objects that bend light and magnify its intensity). Finally, in Section 4, we discuss our results.

2. Data and Light Curves

The radio dataset from 22 January 2009 to 27 June 2020 was obtained from the OVRO database. OVRO consists of a $40\mathrm{m}$ telescope with a cryogenic receiver at a central $15\mathrm{GHz}$ frequency, a $3\mathrm{GHz}$ bandwidth and two symmetric off-axis beams. This observatory has been monitoring blazars since 2008 [29], and one of its main targets is the search for QPOs and correlations between radio and $\gamma$-rays in blazars [33,34]. The signal-to-noise level reported by the OVRO database is such that it produces a systematic flux uncertainty of about 5%.

The optical AAVSO database is public, and it offers long-term datasets. The institution is an international organization of variable star observers who participate in scientific discovery through variable-star astronomy. It was founded in 1910, and its observations of variable stars are collected and archived for worldwide access in collaboration with amateur and professional astronomers. Observations with errors of ≥1.0 magnitudes are rejected by the AAVSO community. An error of 1.0 magnitude represents a signal-to-noise ratio of 1, making it statistically insignificant. The light curve of Mrk 501 was built using the database from 24 June 1998 to 12 September 2021.

The optical and radio data reductions were processed via AAVSO and OVRO, respectively. The data were obtained by consulting the public databases of both observatories. For details of the reduction, calibration, correlation processes, etc., of the radio database, see Richards et al. [29], and for AAVSO, see Kinne [35].

For the X-ray light curve, we used the Swift database from 2 October 2008 to 8 October 2021. This dataset contains energies in the range of 0.3–10 keV. Swift has an X-ray telescope (XRT) with two important characteristics that makes it important for observations: a low background and a constant point spread function across the field of view [36].

The bin size used for obtaining the data was the default time of 5 s, and the command line interface used was xselect, which works in conjunction with the Fermitools software version v11r5p3, specifically designed for astrophysical X-ray analysis. It offers convenient functions to organize input data using the observation catalog, applies various filters to the event data (such as intensity, phase, region, etc.) and creates good time intervals based on a chosen selection criteria. It serves as a valuable tool to work with X-ray data, providing an efficient and customizable way to manage and analyze astrophysical data.

In addition to its major functions, xselect also provides three commands that correspond to different time systems: universal time (UT), modified Julian days (MJD), and the SpaceCraft Clock. We constructed the X-ray light curve using MJD.

The Fermi public database of $\gamma$-rays has fluxes in the range of $100\mathrm{MeV}$–$300\mathrm{GeV}$. For Mrk 501, it covers a time interval from 4 August 2008 up to 13 March 2020. The analysis was performed using the public Fermi/LAT data corresponding to the P8R3 SOURCE $\gamma$-ray event selection within a 15-degree range of search. To ensure data quality, only events corresponding to good time intervals with DATA_QUAL > 0 and LAT_CONFIG == 1 were retained, and a maximal telescope zenith angle of 90 degrees was applied. Data reduction was performed using the Fermitools package v2.0.8. Galactic and extragalactic diffuse emission was taken into account using the gll_iem_v07.fits and iso_P8R3_SOURCE_V2_v1.txt models, respectively. The spectral shape was assumed to follow a LogParabola model (https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/source_models.html, accessed on 13 March 2020).

The data from the Swift and Fermi observatories were reduced by processing the files corresponding to both the photons and the position and the orientation of the satellite on Flexible Image Transport System-format files (FITS) with HEASoft version 6.26.1 (https://heasarc.gsfc.nasa.gov/docs/software/lheasoft/, accessed on 13 March 2020) and Fermitools version v11r5p3 (https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/, accessed on 13 March 2020) software. In both cases, it is necessary to know the relevant parameters of the object in question, such as the right ascension, declination, energy interval, and starting and ending times of the data processing (see, e.g., [37]). The satellite information was reviewed to make various time and position corrections, and with this, the analysis of the maximum likelihood estimation was carried out, in which the FITS files were used (see, e.g., [38]). Finally, a table was built that contains the light curve information, i.e., a file that contains discrete data of time, luminosity and its uncertainty with a signal-to-noise ratio <2$\sigma$.

The multifrequency light curves in Figure 1 are the result of processing the obtained data described above at a $3\sigma$ confidence level; thus, the accepted total fraction of bins is 99.7%.

Table 1. The table shows each electromagnetic band studied for Mrk 501: the number of observations, the duration in years, months and days (y, m, d), and the total days for which at least a measurement was obtained (d).

Band	Observations	Duration	Total
		(y, m, d)	(d)
Radio	615	$11,5,5$	4174
Optical	11,849	$23,2,17$	8441
X-rays	28,000	$12,7,9$	4607
$\gamma$-rays	4199	$11,7,9$	4239

shows the processed data in a synthesized way for all the electromagnetic frequencies studied for Mrk 501.

3. Methods

3.1. Periodograms and Window Functions

With the data processed, the search for periodicities was carried out. For this, the R-package RobPer version 1.2.3 was useful to analyze the associated periodograms on scales of a time of several years. This package was built for the analysis of time series in astrophysics, in which there is the possibility of having data that are not necessarily equidistant in time [39]. It is based on robustly fitting periodic functions of the time-series (light curves) in question and calculates periodograms with irregular (non-equidistant) observations on the time scale. The RobPer routines report significant peaks on the periodograms based on different statistical techniques (see [39] and references therein), with regression options such as least squares, least absolute deviations, least trimmed squares, and M-, S- and $\tau$-regression while accounting for measurement accuracies with weights. RobPer covers previous approaches and introduces model-regression combinations. Instead of relying on fixed critical values, it employs an outlier search on the periodogram to identify valid periods, considering that the traditional assumptions might not hold. RobPer, being a robust technique, is less sensitive to drastic changes in small portions of the data, making it more robust to outliers.

To support the RobPer periodogram analysis, a Lomb–Scargle (L-S) periodogram was also used, applying the lomb statistical routine in R and the LombScargle class of the astropy.timeseries routine in Python. This algorithm was proposed by Lomb [40] and Scargle [41], and to date, it is the most-used method for detecting periodicities in unevenly sampled light curves. The L-S technique is based on Fourier methods, phase-folding methods, least squares methods and Bayesian approaches [42]. Despite the success of the L-S routine, it is well known that it fails in some cases where eclipses are to be detected as periodicities (see, e.g., [43]). This means that one has to proceed with care when studying the L-S periodogram. In any case, not a single statistical method presented in this article can, in general, serve as a unique tool for the detection of periods on a time series. This is the reason why we studied the problem using different statistical techniques.

To avoid false positives in the detected periods with the RobPer and L-S routines, we built a windowing program following the work of Dawson and Fabrycky [44] and applied it to all the electromagnetic bands studied.

This method effectively distinguishes between aliases and true frequencies. Aliases occur at $|f\pm f_s|$, where f represents the frequency, and $f_s$ is a feature in the window function. By comparing the phase and amplitude of predicted aliases, derived from a sinusoidal waveform with the candidate frequency and the actual data, it is possible to determine whether the predicted aliases align with the data. When all aliases match in terms of amplitude, phase and pattern, we can confidently conclude that the true period has been found.

Thus, the spectral window function of an evenly sampled time series is

$$W\left(\nu \right)=\frac{1}{N}\sum _{r=1}^N{e}^{-2\pi i\nu t_r}$$

(1)

with N as the number of data points, m is an integer and $t_r$ is the time of the rth data point. The peak in the spectral window function occurs at $m{f}_s$.

Under this method, false positives are ruled out via a direct comparison of a $W\left(\nu \right)$ vs. $\nu$ diagram and the periodograms obtained with RobPer and L-S. In this way, if there is a false positive periodicity in the data collection (such as, e.g., an annual period of time due to vacation days in which no observations were taken), it will be presented as a peak in the window function (WF). The peaks that result as true positives in a periodogram are, thus, compared with their respective WF troughs in all the electromagnetic bands analyzed. Figure 2 shows the periodograms obtained with RobPer and L-S with the corresponding time, WF. The true mean peak in each periodogram is shown with a vertical dotted line. These mean periodicities are reinforced by the fact that the true mean peaks do not coincide with any peak in the WF. Figure 2 only shows a zoom into the relevant part of the periodogram, where a peak on it approximately coincides in all studied frequencies.

Figure 2. RobPer (magenta) and L-S (teal) periodograms (represented by their normalized coefficient of determination -NCoD), together with their corresponding window function (blue) for radio-, optical-, X- and $\gamma$-ray observations of Mrk 501 are shown from left to right, top to bottom panels. The black dotted vertical line shows the mean periodicity between the RobPer and L-S peaks that are common in all frequencies (radio: $228.03\mathrm{d}$; optical: $226.77\mathrm{d}$; X-rays: $223.20\mathrm{d}$; and $\gamma$-rays: $238.90\mathrm{d}$).

Table 2. The table shows the time in days (d) for the peaks and their mean values for the RobPer and Lomb–Scargle algorithms in each electromagnetic band studied for Mrk 501.

Band	RobPer	L-S	Mean
	(d)	(d)	(d)
Radio	228	$228.06$	$228.03$
Optical	228	$225.54$	$226.77$
X-rays	223	$223.40$	$223.20$
$\gamma$-rays	237	$240.80$	$238.90$

shows the periods obtained for the periodograms presented in the figure. The fact that these mean periodicities do not coincide with peaks in the window function reinforces their true periodic character.

Figure 2. RobPer (magenta) and L-S (teal) periodograms (represented by their normalized coefficient of determination -NCoD), together with their corresponding window function (blue) for radio-, optical-, X- and $\gamma$-ray observations of Mrk 501 are shown from left to right, top to bottom panels. The black dotted vertical line shows the mean periodicity between the RobPer and L-S peaks that are common in all frequencies (radio: $228.03\mathrm{d}$; optical: $226.77\mathrm{d}$; X-rays: $223.20\mathrm{d}$; and $\gamma$-rays: $238.90\mathrm{d}$). Table 2 shows the periods obtained for the periodograms presented in the figure. The fact that these mean periodicities do not coincide with peaks in the window function reinforces their true periodic character.

3.2. VARTOOLS

VARTOOLS (The VARTOOLS free GNU General Public License software is available at: https://www.astro.princeton.edu/~jhartman/vartools.html, accessed on 6 August 2023) is an open-source command-line utility for analyzing light curves developed by Hartman and Bakos [45]. It is written in the C programming language, and it provides a set of tools for processing, manipulating and studying light curves. The routines implemented in VARTOOLS for studying Mrk 501 were the L-S periodogram, the box-least squares (bls), the analysis of variance (aov) and the discrete Fourier transform (DFT).

AoV is a VARTOOLS routine for the detection of sharp (or statistically significant) periodic signals based on the code developed by Devor [46]. This method consists of folding and binning data with a trial period. It has also been used by Schwarzenberg-Czerny [47] for identifying and studying eclipsing binaries within large data sets of light curves. This routine customizes the process with a minimum period (minp option) and a maximum period (maxp option) to define the search range. The initial search uses a frequency resolution of the subsample, and the refined peak periods utilize a finetuned resolution. The program outputs the highest peaks in the periodogram.

The AoV-h VARTOOLS routine consists of the analysis variance using a multi-harmonic model. This method uses periodic orthogonal polynomials to fit the observations and the analysis of the statistical variance to evaluate the quality of the fit [48]. The parameters and behavior are as for the AOV routine.

The plots of the AoV and AoV-h of Mrk 501 are shown in Figure 3, which result in an average periodicity of $227.55\mathrm{d}$ for the AoV and $224.64\mathrm{d}$ for the AoV-h. The gray vertical zone represents the statistical range of these peaks.

The VARTOOLS BLS routine is commonly used in studies of stellar photometric time series in the search for periodic transits of exoplanets [49]. This algorithm uses the minimum and maximum values (the fraction of orbit in transit), and the search for periodicities operates within specified minimum and maximum periods (minper and maxper commands, respectively) and uses a specified number of trial frequencies. Additionally, the command allows the number of peaks to be reported. Using this routine for Mrk 501, we found a mean periodicity in all frequencies of $229.484\mathrm{d}$—see Figure 4.

The VARTOOLS Discrete Fourier Transform (DFT) algorithm calculates the power spectrum of the time series [50]. The routine to compute the DFT was developed by Kurtz [51]. The computation of the light curves is sampled using points per frequency and set to the maximum frequency (with maxfreq option), with a default value based on the minimum time separation. The DFT also allows researchers to find the highest peaks in the clean power spectrum. The average periodicity value obtained with this tool is $229.95\mathrm{d}$ for Mrk 501 in all studied frequencies—see Figure 4.

Table 3. The table shows the mean periodicity value of the AoV, AoV-h, BLS and DFT obtained with VARTOOLS for Mrk 501.

Band	AoV	AoV-h	BLS	DFT
	(d)	(d)	(d)	(d)
Radio	$227.20$	$228.06$	$220.960$	$228.737$
Optical	$226.73$	$227.40$	$224.141$	$225.909$
X-rays	$227.10$	$223.40$	$220.960$	$222.374$
$\gamma$-rays	$229.20$	$219.70$	$240.404$	$242.818$
Averages	$227.55$	$224.64$	$226.616$	$229.959$

summarizes the the periodicity results of VARTOOLS for Mrk 501.

3.3. Power Spectrum Density, Detrended Fluctuation Analysis, and the Colors of Noise

The study of time series or light curves, $f\left(t\right)$, in the study of this article, with sample length n, can be analyzed using the Fourier transform $P\left(\nu \right)$ given via (see, e.g., [52]):

$$P\left(\nu \right)=\frac{1}{n}\sum _{t=0}^{n-1}f\left(t\right)e^{(-\frac{2\pi i\nu t}{n})}.$$

(2)

With this transformation, it is possible to obtain specific periodic pulses in the light curve, even with excessive noise (see, e.g., [53]). If there are multiple periodic fluctuations present in the light curve, then the power spectral density (PSD) will reveal each one with a peak in the power spectrum at the particular periodicity or frequency.

The PSD function approaches a power law at a particular frequency, $\nu$, given via the following:

$$P\left(\nu \right)={\nu }^{-\alpha },$$

(3)

for a fixed exponent, $\alpha$. The spectral noise of the signal can be represented with its color of noise [53], depending on the value of $\alpha$. For astrophysical phenomena, it is commonly denoted as white, pink and Brownian according to the values of

Table 4. The table shows the intervals of the power spectrum density (PSD) exponent $\alpha$ of Equation (3) and its associated color of noise.

PSD Exponent	Colors of Noise
$0.0\lesssim \alpha \lesssim 0.5$	white
$0.5\lesssim \alpha \lesssim 1.5$	pink
$1.5\lesssim \alpha \lesssim 2.5$	Brownian

(see, e.g., [54] and references therein).

A given color of noise has statistical and correlation characteristics. When white noise predominates in the signal, it means that there is no temporal correlation in a specific time series, and the case of Brownian noise is obtained when a temporal correlation in the signal is significant [55]. Pink noise corresponds to the statistical case of phase change, in which there is a transition from a random process to a predictive one [56]. The colors of noise analysis has been performed in the literature for time series of solar bursts, magnetospheric physics, binary black holes and other astrophysical phenomena [54].

To support the PSD results, a detrended fluctuation analysis (DFA) method [57] with the time series can be computed to determine the statistical self-affinity in each frequency. DFA is also useful in the analysis of time series that appear to show long-memory processes. It is useful in the study of chaos theory, stochastic processes and time series analysis. The fluctuation, $F\left(n\right)$, is calculated for different window sizes, n, as follows:

$$F\left(n\right)=\sqrt{\frac{1}{N}\sum _{t=1}^N{(X_t-Y_t)}^2},$$

(4)

where N is time series length, $X_t$ is the cumulative sum or profile of the time series and $Y_t$ is the resulting piecewise sequence of straight-line fits (see, e.g., [58]). The resulting $logF\left(n\right)$ vs. $logn$ measures the statistical self-affinity of the time series, expressed through a fixed $\beta$ exponent from the following relation:

$$F\left(n\right)\propto n^{\beta }.$$

(5)

Table 5. The table shows the obtained values for the exponents $\alpha$ and $\beta$ resulting from the PSD and DFA color of noise analysis to the multi-frequency observations of Mrk 501. In all cases, the resulting color of noise is pink, according to the classification presented in Table 4, which implies that the light curves present temporal correlations with random fluctuations.

Band	$\mathit{\alpha }$	$\mathit{\beta }$
Radio	$1.24$ ± $0.020$	$1.287$ ± $0.0035$
Optical	$1.16$ ± $0.030$	$0.869$ ± $0.0045$
X-rays	$0.61$ ± $0.005$	$1.019$ ± $0.0045$
$\gamma$-rays	$0.87$ ± $0.035$	$0.853$ ± $0.0015$

shows that the application of a PSD and a DFA to the multi-frequency light curves of Mrk 501 results in pink noise. For completeness, Figure 5 shows the PSD plots in all studied wavelengths. The PSD and DFA analyses were performed using VARTOOLS, and the fitting of the corresponding $\alpha$ and $\beta$ exponents were calculated using our own C program with the aid of the GNU Scientific Library (https://www.gnu.org/software/gsl, accessed on 18 January 2024) routines.

The power spectrum density and the detrended fluctuation analysis complement the methods used to search for periodicities. Nevertheless, the variations in the emission of blazars are more complex. For example, Vaughan et al. [59] mention that periodic variations in quasar light curves have raised the possibility of close binary supermassive black holes. However, quasars typically display stochastic variability, making it challenging to identify true periodic candidates. Methods like Bayesian analysis light curves favors a stochastic process over sinusoidal variation, though simulations of the quasar PG 1302-102 light curves demonstrate the occurrence of periodicity in stochastic processes, which is crucial to calibrate false positive rates when searching for periodic signals among numerous targets. The studied X-ray observations of PG 1211+143 [60] showed time lags where strong correlations were observed, including complex variability patterns. This phenomenon was studied as well by Vaughan [61], analyzing X-ray variability in black hole sources including AGN, X-ray binaries and UltraLuminous X-ray objects, finding the accretion flow as the source of the observed X-ray variations.

3.4. Data Fitting via a Jacobi Elliptical Function

The achromatic periodicity, from radio to $\gamma$-rays, of ≲$229\mathrm{d}$ detected in the previous sections, combined with the pink noise color, implies that a correlated phenomenon is producing it. The simplest way in which such a phenomenon may occur is with an eclipse, as shown in Figure 6, since an orbiting massive object partially covering the radiation from the central parts of Mrk 501 will do so in a periodic way. In what follows, we denote the central supermassive black hole as the primary black hole and the eclipsing massive object as the secondary object. When light is obscured due to a non-relativistic massive object, it gets dimmed, as happens in a standard eclipse of, say, an exoplanet when it passes through the central star of its planetary system [62]. However, light is magnified if a relativistic object passes through light beams produced by a primary source of light (see, e.g., [63] and references therein). The following eclipse function, $e\left(t\right)$, serves as a good empirical way to model an eclipse:

$$e\left(t\right)=\pm \left\{\Theta \left(t\right)-\Theta (t-\pi )\right\}{\mathrm{sn}}^2\left(\frac{2K\left(m\right)t}{\pi }\right)$$

(6)

where $\Theta$ represents the heaviside step function, $\mathrm{sn}$ is the elliptic sine or sinus amplitudinis Jacobi function with module m, such that $0\le m\le 1$, and $K\left(m\right)$ is the complete elliptic integral of the first kind given via (see, e.g., [64]):

$$K\left(m\right)={\int }_0^{\pi /2}\frac{\mathrm{d}\zeta }{\sqrt{1-m{sin}^2\zeta }}.$$

(7)

The ± sign on the right-hand side of Equation (6) represents a standard non-relativistic eclipse, i.e., an eclipse that diminishes the amount of received radiation, for the minus sign and a magnification relativistic eclipse for the plus sign. In the limit that occurs when the module $m\to 1$, a single squared pulse is obtained, and when $m\to 0$ occurs, a sinusoidal curved profile is obtained.

To model the obtained periodicity on the light curves of Mrk 501 shown in the previous sections, we built a C program capable of dealing with the expectations of an eclipse fingerprint produced on a long-term light curve based on a modification of Equation (6). The program has as free parameters, the amplitude, A, of the eclipse, the duration time, $t_{\mathrm{e}}$, of the eclipse, the quiescent duration time, $t_{\mathrm{q}}$, when the eclipse is not occurring and the number of occurrences of the eclipse, n. We chose the statistical average $\langle f\rangle$ as the zero or quiescent flux of the light curve. A set of artificial examples of these eclipses is presented in Figure 7.

The data fitting for Mrk 501 was performed using a Monte Carlo method, based on the program mentioned in the previous paragraph. Random seeds were used for the unknown parameters of the program, using ${10}^6$ attempts for each light curve for Gaussian noise.

Table 6. The table shows the best fit of the eclipse model using a Monte Carlo method for the multiwavelength data of Mrk 501. The columns represent electromagnetic wavebands, the eclipse amplitude, A, the elliptic Jacobi function module, m, the duration time of the eclipse, $t_{\mathrm{e}}$, the interval of time where the eclipse is not occurring, i.e., the quiescent time, $t_{\mathrm{q}}$, the number of times, n, the eclipse occurred, the total periodicity, $p:=t_{\mathrm{e}}+t_{\mathrm{q}}$, the mean value of the flux, $\langle f\rangle$, and the dimensionless brightness magnification produced via the eclipse $A/\langle f\rangle$.

Band	A	m	${\mathit{t}}_{\mathbf{e}}$	${\mathit{t}}_{\mathbf{q}}$	n	p	$\langle \mathit{f}\rangle$	$\mathit{A}/\langle \mathit{f}\rangle$
			Days	Days		Days
Radio	$0.09 \pm 0.04$	$0.85 \pm 0.2$	$88.33 \pm 2.83$	130.20 ± 12.77	$16.43 \pm 1.80$	$218.53 \pm 15.6$	$1.11 \pm 0.022$	$0.08 \pm 0.04$
Optical	$1.18 \pm 0.02$	$0.14 \pm 0.21$	$89.76 \pm 3.60$	$133.18 \pm 0.09$	$35.03 \pm 0.43$	$222.94 \pm 3.69$	$13.48 \pm 1.34$	$0.08 \pm 0.01$
X-rays	$0.99 \pm 0.44$	$0.82 \pm 0.02$	$89.37 \pm 0.27$	$137.05 \pm 3.61$	$19.96 \pm 0.51$	$226.43 \pm 0.41$	$2.46 \pm 0.68$	$0.40 \pm 0.29$
$\gamma$-rays	$0.99 \pm 0.07$	$0.20 \pm 0.17$	$86.25 \pm 0.32$	$142.10 \pm 4.45$	$17.69 \pm 0.41$	$228.36 \pm 4.77$	$0.13 \pm 0.48$	$7.30 \pm 26.65$

shows the results of the fit for each waveband.

4. Discussion

A different analysis of the long-term multiwavelength (from radio to $\gamma$-rays) light curves of Mrk 501 show a common achromatic periodicity of ≲$229\mathrm{d}$. The Monte Carlo fitting technique, assuming a relativistic eclipse as the cause of this periodicity, shows a periodicity of ∼$224\mathrm{d}$. The eclipse is produced by a secondary supermassive black hole orbiting the primary supermassive black hole, i.e., about the central engine of Mrk 501. The reasoning for this conclusion is summarized in the following paragraphs.

The periodicities found with the RobPer and L-S algorithms described in Section 3.1 are all consistent in all frequencies of Mrk 501 with average values in radio, optical, X-rays and $\gamma$-rays, respectively, given via $228.03$, $226.77$, $223.20$ and $238.90$ days, according to the results of Table 2. The mean value for these periodicities is $229.225\mathrm{d}$.

With the use of the VARTOOLS software described in Section 3.2, we found that, for the AoV, AoV-h, BLS and DFT routines, the periodicity value lies in the intervals 226.73–229.2, 219.7–228, 220.960–240.404 and 222.374–242.818 days, respectively, according to the results presented in Figure 3 and Figure 4.

The fact that the color of the signal noise in the light curves of Mrk 501 presented in Section 3.3 is pink means that, for each particular waveband, there is a robust oscillation (Brownian noise color) with a periodicity accompanied by a random signal (white noise color).

Due to the achromatic nature of the found periodicity of ≲$229\mathrm{d}$, we modeled this periodicity as a relativistic eclipse caused by an orbiting supermassive black hole about the central engine of Mrk 501. The results of Table 6 show that this model is quite coherent in the fitting of the long-term multi-frequency light curves. The only small inconsistency found is with the dimensionless brightness magnification $A/\langle f\rangle$ presented in X-rays and more prominent in $\gamma$-rays. This is most probably due to the large errors reported in $\gamma$-rays and the large gaps that appear in the X-ray light curve. Figure 8 shows a time-folding of the multi-frequency light curves with the corresponding eclipse function using the results of Table 6. The shaded region represents the time duration, $t_{\mathrm{e}}$, of the eclipse. The dotted horizontal lines represent a $1\sigma$ significance level. It is important to note that the radio, optical and X-ray panels in Figure 8 do not contain all the used data points at a $3\sigma$ confidence level. In other words, they represent zooming in to the light curves in order to emphasize the detected eclipse. The large error bars in the $\gamma$-ray make the bump on the eclipse function cross the $3\sigma$ confidence level.

As mentioned in the introduction, shockwave interactions and motions inside the jet could cause periodic oscillations on the light curves, and although the periodicity in this case may appear, in principle, as an achromatic effect, it is extremely unlikely that this is the case for Mrk 501 since the reported amplitudes, A, in Table 6 are quite similar to one another, a property unlikely to occur in periodic shock wave formation and interactions due to the complexity of the radiation produced at each particular frequency.

In this article, we found an achromatic periodicity of ∼224 d in the radio, optical, X- and $\gamma$-rays light curves of Mrk 501, which we modeled as an eclipsing event produced by a massive (secondary) supermassive black hole orbiting the (primary) supermassive black hole of Mrk 501. Since the eclipse is produced by a relativistic object (the secondary black hole), the radiation brightness is magnified (contrary to what occurs in, e.g., a solar eclipse on Earth).

Using the results of our eclipse model in Table 6 and Kepler’s third law for a circular orbit given via

$$v=\frac{r_{\mathrm{orbit}}}{p}=\sqrt{\frac{G{M}_{\mathrm{prim}}}{r_{\mathrm{orbit}}}},$$

(8)

where v is the velocity of the orbiting test mass (secondary black hole), $r_{\mathrm{orbit}}$ and p are the radius and period of the orbit, $M_{\mathrm{prim}}$ is the mass of the central supermassive black hole and G represents Newton’s constant of gravitation, the following conclusions can be drawn about the results obtained in this article:

• The mass of the central (primary) black hole in Mrk 501 is $M_{\mathrm{prim}}\sim$${10}^9{\mathrm{M}}_{\odot }$ [65], which means that its gravitational radius is $r_{\mathrm{g}-\mathrm{prim}}\sim$$20\mathrm{au}$.
• The radius of the orbit of the eclipsing (secondary) binary black hole is $r_{\mathrm{orbit}}\sim$$200\mathrm{au}\sim 10r_{\mathrm{g}-\mathrm{prim}}$.
  Using a full relativistic approach for Schwarzschild’s space–time, Equation (8) can be written as follows [66]:

  $$v^2=\frac{1}{2}{c}^2\left(r_{\mathrm{g}-\mathrm{prim}}/r_{\mathrm{orbit}}\right)\left(1-3r_{\mathrm{g}-\mathrm{prim}}/2r_{\mathrm{orbit}}\right),$$

  (9)
  The above equation is cubic for the radius $r_{\mathrm{orbit}}$, for which the only real solution is given via the following:

  $$r_{\mathrm{g}-\mathrm{prim}}/r={\left(A\sqrt{1+A}\right)}^{1/3}-A/{\left(A\sqrt{1+A}\right)}^{1/3},$$

  (10)

  where

  $$A:=r_{\mathrm{g}-\mathrm{prim}}^2/c^2{p}^2.$$

  (11)
  By using the same input values as the ones we used for the non-relativistic Kepler formula, we obtain $r_{\mathrm{orbit}}\sim$$5r_{\mathrm{g}-\mathrm{prim}}$, which is of the same order of magnitude as the $10r_{\mathrm{g}-\mathrm{prim}}$ that we obtained with the simple Kepler relation.
• The orbital period of the eclipsing binary black hole is ∼$224\mathrm{d}$.
• The orbital velocity of the eclipsing binary black hole is ∼$0.3\%$ of the speed of light, i.e., ∼$3\times {10}^6\mathrm{km}/\mathrm{h}$.
• The brightness magnification of the radiation produced by an eclipse due to the secondary black hole is ≳$10\%$.
• The coalescence time of the binary system is approximately given via the following [67]:

  $$t_{\mathrm{coal}}\approx \frac{1}{50}\frac{c^5}{G^3}\frac{r_{\mathrm{orbit}}}{M_{\mathrm{prim}}^2{M}_{\sec }}\sim \frac{{10}^8\mathrm{yrs}}{\left(M_{\sec }/M_{\odot }\right)},$$

  (12)

  where $M_{\sec }(\ll M_{\mathrm{prim}})$ is the mass of the secondary or eclipsing black hole. Since the periodicity found corresponds to a few decades of observation—meaning that the orbit is sufficiently stable—it is reasonable to assume that $t_{\mathrm{coal}}\gtrsim {10}^3$, which means that the mass of the secondary black hole $M_{\sec }\lesssim {10}^5{\mathrm{M}}_{\odot }$.
• The results of (9) imply that the orbit of the secondary black hole is half the value obtained from the Newtonian calculation, and so its orbital frequency is about $2^{3/2}=2.83$∼3 its Newtonian value. Since the frequency, f, of the of the binary’s quadrupolar gravitational waves is about a factor of 2 larger than the orbital frequency (see, e.g., [68]), then

  $$f\sim 6\sqrt{GM/r_{\mathrm{orbit}}^3},$$

  (13)

  where the total mass of the binary system is given via $M=M_{\mathrm{prim}}+M_{\sec }\approx M_{\mathrm{prim}}$. In other words,

  $$f\sim 6\times {10}^5\left(\frac{M_{\odot }}{M_{\mathrm{prim}}}\right)\mathrm{Hz},$$

  (14)

  which, for a value of $M_{\mathrm{prim}}$∼${10}^9{M}_{\odot }$, yields f∼$6\times {10}^{-4}\mathrm{Hz}$. This suggests that the frequency of the emitted gravitational waves is just above $0.1\mathrm{mHz}$, which is within the range of the forthcoming LISA space-based interferomenter, and as such, the prediction of a binary black hole system in Mrk 501 should be considered a science target for future LISA observations.

The magnification $\mu$ of the radiation is the ratio of the Einstein radius, ${\theta }_{\mathrm{E}}$, to the distance, $\beta$, from the line of sight connecting the observer to the lens [69]. Since the secondary eclipsing black hole mass ∼${10}^5{M}_{\odot }$, then ${\theta }_{\mathrm{E}}$∼${10}^{-8}$, and so, in order to get a magnification of $\mu =1.1$, then $\beta ⪆{10}^{-8}$. In other words, it would seem that the primary and secondary black holes should be quite well aligned with our line of sight. However, since Mrk 501 is a blazar, we are observing the source within an angle of no more than $30°$ from the emitted jet (in fact, for Mrk 501, the observing angle is between $15°$ and $25°$, as reported by Giroletti et al. [70]), and so the radiation would be a combination of the radiation processes occurring close to the central engine plus the relativistic jet. This extended emission from the relativistic jet means that the chances of having an 10% amplification increase.