DETECTION OF LOW-MASS-RATIO STELLAR BINARY SYSTEMS

O- and B-type stars are often found in binary or multiple systems: Mason et al. (2009) estimate that at least 57% of O stars are in spectroscopic multiple systems, and at least 75% are in any type of binary or multiple system. Yet the multiplicity fraction of high-mass stars may be underestimated due to the difficulty of detecting low-mass secondary stars (Sana & Evans 2011). While the mass-ratio distribution is reasonably well known for high-mass binaries with mass ratio q ≡ M_s/M_p > 0.2, there are almost no constraints for low-mass-ratio binaries. However, binaries of low mass ratio are important probes of star formation since they may have formed in a different way than approximately equal-mass binaries. Here we define low mass ratios to be those with q < 0.2, where M_s is the mass of the secondary (lower-mass binary component), and M_p is the mass of the primary (higher-mass binary component). We also use the term "low mass" to describe star with M ≲ 1 M_☉, and "high mass" to describe stars with M ≳ 5 M_☉.

1.1. Binary Formation in High-mass Stars

1.2. Observing Low Mass-ratio Binaries

We describe here our method for detecting the emission from an approximately solar-mass star orbiting an early B-type star, which we will hereafter call the direct spectral detection method. The basis of this method is to cross-correlate a high signal-to-noise ratio (S/N) B-star spectrum with a synthetic F, G, or K star spectrum. If a low-mass star with such a spectrum is orbiting the B star, we expect to find a peak in the cross-correlation function at the radial velocity corresponding to the low-mass star's motion. A peak in the cross-correlation function should appear even if the flux from the low-mass star is comparable to or even slightly less than the noise level in the spectrum. Figure 1 illustrates the approximately limiting case for the flux ratio. The top panel shows a fully reduced CRIRES spectrum of HIP 108975 (see Section 4) with a model spectrum for a 0.9 M_☉ star at a realistic flux ratio below it. We used evolutionary tracks published by Landin et al. (2008) to evolve the secondary star to 50.1 Myr, the age of the system (Tetzlaff et al. 2010), in order to determine the flux ratio between the primary and secondary. The secondary star model was then added to the telluric-corrected B-star spectrum at two different velocities. It is clear that the model spectrum has an amplitude much smaller than the noise. Nonetheless, the bottom two panels show that a cross-correlation will have a peak with high significance at the velocity of the secondary star.

Zoom In Zoom Out Reset image size

A careful choice of the wavelength region is critical for the direct spectral detection method. First, we want a wavelength region where the B-star spectrum is mostly continuum (i.e., very few spectral lines). Since B stars have few spectral lines, it is easy to find such a spectral region. Second, we want a region where the low-mass star would have many closely spaced, strong lines. The more lines there are in the low-mass star, the stronger the peak will be in the cross-correlation function. Finally, we want a spectral region where the flux ratio between the low-mass and the high-mass star is maximized. It is not helpful to go much redder than a few microns for companions with T > 4000 K, because both the high-mass and the low-mass star are firmly in the Rayleigh–Jeans limit by this point, where the flux ratio is approximately constant. For this project, we choose wavelengths from 2300 to 2400 nm, which is the CO Δν = 0–2 bandhead in the low-mass star.

There is both a lower and upper mass detection limit. Secondary stars that are too cool will be too faint, and any signal will be lost in the noise. Additionally, a more massive (and hotter) primary star will decrease the flux ratio, and push the lower-mass limit up. On the other end, secondary stars that are too hot will dissociate CO, destroying the bandhead that we are looking for. The temperature and size of the secondary star will depend on its age as well as its mass since it will still be evolving toward the main sequence during the lifetime of the B star. The exact mass sensitivity will thus depend on the age, primary star mass, and S/N of the system being observed.

