DIRECT SPECTRAL DETECTION: AN EFFICIENT METHOD TO DETECT AND CHARACTERIZE BINARY SYSTEMS

There has recently been a revival of interest in the multiplicity properties of intermediate-mass stars, spurred largely by the detection of planets orbiting nearby $\sim 2{M}_{\odot }\;$ stars on both wide (e.g., Marois et al. 2008; Lagrange et al. 2010) and close (Johnson et al. 2011) orbits. In this context, binary companions are contaminants; companions complicate radial velocity planet searches because they necessitate simultaneous modeling of both stellar motions (e.g., Bergmann et al. 2015). Likewise, companions complicate direct imaging planet searches by requiring either extremely high-contrast instrumentation (Thalmann et al. 2014) or specialized coronagraphs (Crepp et al. 2010).

However, known binary stars are typically avoided in planet search programs for a more fundamental reason: the binary companion depletes or destroys the planet-forming disk. By combining a binary census of the ∼2 Myr Taurus–Auriga star-forming region with a disk census of the same, Kraus et al. (2012) showed that close ($\lesssim 40\;{\rm{AU}}$) binaries are about 2–3 times less likely to host a protoplanetary disk, and so hasten disk dispersal. Even if a disk survives, it tends to be depleted in mass by a factor of ∼25 for binary separations $\lesssim 30\;{\rm{AU}}$ (Harris et al. 2012).

Multiplicity is an inevitable outcome of star formation, especially for more massive stars where multiplicity is more common (Zinnecker & Yorke 2007), and so is important to study in its own right. Beyond the overall multiplicity rate, the mass ratio, period, and eccentricity distributions of a binary star population can place important constraints on the mode of binary star formation. Specifically, a mass ratio distribution that changes with physical separation could indicate that secondary stars are either forming through disk fragmentation (e.g., Kratter & Matzner 2006; Stamatellos & Whitworth 2011) or are accreting a significant amount of their mass from the primary star disk. De Rosa et al. (2014) recently completed an imaging survey of nearby A-type stars, and found preliminary evidence that the mass ratio distribution does in fact become flatter for closer companions. There is no evidence of such a change for solar-type or later stars (Reggiani & Meyer 2013).

It is difficult to detect low-mass companions very near an intrinsically bright primary star and even more difficult to characterize the companion. There are three commonly used techniques for binary star searches: direct imaging with adaptive optics systems, interferometry, and radial velocity monitoring. Imaging can easily detect low-mass companions at wide apparent separations, but loses sensitivity as the on-sky distance from the primary star decreases (see De Rosa et al. (2014), for typical sensitivity curves). Interferometry can usually achieve smaller working angles than imaging, but cannot achieve as high cont1ast (see e.g., Aldoretta et al. 2015).

Radial velocity monitoring can easily find companions on very short-period orbits, but its sensitivity to low-mass companions drops as the physical separation increases. This is especially true for A- and B-type stars, where radial velocity precision is typically limited to $\sim 1\;\mathrm{km}\;{{\rm{s}}}^{-1}$ by their rotationally broadened lines. Worse, radial velocity monitoring techniques cannot characterize the companion unless the inclination is known or if the companion spectral lines are also visible. All three techniques have separation-dependent sensitivity, which introduces observational bias in any search for a parameter that changes with physical separation.

One technique that is separation independent is to search directly for the composite spectrum of two stars. Burgasser (2007) used single-epoch low-resolution spectroscopy to identify and characterize a brown dwarf binary system by fitting both spectra simultaneously. This method only works if the stars have a similar brightness but very different spectra, such that the spectral features from both components are visible and distinguishable. If, as in the case of binary systems with very large flux ratios, the companion spectrum is buried within the noise of the primary star, a different method is needed.

The direct spectral detection method (hereafter referred to as the DSD method), and variations thereof, has been used to search for binary companions to early B-stars (Gullikson & Endl 2013), main sequence FGK-stars (Kolbl et al. 2015), young K-M stars (Prato et al. 2002), and even "Hot Jupiter" type planets (Snellen et al. 2010; Brogi et al. 2012; de Kok et al. 2013) orbiting FGK-stars. The method relies on the cross-correlation function (CCF) of a high-spectral-resolution spectrum of the primary star with a model spectrum for the expected companion. The CCF uses every pixel in the spectrum, and more importantly every spectral line in the secondary spectrum. A simple experiment with synthetic spectra containing increasing numbers of spectral lines (N) in noisy data shows that the CCF peak significance increases as $\sim \sqrt{N}$. For high-resolution cross-dispersed échelle spectra, this amplification can reach several factors of 10, allowing the detection and characterization of a secondary spectrum where the individual lines are completely buried in noise. Since the DSD method uses a seeing-limited spectrum of the primary star, its sensitivity is independent of separation inside ∼1'' and can make use of small telescopes to detect high-contrast companions.

In this paper, we describe the DSD method in detail, and use it to detect the secondary star in nine known binary systems. We describe the method in Section 2. We describe the observations and data reduction in Section 3, then use the observations to estimate the accuracy with which we can measure the companion temperature in Section 4 and the sensitivity of the method in Section 5. In Section 6 we use the DSD method to search for known companions, and discuss the results in Section 7.

All implementations of the DSD method use high-spectral-resolution and high signal-to-noise spectra, and search for companions with extreme flux ratios by cross-correlating the observed spectra with models for the expected companion. The main differences between the various implementations are the primary star and telluric line removal processes. The "Hot Jupiter" searches (e.g., Snellen et al. 2010) use known orbital phase information and a high degree of phase coverage to simultaneously estimate an empirical stellar and telluric spectrum with minimal contamination from the planet, while Kolbl et al. (2015) subtract a best-fit model spectrum for the primary star. Since Kolbl et al. (2015) use optical data, they do not attempt any telluric correction. Note that the approach of Kolbl et al. (2015) is conceptually similar to the todcor code (Mazeh & Zucker 1994), which is widely used to search for double-lined spectroscopic binary systems.

Unlike most previous work, we focus on not only detecting but accurately characterizing the companion. We additionally optimize our technique for detecting cool companions to rapidly rotating early-type stars, for which it is very difficult to detect the reflex motion of the primary star. We fit and remove the telluric absorption using the TelFit code (Gullikson et al. 2014), and estimate an empirical primary star spectrum with a Gaussian smoothing filter applied to the telluric-corrected data. We chose to use a smoothing filter over subtracting model spectra for two reasons: first, the model spectra are a poor representation of the data, especially at the high signal-to-noise ratios that we use, and so leave very large-scale features in the residual spectrum. Second, the smoothing filter will also remove any large-scale instrumental systematics in the spectrum. We use a smoothing filter with a window size (w) set by

