Constraining the Orbit of Supermassive Black Hole Binary 0402+379

TL;DR

This study uses multi-frequency VLBI observations over 12 years to detect relative motion in the supermassive black hole binary 0402+379, constraining its orbit and mass, and marking it as the first resolved black hole binary system.

Contribution

It provides the first direct measurement of relative motion in a resolved supermassive black hole binary, constraining its orbital parameters and mass.

Findings

Detected significant core motion at 0.0054c

Estimated orbital period of approximately 30,000 years

Constrained SMBHB mass to about 15 billion solar masses

Abstract

The radio galaxy 0402+379 is believed to host a supermassive black hole binary (SMBHB). The two compact core sources are separated by a projected distance of 7.3 pc, making it the most (spatially) compact resolved SMBHB known. We present new multi-frequency VLBI observations of 0402+379 at 5, 8, 15 and 22 GHz, and combine with previous observations spanning 12 years. A strong frequency dependent core shift is evident, which we use to infer magnetic fields near the jet base. After correcting for these shifts we detect significant relative motion of the two cores at $\beta=v/c=0.0054 \pm 0.0003$ at $PA= -34.4^\circ$. With some assumptions about the orbit, we use this measurement to constrain the orbital period $P\approx 3 \times 10^4$ y and SMBHB mass $M \approx 15 \times 10^9\ M_\odot$. While additional observations are needed to confirm this motion and obtain a precise orbit, this is apparently the first black hole system resolved as a visual binary.

Tables

Frequency	Date	Integration	BW	Polarization	IF	Reference
(GHz)		time (min)	(MHz)
4.98	01/24/2005	69	8	2	4	Rodriguez et al. (2006)
4.98	12/28/2009	286	32	4	4	This paper
4.98	06/20/2015	70	32	4	8	This paper
8.15	06/13/2005	69	8	2	4	Rodriguez et al. (2006)
8.15	12/28/2009	261	32	4	8	This paper
8.15	06/20/2015	261	32	4	8	This paper
15.35	03/02/2003	478	16	2	4	Maness et al. (2004)
15.35	01/24/2005	122	8	2	4	Rodriguez et al. (2006)
15.35	12/28/2009	292	32	4	8	This paper
15.35	06/20/2015	286	32	4	8	This paper
22.22	06/13/2005	251	8	2	4	Rodriguez et al. (2006)
22.22	12/28/2009	325	32	4	8	This paper
22.22	06/20/2015	334	32	4	8	This paper

Frequency	$a$(C1)	$b/a$(C1)	$\varphi$(C1)	$a$(C2)	$b/a$(C2)	$\varphi$(C2)
(GHz)	(mas)		(o)	(mas)		(o)
5	0.563	0.000	82.80	1.270	0.130	6.60
8	0.451	0.420	74.00	0.420	0.490	8.60
15	0.249	0.360	77.00	0.230	0.000	21.40
22	0.218	0.160	78.90	0.170	0.390	27.80

Epoch	Frequency	$S_{\nu }\left(C1\right)$	$S_{\nu }\left(C2\right)$	$r$	$\theta$
	(GHz)	(Jy)	(Jy)	(mas)	(o)
2005.07	5	0.057 $\pm$ 0.005	0.014 $\pm$ 0.001	6.942	-75.70
2009.99	5	0.058 $\pm$ 0.005	0.016 $\pm$ 0.001	6.841	-75.93
2015.43	5	0.060 $\pm$ 0.006	0.013 $\pm$ 0.001	6.884	-75.79
2005.45	8	0.067 $\pm$ 0.006	0.016 $\pm$ 0.002	6.913	-76.46
2009.99	8	0.052 $\pm$ 0.004	0.017 $\pm$ 0.002	6.920	-76.42
2015.43	8	0.083 $\pm$ 0.008	0.018 $\pm$ 0.002	6.913	-76.33
2003.17	15	0.070 $\pm$ 0.007	0.020 $\pm$ 0.002	6.929	-76.96
2005.07	15	0.054 $\pm$ 0.005	0.016 $\pm$ 0.002	6.959	-76.81
2009.99	15	0.029 $\pm$ 0.003	0.012 $\pm$ 0.001	6.985	-76.96
2015.43	15	0.058 $\pm$ 0.005	0.015 $\pm$ 0.001	6.956	-76.77
2005.45	22	0.037 $\pm$ 0.003	0.011 $\pm$ 0.001	6.950	-77.08
2009.99	22	0.020 $\pm$ 0.002	0.011 $\pm$ 0.001	6.984	-77.16
2015.43	22	0.040 $\pm$ 0.003	0.012 $\pm$ 0.001	6.969	-77.04

		Technique	${\chi }^2$