The detector resolution is also important for the direct spectral detection method. Deeper lines, providing more contrast from the continuum, are easier to detect than broad, shallow lines. In addition, narrow spectral lines will result in a stronger, narrower peak in the cross-correlation function, which is most sensitive to the steep line edges. Therefore, we want the spectral lines in the low-mass companion to be as deep and narrow as possible. The intrinsic width of CO bandhead lines is roughly 5–7 km s⁻¹. In order for the observed line width to be this small, we need the resolution of the instrument to be R = λ/Δλ ≳ 50, 000.

There are two main difficulties with the direct spectral detection method: telluric line removal and the low flux ratio between the primary and secondary star. Figure 2 shows the transmittance through Earth's atmosphere (the telluric spectrum) in the wavelength range we are interested in. Most of the spectral lines are from methane, with a few deep water lines toward the red end of the range shown (see Section 4 for details on the telluric line removal). The low flux ratio makes the telluric contamination especially troublesome, since the telluric lines are stronger than the lines in the companion star spectrum. The flux ratio of F_s/F_p ∼ 10⁻² effectively sets a lower limit on the S/N for which the direct spectral detection method is possible. Any flux coming from a low-mass star will be completely buried in the Poisson noise for spectra with S/N ≪100. Removal of the telluric contamination will add more noise to the spectrum, so a spectrum should have S/N of a few hundred before telluric line removal to have a good chance of detecting a companion.

Zoom In Zoom Out Reset image size

B-type stars are commonly used in the near-IR as telluric standard stars. Astronomers will observe their science targets, and then move to a B-type star. Since B-type stars have few spectral lines relative to cooler stars, most of the observed spectral lines will be from the absorption of Earth's atmosphere (telluric absorption). Therefore, these stars provide an empirical estimate of the telluric spectrum; division of the science spectrum by the normalized standard star spectrum will mostly remove the telluric lines.

Since B-type stars are commonly used as above, there are many high signal-to-noise ratio (S/N ≳ 100) observations of such stars in archived data. We used the VLT/CRIRES archive in this project. CRIRES is a high-resolution (R = λ/Δλ ≈ 100, 000) infrared spectrograph on the VLT at Paranal Observatory. The detector consists of four 1024×512 CCD chips that are mosaicked end to end, and the spectrum falls across them. There are several wavelength settings available, which determine what parts of the spectrum fall on each chip. For wavelength settings in the CO bandhead near 2300 nm, each chip will hold roughly 10 nm of spectrum with roughly 1–2 nm gaps between the chips.

To generate the sample, we started with all single, main-sequence B0–B5 stars with CRIRES observations from 2300 to 2400 nm. We then excluded any shell stars, which have circumstellar disks (Porter & Rivinius 2003) that may create false positives. Table 1 shows the complete sample used in this project. The spectral types and ages were obtained from a catalog of nearby young stars (Tetzlaff et al. 2010). The one exception is HIP 97611, which was not in Tetzlaff et al. (2010). For this star, the age was taken from Westin (1985) and the spectral type from the Simbad database.¹ The distances to all stars were determined from parallaxes given in the Simbad database. The maximum separation column estimates the approximate maximum separation of the binary orbit we are sensitive to, assuming a seeing of 08 which is typical of Paranal Observatory. The median of the maximum separations to which we are sensitive is 124 AU. The final column gives the number of distinct observations of the star. We count all nodding positions taken on a given night with the same detector wavelength setting as one observation.

The data reduction was done using standard methods in IRAF.² All observations were taken in an AB or ABBA nodding pattern. For each set of AB nods, A–B and B–A frames were made to remove any atmospheric emission lines and dark current. The resulting difference images were then treated to a quadratic nonlinearity correction, using coefficients made available by the CRIRES team.³ The corrected frames were then divided by a normalized flat field. Due to the slit curvature, the spectrum can shift by up to a pixel in the dispersion direction between the A and B nod positions. Therefore, combining the two-dimensional frames before extraction can reduce the spectral resolution and affect the line shapes. For this reason, we combined the nodding positions only after the wavelength calibration and telluric correction. Each nod position was extracted using the optimal algorithm in the apall task in IRAF. The spectra were wavelength calibrated using a model telluric spectrum generated with the atmospheric modeling code LBLRTM⁴ (Clough et al. 2005).