$$\begin{eqnarray}&&w=\displaystyle \frac{v\mathrm{sin}i\cdot f}{c}\cdot \displaystyle \frac{{\lambda }_{0}}{{\rm{\Delta }}\lambda }\end{eqnarray} \tag{ 1 }$$

where $v\mathrm{sin}i$ is the literature rotational velocity of the star, ${\lambda }_{0}$ is the central wavelength of the échelle order, ${\rm{\Delta }}\lambda$ is the wavelength spacing per pixel of the order, c is the speed of light, and f = 0.25 is an empirically determined parameter to give a visually adequate fit. Typical window sizes ranged from 50–100 pixels.

We use the following subset of the Phoenix library of model spectra prepared by Husser et al. (2013) throughout this work:

• 1.  
  ${T}_{\mathrm{eff}}=3000-7000\;{\rm{K}}$ ⁴ , in steps of 100 K
• 2.  
  [Fe/H] = −0.5, 0.0, +0.5
• 3.  
  $v\mathrm{sin}i=1,5,10,20,30\ \mathrm{km}\ {{\rm{s}}}^{-1}.$

Here, the $v\mathrm{sin}i$ is the rotational velocity of the secondary star. We account for the small influence that the smoothing kernel has on the companion spectrum by convolving the model with the same smoothing kernel used for the data, and subtracting the convolved model from the original. Treating the model spectrum in this way is more commonly known as unsharp masking. Finally, we cross-correlate every échelle order that does not have strong telluric residuals against the corresponding model spectrum, and combine the CCFs for each order using a simple average. The method is summarized below.

• 1.  
  Smooth the observed spectrum with a Gaussian smoothing kernel with width given by Equation (1), and subtract the smoothed spectrum from the original.
• 2.  
  Rotationally broaden the Phoenix model spectrum to the requested companion $v\mathrm{sin}i.$
• 3.  
  Smooth the broadened spectrum to the instrumental resolution by convolving it with a Gaussian kernel of appropriate width.
• 4.  
  Unsharp mask the broadened model spectrum.
• 5.  
  Resample the processed model spectrum to the same wavelength spacing per pixel as the observed spectrum.
• 6.  
  Cross-correlate each échelle order against the corresponding processed model spectrum, and combine using a simple average.

We use three separate samples in this work. The first set, given in Table 1, contains A- and B-type stars with the following published properties:

• 1.  
  Spectral type B0V–A9V (only main sequence)
• 2.  
  $V\lt 6$
• 3.  
  $v\mathrm{sin}i\gt 80\mathrm{km}\ {{\rm{s}}}^{-1}$
• 4.  
  No known companions within 3'', and no sign of a companion in our data.

The lower limit on $v\mathrm{sin}i$ in our sample ensures that the empirical primary star template is accurate and only minimally affects any companions.

The second data set (Table 2) contains F-M type stars which have a high-quality temperature estimate in the literature. We use the first two samples in Sections 4 and 5 to assess the accuracy of the temperature estimation using the DSD method and the sensitivity to companions of various temperatures.

Table 2.  Late Type Calibration Stars