Parameters	Value	Bootstrap (95% )	$1\sigma$	$2\sigma$
$\Delta {\mathrm{RA}}_0$(mas)	-6.863	-6.892, -6.855	-6.859, -6.868	-6.858, -6.869
$\Delta {\mathrm{DEC}}_0$(mas)	1.474	1.448, 1.478	1.470, 1.478	1.468, 1.480
${\mu }_{RA}$($\mu$as/y)	-0.887	-0.970, -0.831	-1.245, -0.549	-1.389, -0.405
${\mu }_{DEC}$($\mu$as/y)	1.286	1.200, 1.401	0.878, 1.672	0.713, 1.836
a (mas)	0.756	0.700, 0.818	0.739, 0.777	0.731, 0.785
k	1.591	1.556, 1.823	1.565, 1.617	1.555, 1.628

Redshift	Luminosity	Half opening	Bulk Velocity	Inclination	Lorentz	Doppler Factor
$z$	Distance $D_L$ (Mpc)	angle $\Theta$ (∘)	factor ${\beta }_{app}$	angle $i$ (∘)	factor ${\Gamma }_j$	$\delta$
0.055	242.2 aafootnotemark:	4.29 bbfootnotemark:	0.4 ccfootnotemark:	71.3 ddfootnotemark:	1.077 eefootnotemark:	1.11 fffootnotemark:

Equations

$$B_F(h)=\frac{3.35\times 10^{-11}{D}_L[\mathrm{pc}]\delta \Delta \theta [\mathrm{mas}{]}^5\tan {\Theta }^2}{h[\mathrm{pc}]({\nu }_1^{-1}-{\nu }_2^{-1}{)}^5[(1+z)\sin i{]}^3{F}_v[\mathrm{Jy}{]}^2}$$

$$v=(v_x^2-v_x{v}_{y}x/y{)}^{1/2}$$

$$a=-(x^2-xy{v}_x/v_y{)}^{1/2}$$

$$\cos (i)=[-y{v}_y/(x{v}_x){]}^{1/2}$$

$$\theta =\pi +\mathrm{atan}\left([-y{v}_x/(x{v}_y){]}^{1/2}-[1-y{v}_x/(x{v}_y){]}^{1/2}\right)$$

Full text

1 Introduction

It is commonly believed that the later stages of galaxy evolution are governed by mergers. It is very common for galaxies to collide and interact with each other. Considering that most galaxies in the universe harbor supermassive black holes (SMBH) at their centers (Richstone et al., 1998), it can be inferred that massive black hole pairs should be the outcome of such mergers through the hierarchical formation of galaxies (Begelman et al., 1980). This implies that SuperMassive Black Hole Binaries (SMBHB) should be relatively common in the universe. However, despite very extensive searches, very few such systems have been observed (Burke-Spolaor, 2011; Tremblay et al., 2016). The reason for this could be that black holes in a binary system either merge quickly, or that one of them escapes the system (Merritt et al., 2005). Hence, understanding these SMBHB systems is important to understand a variety of processes ranging from galaxy evolution to active galactic nuclei (AGN) to black hole growth.

There are two types of galaxy mergers, major and minor. Major mergers result when the interacting galaxies are of similar sizes (mass ratio less than 3:1 (Stewart et al. (2009), and references therein)), whereas in the case of minor mergers one galaxy is significantly larger than the other. A crucial expectation related to galaxy mergers is the emission of gravitational waves. When galaxies merge, due to the dynamical friction between them, the black holes at their corresponding centers sink towards a common center. This leads to the formation of a binary system, such that its orbit decays due to the interaction between the stars, gas, and dust of both galaxies. The two black holes may reach a small enough separation that energy losses from gravitational waves allow the binary to coalesce into a single black hole (Begelman et al., 1980; Milosavljević & Merritt, 2003).

