MINING PLANET SEARCH DATA FOR BINARY STARS: THE ψ¹ DRACONIS SYSTEM

Several groups (e.g., Wittenmyer et al. 2006; Fischer et al. 2009; Pepe et al. 2011) have used the radial velocity method to search for planets around nearby stars for well over a decade, and have collectively uncovered several hundred planets to date. Close binary stars are usually cut from the star sample because they complicate the detection method (e.g., Bergmann et al. 2015), and because they have long been suspected to inhibit planet formation by quickly destroying (Kraus et al. 2012) or depleting (Harris et al. 2012) the planet-forming disk.

Previously unknown stellar binary companions are nonetheless still uncovered in planet-search data through large-amplitude linear trends or even full long-period orbits, but may be ignored since the goal is to find planet-mass companions. Since binary stars are usually excluded in the star sample, companions that are found tend to have extreme flux- and mass ratios. Binary stars with extreme mass-ratios on orbits with ∼10 year timescales are precisely the ones that are most difficult to detect and characterize with imaging techniques, and so they should not be ignored.

Several groups have recently worked toward using high-resolution spectroscopy to search for very faint companions to nearby stars, both in the context of detecting emission (Snellen et al. 2010; Gullikson & Endl 2013) or reflection (Martins et al. 2013) from "Hot Jupiter" planets, and in the context of detecting stellar binary systems with high contrast ratios (e.g., Gullikson & Dodson-Robinson 2013; Kolbl et al. 2015). Those groups all use a cross-correlation analysis to search directly for the spectral lines of the faint companion and mainly differ in their treatment of the primary star and telluric lines.

In this paper we use data from the McDonald Observatory Planet Search team to examine the ψ¹ Draconis system, which consists of an F5IV-V star (ψ¹ Dra A) orbited by a G0V star (ψ¹ Dra B) with angular separation 30'' (Mason et al. 2014). Tokovinin & Smekhov (2002) searched for signs of a spectroscopic companion to ψ¹ Dra A from 1991 to 1995, but found no radial velocity variation. More recently Toyota et al. (2009) noted a linear trend in their radial velocity measurements, and predicted a companion with M > 50M_J. Our data have a much longer time baseline than either of the previous studies, and show a significant fraction of the orbit which has recently reached quadrature. Furthermore, Endl et al. (2015) use adaptive-optics imaging to detect a ∼4500 K companion 155 mas from ψ¹ Dra A, which they hypothesize is the source of the orbital motion seen in the primary-star radial-velocity measurements.

Here, we use all of our spectra of ψ¹ Dra A to search directly for the spectral lines of the companion and measure the system mass ratio. We describe the observations and data reduction in Section 2, and the method we use to search for the companion in Section 3. Finally, we estimate the mass ratio of the system and give the parameters for the companion in Section 4.

All data were taken at the 2.7 m Harlan J. Smith Telescope at McDonald observatory using the 2dcoudé échelle spectrograph (Tull et al. 1995) at a resolving power $R\equiv \frac{\lambda }{{\rm{\Delta }}\lambda }=60,000.$ The starlight was filtered through a temperature-stabilized I₂ cell to imprint many sharp absorption lines on each spectrum to use for both a precise velocity metric (Butler et al. 1996) and to model the instrument profile (Endl et al. 2000). The raw CCD data were reduced with standard IRAF¹ tasks, and include steps for overscan trimming, bad-pixel processing, bias-frame subtraction, scattered-light removal, flat-field division, order extraction, and wavelength solution fitting using a Th–Ar calibration lamp spectrum. Particularly strong cosmic-ray hits were removed manually by interpolating across nearby pixels.

We used the Austral code (Endl et al. 2000) to measure the differential radial velocity of ψ¹ Dra A at each observation by comparing each spectrum to a high signal-to-noise ratio template spectrum of the same star. We provide the raw velocity measurements in Table 1, as well as the velocities shifted into the system velocity rest frame. The velocity shift necessary to convert from the differential radial velocities to that frame is found in Section 4. Table 1 also gives the measurements of the companion radial velocity (described in the next section).

We use a cross-correlation analysis inspired by recent work attempting to detect light from planetary companions around late-type stars (Gullikson & Endl 2013; Martins et al. 2013) to search for the companion (ψ¹ Dra C). We start by dividing all spectra by the blaze function of the spectrograph, and further divide them by an empirical I₂ cell absorption spectrum in the spectral orders with 500 < λ < 640 nm. The blaze function is derived by fitting a high-order polynomial to the extracted spectrum of an incandescent light source (a flat lamp), and the empirical I₂ spectrum is the spectrum of a flat lamp with the I₂ cell inserted in the light path. Both the flat lamp and I₂ spectra are observed each day of each observing run. We use the Telfit code (Gullikson et al. 2014) to fit and remove the unsaturated telluric absorption lines in the spectrum, and cross-correlate each residual spectrum against a Phoenix model spectrum (Husser & Ulbrich 2013) with the following parameters:

• 1.  
  T_eff = 4400 K;