For telluric correction, we used a similar procedure to the one outlined by Seifahrt et al. (2010). The atmosphere modeling code LBLRTM was used to generate a synthetic telluric absorption spectrum. The abundances of water, methane, and carbon monoxide were fit using a Python implementation of a Levenberg–Marquardt nonlinear least-squares fitting algorithm. The Levenberg–Marquardt fit also refined the wavelength solution to the telluric model, fit the continuum, and fit the resolution of the spectrograph with a Gaussian profile. The FWHM of the profile was the only free parameter in the resolution fit.

The LBLRTM code expects a model atmosphere, which contains the temperature, pressure, and abundance of 30 molecules as a function of atmospheric height. For the majority of molecular species, we used a mid-latitude nighttime MIPAS⁵ profile, which provides the temperature, pressure, and abundances of various molecules in 1 km intervals from sea level to 120 km. The low-altitude (z ≲ 30 km) temperature, pressure (Kerber et al. 2010), and humidity (Chacón et al. 2010) profiles were obtained from radiosonde data taken from Paranal Observatory.

The LBLRTM atmospheric modeling code comes with a molecular line list based on the HITRAN 2008 database (Rothman et al. 2009), with a few molecules individually updated. Since none of these updates were relevant for the wavelength range from 2300 to 2400 nm, we in essence used the stock HITRAN 2008 database. However, in the process of modeling, we found several water and methane lines that were consistently underfit or overfit. For these cases, we manually adjusted the line strengths in the database. The line strengths were fit visually and should not be considered rigorous new line strengths. Table 2 summarizes these changes.

In some of the 2007 data, the first chip was not well illuminated by the flat-field lamp. This introduced an unphysical continuum shape in the data and made the resulting model fit very poor. For these cases, we ignored the first chip in further analysis. In addition, the fourth chip has several bad pixels on the left edge and a streak down the middle. None of the telluric model fits were very good on this chip, and so we have ignored it completely in our analysis.

After the observed spectrum was fit, we found that the residuals still contained large spikes, even on the good detector chips. These spikes can come from a variety of sources. For the deepest lines, simple Poisson noise can create large residuals when dividing by the telluric model. In addition, a poorly fit continuum may cause the model to overestimate or underestimate the abundance of a given molecule. This can be especially troublesome for water lines, for which only a few exist in the wavelength region we are investigating. If a strong water line is near the edge of the chip, where the continuum is usually least certain, the best-fit water abundance may be skewed and cause none of the water lines to be well fit. We do know that large residuals are not coming from the spectrum of a low-mass star, due to the expected flux ratio between the primary and secondary, F_s/F_p ∼ 10⁻². Any residuals with amplitude greater than 1% of the continuum level come from uncorrected telluric lines, cosmic rays, or bad pixels.

In order to minimize these spikes, we performed a second fit to any residuals significantly above the continuum noise level. To make sure we were not fitting away any low-mass star lines, we only corrected spikes whose amplitude was greater than 5% of the continuum level. In this second fit, we first attempted to fit a Gaussian to each spike. If the spike was well fit by a Gaussian, we divided the residuals by the fit. If not, we simply masked out the line core, so that it would not affect the cross-correlation in later analysis (see Section 2). Figure 3 shows the steps involved in the telluric correction. Note that the secondary correction removes the large residuals, while leaving the rest of the spectrum unaffected.

Zoom In Zoom Out Reset image size