Numerous simulations have been performed to study these SMBH mergers. These simulations deal with various aspects such as black hole mass ratio, self or non-gravitating circum-binary discs, orbital spin, black hole spins, gas or stellar dynamics etc. (Barnes, 2002; Merritt et al., 2005; Escala et al., 2004, 2005; Sesana et al., 2006; Dotti et al., 2007; Callegari et al., 2009, 2011; Khan et al., 2011; Schnittman, 2013). In spite of several attempts, attaining the required resolution (last parsec problem $\sim 0.01$ pc) to study the black hole coalescence has been challenging. The fate of the merger at the parsec scale depends on the amount of surrounding stars and gas, and their interaction with the binary. In the case of a stellar background, due to 3-body scattering, formation of loss cone takes place (Sesana et al., 2007), whereas, in the case of gas rich mergers tidal forces inhibit the gas from falling onto the binary, hence creating a gap (Dotti et al. (2012), and references within). The loss cone causes a decay period longer than the Hubble time (Milosavljević & Merritt, 2003; Merritt et al., 2007), and hardening can be attained for a triaxial stellar remnant with the loss cone being replenished (Merritt et al., 2011) and an expected coalescence time of $\sim 10^{8}$ years (Khan et al., 2011). For gas-rich mergers, the gap doesn’t inhibit gas flow (Roedig et al., 2012) and a massive circumbinary disc around the binary promotes the decay leading to a timescale for an equal mass binary $\sim 10^{7}M_{\odot}$ that is less than the age of the Universe (Hayasaki, 2009). For higher mass black holes ($\sim 10^{8-9}M_{\odot}$) the timescales are greater. In a recent study by Khan et al. (2016), they report that for massive galaxies at high redshifts ($z>2$) it takes about few million years for black holes to coalesce once they form a binary, whereas, at lower redshifts where nuclear density of host is lower, it may take longer time of order a Gyr.

Gravitational waves from merging black holes are expected as a result of Einstein’s General Theory of Relativity (Einstein, 1916, 1997, 1918, 2002)111http://einsteinpapers.press.princeton.edu/. In September 2015, the Laser Interferometer Gravitational Wave Observatory (LIGO) discovered a gravitational wave (GW) source GW150914, and identified it as a merging binary black hole (BBH) (Abbott et al. 2016a). Although the masses of the two BHs are much smaller ($\thicksim$30 $M_{\odot}$) in comparison to SMBHBs (\thicksim$$10^{7}-10^{10}M_{\odot}), this discovery provides the first observational evidence for the existence of binary BH systems that inspiral and merge within the age of the universe. It motivates further studies of binary-BH formation astrophysics, and with the upcoming detectors such as evolving Laser Interferometer Space Antenna (eLISA, Amaro-Seoane et al. (2013)), it will be possible to detect low-frequency GW (around one mHz), emitted from the inspiral of massive black holes (Klein et al., 2016). While mergers of SMBHB’s are expected to be common emitters of GW radiation, modulating pulsar timing observations have not yet detected any evidence for a GW signal (Arzoumanian et al., 2016). Pulsar timing observations, unlike LIGO, should be more sensitive to SMBHB mergers (Shannon et al., 2015). 0402+379, with a separation of 7.3 pc between its core components, is one of the most important precursors of GW sources, and is important to understand the reason behind the low incidence of such systems. From the elliptical morphology of the 0402+379 host galaxy (Andrade-Santos et al., 2016), we believe this object to be the result of a major merger.

The radio galaxy 0402+379 was first observed by Xu et al. (1995) as a part of the first Caltech Jodrell Bank Survey (CJ1), although at that time it was not identified as a SMBHB. This source first acquired attention as a Compact Symmetric Object (CSO) candidate (small AGN with jets oriented close to the plane of the sky such that the radio emission from the jets is detected on both sides of the core) in 2003, in the full polarimetry analysis by Pollack et al. (2003). Subsequently, Maness et al. (2004) studied this source at multiple frequencies using the VLBA222The National Radio Astronomy Observatory is operated by Associated Universities, Inc., under cooperative agreement with the National Science Foundation. (Very Long Baseline Array), and on the basis of its properties, they classified it to be an unusual CSO. Rodriguez et al. (2006) studied this source in more detail and arrived at the conclusion that it is a SMBHB. This source contains two central, compact, flat spectrum and variable components (designated C1 and C2, see Fig. 1), a feature which has not been observed in any other compact source. This is one of the only spatially resolved SMBHB candidates (Gitti et al., 2013; Deane et al., 2014). The milliarcsecond scale separation requires high resolving power available only with a telescope such as the VLBA. Although other systems like RBS 797 and J1502+1115 (a triple system) with a separation of about 100 pc have been detected, no system other than 0402+379 has been resolved at parsec scales. We believe that this SMBHB is in the process of merging.

Rodriguez et al. (2006) imaged this source at multiple frequencies, studying the component motion at 5 GHz, but finding no significant detection of core displacement. In this paper, we incorporate new 2009 and 2015 epochs of 0402+379 VLBA observations at 5, 8, 15 and 22 GHz, while re-analyzing the 2003 and 2005 observations. These data show strong evidence for a frequency dependent core-shift effect (Lobanov, 1998; Sokolovsky et al., 2011). After accounting for this effect, our data set allows a detection of the relative motion of the two cores, making this the first visual SMBHB. We comment on the implications for the orbital period and masses. Throughout this discussion, we assume $H_{0}=71$ kms*-1Mpc-1* so that 1 mas = 1.06 pc.

