A Large Systematic Search for Close Supermassive Binary and Rapidly Recoiling Black Holes - III. Radial Velocity Variations

TL;DR

This study monitors 88 quasars over up to 12 years to identify supermassive black hole binaries through radial velocity variations, finding a few promising candidates and establishing mass limits to evaluate the binary hypothesis.

Contribution

It introduces a long-term spectroscopic monitoring method and a new approach to constrain the mass and period of potential supermassive black hole binaries.

Findings

Three quasars show systematic velocity changes suggesting binary motion.

Lower mass limits for candidates range from 47,000 to 380 million solar masses.

Most candidates do not yet rule out the binary hypothesis due to observational constraints.

Abstract

We have been spectroscopically monitoring 88 quasars selected to have broad H$\beta$ emission lines offset from their systemic redshift by thousands of km s$^{-1}$. By analogy with single-lined spectroscopic binary stars, we consider these quasars to be candidates for hosting supermassive black hole binaries (SBHBs). In this work we present new radial velocity measurements, typically 3-4 per object over a time period of up to 12 years in the observer's frame. In 29/88 of the SBHB candidates no variability of the shape of the broad H$\beta$ profile is observed, which allows us to make reliable measurements of radial velocity changes. Among these, we identify three objects that have displayed systematic and monotonic velocity changes by several hundred km s$^{-1}$ and are prime targets for further monitoring. Because the periods of the hypothetical binaries are expected to be long, we cannot hope to observe many orbital cycles during our lifetimes. Instead, we seek to evaluate the credentials of the SBHB candidates by attempting to rule out the SBHB hypothesis. In this spirit, we present a method for placing a lower limit on the period, and thus the mass, of the SBHBs under the assumption that the velocity changes we observe are due to orbital motion. Given the duration of our monitoring campaign and the uncertainties in the radial velocities, we were able to place a lower limit on the total mass in the range $4.7\times10^4-3.8\times10^8$ $M_{\scriptscriptstyle \odot}$, which does not yet allow us to rule out the SBHB hypothesis for any candidates.

Tables

Table 1. Table 1: Results of Cross-correlation Analysis of Broad H β 𝛽 \beta Lines (Abridged) †

			H$\beta$ Offset in	H$\beta$ Shift Relative
Object	Date	$\Delta t_{\mathrm{rest}}$111Time intervals are given in the rest frame.	SDSS Spectrum	to SDSS Spectrum222Numbers without brackets denote a statistically significant shift, while brackets indicate a shift measurement that is consistent with zero at 99% confidence. A positive shift means that the post-SDSS spectrum is shifted towards longer wavelengths.	$2.6\sigma$ Error333Uncertainties on the broad H$\beta$ shift at 90% and 68% confidence ($2.6\sigma$ and $1\sigma$, respectively). These are used in the test described in Section 4 to determine minimum periods that are consistent with the radial velocity curves.	$1\sigma$ Error333Uncertainties on the broad H$\beta$ shift at 90% and 68% confidence ($2.6\sigma$ and $1\sigma$, respectively). These are used in the test described in Section 4 to determine minimum periods that are consistent with the radial velocity curves.
(SDSS J)	(UT)	(years)	(km s-1)	(km s-1)	(km s-1)	(km s-1)
001224	2001/08/20		$-$1870$\pm$800
	2009/12/16	07.2		$-120$	$+50/-40$	$\pm 20$
	2010/10/12	07.9		$-170$	$+30/-40$	$\pm 20$
	2014/08/28	11.3		$-120$	$+40/-70$	$+30/-40$
002444	2000/12/22		00$-$700$\pm$100
	2009/12/18	07.5		$\left(100\right)$	$\pm 200$	$+100/-90$
	2014/08/29	11.5				ddfootnotemark:
015530	2001/09/16		1400$\pm$100
	2009/12/17	07.7		$\left(100\right)$	$\pm 200$	$\pm 100$
	2011/08/31	09.2		$\left(-170\right)$eefootnotemark:	$+100/-80$	$\pm 40$
	2014/08/28	12.0		$-390$	$+60/-100$	$+20/-30$

Table 2. Table 2: Lower limits on binary orbital parameters 1 1 1 Orbital parameters are determined assuming that the observed radial velocities refer to the secondary BH, which is in a circular orbit. See Equation ( 3 ) in Section  4 and associated discussion.

	Minimum		Minimum	Minimum	Minimum
Object	Period222Expressed in the rest frame.	$V_2\sin i$333The time of the next inferior conjunction of the secondary corresponding to the minimum period.	$t_0$333The time of the next inferior conjunction of the secondary corresponding to the minimum period.	Total Mass444The broad H$\beta$ profile varied substantially between observations, thus a meaningful shift could not be determined. These variations may include large changes in the profile shape or changes in the width that cause the two wings of the line to move in opposite directions.	Orbital Separation444Computed using the the values of the rest-frame minimum period and $V_2\sin i$ reported in this table.
(SDSS J)	(years)	(km s-1)	(year)	(M⊙)	(${10}^{-2}$ pc)
001224	180	2000	2033	$2.5\times {10}^7$	12.6
015530	150	1400	2050	$5.3\times {10}^6$	11.1
031715	170	1500	2073	$9.5\times {10}^6$	11.7
074007	190	1700	2084	$1.1\times {10}^6$	11.0
082150	150	1300	2061	$4.7\times {10}^6$	11.1
091928	120	4300	2015	$5.9\times {10}^7$	11.3
092712	130	2600	2022	$2.2\times {10}^7$	11.3
093844	110	1910	2044	$3.2\times {10}^6$	11.6
094603	150	2000	2050	$1.5\times {10}^7$	11.5
094620	150	1800	2026	$1.1\times {10}^7$	11.5
095036	120	1700	2012	$3.7\times {10}^6$	10.6
095539	190	1640	2079	$9.2\times {10}^5$	10.9
105041	140	3600	2052	$2.5\times {10}^8$	18.1
110556	170	1200	2028	$4.2\times {10}^6$	11.3
111537	100	2200	2044	$3.7\times {10}^7$	13.3
113330	190	1200	2052	$6.0\times {10}^6$	11.7
113904	130	1500	2156	$1.8\times {10}^7$	13.2
115158	120	1390	2025	$4.7\times {10}^4$	10.1
120924	150	1800	2024	$1.1\times {10}^7$	11.4
125142	110	1700	2043	$1.9\times {10}^7$	12.9
131945	140	1000	2019	$1.5\times {10}^6$	10.6
140251	130	1980	2109	$4.8\times {10}^6$	12.1
151132	190	1100	2100	$5.4\times {10}^6$	11.7
153636	240	3500	2092	$3.8\times {10}^8$	13.4
153644	150	1300	2067	$4.9\times {10}^6$	11.2
154637	130	2600	2019	$1.8\times {10}^7$	11.1
155654	100	3900	2037	$1.9\times {10}^8$	15.7
161911	120	1510	2016	$1.2\times {10}^5$	10.2
162914	140	1930	2051	$1.4\times {10}^6$	10.7

Table 3. Table 3: Parameters of prior distributions a

Parameter	Distribution	Lower Bound	Upper Bound
$P$ (years)……………	Jeffreys	$2T$	$2.652\times {10}^6/u_{2,3}^{\mathrm{obs}}\left(t_1\right)$
$V_2\sin i$ (km s-1)…	Modified Jeffreys	0	$c$
$t_0$ (years)……………	Jeffreys/Uniform222We do not have strong physical constraints on the prior distribution of the time of the next conjunction, $t_0$. Therefore, we try two fairly different distributions.	$t_1$	$t_1+T+P$
$u_2\left(t_1\right)$(km s-1)333This parameter is effectively not free because the prior distribution is constrained by measurements of the SDSS data point.….	Gaussian	All reals

