Radio jet precession in M81*

We report four novel position angle measurements of the core region of M81* at 5GHz and 8GHz, which conﬁrm the presence of sinusoidal jet precession of the M81 jet region as suggested by Martí-Vidal et al. (2011). The model makes three testable predictions on the evolution of the jet precession, which we test in our data with observations in 2017, 2018


Introduction
The galaxy M 81 appears as a bright radio source and is located at a distance of 3.36 ± 0.34 Mpc (Freedman et al. 1994).The black hole in the center of M 81 belongs to the class of low-luminosity active galactic nuclei (LLAGN) and exhibits relatively weak radio emission (F ν=4.8 GHz ∼ 150 mJy; e.g., Brunthaler et al. 2006).It is the closest LLAGN to Earth.Due to its high apparent luminosity, it serves as an optimal test bed for characterizing this class of accreting black hole (Markoff et al. 2008).As such, it may be the best candidate for bridging the accretion processes of LLAGN to that of Sgr A* (GRAVITY Collaboration 2020).M 81 shows slow changes in radio flux density on yearly timescales (e.g., Ho et al. 1999).Further, it shows fast intra-day flare-like variability at millimeter wavelengths (Sakamoto et al. 2001).This radio and millimeter behavior is consistent with the van der Laan expanding-blob scenario (van der Laan 1966).Here, the millimeter variability is created by blobs that expand as they move along the jet, where they become observable in the radio (Ho et al. 1999;Sakamoto et al. 2001).In the X-ray, M 81 is detected with a luminosity of ∼10 40 erg s −1 and flux changes on the order of a few tens of percent (Ishisaki et al. 1996).The overall spectral energy distribution was modeled with a jet-dominated ADAF model by Markoff et al. (2008) using a set of simultaneous multiwavelength observations of M 81 and shows remarkable similarity to that of Sgr A* (e.g., von Fellenberg et al. 2018).
M 81*, the core region of M 81, is a regular target for global very-long-baseline interferometry (VLBI) observations, in particular since the explosion of the radio-luminous supernova SN 1993J (Ripero et al. 1993;Weiler et al. 2007).This supernova was located at a close angular separation from M 81* and was thus frequently used as a calibration source.In VLBI observations, M 81* is typically marginally resolved and shows a jet in the northeastern direction.Bietenholz et al. (2000) and Bietenholz et al. (2004) determined that M 81 is at rest with respect to SN1993J.Further, they found a frequencydependent shift in the peak brightness of M 81*.This is consistent with the known core-shift effect of many active galactic nucleus jets (Marcaide & Shapiro 1984;Lobanov et al. 1998;Kovalev et al. 2008), which is believed to be caused by the synchrotron self-absorption of photons in the jet plasma.In this picture, the apparent shift in the luminous component results from the increasing opacity (and thus increasing luminosity) as a function of wavelength (Blandford & Königl 1979;Konigl 1981;Falcke & Biermann 1995;Marscher & Travis 1996;Davelaar et al. 2018).M 81* shows a decreasing core size (Θ) as a function of frequency, with an almost linear relationship: Θ ∝ λ ∼0.9 (Bartel et al. 1982;Kellermann et al. 1976;Bietenholz et al. 2000;Markoff et al. 2008).This relationship holds down to millimeter frequencies, with a confirmation at 43 GHz (Ros & Pérez-Torres 2012) and one at 87 GHz (Jiang et al. 2018); the latter authors report a relation of Θ ∝ λ 0.89±0.03 .Martí-Vidal et al. (2011) confirmed the basic findings reported in earlier works.They further showed that the M 81* intensity peak is shifted as a function of frequency along the direction of the jet using VLBI measurements from 1993 to 2005.Its size increases with decreasing frequency.The authors determined the location of the jet base to within 20, µas of the black hole, and they constrained the black hole mass to ∼2 × 10 7 M based on the strongly magnetized accretion flow scenario (Kardashev 1995).Alberdi et al. (2013) extended the temporal baseline of observations to the year 2012 and confirmed the basic findings reported in Martí-Vidal et al. (2011).By fitting elliptical Gaussian models to the central source and, if detected, the jet component, they derived a sinusoidal modulation of the jet position angle.They found a precession period of ∼7 yr on top of a linear increase in the position angle of ∼0.5 • /yr.They found this precession to be present both in their 5 GHz observations and in their 1.7 GHz observations.However, they found a lag of 1.9 ± 0.4 years for the precession between the frequency bands, which they interpreted as a core-shift effect.In this picture, the jet shows a differential corkscrew-like precession, and different frequencies probe different regions along the jet.Lastly, they connected the observed jet precession with a four-year flare and argued that the increase in flux is caused by Doppler boosting of the jet along the line of sight.
Their model therefore provides three testable hypotheses: (1) a prediction of the position angle of the core component as a function of time; (2) a prediction of this modulation at different frequencies; and (3) a prediction of the expected flux level at a given time point.
In this Letter we investigate whether the three predictions hold against new VLBI observations of M 81* obtained in the years 2017, 2018, and 2019.