2 Observations and Data Reduction

2.1 VLBA Observations

Observations were conducted on December 28, 2009 and June 20, 2015 with the VLBA at 4.98, 8.41, 15.35, and 22.22 GHz. For the 2015 observations, the total time on source was 70 min at 5 GHz, 260 min at 8 GHz, 290 min at 15 GHz, and 330 min at 22 GHz. 3C84 and 3C111 were observed for bandpass and gain calibration, respectively. The data recording rate was 2048 Mbps with two bit sampling. Each frequency was measured over eight intermediate frequencies (IFs) such that every IF consisted of a bandwidth of 32 MHz across 64 channels in both circular and their respective cross polarizations.

Standard data reduction steps including flagging, instrumental time delay, bandpass corrections, and frequency averaging were performed with the NRAO Astronomical Image Processing System (AIPS) (van Moorsel et al., 1996; Ulvestad et al., 2001). For all iterative self-calibration methods the initial model was a point source. Further cleaning, phase and amplitude self-calibration were executed manually using Difmap (Shepherd et al., 1995). The source structure was later model-fitted in the visibility $(u,v)$ plane with Difmap using circular and elliptical Gaussian components.

Fully calibrated VLBI archival data from the 2003 epoch (Maness et al., 2004) at 15 GHz, and the 2005 (Rodriguez et al., 2006) epochs at 5, 8, 15 and 22 GHz, have been included to study the core component motion with frequency and time. The visibilities were imaged and model fitted in Difmap, as with the 2009 and 2015 data, to obtain the core positions. The calibrator 3C111 has been observed in the same configuration across all four epochs. Details of the observations can be found in Table 1.

3 Measurement and Fits

The new 2015 data set provides the most sensitive measurement of the source geometry. Using model fitting in Difmap to the visibilities we determine the Gaussian size, axis ratio and position angle for each core component. These measurements are listed in Table 2. Additional Gaussian components are included in each model to account for the extended structure. We then follow Rodriguez et al. (2006) and Maness et al. (2004) in fixing these parameters in fits to the other epochs, while allowing only the positions and fluxes of C1 and C2 to vary. These are listed in Table 3, with the relative positions listed as angular separation $r$ and position angle $\theta$ measured north through east. The reported flux density errors combine the map rms $\sigma_{\rm rms}$ and an estimated systematic error in quadrature: $\sigma_{S}=[(0.1S_{\nu})^{2}+\sigma_{\rm rms}]^{1/2}$.

The effective position errors are more difficult to estimate. The statistical errors $\approx a/2(S/N)$ (Fomalont, 1999) are very small ($\sim 2\ \mu$as), and systematic effects certainly dominate. We made an initial check on these errors, by re-doing the previous Stokes I map analysis in Stokes LL and RR. 0402+379 is an unpolarized source so we expect data from both Stokes RR and LL to yield similar distance measurements. Another way to obtain these error estimates is to split the data in time or frequency. However, these maps would have a different (u,v) coverage, thus making it difficult to make a comparison between them. If the polarization-dependent structure differences are small as expected, then any measured differences can be attributed to systematic errors. This technique has been discussed previously by Roberts et al. (1991) and McGary et al. (2001). Decomposing the relative positions into RA($x$) and DEC($y$), we find that the median shifts in the core centroid positions are $\sigma_{x}=7\ \mu$as and $\sigma_{y}=8\ \mu$as for the higher frequencies, and $\sigma_{x}=31\ \mu$as and $\sigma_{y}=34\ \mu$as at 5 GHz.

3.1 Analysis

The core component flux density arises from the surface at which the self-absorption optical depth is unity (Blandford et al., 1979). Since this is strongly frequency dependent, we expect an asymmetric extended structure, such as the jet base which defines the core, to have a frequency-dependent centroid. This effect is very obvious in the raw positions (Fig. 2), where the lower frequency centroids are shifted to the NE along the larger scale outflow position angle (Fig. 1a). Similar shifts have been detected in a number of AGN (Lobanov, 1998; Sokolovsky et al., 2011). Since we are measuring position relative to C1, any extension to that source may also contribute to the relative core shift; this need not be at the same position angle. However the combined shift appears to be dominated by C2 and we indeed find that the frequency-dependent shift is along the $47^{\circ}$ (N-E) position angle of the C2 outflow. According to Lobanov (1998) the shift can be parameterized as $r_{c}=a\nu^{-1/k}$, where $a$ is the shift amplitude and $k$ depends on the jet geometry, particle distribution and magnetic field. For example a conical jet with a synchrotron self-absorbed spectrum gives $k=1$ (Lobanov, 1998).