Equations

$$u_2(t)=(V_2\sin i)\sin \left[\frac{2\pi }{P}(t-t_0)\right],$$

$$\Delta u_2(t)$$

$$M>3.8\times 10^5(1+q{)}^3\left(\frac{P_{\min }}{10\text{yr}}\right){\left(\frac{V_2\text{sin}i}{10^3{\text{km s}}^{-1}}\right)}^3{\mathrm{M}}_{\odot },$$

Full text

††pagerange: A Large Systematic Search for Close Supermassive Binary and Rapidly Recoiling Black Holes - III. Radial Velocity Variations–A Large Systematic Search for Close Supermassive Binary and Rapidly Recoiling Black Holes - III. Radial Velocity Variations††pubyear: 2016

1 Introduction

Supermassive black hole binaries (SBHBs) are thought to be a common, if not inevitable, outcome of hierarchical galaxy evolution scenarios (e.g., Menou et al., 2001; Volonteri et al., 2003; Hopkins et al., 2006) in which most massive galaxies host a central supermassive black hole (BH; e.g., Kormendy & Richstone, 1995). We now know of many examples of dual active galactic nuclei (AGNs) in interacting galaxies, with separations of a few kpc (Comerford et al., 2009; Wang et al., 2009; Smith et al., 2010; Liu et al., 2010b, a; Shen et al., 2011; Fu et al., 2011; Fu et al., 2012; Barrows et al., 2012; Comerford et al., 2012; Comerford et al., 2013, 2015; Müller-Sánchez et al., 2015; McGurk et al., 2015; Koss et al., 2016). These are thought to represent the initial stage of the evolution of SBHBs, as described early on by Begelman et al. (1980). However, there are no known unambiguous cases of close, bound SBHBs, with separations of 1 pc or less, which would represent the late stages of the evolution of such systems. A great deal of effort has recently been devoted to finding close SBHBs because they afford tests of supermassive black hole and galaxy co-evolution scenarios and because they are the progenitors of low-frequency gravitational waves detectable by future space-based missions (Amaro-Seoane et al., 2012; Arzoumanian et al., 2014; Sesana, 2015).

Begelman et al. (1980) describe in detail the evolution of such a SBHB via processes that extract energy and angular moment from the system. Initially, the supermassive black holes sink to the centre of the gravitational potential well in the merger remnant and the binary orbit decays due to dynamical friction. The duration of this stage is relatively short, $\sim 10^{6}$ years. Once the orbital velocity becomes comparable to the stellar velocity dispersion the SBHB is considered “hard.” Scattering events with stars in the nucleus are now the dominant angular momentum loss mechanism for the SBHB. If scattered stars are not replaced, the SBHB evolves slowly to separations of $\sim 1$ pc. At separations on the order of $10^{-2}$ or $10^{-3}$ pc, the emission of gravitational waves takes over and the SBHB evolves to coalescence on a time scale of $\sim$a few$\times 10^{8}\;$years. It is not clear whether the process of scattering stars alone can bring the SBHB to the point where gravitational wave emission becomes efficient within a Hubble time. The 3-body scattering process requires a sufficient supply of stars that occupy a certain region of phase space and as these are depleted, the binary separation may stall at $\sim 1\;$pc. Solutions to this “last parsec problem” include replenishing the population of stars available for scattering (e.g., Yu et al., 2005), non-spherical potentials (e.g., Yu, 2002; Merritt & Poon, 2004; Berczik et al., 2006; Khan et al., 2013; Vasiliev & Merritt, 2013; Vasiliev et al., 2014, 2015; Holley-Bockelmann & Khan, 2015), or interactions between the SBHB and a reservoir of gas (e.g., Armitage & Natarajan, 2002; Escala et al., 2004; Dotti et al., 2007, 2009; Cuadra et al., 2009; Lodato et al., 2009).

In light of this evolutionary scenario, it is particularly interesting to search for SBHBs with sub-parsec orbital separations. The chances of finding systems in this regime may be high because, for likely evolutionary pathways, the residence time is relatively long (\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{>}}}10^{8} years). However, the details of the evolution of systems through these separations are not known, so the detection of such SBHBs would also contribute much needed constraints on relevant physical models. Here, we implicitly consider a model that involves interaction of the SBHB with a gaseous circumbinary disc (e.g. Cuadra et al., 2009). This disc acts as a source of fuel for the BHs and as a means of removing angular momentum from the system (see, for example, the models of Artymowicz & Lubow, 1996; Hayasaki et al., 2007; Cuadra et al., 2009). Thus, accretion on to one or both of the BHs creates an observational signature that one can exploit to search for such systems. Depending on the mass ratio of the two BHs, either one or both of them may accrete at a substantial rate. If the mass ratio differs significantly from unity, the secondary (less massive) BH has easier access to the gas in the circumbinary disc, while the primary resides close to the centre of mass of the system.

Two methods have been widely used to search for SBHBs, under the assumption that at least one of the BHs is actively accreting. One approach is to search for variability in the light curves of AGN, assuming that the accretion is modulated on time-scales of the order of the orbital period of the SBHB. Thus far, these searches have identified candidates with periods of order several years or less (e.g., Graham et al., 2015b, a; Liu et al., 2015; Charisi et al., 2016). If these are SBHBs, they likely correspond to the last phase of binary evolution where separations are small enough that the decay of the orbits is driven by gravitational wave radiation. A significant uncertainty in this approach arises from the fact that regular quasar light curves exhibit red noise (e.g., MacLeod et al., 2010, and references therein), blurring the distinction between a periodic signal and normal quasar variability when only a few oscillations of the light curve are observed (see discussion by Vaughan et al., 2016). The second approach is based on an analogy to spectroscopic binary stars and seeks to identify periodic variability in the offsets of quasar broad emission lines (e.g., Gaskell, 1996) indicative of bulk motion through space of one or both BHs and associated broad-line regions (BLRs). The double-line spectroscopic binary case has been tested and rejected for quasars with double-peaked broad emission lines (Eracleous et al., 1997; Liu et al., 2016). Quasars whose broad emission lines are displaced by a substantial amount represent the single broad line spectroscopic binary case, which remains a viable way to look for the signature of orbital motion using quasar spectra. To be sure, this method does have drawbacks and is subject to several uncertainties, which we discuss in later sections of this paper. A “hybrid” approach has also been used to look for SBHBs in nearby Seyfert galaxies with long and well-sampled photometric and spectroscopic time series (e.g., NGC 4151 and NGC 5548; see Bon et al., 2012, 2016; Li et al., 2016). In this approach one looks for periodic behavior in the light curves and also tries to find regular variability of the broad Balmer line profiles on the same period. This approach is subject to the same uncertainties resulting from red noise in the light curves. Moreover, interpretation of the variability of the line profiles involves fitting them with linear combinations of several Gaussian components and assigning physical meaning to these components, which is highly subjective at best.

