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

Analyzing single-dish and very long baseline interferometry radio, as well as Fermi Large Area Telescope γ-ray observations, we explained the three major flares in the γ-ray light curve of FSRQ J1048+7143 with the spin–orbit precession of the dominant mass black hole in a supermassive black hole binary system. Here, we report on the detection of a fourth γ-ray flare from J1048+7143, appearing in the time interval that was predicted in our previous work. Including this new flare, we constrained the mass ratio into a narrow range of 0.062 < q < 0.088, and consequently we were able to further constrain the parameters of the hypothetical supermassive binary black hole at the heart of J1048+7143. We predict the occurrence of the fifth major γ-ray flare that would appear only if the jet will still lay close to our line of sight. The fourth major γ-ray flare also shows the two-subflare structure, further strengthening our scenario in which the occurrence of the subflares is the signature of the precession of a spine–sheath jet structure that quasiperiodically interacts with a proton target, e.g., clouds in the broad-line region.


INTRODUCTION
In the recent years, multimessenger astronomy became one of the most rapidly evolving fields in astronomy.As being the synergy of coordinated observations Corresponding author: Emma Kun emma@tp4.ruhr-uni-bochum.detargeting 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 highenergy 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,b).Patterns of high-energy particles 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 yet (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 (GZK) horizon (e.g.Stanev et al. 2000).High-energy γ-photons lose also their energy on their way to Earth in interactions with the photons of the Cosmic Microwave Background (CMB).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 gravitational wave 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 a SMBBH (Gergely & Biermann 2009) in the heart of the galaxy.This active galactic nucleus (AGN) indeed belongs to a group of blazars which 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 multi-epoch 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 a 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, therefore its uncertainty (90% containment) is a relatively large ∼ 3.58 • (the neutrino arrived with the energy of ∼ 170 TeV and the event has a signalness of 42%, Blaufuss et al. 2022).In this paper, we extend the γray light curve of J1048+743 with 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.

THE UPDATED γ-RAY LIGHT CURVE OF J1048+7143
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 right ascension RA J2000 = 162.• 1067 and declination DEC J2000 = 71.• 7297.The region of interest (ROI) 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)  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 5% threshold of Table 1.Exponential fitting of the γ-ray light curve of 4FGL J1048.4+7143.Column (1) sets the fitted quantity (a measures the height of the respective flare, b the time location of its peak, and c its slope).The numbered columns (2)-( 7) contain the values of the parameters of the exponential functions fitted to the respective flares, where Fi,j means the j-th subflare (j=1 or j=2) of the i-th main flare (from 1 to 4).The fitted baseline is a0 = (0.12 ± 0.02) × 10 −7 ph cm −2 s −1 .3.20 ± 0.12 3.07 ± 0.08 3.06 ± 0.17 relative flux error.In the resulting light curve, shown in Fig. 1, we see a new major flare starting in the beginning of 2021 (after MJD 59200), in addition to the previously reported three flares (Kun et al. 2022b).This 4th major flare also shows a double sub-flare 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

Parameter
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 0 ).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 Fig. 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 (MCMC) 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).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).

JET KINEMATICS AT 8.6 GHZ, UTILIZING ARCHIVAL VLBI DATA
To derive the structural and kinematic properties of the mas-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 database1 .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-tomajor 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): where S ν is the flux density (in Jy), d 1 and d 2 are the major and minor axes of the core (in mas), 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 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 and then by fitting accelerating proper motions as where µ is the linear proper motion measured in mas yr −1 , t ej the ejection time, μ 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 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 Fig. 2. The resulted best-fit proper motions, as well as the derived apparent speeds are shown in Table 3.As it can be seen from the reduced χ 2 values in Table 3 and from Fig. 2, the accelerated proper motion fits give better results.
The apparent speed of the jet components is related to the bulk jet speed β and their inclination angle ι as (e.g.Urry & Padovani 1995) where the jet speed is related to the Lorentz factor Γ as The maximal β app occurs for the critical inclination cos ι c = β, sin ι c = Γ −1 , having the value β app,max = (Γ 2 − 1) 1/2 .In practice, it is unlikely that a jet component will exactly move in this direction, therefore the maximal speed β obs app,max ≤ β app,max seen in the VLBI jet gives a lower limit on the Lorentz factor of the jet Taking the accelerated proper motion of C1, µ max = (0.113 ± 0.002) mas yr −1 , that is the fastest jet component seen in J1048+7143 before the flaring behavior starts, we get the maximum observed apparent speed as β obs app,max = 6.56 ± 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 the light).Then the related critical inclination angle is ι c = arcsin(Γ −1 low ) ≈ 10.9 • and the half-opening angle of the jet is estimated as ζ ≈ 1/Γ ≈ 8.3 • (Boettcher et al. 2012).