By correcting to an infinite frequency we can mitigate the effect of this core shift on the position of the nucleus (Fig. 3). We are of course especially interested in the relative motion of C1 and C2 and so our model includes a fiducial relative position (at epoch 2000.0) as well as relative proper motion in RA and DEC. Thus our model has six fit parameters (if we include the core shift position angle as a fit parameter, we do indeed obtain $46\pm 1^{\circ}$, but prefer to fix this via the larger scale jet axis). Our data set are the 13 $r$, $\theta$ ($x$, $y$) position pairs over four epochs and four frequencies, so the fit has $26-6=20$ degrees of freedom (DoF).

Using our measurements and the position error estimates above, we performed a $\chi^{2}$ minimization to determine the model parameters. These are listed in Table 4. The parameter error estimates are somewhat subtle. Since the fit minimum has $\chi^{2}$/DoF = 2.78, we must have systematic errors beyond those estimated above. The conventional approach is to uniformly inflate all errors until the effective $\chi^{2}$/DoF=1. This is equivalent to estimating errors using the $\chi^{2}$ surface with increases of $+1,\,+2\,...\times\chi^{2}$/DoF. We list these “1 $\sigma$” and “2 $\sigma$” confidence intervals in Table 4.

However, this uniform inflation assumes that all errors have a Gaussian distribution and are equally underestimated. This is unlikely to be true. An alternative approach is to estimate errors via a bootstrap analysis (Efron 1987). This has the virtue of using only the actual data values (not the error estimates), but does pre-suppose that the observed data values are an unbiased draw from an (unknown) error distribution about the true values. Although our set of 13 position pairs is somewhat small for a robust bootstrap, we generated 10,000 re-sampled realizations of the data set, replacing five pairs in each realization with random draws from the remaining pairs. Each realization was subject to the full least-squares fit for all model parameters. The histograms of the fit values for the individual parameters were used to extract 95% confidence intervals for each quatity. These confidence intervals are listed in Table 4. They accord fairly well with the inflated $\chi^{2}$ estimates.

In general, the parameters appear well-constrained. The core-shift coefficients $a$ and $k$ are estimated to $\sim$3 % accuracy (Table 4). The epoch position range is somewhat larger in the bootstrap error analysis, evidently as a result of the substantial offset of the 2009.9 position from the general trend.

The coefficient $k$ depends on the shape of electron energy spectrum, magnetic strength and particle density distribution (Lobanov, 1998). If $k$=1, it implies that the jet has a conical shape, where synchrotron self-absorption is the dominant absorption mechanism. We have obtained $k=1.591\pm.232$ (via the bootstrap technique), in accordance with the literature (Sokolovsky et al. (2011)) where the highest reported $k$ value is $\thicksim$ 1.5. This implies that our observations are consistent with the synchrotron self-absorption mechanism.

The core-shift depends mainly on the frequency as well as the magnetic field and spectral index (Lobanov, 1998). All our measurements are derived from the same frequencies over time, however, there is a possibility of variation with magnetic field and spectral index. We calculated the spectral index for three epochs (2005, 2009 and 2015) using the peak intensities, and we have found these values to range from $-0.58$ to $-0.98$ and $-0.43$ to $-0.50$ for C1 and C2 respectively. From these ranges of spectral index, we can see that the variation over 10 years is quite small ($\thicksim 0.4$). To our knowledge, no time varying core-shift offsets have been reported in the literature. For our analysis, we assume that the core-shift is constant over time.

In Fig. 3 we plot the relative C2 core position, shifted to infinite frequency according to the best-fit $a$ and $k$, for each of the 13 observation frequencies and epochs. The plotted error ellipses for the 8, 15 and 22 GHz observations are the formal $\sigma_{x}$ and $\sigma_{y}$ from the polarization analysis. There are large outliers, especially the 2009.9 epoch. However, the overall shift is quite significant with motion detected in both RA and DEC, at $>3\ \sigma$ in the $\chi^{2}$ analysis and at well over $95$% confidence in the bootstrap analysis. Additional epochs, especially at high frequency will, however, be needed to make the motion visually clear.

To compare these above results, we also studied the motion of jet components. We find that the bright southern jet components continue to move away from the core, consistent with the previous results (Rodriguez et al., 2006). For the weaker northern hotspot, the agreement is not as good, as it seems to exhibit inconsistent motion for some frequencies. However, this may be the result of errors in the previous measurements based on 5 GHz observations.