Starting with Boroson & Lauer (2010), a significant effort has been invested in observational searches for close, SBHBs via the single broad line spectroscopic binary test. The effort has been directed towards both the population of quasars with broad lines located at their systemic velocities (in the context of the SBHB hypothesis such cases would correspond to binaries in conjunction; see Shen et al., 2013; Ju et al., 2013; Wang et al., 2016) and those where the broad emission lines are offset from the rest frame by thousands of kilometers per second (Tsalmantza et al., 2011; Eracleous et al., 2012; Decarli et al., 2013; Liu et al., 2014). So far, these searches have sought to test the SBHB hypothesis for candidates using multiple (generally two) spectra taken years apart to measure or constrain the projected acceleration. The main caveat to this approach is that normal quasar spectral variability, that can mimic the sought after acceleration, has not been characterized on these time-scales. Thus, in order to gain confidence in these results a long-term spectroscopic monitoring campaign is necessary. Long-term spectroscopic monitoring would also help test the recoiling BH hypothesis for the offset of the broad emission lines of the quasars. Recoils of up to several thousand km s*-1* can result from the merger of two BHs and can displace the product of the merger from the centre of the host galaxy (e.g., Campanelli et al., 2007a, b; Baker et al., 2008). These BHs can accrete either from a reservoir of gas they carry with them or from gas that they encounter after they are ejected (e.g. Loeb, 2010; Guedes et al., 2011; Blecha et al., 2011; Blecha et al., 2016). On time scales of order decades we do not expect substantial and monotonic velocity variations in the emission lines from such systems. Such is the case with QSO E1821+643, where the lack of velocity variations over a period of 24 years has been used to argue in favor of the recoiling BH over the SBHB scenario (Shapovalova et al., 2016), although it is possible to explain the observations without invoking either.

As part of a systematic, long-term effort we have been spectroscopically monitoring the 88 SBHB candidates of Eracleous et al. (2012, hereafter Paper I) since 2009 in order to test the SBHB and recoiling BH hypotheses. In Runnoe et al. (2015, hereafter Paper II) we presented the followup spectra, measurements of spectroscopic properties from all the spectra, and an analysis of the variability of the integrated fluxes of the broad H$\beta$ lines. We found that the continuum and broad H$\beta$ lines of SBHB candidates have the same variability properties as those of typical quasars of the same luminosity and redshift. We concluded that, if the line variability is driven by continuum variability, then the BLRs of the SBHB candidates have sizes comparable to those of typical quasars. The goal of this paper, the third in the series, is to present the radial velocity curves derived from the spectra we presented in Paper II. As we noted in section 2 of Paper I, one may expect radial velocity changes of up to a few hundred km s*-1* over the course of our monitoring period, assuming that the secondary BH is active (see the physical picture described in Paper I and re-iterated above). In Section 2 we describe how the sample was identified and its properties. The shift measurements and methodology are presented in Section 3.1 with the most significant caveats to this approach discussed in Section 3.2. We analyse the radial velocity curves in the context of the SBHB hypothesis in Section 4 and discuss them in the context of previous work in Section 5. Our results are summarized in Section 6. Throughout this work, we adopt a cosmology with $H_{0}=73$ km s*-1* Mpc*-1*, $\Omega_{\Lambda}=0.73$, and $\Omega_{m}=027$.

2 Sample and Data

The SBHB candidates were selected from SDSS DR7 quasars with $z<0.7$ to have broad H$\beta$ emission lines that were offset by \lower 3.00003pt\hbox{,\buildrel>\over{\sim},}1,000\;{\rm km\;s}^{-1} relative to the rest frame set by the [O$\,\textsc{iii}]$ $\lambda$5007 narrow line. We note that some redshifts have been refined since that data release. Most notably, the redshift of J092712 was listed as 0.697 in DR7 putting this object in the redshift range of our sample, while its new redshift is 0.713. Potential candidates were initially identified via the spectral principal component analysis described by Boroson & Lauer (2010) and were added to the sample after a single offset peak in the broad line was verified by visual inspection. The final sample includes 88 candidates, with redshifts of $0.077<z<0.713$ and absolute V-band magnitudes of $-21.25<V<-26.24$ after Galactic extinction corrections but not K corrections. Additional details on the sample selection and properties are available from Paper I. In comparison to other similar surveys, the objects in our sample and the sample of Decarli et al. (2013) are at the high end of the distribution of broad H$\beta$ velocity offsets. The objects in the samples of Shen et al. (2013) and Ju et al. (2013) have velocity offsets close to zero by selection. The objects in the sample of Liu et al. (2014) have intermediate offsets (typically less than 1,000  km s*-1*)111In the context of the single broad line spectroscopic binary analogy, described in §1, and assuming circular orbits and a relatively narrow range of masses, periods, and orbital inclinations the objects with the highest velocity offsets may correspond to orbital phases close to quadrature, where the residence times are large while the accelerations are small. Here, we are referring to the velocity and acceleration projected along the line of sight (for circular orbits the true velocity and acceleration have a constant magnitude and vary in direction with orbital phase). In contrast, the objects with velocity offsets close to zero may be close to conjunction, where the residence times are small and the accelerations are large. Objects with intermediate velocity offsets may represent intermediate orbital phases..

We have been monitoring the candidates spectroscopically with a variety of ground-based facilities since late 2009. In addition to the original 88 SDSS spectra, we have obtained 320 follow-up spectra spanning time intervals from a few weeks to 12 years in the observer’s frame. Each object has been observed at least twice and many have three or more observations. All of the data were presented in Papers I and II; we refer the reader to those papers for details of our observations and data reduction.

The wavelength calibration is critical for measuring radial velocity shifts so we summarize it here. For each observing run, we derived a low-order (order$\;<5$) polynomial wavelength solution based on between 20 and 60 arc emission lines. We required that the standard deviation of the residuals from this fit be less than 0.1 pixel, which set the uncertainty in the relative velocity scale for each spectrum at 6400 Å, the median wavelength of H$\beta$ in our sample, to be between 5 and 11 km s*-1*, depending on the telescope and instrument. The absolute wavelength scale was set with the help of the [O$\,\textsc{i}]$ $\lambda$5577 night sky line, which was recorded simultaneously with the quasar spectrum. After applying heliocentric velocity corrections and converting the wavelengths scale to vacuum values, we aligned the spectra of the same object by comparing the wavelengths of [O$\,\textsc{iii}]$ $\lambda\lambda$4959,5007 doublet lines and using the cross-correlation method described in Paper I. The uncertainty in this alignment translates directly into how well we can determine the relative velocity of the broad H$\beta$ line between epochs. This uncertainty is smaller than 55 km s*-1* for all objects where we detect a significant velocity shift and smaller than 35 km s*-1* for 80 per cent of these objects.

The first-epoch velocity offset in the broad H$\beta$ emission line, which anchors the radial velocity curve, was measured from the SDSS spectra of our objects and reported in Paper I. In summary, the peak was located by fitting the region around it with a Gaussian. To assess the uncertainties resulting from the finite signal-to-noise ratio (S/N) and possible asymmetries of the profile near the peak, we repeated the fitting exercise several times, changing the range of pixels that we fitted each time. In the end we took the average peak wavelength from our trials to be the best estimate of the true value and we adopted the full range of peak wavelengths from those trials as a measure of the uncertainty. As such, the error bars on the peak velocity measured from the SDSS spectrum effectively indicate the 99 per cent confidence limits. Velocity changes are then measured relative to the first-epoch spectrum, as we describe in the next section.

3 Broad H$\beta$ radial velocity variations

3.1 Measurement of velocity variations