The telluric correction described above usually reduced any telluric lines to near the Poisson noise level in the spectrum, which is the best a fitting routine can do. To search for any systematic errors in the telluric correction, we added all spectra of each wavelength setting together to make a series of master telluric residual spectra. These master spectra had less random noise than any individual observation, and therefore we were immediately able to see whether some telluric lines are systematically underfit or overfit. We found that for wavelength settings with at least 10 spectra in our sample, dividing the corrected spectra by this systematic telluric residual template increased the sensitivity to companions. We did not make this final correction for wavelength settings with fewer than 10 individual spectra in our sample.

Each telluric residual spectrum was cross-correlated against a suite of model atmospheres generated by the Phoenix stellar atmosphere code (Hauschildt et al. 1999). All model spectra had solar metallicity. The effective temperatures ranged from 3000 to 7200 K, in 100 K intervals. We used several surface gravities based on the stellar temperature. For the model secondary stars with 3000 < T_eff < 3600, which would have to be very young (and large) to be detectable, we used a log (g) = 3.5. For the secondaries with 3600 < T_eff < 6500, we used log (g) = 4.0. Finally, we used log (g) = 4.5 for T_eff > 6500, which can be detected closer to the main sequence. We found that the surface gravity has only a very small effect on the cross-correlation, which is more sensitive to the line position than its precise width or depth. We compiled a list of all cross-correlations that show a single peak with at least 3σ significance. For a given telluric-corrected residual spectrum, several model atmospheres may generate a significant peak at the same velocity. This is because the model spectra of two stars differing by only a few hundred kelvin are not very different. To keep from counting peaks twice, we only counted the cross-correlation that resulted in the most significant peak at a given velocity. We then attempted to reject spurious peaks caused by the noise or incomplete telluric line removal in a multi-stage process.

The first rejection stage was done by identifying peaks in the cross-correlation caused by telluric residuals. To do this, we cross-correlated a spectrum uncorrected for telluric absorption with the same suite of model atmospheres as we used for the corrected spectra (see above). We did these cross-correlations for one observation of each wavelength setting. The cross-correlation of uncorrected spectra with model-atmosphere spectra generated a series of cross-correlations with peaks arising exclusively from telluric lines. We visually compared all of the binary candidate signals with these telluric cross-correlations. If the dominant cross-correlation peak was at the same velocity and had a similar width as a peak in the telluric cross-correlation function corresponding to the same wavelength setting and secondary model temperature, we assumed that the peak was caused by incomplete telluric removal and rejected the candidate. There were several cross-correlations with peaks at the same location as a telluric peak, but with a different width. In these cases, we marked the candidate as probably coming from incomplete telluric correction, but did not reject the candidate.

Next, we determined whether the S/N and telluric line removal in a given observation would allow us to detect the candidate companion star. To do this, we added a model atmosphere with the same temperature as the candidate to the telluric-corrected spectrum at 17 different radial velocities ranging from −400 to 400 km s⁻¹. We do not expect to see any peaks from real companions with |v| > 400 km s⁻¹, the approximate radial velocity of a 1 M_☉ star orbiting a 10 M_☉ star such that the stellar surfaces are in contact. The flux ratio of the model atmosphere to the primary was obtained by interpolating pre-main-sequence evolutionary tracks from Landin et al. (2008) at the age of the system, as well at age±σ_age. The primary star ages for our sample are given in Table 1. We cross-correlated each of these semi-synthetic spectra against the model spectrum; if the largest peak was at the correct velocity, we counted the star as detected. If the star was not detected in the sensitivity analysis at least 50% of the time, we rejected the candidate. The significance of the correct peak in the cross-correlation function can vary greatly, depending on where the stellar spectrum falls in relation to the telluric line residuals. Therefore, we cannot usually reject a peak based solely on its significance.

We then visually inspected the remaining candidate cross-correlations, picking out those with a single dominant peak with |v_r| < v_max where v_max is the maximum possible radial velocity for a star of temperature T_sec to be orbiting a hotter star of temperature T_prim with semimajor axis a. Assuming a circular orbit, v_max is given by

