Constraining the Orbit of the Supermassive Black Hole Binary 0402+379

Observations were conducted on 2009 December 28 and 2015 June 20 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 minutes at 5 GHz, 260 minutes at 8 GHz, 290 minutes at 15 GHz, and 330 minutes 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 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) epoch 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 was done for 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.

The core-component flux density arises from the surface at which the self-absorption optical depth is unity (Blandford & Königl 1979). Since this is strongly frequency-dependent, we expect an asymmetric extended structure, such as the jet base that defines the core, to have a frequency-dependent centroid. This effect is very obvious in the raw positions (Figure 2), where the lower-frequency centroids are shifted to the northeast along the larger-scale outflow position angle (Figure 1(a)). Similar shifts have been detected in a number of AGNs (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° (north to east) 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).

Zoom In Zoom Out Reset image size

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

Zoom In Zoom Out Reset image size

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σ" and "2σ" confidence intervals in Table 4.

Table 4.  Fitting Parameters

		Technique	${\chi }^{2}$
Parameters	Value	Bootstrap (95% )	$1\sigma$	$2\sigma$
${\rm{\Delta }}{\rm{R}}.{\rm{A}}{.}_{0}$(mas)	−6.863	−6.892, −6.855	−6.859, −6.868	−6.858, −6.869
${\rm{\Delta }}\mathrm{decl}{.}_{0}$(mas)	1.474	1.448, 1.478	1.470, 1.478	1.468, 1.480
${\mu }_{{\rm{R}}.{\rm{A}}.}$(μas yr−1)	−0.887	−0.970, −0.831	−1.245, −0.549	−1.389, −0.405
${\mu }_{\mathrm{decl}.}$(μas yr−1)	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

Note. Fitted parameter values (Column 2) and their corresponding confidence intervals obtained from two different techniques: bootstrap analysis (Column 3) and ${\chi }^{2}$ minimization (Column 5 and 6). ${\rm{\Delta }}{\rm{R}}.{\rm{A}}{.}_{0}$ and ${\rm{\Delta }}\mathrm{decl}{.}_{0}$ are infinite-frequency core offsets at epoch 2000.0; ${\mu }_{{\rm{R}}.{\rm{A}}.}$ and ${\mu }_{\mathrm{decl}.}$ are proper motion estimates; a and k are core shift fitting parameters (${r}_{c}=a\ {\nu }^{(-1/k)}$).

Download table as:  ASCIITypeset image

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 quantity. 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 ∼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 the 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 $\sim$ 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 ($\sim 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 Figure 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 R.A. and decl., at $\gt 3\sigma$ in the ${\chi }^{2}$ analysis and at well over 95% confidence in the bootstrap analysis. However, additional epochs, especially at high frequency, will be needed to make the motion visually clear.

To compare these above results, we also studied the motions 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 }_{{\rm{R}}.{\rm{A}}.}$ and ${\mu }_{\mathrm{decl}.}$, we shift the raw data points (Table 3) to epoch 2000.0. For each frequency, we obtain an average relative R.A. and decl., 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 Figure 4. The plotted errors have been obtained by error propagation using the above stated errors (${\sigma }_{x}$ and ${\sigma }_{y}$).

Zoom In Zoom Out Reset image size

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, provided 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

$$\begin{eqnarray}&&{B}_{F}(h)=\displaystyle \frac{3.35\times {10}^{-11}\ {D}_{L}[\mathrm{pc}]\ \delta \ {\rm{\Delta }}\theta {[\mathrm{mas}]}^{5}\ \tan {{\rm{\Theta }}}^{2}}{h[\mathrm{pc}]\ {({\nu }_{1}^{-1}-{\nu }_{2}^{-1})}^{5}\ {[(1+z)\sin i]}^{3}\ {F}_{v}{[\mathrm{Jy}]}^{2}},\end{eqnarray} \tag{ 1 }$$