The first step in measuring the velocity shifts of the broad H$\beta$ emission line from the initial offset observed in the SDSS spectrum is to apply the $\chi^{2}$ cross-correlation method of Paper I. In this approach, we shift the follow-up spectrum past the original SDSS spectrum in small steps. At each step, we rebin the observed follow-up spectrum to the wavelength scale of the SDSS spectrum and renormalize it to match the flux scale of the SDSS spectrum, as described by equation (6) of Paper I. This normalization effectively highlights changes in the shape of the broad-line profile, but it hides changes in the broad-line flux and continuum. These can still be identified from the narrow lines in Figures 1 and 2; an apparent mismatch in the strength of the narrow lines indicates that the broad-line flux has changed between observations. Finally, we measure the $\chi^{2}$ statistic in a window that includes emission from the broad H$\beta$ line and avoids the narrow lines, effectively excluding their contribution from the cross-correlation signal. In this way, we build a $\chi^{2}$ curve as a function of shift and can identify the best shift and uncertainties at the 68 per cent and 99 per cent confidence levels (‘$1\sigma$’ and ‘$2.6\sigma$’ respectively). Since we make this measurement on the observed spectra, the resulting shifts do not depend on any parametric model for the emission-line profiles. We also repeat the measurement with the SDSS spectrum rebinned to the wavelength scale of the follow-up spectrum, so that the shift is measured twice per follow-up spectrum. An example of a diagnostic plot for one velocity shift measurement is shown in Figure 1. The measured shifts are reported in Table 1 where we give the average between the two measurements calculated for each follow-up spectrum and the ‘$1\sigma$’ and ‘$2.6\sigma$’ uncertainties.

Following Paper I, we carried out simulations to check that the uncertainties resulting from the cross-correlation method were not underestimated. In these simulations we cloned the higher-S/N spectrum in each pair and injected in it synthetic noise to match that of the lower-S/N spectrum. We then cross-correlated the two spectra. After 1,000 trials we examined the distribution of shifts about zero to verify that the range encompassing 99 per cent of the trials agreed with the 99 per cent confidence limits derived from the $\chi^{2}$ cross-correlation method. Indeed, the agreement was excellent.

To guard against profile shape variability at least two authors visually inspect all the diagnostic plots for each shift measurement, and a third performs this task if the first two provide differing visual classifications. Visual inspection is used in conjunction with the minimum value of $\chi^{2}$ from the cross-correlation procedure to determine whether profile shape variability has influenced the shift measurement. Moreover, we require that the two shift measurements that result from rebinning the SDSS or follow-up spectrum are consistent with each other in magnitude and have opposite signs. Thus, measurements that appear to be affected by profile shape variability are deemed unreliable, while those that appear not to be affected by profile variations are hereafter referred to as reliable.

Following Paper I, we place each shift measurement into one of five categories, noted in Table 1: (1) we reliably detect a statistically significant shift, (2) we reliably determine that the shift is consistent with zero within measurement uncertainties, (3) the S/N in the spectrum is too poor to obtain a meaningful shift measurement (the cross-correlation algorithm is fooled by noise spikes), (4) we obtain a shift measurement, albeit with larger uncertainties because of slight profile shape variability, (5) the profile varied substantially between observations and no shift measurement can be made. An example of profiles yielding a reliable shift measurement (category 1) is shown in Figure 1. Examples of profiles from categories 2–5 are shown in Figure 2. It is specifically the measurements from categories 1 and 2 that were deemed “reliable.” We find that 29/88 objects have at least three reliable measurements. Ignoring the low S/N measurements, we can build reliable radial velocity curves for these objects that are useful for further analysis. We show these radial velocity curves in Figure 3. It is also noteworthy that 39/88 objects exhibit substantial profile shape variability in at least one spectrum (i.e., category 5 above). While we cannot rule out the SBHB hypothesis for objects whose broad H$\beta$ profiles have varied substantially, we are unable to measure the radial velocity changes in these cases.

3.2 Caveats and Uncertainties

The most significant caveat to this method of searching for SBHBs is that variability of the broad H$\beta$ profiles may mimic the radial velocity shifts that we hope to measure. Indeed, based on our visual inspection, 44 per cent (39/88) of the sample showed substantial profile shape variability over the course of our campaign. Profile shape variability at smaller amplitudes would not be apparent in a visual inspection but could be mistaken for a velocity shift. There are multiple mechanisms that can produce profile shape variability on time-scales of the order of the span of our observations, which we describe below.

On time-scales comparable to the variability time-scale of the continuum source, the broad-line centroids or peaks can display velocity shifts as a result of changes in the illumination of the BLR that are similar to those we expect from orbital motion. Barth et al. (2015) highlight this issue in the context of searches for SBHBs and demonstrate the danger based on high-quality reverberation mapping data for several well-studied AGN. If one side of the broad-line profile responds to variations in the continuum source before the other, the centroid of the line profile can shift by up to a few hundred km s*-1*, even when measured using cross-correlation techniques. Such an asymmetric response can result from an non-axisymmetric distribution of matter or the presence of radial motions in the BLR. Examples of BLRs with such properties have been found through reverberation mapping campaigns with sufficiently high-S/N spectra (e.g., Bentz et al., 2010; Grier et al., 2013; Skielboe et al., 2015). Searches for SBHBs based on radial velocity changes of line profiles that are close to their systemic velocity are particularly susceptible to this effect. We cannot rule out the possibility that this effect is also at work in our targets since we do not monitor the continuum variations of our targets. Even in a SBHB, this type of variability might be expected and will produce jitter on top of a systematic velocity change. Thus, long-term monitoring is crucial in order to measure long-term radial velocity variations.

On time-scales comparable to the dynamical time-scale in the BLR, of order years or longer, the broad-line profile may vary independently of the continuum source as a result of structural changes in the BLR (e.g., redistribution of the gas in position and/or velocity space). The variations of the line profiles resemble a “see saw” pattern and can be quantified by following the magnitude and sign of the skewness. This behavior has been observed in both quasars (e.g., Kassebaum et al., 1997; Romano et al., 1998; de Diego et al., 1998) and Seyfert galaxies (e.g., Peterson et al., 1999; Kollatschny et al., 2000; Wang & Li, 2011; Doroshenko et al., 2012; Schimoia et al., 2012). The same behavior is also fairly common in the profiles of broad, double-peaked Balmer lines found in some radio-loud AGNs (e.g., Eracleous et al., 1997; Sergeev et al., 2007; Gezari et al., 2007; Lewis et al., 2010; Popović et al., 2014). Thus, we performed simulations to investigate whether profile shape variability can masquerade as the radial velocity shifts we hope to measure. We describe the simulations in detail in Appendix A. In summary, we induced both shape variability and a shift of the broad H$\beta$ profile in the observed spectra of the SBHB candidates, and applied the same cross-correlation method to measure shifts. We found that it is much easier to mistake profile shape variability for a bulk shift of the broad line in broad boxy profiles, while in cuspy profiles it is easier to notice in a visual inspection that the profile has changed. However, we were able to find examples of both boxy and cuspy profiles where the profile shape variability and bulk shift produced visually indistinguishable shift measurements, as we illustrate in Appendix A. We conclude that it is possible to be misled by this type of profile shape variability into thinking that we are measuring radial velocity shifts. Again, long-term spectroscopic monitoring is the best way to distinguish between profile shape variability and real shifts. Moreover, independent tests of the nature of the SBHB candidates can bolster the radial velocity results.

4 Interpretation of radial velocity variations in the context of the binary hypothesis

It is useful to consider whether the apparent radial velocity variations that we have measured in our campaign so far are consistent with those expected from orbital motion. To this end, we present here a methodology for testing whether the observed radial velocity curves are consistent with sinusoids and obtaining constraints on the orbital parameters. In view of the long periods that we expect, the main attainable goal of this exercise is to set restrictive lower limits on the BH masses that make the SBHB hypothesis untenable for some systems. Even objects in which no significant variations are observed can lead to interesting limits on the properties of the hypothesized SBHBs so we consider all objects with reliable radial velocity curves, even if their radial velocity variations are consistent with zero within uncertainties. We follow a variant of the methodology described by Eracleous et al. (1997) (see also Halpern & Filippenko, 1988, 1992; Liu et al., 2016) to identify the minimum period that is consistent with the observed radial velocity curve. We assume that both members of the binary follow circular orbits around the centre of mass and only the secondary BH is active.