Star	R.A.	Decl.	V	K	${T}_{\mathrm{eff}}$ (K)	Instrument	Date	Exp. Time (minute)
HD 33793	05:11:40.5	−45:01:06.2	8.85	5.05	3570 ± 1601	CHIRON	2015 Jan 13	60.00
HD 36379	05:30:59.9	−10:04:51.9	6.91	5.56	6030 ± 142	CHIRON	2015 Jan 14	9.58
HD 38858	05:48:34.9	−04:05:40.7	5.97	4.41	5646 ± 453	CHIRON	2015 Jan 14	5.31
HD 42581	06:10:34.6	−21:51:52.7	8.12	4.17	3814 ± 1134	CHIRON	2015 Jan 14	30.62
HD 45184	06:24:43.8	−28:46:48.4	6.39	4.87	5869 ± 142	CHIRON	2015 Jan 14	5.30
HD 50806	06:53:33.9	−28:32:23.2	6.04	4.33	5633 ± 152	CHIRON	2015 Jan 14	3.99
HD 61421	07:39:18.1	+05:13:29.9	0.37	−0.65	6582 ± 163	CHIRON	2015 Jan 16	0.05
HD 69830	08:18:23.9	−12:37:55.8	5.95	4.16	5402 ± 282	CHIRON	2015 Jan 14	5.03
HD 102634	11:49:01.2	−00:19:07.2	6.15	4.92	6215 ± 445	CHIRON	2015 Jan 17	5.18
GJ 465	12:24:52.5	−18:14:32.2	11.27	6.95	3472 ± 1106	CHIRON	2015 Jan 17	65.00
HD 115617	13:18:24.3	−18:18:40.3	4.74	2.96	5558 ± 192	CHIRON	2015 Jan 17	1.32
HD 125072	14:19:04.8	−59:22:44.5	6.66	4.33	4903 ± 445	CHIRON	2015 Feb 11	9.14
HD 128621	14:39:35.0	−60:50:15.0	1.33	−0.60	5232 ± 83	CHIRON	2015 Feb 06	0.03
HD 154363	17:05:03.3	−05:03:59.4	7.71	4.73	4723 ± 892	CHIRON	2015 Mar 12	26.27
HD 157881	17:25:45.2	+02:06:41.1	7.56	4.14	4124 ± 607	CHIRON	2015 Mar 13	25.86
HD 165222	18:05:07.5	−03:01:52.7	9.36	5.31	3416 ± 407	CHIRON	2015 Feb 11	3.83
HD 225239	00:04:53.7	+34:39:35.2	6.11	4.44	5699 ± 808	HRS	2002 Sep 18	8.00
HD 3651	00:39:21.8	+21:15:01.7	5.88	4.00	5046 ± 863	HRS	2005 Jul 30	3.00
HD 16895	02:44:11.9	+49:13:42.4	4.11	2.78	6344 ± 445	HRS	2006 Dec 02	0.11
HD 38529	05:46:34.9	+01:10:05.4	5.94	4.21	5697 ± 445	HRS	2004 Dec 02	0.55
GJ 270	07:19:31.2	+32:49:48.3	10.05	6.38	3668 ± 549	HRS	2002 Dec 11	20.00
HD 58855	07:29:55.9	+49:40:20.8	5.36	4.18	6398 ± 808	HRS	2006 Mar 12	0.50
GJ 281	07:39:23.0	+02:11:01.1	9.59	5.87	3776 ± 14510	HRS	2003 Jan 19	20.00
HD 69056	08:15:33.2	+11:25:51.4	7.70	6.06	5635 ± 558	HRS	2003 Dec 02	13.00
HD 73732	08:52:35.8	+28:19:50.9	5.95	4.01	5235 ± 445	HRS	2003 Oct 15	3.33
GJ 328	08:55:07.5	+01:32:56.4	9.98	6.35	3828 ± 16810	HRS	2003 Jan 14	20.00
HD 79969	09:17:53.4	+28:33:37.8	7.21	4.77	4825 ± 811	HRS	2003 Dec 02	10.00
HIP 53070	10:51:28.1	+20:16:38.9	8.22	6.83	6110 ± 768	HRS	2009 Feb 14	20.00
HIP 53169	10:52:36.4	−02:06:33.5	9.82	7.05	4525 ± 4712	HRS	2009 Jan 09	15.00
GJ 411	11:03:20.1	+35:58:11.5	7.52	3.34	3464 ± 1513	HRS	2001 Dec 27	5.00
HD 114783	13:12:43.7	−02:15:54.1	7.55	5.47	5135 ± 445	HRS	2005 Jan 08	7.08
GJ 525	13:45:05.0	+17:47:07.5	9.75	6.22	3680 ± 15014	HRS	2008 Apr 21	15.00
GJ 535	13:59:19.4	+22:52:11.1	9.04	6.24	4580 ± 711	HRS	2002 Apr 29	12.16
HD 142267	15:53:12.0	+13:11:47.8	6.12	4.53	5756 ± 445	HRS	2002 Aug 11	2.08
GJ 687	17:36:25.8	+68:20:20.9	9.15	4.55	3413 ± 2813	HRS	2002 Apr 30	12.50
GJ 699	17:57:48.4	+04:41:36.2	9.51	4.52	3222 ± 1013	HRS	2002 May 25	15.00
GJ 699	17:57:48.4	+04:41:36.2	9.51	4.52	3222 ± 1013	HRS	2002 May 25	35.00
GL 15A	00:18:22.8	+44:01:22.6	8.13	4.02	3567 ± 1113	IGRINS	2014 Nov 23	4.00
GL 15B	00:18:25.4	+44:01:37.6	11.04	5.95	3218 ± 607	IGRINS	2014 Nov 23	8.00
HD 1835	00:22:51.7	−12:12:33.9	6.39	4.86	5837 ± 445	IGRINS	2014 Dec 07	8.00
HD 4614	00:49:06.2	+57:48:54.6	3.44	1.99	5973 ± 83	IGRINS	2014 Dec 06	0.67
HD 10476	01:42:29.7	+20:16:06.6	5.24	3.25	5242 ± 123	IGRINS	2014 Nov 18	2.67
GL 1094	07:02:42.9	−06:47:57.2	8.35	5.76	4698 ± 912	IGRINS	2014 Nov 24	6.00
HD 58946	07:29:06.7	+31:47:04.3	4.18	2.98	6597 ± 183	IGRINS	2015 Jan 20	0.83
HD 67767	08:10:27.1	+25:30:26.4	5.73	3.84	5344 ± 445	IGRINS	2015 Jan 20	1.50
HD 71148	08:27:36.7	+45:39:10.7	6.30	4.83	5818 ± 445	IGRINS	2015 Jan 20	4.67
HD 76151	08:54:17.9	−05:26:04.0	6.00	4.46	5788 ± 232	IGRINS	2014 Nov 23	8.00
HD 87141	10:04:36.3	+53:53:30.1	5.72	4.50	6401 ± 808	IGRINS	2015 Jan 23	4.67
HD 87822	10:08:15.8	+31:36:14.5	6.24	5.13	6586 ± 808	IGRINS	2015 Jan 23	26.67
HD 91752	10:36:21.4	+36:19:36.9	6.30	5.20	6543 ± 808	IGRINS	2015 Jan 20	24.00
HD 95128	10:59:27.9	+40:25:48.9	5.04	3.75	5882 ± 445	IGRINS	2015 Jan 23	4.00
HD 95735	11:03:20.1	+35:58:11.5	7.52	3.34	3464 ± 1513	IGRINS	2015 Jan 23	2.00
BS 5019	13:18:24.3	−18:18:40.3	4.74	2.96	5558 ± 192	IGRINS	2015 Jan 06	6.00
HD 119850	13:45:43.7	+14:53:29.4	8.50	4.41	3618 ± 3113	IGRINS	2015 Jan 27	2.00
HD 122120	13:59:19.4	+22:52:11.1	9.04	6.24	4580 ± 711	IGRINS	2015 Jan 27	6.00
HD 122652	14:02:31.6	+31:39:39.0	7.15	5.88	6093 ± 445	IGRINS	2015 Jan 27	4.00
GJ 570A	14:57:28.0	−21:24:55.7	5.72	3.10	4507 ± 5813	IGRINS	2014 May 27	1.33
GJ 576	15:04:53.5	+05:38:17.1	9.81	6.47	4450 ± 10015	IGRINS	2015 Jan 27	6.00
GJ 758	19:23:34.0	+33:13:19.0	6.36	4.49	5453 ± 445	IGRINS	2014 Oct 10	3.00
GJ 820 A	21:06:53.9	+38:44:57.9	5.21	2.68	4361 ± 1713	IGRINS	2014 Dec 05	0.67
GJ 820 B	21:06:55.2	+38:44:31.4	6.03	2.32	3932 ± 2513	IGRINS	2014 Dec 05	0.67
HD 220339	23:23:04.8	−10:45:51.2	7.80	5.59	5029 ± 522	IGRINS	2014 Dec 07	13.00
HIP 117473	23:49:12.5	+02:24:04.4	8.99	5.04	3646 ± 607	IGRINS	2014 Nov 24	4.00
HD 4614	00:49:06.2	+57:48:54.6	3.44	1.99	5973 ± 83	TS23	1998 Jul 16	2.50
HD 10700	01:44:04.0	−15:56:14.9	3.50	1.68	5290 ± 393	TS23	1998 Jul 16	3.00
GJ 74	01:46:38.7	+12:24:42.3	8.89	6.32	4638 ± 7212	TS23	2008 Apr 12	20.00
HD 22049	03:32:55.8	−09:27:29.7	3.73	1.67	5077 ± 3513	TS23	2000 Sep 22	1.67
HR 1287	04:10:49.8	+26:28:51.4	5.40	4.48	6912 ± 808	TS23	2008 Mar 30	5.00
HD 30652	04:49:50.4	+06:57:40.5	3.19	2.05	6414 ± 193	TS23	1998 Nov 03	1.00
HD 40590	05:59:51.5	+00:03:21.4	8.07	6.91	6528 ± 758	TS23	2004 Feb 03	21.67
HR 3538	08:54:17.9	−05:26:04.0	6.00	4.46	5788 ± 232	TS23	2000 Jan 15	15.00
GJ 380	10:11:22.1	+49:27:15.2	6.61	3.26	4085 ± 1413	TS23	2012 Oct 02	13.33
GJ 411	11:03:20.1	+35:58:11.5	7.52	3.34	3464 ± 1513	TS23	2008 Mar 27	10.00
61 Vir	13:18:24.3	−18:18:40.3	4.74	2.96	5558 ± 192	TS23	2000 Jan 12	12.00
70 Vir	13:28:25.8	+13:46:43.6	4.97	3.24	5406 ± 643	TS23	1998 Jul 14	8.00
HD 142860	15:56:27.1	+15:39:41.8	3.84	2.62	6222 ± 133	TS23	1998 Jul 14	2.50
GJ 699	17:57:48.4	+04:41:36.2	9.51	4.52	3222 ± 1013	TS23	2000 May 24	35.00
70 Oph A	18:05:27.3	+02:29:59.3	4.20	1.79	5407 ± 5213	TS23	1998 Jul 14	3.00
16 Cyg A	19:41:48.9	+50:31:30.2	5.95	4.43	5750 ± 573	TS23	2005 Oct 12	6.67
16 Cyg B	19:41:51.9	+50:31:03.0	6.20	4.65	5678 ± 663	TS23	2002 Sep 21	13.33
61 Cyg B	21:06:55.2	+38:44:31.4	6.03	2.32	3932 ± 2513	TS23	1998 Jul 14	10.00
GJ 864	22:36:09.6	−00:50:30.0	9.92	6.16	3916 ± 617	TS23	2002 Nov 22	25.00
HD 216625	22:54:07.4	+19:53:31.3	7.02	5.73	6212 ± 445	TS23	2001 Jul 25	20.00

