Follow-up on the Supermassive Black Hole Binary Candidate J1048+7143: Successful Prediction of the Next Gamma-Ray Flare and Refined Binary Parameters in the Framework of the Jet Precession Model

In recent years, multimessenger astronomy has become one of the most rapidly evolving fields in astronomy. Being the synergy of coordinated observations targeting all of the four extragalactic messengers, such as the electromagnetic (EM) radiation, cosmic rays, neutrinos, and gravitational waves (GWs), multimessenger astronomy is a very useful tool for studying the most energetic phenomena of the Universe, such as the merging of supermassive black holes (SMBHs), see, e.g., Becker (2008).

In 2011, the IceCube Collaboration detected the high-energy neutrino flux with cosmic origin for the first time ever (IceCube Collaboration 2013). In 2016, the LIGO Scientific Collaboration and Virgo Collaboration announced the twofold direct detection of gravitational waves (GWs) in both Advanced LIGO observatories in 2015 (Abbott et al. 2016). In 2017, the GW signal of two colliding neutron stars, as well as a short γ-ray burst signal have been also detected (Abbott et al. 2017a, 2017b). Patterns of high-energy particle sources are thoroughly discussed in the literature (Aartsen et al. 2020; Franckowiak et al. 2020; Kun et al. 2022a; Novikova et al. 2023); however, there are only two direct neutrino-source associations to date (IceCube Collaboration et al. 2018, 2022). The simultaneous observation of GWs and high-energy neutrinos still lies ahead.

Since neutrinos interact only weakly with matter, they are able to pinpoint cosmic high-energy particle accelerators in the sky that would be otherwise hidden because ultra-high-energy cosmic rays (UHECRs) scatter due to the Galactic and intergalactic fields (Dermer et al. 2009), while their energy is strongly attenuated beyond the Greisen–Zatsepin–Kuzmin horizon (e.g., Stanev et al. 2000). High-energy γ-photons also lose their energy on their way to Earth in interactions with the photons of the cosmic microwave background. The cosmic neutrinos are indeed the key messengers to explore the Universe where it is opaque to the cosmic rays and photons. Violent phenomena in the distant Universe, for example merging SMBHs in the young Universe, can be observed through the observation window of GWs and/or neutrinos only. Based on the potential detection of the GW background by the NANOGrav Collaboration (Agazie et al. 2023), supermassive binary black hole (SMBBH) mergers should happen more often than previously assumed. Interpreting the X-shaped jets of radio galaxies as relics of jet precession, which is a plausible assumption, e.g., based on their mass excess (see Gopal-Krishna et al. 2012), one can argue that there are more such mergers potentially observed. The frequency at which the SMBH mergers happen is not quite clear though, and potential periodic neutrino sources might help understand the mergers and their rate better.

In Kun et al. (2022b), we analyzed and modeled the Fermi Large Area Telescope (LAT) light curve of the blazar J1048+7143, and concluded that the multiwavelength behavior of this source is compatible with the spin–orbit precession of an SMBBH (Gergely & Biermann 2009) in the heart of the galaxy. This active galactic nucleus (AGN) indeed belongs to a group of blazars that reveal precession-induced variability (see Britzen et al. 2023, and references therein). Generalizing the model of de Bruijn et al. (2020), we predicted the GW signal of this blazar. We also predicted the appearance time of a fourth γ-ray flare. This model was already successfully applied to the multiepoch neutrino observations from the direction of TXS 0506+056 (IceCube Collaboration et al. 2018; IceCube Collaboration 2018) by Becker Tjus et al. (2022). Britzen et al. (2019) also proposed an SMBBH for the same source. de Bruijn et al. (2020) predicted a neutrino episode for TXS 0506+056, which was later confirmed by the detection of a high-energy neutrino by the Antarctic IceCube Neutrino Detector from the approximate direction of this blazar. We note that this neutrino only skimmed the edge of the detector, and therefore its uncertainty (90% containment) is a relatively large ∼358 (the neutrino arrived with the energy of ∼170 TeV and the event has a signalness of 42%; Blaufuss et al. 2022). In this Letter, we extend the γ-ray light curve of J1048+743 by analyzing the new Fermi-LAT data, and indeed find the predicted fourth flare. Using the extended γ-ray light curve, we fine-tune the binary parameters that are consistent with the observations within our model.