Using our best-fit $\mu_{RA}$ and $\mu_{Dec}$, we shift the raw data points (Table 3) to epoch 2000.0. For each frequency, we obtain an average relative RA and DEC, which are subsequently subtracted from the fiducial 2000.0 point (Table 4) to obtain the distance from the core ($r_{c}$). This has been plotted to demonstrate our fitted frequency dependence, shown in Fig. 4. The plotted errors have been obtained by error propagation using the above stated errors ($\sigma_{x}$ and $\sigma_{y}$).

3.2 Magnetic Field Estimate

The core-shift effect is useful to deduce various jet related physical parameters, including the magnetic field strength. Lobanov (1998) and Hirotani (2005), for example provide a derivation that assumes equipartition between the particle and magnetic field energy densities in the jet. An alternate formulation by Zdziarski et al. (2015) avoids the equipartition assumption, using the flux density $F_{v}$ at a jet axis distance $h$ to estimate the magnetic field strength as

[TABLE]

where $z$ is redshift, $D_{L}$ is the luminosity distance in pc, $\delta=[\Gamma_{j}(1\ -\ \beta_{j}cos\ i)]^{-1}$ is the Doppler factor, $\Gamma_{j}$ is the minimum Lorentz factor, $\beta_{j}$ is the jet bulk velocity factor (obtained from Rodriguez et al. (2009)), $i$ is the inclination angle, $\Delta\theta$ is the observed angular core shift (Lobanov 1998), and $\Theta=$ arctan $\frac{\sqrt{d^{2}-b_{\phi}^{2}}}{2r}$ is the jet half opening angle (Pushkarev et al., 2012). From the latest 8 GHz map, we have obtained the full width at half maximum of a Gaussian fitted to the transverse jet brightness component, $d=4.130\pm 0.017$ mas (minor axis of jet); the beam size along the jet direction, $b_{\phi}=1.26$ mas; and the distance to the core along the jet axis, $r=26.320\pm 0.017$ mas. Since the extended jet components are not readily detected at 15 and 22 GHz, we assume the same opening angle for all frequencies. Table 5 gives our estimated values for these parameters, with the origins in the footnotes. The numerical constant in the above equation has been obtained for $p=2$, where $p$ is index of the electron power law (see Zdziarski et al. (2012, 2015)).

From a weighted linear fit of $\delta\theta$ against $\nu_{1}^{-1}-\nu_{2}^{-1}$, we obtain a slope = $1.128\pm 0.152$ and intercept = $0.008\pm 0.010$. Instead of calculating $\frac{\delta\theta}{\nu_{1}^{-1}-\nu_{2}^{-1}}$ for each frequency separately, its slope has been used in calculating the magnetic field strength. Our fit indicates a magnetic field strength $0.71\pm 0.25$ G at $h=1$ pc, similar to that for other jets (O’Sullivan el al., 2009).

3.3 Orbital Models

Our measured proper motion $\mu_{\rm RA}=-0.89\pm 0.07\,\mu$as/y, $\mu_{\rm Dec}=1.29\pm 0.10\,\mu$as/y (symmetric width of the 95 % CL bootstrap range) corresponds to a proper motion $\mu$ of $1.57\pm 0.08\,\mu$as/y at $PA_{\mu}=-34.6\pm 2.9^{\circ}$ (if we use the “1 $\sigma$” $\chi^{2}$ errors, the amplitude uncertainty is $\pm 0.38\,\mu$as/y). Thus this is at least a $4\ \sigma$ detection. It is consistent with the non-detection of a proper motion in Rodriguez et al. (2006), where 15 years of 5 GHz data (1990-2005) were used to estimate $\mu=6.7\pm 9.4\,\mu$as/y; our higher frequency data and core-shift correction are essential for measuring the much smaller motion.

We now ask if this proper motion is consistent with a shift due to the relative orbits of the two BH. At $z=0.055$, it corresponds to a projected space velocity of $\beta=v/c=0.0054\pm 0.0003$, so a Keplerian analysis suffices. First, the ratio $2\pi r/\mu=2\pi\times 7.02\ {\rm mas}/0.00157$ mas/y = 28,000 y gives a characteristic orbital timescale. Thus over our 12 y baseline, the core position PA has rotated by less than a degree. This does not allow us to fit for precise orbital parameters. In particular, we have four measurements from the VLBI analysis (relative position and proper motion) while we need six parameters to define the relative orbit.

We note that the above derived $\sim 28000$ y period is rather close to the Earth’s spin axis precession period of $\sim 26000$ y, we believe this to be a coincidence. The differential astrometry performed here should not be expected by precession as that term has been removed with the correlator model and affects both the sources in an identical way.

