VERY LOW-MASS STELLAR AND SUBSTELLAR COMPANIONS TO SOLAR-LIKE STARS FROM MARVELS. VI. A GIANT PLANET AND A BROWN DWARF CANDIDATE IN A CLOSE BINARY SYSTEM HD 87646

Subsequent observations were performed using the ET instrument at KPNO (Ge et al. 2006a). Initial follow-up was performed in 2007 November, which confirmed the variability seen in the KeckET data. Additional data points were obtained at KPNO in 2008 January, February, and May. The integration time was 35–40 minutes in 2007 November and 20 minutes in 2008 January, February, and May.

The data were reduced using software described in Mahadevan et al. (2008) and references therein. See van Eyken et al. (2010) for the theory behind the technique. A total of 40 data points were obtained from 2007 November to 2008 May and are also listed in Table 1. The observations confirmed the linear trend shown in the KeckET data, which will be found later due to another substellar companion.

Follow-up observations of HD 87646 were conducted with the fiber-fed High Resolution Spectrograph (HRS, Tull 1998) of the Hobby Eberley telescope (HET, Ramsey et al. 1998). The observations were executed in queue scheduled mode (Shetrone et al. 2007) and used a 2 arcsec fiber, with the HRS slit set, to yield a spectral resolution of $R\sim {\rm{60,000}}$. A total of 29 data points were obtained between 2007 December and 2008 March. An iodine cell was inserted into the beam path to yield a fiducial velocity reference. The radial velocities were obtained using the procedure and analysis techniques described in Marcy & Butler (1992). The HRS spectra consisted of 46 echelle orders recorded on the blue CCD (407−592 nm) and 24 orders on the red one (602−784 nm). The spectral data used for RV measurements were extracted from the 17 orders (505−592 nm) in which the I₂ cell superimposed strong absorption lines. The radial velocities obtained are also provided in Table 1.

HD 87646 was selected as an RV survey target by the MARVELS preselection criterion (Paegert et al. 2015). The star has been monitored at 23 epochs using the MARVELS instrument mounted on the SDSS 2.5 m Telescope at APO between 2009 May and 2011 December (Ge et al. 2008, Ma et al. 2013). The MARVELS instrument is a fiber-fed dispersed fixed-delay interferometer instrument capable of observing 60 objects simultaneously and covers a wavelength range of 5000–5700 Å with a resolution of R ∼ 12,000. The data processing and error estimation algorithm have been described in detail by Thomas et al. (2016). The final differential RV products are included in the SDSS Data Release 12 (Alam et al. 2015) and are presented in Table 1.

To investigate the nature of the linear RV trend found in previous RV data, we have obtained additional observations of HD 87646 with a fiber-fed echelle spectrograph situated at the 2 m Automatic Spectroscopic Telescope (AST) in the Fairborn Observatory (Eaton & Williamson 2004, 2007). The robotic nature of the AST allowed for high cadence observations, which removed orbital period degeneracies and helped solidify the longer-period companion's orbit. Through 2011 June, the detector was a 2048 × 4096 SITe ST-002A CCD with 15 μm pixels. The AST echelle spectrograph has 21 orders that cover the wavelength range of 4920−7100 Å, and has an average resolution of 0.17 Å. Beginning in 2010 January, several upgrades were made to increase the throughput, sensitivity, and flexibility of the AST (Fekel et al. 2013). In the summer of 2011, the SITe CCD detector and dewar were replaced with a Fairchild 486 CCD having 4K × 4K 15 μm pixels, which required a new readout electronics package, and a new dewar with a Cryotiger refrigeration system. The echelle spectrograms that were obtained with this new detector have 48 orders, covering the wavelength range of 3800−8260 Å. The data reduction and RV measurements are discussed in Eaton & Williamson (2007). A total of 135 data points were obtained from 2009 March through 2013 October and are listed in Table 1. With these additional RV data, we were able to detect the turnaround of the previous identified linear RV trend and start to uncover the second substellar companion's orbit.

A total of nine high-resolution spectra (R = 164,000) of HD 87646 were obtained with the SARG spectrograph on the 3.5 m TNG telescope at La Palma on 2008 March 21, 22, and 28, and April 03 and 11. These data were used to monitor line bisector variations to determine stellar properties (metallicity, log g, ${T}_{\mathrm{eff}}$ and $v\sin i$), and to search for evidence of a second set of lines in the system. The typical signal-to-noise ratio (S/N) for each spectrum is about 150 per resolution element around $5500\,{\mathring{\rm{A}}}$.

We also obtained spectroscopic observations of HD 87646 from the 2.1 m telescope at the Kitt Peak National Observatory using the R = 30,000 Direct Echelle Mode of the EXPERT spectrograph (Ge et al. 2010). A total of seven EXPERT spectra were acquired between 2014 February to 2014 June. The exposure time for each observation ranged from 20 to 40 minutes, yielding an ${\rm{S}}/{\rm{N}}\sim 250$ per resolution element around 5500 Å.