Observations and data reduction
In total, we analyzed four sets of observations of M 81*.The first set, obtained by the European VLBI Network (EVN) at 5 GHz in June 2017, was a dedicated observation to measure the position angle of M 81* (Prog.ID ED042, PI Jordy Davelaar).The second set consists of an observation at 8 GHz obtained in June 2018 at the Very Long Baseline Array (VLBA; Prog.ID BJ090, PI Wu Jiang), with the intent to obtain phase-referenced observations of the core-shift effect in M 81*.Lastly, we used 5 GHz and 8 GHz VLBA observations obtained in 2019 for the purpose of studying the jet components in M 81* (Prog.ID BJ099, PI Wu Jiang).All data were reduced using the rPicard1 VLBI pipeline version v7.1.5(Janssen et al. 2019), which makes use of CASA v6.5 (THE CASA TEAM 2022) and the latest VLBI features (van Bemmel et al. 2022)2 .To derive the position angle, we fit all data with a single Gaussian model using Difmap (Shepherd et al. 1997).In all cases, the source is marginally resolved; however, we opted for a simple description by a singular Gaussian component to derive a robust measurement of the position angle.Table 1 reports the dates, frequency bands, derived values, and the respective imaging χ 2 reduced values.The corresponding maps and models are shown in Appendix D.
In order to obtain an as complete picture of M 81* as possible, we searched the VLBA and EVN archives for available observations since 2012.Several observations at higher frequencies exist, which we do not study in this Letter (e.g., BJ086  and BB303).Three more observations in the L, S , or X band exist: RP023A (EVN), RP023B (EVN), and BD185 (VLBA).They, unfortunately, do not allow a determination of the position angle as the data quality is not sufficient.For RP023A, no long baseline stations participated in the observations.For RP023B, several telescopes suffered from sensitivity losses at the start of the scan, possibly due to being late on source.When the bad measurements were flagged, the scan durations were no longer long enough to obtain good fringe solutions.For BD185, we find large instrumental delay corrections of ∼100 ns from the calibrator sources and only a few robust fringe detections on M 81 over the full Nyquist search window.These issues possibly originate from an error in the clock search at the correlator.
We further validated that we can reproduce the values published in Alberdi et al. (2013) with one example (BB293 B), which gives consistent results.

EVN observation
We analyzed a set of 5 GHz VLBI observations carried out by the EVN observatory, which were obtained on 21 June 2017 and which lasted for about 12 h.The IR, YS, TR, SH, NT, MC, WB, JB, and EF3 stations participated in the observations, and apart from a full loss of the right hand polarization at the EF station, no major technical difficulties occurred.The data were calibrated using observations of J0958+6533.

VLBA BJ090 observation
We analyzed the 8 GHz subset of the observations carried out in the VLBA BJ090 observation campaign from 10 June 2018.The data included the BR, FD, HN, KP, LA, MK, NL, OV, PT, and SC stations and used OJ287, J0954+658, and J1331+305 as calibration sources.No major technical difficulties occurred during the observations, with overall good data quality.The data were reduced in the same fashion as the EVN observations.In addition to the 8 GHz observations (X-band), the BJ090 observations included K-, Q-, and W-band observations, which we did not include in the analysis.

VLBA BJ099 observation
We analyzed the 5 GHz and 8 GHz subsets of the observations obtained as part of the VLBA BJ099 observation campaign, which was carried out on 4 November 2019.The data included the BR, FD, HN, KP, LA, MK, NL, OV, PT, and SC stations and used OJ287 and J0954+658 as calibration sources.While the 5 GHz observations showed overall agreeable data quality, the 8 GHz data set suffered from poor observations by the PT and SC stations.In order to derive a position angle measurement in  the 8 GHz band, we excluded the baselines to the PT and SC stations, the latter of which contributes the longest baselines.For both observations, the fit is relatively poor, with reduced χ 2 values of between four and five.Apart from the 5 GHz and 8 GHz observations (C-band and X-band), the BJ099 observations included K-, Q-, and W-band observations, which we did not include in the analysis.