If we assume circular motion ($e=0$), then we can determine the relative orbit in terms of one additional free parameter. In practice it is easiest to select the PA of the projected orbit normal (measured N through E) and then resolve the relative positions $x$ and $y$ (in mas) and relative velocities $v_{x}$ and $v_{y}$ (in mas/y) in this rotated coordinate system. Then the orbit parameters are

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $a$ and $v$ are the relative orbit radius and velocity, $i$ is the inclination and $\theta$ gives the phase at our observation epoch. In fact it is more interesting to plot the total mass $M=v^{2}a/G$ and period $P=2\pi a/v$ against the orbit inclination $i$ (Fig. 5). Note that with our assumption of a circular orbit only fairly large inclinations are consistent with our C1-C2 offset and relative motion (Fig. 3). Typical orbital periods are indeed 20-30ky, but the masses required by our apparent velocity are quite large $\geq 15\times 10^{9}\ M_{\odot}$. With our nominal fit errors, the minimum mass is $M_{9}=15.4\pm 1.3$. If one relaxes the $e=0$ assumption, smaller masses are allowed, but then the solution is nearly unconstrained.

The orbital eccentricity grows as the orbit shrinks since both stars and gas extract energy and angular momentum from the binary. For a stellar background, it depends on the mass ratio, with equal mass binaries producing orbits that are usually circular or slightly eccentric with $e<0.2$ (Merritt et al., 2007). If a pair during the binary formation starts out with a non-zero eccentricity it may never become circular instead it tends to become more eccentric (Matsubayashi et al., 2011). In the case of gas driven mergers, it depends on the disc thickness and the SMBBH’s location inside the disc. The critical value of e is reported to be $\sim 0.6$, such that system with high eccentricities tend to shrink to this value (Armitage & Natarajan, 2005; Cuadra et al., 2009; Roedig et al., 2012). In the case of 0402+379, it has been found to be embedded in cluster gas (Andrade-Santos et al., 2016), which makes it likely to have a non-zero eccentricity.

In Rodriguez et al. (2009) HI absorption measurements were used to infer kinematic motion about an axis inclined $\sim 75^{\circ}$ to the Earth line of sight, passing through C2, the origin of the kpc-scale jets. If we look at the solution derived here we see that the $PA=47^{\circ}$ axis corresponds to $i=71.3^{\circ}$ (Fig. 5 & Fig. 6 ), in reasonable agreement with the HI estimate. The binary spin can be different from the orbital angular momentum, however, it’s been found that if the amount of gas accreted is high (1-10% of the black hole) on the timescales of binary evolution, it can change according to the orbital axis (Schnittman (2013), and references therein). Binary orbital axis and individual black hole spins tend to realign due to interaction with external gas except when the mass ratios are extreme ($>>1$) whereas, torques from stars can cause misalignments of the binary orbit orientation from the disc (Coleman Miller and Krolik, 2013). For this fit we have $P=49$ ky and $M=16.5\times 10^{9}\ M_{\odot}$.

Because we find a large proper motion $\mu$ we expect the orbital motion to induce a substantial radial velocity in the relative orbit. Some values are given in Fig. 6 and for $PA=47^{\circ}$ we expect a relative $v_{r}=700$ km/s. While the HI measurements do show velocity differences of this order, we do not see such large velocities in the optical line peaks. Examining the Keck spectra in Romani et al. (2014) we see that the stellar features of the elliptical host center on $16,618\pm 53$ km/s, while the Seyfert I-type narrow-line core emission centers on $16,490$ km/s. Narrow line emission extends several arcsec from the core spanning $\sim 300$ km/s while in the unresolved kpc core the velocity dispersion is $\sim 750$ km/s. Thus, while at least $2\times 10^{10}\ M_{\odot}$ lies within the central kpc, we do not see multiple components shifted by $>500$ km/s. However, the full line width does accommodate such velocities and the wings of the H$\alpha$ complex are centered at $\sim 17,020$ km/s suggesting that fainter broad line emission might include components spread over $>1000$ km/s. Further $v_{r}$ above is the relative velocity; if only the heavier component has bright optical emission, then the broad line velocity shift from the background galactic velocity (arguably near the center of mass velocity) will be reduced to $mv_{r}/M_{Tot}$. Indeed, Rodriguez et al. (2009) assume that the jet-producing C2 core is the dominant mass and the center of rotation. If this core also dominates the broad line emission, then we expect that our VLBI relative velocity is dominated by the motion of C1 and the optical radial velocity shift from the host velocity may be small.