Our assumption of circular orbits is motivated by both practical and theoretical considerations. From a practical point of view, the radial velocity curves shown in Figure 3 often include three or four measurements, which are not enough to constrain the free parameters of an elliptical orbit. In some cases the radial velocity curves include five to seven points but in these cases there are closely separated pairs of measurements, which hampers our ability to get meaningful constraints on the parameters of an elliptical orbit. From a theoretical point of view, predictions for the orbital eccentricity of SBHBs vary substantially and depend on the formation channel of the system and the mass of the circumbinary disc. Simulations of the evolution of a SBHB through interactions with stars suggest that the eccentricity is damped very quickly (Khan et al., 2013; Holley-Bockelmann & Khan, 2015). However, the question of eccentricity evolution in these models is not settled because it depends sensitively on the initial conditions, the SBHB mass ratio, and other properties of the system (e.g., Quinlan, 1996; Sesana et al., 2006; Wang et al., 2014; Vasiliev et al., 2015), and can be high (Vecchio et al., 1994; Aarseth, 2003). Eccentricity evolution in models where the SBHB interacts with a gaseous reservoir is also uncertain. In models where the SBHB initially evolves in a large, smooth, gaseous disc (e.g. Dotti et al., 2006), the eccentricity is damped very quickly. However, if the disc includes massive clumps, the interaction of the BHs with these clumps imparts a substantial eccentricity to their orbits (e.g. Fiacconi et al., 2013; del Valle et al., 2015). At the late stages of the evolution of the SBHB, after it has opened a gap in the gaseous disc, the evolution of the eccentricity depends on the mass of the disc exterior to the binary orbit. If the disc is sufficiently massive, the combination of torques from the disc and accretion streams will give rise to an appreciable eccentricity, even if the orbit is initially circular (see Roedig et al., 2011, and references therein). Otherwise, the eccentricity will not be affected by interactions of the SBHB with the circumbinary disc (e.g. Hayasaki et al., 2007; Farris et al., 2014). Thus, there is a wide variety of physical scenarios in which SBHBs have circular orbits.

Under the assumption of circular orbits, we write the model for the radial velocity curve of the secondary as

[TABLE]

where $u_{2}$ is the observed velocity of the secondary, $V_{2}$ is the (true, tangential) orbital velocity of the secondary, $i$ is the inclination angle between the normal to the orbital plane and the line of sight, $P$ is the orbital period, and $t_{0}$ is the time of the next inferior conjunction of the secondary (i.e., the time at which the secondary starts to recede from the observer and its projected position is between the primary and the observer; if the orbital plane were viewed edge on, this would be the time of eclipse of the primary by the secondary). Moreover, we denote the time of the $n^{th}$ observation in our time series as $t_{n}$ and the observed radial velocity at that time as $u_{2}^{\rm obs}(t_{n})$ (the superscript “obs” identifies measured radial velocities). In this convention, $u_{2}^{\rm obs}(t_{1})$ is the radial velocity measured from the first spectrum in the series, taken as part of SDSS.

We have to treat the uncertainty in the first radial velocity measurement differently from the uncertainties in subsequent measurements because the first measurement is an absolute one while subsequent measurements are differential (relative to the first one). Therefore, we re-cast the model for the radial velocity curve as

[TABLE]

where $u_{2}(t_{1})$ is the initial velocity offset in the time series and $\Delta u_{2}(t)$ represents a velocity change relative to the initial value. Following the convention of the previous paragraph, we use $\Delta u_{2}^{\rm obs}(t_{n})$ to identify differential velocity measurements. By writing the model in terms of $\Delta u_{2}(t)$, we can use the uncertainties in relative velocity changes determined as described in §3.1.

In order to test this model, we ran the Markov Chain Monte Carlo (MCMC) sampler emcee (Foreman-Mackey et al., 2013), which we describe in detail in Appendix B. In summary, the simulations determine the likelihood of our data given the above model, and produce posterior distributions for all model parameters and an accurate accounting of their covariance. We then generate the cumulative posterior distribution of $P$, $\psi(>\!P)$, from which we can determine the 68th, 90th, and 99th percentile minimum periods and calculate the corresponding orbital models. In Figure 4 we show two examples of the ultimate output of the MCMC simulations, one for an object with a flat radial velocity curve and one where jitter has likely been interpreted as orbital motion. The output consists of the $\psi(>\!P)$ distributions and the radial velocity curve with the orbital models corresponding to the 68th, 90th, and 99th percentile minimum periods superposed. These models are illustrated for all objects in Figure 3.

We take the period corresponding to the 99th percentile of $\psi(>\!P)$ to be the minimum period that is consistent with the observations, $P_{\rm min}$. We combine this period with the corresponding velocity normalization of the sinusoid to place a lower limit on the total mass and orbital separation of the SBHB. In the framework described in §2 of Paper I,

[TABLE]

where we have defined $M\equiv M_{1}+M_{2}$ as the total mass of the SBHB, $q\equiv M_{2}/M_{1}\leq 1$ as the mass ratio, and we note that $(1+q)^{3}>1$. If, instead, we were observing the radial velocity variations of the primary, we would replace $V_{2}$ by $V_{1}$ and $(1+q)^{3}\geq 1$ with $(1+q)^{3}/q^{3}\geq 8$ in Equation (3). Thus, any limits on the total mass would be multiplied by 8. The limits on the periods and total masses, derived using Equation (3), for objects with reliable radial velocity curves are given in Table 1. The physical models corresponding to the 68th, 90th, and 99th percentile periods are shown in Figure 3.

With the data available so far, none of the limits on the total mass of the SBHB are sufficiently stringent to rule out the binary hypothesis. According to Equation (3), for velocity offsets of order a few$\times 10^{3}$ km s*-1* we would need a very high minimum period, $\sim 10^{4}$ yrs, to obtain a limiting mass $\sim 10^{10}$ M*⊙. With only 3–4 observations, it is rare that the radial velocity curve is actually flat over the duration of our campaign. In other words, velocity jitter and finite measurement uncertainties can make a flat radial velocity curve appear curved if only a few measurements are available. Specifically, in 10/88 cases, we measured statistically significant radial velocity variations, which introduce curvature into the radial velocity curves and lead to less restrictive constraints on the periods. In 19/88 cases we found that all of our reliable radial velocity measurements were consistent with the original SDSS measurement at the 99 per cent confidence level, however the sampling pattern and measurement errors allow for substantial curvature in the orbital models. For example J153636, which has a stable broad line profile observed over six years but substantial jitter, yields the largest mass limit at $M>3.8\times 10^{8}$ M⊙*. In the case of J105041, the 2012 measurement with its large uncertainty does not constrain the model as much as the 2010–2011 measurements, with the result that shorter periods are possible. From these results, we conclude that five or more radial velocity measurements with small measurement uncertainties are needed to get better period constraints than what we have gotten so far.

5 Discussion