We obtained archival data taken with the LAT instrument onboard the Fermi Gamma-ray Space Telescope to measure the γ-ray flux of 4FGL J1048.4+7143 (associated with J1048+7143; Abdollahi et al. 2020), located at the sky coordinates of R. A._J2000 = 1621067 and decl._J2000 = 717297. The region of interest with a 15° radius was centered at the J2000 sky coordinates of 4FGL J1048.4+7143. We extended our previous analysis published in Kun et al. (2022b) by adding Fermi-LAT data in the time range from 2022 March 14 to 2023 June 4 (MJD 59635–60099) in the energy range 100 MeV–800 GeV. Now our analysis covers almost 15 yr between 2008 August 4 and 2023 June 4.

We repeated the same binned likelihood analysis of Fermi-LAT data for the time range MJD 59635–60099 as we did earlier (Kun et al. 2022b). Technical details of the analysis chain can be found in that work. The binning of the light curve and the optimum energy were constrained using adaptive binning with a 5% threshold of relative flux error. In the resulting light curve, shown in Figure 1, we see a new major flare starting at the beginning of 2021 (after MJD 59200), in addition to the previously reported three flares (Kun et al. 2022b). This fourth major flare also shows a double subflare structure similar to the three major flares already reported in Kun et al. (2022b). We refitted the light curve with two-sided exponential functions in the form of

$$\begin{eqnarray}&&{F}_{i,j}(t)={a}_{0}+{a}_{i,j}\times \exp \left[\displaystyle \frac{-| t-{b}_{i,j}| }{{c}_{i,j}}\right],\end{eqnarray} \tag{ 1 }$$

where t is time, the index i runs from 1 to 4 for the four main flares, j is 1 or 2 for the subflares of the ith main flare, a_i,j the height, b_i,j the time location of the peak, and c_i,j the slope of the exponential function. Since in the quiescent phase the γ-ray flux did not vanish, we also fitted a constant baseline (represented by a₀). We present the resulting fit parameters in Table 1. We derived the center time of each of the four major flares in the extended light curve using the centroid method described in Kun et al. (2022b). We present the flare characteristics using the centroid method in Table 2 and show the centroids in Figure 1 as the cross of the green dashed lines. The intersection of the horizontal green dashed line with the exponential fit (blue continuous line) of the respective main flare indicates the start and end of the flare. The duration is thus seen as the gray shaded area and given with error bars in Table 1. Its uncertainty is determined using the Markov Chain Monte Carlo method using 10,000 iterations for each flare. In addition, the vertical green dashed line indicates the date of the main flare center (also listed in Table 2 with its uncertainty).

Zoom In Zoom Out Reset image size

The periods between the flares are computed as the time elapsed between these main flare centers. As a result, P_1→2 = (3.20 ± 0.13) yr is obtained for the time between main flare one and two, P_2→3 = (3.07 ± 0.09) yr between two and three, and P_3→4 = (3.06 ± 0.17) yr between three and four (see Table 2).

To derive the structural and kinematic properties of the milliarcsecond-scale jet of J1048+7143, we used archival calibrated 8.6 GHz very long baseline interferometry (VLBI) visibility data taken with the Very Long Baseline Array (VLBA). The observations span more than 26 yr at 69 epochs (between 1994.61 and 2020.73), and the calibrated visibilities are available in the Astrogeo database. ¹² We introduced these observations already in Kun et al. (2022b).

We model the brightness distribution of the source visibility data by fitting elliptical Gaussian components (Pearson 1995) to the core and circular Gaussians to the jet components using the Difmap (Shepherd et al. 1994) software. The fitted parameters are the component flux density, the position, the component width, and for only the core the position angle of the major axis (measured from north through east) and the minor-to-major axial ratio of this elliptical component. We followed Kun et al. (2014) for the error estimation of the fitted parameters.