$$\begin{equation} v_{\rm max} = \sqrt{\frac{2G(M_{\rm prim} + M_{\rm sec})}{R_{\rm prim}}} \cdot \frac{T_{\rm sec}}{T_{\rm prim}}, \end{equation} \tag{ 1 }$$

where M_prim and M_sec are the masses of the primary and secondary stars, respectively, and R_prim is the radius of the primary star. The primary masses are given in Tetzlaff et al. (2010), while the radii and primary star temperatures were estimated from spectral type relations given in Carroll & Ostlie (2006). An eccentric orbit could have a larger maximum velocity than that estimated by Equation (1), if the orbit was oriented such that the secondary star was moving directly toward or away from earth at or near periastron. A star in an eccentric orbit cannot get so close to the primary that they touch, and its average distance must still be far enough to allow for the observed secondary star temperature. These conditions lead to a maximum eccentricity, given by

$$\begin{equation} e_{\rm max} = 1 - \frac{v_{\rm max}^2 (R_{\rm prim} + R_{\rm sec})}{G(M_{\rm prim} + M_{\rm sec})}, \end{equation} \tag{ 2 }$$

where R_prim and R_sec are the radii of the primary and secondary stars, respectively, and v_max is the maximum circular velocity given by Equation (1). Typical values give e_max ≈ 0.6. Eccentricities near this value could allow for velocities significantly greater than v_max given by Equation (1).

The above analysis is summarized in Table 3. There are two binary candidate systems that we have not been able to reject. For both of these, we checked what other observations the candidate star had within the CO bandhead spectral region (2300–2400 nm). The analysis of each of these stars is done separately below.

5.1. HIP 26713

The cross-correlation for this candidate is shown in Figure 4. The strong peak at −220 km s⁻¹ has a significance just over 4σ, and corresponds to a 5600 K star model. A sensitivity analysis (Table 3, step 2) gives a median peak significance of ∼6σ with a large (∼2σ) spread. Due to the large spread, we cannot reject the peak based on the observed significance.

Zoom In Zoom Out Reset image size

The candidate radial velocity amplitude of 220 km s⁻¹ is very near the upper limit given by Equation (1) of 256 km s⁻¹. If this candidate is a real binary companion, it must have been observed very near its radial velocity maximum and the orbital inclination must be very near edge-on. Assuming a circular orbit and equal probabilities of observing any given phase or inclination, the probability of this occurring for this system is ∼0.01. However, the probability for one system in our entire sample to be caught with this chance alignment rises to 0.14. In addition, we cannot discount the possibility of an eccentric orbit leading to a higher value of v_max than estimated by Equation (1). Without more data, we can neither confirm nor disprove the existence of a companion star orbiting HIP 26713.

If HIP 26713 is a true binary system, evolutionary tracks by Landin et al. (2008) give a secondary star mass of 1.6 ± 0.2 M_☉. The mass of the primary star is 9.4 ± 0.2 M_☉ (Tetzlaff et al. 2010), giving a mass ratio of q = 0.17 ± 0.02. For a circular orbit, the mass ratio and observed velocity correspond to an orbital period of 10.0 days. If the companion star is on an eccentric orbit, its true period would be longer than this.

5.2. HIP 92855

There were seven observations of HIP 92855 on different dates, all with the 2336 nm wavelength setting. The cross-correlations for four of the observation dates show a single strong peak when using a 6100 K model star as template. Figure 5 shows the cross-correlations of the telluric-corrected spectra with a 6100 K model for all of the observations. Table 4 summarizes the sensitivity and cross-correlation significance. The detection rate is the fraction of the 17 radial velocities between −400 and 400 km s⁻¹ that were correctly detected in the sensitivity analysis (Table 3, step 2). The variation in detection rate is due to the different S/N levels and telluric line corrections in the observations at different dates. The expected significance is the median significance of the radial velocities which were detected in the sensitivity analysis, in units of the standard deviation of the cross-correlation function. The observed significance and velocity are for the observed peaks. The velocities in Table 4 are corrected for the barycentric motion and the known systematic radial velocity of HIP 92855, while those in Figure 5 are not.