Note. The temperatures come from the following sources, and are labeled as superscripts after the temperature. (1) Woolf & Wallerstein (2005), (2) Sousa et al. (2008), (3) Boyajian et al. (2013), (4) Alonso et al. (1996), (5) Valenti & Fischer (2005), (6) Neves et al. (2014), (7) Mann et al. (2015), (8) Casagrande et al. (2011), (9) Ramírez and Meléndez (2005), (10) Casagrande et al. (2008), (11) Mishenina et al. (2012), (12) Casagrande et al. (2010), (13) Boyajian et al. (2012), (14) Pecaut & Mamajek (2013), (15) Zboril & Byrne (1998).

Download table as:  ASCIITypeset images: 1 2

Finally, we use the third data set (Table 3) to search for the spectral lines of the companion in several known binary systems. The third data set has the same properties as the first, except that they have one and only one known companion within 1''. We further require that the literature data either puts no constraints on the companion temperature (as in the case of single-lined spectroscopic binaries) or that the companion has ${T}_{{\rm{eff}}}\lt 7500\;{\rm{K}}$.

We estimate the expected companion temperature depending on whether it is part of a spectroscopic (Table 4) or visual (Table 5) binary system. In the case of double-lined spectroscopic binaries, we use the ratio of the semi-amplitudes given in the 9th catalog of spectroscopic binary orbits (SB9, Pourbaix et al. 2004) to estimate the mass ratio of the system. We convert the primary star spectral type from the Simbad database (Wenger et al. 2000) to mass by interpolating Table 5 of Pecaut & Mamajek (2013). The mass ratio and primary mass gives an estimate of the companion mass, which we convert to temperature by interpolating the same table. Most of the directly imaged binary systems do not have orbital information, so we use the magnitude difference published in the Washington Double Star catalog (WDS, Mason et al. 2014). We use the Simbad spectral type of the primary star and Table 5 of Pecaut & Mamajek (2013) to estimate the primary star temperature (T₁) and radius (R₁). We then find the companion temperature that minimizes the following function for the companion temperature (T₂), given the observed magnitude difference (${\rm{\Delta }}{m}_{\mathrm{obs}}$)

$$\begin{eqnarray}&&Q={(m({T}_{2},{R}_{2})-m({T}_{1},{R}_{1})-{\rm{\Delta }}{m}_{\mathrm{obs}})}^{2}\end{eqnarray} \tag{ 2 }$$

where $m(T,R)$ is the Vega magnitude of a star with temperature T and radius R. We use the pysynphot package⁵ and a Kurucz model grid (Castelli & Kurucz 2003) to calculate $m(T,R)$, and assume the companion is on the main sequence to estimate its radius (R₂). For both spectroscopic and visual binary systems, we assume spectral type uncertainties of ±1 subtype on the primary stars, and propagate the uncertainties into uncertainty in the companion temperature. We include the binary system in the sample if the expected companion temperature is $\lt 7500\;{\rm{K}}$.

Table 4.  Literature Spectroscopic Data

	K1	K2	Period
Star	($\mathrm{km}\;{{\rm{s}}}^{-1}$)	($\mathrm{km}\;{{\rm{s}}}^{-1}$)	(days)
HIP 33004	11.90	⋯	940.20
HIP 127194	8.80	⋯	490.00
HIP 131655	24.80	⋯	3.85
HIP 153386	20.00	⋯	36.50
HIP 175637	26.80	⋯	1.69
HIP 228407	24.50	⋯	24.10
HIP 2606311	13.50	⋯	119.09
HIP 2656312	28.60	⋯	445.74
HIP 2869113	12.22	⋯	4741.10
HIP 4412712	6.00	⋯	4028.00
HIP 5859012	26.20	⋯	282.69
HIP 7626714	35.40	99.00	17.36
HIP 7785815	32.90	⋯	1.92
HIP 7940415	31.50	⋯	5.78
HIP 853857	17.10	⋯	8.96
HIP 9272816	39.70	⋯	88.35
HIP 10022118	49.70	⋯	5.30
HIP 10678612	11.30	⋯	8016.00
HIP 1166113	25.19	⋯	0.50