For VLBI jets, seen in small viewing or inclination angle with respect to the line of sight, the apparent brightness temperature T_b usually exceeds the limiting intrinsic brightness temperature T_int due to strong Doppler boosting. The Doppler factor δ connects them as T_b = δ T_int. The brightness temperature of the VLBI components is estimated as (e.g., Condon et al. 1982)

$$\begin{eqnarray}&&{T}_{{\rm{b}},\mathrm{VLBI}}=1.22\cdot {10}^{12}\times (1+z)\displaystyle \frac{{S}_{\nu }}{{d}_{1}{d}_{2}{f}^{2}}\text{}({\rm{K}}),\end{eqnarray} \tag{ 2 }$$

where S_ν is the flux density (in Jy), d₁ and d₂ are the major and minor axes of the core (in milliarcseconds), respectively, f is the observing frequency (in GHz), and z = 1.15 (Polatidis et al. 1995) is the redshift of the source. Assuming that the intrinsic brightness temperature is slightly lower than the equipartition brightness temperature, T_int ≈ 3 × 10¹⁰ K (Homan et al. 2006), we obtained the median value of δ ≈ 5.1 for the core at 8.6 GHz.

We were able to identify two moving components, C1 and C2. We calculated the apparent proper motion of the components first by fitting linear proper motions as

$$\begin{eqnarray}&&{r}_{{\rm{L}}}=\mu (t-{t}_{\mathrm{ej}}),\end{eqnarray} \tag{ 3 }$$

and then by fitting accelerating proper motions as

$$\begin{eqnarray}&&{r}_{{\rm{A}}}={r}_{{\rm{L}}}+\displaystyle \frac{\dot{\mu }}{2}{\left(t-{t}_{\mathrm{mid}}\right)}^{2},\end{eqnarray} \tag{ 4 }$$

where μ is the linear proper motion measured in mas yr⁻¹, t_ej the ejection time, $\dot{\mu }$ the angular acceleration, and t_mid the half of the sum of the maximum and the minimum epochs when the respective components are detected.

We fitted the core separation of C1 and C2 as a function of time before the flaring behavior of J1048+7143 started according to its single-dish radio flux density curve (see in Kun et al. 2022b). We show the proper motion fits in Figure 2. The resulting best-fit proper motions, as well as the derived apparent speeds are shown in Table 3. As can be seen from the reduced χ² values in Table 3 and from Figure 2, the accelerated proper motion fits give better results.

Zoom In Zoom Out Reset image size

Table 3. Component Speed Fits

ID	μL	tej,L	χ2	μA	${\dot{\mu }}_{{\rm{A}}}$	tej,A	χ2
	(mas yr−1)	(yr)		(mas yr−1)	(mas yr−2)	(yr)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)
C1	0.118 ± 0.006	1996.160 ± 0.250	3.503	0.113 ± 0.002	−0.033 ± 0.004	1995.676 ± 0.139	0.578
C2	0.074 ± 0.016	1993.811 ± 3.182	5.034	0.073 ± 0.014	0.020 ± 0.011	1994.565 ± 2.606	3.006

Note. (1) jet component identifier. From linear fit (L): (2) linear proper motion, (3) ejection time, (4) reduced χ². From accelerated fit (A): (5) nonlinear proper motion, (6) acceleration rate, (7) ejection time, (8) reduced χ².

Download table as:  ASCIITypeset image

The apparent speed of the jet components is related to the bulk jet speed β and their inclination angle ι as (e.g., Urry & Padovani 1995)

$$\begin{eqnarray}&&{\beta }_{\mathrm{app}}=\displaystyle \frac{\beta \sin \iota }{1-\beta \cos \iota },\end{eqnarray} \tag{ 5 }$$

where the jet speed β is related to the Lorentz factor Γ as