On 2008 May 29, high-angular-resolution lucky images of HD 87646 were obtained with the FastCam instrument (Oscoz et al. 2008) on the 2.5 m Nordic Optical Telescope at the Roque de los Muchachos Observatory in La Palma (Spain). Five data cubes of 1000 images each were obtained using the I-band filter. Individual exposure times were 30 ms for each image. High-spatial-resolution was obtained by combining the best 1% of the images. Figure 1 shows the processed image. The image scale was 30.95 ± 0.05 mas pixel⁻¹. In this figure, the secondary star HD 87646B in this binary system is visible. The point-spread function (PSF) of the star has a full width half maximum (FWHM) value of 011.

Zoom In Zoom Out Reset image size

On 2009 June 4, we acquired high-resolution AO images of the binary star system HD 87646 from PALAO on the 200 inch Hale Telescope at Palomar. Data sets were taken in the J and K bands. The AO system was running at 500 Hz. The seeing was roughly 13 in the K band. We also observed a calibrator star and subtracted the K band's PSF, which improved sensitivity by a factor of a few. The utility of PSF subtraction is limited, in this case, by the difference in stellar spectral types, since the filters are broad.

Images were flat-fielded, background subtracted, and cleaned; the final images are displayed on a logarithmic scale in Figure 2 with a scale of 25 mas pixel⁻¹. The binary system is well-resolved in both bands. The angular separation is measured to be 401 ± 12 mas, nearly twice that quoted from the Hipparcos and Tycho Catalogs (Perryman 1997). The PA is $69\buildrel{\circ}\over{.} 8\pm 0\buildrel{\circ}\over{.} 5$. There is no evidence in the high-resolution images for a tertiary (stellar) companion. The J- and K-band brightness ratios of the two stars are 6.18 ± 0.12 and 5.82 ± 0.10, respectively. We cannot use these ratios and their error bars to put a meaningful constraint on the optical band flux ratio because the spectral energy distribution (SED) curve slope is basically flat in the J and K bands, but very sharp in the optical band.

Zoom In Zoom Out Reset image size

Photometry of HD 87646 was obtained in the Stromgren b and y bands between 2007 December and 2015 June with the T12 0.8 m APT at Fairborn Observatory in Arizona. Our primary goal with photometry is to detect if the companions transit the primary star. The data were processed using software described in Henry (1999). The measurements have a typical accuracy of ∼0.001 mag.

The SARG spectra taken at TNG without the iodine cell were used to derive the stellar parameters. HD 87646 is flagged as a binary in the Hipparcos and Tycho Catalogs (Perryman 1997), with a Hipparcos magnitude (a broadband V filter) difference between the primary and secondary to be 2.66 ± 0.96 mag. Taking into account the binary nature of the object, we explored possible combinations of stellar parameters for the primary and secondary. The SED and the colors change slightly due to the secondary contribution; however, the normalized spectra show minor changes, only affecting the wings of strong lines (e.g., the Mgb lines, Figure 3). The equivalent widths of most weak lines are essentially unchanged, so we used Fe i and Fe ii lines with equivalent widths below 140 mÅ and performed a traditional spectroscopic analysis.

Zoom In Zoom Out Reset image size

We use the latest MARCS model atmospheres (Gustafsson et al. 2008) for the analysis. Generation of synthetic spectra and the line analysis were performed using the turbospectrum code (Alvarez & Plez 1998), which employs line broadening according to the prescription of Barklem & O'Mara (1998). The line lists used are drawn from a variety of sources. Atomic lines are taken mainly from the VALD database (Kupka et al. 1999). The molecular species CH, CN, OH, CaH, and TiO are provided by B. Plez (see Plez & Cohen 2005), while the NH, MgH, and C₂ molecules are from the Kurucz line lists. The solar abundances used here are the same as Asplund (2005). We used the Fe i excitation equilibrium and derived an effective temperature ${T}_{\mathrm{eff}}$  = 5770 ±  80 K, which is slightly lower than the effective temperature derived from photometry (Feltzing & Gustafsson 1998). The Hα and Hβ wings also agree better for this lower ${T}_{\mathrm{eff}}$ value. We find log $g=4.1\pm 0.1$, based on the ionization equilibrium of Fe i and Fe ii lines and by fitting the wings of the Mgb lines, which agrees with the previous estimates. A microturbulence value of 1.8 km s⁻¹ is derived by forcing weak and strong Fe i lines to yield the same abundances. We are not able to confirm the super solar metallicity of this object (Feltzing & Gustafsson 1998); we derived [Fe/H] = −0.17 ± 0.08. When we adopt the same ${T}_{\mathrm{eff}}$ and microturbulence as Feltzing & Gustafsson (1998), we obtain the same metallicity value of [Fe/H] = 0.3, but with a large slope in the excitation potential versus the Fe i abundance and the reduced equivalent width versus the Fe i abundance. The derived stellar parameters are summarized in Table 2.