Note. Known binary stars with spectroscopic orbit solutions. The orbital data is from the SB9 database (Pourbaix et al. 2004), and the original references are provided as superscripts after the star names: (1) Hill et al. (1971), (2) Lloyd (1981), (3) Rucinski et al. (2005), (4) Abt & Levy (1978) (5) Pourbaix et al. (2004), (6) Morrell & Abt (1992), (7) Abt et al. (1990), (8) Fekel & Tomkin (1982), (9) Lucy & Sweeney (1971), (10) Pogo (1928), (11) Duerbeck (1975), (12) Abt (1965), (13) Scarfe et al. (2000), (14) Tomkin & Popper (1986), (15) Levato et al. (1987), (16) Richardson & McKellar (1957), (17) Leone & Catanzaro (1999), (18) Hube (1973), (19) Pearce (1936).

Download table as:  ASCIITypeset image

Table 5.  Literature Imaging Data

	Separation		Wavelength
Star	('')	${\rm{\Delta }}m$	(nm)
HIP 13661	0.06	⋯	549
HIP 229583	0.65	4.15 ± 0.14	511
HIP 249023	0.38	2.97 ± 0.06	511
HIP 333723	0.75	3.27 ± 0.04	511
HIP 654774	1.11	5.18 ± 0.07	4770
HIP 775165	0.20	1.70 ± 0.05	780
HIP 791996	1.12	4.62 ± 0.12	1250
HIP 816417	0.04	1.90 ± 0.00	551
HIP 846063	0.84	4.02 ± 0.08	511
HIP 882903	0.58	2.95 ± 0.04	511
HIP 911188	0.16	⋯	549
HIP 920279	0.17	1.53 ± 0.00	550
HIP 9805510	0.10	0.51 ± 0.00	550
HIP 1137883	0.39	2.17 ± 0.02	511
HIP 11624711	0.84	2.43 ± 0.15	541
HIP 11661112	0.95	5.93 ± 0.09	2169

Notes. Known binary stars detected through either high-contrast imaging or interferometry. The imaging data comes from the Washington Double Star Catalog (Mason et al. 2014), and the most recent measurements are given as superscripts to the star name: (1) McAlister et al. (1989), (2) Roberts et al. (2007), (3) ESA (1997), (4) Mamajek et al. (2010), (5) Drummond (2014), (6) Shatsky & Tokovinin (2002), (7) Tokovinin et al. (2010), (8) McAlister et al. (1987), (9) Horch et al. (2010), (10) Horch et al. (2008), (11) Horch et al. (2001), (12) De Rosa et al. (2012).

Download table as:  ASCIITypeset image

We use the same set of instruments and settings for all observations throughout the three data sets. We use the CHIRON spectrograph (Tokovinin et al. 2013) on the 1.5 m telescope at Cerro Tololo Inter-American Observatory for most southern targets. This spectrograph is an $R\equiv \lambda /{\rm{\Delta }}\lambda =80000$ cross-dispersed échelle spectrograph with wavelength coverage from 450 to 850 nm, and is fed by a 27 optical fiber. The data are automatically reduced with a standard CHIRON data reduction pipeline, but the pipeline leaves residuals of strong lines in adjacent orders. We therefore bias-correct, flat-field and extract the spectra with the optimum extraction technique (Horne 1986) using standard IRAF⁶ tasks, and use the wavelength calibration from the pipeline reduced spectra.

For the northern targets, we use a combination of the High Resolution Spectrograph (HRS, Tull 1998) on the Hobby Eberly Telescope, and the Tull coudé (TS23, Tull et al. 1995) and IGRINS (Park et al. 2014) spectrographs, both on the 2.7 m Harlan J. Smith Telescope. All three northern instruments are at McDonald Observatory. For the HRS, we use the R = 60000 setting with a 2'' fiber, and with wavelength coverage from 410 to 780 nm. We bias-correct, flat-field, and extract the spectra using an IRAF pipeline. The HRS spectra are wavelength-calibrated using a Th–Ar lamp observed immediately before or after the science observations.

For the TS23, we use a 12 slit in combination with the E2 échelle grating (53 grooves/mm, blaze angle 65^◦), yielding a resolving power of R = 60000 and a wavelength coverage from 375 to 1020 nm. We reduce the data using an IRAF pipeline very similar to the one we use for the HRS, and wavelength calibrate using a Th–Ar lamp observed immediately before the science observations.

IGRINS only has one setting with R = 40000. It has complete wavelength coverage from 1475 to 2480 nm, except for where telluric absorption is almost 100% from 1810 to 1930 nm. Each star is observed in an ABBA nodding mode, and reduced using the standard IGRINS pipeline (Lee 2015). The standard pipeline uses atmospheric OH emission lines as well as a Th–Ar calibration frame to calibrate the wavelengths; we further refine the wavelength solution using telluric absorption lines and the TelFit code.

In the absence of noise, the CCF of an observed spectrum with a perfect model will have a value of 1 at the radial velocity of the star. As the model becomes a worse representation of the data, the peak height of the resulting CCF will decrease. Thus the CCFs act in a similar way as a ${\chi }^{2}$ map of the parameter space, allowing us to measure the effective temperature, metallicity, and rotational broadening of the secondary star. However, the presence of noise and the imperfections in the model spectra cause the measured values to deviate from the true parameters of the secondary star.

To measure the impact of both random and systematic noise on the parameter estimation, we created several hundred synthetic binary systems for each instrument used in our program. We made the synthetic binary systems by combining the early-type star spectra from Table 1 with those of the late-type stars in Table 2 in every possible combination, provided both observations came from the same instrument. By combining actual observations of early-type and late-type stars, our synthetic binary observations retain any instrument-specific effects that may impact the temperature estimation. We scaled the flux of the late-type star such that the flux ratio (${F}_{{\rm{secondary}}}/{F}_{{\rm{primary}}}$) is ten times larger than the expected flux ratio for main sequence components. The artificial brightening relative to a real binary system is to ensure that the temperature estimation uncertainties are separate from the overall sensitivity of the method, which we discuss in Section 5. We estimate the main sequence flux ratio from the published temperature of the late type star (given in Table 2) and the published spectral types of the primaries available on Simbad (Wenger et al. 2000), and convert to temperature and luminosity by using Table 5 of Pecaut & Mamajek (2013).