$$\begin{eqnarray}&&\beta =\sqrt{1-\displaystyle \frac{1}{{{\rm{\Gamma }}}^{2}}}.\end{eqnarray} \tag{ 6 }$$

The maximal β_app occurs for the critical inclination $\cos {\iota }_{c}=\beta$, $\sin {\iota }_{c}={{\rm{\Gamma }}}^{-1}$, with the value ${\beta }_{\mathrm{app},\max }={({{\rm{\Gamma }}}^{2}-1)}^{1/2}$.

In practice, it is unlikely that a jet component will exactly move in this direction; therefore, the maximal speed ${\beta }_{\mathrm{app},\max }^{\mathrm{obs}}\leqslant {\beta }_{\mathrm{app},\max }$ seen in the VLBI jet gives a lower limit on the Lorentz factor of the jet,

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{low}}=\sqrt{1+{\left({\beta }_{\mathrm{app},\max }^{\mathrm{obs}}\right)}^{2}}.\end{eqnarray} \tag{ 7 }$$

Taking the accelerated proper motion of C1, ${\mu }_{\max }=(0.113\pm 0.002)$ mas yr⁻¹, that is the fastest jet component seen in J1048+7143 before the flaring behavior starts, we get the maximum observed apparent speed as ${\beta }_{\mathrm{app},\max }^{\mathrm{obs}}=6.56\pm 0.12$. Then, the limiting Lorentz factor emerges as Γ_low = 6.87 and the minimum intrinsic jet speed is β_low = 0.989 (in the units of the speed of light). Then, the related critical inclination angle is ${\iota }_{{\rm{c}}}=\arcsin ({{\rm{\Gamma }}}_{\mathrm{low}}^{-1})\approx 10\buildrel{\circ}\over{.} 9$ and the half-opening angle of the jet is estimated as ζ ≈ 1/Γ ≈ 83 (Boettcher et al. 2012).

In Kun et al. (2022b) we have shown that the γ-ray light curve of J1048+7143 (E > 168 MeV) and its single-dish radio flux density curve obtained at 4.8 GHz (see bottom panel of Figure 1) are consistent with the jet precession driven by the spin–orbit precession, occurring in an SMBBH with a total mass of m = 10^9.16 M_⊙. We were able to constrain the mass ratio only from the above, such that q < 0.2, where q = m₂/m₁ (m = m₁ + m₂, m₁ > m₂). We note the ∼0.1 Jy difference in flux density seen between the single dish radio flux density curves in Figure 1 is most probably caused by the difference in the resolving power between the RATAN600 and Nanshan (Urumqi) radio telescopes at 4.8 GHz.

In the inspiral phase of the merger of an SMBBH, first the dynamical friction, then the gravitational radiation shrinks the orbit of the black holes. In the latter case, the binary goes through three phases, the inspiral, the plunge, and the ring-down. In the inspiral phase, the equation of motion can be approximated by a series expansion in terms of a small parameter, the post-Newtonian (PN) parameter $\varepsilon ={Gm}{({c}^{2}r)}^{-1}$, where G is the gravitational constant, c is the speed of the light, m is the total mass, and r is the separation of the binary.