Our main observational result is the measurement of radial velocities of the SBHB candidates in our sample at several epochs. We are able to produce reliable radial velocity curves for 29/88 candidates. As we noted in §3.2, variability of the line profiles is the most serious limitation in our ability to put together long time series of radial velocity measurements. In approximately half of the objects in our sample we have observed substantial changes in the profiles of the broad H$\beta$ lines that make the measurements of radial velocity variations impossible. In objects where the broad H$\beta$ profiles appear stable, there is still a danger that profile shape variability causes radial velocity jitter beyond the measurement errors and complicates the interpretation of the radial velocity curves. The potential effects of profile shape variability were pointed out and discussed by Decarli et al. (2013), Shen et al. (2013), and Liu et al. (2014) who also noted that broad lines with larger velocity offsets are more likely to display profile shape variability. Because the dynamical time of the broad-line region of any one of the two BHs is shorter than the orbital period of the SBHB, we cannot expect that the profiles of the broad lines will remain stable over the course of an orbital cycle. This means that we will typically be able to measure only segments of the radial velocity curve during which the line profiles are stable. In order to put together long-term time series, we will need to devise methods of combining such short segments together and overcoming the problem caused by the changing velocity zero point that is brought about by profile variations.

The above cautions notwithstanding, we would like to draw attention to three objects in our sample, J093844, J095036, and J161911, which have shown systematic and monotonic velocity changes by several hundred kms s*-1*. These represent the most promising cases in our sample and deserve further attention. They are high-priority targets for continued spectroscopic monitoring as well as additional observational tests, such as high-resolution imaging of the host galaxies in search of morphological features associated with a merger.

In addition to making new measurements, we have also developed the methodology for setting a lower limit on the period, hence the mass, of a hypothesized SBHB using a radial velocity curve of finite duration. The method uses an MCMC algorithm to explore the space of parameters of circular orbits and allows one to determine the shortest orbital period consistent with the observed radial velocities at the 99 per cent confidence level. Applying this method to the available data, we find that we are not yet able to put interesting constraints on the SBHB candidates. The limits we give in Table 1 reflect the properties of our radial velocity curves (duration, uncertainties, jitter, etc.) rather than physically interesting constraints on the properties of the hypothesized binaries. With more measurements over a longer baseline, we can obtain more stringent constraints and also attempt to constrain models of elliptical orbits, which is not possible with the data available so far. Another noteworthy point is that the value of $V_{2}\sin i$, derived from the model fits, influences the mass constraints considerably because it enters in the third power in Equation (3). For objects with small observed radial velocities, $V_{2}\sin i$ is likely to be small, therefore a combination of a longer temporal baseline and smaller error bars is needed in order to obtain interesting constraints on the mass of the hypothesized SBHB.

If the orbital periods are of order centuries, we cannot hope to confirm in our lifetimes that any specific objects in our sample are SBHBs based on radial velocity variations. However, we will be able to place more stringent constraints on the properties of the hypothesized binary orbits with continued monitoring. In the spirit of Decarli et al. (2013) we ask how much longer we need to monitor the SBHB candidates before we obtain interesting constraints on their orbital parameters. With simulations, we explore how long it will take to obtain a restrictive mass constraint in the event that the observed radial velocity curve is constant, within the uncertainties and allowing for jitter. We use the observed radial velocity curve of J153636 as the starting point for this experiment, because it has a relatively large velocity amplitude and has remained flat over our monitoring period. We simulate 50 additional observations using the measurement uncertainties and jitter properties determined from the available observations. The simulated observations are made every two years with several gaps to simulate the effects of bad weather. The simulated data for a century of monitoring are shown in Figure 5. After analyzing these data with the procedure of §4, we find that a monitoring period of of $\sim 100$ years yields an interesting mass limit of $M>9\times 10^{9}$ M*⊙(see Figure 5), provided that the time series is not sparsely sampled, i.e., as long as we have an observation every couple of years. For reference, we note that the largest BH masses measured so far (in brightest cluster galaxies, e.g., McConnell et al., 2012; Thomas et al., 2016, and references therein) are approximately 2\times 10^{10}\;{\rm M}_{\hbox{\odot}} while arguments have been made that BH with masses sevaral times higher than this value could exist (see Inayoshi & Haiman, 2016, and references therein). The velocity jitter, which is evident in the observed radial velocity curve for J153636 and which we reproduce in our simulations, significantly hinders our ability to get an interesting mass limit by allowing orbital models with short orbital periods to describe the data adequately. Carrying out an analogous simulation with measurement errors only and no jitter yields a more restrictive mass constraint of $M>1\times 10^{10}$ M⊙*. It may be possible to obtain such interesting constraints within a few decades by introducing a model for the jitter in the MCMC simulation. This is beyond the scope of this work and, notably, would require several more observations of radial velocity changes before a reliable model for the jitter can be derived. Thus, patience and continued monitoring are required for the radial velocity test to achieve its maximum discriminating power.

With additional monitoring we can possibly carry out additional tests of the SBHB hypothesis by looking for radial velocity variations that are inconsistent with orbital motion. The double-peaked emitter 3C 390.3 is a good case in point: the double-line spectroscopic binary scenario was tested and rejected for this object after three decades of monitoring (Liu et al., 2016, and references therein). The velocity offsets of the peaks of the broad Balmer lines abruptly departed from their long-term monotonic trend (positive acceleration) and remained steady for more than a decade before the acceleration changed sign without the velocity passing through zero (e.g., see Fig. 5 of Eracleous et al. 1997 and Fig. 10 of Shapovalova et al. 2001). This behavior is inconsistent with orbital motion. Additional monitoring can also constrain the properties of the BLR, if we assume that our candidates are true SBHBs. As we noted in Paper II and illustrated in Figure 14 therein, constraints on the orbital period of an SBHB translate into constraints on the extent of the BLR around the secondary (active) BH. Thus, if further monitoring increases the lower limit on the orbital periods to values of 500 years or greater, the scenario of a tidally truncated BLR will be favoured over the scenario of a BLR that fills the Roche lobe of the secondary BH.

The shifted broad lines that were used to select the SBHB candidates may also be a signpost of recoiling BHs. As discussed in detail in Paper I, this scenario seems unlikely to produce the observed properties of the objects in our sample. Theoretical models of recoiling BHs suggest that velocity offsets in excess of 400 km s*-1* are rare (e.g., Dotti et al., 2010). Additionally, the illumination of the narrow-line region and the obscuring torus from outside the host galaxy should leave its imprint on the relative strengths and widths of the line profiles and the shape of the spectral energy distribution in the infrared (Bonning et al., 2007; Shields et al., 2009a; Shields et al., 2009b; Lusso et al., 2014), which is not observed in these objects. In addition to these tests, which were carried out after the initial selection of the sample, recoiling BHs are expected to exhibit stable radial velocity curves. The 19 objects with reliable, non-variable radial velocity curves are therefore the best candidates for this interpretation, notwithstanding the disagreement of the high velocity offsets, narrow-line ratios, and IR spectral energy distributions with expectations. Of these, J092712 is particularly noteworthy because recent observations by Decarli et al. (2014) find molecular gas in the host galaxy at the same redshift as the broad emission lines from the quasar. This result strongly disfavors the recoiling BH interpretation for this object. The 10 objects with variable radial velocity curves are unlikely to represent recoiling BHs unless the radial velocity variations can be associated with variability of a single BLR. Based on our spectroscopic monitoring, of these 10 objects there are three (J093844, J095036, J161911) in which the radial velocity curves show monotonic variations by many hundreds of kms s*-1*, incompatible with a single BLR bound to a recoiling BH. It is difficult to evaluate the remaining 7/10 of the objects with reliable variable radial velocity curves because their radial velocities are erratic and consistent with jitter.

6 Summary and Conclusions