We analyzed each synthetic binary star system using the method described above, and measured the temperature (T_m) and variance (${\sigma }_{T}^{2}$) as a weighted sum near the grid point with the highest CCF peak value, weighting by the peak CCF height at each temperature (C_i):

$$\begin{eqnarray}&&{T}_{m}=\displaystyle \sum _{i}{C}_{i}{T}_{i}/\displaystyle \sum _{i}{C}_{i},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{\sigma }_{T}^{2}=\displaystyle \frac{{\displaystyle \sum }_{i}{C}_{i}{({T}_{i}-{T}_{m})}^{2}}{{\displaystyle \sum }_{i}{C}_{i}-{\displaystyle \sum }_{i}{C}_{i}^{2}/{\displaystyle \sum }_{i}{C}_{i}}.\end{eqnarray} \tag{ 4 }$$

Each synthetic binary observation contributes a pair of measured and actual (literature) companion temperatures, and so each late-type star in Table 2 has many independent temperature measurements made with the DSD method. To determine the correspondence between measured and actual temperature, we perform a Markov Chain Monte Carlo (MCMC) fit to a straight line using the emcee code (Foreman-Mackey et al. 2013). We plot the mean and standard deviation of the measured temperatures in Figure 1, along with 300 MCMC samples for the linear fit and the line of unity slope. The MCMC samples give posterior probability distributions for the parameters a and b relating the actual temperature (T_a) to the measured temperature (T_m) through

$$\begin{eqnarray}&&{T}_{m}=a+{{bT}}_{a}.\end{eqnarray} \tag{ 5 }$$

In Section 6 we invert this relation to determine the actual companion temperature, given the measured temperature from Equation (3).

Zoom In Zoom Out Reset image size

The typical temperature uncertainty with the DSD method is ∼150–200 K, but the optical instruments systematically overestimate the companion temperature. The error analysis is therefore not just important to measure the parameter uncertainties, but also to get the correct answer. We suspect the systematic biases come from a mismatch between the Phoenix model spectrum template and the real spectrum of a late-type star. The biases are different for each instrument because the instruments have different wavelength ranges, and so the spectral lines that contribute most to the CCF are different.

The detectability of a companion decreases primarily as the contrast between it and the primary star increases. Rotation plays an important role in the detection rate as well, since the CCF derives most of its power from narrow spectral features. We follow a similar strategy as above to estimate the detection rate as a function of temperature and rotational velocity for each star, with the key differences that we scale the model spectra to replicate a binary star observation with main sequence observations (rather than scaling the companion to ten times main sequence), and that we add Phoenix model spectra for late-type stars to the data instead of real spectra. We use synthetic spectra so that we can use a finer grid of temperatures and rotational broadening and not be limited by the temperatures or the temperature estimation uncertainties of real late-type stars. However, since we are comparing models to models any mismatch between the model spectrum and the real spectrum of a star of that temperature will tend to make the sensitivity calculations somewhat optimistic. This will have the largest impact for very cool stars, where the difficult to model molecular absorption is more important.

For each observed early-type star, we generate several synthetic binary star observations by adding model spectra for stars with ${T}_{\mathrm{eff}}=3000-7000\;{\rm{K}}$ in steps of 100 K and rotational velocities $v\mathrm{sin}i=0-50$ $\mathrm{km}\;{{\rm{s}}}^{-1}$ in steps of $10\;\mathrm{km}\;{{\rm{s}}}^{-1}$. For each temperature and $v\mathrm{sin}i$ combination, we make 17 independent synthetic observations by adding the model to the data with a radial velocity shift between −400 and 400 $\mathrm{km}\;{{\rm{s}}}^{-1}$ in steps of 50 $\mathrm{km}\;{{\rm{s}}}^{-1}$. We label a companion as detected if the highest peak in the CCF of the synthetic data with the model spectrum of the same temperature is within $5\;\mathrm{km}\;{{\rm{s}}}^{-1}$ (the approximate instrumental broadening) of the correct velocity.

The median detection rate for targets in Table 3 is shown in Figure 2.⁴ We can usually detect very cool stars if they are slowly rotating, but the sensitivity quickly degrades as the companion $v\mathrm{sin}i$ increases. Cool stars spin down as they age (Barnes 2003) so the rotation speed dependence is equivalent to an age dependence. We estimate the impact of rotation on our detection method by using the gyrochronology relation given in Barnes (2010):

$$\begin{eqnarray}&&\displaystyle \frac{{k}_{C}t}{\tau }=\mathrm{ln}(\displaystyle \frac{P}{{P}_{0}})\displaystyle \frac{{k}_{I}{k}_{C}}{2{\tau }^{2}}({P}^{2}-{P}_{0}^{2}).\end{eqnarray} \tag{ 6 }$$

Zoom In Zoom Out Reset image size

In Equation (6), k_C and k_I are constants fit to data with known ages and rotation periods, P and P₀ are respectively the current and zero-age main sequence (ZAMS) rotation periods, τ is the convective turnover time scale and t is the current age of the star. We use the same values that Barnes (2010) use for the constants:

• 1.  
  k_C = 0.646 day/Myr
• 2.  
  ${k}_{I}=452\;{\rm{Myr}}$/day.

We use Equation (6) to estimate the expected rotation period for a companion star of given temperature and age as follows: First, we convert from temperature to convective timescale (τ) by interpolating Table 1 in Barnes & Kim (2010). Next we sample an appropriate probability density function (PDF) for the age of the binary system; if the primary star was analyzed in David & Hillenbrand (2015), we use their posterior age PDFs. Otherwise, we use a uniform PDF from the Zero Age Main Sequence (ZAMS) age of the primary star to its main sequence lifetime (typically 10–200 Myr for our sample). Following the discussion in Barnes (2010), we uniformly sample initial rotation periods from 0.2 to 5 days for all stars. We estimate the current rotation period for each pair of age and initial rotation period samples using Equation (6) to build up a PDF of current rotation periods. We transform the period distribution into a PDF for $v\mathrm{sin}i$ using the main sequence radius of a star of the given temperature, obtained by interpolating Table 1 of Barnes & Kim (2010), and a uniform sampling of inclinations ($\mathrm{sin}i$). Figure 3 shows a typical $v\mathrm{sin}i$ distribution, which peaks near ∼5–10 $\mathrm{km}\ {{\rm{s}}}^{-1}$ and has a long tail extending to $\sim 40-50\ \mathrm{km}\;{{\rm{s}}}^{-1}$.

Zoom In Zoom Out Reset image size