• 2.  
  log g = 4.5 (cgs units);
• 3.  
  [Fe/H] = 0.0.

The model temperature was chosen on the basis of high-contrast imaging in Endl et al. (2015), which finds a companion with approximately that temperature. We shift each CCF so that the dominant peak, which signifies the match of the M-star model template with the F-type primary star, falls at v = 0. That effectively puts the cross-correlation functions in the rest frame of the primary star, although there is a constant velocity offset caused by small errors in the vacuum to air wavelength conversion and spectrograph wavelength drift throughout the night. We denote this shift as Δv₂ in later sections of this paper.

We normalize each CCF by subtracting a quadratic function that we fit well away from the peak, and then dividing by the height of the CCF at v = 0 (the peak). The average of the shifted CCFs is a close estimate for the cross-correlation function of the M-star template with the F5 primary star, since the contribution from the companion is diluted by shifting the CCFs to the primary star rest frame. We remove the contribution from the primary star by subtracting the average from each CCF. The result is a series of residual cross-correlation functions that are estimates for the CCF of the companion spectrum against the 4400 K model spectrum template, with significant noise. We show the residual CCFs in Figure 1; the trace of the companion star is easily visible as the dark curve near the top middle. We are unable to recover the companion signal at early dates when the two stars were close to one another in velocity space.

Zoom In Zoom Out Reset image size

We measure the radial velocity of the companion at each epoch by finding the maximum and FWHM of the residual CCF. On the basis of Figure 1, we use only the portion of the CCFs with −50 < v < 10 km s⁻¹; we show a typical residual CCF in Figure 2. Since the CCFs were shifted to subtract the contribution from the primary star, the measured velocities (v_m,2) are related to the true barycentric velocities (v₁ and v₂ for the primary and secondary, respectively) and the constant shift described above (Δv₂) through

$$\begin{eqnarray}&&{v}_{{\rm{m}},2}(t)={v}_{2}(t)-{v}_{1}(t)+{\rm{\Delta }}{v}_{2}.\end{eqnarray} \tag{ 1 }$$

Zoom In Zoom Out Reset image size

We give the measured companion velocities in Table 1 (column 5). We additionally provide the velocities in the system velocity rest-frame (v₂) by using Equation (1) and the results of the analysis described below. The uncertainties given in Table 1 are determined from the CCF peak width and the scaling factor (f) derived below. The shifted primary and secondary velocities given in Table 1 are for the reader's convenience since they are in the same reference frame; we use the raw measurements in the orbital fit.

We now use the radial-velocity measurements to find the best orbital parameters to describe the orbit, as well as some data scaling and shifting factors. The orbit is described by the semi-amplitudes for both the primary and secondary stars (K₁ and K₂, respectively), the longitude of pericenter (ω), the eccentricity (e), the period (P), and the periastron-passage epoch (T₀). We cannot measure the system radial velocity, which is usually the final orbital element, because the measured primary-star velocities are differential and the secondary-star velocities are measured relative to the primary star.

Since the primary-star radial velocity measurements are differential measurements, we must also fit a constant shift (Δv₁) to account for the absolute radial velocity of the primary at the time at which our template spectrum was observed. We include an rv-jitter term (σ_J) to the fit to account for radial-velocity variations not encompassed by the orbital solution, and add the value in quadrature with the formal uncertainties on the primary star-velocity measurements.

The companion radial velocities are measured relative to the primary star plus a small velocity shift (Δv₂) caused by slight inaccuracies in the vacuum-to-air wavelength conversion in the model spectrum and spectrograph wavelength drift throughout the night. Finally, the CCF peak full-width at half maximum vastly over-estimates the velocity uncertainty and so we fit a scale factor (f) to apply to the companion velocity uncertainties. The uncertainties given in Table 1 are already scaled by that factor. The full log-likelihood function $({\mathcal{L}})$ is then given by

$$\begin{eqnarray*}{s}_{1} & = & \displaystyle \sum _{{t}_{{\rm{m}}}}\displaystyle \frac{{({v}_{1,{\rm{m}}}-{v}_{1}({t}_{{\rm{m}}})-{\rm{\Delta }}{v}_{1})}^{2}}{{\sigma }_{{v}_{1}}^{2}+{\sigma }_{{\rm{J}}}^{2}}+\mathrm{ln}2\pi ({\sigma }_{{v}_{1}}^{2}+{\sigma }_{{\rm{J}}}^{2})\\ {s}_{2} & = & \displaystyle \sum _{{t}_{{\rm{m}}}}\displaystyle \frac{{({v}_{2,{\rm{m}}}-{v}_{2}({t}_{{\rm{m}}})(1+\displaystyle \frac{{K}_{1}}{{K}_{2}})-{\rm{\Delta }}{v}_{2})}^{2}}{f{\sigma }_{{v}_{2}}^{2}}+\mathrm{ln}2\pi (f{\sigma }_{{v}_{2}}^{2})\\ {\mathcal{L}} & = & -0.5({s}_{1}+{s}_{2}),\end{eqnarray*}$$