In Paper I we identified a sample of 88 SBHB candidates selected based on their broad H$\beta$ emission lines, which are offset by thousands of kilometers per second relative to the frame of the narrow lines. As part of a long-term effort we have been monitoring the candidates spectroscopically since 2009, with the latest installment of data from our campaign, typically 3–4 spectra per object spanning an interval of 12 years in the observer’s frame, reported in Paper II. The main observational result of this paper, the third of the series, is the measurement of radial velocity variations of the broad H$\beta$ lines from these spectra and the presentation of the radial velocity curves of the SBHB candidates. Furthermore, we establish a methodology for evaluating the binary hypothesis by setting a lower limit on the period, and thus the mass, of the hypothetical SBHBs. Our results are:

We identify 29/88 of SBHB candidates where, based on visual inspection, no profile shape variability has occurred and we can make reliable measurements of the radial velocities. This number includes cases where we detect statistically significant radial velocity changes as well as those where the radial velocities are reliably steady. For these objects, we construct and present radial velocity curves. 2. 2.

We identify 3 objects (J093844, J095036, and J161911) in the sample that show systematic and monotonic velocity changes consistent with the binary hypothesis. These sources are high-priority targets for further spectroscopic monitoring and additional observational tests. 3. 3.

We observe substantial profile shape variability in at least one spectrum in 39/88 of the SBHB candidates. This highlights the main caveat associated with searching for SBHBs with our method: variability of the broad-line profile shape can mimic radial velocity variations. Through simulations of profile shape variability, we have demonstrated that subtle profile shape variability may elude visual inspection and masquerade as a radial velocity change. Broad lines with wide, boxy profiles are more susceptible to this effect than cuspy profiles, where it is easier to discern that the profile shape has changed. As a result, it is critical to characterize the long-term variability properties of the broad emission line profiles of regular quasars, and to continue to observe the candidates on time-scales where orbital motion can be distinguished. 4. 4.

Under the assumption that we can interpret the radial velocity variations of the broad H$\beta$ emission line as orbital motion, we develop a methodology for finding the minimum period that is consistent with the observed radial velocity curves at the 99 per cent level. This allows us to calculate the minimum total mass of the SBHB. A sufficiently high mass limit can lead to a rejection of the SBHB hypothesis. The short duration of our campaign so far leads to mass limits between a few$\times 10^{6}$ and a few$\times 10^{8}\;{\rm M}_{\odot}$, which are not restrictive enough.

Continued spectroscopic monitoring will enable substantial gains in several ways. The addition of points to the radial velocity curves will allow us to test more complex physical models, such as eccentric orbits. By characterizing the jitter and noise properties in the radial velocity curves we may be able to uncover the underlying trends and obtain more restrictive mass limits within a few decades of observing. With extended monitoring we also open the door to observing radial velocity variations that are inconsistent with orbital motion. Given that the periods of the SBHBs are expected to be long, abrupt departures from established monotonic changes, as has been observed in double-peaked emitters, are not expected. We therefore conclude that continued spectroscopic monitoring is required to further evaluate the credentials of the SBHB candidates.

7 Acknowledgements

We thank Anh N. D. Doan for his help in testing our local implementation of the emcee code. JCR would like to thank Becky Nevin and Alex Hagen for helpful discussions during the preparation of this work. The authors also wish to thank the referee, Roberto Decarli, for a careful reading of the manuscript and suggestions that improved the paper.

This work was supported by grant AST-1211756 from the National Science Foundation and an associated REU supplement. ME thanks the members of the Center for Relativistic Astrophysics at Georgia Tech and the Department of Astronomy at the University of Washington, where he was based during some phases of this work, for their warm hospitality. SS thanks the Aspen Center for Physics for hospitality and the ACP is supported by NSF grant PHY-1066293. TB acknowledges support from the National Aeronautics and Space Administration under Grant No. NNX15AK84G issued through the Astrophysics Theory Program and support from the Research Corporation for Science Advancement through a Cottrell Scholar Award. We thank the staff at Kitt Peak National Observatory, Apache Point Observatory, and the Hobby-Eberly Telescope for their expert help in carrying out the observations.

This work is based on observations obtained with the Apache Point Observatory 3.5-meter telescope, which is owned and operated by the Astrophysical Research Consortium.

This work is based on observations at Kitt Peak National Observatory, National Optical Astronomy Observatory (NOAO Prop. ID: 2014A-0098; PI: Runnoe), which is operated by the Association of Universities for Research in Astronomy (AURA) under cooperative agreement with the National Science Foundation.

The Hobby-Eberly Telescope (HET) is a joint project of the University of Texas at Austin, the Pennsylvania State University, Stanford University, Ludwig-‘Maximillians-Universität München, and Georg-August-Universität Göttingen. The HET is named in honor of its principal benefactors, William P. Hobby and Robert E. Eberly.

The Marcario Low-Resolution Spectrograph is named for Mike Marcario of High Lonesome Optics, who fabricated several optics for the instrument but died before its completion; it is a joint project of the Hobby-Eberly Telescope partnership and the Instituto de Astronomía de la Universidad Nacional Autónoma de México.

This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

Appendix A Simulated profile shape variability

Variability of the profile shape is an important impediment in our ability to measure velocity shifts of the broad H$\beta$ emission line. In making the radial velocity measurements, profile shape variability was diagnosed by eye and with the help of the $\chi^{2}$ test so that shifts were only measured in cases where the profile was stable. Here, we investigate whether changes in the profile shape can cause us to measure radial velocity shifts and whether we can identify these cases based on visual inspection. That is, we ask if we can actually distinguish between true velocity shifts in the broad-line profile and subtle changes in the profile shape that may occur over the course of several years. In order to address this question, we simulate both profile shape variability and real shifts in the observed spectra of the SBHB candidates and consider the results of our shift measurements.

We specifically simulate different types of “see saw” variability of the line profiles (i.e., the shapes) because this type of variability is often observed (see discussion in §3.2) and because it has the potential of mimicking the effects of a velocity shift. We take advantage of the decomposition of the spectra into continuum and broad and narrow lines, carried out in Paper II, in order to isolate the profile of the broad H$\beta$ line. Thus, we extract the smooth parametrization of the broad H$\beta$ profile, perturb it in a variety of ways, and then put it back in the spectrum.

The perturbations we use take the form of a suppression of one shoulder/wing of the broad H$\beta$ line and a boosting of the other. In practice, we multiply the parametric model for the broad H$\beta$ profile with three families of odd functions: a tilted line of the form $f(x)=1+ax$, a sinusoid of the form $f(x)=1+a\,\sin(b\,x)$, and a modified Gaussian of the form $f(x)=1+a\,x\,e^{-b\,x^{2}}$. In all cases $x=\lambda-\lambda_{peak}\propto\Delta v$ with $\lambda_{peak}$ the wavelength of the peak of the line profile and $\Delta v$ the velocity separation relative to the peak of the profile. The parameter $a$ controls the amplitude of the perturbations and the parameter $b$ controls how the perturbation evolves with velocity relative to the peak of the profile. These functions are illustrated in Figure 6. In the case of the tilted line, the perturbation increases steadily with increasing $\Delta v$. In the case of the sinusoid, the parameter $b$ has different values on the blue and red sides of the peak so that the variability reaches zero where the line profile falls below 1 per cent of the peak value. In the case of the modified Gaussian, the parameter $b$ also has different values on the blue and red sides of the peak so that the maximal perturbation occurs at 40 per cent of the full-width at half-maximum (see references in §3.2). We find that the modified Gaussian creates the most realistic profile shape variability, particularly for boxy line profiles. The tilted line and sinusoid induce unrealistically large variability in the far wings of broad line profiles.