where z is redshift, D_L is the luminosity distance in parsecs, $\delta ={[{{\rm{\Gamma }}}_{j}(1-{\beta }_{j}\cos i)]}^{-1}$ is the Doppler factor, ${{\rm{\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, ${\rm{\Delta }}\theta$ is the observed angular core shift (Lobanov 1998), and ${\rm{\Theta }}\,=$ arctan $\tfrac{\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 FWHM 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 the index of the electron power law (see Zdziarski et al. 2012, 2015).

Table 5.  C2 Jet Parameters

Redshift	Luminosity	Half-opening	Bulk Velocity	Inclination	Lorentz	Doppler Factor
z	Distance DL (Mpc)	Angle Θ (${}^{^\circ }$)	Factor ${\beta }_{\mathrm{app}}$	Angle i (${}^{^\circ }$)	Factor ${{\rm{\Gamma }}}_{j}$	δ
0.055	242.2a	4.29b	0.4c	71.3d	1.077e	1.11f

Notes.

^aLuminosity distance was obtained for the cosmological model : ${H}_{0}=71$ km s⁻¹ Mpc⁻¹, ${{\rm{\Omega }}}_{{\rm{\Lambda }}}$ = 0.73, ${{\rm{\Omega }}}_{M}$ = 0.27. ^bPushkarev et al. (2012). ^cRodriguez et al. (2009). ^dSection 3.1. ^eLorentz factor, ${{\rm{\Gamma }}}_{j}={(1+{\beta }_{\mathrm{app}}^{2})}^{\tfrac{1}{2}}$ (Zdziarski et al. 2015). ^fDoppler factor, $\delta ={[{{\rm{\Gamma }}}_{j}(1-{\beta }_{j}\cos i)]}^{-1}$ (Zdziarski et al. 2015).

Download table as:  ASCIITypeset image

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

Our measured proper motion ${\mu }_{{\rm{R}}.{\rm{A}}.}=-0.89\pm 0.07\,\mu \mathrm{as}\,{\mathrm{yr}}^{-1}$, ${\mu }_{\mathrm{decl}.}=1.29\pm 0.10\,\mu \mathrm{as}\,{\mathrm{yr}}^{-1}$ (symmetric width of the 95% CL bootstrap range) corresponds to a proper motion μ of $1.57\pm 0.08\,\mu \mathrm{as}\,{\mathrm{yr}}^{-1}$ at ${\mathrm{PA}}_{\mu }=-34\buildrel{\circ}\over{.} 6\pm 2\buildrel{\circ}\over{.} 9$ (if we use the "1σ" ${\chi }^{2}$ errors, the amplitude uncertainty is $\pm 0.38\,\mu \mathrm{as}\,{\mathrm{yr}}^{-1}$). Thus, this is at least a 4σ 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 \mathrm{as}\,{\mathrm{yr}}^{-1}$; 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 BHs. 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\ \mathrm{mas}/0.00157$ mas yr⁻¹ = 28,000 year gives a characteristic orbital timescale. Thus, over our 12-year 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 ∼28,000-year period is rather close to the Earth's spin axis precession period of ∼26,000 years, but 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 yr⁻¹) in this rotated coordinate system. Then the orbit parameters are

$$\begin{eqnarray*}&&v={({v}_{x}^{2}-{v}_{x}{v}_{y}x/y)}^{1/2},\end{eqnarray*}$$

$$\begin{eqnarray*}&&a=-{({x}^{2}-xy{v}_{x}/{v}_{y})}^{1/2},\end{eqnarray*}$$

$$\begin{eqnarray*}&&\cos (i)={[-y{v}_{y}/(x{v}_{x})]}^{1/2},\end{eqnarray*}$$

$$\begin{eqnarray*}&&\theta =\pi +a\tan ({[-y{v}_{x}/(x{v}_{y})]}^{1/2}-{[1-y{v}_{x}/(x{v}_{y})]}^{1/2}),\end{eqnarray*}$$

where a and v are the relative orbit radius and velocity, i is the inclination, and θ 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 (Figure 5). Note that with our assumption of a circular orbit, only fairly large inclinations are consistent with our C1–C2 offset and relative motion (Figure 3). Typical orbital periods are indeed 20–30 ky, but the masses required by our apparent velocity are quite large $\geqslant 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.

Zoom In Zoom Out Reset image size

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\lt 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. 2007). In the case of gas-driven mergers, it depends on the disk thickness and the SMBBHs location inside the disk. The critical value of e is reported to be ∼0.6, such that systems 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) H i 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 kiloparsec-scale jets. If we look at the solution derived here we see that the $\mathrm{PA}=47^\circ$ axis corresponds to $i=71.3$° (Figures 5 and 6), in reasonable agreement with the H i estimate. The binary spin can be different from the orbital angular momentum; however, it has been found that if the amount of gas accreted is high (1%–10% of the black hole), on  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 for when the mass ratios are extreme ($\gg 1$), whereas torques from stars can cause misalignments of the binary orbit orientation from the disk (Coleman Miller & Krolik 2013). For this fit we have P = 49 ky and $M=16.5\times {10}^{9}\,{M}_{\odot }$.

Zoom In Zoom Out Reset image size

Because we find a large proper motion μ we expect the orbital motion to induce a substantial radial velocity in the relative orbit. Some values are given in Figure 6 and for $\mathrm{PA}=47^\circ$ we expect a relative v_r = 700 km s⁻¹. While the H i 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 arcseconds from the core spanning ∼300 km s⁻¹,while in the unresolved kiloparsec core the velocity dispersion is ∼750 km s⁻¹. Thus, while at least $2\times {10}^{10}\,{M}_{\odot }$ lies within the central kiloparsecs, we do not see multiple components shifted by $\gt 500$ km s⁻¹. However, the full line width does accommodate such velocities and the wings of the Hα complex are centered at ∼17,020 km s⁻¹, suggesting that fainter broad-line emission might include components spread over $\gt 1000$ km s⁻¹. Furthermore, 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}_{\mathrm{Tot}}$. Indeed, Rodriguez et al. (2009) assumed 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 kiloparsecs. H i absorption velocities require $\gt 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.

The process of hierarchical merging should make close SMBHB common, but to date few candidates at sub-kiloparsec 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. Thus, such environments seem to be a good place to look for additional systems. Another path to discovering multi-BH nuclei was described by Deane et al. (2014), who found 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 parsec scales seems very small, with 0402+379 remaining as the only clear example out of several thousand mapped sources (Burke-Spolaor 2011; Tremblay et al. 2016).

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 for pulsar-timing studies (Ravi et al. 2015). At present, we rely on arguments about the evolutions of the wider systems to infer the existence of merging SMBHBs. If such evolutions occur we might hope for a discovery of an intermediate $r\sim 0.1$ pc-scale massive $\gt {10}^{9}\,{M}_{\odot }$ system at low z, where a sufficient resolution for a kinematic binary study is possible with high-frequency VLBI. Such a binary would have $P\lt {10}^{3}$ yr and a well-constrained visual orbit would be achievable, making possible a precision test of its 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 SMBHBs may not always be smooth.