Results
In the following sections we test the three hypotheses presented in Martí-Vidal et al. (2011) and Alberdi et al. (2013).

Sinusoidal jet procession
We adapted the model proposed in Martí-Vidal et al. (2011): where θ(t) is the position angle as a function of time, θ 0 is the initial position angle value at time 1996 − t 0 , T is the oscillation period in years, and β is the linear drift slope in degrees per year.
We extracted the measurements presented by Martí-Vidal et al. (2011) and fit the data up to the year 2012.The second column of Unfortunately, it is not straightforward to determine the uncertainty of the position angle accuracy.We estimated the uncertainty of each data point by fitting the evolution of the data with a nonparametric Gaussian process4 .The uncertainty of the Gaussian process model was then used as a proxy of the uncertainty of a datum at a given time.This allowed us to do a model selection test against a simpler baseline model: a simple linear drift (Eq.(A.1), Fig. A.2). We report the best-fit values of this model in the fourth column of Table 2.The reduced χ 2 of the linear drift model is χ 2 linear ≈ 38.5, whereas the sinusoidal model has χ 2 sinusoidal ≈ 4.5.The sinusoidal model is favored over a simple linear drift.

Frequency-dependent jet modulation
Martí-Vidal et al. ( 2011) reported a sinusoidal modulation of the M 81* core position angle at both 5 GHz and 1.7 GHz and found a (1.9 ± 0.4) yr time lag between the frequency bands.This was interpreted as differential corkscrew-like precession, where the core-shift effect (Konigl 1981;Marscher & Travis 1996) leads to longer wavelengths probing regions at a larger separation from the black hole.The model therefore makes a prediction for our higher frequency observation at 8 GHz.Typically, a linear relation (r c ∝ ν −1 ) for the core shift is found (e.g., Sokolovsky et al. 2011), and we thus expect a similar lag of 1.9 yr in the 8 GHz observations, but in the opposite direction.Figure 2 shows that while the observations are consistent with the period and amplitude of the oscillation, the data favor a larger time lag of ∆P 8 GHz ∼ 3.5 yr.This result is driven by the 2019 observation, which suffered from poor data quality.If one disregards the 2019 data, the observations are consistent with a shift of 1.9 years (see Fig. B.1).This illustrates that the available data set is insufficient to test this hypothesis, and future observations are required for a more definitive statement.2011) and Alberdi et al. (2013), and the orange pentagon denotes our newest measurement.The gray line shows a best-fit Gaussian model, which approximates the observed flare.The vertical line indicates the peak of the Gaussian flare model.The Gaussian flare shape is shifted by the precession period (6.9 yr).

Martí
the direction of the observer (see also Ros & Pérez-Torres 2012).The top panel of Fig. 3 shows the best-fit model and the position angle measurements.We include the two newest 5 GHz flux density measurements in the bottom panel.We fit a simple Gaussian + flux offset model to the flare data.The best-fit model is plotted in gray in Fig. 3.If the flare is caused by the precession of the jet, one expects the reoccurrence of a similar flare shape shifted by the period.We illustrate this by plotting shifted realizations of the best-fit Gaussian flare model, shifted by multiples of the precession period.It is clear that the low flux measurements in the years after 2002 are inconsistent with such a scenario.However, given the sparse sampling of the light curve, it is possible that the observations miss flaring activity if the correlation to the precession period is not very tight.Britzen et al. (2018) analyzed the precession and nutation of the resolved jet components of OJ287.In this section we follow their  2000), has seven parameters, which are given in Table 3.While some of the older observations resolved the core structure, our prime observable is the position angle of the core on sky, η(t).Further, we do not detect a strong deviation from a linear trend of the core-position angle.Thus, we cannot constrain the periodicity of the jet precession, P p , directly; we can only constrain its value from the observed nutation and the linear trend.We fit the temporal evolution of the position angle using Eq. ( 8) of Britzen et al. (2018) and refer the reader to Appendix C for the mathematical details.We determined the posterior parameters of the model using dynesty (Skilling 2006;Feroz et al. 2009;Skilling et al. 2004;Speagle 2020).Figure 4 shows the posterior distributions of the precession period, P p , the nutation period, P n , and the precession cone half-angle, Ω p .As argued before, the nutation is driving the sinusoidal modulation of the position angle with a period P n ∼ 7 years, the linear trend is caused by a large-scale precession of the jet with a largely unconstrained period (i.e., more than 200 years and less than 1800 years), and the precession cone angle is constrained to be positive.