Table 2.  Parameters of the Star HD 87646A

Parameter	Value
${T}_{\mathrm{eff}}$	5770 ± 80 K
$\mathrm{log}(g)$	4.1 ± 0.1
$[\mathrm{Fe}/{\rm{H}}]$	−0.17 ± 0.08
$V\sin i$	7.5 km s−1
${\xi }_{t}$	1.8 km s−1
Mass	1.12 ± 0.09 M⊙
Radius	1.55 ± 0.22 R⊙

Download table as:  ASCIITypeset image

We attempted to place constraints on the secondary star by fitting the Balmer line and Mgb line wings. Based on the Hipparcos data, which suggest a Hipparcos magnitude difference between the primary and secondary 2.66 ± 0.96 mag, the secondary has a flux contribution of 8%–10% with respect to the primary. We synthesized binary spectra with a G-dwarf primary and a K-dwarf secondary, with 10% and 50% flux contributions from the secondary. The binary spectrum synthesis results are consistent with the Hipparcos data with a 10% contribution from the secondary. We cannot place better estimates based on the spectral line profiles because of the large errors in the continuum normalization of these lines, since they spread over the entire echelle order.

We use the empirical polynomial relations of Torres et al. (2010) to estimate the mass and radius of the primary star, HD 87646A, from ${T}_{\mathrm{eff}}$, $\mathrm{log}(g)$, and [Fe/H]. These relations were derived from a sample of eclipsing binaries with precisely measured masses and radii. We estimate the uncertainties in M_* and R_* by propagating the uncertainties in ${T}_{\mathrm{eff}}$, $\mathrm{log}(g)$, and [Fe/H] using the covariance matrices of the Torres et al. (2010) relations (kindly provided by G. Torres). Since the polynomial relations of Torres et al. (2010) were derived empirically, the relations were subject to some intrinsic scatter, which we add in quadrature to the uncertainties propagated from the stellar parameter measurements (${\sigma }_{\mathrm{log}m}=0.027$ and ${\sigma }_{\mathrm{log}r}=0.014;$ Torres et al. 2010). The final stellar mass and radius values obtained are ${M}_{* }=1.12\pm 0.09\,{M}_{\odot }$ and ${R}_{* }=1.55\pm 0.22\,{R}_{\odot }$.

HD 87646 is a binary system, and contamination of the primary star's spectrum by the secondary star leads to an increased RV jitter because it interferes with the analysis pipeline. In this section, we investigate the possible systematic RV errors caused by the blended binary spectra using simulations. Since our RV observations were produced by two different kinds of spectrographs, we decided to perform two simulations, one for the DFDI instruments, including KeckET, KPNO ET, and MARVELS, and the other for traditional echelle spectrographs, including HRS at HET and the AST fiber-fed echelle spectrograph at Fairborn. In both simulations, we first create a set of stellar spectra by combining a G-type star (for the primary) and a K-type (for the secondary) star spectra with varying radial velocities for both stars. Then we calculate the differential radial velocities for the G star from the simulated spectra. The differences between the output G-star RVs and input G-star RVs are the RV errors caused by secondary star spectra contamination. Both simulations yield similar RV errors on the order of 200 m s⁻¹. We expect to see this level of systematic error and will include it in the RV "jitter" term when we perform the RV curve fitting in the next section. Traditionally, the "jitter" term is used to denote any RV noise caused by stellar activity; our "jitter" term also contains the RV noise caused by blended binary spectra.

We have performed a Markov Chain Monte Carlo (MCMC) analysis of the combined RV observations from KeckET, ET, HET, MARVELS, and Fairborn instruments. In this analysis, we initially used a one planet RV model to fit our RV observations, and later found that there is another strong periodic RV signal in the RV residuals. We then adopted a two-object (a planet and a BD) RV model to fit our RV data. The RV model details are presented in Section 2 of Gregory (2007). We have attempted to add in another planet to fit our RV data around the third peak in the periodogram. The addition of another planet did not significantly improve our RV fitting. Furthermore, the fact that the newly added planet has a period half of the first giant planet and an eccentricity of 0.99, we consider it to be an alias and over-fit of the noise. Thus we have rejected the RV model with three substellar companions. Throughout the paper, we only present and discuss the RV model with two substellar companions.

Each state in the Markov chain is described by the parameter set

$$\begin{eqnarray}&&{\boldsymbol{\theta }}=\{{P}_{1},{K}_{1},{e}_{1},{\omega }_{1},{M}_{1},{P}_{2},{K}_{2},{e}_{2},{\omega }_{2},{M}_{2},{C}_{i},{\sigma }_{\mathrm{jitter}}\},\end{eqnarray} \tag{ 4 }$$