Zoom In Zoom Out Reset image size

Table 4. Summary of HIP 92855 Observations

Date	Detection Rate	Expected Significance	Observed Significance	Velocity
				(km s−1)
2007 May 9	0.35	3.3σ	4.1σ	−74
2007 Jun 9	0.18	3.6σ	N/A	N/A
2007 Jul 25	0.35	3.7σ	N/A	N/A
2007 Aug 2	0.47	3.6σ	3.5σ	−234
2007 Sep 16	0.82	4.2σ	N/A	N/A
2008 Sep 19	0.71	3.4σ	4.1σ	−126
2008 Sep 10	0.59	3.3σ	4.0σ	116

Notes. Significance is in units of the standard deviation of the cross-correlation function. The detection rate and expected significances are from the sensitivity analysis (Table 3, step 2).

Download table as:  ASCIITypeset image

As Table 4 shows, a 6100 K star orbiting HIP 92855 is at the limit of detectability with the direct spectral detection method. With the exception of the observation on 2007 September 16, the cross-correlations for the dates with the highest detection rates have a single large peak. We consider this an excellent candidate for follow-up observations.

If HIP 92855 is a true binary system, evolutionary tracks by Landin et al. (2008) give a secondary star mass of 1.2 ± 0.2 M_☉. The mass of the primary star is 7.8 ± 0.2 M_☉ (Tetzlaff et al. 2010), giving a mass ratio of q = 0.15 ± 0.04. The maximum radial velocity, observed on 2007 August 2, was 234 km s⁻¹. If this is the radial velocity semi-amplitude and assuming a circular orbit, the binary orbit would have a period of 6.8 days. If the true velocity semi-amplitude is larger either because no observation was taken when the companion star was at quadrature or because the orbit is inclined, the period would be shorter than this. The period may be longer than 6.8 days if the orbit is eccentric.

We now estimate the completeness of the direct spectral detection method applied to this data set. For each telluric-corrected observation, we created a series of synthetic binary-star spectra by adding stellar models to the data at various flux ratios, temperatures, and radial velocities. We used evolutionary tracks from Landin et al. (2008) to find the luminosity of model companion stars with temperatures ranging from 3000 K to 7000 K in steps of 500 K. To find the model flux ratio F_s/F_p, we used the best-fit age for each star quoted in Table 1, as well as the best-fit age±σ_age. Model secondary spectra generated with the Phoenix code (Hauschildt et al. 1999) were added to the telluric-corrected observations at 17 different radial velocities ranging from −400 to +400 km s⁻¹. Changing the radial velocity of the model spectrum changes where the companion spectral lines fall with respect to the telluric lines. Finally, we cross-correlated each synthetic spectrum with its corresponding model secondary star spectrum, and examined the cross-correlation function.

If the highest peak in the cross-correlation function was at the correct velocity, the companion was considered detected. We then tabulated how many times the companion star was detected in the 17 radial velocity trials. Figure 6 shows the fraction of trials that detected the companion for all of the model radial velocities, as a function of primary (B-star) mass and the binary mass ratio. The points correspond to individual spectra, with their sizes indicating the S/N in the spectrum, and the contours are drawn by interpolating between the points. The four panels are for different companion star temperatures. Companion stars with effective temperatures from 4600  to 5400 K have large regions with a very high detection rate. Stars cooler than about 4600 K are too dim to detect without much higher S/Ns than present in our data set, and stars hotter than about 5400 K do not have a strong CO bandhead and so the cross-correlation function is not as sensitive.

Zoom In Zoom Out Reset image size

Figure 6 shows that the direct spectral detection method is able to find companion stars with a mass ratio of q ≈ 0.1–0.2, for a range of effective temperatures. The regions with a detection rate near 1 are completely sampled, and a companion star in that region would be detected. For primary stars with M < 10 M_☉, we are sensitive to almost all companions with 4600 K < T_eff < 5400 K.