Astrophysical origin of precession
Jet precession is observed for many jetted active galactic nuclei, and its origin is typically attributed to being of a purely stochastic nature (i.e., without true periodicity), induced by disk instabilities in tilted accretion disks, or induced gravitationally by a pacemaker companion.Following the arguments presented by Vaughan et al. (2016), we cannot rule out a stochastic nature for the observed modulation in the M 81* core region position angle.Nevertheless, the apparently sinusoidal modulation of the precession angle proposed by Martí-Vidal et al. (2011) has so far withstood both the extrapolation to the future (albeit with a slightly increased linear slope) and the extrapolation to different frequencies (precession in the X band leading the C band).The latter is, however, not proof of periodicity, as one would also expect a frequency-dependent lag of the position angle for a purely stochastic jet modulation.Assuming that the observed modulation is truly periodic, we discuss two possible scenarios for its astrophysical origin.For an accretion disk that is not aligned with the black hole's spin axis, Lense-Thirring precession (Thirring 1918)  a nodal precession of test particle orbits.The strength of this effect is frequency dependent and thus results in a precessing warp of the accretion disk.Fragile et al. (2007) demonstrated that this can lead to a constant period precession of the disk around the black hole.Such a disk precession is thought to be the for the Type-C quasi-periodic oscillations observed in X-ray binaries (e.g., Stella & Vietri 1998) and can be used to estimate the black hole mass and spin using certain assumptions (e.g., Ingram & Motta 2014).Building on the simulations by Fragile et al. (2007), Liska et al. (2018) demonstrated that titled accretion flows are able to launch jets and that the jet precession is aligned with the precession of the disk.In this set of simulations, the amplitude of the jet precession depends on the separation from the black hole, dropping from 100 • at a separation of R ∼ 10R g to ∼50 • at 100R g .Further, Fragile et al. (2007) found an empirical relation between the precession period of the disk, T p ∼ 0.3(m/M ) s, which corresponds to roughly 0.2 years assuming a black hole mass of 2 × 10 7 M .Both the amplitude of the oscillation and the precession period are in slight tension with the observed values (2Θ p ∼ 7 • , P p ∼ 7 yr).Assuming a precession-nutation model, the amplitude of precession is more similar to the expected value (∼80 • ); however, the much longer period in this case is even harder to contextualize.In this simple argumentation, however, we have not accounted for projection effects, none of the simulations have been tailored to M 81* (i.e., the unknown spin and the accretion-disk tilt), and we have ignored the fact that the disk precession timescales are dependent on the initial disk mass.We therefore suggest that the Lense-Thirring precession scenario may well be applicable in the case of M 81*.
Finally, a gravitational pacemaker may explain the observed precession of the M 81* position angle.Such a scenario has been found in the binary system SS433 (e.g., Stephenson & Sanduleak 1977;Clark & Murdin 1978), where two large-scale precessing jets create a corkscrew-like structure.For this system, the precession is thought to originate from the gravitational torque of the donor star on the accretion disk of a compact object, which is either a black hole or a neutron star (slaved disk model; Roberts 1974;Waisberg et al. 2019).A similar scenario has been proposed for OJ 287 and 3C 345, where an orbital modulation is present in both the light curve and the jet components (Lobanov & Roland 2005;Britzen et al. 2018).Further, Caproni et al. (2013) found a 12.1 yr precession period in BL Lacertae.Britzen et al. (2018) found a precession period of ∼23 yr on top of a much faster ∼1 yr nutation period and derived a binary separation of 0.001 pc < d < 0.1 pc.For 3C 345, Lobanov & Roland (2005) found a precession period of ∼10 yr for the position angle, which exhibited a linear trend; this is very similar to the behavior of M 81* (a short, ∼7 yr, nutation period and a linear trend, which we interpret as a long-period precession of ∼800 yr).
If we interpret the inferred seven-year period as the orbiting period of the companion object, we can derive the binary semimajor axis as well as the gravitational binary merger time.An orbiting secondary black hole induces gravitational torques on the accretion disk, which results in the precessing motion in the opposite direction of the disk rotation.The precessing motion is accompanied by the short-term nutation motion caused by the torque of a similar magnitude.However, the amplitude is smaller than for the precession by the ratio of the precession and the orbital frequencies.Following Katz et al. (1982, see also Caproni et al. 2013), the nutation angular frequency is equal to twice the difference of the orbital and precession angular frequencies, where ω p is negative because of the opposite direction with respect to the orbital motion.Using Eq. ( 2), the orbital period can be expressed as The putative supermassive black hole binary in M 81 has an orbital period of P orb ∼ 14 to 15 years for P p ≥ P n , which seems to be the case based on the long-term linear drift.The semimajor axis of the binary system can be estimated from the third Kepler law, which corresponds to a bin ∼ 8000R g (R g = GM tot /c 2 ) in gravitational radii if the primary supermassive black hole dominates the total mass.Using a primary mass fraction of x p = m 1 /M tot ∼ 0.9 (the secondary mass fraction is x s ∼ 0.1), the binary merger timescale is For an even more extreme ratio, x p ∼ 1, x s ∼ 10 −2 , we obtain τ merge ∼ 25 Gyr, and hence the system would generally be long-lived.