SPIN-ORBIT PRECESSION AND PREDICTION OF THE NEXT γ-RAY FLARE
In Kun et al. (2022b) we have shown that the γ-ray light curve of J1048+7143 (E > 168 MeV) and its singledish radio flux density curve obtained at 4.8 GHz (see bottom panel of Fig. 1) are consistent with the jet precession driven by the spin-orbit precession, occurring in a 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 2 /m 1 (m = m 1 + m 2 , m 1 > m 2 ).We note, the ∼ 0.1 Jy difference in flux density seen between the single dish radio flux density curves in Fig. 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 a 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 ε = 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 1 + S 2 + L N are conserved during the motion.If the S 1 and S 2 the spins of the two SMBHs (m 1 > m 2 ) are not parallel with the Newtonian orbital momentum L N , the spins begin to precess (Barker & O'Connell 1975, 1979): 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 till the final coalescence of the compact binary as:   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: 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 4th major γ-ray flare, the mass ratio of the SMBBH is now constrained as 0.062 < q < 0.088 (see Fig. 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 gravitational waves leads to shrinking orbits, the Newtonian orbital angular momentum L N and the direction of the dominant spin S 1 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 1 , 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 1 , we do not expect a flaring behavior anymore due to the weak Doppler boosting.Given the jet still boosted in the close future, a 5th γray flare is expected to arrive between MJD 60379 (2024 March 10) and MJD 61350 (2026 Nov 6).

DISCUSSION AND CONCLUSIONS
Periodicites 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. 2020Kun et al. , 2023;;Britzen 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 depend on the complexity of the observed behavior: while models relying on e.g.single black holes or the precession of 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, quasi-periodic oscillations (QPOs) 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 Fig. 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 3 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 Fig. 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 Fig. 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 Fig. 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 (Fig. 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 (AAVSO), also show 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.
Due to gravitational radiation, the dominant spin S 1 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 1 is narrowing.The ∼ 90 • misalignment between the pc-scale and kpc-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 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 highenergy neutrinos is expected from p − p interactions, J1048+7143 might be an interesting source to search for high-energy neutrinos in archival neutrino datasets.We encourage further multiwavelength and multimessenger studies of this source.
by adding Fermi-LAT data in the time range from 2022 Mar 14 to 2023 Jun 4 (MJD 59635-60099) in the energy range 100 MeV − 800 GeV.Now our analysis covers almost 15 years between 2008 Aug 4 and to 2023 Jun 4.

Figure 1 .
Figure 1.Gamma-ray and radio flares of J1048+7143.Top: The new exponential fit of the γ-ray light curve of J1048+7143 with 5% adaptive binning (best-fit parameters are shown in Table 1).The sum of the contributions from individual sub-flares is displayed as the blue continuous line.Application of the centroid method is indicated via the red crosses.Their weighted average determine the centers of main flares as green dashed line.Intersections of these lines with the blue lines indicate the duration of the flares, illustrated as gray areas.The red continuous line shows the latest data point available in the Fermi -LAT Light Curve Repository (2023 September 22 or MJD 60209).Bottom: The extended radio flux density curve of J1048+7143, obtained with the RATAN600 and Nanshan (Urumqi) radio telescopes at 4.8 GHz (see details in Kun et al. 2022b).

Figure 2 .
Figure 2. Linear (left) and accelerated (right) proper motion fits for jet components C1 and C2 of J1048+7143 at 8.6 GHz.Only observations before the flaring behavior of J1048+7143 started are used to make these plots.

Figure 3 .
Figure 3. Prediction plot for γ-ray flares from J1048+7143.Gray areas show the duration's of the four flares determined with the centroid method in Table 2.The blue hashed bands represent an overlap of the predictions using the period Pn→n+1, where n = 1 or 2. Invalid mass ratios are displayed via the red crossed area.This plot is based on the duration of the second flare and the jet half-opening angle ζ = 8.3 • .The latter was derived from the pc-scale jet kinematics of J1048+7143 seen at the f = 8.6 GHz observing frequency with VLBA.

Figure 4 .
Figure4.Geometry of the jet-precession-induced light curve variation.The jet, driven by the spin-orbit precession in a SMBBH, goes through a target that supplies target protons for the hadronic processes between them and highenergy cosmic rays accelerated in spine-sheath structured jet.
E.K. thanks the Alexander von Humboldt Foundation for its Fellowship.J.T. and I.J. acknowledges 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 number 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 https://fermi.gsfc.nasa.gov/ssc/data/access/Fermi Science Support Center.On behalf of Project 'fermi-agn', we thank for the usage of the HUN-REN Cloud that significantly helped us achieving 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.

Table 2 .
Flare characteristics using the centroid method.