Observational Evidence of a Centi-parsec Supermassive Black Hole Binary Existing in the Nearby Galaxy M81

Studying a centi-parsec supermassive black hole binary (SMBHB) would allow us to explore a new parameter space in active galactic nuclei, and these objects are also potential sources of gravitational waves. We report evidence that an SMBHB with an orbital period of ∼30 yr may be resident in the nearby galactic nucleus M81. This orbital period and the known mass of M81 imply an orbital separation of ∼0.02 pc. The jet emanating from the primary black hole showed a short period of jet wobbling at ∼16.7 yr, superposing a long-term precession at a timescale of several hundred years. Periodic radio and X-ray outbursts were also found two times per orbital period, which could be explained by a double-peaked mass accretion rate variation per binary orbit. If confirmed, M81 would be one of the closest SMBHB candidates, providing a rare opportunity to study the final parsec problem.


INTRODUCTION
The nearby ∼3.96 Mpc (Bartel et al. 2007) spiral galaxy M81 (NGC 3031), harbors one of the nearest low luminosity active galactic nuclei (Nagar et al. 2002;Ho 2008).The almost face-on aspect with an inclination angle to the line of sight (LOS) 14 ± 2 degrees of the host galaxy allows a clear and unobscured view of the nucleus.With a black hole mass M ∼ 7 × 10 7 M ⊙ determined spectroscopically (Devereux et al. 2003) with the Hubble Space Telescope (HST), an angular resolution of 1 micro-arcsecond (µas) corresponds to a linear size of 6 gravitational radius (r g = GM c 2 ).As seen with very long baseline interferometry (VLBI), the AGN in M81 (i.e., M81*) was identified with a compact core and a one-sided jet to the northeast (Bietenholz et al. 1996(Bietenholz et al. , 2000;;Markoff et al. 2008).The core-jet structure was confirmed by locating the core at multiple frequencies (Bietenholz et al. 2004), and the frequency-size relation of the core was measured up to 88 GHz (Ros & Pérez-Torres 2012;Jiang et al. 2018).The first series of VLBI observations on M81* were carried out in 1980s (Bartel et al. 1982).Since the early 1990's, M81* has been frequently observed, mainly to monitor the evolution of the adjacent supernovae SN1993J (Marcaide et al. 1995(Marcaide et al. , 1997(Marcaide et al. , 2009;;Bartel et al. 2000Bartel et al. , 2002Bartel et al. , 2007;;Bietenholz et al. 2001Bietenholz et al. , 2003) ) and SN2008iz (Kimani et al. 2016), as well as to track the relative motions between M81 and M82.A long-term VLBI monitoring of M81* showed strong evidence of jet precession with a period of ∼7.27 years, it was also found a flare in the peak flux densities, which lasted from mid 1997 to mid 2001 (Martí-Vidal et al. 2011).But the long-duration flare did not repeat in a late observation campaign in 2005 (Markoff et al. 2008).The jet precession and its position angle (PA) drift were confirmed with new data in 2017-2019(von Fellenberg et al. 2023)), while it could not constrain the model well due to the sparse sampling in time domain.
In this paper, we collected all the archived VLA/VLBI and also our newly VLBI observational data from 1980 to 2022 at multifrequency bands, as well as the longterm X-ray monitoring data covering the same period.We found M81* underwent a periodic change of its jet PAs, with three main outbursts in both the radio and the X-ray.We propose a scenario involving a supermassive black hole binary in an orbit misaligned with the accretion disk that exists in M81*.The binary model fits the jet PAs well, and can explain the repeated outbursts.In Section 2 we summarize the observations and data reduction of M81*, we describe our results in Section 3 and the fitting by the proposed binary model in Section 4. Our discussions are presented in Section 5, followed conclusions in Section 6.  1.
The VLBI data reduction followed standard procedures to the amplitude and phase calibrations in AIPS and Difmap, described in references (Martí-Vidal et al. 2011;Jiang et al. 2018).Several iterations of phase selfcalibration were executed in Difmap to make the final image, prior to an amplitude self-calibration with a solution interval of 10 minutes.The VLBI flux density measurements were estimated from the source core.Since the VLBA is an homogeneous array, the uncertainty of the amplitude calibration is ∼5% (Kovalev et al. 2005).Hence, the VLBA data can be used as a template to constrain the initial source model in joint (i.e., global) observations.In the case of KaVA, we observed the amplitude calibrator 4C 39.25, whose flux density was ∼4.84 Jy at 22 GHz and ∼1.97 Jy at 43 GHz in our observations, almost at the same level as that in previous observations (Niinuma et al. 2014).While M81* showed a flux density of at least two times higher than its normal stage, the flux at 43 GHz was even higher than that at 22 GHz, indicating radio outbursts during these periods.
VLA Observations The archived observational data of the interferometer VLA at 5.0, 8.4 and 15 GHz from 1980 to 1991 were re-analyzed with AIPS and identified an outburst around the middle of 1985.The VLA data around 1999 and 2015 were added for increasing the cadences during outburst phases.The 2014-2015 data was utilizing the CASA software package (CASA Team et al. 2022), version 5.5.0.The origin of CASA lies in the AIPS++ project and it is better suited for processing the datasets of 2014-2015 observations.The visibilities were examined for radio frequency interference, which primarily affected the low frequency bands.Standard reduction techniques were used to transfer the flux density scale from the flux density calibrator to phase calibrator and M81*, and then bandpass, gain calibrations, and primary beam correction (if necessary) were done.A Gaussian centered at the phase center was then fit to the visibilities.A 10% error was added to the statistical fit error, to account for the uncertainty in the absolute flux density.The results are presented in Figure 1B and Table 1.

X-ray
We analyzed the soft X-ray emission from the nuclear source of M81* from 2005 to 2022, using the archival observations of the X-ray Telescope (XRT) on board the Neil Gehrels Swift Observatory.The cleaned events were reduced by the tool xrtpipeline.The source spectra were extracted from a circle region with a radius of 30 pixels, for photon counting (PC) mode, and a box region with a width of 60 pixels and a height of 20 pixels, for window timing (WT) mode.For the observations suffering pile-up effects (count rate >0.6 c/s for PC mode), we excluded the central pixels determined by using the updated XRT point spread function (PSF)1 .The background spectra were extracted from an annulus region with inner radius of 180 pixels and outer radius of 200 pixels, to avoid other point sources in the Swift/XRT filed of view.
The X-ray spectra of Swift/XRT were fitted with a model tbabs*(powerlaw) in the energy range of 2-10 keV, where the column density n H is fixed as the Galactic value 7.22×10 20 cm −2 towards to the direction 2 .The X-ray spectral fitting was performed with Bayesian Xray Analysis (BXA, Buchner et al. 2014) within PyXspec (Gordon & Arnaud 2021).Then, we used the convolution model cflux to calculate the flux in the energy range 2-10 keV.We found a negative correlation between the photon index and X-ray flux, which was fitted with a function as Γ = −0.55 × log(F 2−10keV ) − 4.25, implemented in the software package UltraNest (Buchner 2016(Buchner , 2019(Buchner , 2021)).
In order to investigate the long term variability of M81*, we also used the flux measurements of other X-ray instruments in the literature.Elvis & van Speybroeck (1982) reported the 2-10 keV X-ray flux of two Einstein observations during 1978-1979 (Figure 1).Petre et al. (1993) reported the X-ray flux in 2-10 keV of Einstein, EXOSAT , Ginga, ROSAT and BBXRT (Figure 1).La Parola et al. (2004) reported the X-ray flux in 0.5 -2.4 keV for more than twenty years, including the X-ray mission ROSAT,ASCA, Chandra and XMM-Newton.In order to convert the X-ray flux into 2-10 keV, we used the photon index 1.79 for ROSAT and ASCA, and the reported photon index of Chandra and XMM-Newton in La Parola et al. (2004).

Position angle
The compact core region of M81* at all frequencies is elongated in the direction of the jet extension and shows a flat or inverted spectrum (Bietenholz et al. 1996(Bietenholz et al. , 2000;;Martí-Vidal et al. 2011;Jiang et al. 2018).Hence, we fitted the VLBI calibrated visibility data with elliptical Gaussian models.The radio emission from the core could be well fitted by a single Gaussian (Bietenholz et al. 2000;Martí-Vidal et al. 2011;Jiang et al. 2018).The jet orientation in this model corresponds to the orientation of the Gaussian main axis.A similar approach has been applied by the team of Monitoring of Jets in Active Galactic Nuclei with VLBA Experiments and demonstrated to be a robust way to derive the core parameters (Lister & Homan 2005;Kovalev et al. 2005).Besides this model of an elliptical Gaussian (plus an extended circular Gaussian, to account for the more distant jet emission), we also fitted the data using two close-by circular Gaussians, which resulted in fits with a similar quality.For this alternative model, the PA of the jet can be estimated from the orientation between the two Gaussian components.The PAs derived from the elliptical-Gaussian fitting are consistent with those from the two circular-Gaussian components.We adopted the PAs estimated from the elliptical Gaussians from 1993 to 2005 (Bietenholz et al. 2000;Martí-Vidal et al. 2011), to cover another radio and X-ray outbursts.While the practical PA errors are three times of the norminal errors in their papers.In the simultaneous multi-frequency observations, the PAs at 5.0, 15, 22, and 43 GHz were different to those observed at 8.4 GHz by 1.9 ± 3.3, −1.1 ± 1.6, −1.3 ± 0.7, and −3.1 ± 1.4 degrees, respectively, implying a slight frequency-dependent jet orientation, likely related to opacity effects (Bietenholz et al. 1996;Martí-Vidal et al. 2011).These slight PA differences were considered and interpolated to the epochs without observation at 8.4 GHz.We assembled all the PAs of VLBI observations from 1981 to 2022 at different frequencies to that of 8.4 GHz, and show the results in Figure 1A.Six VLBI epochs at 8.4 GHZ are also shown in Fig. 1D, where the jet wobbling can be seen around a projected on-sky axis (shown in red).There is a hint of periodic behavior, which can be fitted using a sinusoidal model.

Light curves
The radio light curve is assembled to the flux densities at 8.4 GHz or interpolated at 8.4 GHz with a spectral index if the flux density measurements other than that at 8.4 GHz are available.The averaged spectral index α of M81 core is −0.11,F υ (υ) ∝ υ −α .It appears that there are repeated outbursts in the radio and Xray light curves (Figure 1).We search the periodicity in the time series of PA by using the Lomb-Scargle periodogram (LSP) calculated by astropy (Astropy Collaboration et al. 2013Collaboration et al. , 2018Collaboration et al. , 2022)).There is an obvious peak at 15.2 years (Figure 2C), which indicates that the PA shows a periodic variation with ∼ 15 years.We also detect a peak at ∼ 15.2 years in the LSP of X-ray light curve (Figure 2A), which is consistent with that of PA.Meanwhile, there are two comparable peaks at 27.6 and 12.5 years in the radio LSP (Figure 2B).The differences of periods between radio and X-ray light curves should be due to the sparse sampling in the 1980s.The measuring of the periodicity would be improved with simultaneously intensive VLBI and high-energy monitoring in the future.

SUPERMASSIVE BLACK HOLE BINARY MODEL
In a binary black hole system, where the disk plane is misaligned with respect to the orbital plane, the jet of the primary black hole is expected to undergo jet wobbling and nodding besides of the jet precession (see Fig. 3), due to the oscillating torque from an orbital companion (Katz et al. 1982;Bate et al. 2000).The wobble is along the precession direction and its rate ω b is approximately twice of the orbital rate (ω o ), while the amplitude in radians is roughly Ω p /2ω o , where Ω p is the precession rate Ω p ≈ 0.05ω o .The nodding period is the same but the amplitude is multiplied by a factor of tan θ p , where θ p is the inclined angle of the orbit plane to the disk plane.The secondary is orbiting around the primary black hole with an angular velocity of ωo and a separation a.The orbital plane is misaligned θp to the disk plane of primary black hole.ξ0 and ϕo are the projection of the orbital axis z orb in the sky and its inclined angle to the LOS, respectively.The disk-jet of primary black hole undergoes short term (2π/ω b ) nodding and wobbling (smaller purple ellipse), superimposing a long-term (2π/Ωp) precession (blue ellipse).

Observer
A series of rotations are used to transform the motion of the jet vector n j = (0, 0, 1) to the observer's frame.All the rotations are in right handed coordinate systems.The first step is to be parallel to the orbital axis.n orb in the orbital frame can be calculated using the following equations where t 0 is set to 1980 as the reference time here.Two more rotations are applied to transform to the observer's frame.
where ξ 0 and ϕ o is the projection of the orbital axis in the sky and its inclined angle to the LOS, respectively.

Bayesian analysis
The Bayesian analysis based on Markov Chain Monte Carlo (MCMC) is one of the most effective approaches to infer posterior distributions of model parameters, which describe the data while marginalizing over nuisance parameters and taking into account systematic uncertainties.It becomes a ubiquitous technique in astronomy (Sharma 2017).A Bayesian analysis using an MCMC algorithm (Foreman-Mackey et al. 2013) is performed with six parameters (Ω p , θ p , ω b , η p , ξ 0 , and ϕ o ) in the binary black hole model.The PA likelihood is assumed to be Gaussian, − ln L P A = i . In addition to the PA measurements, we consider joint constraints with the jet viewing angles θ v,2011 < 56 • in 2011 (King et al. 2016) and θ v,2018 < 15 • in 2018 (Wang et al. in prep).The corner plot of MCMC fitting results is presented in Fig. 4. All the results (given with 1σ uncertainties) are presented here.The jet wobble and precession model of a SMBHB fits the evolution of PA well with a reduced χ 2 equal to 1.21 (the red line in Figure 1A).The MCMC results tell the jet is wobbling every ∼ 16.7 years, and precessing with an additional long period ∼ 785 years.The binary model prefers an orbit highly inclined or misaligned to the inner accretion disk with θ p = −54 +10 −15 • .

Repeated outbursts
In the relativistic beaming theory, the observed and rest-frame flux densities are linked by F υ (υ) = δ 2+α F ′ υ ′ (υ), where the Doppler factor δ = [γ(1 − β cos θ v )] −1 and the Lorentz factor γ = (1 − β 2 ) −1/2 (Raiteri et al. 2017).We use a bulk velocity β = 0.71 in the jet core given the apparent speed of 0.51c (King et al. 2016) and the viewing angle 12.8 • from fitting results.The averaged spectral index and the Doppler factor incorporating with the fitted viewing angle θ v are applied for the radio light curve.The X-ray light curve (Fig 1C) is also compensated the beaming effect due to the changes of viewing angles.The detection of relatively slow apparent speeds of knots (∼ 0.5 c near the core in 2011 and ∼ 0.1 c down the jet in 2018) indicates the Doppler beaming effect in the jet should be small.The fitted descending viewing angles are not responsible for the repeated outbursts in the light curves.The LSP results of original light curves and those corrected the beaming effect are almost identical.We suppose a compact object circling the primary black hole in a Keplerian orbit.When the compact object goes through the accretion disk, the mass accretion rates are raised and therefore cause outburst periodically.A binary system can enhance the mass accretion rate twice per binary orbit (Hayasaki et al. 2013;Gold et al. 2014).

Implications of the Observational Results
Implied from the SMBHB model as shown in Fig. 3, the troughs of fitting viewing angles in each wobbling period occur when the secondary black hole goes through the accretion disk.The scenario is similar to the most promising SMBHB candidate OJ 287, although the origin of flares is proposed to be impact flashes (Valtonen et al. 2008) or precession-induced (Britzen et al. 2023).The occurence time of troughs in viewing angles are consistent with those of flux peaks in light curves.Interestingly, the intervals among three peaks in the light curves are similar but not identical.That means the two nodes in an orbit are not symmetric and the orbit period should not be exactly twice of the jet wobbling period.We adopt the time interval of 29.7 ± 0.4 yr between the peaks in 1985 and 2015 as the orbit period, which corresponds to an orbit separation a ∼ 0.02 pc (∼ 1 mas at the distance ∼ 4 Mpc) via the Kepler's third law.The post-Newtonian parameter (Will 2011), ε = GM c −2 r −1 is about 1.5×10 −4 .A value of 10 −8 < ε < 10 −3 confirms that it is a bounded wide binary system where dynamical friction is dominant at this stage.A mass ratio of the companion to the primary black hole, 0.07 < q < 0.10, could be constrained through its relations to the outer disk radius and the precession period, the equations ( 14) and ( 16) in Britzen et al. (2018), respectively.A small mass ratio indicates the companion is probably radio quiet due to no normal accretion disk at most of its orbital phase.The gravitational wave (GW) strain h, as expressed by the equation ( 6) in Yan et al. (2020), is about 3 × 10 −17 , which is more than one order lower below the detection limit of the strain-sensitive GW detector, international pulsar timing arrays (Antoniadis et al. 2022;Agazie et al. 2023;EPTA Collaboration, & InPTA Collaboration 2023;Reardon et al. 2023;Xu et al. 2023).The parameters of the SMBHB model are listed in Table 2.We estimate the separation when gravitational radiation drives the binary to merger, it is about 0.004 pc.The coalescence timescale in the GW-dominated regime is ∼3.9 million years (Loeb 2010).While the inspiral from 0.02 pc to 0.004 pc could be driven by the gasassisted orbital contraction, which depends on the accretion dynamics of circumbinary disk and the binary system itself.This is related to the interesting "finalparsec problem" (Begelman et al. 1980;Lai & Muñoz 2023) and should be investigated with further studies.
Regarding the jet launching of the primary black hole, two widely adopted scenarios have been proposed.The outer disk radius of inner accretion disk is ∼ 1000 r g , far less than the orbital separation ∼ 6000 r g , which suggests that the secondary black hole is not significantly affecting the Blandford-Payne (BP) jet launched by the primary disk (Blandford & Payne 1982).Blandford-Znajek (BZ) jet is coupled to the spin of black hole (Blandford & Znajek 1977), it is thus more independent to the companion black hole.Current data can not distinguish the scenario of BP or BZ.In the case of mass ratios ≪1, the accretion and variability are dominated by the accretion flow onto the primary (Gold et al. 2014).As proposed in Section 4, the jet precession of primary black hole is due to the misaligned disk to orbital plane (Katz et al. 1982;Bate et al. 2000;Lai & Muñoz 2023).We can not distinguish the scenario of disk-jet precession or precession of disk-jet and its central black hole as a whole either.

Other Possible Scenarios
A jet launched by a tilted accretion disk of spinning SMBH undergoes Lense-Thirring (LT) precession (Fragile et al. 2007;Liska et al. 2018).The sinusoidal changes of PA could be related to a precessing jet.The longterm precession modulation (Ω p ) could be contributed by an outer elliptical-disk (Bower et al. 1996;Storchi-Bergmann et al. 2003).Including Ω p , there are six parameters in the LT precession model (ω LT , θ LT , Ω p , η 0 , t 0 , and Φ o ).We also perform a Bayesian analysis using an MCMC algorithm.We consider joint constraints of the PA measurements together with the jet viewing angles θ v,2011 < 56 • in 2011 and θ v,2018 < 15 • in 2018, as well as a further constraint is given by the inclination angle 14 ± 2 • of the LOS to the normal direction of the gaseous disk (Devereux et al. 2003), Φ o = 194±2 • as the LOS is in the opposite direction to the observer.The fitting results indicate a 17.5 ± 1.5 year periodic LT precession (2π/ω LT ), superimposed to a long-term rotation period (2π/Ω p ) of 628.3 ± 63.0 years.The tilted angle of the inner thick disk is θ LT = 1.7 ± 0.6 degrees.
A small tilt angle implies the disk precession by LT should have entered a stable stage.According to the equation ( 43) in Fragile et al. (2007), the short period of LT precession indicates a truncated inner disk at 160 − 170r g , incorporating a black hole mass of 7×10 7 M ⊙ and a spinning parameter a • = 0.9.However, the radius of truncated inner disk by LT is larger than a typical value of several tens r g near the black hole horizon.And forming an outer elliptical-disk implies a possible companion.Additionally, the repeated outbursts could not be explained by the Doppler boosting of a reduced viewing angle in the case of LT precession.Using the fitted viewing angles, the peaks in the light curves (Fig. 1B) should be produced with a moderate Lorentz factor Γ ∼ 12.9, a relativistic velocity β = 0.997 of the emitting source in units of the speed of light.This is inconsistent with the slow motions of knots detected in the jet of M81*.
Lense-Thirring precession of the SMBH itself has a normal period of millions of years.Otherwise, the spinning parameter of M81* black hole would have to be very small (∼ 10 −12 for a timescale of the order of a decade), according to equation (10) in Britzen et al. (2018) with a viscosity parameter α vis = 0.1.The scenario of the SMBH precession is thus rather implausible.

CONCLUSIONS
The jet precession in M81* is investigated with fourdecade data thoroughly both in radio and X-ray observations.Owning to the proximity of M81*, the high resolution VLBI approach the inner core region adjacent to the black hole, and reveals that the M81* jet is undergoing a short period of jet wobbling ∼ 16.7 years, superposing a long-term precession at a time-scale of several hundred years.Combining the PA evolution and the light curves, this can be best explained by oscillations of a secondary black hole in a binary system, rather than black-hole precession or tilted disk-jet precession due to LT effect.The binary system has an orbital separation of ∼ 0.02 pc and an orbit period of ∼ 30 yr.This makes M81* as one of the closest SMBHB candidate and a potential target to resolve the "final parsec problem".Nevertheless, we notice that similar claims about existing binary black holes have been made before but turn out to be challenging, as evidenced by the detection of the offset of VLBI component not associated with the VLBI core (Roland et al. 2013) or possibly dual inverted-spectrum radio cores (Kharb et al. 2017), as well as the extensive literature on OJ287 and its multiple unsuccessful predictions.Some recent evidences even contradict the existence of supermassive binary black holes, e.g. the SMBHB candidate J1502SE/SW are double hotspots (Wrobel et al. 2014), PSO J334 is most likely a jetted active galactic nucleus with a single SMBH (Benke et al. 2023), against the binary scenario evidenced by periodic variations in its optical light curve.Future intensive contemporary VLBI and high-energy monitoring campaigns are necessary to further verify the SMBHB scenario in M81*.
The authors thank the anonymous referee for very constructive suggestions in improving the manuscript.WJ thanks Fangzheng Shi, Zhiyuan Li for useful discussions.This work was supported in part by the National Natural Science Foundation of China (Grant No. 12173074, 11803071).IMV thanks the Generalitat Valenciana for funding in the frame of the CIDEGENT/2018/021 and ASFAE/2022/018 Projects, the MICINN Research Project PID2019-108995GB-C22, and the Astrophysics and High Energy Physics programme by MCIN, with funding from European Union NextGenerationEU (PRTR-C17I1).ZY is supported by the Natural Science Foundation of China (grants U1838203,U1938114), the Youth Innovation Promotion Association of CAS (id 2020265) and funds for key programs of Shanghai astronomical observatory.YL is supported in part by the Natural Science Foundation of Shanghai (grant NO. 23ZR1473700).We are grateful to all staff members at Effelsberg-100m, EVN, GBT, KaVA, Sheshan & Tianma-65m, VLA, and VLBA, who helped to operate the array and to correlate the data.VLBA, VLA, GBT, and the data archive system are operated by the National Radio Astronomy Observatory, which is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.The European VLBI Network is a joint facility of independent European, African, Asian, and North American radio astronomy institutes.The KVN is a facility operated by the Korea Astronomy and Space Science Institute.VERA is a facility operated by Na-tional Astronomical Observatory of Japan in collaboration with Japanese universities.

Figure 1 .
Figure 1.Radio and X-ray variations of M81* from 1980 to 2022.a) Jet position angles observed at multifrequencies and assembled at 8.4 GHz.The solid curves are the Bayesian fitting to the SMBHB model and results of evolution of PA (red) and viewing angle (blue).b) Radio light curve formed assembled at 8.4 GHz (grey: original; blue: Doppler beaming corrected).c) X-ray (2-10 keV) light curve (grey: original; dark: Doppler beaming corrected).d) 6 epoch VLBI images at 8.4 GHz from 2005 to 2022, all the uniform-weighted interferometry beams shown in bottom left corner are restoredto a normal beam size (0.5 mas×0.5 mas), which is the intercontinental interferometry including VLBA and EF or SH at 8.4 GHz.The contour levels are 0.3×(-1, 5, 10, 20, 40,  80,160, 240)  mJy beam −1 .The lowest level is about three times root-mean-square of noise in the image.The red line indicates the projection on the sky of the rotation axis of M81 galaxy at 62 •(Bartel et al. 1982); the dashed green lines indicate the position angles of jet in the center region from the results of model-fitting.