The compact binary dynamics is conservative up to the second PN order, such that the total energy E and the total momentum J = S ₁ + S ₂ + L _N are conserved during the motion. If the S ₁ and S ₂ the spins of the two SMBHs (m₁ > m₂) are not parallel with the Newtonian orbital momentum L _N , the spins begin to precess (Barker & O'Connell 1975, 1979):

$$\begin{eqnarray}&&{\dot{{\boldsymbol{S}}}}_{i}={{\boldsymbol{\Omega }}}_{i}\times {{\boldsymbol{S}}}_{i},\end{eqnarray} \tag{ 8 }$$

where Ω_i is the angular momentum of the ith spin and it contains spin–orbit (1.5PN), spin–spin (2PN), and quadrupole–monopole contributions (2PN) up to 2PN order. The gravitational radiation circularizes the orbits (Peters 1964); therefore, from here we take only circular orbits into account.

Expanding the work of de Bruijn et al. (2020), in Kun et al. (2022b) we determined the direction angle of the dominant spin (ϕ) as a function until the final coalescence of the compact binary as

$$\begin{eqnarray}&&\begin{array}{l}\phi ({\rm{\Delta }}{T}_{\mathrm{GW}},q)=-\displaystyle \frac{2\,(4+3q)}{{\left(1+q\right)}^{2}}\times \\ \,\times {\left(\displaystyle \frac{5\,c}{32\,{G}^{1/3}{m}^{1/3}}\cdot \displaystyle \frac{{\left(1+q\right)}^{2}}{q}\right)}^{3/4}{\left({\rm{\Delta }}{T}_{\mathrm{GW}}\right)}^{1/4}+\psi ,\end{array}\end{eqnarray} \tag{ 9 }$$

where ψ is an integration constant.

Assuming the dominant spin makes a complete 360° turn between the subsequent major flares in the γ-ray light curve, we can establish the connection between two flares as

$$\begin{eqnarray}&&\phi ({\rm{\Delta }}{T}_{\mathrm{GW}},q)=\phi ({\rm{\Delta }}{T}_{\mathrm{GW}}-{P}_{\mathrm{jet}},q)\pm \zeta \,,\end{eqnarray} \tag{ 10 }$$

where P_jet is the jet precession period, which is the time that elapsed between two subsequent major flares in the γ-ray light curve.

Updating the light curve with the fourth major γ-ray flare, the mass ratio of the SMBBH is now constrained as 0.062 < q < 0.088 (see Figure 3). In our previous work, we were able to constrain the mass ratio as q < 1/4. This means that including the new flare, now we can constrain the mass ratio into a significantly narrower range, 0.062 < q < 0.088, and consequently we can further constrain the parameters of the hypothetical SMBBH at the heart of J1048+7143. Since the emission of GWs leads to shrinking orbits, the Newtonian orbital angular momentum L _N and the direction of the dominant spin S ₁ also change: the former is moving away, while the latter is moving closer to the total angular momentum J . Therefore, the jet, attached to S ₁, also shows a slow, secular change in its direction. We see the Doppler boosting governed by the spin–orbit precession while the jet lies close to our line of sight. Once the jet is far away from our line of sight because of the secular change in the direction of S ₁, we do not expect a flaring behavior anymore due to the weak Doppler boosting. Given the jet still boosted in the close future, a fifth γ-ray flare is expected to arrive between MJD 60379 (2024 March 10) and MJD 61350 (2026 November 6).

Zoom In Zoom Out Reset image size

Periodicities in AGN can usually be explained by a number of model families, relying on single (e.g., Lai 2003) or binary black holes (e.g., Roos et al. 1993; Bon et al. 2012; Kun et al. 2018; Jaiswal et al. 2019; Kun et al. 2020; Britzen et al. 2023; Kun et al. 2023), on the accretion disk (Bardeen & Petterson 1975; Caproni et al. 2006), or instabilities within the jet (e.g., Hardee et al. 1994; Perucho et al. 2006). Application of these models depends on the complexity of the observed behavior: while models relying on, e.g., single black holes or the precession of an accretion disk can usually explain a stable period, for cases when the time duration between flares changes, inclusions are needed. We note that in quasars, quasiperiodic oscillations can be attributed to stochastic changes in the accretion disk (e.g., Kuźmicz et al. 2022), or disk disturbances due to companion black holes (e.g., Caproni & Abraham 2004), etc. However, J1048+7143 belongs to the blazar subclass of radio-loud AGN as its jet is seen under a small viewing angle. For such sources, relativistic beaming dominates the electromagnetic emission as the boosted jet can overshine any intrinsic signal.

Here we offer a qualitative scenario to explain the multiwavelength behavior of J1048+7143 (see the γ-ray light curve and single-dish radio flux density curve in Figure 1). We note that any modeling targeting this source should explain the decreasing nature of the time duration between subsequent flares and their double structure in γ-rays. The scenario we found to be the most compatible with these is the spin–orbit precession of the supermassive black hole that launches a spine–sheath jet and drives the multiwavelength appearance of the source. Due to the precession, an (e⁻, e⁺) spine–(p, e⁻) sheath jet periodically goes through our line of sight. The precessing jet approaches the target (e.g., clouds in the broad-line region) close to our line of sight. Their interaction imprints three phases in the γ-ray and radio regimes. First, the "upper" sheath with accelerated cosmic rays meets the target (that gives the target protons), which leads to the first (pionic) γ-ray subflare. Meanwhile, the radio flux density goes up as relativistic boosting gets stronger due to the decreasing viewing angle (Phase I in Figure 4). Then, the spine meets the target leading to the minimum in the pionic γ-rays between the subflares since the proton density in the spine is negligible, and therefore there are no hadronic interactions between the jet and the target (Phase II in Figure 4). Meanwhile, the radio flux density curve reaches its maximum as we see the jet at the smallest viewing angle. This boosted radio emission is expected to come dominantly from the spine of the jet. Then, the "lower" sheath meets the target that leads to the second γ-ray subflare, while the radio flux density goes down as relativistic boosting gets weaker due to the increasing viewing angle of the jet (Phase III in Figure 4). Then, the jet leaves the target and moves more away from our line of sight. Since the VLBI jet is visible between the major flares, hinting that considerable boosting is still present, we can conclude that the jet moves not far from our line of sight. This is when the quiescence phase happens between the major flares. After one spin–orbit period, the jet approaches again, triggering the next sequence of I/II/III phases (Figure 4). We note that the optical light curve of J1048+7143, constructed based on visual optical observations sent to the American Association of Variable Star Observers, also shows the double flaring structure, similarly to the γ-ray light curve. A follow-up paper elaborating on the optical light curves of J1048+7143 is in preparation.

Zoom In Zoom Out Reset image size

Due to gravitational radiation, the dominant spin S ₁ slowly approaches the total angular momentum J ; therefore, the jet meets the target at different angles as the binary merger is progressing. After several spin–orbit periods, the jet just does not meet with the target again, as the precession cone of S ₁ is narrowing. The ∼90° misalignment between the parsec-scale and kiloparsec-scale jets of J1048+7143 is indeed compatible with a secular change in the jet direction (Kun et al. 2022b). Qualitatively it explains why we do not see a flaring radio behavior before 2008—the jet is circling somewhere away from our line of sight. The peak flux of the γ-ray subflares is first increasing and then decreasing. It can also be explained in the framework of this scenario, with the combination of changing characteristic relativistic boosting of the jet and column density changes in the target.

Since in Phases I and III the emission of high-energy neutrinos is expected from p–p interactions, J1048+7143 might be an interesting source to search for high-energy neutrinos in archival neutrino data sets. We encourage further multiwavelength and multimessenger studies of this source.

E.K. thanks the Alexander von Humboldt Foundation for its Fellowship. J.T. and I.J. acknowledge support from the German Science Foundation DFG, via the Collaborative Research Center SFB1491: Cosmic Interacting Matters—from Source to Signal. Support from the Hungarian National Research, Development and Innovation Office (NKFIH) is acknowledged (grant No. OTKA K134213). This project has received funding from the HUN-REN Hungarian Research Network. L.C. was supported by the CAS "Light of West China" Program (No. 2021-XBQNXZ-005) and the Xinjiang Tianshan Talent Program. This paper makes use of publicly available Fermi-LAT data provided online by the Fermi Science Support Center (http://fermi.gsfc.nasa.gov/ssc/data/access/). On behalf of Project "fermi-agn," we are grateful for the usage of the HUN-REN Cloud that significantly helped us achieve the results published in this paper. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under the cooperative agreement by Associated Universities, Inc. The XAO-NSRT is operated by the Urumqi Nanshan Astronomy and Deep Space Exploration Observation and Research Station of Xinjiang (XJYWZ2303). We acknowledge the use of data from the Astrogeo Center database maintained by Leonid Petrov.