where ${v}_{\mathrm{1,2}}(t)=v({T}_{0},P,e,{K}_{\mathrm{1,2}},\omega ,t)$ is the velocity at time t given by the orbital elements T₀, P, e, K, and ω.

We use the affine invariant sampler provided in the emcee code (Foreman-Mackey et al. 2013) to perform a Markov Chain Monte Carlo (MCMC) fit to all of the parameters described above. We use flat priors in all variables except for the rv-jitter and companion rv uncertainty scale factors (σ_J and f, respectively), for which we use log-uniform priors to allow for a large range of values. We give the median value and uncertainty for each parameter in Table 2. The uncertainties are estimated from the posterior probability distribution samples such that the lower and upper bounds give the 16th and 84th percentile (i.e., they are 1σ credibility intervals). We plot the best-fit orbit with the data in Figure 3, with the uncertainties on the companion velocities scaled and the velocities shifted by Δv₁ and Δv₂.

Zoom In Zoom Out Reset image size

Table 2.  Orbital Parameters for the ψ¹ Draconis A Subsystem

Parameters From Endl et al. (2015)
Teff,1 (K)	6544 ± 42
log g	3.90 ± 0.11
[Fe/H]	−0.10 ± 0.05
Parameters Derived in This Work
K1 (km s−1)	${5.18}_{-0.03}^{+0.04}$
K2 (km s−1)	11.1 ± 0.2
Period (days)	${6774}_{-167}^{+271}$
Periastron passage time (JD)	${2450388}_{-273}^{+169}$
ω (°)	32.6 ± 0.7
e	${0.679}_{-0.004}^{+0.006}$
Δv1 (km s−1)	${4.10}_{-0.09}^{+0.06}$
Δv2 (km s−1)	$-{5.4}_{-0.2}^{+0.3}$
Companion uncertainty scale factor (f)	0.17 ± 0.02
rv jitter (m s−1)	${72}_{-5}^{+7}$
reduced χ2	0.41
q	0.466 ± 0.008
M1 (M⊙)	${1.38}_{-0.08}^{+0.15}$
M2 (M⊙)	0.70 ± 0.07
Teff,2 (K)	4400 ± 300
i (°)	31 ± 1
a (AU)	${9.1}_{-0.3}^{+0.4}$

Note. The primary mass is derived using the spectroscopic T_eff, log g, and [Fe/H] and interpolating Dartmouth isochrones. The companion temperature is likewise derived from the companion mass using Dartmouth isochrones of the same metallicity.

Download table as:  ASCIITypeset image

Next, we calculate a series of derived quantities to characterize the companion and report them in Table 2. The mass ratio of the system is the ratio K₁/K₂ = 0.47. We estimate the primary star mass by interpolating Dartmouth isochrones (Dotter et al. 2008) with the "isochrones" code (described in Montet et al. 2015), and using spectroscopic parameters derived in Endl et al. (2015). The secondary mass is ${M}_{2}={{qM}}_{1}\sim 0.70\ {M}_{\odot };$ assuming the same age and metallicity as the primary, the Dartmouth isochrones give an expected temperature of ∼4400 K. That temperature is in excellent agreement with the high-contrast-imaging data, which support a companion of ∼4400 K with large uncertainty. With both the primary and secondary star mass, we calculate the orbital inclination and semimajor axis and report them in Table 2.

We use nearly 15 years of time-series spectra of the star ψ¹ Draconis A to search for the spectral lines of a companion identified by a large-amplitude trend in the primary-star radial velocities and later by direct imaging. We cross-correlate each spectrum against a Phoenix model spectrum of a 4400 K star and subtract the average CCF. The residual CCFs clearly show the template match with the companion (Figure 1), and we are able to measure the companion radial velocities for most dates.

We use the radial-velocity measurements for both the primary and secondary stars to find an orbital solution for the now double-lined spectroscopic binary. The summary values of the fitted parameters are given in Table 2. Finally, we report the mass and expected temperature of the companion as well as the orbital inclination and semimajor axis. The temperature agrees well with high-contrast imaging, validating our method.

The ψ¹ Draconis system is therefore a hierarchical multiple system with the component parameters given in Table 2. ψ¹ Dra A and B are separated by ∼680 AU and have a mass ratio q = 0.82, while A and C (the new companion) have a much closer orbit with with a = 9.1 AU and q = 0.47.

This method could be used to search for the spectral lines of stellar companions to other stars observed with high-precision radial-velocity surveys. To that end, and in the goal of open science, we make the source code used for the analysis and generating the plots for this paper available at https://github.com/kgullikson88/Companion-Finder. The raw radial-velocity measurements for both the primary and secondary star, as well as the MCMC chains, are available at the same url.

This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France, and of Astropy, a community-developed core Python package for Astronomy (Robitaille et al. 2013). It was supported by a start-up grant to Adam Kraus from the University of Texas. The McDonald Observatory planet search is supported by the National Science Foundation under grant AST-1313075. We would like to thank the referee for various suggestions that improved this paper.