Figure 2 .
Figure2.Lomb-Scargle periodograms (LSP) of time series of X-ray (upper panel), radio fluxes (middle panel) and PAs (lower panel) for M81*.A peak at ∼15.2 years can be detected in the LSP of PAs and the X-ray light curve, respectively.There are two comparable peaks at 27.6 and 12.5 years in the radio LSP.All these time series measurements suggest a periodic mechanism in the central region.

Figure 4 .
Figure 4. Corner plot of the MCMC fitting to a binary black-hole model in M81*.The six parameters in the binary black hole model Ωp, θp, ω b , ηp, ξ0, and ϕo is the precession rate, the inclined angle of the orbit plane to the disk plane, the wobbling and nodding rate, the initial projection angle of jet on the orbital plane at t0, the projection of the orbital axis in the sky and its inclined angle to the LOS, respectively.
Radio VLBI Observations M81* has a flat or inverted spectrum at a flux density level of ∼ 100 mJy.A high declination of ∼ 70 • makes it circumpolar for most VLBI observatories in the Northern hemisphere, hence visible in common with all telescopes.We observed M81* with the Korean VLBI Network and VLBI Exploration of Radio Astrometry array (KaVA) in 2015 at 22 and 43 GHz, the Very Long Baseline Array (VLBA) in 2016, 2018, 2019, 2021 and 2022 (the latest epoch had VLBA plus Tianma-65m and Sheshan-25m telescope in Shanghai) at multiple frequencies covering 5.0, 8.4, 22 and 43 GHz.We also analyzed the archived VLBI observational data since 2005, which mainly included the VLBA and Effelsberg-100 m (EF) telescope, some epochs had stations of European VLBI network, Greenbank (GB) Telescope and Very Large Array (VLA) joining.The VLBA, with additional inter-continental baselines, resolved M81* up to the ground-based highest resolution at each frequency.The usual 8 or 12 hour tracks for each epoch made the spatial-frequency sampling, known as UV coverage, rather complete for the image reconstruction.The published VLBI results of M81* in 1981 and from 1993 to 2005, with similar array configurations, were assembled for a complete understanding the periodic phenomenon.A summary of observations since 2005 is presented in Table

Table 2 .
Parameters of SMBHB system in M81*