where P₁ and P₂ are orbital periods, K₁ and K₂ are the RV semi-amplitudes, e₁ and e₂ are the orbital eccentricities, ${\omega }_{1}$ and ${\omega }_{2}$ are the arguments of periastron, M₁ and M₂ are the mean anomalies at chosen epoch (τ), C_i is the constant velocity offset between the differential RV data shown in Table 1 and the zero-point of the Keplerian RV model (i = 1 for KeckET observations, i = 2 for KPNO ET observations, i = 3 for HET observations, i = 4 for MARVELS observations, and i = 5 for Fairborn observations), and ${\sigma }_{\mathrm{jitter}}$ is the "jitter" parameter. The jitter parameter describes any excess noise, including both astrophysical noise (e.g., stellar oscillation, stellar spots; Wright 2005), any instrument noise not accounted for in the quoted measurement uncertainties and systematic RV errors from analyzing blended binary spectra discussed in the last section. We use standard priors for each parameter (see Gregory 2007). The prior is uniform in the logarithm of the orbital period (P₁ and P₂) from 1 to 5000 days. For K₁, K₂, and ${\sigma }_{\mathrm{jitter}}$, we use a modified Jefferys prior, which takes the form of $p{(x)=(x+{x}_{o})}^{-1}{[\mathrm{log}(1+{x}_{\max }/{x}_{o}]}^{-1}$, where x_o = 0.1 ${\rm{m}}\,{{\rm{s}}}^{-1}$ and ${x}_{\max }=2128$ ${\rm{m}}\,{{\rm{s}}}^{-1}$ (Gregory 2005). Priors for e₁ and e₂ are uniform between zero and unity. Priors for ${\omega }_{1}$, ${\omega }_{2}$, M₁, and M₂ are uniform between zero and $2\pi$. For C_i, the priors are uniform between min(v_i) − 5 $\mathrm{km}\,{{\rm{s}}}^{-1}$ and max(v_i) + 5 $\mathrm{km}\,{{\rm{s}}}^{-1}$, where v_i is the set of radial velocities obtained from each of the four instruments (i = 1 for KeckET observations, i = 2 for KPNO ET observations, i = 3 for HET observations, i = 4 for MARVELS observation, and i = 5 for Fairborn observations). We verified that the chains did not approach the limiting values of P₁, P₂, K₁, K₂, and ${\sigma }_{\mathrm{jitter}}$.

Following Ford (2006), we adopt a likelihood (i.e., conditional probability of making the specified measurements given a particular set of model parameters) of

$$\begin{eqnarray}&&p(v| {\boldsymbol{\theta }},M)\propto \displaystyle \prod _{k}\displaystyle \frac{\exp [-{({v}_{k,\theta }-{v}_{k})}^{2}/2({\sigma }_{k,\mathrm{obs}}^{2}+{\sigma }_{\mathrm{jitter}}^{2})]}{\sqrt{{{\sigma }_{k,\mathrm{obs}}}^{2}+{{\sigma }_{\mathrm{jitter}}}^{2}}},\end{eqnarray} \tag{ 5 }$$

where v_k is observed RV at time t_k, ${v}_{k,\theta }$ is the model velocity at time t_k given the model parameters ${\boldsymbol{\theta }}$, and ${\sigma }_{k,\mathrm{obs}}$ is the measurement uncertainty for the RV observation at time t_k.

We combine the Markov chains described above to estimate the joint posterior probability distribution for the orbital model for HD 87646. In Table 3, we report the median value and an uncertainty estimate for each model parameter based on the marginal posterior probability distributions. The uncertainties are calculated as the standard deviation about the mean value from the combined posterior sample. Since the shape of the marginal posterior distribution is roughly similar to a multivariate normal distribution, the median value plus or minus the reported uncertainty roughly corresponds to a 68.3% credible interval. In the same table, we also reported the rms of the RV fitting residuals for the five different RV instruments used here, which are ${\mathrm{rms}}_{1}=245$ m s⁻¹ for KeckET, ${\mathrm{rms}}_{2}=248$ m s⁻¹ for KPNO ET, ${\mathrm{rms}}_{3}=261$ m s⁻¹ for HET, ${\mathrm{rms}}_{4}=270$ m s⁻¹ for MARVELS, and ${\mathrm{rms}}_{5}=312$ m s⁻¹ for Fairborn RV observations.

Table 3.  Orbital Parameters for HD 87646b and HD 87646c

Parameter	HD 87646b	HD 87646c
Minimum mass	12.4 ± 0.7 MJup	57.0 ± 3.7 MJup
a	0.117 ± 0.003 au	1.58 ± 0.04 au
K	956 ± 25 m s−1	1370 ± 54 m s−1
P	13.481 ± 0.001 day	674 ± 4 day
e	0.05 ± 0.02	0.50 ± 0.02
ω (radians)	5.20 ± 0.52	1.95 ± 0.06
${T}_{\mathrm{prediction}\mathrm{for}\mathrm{transit}}$ (JD${}_{\mathrm{UTC}}$)	2454093.85 ± 0.12 day
${T}_{\mathrm{periastron}}$ (JD${}_{\mathrm{UTC}}$)	2454088.3 ± 1.1 day	2453707 ± 9 day
${\sigma }_{\mathrm{jitter}}$	240 ± 12 m s−1
C1	20.878 ± 0.050 km s−1
C2	20.777 ± 0.090 km s−1
C3	0.908 ± 0.080 km s−1
C4	−0.306 ± 0.050 km s−1
C5	20.786 ± 0.070 km s−1
rms1	245 m s−1
rms2	248 m s−1
rms3	261 m s−1
rms4	270 m s−1
rms5	312 m s−1

Download table as:  ASCIITypeset image

This two-Keplerian orbital solution is shown in Figures 4 and 5 together with the KeckET, HET, KPNO ET, Fairborn, and MARVELS RV data. The residuals in these two plots cannot be explained only by the errors in our RV data. A stellar jitter term ${\sigma }_{\mathrm{jitter}}=240\pm 12$ m s⁻¹ is required in our fitting to explain these residuals. As discussed in the last section, most of the "jitter" noise arises from our data pipeline handling the blended binary spectra instead of a single star spectra. We also did an MCMC analysis using two "jitter" noise terms, one for DFDI instruments and the other one for echelle spectrographs, and find the orbital parameters for the giant planet candidate and the BD candidate are barely changed within the error bars. Thus, to keep it simple, we choose to use one "jitter" term for all our RV observations from different instruments.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

HD 87646 is a binary system, so we have done another MCMC analysis by including a linear RV trend (${v}_{\mathrm{trend}}\times (t-{t}_{0})$) to the two-object RV model used above. This linear RV trend is used to account for the perturbation of the primary star induced by the gravitational force of the secondary star. We note here that the offsets between different data sets will hinder the modeling of this linear trend as there is expected strong correlation between the offsets and this linear trend. Our new RV fitting yields orbital parameters for the two substellar companions in addition to a linear RV trend of ${v}_{\mathrm{trend}}=-12\pm 18$ m s⁻¹ yr, which are summarized in Table 4. Since all of the main orbital parameters of the two substellar companions are barely changed within their respective error bars and the strong correlation between this RV trend and telescope RV offsets, we decide to keep using the numbers present in Table 3 throughout this paper. This linear trend is not significant, which means it is more likely that either the secondary star is close to its ascending or descending node during 2008–2013, or this binary is on a relatively low-inclination (face on) orbit. It is not possible for us to distinguish these two scenarios using our current data. Future high-precision astrometry observations, such as GAIA, will help to solve this binary orbital problem. We also want to note here that this linear trend is not the exact real RV trend of the primary star induced by the gravitational perturbation of the secondary star because of the secondary star's spectral contamination. It is close to ∼70%–80% of the real trend value assuming the flux ratio of the two stars is ∼10 in the optical band and the mass ratio is ∼2.

Table 4.  Orbital Parameters for HD 87646b and HD 87646c with a Linear RV Trend

Parameter	HD 87646b	HD 87646c
Minimum mass	12.4 ± 0.7 MJup	57.0 ± 3.7 MJup
a	0.117 ± 0.003 au	1.58 ± 0.04 au
K	954 ± 24 m s−1	1370 ± 56 m s−1
P	13.481 ± 0.001 day	673 ± 4 day
e	0.05 ± 0.02	0.50 ± 0.02
ω (radians)	5.18 ± 0.48	1.95 ± 0.06
${T}_{\mathrm{prediction}\mathrm{for}\mathrm{transit}}$ (JD${}_{\mathrm{UTC}}$)	2454093.86 ± 0.14 day
${T}_{\mathrm{periastron}}$ (JD${}_{\mathrm{UTC}}$)	2454088.2 ± 1.0 day	2453709 ± 8 day
${v}_{\mathrm{trend}}$	−12 ± 18 m s−1 yr−1
${\sigma }_{\mathrm{jitter}}$	240 ± 13 m s−1
C1	20.89 ± 0.070 km s−1
C2	20.79 ± 0.10 km s−1
C3	0.92 ± 0.09 km s−1
C4	−0.262 ± 0.06 km s−1
C5	20.838 ± 0.09 km s−1
rms1	246 m s−1
rms2	248 m s−1
rms3	261 m s−1
rms4	269 m s−1
rms5	312 m s−1

Download table as:  ASCIITypeset image

Santos et al. (2002) found small RV variations and line asymmetries for the star HD 41004, which is a visual binary and is unresolved at the spectrograph. It was initially thought to have a planetary companion around the primary star, but from the line bisector analysis they were able to infer a possible BD orbiting the secondary star instead of a planet orbiting the primary star. Their conclusions were subsequently corroborated by Zucker et al. (2003).

Similar to HD 41004, HD 87646 is also a binary system with a small angular separation ($0\buildrel{\prime\prime}\over{.} 4$), which renders the spectrum a blended spectrum of the two stellar components. Following the same philosophy of Santos et al. (2002), we performed a bisector analysis for HD 87646 to determine from which star in the binary system the RV signal was produced. We have analyzed spectra taken at the kitt peak 2 m telescope using EXPERT (Ge et al. 2010). Spectra were reduced using an IDL pipeline modified from an early version described in Wang (2012). Frames were trimmed, bias subtracted, flat-field corrected, aperture-traced, and extracted. Cross-correlation functions (CCFS) are derived by cross-correlation with a spectral mask from the wavelength range of 4900–6300 Å. Then we compute the bisector velocity for 10 different levels of the CCF. The values for the upper (near continuum) and lower bisector points are averaged and subtracted. The resulting quantity (the bisector inverse slope, BIS) can be used to measure the line bisector variations (Queloz et al. 2001). The result of the bisector analysis is presented in Figure 6, which demonstrates that the BIS varies in phase with the RV.

Zoom In Zoom Out Reset image size

We then created two simulations following Santos et al. (2002), one by assuming a giant planet orbiting the primary star of the binary (solid line in Figure 6) and the other by assuming a heavier BD orbiting the secondary star (dotted line in Figure 6). Both scenarios can explain the RV curve seen for HD 87646, but clearly only the one in which the giant planet is orbiting the primary star (solid line in Figure 6) is consistent with the BIS analysis. Our conclusion from the BIS analysis is that the 13.5 day period giant planet is orbiting the primary star HD 87646A.

The top panel of Figure 7 presents all 1077 photometric observations plotted against the latest transit ephemeris of HD 87646a: ${T}_{{\rm{c}}}=2454093.85$ day, $P=13.481$ day. The differential magnitudes are measured against the mean of three comparison stars to improve precision. The standard deviation of these data from their mean is 0.0014 mag. A least-squares sine-curve fit to the phased data yields a full amplitude of 0.000089 ± 0.000065 mag. There is no detectable brightness variation on the short-period RV period; this result supports the interpretation that the observed RV variations are caused by a companion.

Zoom In Zoom Out Reset image size

The bottom panel of Figure 7 is similar to the top panel except that it displays only the data within ±0.1 phase units from the predicted transit time. We also show the predicted central transit, phased at 0.5, for a duration of 0.21 days or ∼0.015 units of phase and a depth of 0.5% or ∼0.005 mag (Kane & von Braun 2008). The $\pm 1\sigma$ uncertainty in the transit window timing is indicated by the two vertical dotted lines. There are 1005 observation that lie outside the predicted transit window, which have a mean of 0.99998 ± 0.00005 mag and 72 observations that fell within the transit window have a mean of 1.0003 ± 0.0002 mag. The difference between these two mean brightnesses is 0.0003 ± 0.0002 mag. Full transits with a predicted depthn of ∼0.005 mag are excluded by the photometry at the predicted transit time.

The mass function is related to the observed period, eccentricity, and RV semi-amplitude as

$$\begin{eqnarray}&&\displaystyle \frac{{(m\sin i)}^{3}}{{({M}_{* }+m)}^{2}}=\displaystyle \frac{P{(1-{e}^{2})}^{\tfrac{3}{2}}{K}^{3}}{2\pi G}\end{eqnarray} \tag{ 6 }$$

where M_* is the mass of the primary and m is the mass of the companions. Since the first companion is known not to transit the star, we cannot break the degeneracy of mass and $\sin i$ with RV observations alone. Using the derived stellar mass (1.12${M}_{\odot }$) for the primary with the orbital parameters determined from the RV (Table 3), we determine that the minimum mass of the inner companion for an edge-on orbit ($\sin i=1$) is 12.4 ± 0.7 M${}_{\mathrm{Jup}}$. This mass is quite close to the deuterium burning limit, and the detected companion is likely burning deuterium, though its minimum mass places it in the giant planet regime. The second companion's minimum mass, when assuming an edge-on orbit, is 57.0 ± 3.7 ${M}_{\mathrm{Jup}}$, which falls right in the BD regime.

Our SDSS MARVELS pilot survey and additional observations at the HET, KPNO 2.1 m telescope, and Fairborn observatory confirm the detection of two massive substellar companions in a close binary system HD 87646. The first companion, HD 87646Ab, has a minimum mass of 12.4 ± 0.7 ${M}_{\mathrm{Jup}}$, a period of 13.481 ± 0.001 days, and an eccentricity of 0.05 ± 0.02. The measured eccentricity is in line with other short-period giant planets in binaries (e.g., Eggenberger et al. 2007). This companion is likely to be a giant planet or a BD, depending on its inclination angle. Our bisector analysis has shown that this companion is in orbit around the primary star. The second companion has a minimum mass of 57 ±  3.7 ${M}_{\mathrm{Jup}}$, a period of 674 ± 4 days, and an eccentricity of 0.50 ± 0.02. This companion is likely to be a BD. Ma & Ge (2014) have found that long-period, high-mass (>42M${}_{\mathrm{Jup}}$) BDs tend to have higher eccentricities. This new BD is consistent with this trend.

This is the eleventh detection of a substellar companion(s) in a binary system with a separation of only about 20 au. The other 10 systems are Gliese 86 (Queloz et al. 2000; Mugrauer & Neuhäuser 2005; Lagrange et al. 2006), γ Cephei (Hatzes et al. 2003), HD 41004 (Zucker et al. 2003, 2004), HD 188753 (Konacki 2005), HD 176051 (Muterspaugh et al. 2010), HD 126614 (Howard et al. 2010), α Centauri (Dumusque et al. 2012), HD 196885 (Correia et al. 2008), OGLE-2013-BLG-0341 (Gould et al. 2014), and HD 59686 (Ortiz et al. 2016). However, Eggenberger et al. (2007) did not confirm the planet in HD 188753, and Rajpaul et al. (2016) suggest that the planet signal discovered from α Centauri B is not from a real planet, but from the observation window function. To the best of our knowledge, HD 87646A is the first multiple planet/BD system detected in such close binaries.

In this section, we will discuss the dynamical stability of the binary system HD 87646. First, we have collected observational data of HD 87646 from the literature (Horch et al. 2008, 2010; Hartkopf & Mason 2009; Balega et al. 2013) and combined them with our AO data to constrain the binary orbit of HD 87646. Our best-fitting binary orbital solution (solid line) and the observational data (black dots) are shown in Figure 8. The best-fitting parameters are $P=51.6\,\mathrm{year}$, e = 0.54, and a = 0.26 arcsec. At a distance of 73.58 ± 9.68 pc, this angular separation corresponds to a binary semimajor axis of 19 ± 2 au. Since the observational data only cover less than half of the binary orbit, these fitting parameters are very preliminary. There are many previous cases for which binary orbital parameters were revised significantly with new and better astrometry measurements, especially when the binary orbit is not yet completely covered by astrometry observations (Hartkopf & Mason 2009). Using the same classification of visual binaries as that in Hartkopf et al. (2001) for the Catalog of Visual Binary, our orbital solution has a grade of 4 (1—definitive, 2—good, 3—reliable, 4—preliminary, 5—indeterminate) and formal errors of the orbital solution were considered to have little meaning.

Zoom In Zoom Out Reset image size

We then performed a numerical simulation of the binary system including the giant planet and BD discovered in this paper. The binary-planetary-BD system of HD 87646 was integrated numerically using the Bulirsch–Stoer integrator of the N-body integration package Mercury (Chambers 1999). For each simulation, we tested the dynamical stability of this system given the semimajor axis (a${}_{{\rm{B}}}$) and eccentricity (e${}_{{\rm{B}}}$) of the binary system up to a million years. We have assumed the giant planet, BD, and binary to be coplanar. The results are shown in Figure 9. There is a stable zone in the binary ${a}_{{\rm{B}}}$–e${}_{{\rm{B}}}$ diagram. We have scaled the error bars of these astrometry data to force the best orbital fit to have a reduced chi-squared ${\chi }_{\mathrm{red}}^{2}=1$ and then over-plotted the binary orbital parameter fit from astrometry data with max ${\chi }_{\mathrm{red}}^{2}=2$ in Figure 9. The big uncertainty of the binary orbital fitting from astrometry data arises from the big error bar (∼0.1 arcsec) of the 1991 Hipparcos data point. We want to point out that there is another caveat of this plot, which is from our scaling of the error bars of all the astrometry data points. These data are from five different previous observation programs and scaling all of them together to make the best fit have a ${\chi }_{\mathrm{red}}^{2}=1$ can be problematic. The main conclusion from the simulation study and astrometry data fitting is, with a large binary semimajor axis (a${}_{{\rm{B}}}\gt 17$ au) and a relatively low binary eccentricity (${e}_{{\rm{B}}}\lt ({a}_{{\rm{B}}}-17)\times 0.57+0.2$), the binary-planet-BD system discovered in this paper is stable.

Zoom In Zoom Out Reset image size

HD 87646 is the first known system to have two massive substellar objects orbiting a star in a close binary system. Interestingly, the masses of these two substellar objects are close to the minimum masses for burning deuterium (∼13 ${M}_{\mathrm{Jup}}$, Spiegel et al. 2011) and hydrogen (∼75M${}_{\mathrm{Jup}}$, Chabrier & Baraffe 2000), which are generally assumed to be the general mass boundaries between the planet and the BD and between the BD and the star, respectively. All of these peculiarities raise a question: how could such a system be formed? Here we briefly discuss this intriguing issue.

The large masses of these two substellar objects suggest that they could be formed as stars with their binary hosts: a large molecular cloud collapsed and fragmented into four pieces; the larger two successfully became stars and formed the HD 87646 binary, and the other smaller ones failed to form stars and became the substellar objects in this system (Chabrier et al. 2014). This scenario might be relevant for the binary stars but seems to be problematic for the two substellar objects on orbits within ∼1 au because it is unclear whether fragmentation on such a small scale can occur (Kratter & Murray-Clay 2011).

Perhaps a more plausible explanation is that the two substellar objects were formed like a giant planet in a protoplanetary disk around HD 87646A. As for giant planet formation, there are currently two main models: core accretion versus disk instability (see the review by D'Angelo et al. 2010; Helled et al. 2014). Recently, many studies have examined the core accretion model's ability to form planets in close binaries with separations of ∼20 au (see the review by Thebault & Haghighipour 2014, and references therein). A commonly recognized issue is that the binary perturbations generally inhibit the growth of planetesimals in the disk (Thebault 2011). Even if their growth could proceed in certain favorable conditions (Xie & Zhou 2008, 2009), it would be significantly slowed, requiring a timescale of 10⁶ years or even longer (Xie et al. 2010). This result raises a problem for the formation of a gaseous giant planet because it would not form a planetary core (via planetesimal growth) to accrete gas before the gas disk dissipation, which takes a timescale as short as 10⁵–10⁶ years for such close binaries (Cieza et al. 2009; Kraus et al. 2012). Following the above logic, Xie et al. (2010) found that Jupiter-like planets are unlikely to form around Alpha Centauri B. As for the case of HD 87646, the formation of the two massive substellar objects via the core accretion model would be more problematic because it requires a more massive disk with a mass larger than 68 ${M}_{\mathrm{Jup}}$. Such a massive disk is seldom observed in close binaries, which indicates that should such a massive disk exist, it would dissipate much faster than a normal lighter one.

Conversely, the disk instability model could circumvent most of the above barriers. First, disk instability usually requires a very high disk mass, which is in line with the masses of the two detected substellar objects. Second, planet formation via disk instability requires a short timescale, which is also consistent with the short disk dissipation timescale observed in close binaries. In addition, the model of disk instability is recently advocated by Duchêne (2010), who argued that planet formation might be dominated by disk instability in close binaries based on the fact that exoplanets within close binaries (separation < 100 au) are significantly more massive than those within wide binaries or single stars. We could use the planet and BD mass to estimate the minimum surface density of the primordial disk, and test if such a disk is gravitationally unstable. We adopt the similarity solution of the evolving viscous disk (Hartmann et al. 1998) where the surface density follows R⁻¹ from the star to the disk edge. Since the binary separation is 19 au, the tidal truncation radius for the circumstellar disk around the primary star is ∼6 au (1/3 of the binary separation). Spreading the total mass of the 12.4 + 57 Jupiter mass to this 6 au disk, the disk surface density at 6 au is 2604 g cm⁻². At 6 au, the temperature is around 90 K with 1 solar luminosity. The sound speed is 0.56 km s⁻¹. Then the Toomre Q parameter with a 1.12 solar mass star is 1.5. At 1 au, the temperature is around 220 K. The surface density is 15610 g cm⁻². The sound speed is 0.87 km s⁻¹. The Toomre Q parameter is 5.7. Considering the primordial disk mass should be a lot more massive than the planet and BD mass, the disk is likely to be gravitationally unstable throughout the disk. This is consistent with gravitational instability leading to planet formation. Although several advantages exist for the disk instability model, we consider that such an explanation for the formation of the HD 87646 system should be taken with caution because whether disk instability can be triggered in the present of a close stellar companion remains an issue under debate (Nelson 2000; Mayer et al. 2005; Boss 2006).

Next, we want to discuss how the giant planet, b, moved to its current position with a very low eccentricity. Although the Kozai-Lidov mechanism and subsequent tidal dissipation (Wu & Murray 2003; Fabrycky & Tremaine 2007) have often been invoked to explain the formation of hot Jupiters, it is unlikely that the BD, c, has helped to move b inward because such a process cannot help to form "warm" Jupiters (Antonini et al. 2016). Thus, we prefer a scenario in which b initially formed in the disk and then migrated inward in the disk to its current position, which explains why it has a near-zero eccentricity. As for the BD, c, after it formed in the disk, scattering between c and other objects formed in the disk moved it to a higher eccentricity. During such a scattering process, lower mass objects tend to be ejected out and more massive objects are kicked inward with a higher eccentricity according to the simulation of Chatterjee et al. (2008). The stellar companion, B, cannot excite the eccentricity of c because Kozai oscillations will be destroyed by the presence of other massive substellar objects in the system (in this case, the giant planet b) according to the studies of Wu & Murray (2003) and Fabrycky & Tremaine (2007). However, the presence of the stellar companion, B, can help to enhance the scattering process between c and other objects formed in the protoplanetary disk around A (Marzari et al. 2005). Future astrometry observations, like those from Gaia (Perryman et al. 2001), can provide a better binary orbital solution and even possibly constrain the BD candidate's orbit. These data will help us study the dynamic structure of this complicated system and give more insight into its formation scenario.