By combining the sensitivity calculations described above with the $v\mathrm{sin}i$ samples, we marginalize over the expected rotation periods of the secondary stars to get simpler curves of detection rate as a function of companion star temperature. We show the median and approximate range of the marginalized detection rate in Figure 4. The DSD method can reliably detect companions as cool as 3700 K in most cases, although the primary star spectral type plays a dominant role in setting the coolest detectable companion.

Zoom In Zoom Out Reset image size

Companions with $T\gtrsim 6250\;{\rm{K}}$, the canonical limit at which the convective zone is too small to transfer angular momentum to the stellar wind and spin down the star (Pinsonneault et al. 2001), may have rotational velocities comparable to that of the primary star. In that case, estimating the primary star spectrum with a Gaussian filter may remove much or all of the companion spectrum. Since these are the types of stars with less extreme flux- and mass ratios, they are easier to detect with more conventional methods. However, this shortcoming could be overcome by using model spectra for the primary star as in Kolbl et al. (2015). In this work, we have optimized the method for finding cool companions.

We now use the DSD method to measure the temperatures of several known binary systems (Table 3). We cross-correlate the spectra against the full grid of model spectra enumerated in Section 2, and find the temperature of the companion using Equation (3). We then convert the measured temperature to PDFs of the true companion temperature using the MCMC chains developed in Section 4 (see also Figure 1). For stars with multiple observations, we multiply the PDFs from each detection. Finally, we calculate the companion temperature and confidence interval from the integral of the PDF:

$$\begin{eqnarray}&&f={\displaystyle \int }_{-\infty }^{x}\mathrm{pdf}(T){dT}.\end{eqnarray} \tag{ 7 }$$

We use as the central value the value of x such that f = 0.5 (the median). Likewise, we calculate the $1\sigma$ lower and upper bounds such that f = 0.16 and f = 0.84, respectively. The CCFs for the companions that we detect are shown in Figures 5 and 6. For each star, we show the CCF which has the maximum peak value and annotate the figures with the parameters. Most of the CCFs have very strong peaks. The exception is HIP 22958; however, the detection is strengthened by the fact that we observed this star twice and measured a similar temperature both times. The CCFs for HIP 22958 and HIP 24902 demonstrate the adverse effect a large companion rotational velocity has on the detection significance.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

6.1. Comparison to Literature Data

We use the literature data to predict an expected temperature for each companion in order to directly compare our measurements to previous results. The procedure outlined in Section 3 using the magnitude difference or orbital information alone produces reasonable estimates, but in many cases there is additional information in the literature to refine the estimates. The refined estimates are described below.

HIP 76267 and HIP 84606 are found in the David & Hillenbrand (2015) sample; we use the mass and temperature estimates provided there rather than going through the Simbad spectral type and assuming main-sequence relationships.

Shatsky & Tokovinin (2002) provide a color estimate of the companion star to HIP 79199 ($J-K=0.57\pm 0.12$). We convert this directly into a temperature estimate through Table 5 of Pecaut & Mamajek (2013).

Zorec & Royer (2012) find fundamental parameters for HIP 22958, and determine a temperature slightly cooler and luminosity much greater than the spectral type (B6V) would suggest. Because of this the usual analysis, which uses main sequence relationships, results in a biased answer. We estimate the companion temperature by assuming that the companion does follow the main sequence relationships as described in Section 3, but sample the uncertainty distributions given in Zorec & Royer (2012) for the temperature and radius of the primary star.

We compare our companion temperature measurements from the DSD method to the estimates described above in Figure 7. There is overall excellent agreement between the temperatures, with 5/6 falling within $1\sigma$ of equality. We test for a bias (Δ) between the measured temperatures (T_m) and expected temperatures (T_a) with the equations

$$\begin{eqnarray}&&{\rm{\Delta }}=\displaystyle \sum _{i}({T}_{m,i}-{T}_{a,i})\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{\sigma }_{{\rm{\Delta }}}^{2}=\displaystyle \sum _{i}({\sigma }_{{T}_{m,i}}^{2}+{\sigma }_{{T}_{a,i}}^{2})\end{eqnarray} \tag{ 9 }$$

which results in ${\rm{\Delta }}=-580\pm 770\;{\rm{K}}$. Our temperature measurements are consistent with the expected temperatures.

Zoom In Zoom Out Reset image size

We list our measurements as well as the expected temperatures described above in Table 6. The expected $v\mathrm{sin}i$ values come from application of Equation (6) as described in Section 5. While we do give the measured $v\mathrm{sin}i$ and metallicity for our detections, the accuracy of these parameters is not calibrated and is determined with a coarse grid; the values should only be taken as rough estimates. We do note that most of the measurements have $[\mathrm{Fe}/{\rm{H}}]=-0.5$. This is likely a measurement bias since we do not expect the binary systems to have significantly sub-solar metallicity. As metallicity increases, so do the line depths of most of the lines in the spectrum. Any lines that are poorly modeled will then have a larger negative impact on the resulting CCF; thus the bias toward low metallicity is likely a result of imperfect model atmosphere templates. We do not attempt to identify the poorly modeled lines in this work.

Table 6.  Companion Data

	Measured Values			Expected Values
Primary	${T}_{\mathrm{eff}}$	[Fe/H]	$v\mathrm{sin}i$	${T}_{\mathrm{eff}}$	$v\mathrm{sin}i$
Star	(K)	(dex)	(km s−1)	(K)	(km s−1)
HIP 13165	5770 ± 162	−0.5	5	⋯	⋯
HIP 22958	6070 ± 112	−0.5	30	${6240}_{-409}^{+579}$	${11}_{-8}^{+16}$
HIP 24902	5680 ± 154	0.0	30	${5950}_{-84}^{+34}$	${4}_{-3}^{+3}$
HIP 33372	⋯	⋯	⋯	${6955}_{-550}^{+238}$	${12}_{-8}^{+16}$
HIP 65477	⋯	⋯	⋯	${3861}_{-11}^{+34}$	${2}_{-1}^{+3}$
HIP 76267	5450 ± 158	−0.5	5	${5670}_{-281}^{+193}$	${3}_{-2}^{+4}$
HIP 77516	6820 ± 141	0.5	5	${6600}_{-167}^{+530}$	${12}_{-8}^{+17}$
HIP 79199	4620 ± 158	−0.5	5	${4800}_{-484}^{+484}$	${6}_{-4}^{+7}$
HIP 79404	4770 ± 112	−0.5	10	⋯	⋯
HIP 81641	⋯	⋯	⋯	${6475}_{-96}^{+126}$	${10}_{-7}^{+12}$
HIP 84606	5480 ± 154	0.0	10	${5450}_{-45}^{+71}$	${3}_{-2}^{+3}$
HIP 88290	⋯	⋯	⋯	${5847}_{-24}^{+41}$	${4}_{-3}^{+4}$
HIP 91118	6490 ± 154	−0.5	10	⋯	⋯
HIP 92027	⋯	⋯	⋯	${6752}_{-96}^{+315}$	${12}_{-8}^{+16}$
HIP 98055	⋯	⋯	⋯	${7366}_{-113}^{+391}$	${13}_{-9}^{+18}$
HIP 113788	⋯	⋯	⋯	${6276}_{-92}^{+34}$	${7}_{-4}^{+6}$
HIP 116247	⋯	⋯	⋯	${6351}_{-218}^{+482}$	${8}_{-6}^{+10}$
HIP 116611	⋯	⋯	⋯	${3842}_{-50}^{+103}$	${2}_{-1}^{+4}$