Conclusions
We present novel observations of the core region M 81*.Our observations are consistent with a precession of the jet of this galaxy, which was first proposed by Martí-Vidal et al. (2011).The jet precession period is roughly 7 years.On top of the precession (amplitude ∼7 • ), the jet exhibits a small linear drift of roughly 0.5 • /yr.However, we cannot confirm that the flare observed in the years 1998-2002 is connected to the precession of the jet, for instance, through Doppler boosting.Our sparse flux measurements do not show elevated flux in the years 2017-2019.Further, historic data do not show any flaring activity since the end of the last flare in 2002.We can rule out a self-similar modulation of the light curve with the period of the modulation.However, the sampling of the light curve is too spare to definitely rule out any modulation scenario.Thus, we consider it more likely that the flare observed in the years from 1998 to 2002 was not caused by variable Doppler boosting.Our observations are consistent with either a Lense-Thirring-induced precession of the jet or a binary-induced precession of the jet.They are therefore analogous to observations of other precessing jets, such as in 3C 345 and OJ 287, which both show a fast nutation-like component as well as a slower, but larger-amplitude, precessionlike variation in the jet position angle.If Lense-Thirring precession is responsible for the observed jet precession in all three active galactic nuclei, the underlying coupling of the accretion disk precession to the jet precession would show a remarkable self-similar coupling through vastly different accretion regimes.While both 3C 345 and OJ 287 accrete close to their Eddington limit, M 81* belongs to the class of radiatively inefficient accretion flows.This may hint that the jet physics responsible for the observed precession in these systems is accretion-flow-rate and accretion-flow-state independent.On the other hand, if a binary black hole is responsible for the apparent precession, then it may be the closest supermassive or intermediate-mass binary black hole candidate known to date.However, more data are necessary to confirm the precessing nature of M 81*, and in particular the proposed frequency-dependent time lag.

Notes.Fig. 1 .
Fig. 1.Comparison of the recent observation with the model extrapolation.Black dots show the position angle measurement as a function of time up to the year 2012.The orange pentagons show the 5 GHz data derived in this work.The black line shows the best-fit model to the data up to 2012 by Alberdi et al. (2013), excluding the new measurements, and extrapolated to the year 2021.The gray shaded region indicates the 1σ and 3σ contours.
-Vidal et al. (2011) proposed that the flare observed in the years 1998-2002 was caused by the variable Doppler boosting of the precessing jet as the viewing angle periodically changes in

Fig. 2 .Fig. 3 .
Fig. 2. Multifrequency data of the core position angle of M 81*.The black points indicate the 5 GHz position angle as presented by Martí-Vidal et al. (2011), and the orange pentagons show the new 2018 and 2019 measurements.The black line with gray contours (1σ, 2σ, and 3σ) shows the best fit to all 5 GHz data.The brown points show the 1.7 GHz data; the brown line shows the best fit to the 5 GHz data, shifted by 1.9 yr, as suggested byMartí-Vidal et al. (2011).The green points show the 8 GHz measurement; the green line shows the best-fit 5GHz model shifted by ∼−3.5 yr.

Fig. 4 .
Fig. 4. Posterior of the precession period, P p , and the nutation period, P n , as well as the precession cone half-angle, Ω p .

Table 1 .
Date, observation band, derived position angle, and reduced χ 2 of the observation analyzed in this study.

Table 2 .
Best-fit values of the sinusoidal and linear drift model derived from a χ 2 fit to the observed data.Parameter Sinusoidal model up to 2012 Sinusoidal model up to 2019 Linear drift model up to 2019 Table 2 reports the values derived in this manner.Figure 1 shows the best fit.The shaded regions show the 1σ and 3σ contours derived by sampling 1000 model realizations from the best-fit values and corresponding errors, where we assumed a Gaussian distribution and took the covariance of the best-fit values into account.The model is able to predict the observations in 2017 and 2019 well, with the predicted values within the 3σ contour.Figure A.1 shows the best-fit model when the new data are included; the best-fit parameters are reported in the third column of Table 2.