To explore the effect of the simulated perturbations we apply them to the SDSS spectrum of each SBHB candidate in our sample as described above. For each object and perturbation prescription we generate 50 synthetic spectra with a range in the perturbation amplitude. Two examples of synthetic spectra, generated using the modified Gaussian perturbation, are shown in Figure 7, one for a “cuspy” and one for a “boxy” line profile. In a similar fashion, we also generated spectra where the broad H$\beta$ line was actually shifted in order to carry out a visual comparison with the spectra where we induced profile shape variability.

We measure the shifts of the broad H$\beta$ lines in the synthetic spectra relative to the observed spectrum using the methodology of §3.1. Examples of the results of this exercise are shown in Figures 8 and 9 for the same two objects shown Figure 7. We have chosen these examples to illustrate that shifted and perturbed profiles produce an effectively indistinguishable signal. Visual inspection in such cases does not help us distinguish between shifts and profile variations because the magnitude of the perturbations is small. Such cases are common in the tests that we have carried out; we encounter them in 24 per cent of the objects in our sample. Moreover, we find that genuine shifts are more difficult to distinguish from profile shape variability if the broad line profile is “boxy” e.g., J095539 in Figure 7b and Figure 9. Of course, our ability to distinguish profile shape variability from real shifts using visual inspection and a quantitative comparison of the line profiles via the $\chi^{2}$ test depends critically not only on the shape of the profile but also on the type of perturbation and the S/N of the spectra. Therefore, we cannot make a reliable estimate of the likelihood that any of the shifts that we measure may be false. Nor can we formulate specific criteria for identifying cases that are particularly susceptible to this effect. None the less, the simulations we have carried out do indicate that there is cause for concern. The best way to guard against false positive signals of this type is continued spectroscopic monitoring. By extending the time series of spectra we can see if perturbations in the line profiles grow. Also, the long-term radial velocity curves will show different patterns in cases of shifts and profile variations.

Appendix B Markov Chain Monte Carlo simulations of the radial velocity curves

In order to determine the minimum period consistent with the observed radial velocity curves at the 99% confidence level, we adopt the physical model and measurements described in Section 4, and run an MCMC simulation which we describe here in detail. We use the code emcee222emcee is written in Python and made publicly available via http://dan.iel.fm/emcee/current/ (Foreman-Mackey et al., 2013), which implements an affine invariant MCMC sampler. This code determines the likelihood of the data given the model in Equation 2 and produces posterior probability distributions of all model parameters. As inputs the code takes the data, prior distributions for all model parameters, and an initial distribution of “walkers,” i.e., starting points of trajectories for exploring the model parameter space. In practice, a “burn-in phase” is required before the simulation will yield meaningful results. That is, initial, throw-away iterations of the simulation are needed in order to transform the input parameter and walker distributions to a distribution of walkers in parameter space that reflects the desired posterior distributions. The simulation itself is then started with the walker distribution reached at the end of the burn-in phase and all other input distributions as supplied initially by the user.

This code also allows us to properly treat the uncertainty in the first radial velocity measurement, as follows. We allow $u_{2}(t_{1})$ in Equation (2) to be a free parameter but we constrain it to follow a Gaussian distribution with a mean equal to the measured value of $u_{2}^{\rm obs}(t_{1})$ and a standard deviation equal to the ‘$1\sigma$’ uncertainty on the measurement. For the first measurement in the relative radial velocity curve, we adopt $\Delta u_{2}^{\rm obs}(t_{1})=0\;{\rm km\;s}^{-1}$ and an uncertainty of 1 km s*-1*, much smaller than the other points, effectively forcing the model to follow the $u_{2}(t)$ model parameter. Other inputs to this code (i.e., “priors”) include the distributions of the orbital period, $P$, the time of the next conjunction, $t_{0}$, and the amplitude of the radial velocity variations, $V_{2}\,\sin\,i$. We choose the orbital period distribution to be $\psi(P)\propto 1/P$ (i.e., a Jeffreys prior, flat in $\log P$) and ranging from $2T$ (where $T$ is the duration of the monitoring period for a particular object) to a maximum value, $P_{\rm max}$, chosen as follows. We use equation (1) of Paper I, which relates the orbital period to the measured velocity offset and the mass, phase angle, and inclination angle of the binary orbit. We set the total SBHB mass to $10^{11}\;\rm M_{\odot}$ (see §5) and select the values of other parameters so as to maximize the period ($q\ll 1$, $\sin i=\sin\phi=1$). Thus we get the maximum orbital period corresponding to a mass of $10^{11}\;\rm M_{\odot}$ and compatible with the measured value of $u_{2}^{\rm obs}(t_{1})$, $P_{\rm max}=2.652\times 10^{6}\,[u_{2}^{\rm obs}(t_{1})/10^{3}\;{\rm km\;s}^{-1}]^{-1}\;$yr. With this scheme, the code explores orbital models corresponding to SBHB masses that are at least as large as $10^{11}\;M_{\odot}$. Since the SBHB mass is a derived parameter that depends on the values of $P$ and $V_{2}\,\sin\,i$, the code can actually explore SBHB masses that are larger than $10^{11}\;{\rm M}_{\odot}$ because the value of $V_{2}\,\sin\,i$ obtained from the model fit can be larger than the value of $u_{2}^{\rm obs}(t_{1})$ used to set $P_{\rm max}$. For $t_{0}$ we choose a Jeffreys distribution ranging from $t_{1}$ to $t_{1}+T+P$ (where the upper limit is set during the simulation depending on $P$). We also explore the effect of adopting a flat distribution for $t_{0}$ and find that it yields systematically higher periods by 1 per cent with a scatter of 10 per cent, a negligible difference. For $V_{2}\,\sin\,i$ we choose a modified Jeffreys prior (which is flat below 1 km s*-1*) ranging from 0 to the speed of light. We summarize our selection of prior distributions in Table 3.

We run the emcee code in two stages so as to explore the large parameter space efficiently but also sample the shortest periods fully. First, we perform a high-resolution simulation of the shortest periods. We initialize 500 walkers clustered very close together around an initial guess for each model parameter ($\delta Q/Q\sim 10^{-4}$, where $Q$ is any model parameter). We follow a three-stage burn-in process, where for each stage we run the simulation for $10^{3}$ iterations and then reset the walkers in a cluster around their average value for each parameter. Then a final simulation is run for another $10^{3}$ iterations without resetting the walkers after the last burn-in phase. Second, we perform a lower-resolution simulation that samples the longest periods. Based on the results of the high-resolution simulation, we identify a high-likelihood ridge in the $P-t_{0}$ parameter space. The 500 walkers for the low-resolution simulation are initialized along this ridge according to the priors on $P$ and $t_{0}$, while the walkers for $V_{2}\,\sin\,i$ are distributed in a cluster around the high-resolution average value with dispersion equal to the measured standard deviation. After a burn-in phase of $3\times 10^{3}$ iterations, the low-resolution simulation is run for $10^{3}$ iterations.

Each simulation produces posterior distributions for all model parameters and an accurate accounting of their covariance. We inspect the posterior distributions of $u_{2}(t_{1})$ to confirm that, indeed, the prior constraints were obeyed in all cases. Using the posterior $\psi(P)$ distributions from the high- and low-resolution simulations, we can create a combined, cumulative posterior distribution of $P$, $\psi(>\!P)$. We identify a range of periods where the posterior distributions from the high- and low-resolution simulations overlap and are well populated and we normalize them so that they yield equal probability over this range and then combine them. From the final $\psi(>\!P)$ we can determine the 68th, 90th, and 99th percentile periods and calculate the corresponding orbital models.