We must also compare with other core mass estimates. As noted above the optical lines indicate several $\times 10^{10}\ M_{\odot}$ in the central kpc. HI absorption velocities require $>7\times 10^{8}\ M_{\odot}$ (Rodriguez et al., 2009). And finally the host bulge luminosity also indicates a large hole mass $M_{\bullet}\sim 3\times 10^{9}\ M_{\odot}$ (Romani et al., 2014). All of these suggest substantial hole mass. The very large $M_{9}\sim 15$ masses implied by our fits are not excluded but do stretch the available mass budget.

3.4 Comments on the resolved SMBHB Population

The process of hierarchical merging should make close SMBHB common, but to date few candidates at sub-kpc separations have been seen. The resolved (massive) SMBHB seem to be preferentially in galaxy clusters or their products. For example the SMBHB candidate RBS 797 (Gitti et al., 2013) resides in a cool-core cluster at z=0.35. 0402+379 itself lies in a massive galaxy and dense X-ray halo (likely a fossil cluster) at z = 0.055. So such environments seem to be a good place to look for additional systems. Another path to discovering multi-BH nuclei has been described by Deane et al. (2014) who find J1502+1115 to be a triple system, with a closest pair separation of 140 pc at redshift z=0.39. Compact radio jets in the closest pair of this source exhibit rotationally symmetric helical structure, plausibly due to binary-induced jet precession. However the total number of resolved compact cores at pc scales seems very small with 0402+379 remaining the only clear example out of several thousand mapped sources (Burke-Spolaor, 2011; Tremblay et al., 2016).

Of course systems of even smaller separation are of the greatest interest since at $r\sim 0.01$ pc losses from gravitational radiation will dominate and the merging binaries can be an important signal in pulsar timing studies (Ravi et al., 2015). At present, we rely on arguments about evolution of the wider systems to infer the existence of merging SMBHB. If such evolution occurs we might hope for a discovery of an intermediate $r\sim 0.1$ pc scale massive $>10^{9}\ M_{\odot}$ system at low $z$ where sufficient resolution for a kinematic binary study is possible with high frequency VLBI. Such a binary would have $P<10^{3}$ y and a well-constrained visual orbit should be achievable, making possible a precision test of the SMBHB nature (Taylor, 2014). However, we should note that our study of the galactic halo of 0402+379 (Andrade-Santos et al., 2016) implies that it has stalled at its present 7 pc separation for several Gyr, so the path between resolvable and gravitational radiation-dominated SMBHB may not always be smooth.

4 Conclusion

In this study of 0402+379, we have focused on two aspects: frequency dependent core-shift and secular relative core motion. Both effects are observed, but the measured values present interpretation challenges.

The strong observed core-shift matches well with the large-scale jet axis. It also provides quite typical estimates for the jet base magnetic field of $\sim 0.45-0.95$  G. However the core-shift index 1.591 (1.556-1.823 95% CL range) is somewhat large (expected $k\sim 1$), with the highest reported value in literature is $k\thicksim$ 1.5 (Sokolovsky et al., 2011). This may be an artifact of our constant core-shift fit assumption, as perturbations could arise from the epoch-to-epoch variation in the underlying core component fluxes.

After accounting for this core shift the infinite frequency relative positions of the C1 and C2 cores undergo a statistically significant secular proper motion. The motion corresponds to $\beta=0.0054\pm 0.0003$ and, if orbital, it represents the first direct detection of orbital motion in a SMBHB, and promotes this system to a visual binary. Although we do not have sufficient observables to solve for an orbit, we can find plausible orbits, even assuming $e=0$. Intriguingly such orbits align well with the large scale jet axis and have similar inclination to those estimated with HI absorption VLBI. But the required masses are quite large (highest reported mass is 21 billion $M_{\odot}$, McConnell et al. (2011)). To test our orbital picture, additional VLBI epochs to confirm the consistency of the proper motion will be essential, and further studies of the core dynamics, especially at sub-kpc scales will also be very important. We should not forget that including a finite orbital eccentricity can allow smaller masses, but we will need additional kinematic constraints to motivate such solutions.

Thus discovery of possible orbital motion in 0402+379 presents the exciting prospect of probing a SMBHB’s kinematics. Certainly, extensions to our high frequency VLBI campaign can improve the measurements, but this proper motion is perhaps the most exciting as a spur to searches for tighter, faster and more easily measured examples of resolved SMBHBs.

Acknowledgements:

This research has made use of NASA’s Astrophysics Data System. English translations of Einstein (1916) and Einstein (1918) made available via The Digital Einstein Papers of the Princeton University Press aided the preparation of this paper. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.. Partial support for this work was provided by the National Aeronautics and Space Administration through Chandra Award Number GO4-15121X issued by the Chandra X-ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of the National Aeronautics Space Administration under contract NAS8-03060.