Download table as:  ASCIITypeset image

6.2. Non-Detections

There are many companions in Table 3 that we do not detect. Most of these are single-lined spectroscopic binaries (Table 4), and are likely too cool to detect with our data; very high signal-to-noise spectra with a near-infrared instrument such as IGRINS may uncover them. Several of the remaining un-detected companions have expected temperatures $T\gt 6250\;{\rm{K}}$, and so are likely to be rapid rotators. Since the CCF gets most of its power from sharp spectral features, these rapidly rotating companions are difficult to detect (see Figure 2).

Finally, HIP 88290 is hot enough and expected to be rotating slowly enough that we should be able to easily detect it. In fact, we would expect to be able to directly see the companion in the spectra (the green lines in Figure 8). The fact that we do not see the composite spectrum or see a peak in the CCF implies that the companion must be rotating with $v\mathrm{sin}i\gt 50\;\mathrm{km}\;{{\rm{s}}}^{-1}$, much more quickly than Equation (6) predicts, that the primary is a giant and therefore much brighter than main-sequence relationships suggest, or that the companion fell outside the spectrograph slit. This star is in the David & Hillenbrand (2015) sample and has an effective temperature and mass consistent with main sequence, so we can rule out the giant primary possibility. Additionally, the binary separation is 047 (Tokovinin et al. 2015) and the CHIRON spectrograph has a ∼27 diameter fiber; light from the companion is guaranteed to fall on the slit.

Zoom In Zoom Out Reset image size

We have presented and extensively characterized the direct spectral detection method for finding companions to intermediate-mass stars using high-resolution cross-dispersed échelle spectroscopy. Using a very large number of synthetic but realistic binary star observations, we constrained the uncertainty and systematic errors present in determining the companion temperature with the direct spectral detection method. The typical uncertainties are of the order of 200 K across all instruments used in this study, with a systematic offset of similar magnitude (for the optical instruments). We used the synthetic binary star analysis to calibrate the direct spectral detection method for the four instruments used in this study between the temperatures $3500\;{\rm{K}}\lt T\lt 6500\;{\rm{K}}$.

We also estimated the sensitivity to detection of companions with a range of temperature and $v\mathrm{sin}i$ by creating a second set of synthetic companions. The method can detect companions as late as M0 in most cases, although the lower limit depends on the primary star spectral type, the signal-to-noise ratio achieved, and the instrument used. The median detection limit corresponds to average flux ratios as small as ${F}_{{\rm{sec}}}/{F}_{{\rm{prim}}}\sim {10}^{-3}$ and binary mass ratios ${M}_{{\rm{sec}}}/{M}_{{\rm{prim}}}\sim 0.2$, or a main-sequence M0 star orbiting an A0V primary.

The lowest detectable mass ratio is even more striking for young stars. At 1 Myr, both the A0 star and its companion are still contracting onto the main sequence (Bressan et al. 2012). The flux ratio limit corresponds to a ∼M1 companion, similar to the main sequence case. However the mass ratio in this young system is ${M}_{{\rm{}}\;{\rm{sec}}}/{M}_{{\rm{prim}}}\sim 0.1$, half that of main-sequence components with similar spectral types. The direct spectral detection method is therefore well suited for finding close, low-mass companions to massive young stars.

There is also an upper detection limit near $6500\;{\rm{K}}$ set by rotation. Our method of removing the primary star spectrum can also remove the companion spectrum if it has a similar rotational velocity, which hot companions are likely to have. Subtracting a model atmosphere for the primary star would remove the upper limit, but would reduce the detection rate for cool companions that are most difficult to detect with any other means.

Finally, we applied the direct spectral detection method to a set of known binary systems with close, late-type companions. We detected the companion spectrum in 9 of 34 known binary systems, 3 of which we characterized for the first time. Most of the companions we failed to detect are likely very cool, falling below the sensitivity limit of our data.

The direct spectral detection method is able to detect close binary companions with comparable or better sensitivity than imaging techniques, and does not require large telescopes with extremely competitive time allocation requests. This method is an excellent way to identify and perform initial characterization on new binary systems using smaller telescopes, but care must be taken to calibrate the parameter estimation.

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 (Astropy Collaboration, 2013). It was supported by a start-up grant to Adam Kraus as well as a University of Texas Continuing Fellowship to Kevin Gullikson. J.-E. Lee was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) (grant No. NRF-2015R1A2A2A01004769) and the Korea Astronomy and Space Science Institute under the R&D program (Project No. 2015-1-320-18) supervised by the Ministry of Science, ICT and Future Planning.

This work used the Immersion Grating Infrared Spectrograph (IGRINS) that was developed under a collaboration between the University of Texas at Austin and the Korea Astronomy and Space Science Institute (KASI) with the financial support of the US National Science Foundation under grant AST-1229522, of the University of Texas at Austin, and of the Korean GMT Project of KASI.

The Hobby-Eberly Telescope (HET) is a joint project of the University of Texas at Austin, the Pennsylvania State University, Stanford University, Ludwig-Maximilians-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.

Based on observations at Cerro Tololo Inter-American Observatory, National Optical Astronomy Observatory (NOAO Prop. IDs: 13A-0139, 13B-0112, 2014A-0260, 14A-0260, 15A-0245; PI: Kevin Gullikson), which is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.

We would like to thank Bill Cochran and Mike Endl for observing some of the spectra used in this project. Finally, we would like to thank the anonymous referee for numerous comments that improved the paper.