We have not found any unambiguous low-mass-ratio companions in our sample, though we do have two candidates that require follow-up observations (HIP 92855 and HIP 26713). From the work described in Section 6, we define a range of mass ratios for which the direct spectral detection method is sensitive for each primary (B) star. We can then rule out any companions with mass ratios in that range, for that primary star.

In order to convert these star-by-star limits on the presence of a companion into upper limits on the multiplicity fraction of the parent population, we first count the number of stars that rule out companions in a particular range of mass ratios. We then apply binomial statistics, where the probability P of finding k companions from n samples of a parent population with a true binary fraction p is given by

$$\begin{equation} P(k|p,n) = \frac{n!}{k!(n-k)!}p^k(1-p)^{n-k}. \end{equation} \tag{ 3 }$$

For no detected binary companions (k = 0), the corresponding likelihood function for the binary fraction is

$$\begin{equation} P(p|k=0, n) = (n+1)\cdot (1-p)^n. \end{equation} \tag{ 4 }$$

A 90% upper limit is given by

$$\begin{equation} 0.9 = \int _0^p (n+1)\cdot (1-p^{\prime })^n dp^{\prime } \end{equation} \tag{ 5 }$$

with the solution

$$\begin{equation} p_{90} = 1 - 0.1^{\frac{1}{n+1}}. \end{equation} \tag{ 6 }$$

A similar derivation gives 90% upper limits for one detection (k = 1). Figure 7 shows the 90% upper limits to the binary fraction as a function of mass ratio. To find n in Equation (6), we counted the number of stars that ruled out a companion at the given mass ratio. We only counted companions that were found in all 17 radial velocity trials (see Section 6), and were always found with at least 4σ significance. We also included upper limits assuming that HIP 92855, the more likely of our two candidates, is a true binary system. Figure 7 also shows the lower multiplicity limit set by Evans et al. (2011; blue dotted line) and the intrinsic O-star mass-ratio distribution derived by Sana et al. (2012). Evans et al. (2011) do not split their multiplicity by mass ratio, and so the line shown in Figure 7 is an average value. While we include them for comparison with our results, neither of these studies are directly comparable to our sample. Sana et al. (2012) have only O stars in their sample and are only able to measure mass ratios q ≈ 0.2–1, although they consider mass ratios down to q = 0.1 when deriving the intrinsic distribution. Evans et al. (2011) are sensitive to similar mass ratios as our sample, but they sample late B stars (B4–B9) while we sample early B stars (B0–B5). Our results are almost perfectly complementary to those of Evans et al. (2011). It is encouraging that our upper limits for 0.1 < q < 0.2, where our sample is most complete, lie in between the results of Sana et al. (2012) who measure the binary fraction of more massive primaries, and Evans et al. (2011) who measure the binary fraction of less massive primaries than we do.

Zoom In Zoom Out Reset image size

The above analysis assumes that the sample of B stars given in Table 1 is representative of the B-star population as a whole. Most of the sample stars are field B stars, which have a lower overall multiplicity than cluster or association B stars (Mason et al. 2009). However, the close binaries this method is sensitive to would be difficult to disrupt with dynamical interactions in a cluster environment, and we therefore expect this sample to be representative of both populations. A potential complication is that telluric standard stars are often chosen specifically because they are not known binary systems. While this may introduce a bias for large mass ratios, the direct spectral detection method is sensitive to low mass-ratio companions that may not have been found with other methods such as classical spectroscopy or imaging and so companions with q ≲ 0.2 may be less effected. The amount of bias introduced into the above measurement depends on how carefully the individual observers chose the telluric standard stars, and so is very difficult to assess.

We have described a new technique for finding binary systems with a flux ratio of F_p/F_s ≈ 100, where F_p and F_s are the fluxes from the primary and secondary star, respectively. In this technique, which we call the direct spectral detection technique, we use high S/N, high-resolution spectra of a binary candidate. We remove the contamination from Earth's atmosphere with the telluric modeling code LBLRTM, and cross-correlate the residuals with a library of stellar models for late type stars (F2–M5). A binary detection would appear as a strong peak in the cross-correlation function.

We prove the feasibility of the direct spectral detection method by adding a synthetic signal to real data, and successfully recovering it. We further investigate the completeness of the method in Section 6. This method is sensitive to detecting a range of companion stars, set by the spectral type and age of the primary and the S/N of the observation. Our sample is sensitive to almost all companion stars with 4600 < T < 5400, corresponding to binary mass ratios of 0.1 ≲ q ≲ 0.2.

We have applied this technique to a sample of 34 archived main sequence early B stars (B0–B5) with spectra taken with the CRIRES near-infrared spectrograph. We found no unambiguous companions in our sample, but identify two targets as candidate binary systems: HIP 92855 and HIP 26713. HIP 92855 is B2.5V type star, with a candidate companion star with effective temperature T = 6100 K and mass 1.2 ± 0.2 M_☉. Such a companion star is very near the detection limit of the direct spectral detection technique and deserves further follow-up observations. HIP 26713 is a B1.5V type star, with a candidate companion star with T = 5600 K and M = 1.6 ± 0.2 M_☉. This star was only observed once in our sample, and so may be a series of incompletely removed telluric absorption lines masquerading as a companion star.

We set upper limits on the binary fraction of early B stars as a function of binary mass ratio (see Figure 7). As well as showing the upper limit for no detections in our sample, we also show the binary fraction upper limit assuming the HIP 92855 is a true binary system. Our upper limits are strongest for mass ratios q ≈ 0.1–0.15, and are about 20%. We compare our limits to the intrinsic binary mass-ratio distribution for O-type primaries derived by Sana et al. (2012), as well as the lower limit average binary fraction seen by Evans et al. (2011) for late B stars (B4–B9). Our strongest upper limits (0.1 ⩽ q ⩽ 0.15) fall in between these two previous studies.

Companion stars formed by circumstellar disk instability are expected to have typical mass ratios near q = 0.1 (Kratter & Matzner 2006; Stamatellos & Whitworth 2011), near where our upper limits are strongest. If there was a large population of low mass-ratio companions formed by disk fragmentation, we would expect to see a peak in the mass-ratio distribution near q ≈ 0.1. Since our results agree well with the nearly flat mass-ratio distribution derived by Sana et al. (2012), it is unlikely that such a peak exists. There are several possible interpretations of this result, three of which we give below.

• 1.  
  Stellar companions formed by disk instability have a much lower characteristic mass ratio than q = 0.1, which remain invisible to observations.
• 2.  
  The mass-ratio distribution of companions formed by disk instability is very broad, and so we would not expect a strong peak in the observed mass-ratio distribution. This interpretation may be supported by the flatness of the mass-ratio distribution as well as disk instability simulations that end with massive companions (e.g., Krumholz et al. 2009; Clarke 2009).
• 3.  
  Disk instability is not a dominant formation mechanism for low-mass-ratio binary systems, and molecular core fragmentation alone can generate the nearly flat mass-ratio distribution down to low mass ratios.

It is difficult to distinguish between these interpretations at this time. More observational work, as well as computational work, is required to explain the binary properties of high-mass stars.

We acknowledge Andreas Seifahrt, Rob Robinson, and Daniel Jaffe for their generous help with the telluric modeling procedure and some of the statistical aspects of this work. We also thank the referee for a quick review and many helpful comments. This research has made use of the following online resources: the SIMBAD database and the VizieR catalogue access tool at CDS, Strasbourg, France, and the ESO Science Archive Facility. Funding for this work was provided by a National Science Foundation CAREER award to Sarah Dodson-Robinson (AST-1055910) and by start-up funding from the University of Texas College of Natural Sciences.