MULTIPLE STAR FORMATION TO THE BOTTOM OF THE INITIAL MASS FUNCTION

Nearby star-forming regions have been the target of numerous wide-field photometric imaging surveys to detect new low-mass members (e.g., Luhman 2004, 2006; Slesnick et al. 2006a, 2006b, 2008). These surveys identified candidate members based on their location on an optical or near-infrared color–magnitude diagram, and membership was then confirmed spectroscopically via the detection of lithium absorption, excess Hα emission, or low surface gravity, all of which are indicators of youth. We chose to concentrate on Taurus-Auriga (τ ∼ 1–2 Myr, d ∼ 145 pc; Kraus & Hillenbrand 2009a; Torres et al. 2009) and Upper Scorpius (τ ∼ 5 Myr, d ∼ 145 pc; de Zeeuw et al. 1999; Preibisch et al. 2002) because they are the nearest young associations that are accessible from the northern hemisphere. In our analysis, we also use results from the literature on the Cha-I association (τ ∼ 2–3 Myr, d ∼ 185 pc; Luhman 2004), which is similar to Taurus in age and environment. All regions have very low stellar densities (N ≲ 10 pc⁻³), which is critical for minimizing the complicating role of dynamical interactions.

Our initial observational sample included all late-type members of each association (SpT ⩾ M4) that had been identified by 2006 and that had not been observed at high angular resolution. However, as we describe in Section 3, we lost a significant fraction of our observing time to poor weather and instrument problems. This left our Taurus sample significantly incomplete for members discovered in 2006, plus we were unable to observe three members that had been identified earlier (J1-4423, V410 X-ray 6, and 2MASS J04163049). The effect on Upper Sco was to limit our observed sample to only the latest-type members (SpT ⩾ M6.5). We also were unable to observe 10 Taurus members that had no suitable tip-tilt stars available. This omission introduces a bias against the most reddened members of Taurus since the most heavily extincted stars were least likely to have an optically bright star nearby that could serve as a tip-tilt reference. The density of field stars is very high in Upper Sco, plus the association is almost completely cleared of its primordial molecular material, so we were always able to find a suitable tip-tilt star.

We have supplemented this observational sample with the results of numerous previous multiplicity surveys. Taurus has been a very popular target for multiplicity surveys, and association members in our spectral type range have been observed with speckle interferometry (Ghez et al. 1993; Konopacky et al. 2007), lunar occultations (Simon et al. 1995), Hubble Space Telescope (HST) imaging (Padgett et al. 1999; White & Ghez 2001; Kraus et al. 2006), adaptive optics (AO) imaging (Correia et al. 2006), and aperture-masking interferometry (Kraus et al. 2011). Upper Scorpius has been the subject of several surveys as well, and members have been observed with speckle interferometry (Köhler et al. 2000), HST imaging (Kraus et al. 2005), and AO imaging and aperture-masking interferometry (Kraus et al. 2008). In most cases (and almost certainly in aggregate), the sample members were selected seemingly at random. As a result, we adopt the combined set as an unbiased sample.

Finally, we also observed a small number of other targets that fall outside these selection parameters, but were considered interesting for other reasons. In both associations, we observed a number of candidate wide binary systems that seemed to have unusually low binding energies. We already described the observations for UScoJ1606-1935 in a previous paper (Kraus & Hillenbrand 2007c), and we reported the astrometric measurements for the rest in our paper on wide binary formation (Kraus & Hillenbrand 2009b). In this paper, we report on the search for higher-order multiplicity. We also observed several candidate Taurus members discovered by Slesnick et al. (2006b) that are not part of the young Taurus population, but might represent an older, more widely distributed population of young stars and brown dwarfs. Finally, we observed the known binary V928 Tau because it served as the tip-tilt reference for CFHT-Tau-7 and we typically imaged tip-tilt references for a data quality check. However, we do not report any results for CFHT-Tau-7 because the observing conditions at the time were too marginal for AO to yield any meaningful correction.

In Table 1, we list the young association members that we observed in our study. The K magnitude for each target was taken from Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006), while the R magnitude and distance to each star's tip-tilt reference are from the USNO-B1.0 catalog (Monet et al. 2003). We also list references for the handful of objects which have also been observed in other high-resolution imaging surveys; in several cases, our detection limits for small separations were superceded by the survey by Konopacky et al. (2007), so we adopted those limits where appropriate.

Most of the data that we summarize were obtained in four observing runs, totaling 10 nights, between 2005 December and 2007 March. One source was observed during a time trade in 2006 December. Most of our observations were obtained using laser guide star adaptive optics (LGSAO; Wizinowich et al. 2006) on the Keck-II telescope with NIRC2 (K. Matthews, in preparation), a high spatial resolution near-infrared camera. During some periods of moderate cloud cover that were not suitable for laser operation, we also used natural guide star adaptive optics (NGSAO) to observe sample members with very close and bright tip-tilt stars. In the worst conditions, we also observed some higher-mass stars that did not fall in our sample; most of these observations have been described in our previous papers (Kraus & Hillenbrand 2008, 2009b), so we report the rest here for completeness.

The weather conditions were highly variable over the course of our campaign, with only five nights suitable for laser operations. Several of the remaining nights were also impacted by poor seeing. Despite numerous difficulties, we report 82 observations of young association members with LGSAO and 5 additional observations with NGSAO, encompassing 78 different targets. In Table 2, we summarize the observations for each target. We also list the typical point-spread function (PSF) FWHM for each target; a significant fraction of the targets used bright, on-axis tip-tilt stars that should have achieved diffraction-limited performance (strehl ∼ 30%), but we only achieved this performance on two nights.

All of the images presented here were produced with the narrow camera, which has a field of view (FOV) of 102 and a pixel size of 9.963 mas pixel⁻¹ (Ghez et al. 2008). All targets were observed with the K_p filter; in most cases, we did not obtain observations in other filters because most background stars have J − K and H − K colors that are not sufficiently different from young stars as to allow secure identification. Many results we draw from the literature appear to have been observed with K or K_s filters (though some do not specify), but we treat all three filters equivalently since the color terms are smaller than the typical photometric uncertainties (e.g., Carpenter et al. 2002). During early observing runs, we used a three-point dither pattern that was designed to avoid the bottom-left quadrant, which suffers from high read noise. After 2006 February, we obtained all of our observations in a diagonal two-point dither pattern because experience showed that dithers degrade the AO correction until several exposures have been taken with the low-bandwidth wavefront sensor. The delay before returning to optimal correction represented a significant overhead that we sought to minimize.

Many of the targets are relatively bright in the NIR and require very short integration times to avoid saturation or nonlinearity, so a large fraction of our observations were taken in correlated double-sampling mode, for which the array read noise is 38 electrons read⁻¹. Where possible, we observed targets in multiple correlated double-sampling mode, where multiple reads are taken at the beginning and ending of each exposure; this choice reduces the read noise by approximately the square root of the number of reads. This is doubly significant because the read noise per co-add and the total number of co-adds per exposure are both reduced. In all cases, the read noise is negligible compared to PSF variations from the primary at separations of <1''. However, it always dominated over the sky background in determining our faint-source detection limits at large separations from the science target. In all cases, the images were flat-fielded and dark- and bias-subtracted using standard IRAF procedures.

Source identification in AO imagery is a complicated endeavor. In NGSAO mode, the gross shape of the PSF depends on the target's optical brightness and the seeing, while the fine structure is determined by speckle patterns that continuously change on timescales ranging from seconds to hours. The LGSAO PSF is further complicated by variations in laser return strength, tip-tilt anisoplanatism with respect to off-axis guide stars, and heightened sensitivity to telescope effects like wind shake. Finally, observations in poor weather are further complicated by rapid PSF quality variation due to changing atmospheric conditions. The source detection process can be divided into two regimes: a wide regime where the PSF core is negligible and speckle confusion dominates (projected separation ρ ≳ 2 times the core FWHM), and a close regime where shape and width variations in the PSF core dominate and speckle confusion is negligible. We have adopted a different method in each separation regime.

3.1.1. The Wide, Speckle-dominated Regime

As was summarized by Metchev & Hillenbrand (2009), there are four common methods used to subtract the primary star's flux and identify companions in AO imagery: subtracting a median PSF representing all similar observations, subtracting a 180° rotated version of the same image, high-pass filtering by subtracting a Gaussian-smoothed version of the same image, or subtracting the azimuthally averaged profile. We conducted experiments with these techniques, but we found that the speckle mitigation strategies that are vital for high-strehl NGSAO data are actually only marginally useful for low-strehl LGSAO data. Most of the flux that would be found in discrete speckles in high-strehl data is instead averaged into the seeing-limited halo, rendering the PSF less azimuthally variable at a given separation. The brightest speckles remain distinguishable, but at far lower contrast with respect to the surrounding median flux. Since the noise floor is brighter and the noise ceiling is fainter, there is less to be gained from exceeding the noise ceiling.

The sparsity and relatively low contrast of speckles in LGSAO data suggest that a different strategy is optimal for our data. Speckle mitigation and subtraction of the primary star's flux are observationally expensive, so the preferred strategy should be to characterize the mean and standard deviation of the brightness distribution of the PSF as a function of separation, then set the source detection limit above the expected ceiling for speckle brightness. We characterized the brightness distribution of each target's PSF by measuring the flux through photometric apertures placed at a range of separations and position angles (P.A.s), then measuring the mean and standard deviation for all apertures in a given bin of separation. The apertures were placed on a rectangular grid with spacing of 25 mas in order to ensure that the small number of speckles were detected, and the aperture sizes were matched to the FWHM of the PSF core for the primary. We measured this aperture photometry using the IRAF⁴ task PHOT, which is distributed as part of the DAOPHOT package (Stetson 1987).

In Figure 1, we show the contrast as a function of separation for three stars that span our survey's data quality, as well as the 5σ envelope for each source. We found that there was typically one detection at 4σ–5σ per 2 stars, one detection at 5σ–6σ per 10 stars, and no detections among any of our targets at 6σ–10σ. This indicates that a 6σ clip should be uncontaminated by spurious detections, while a 5σ clip can be adopted if the few remaining speckles can be confidently identified as such. We found that all of the 5σ–6σ candidate detections fell on the PSF's diffraction spikes, were sufficiently short-lived so as to not appear in all observations of a target, or were sufficiently long-lived so as to appear in observations of multiple sequential targets. We therefore suggest that all such 5σ–6σ candidate detections are spurious, and adopt a 5σ clip as our survey's detection limits. All candidate detections that sit well above the 5σ limit appear to be genuine astronomical sources, though not necessarily comoving companions; we will revisit this distinction in Section 5.2.

Zoom In Zoom Out Reset image size

3.1.2. The Close, Core-dominated Regime

For separations near the PSF FWHM, the detection limits are driven by time-dependent and spatially dependent variations of the shape of the PSF core. LGSAO observations seem to be more susceptible to all of the weather effects that can degrade NGSAO observations, so distinguishing genuine companions from PSF artifacts is a significant challenge. The primary effect we see is for wind shake to cause PSF elongation in the direction of the zenith, perhaps because the tip-tilt sensor and low-bandwidth wavefront sensor operate at a lower rate and can not fully sample high-frequency oscillations of the telescope. Tip-tilt anisoplanatism is also significant for observations with off-axis tip-tilt guide stars, an effect that becomes worse in poor seeing because the isokinetic angle becomes smaller. This causes elongation along the P.A. to the tip-tilt star. Finally, significant variations in the AO correction cause the PSF FWHM itself to vary by a factor of ∼3 across our sample; a few of the lowest-quality observations have a PSF FWHM approaching 150 mas.

We have characterized these effects by fitting each science target's PSF core with a bivariate Gaussian distribution. This fit directly yields the PSF elongation (the ratio of the major axis a and minor axis b) and its direction (the P.A. of the major axis). For each of our targets, we report the minor axis FWHM (i.e., the resolution prior to elongation effects) and the fractional elongation in Table 2. If wind shake and tip-tilt anisoplanatism are the dominant sources of PSF asymmetry, then most sources should have a P.A. that is preferentially aligned with either the tip-tilt angle or the zenith angle. A set of selection criteria based on these quantities also lends itself to rigorous characterization of the detection limits, as artificial star tests can be used to determine whether a system with given separation and contrast (and perhaps P.A.) would be detected.

In Figure 2 (top), we plot the PSF elongation and the relative angle between the PSF and the zenith angle (|θ_PSF–θ_zen|) for all science targets which served as their own tip-tilt references. Of the 21 sources which are not independently confirmed binary systems, 9 are aligned to within <10° of the zenith angle and 7 of the remaining 12 are aligned to within 10°–45°. This strong trend indicates that wind-induced elongation was significant across the majority of our sample, even on those nights with moderate winds. Also, all of the targets that are not independently confirmed binary systems have elongations of <30%, which seems to be the ceiling for PSF elongation due to telescope or atmospheric effects.

Zoom In Zoom Out Reset image size

In Figure 2 (bottom), we consider the rest of our sample in plotting the PSF elongation versus the minimum of the relative angle either between PSF and zenith (|θ_PSF–θ_zen|) or between PSF and tip-tilt (|θ_PSF–θ_TT|). These targets also show a pronounced tendency to align with either the zenith or tip-tilt, though the result is more complicated because many targets have a net elongation intermediate between the two directions. There are fewer confirmed binary systems for comparison, but most of the targets fall below elongations of 40%, suggesting that this is the ceiling for combined wind and tip-tilt effects. Many of the targets with elongation >40% fall among our poorest sample and have significantly different elongation angles and magnitudes in each exposure, while all of the remainder present consistent and apparently double-peaked PSFs.

In light of these trends, we have adopted two criteria for identifying a source as a candidate binary system. First, it must have a PSF elongation of >40%, which appears to be the ceiling for any weather-based effect in all but the worst data. Second, similar PSF elongation must not be seen for other sources in the science target FOV, in our preliminary image of the tip-tilt guide star, or in observations of the previous or subsequent science target. Finally, the astrometry and photometry for a fit of two-point sources must be consistent across all exposures; the poor data with elongations >40% tend to vary its PSF shape on extremely short timescales, yielding extremely inconsistent fits across the full data set. We have inferred the source detection limit for each of our targets by measuring its minor axis FWHM, then using artificial star tests to determine what ranges of companion separation and companion brightness would have elongated a circular PSF with that FWHM to >40%.

Finally, there are also some cases where companions can be confidently studied below our survey's detection limit, such as if the companions were previously identified in another survey (i.e., V410 X-ray 3) or if a third bright star can be used as an independent PSF calibrator (UScoJ1607-2019). We have used PSF-fitting techniques (Section 3.2) to recover the photometry and astrometry for these close binary pairs, though we generally cannot include them in our statistics if their detection relied on a special feature of the system like high-order multiplicity. We also note that we were unable to recover accurate astrometry and photometry for MHO-Tau-8, which suggests that its orbital motion might have carried it inward from its last-known projected separation (∼40 mas; Kraus et al. 2006).

We measured relative photometry and astrometry for candidate companions using the IRAF package DAOPHOT (Stetson 1987). For source pairs with small separations, where the two PSF cores were not unambiguously resolved, we used the PSF-fitting ALLSTAR routine. For pairs with wider separations, we used the aperture photometry package PHOT. We analyzed each frame separately in order to estimate the uncertainty from the scatter between all frames; this also allowed us to reject some frames with subpar AO correction. Our final results represent the mean and standard deviation for all observations.

For the close binaries that we analyzed with ALLSTAR, we reconstructed the single-star PSF out of the merged binary PSF using the algorithm described in Kraus & Hillenbrand (2007c), which iteratively fits a template PSF to the primary and then subtracts the secondary to fit an improved estimate of the primary. For one triple system, UScoJ1607-2019, we chose instead to use the single secondary as a PSF template for fitting the close pair constituting the primary. This choice allowed us to clearly distinguish the close pair despite a separation (∼50 mas) that was significantly lower than the PSF FWHM (∼70 mas).

Our relative astrometric measurements were distortion-corrected using a new high-order distortion solution (Cameron 2008) that delivers a significant performance improvement as compared to the solution presented in the NIRC2 pre-ship manual.⁵ This distortion solution was derived from observations of a pinhole mask in the NIRC2 filter wheel, so it does not include any distortions introduced upstream of this point. The remaining residuals due to these uncorrected distortions are ∼5 mas for positions separated by ∼5–10'' (J. Lu 2011, private communication). We calibrated our photometry using the known 2MASS K_s magnitudes for each of our science targets; these absolute magnitudes are uncertain by ∼0.1–0.2 mag due to the intrinsic variability of young stars (resulting from accretion or rotation).

Stellar properties can be difficult to estimate, particularly for young stars, since pre-main-sequence stellar evolutionary models are not well calibrated. The model-predicted masses of young stars could be systematically uncertain by as much as 20% (e.g., Hillenbrand & White 2004), and estimates for individual stars could be uncertain by factors of two or more if their observed luminosities are biased by unresolved multiplicity or the intrinsic variability that young stars often display (from accretion or from rotational modulation of star spots). This suggests that any prescription for determining stellar properties should be treated with caution.

We estimated the properties of all of our sample members using the methods described in Kraus & Hillenbrand (2007a). This procedure combines the 2 or 5 Myr isochrones of Baraffe et al. (1998) and Chabrier et al. (2000) with the temperature scales of Schmidt-Kaler (1982) and (Luhman et al. 2003) to directly convert observed spectral types to masses. Relative properties (mass ratios q) for all binaries in our sample were calculated by combining these isochrones and temperature scales with the empirical NIR colors and K-band bolometric corrections of Kraus & Hillenbrand (2007b) to estimate q from the observed flux ratio ΔK'.

For all binary systems, we have adopted the previously measured(unresolved) spectral type for the brightest component and inferred its properties from that spectral type. This should be a robust assumption since equal-flux binary components will have similar spectral types and significantly fainter components would not have contributed significant flux to the original discovery spectrum. Projected spatial separations are calculated assuming the mean distance to the associations (∼145 pc for both Upper Sco and Taurus; de Zeeuw et al. 1999; Torres et al. 2009). If the total radial depth of each association is equal to its angular extent (±8° or ±20 pc), then the unknown depth of each system within the association implies an uncertainty in the projected spatial separation of ±14%. The systematic uncertainty due to the uncertainty in the mean distance of each association is negligible in comparison (≲2%).

Our search for sources in the speckle-dominated regime (at separations ≳1.5 times the PSF FWHM and extending to the edge of the detector, Section 3.1.1) yielded 45 candidate companions among the 78 young stars and brown dwarfs in our observed sample. All candidates within ≲1'' of the target sit well above the 5σ detection limit, so they all represent secure detections and do not seem to be spurious structures in the primary star's PSF. We also found numerous possible detections with significance levels of 5σ–6σ in this regime, but as we described in Section 3.1, all of these possible detections appear to be spurious. Our corresponding search for sources in the core-dominated regime (at separations of order the PSF FWHM, Section 3.1.2) yielded nine targets with PSF cores consistently elongated by >40%. Many targets had PSF elongations below this limit, but as we discussed above, most appear to be distorted due to observational effects and not the presence of a companion. We address the membership probabilities of these candidates in Section 4.2.

In Table 3, we list our survey's candidate companions and report their flux ratios, separations, and P.A.s. We also plot the flux ratio ΔK' and the candidate companion brightness K' as a function of separation in Figures 3 and 4, respectively. Finally, in Figure 5, we show contour plots for the 12 candidate companions with separations <1'' and flux ratios of ΔK' < 4; we do not show contour plots for wider systems because they are completely resolved, and we do not show any faint companions because few are likely to be bound companions. The vast majority of our wide-separation candidates are near the detection limits of our survey, in the brightness range where a real companion would fall below the deuterium-burning limit (∼13 M_Jup), so we expect that almost all are unassociated field stars. However, a handful of planetary–mass companions have been identified around young stars and brown dwarfs (e.g., 2M1207b; Chauvin et al. 2004), so we must consider the possibility that some of the candidates in our survey are also extremely low mass companions. We will present second-epoch observations that test the membership of these candidates in a future paper. In all cases, the high astrometric precision of Keck/NIRC2 (σ ∼ 1–2 mas) and the relatively large proper motions (μ ∼ 20–30 mas yr⁻¹) allow for a test of common proper motion using observations from consecutive observing seasons.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Of the 54 combined sources from the speckle-dominated and core-dominated regimes, 11 had already been identified as candidate binary companions by past survey efforts, so we note these past identifications in Table 3. Three of our candidate companions merit special attention. V410 X-ray3 AB was identified as a candidate binary in our previous HST/ACS multiplicity survey based on a marginal elongation of its PSF, but the best-fit separation was well inside the HST diffraction limit for the i' and z' filters. By contrast, the system is almost resolved in our K' observations, and corresponding JH images (Kraus 2009) clearly reveal V410 X-ray3 to be a genuine binary system. Two wide binary systems, USco 80 AB and UScoJ1607-2019 AB, were resolved to be hierarchical triples. USco 80 A is clearly resolved to be a close pair, while UScoJ1607-2019 A is marginally resolved.

The other 26 members of our sample have no resolved neighbors within our detection limits. These detection limits, which we derived using the methods described in Section 3.1 and list in Table 4, are extremely heterogeneous due to the wide range in observing conditions during our observing campaign. Some observations nearly achieve the expected limits for diffraction-limited images, while images from most of the other nights achieved significantly poorer conditions. We also note that our nominal detection limits for bound companions at wide separations was constrained by our follow-up efforts. We found many more faint candidates than we were able to follow up with second-epoch imaging, so we cannot claim completeness beyond the maximum separation or flux ratio at which we have identified all of the field stars by testing for common proper motion.

Companion searches must address the prospect of chance alignments with background stars, especially surveys with very deep detection limits. We are obtaining multi-epoch astrometry for several candidates in order to test for common proper motion, but we can use statistical arguments to determine which candidate companions require those follow-up observations. As we describe more fully in Appendix A, we have updated the Milky Way model of Bahcall & Soneira (1980) to predict star counts as a function of magnitude for the line of sight toward each of our targets; these models allow us to predict the field star contamination rate and thereby determine which companions have a significant probability of being background stars.

In Figure 4, we plot the joint magnitude–separation limits at which our models predict we should find 1, 3, 5, or 10 background stars among all the targets observed in that association. In both associations, our models predict that we should find <1 background star with K ≲ 15 within <5'', which suggests that all of the bright sources we observe well inside this limit are genuine companions. This limit agrees with our estimate based on 2MASS source counts in the direction of Taurus and Upper Sco (Kraus & Hillenbrand 2007a), which found that all neighbors down to the 2MASS 10σ limit (K = 14.3) could be assumed to be bound association members out to separations of 3''–5''.

The status of our fainter candidates is not as clear. Our models predict that we should find only one background star with separation ρ < 2'' and brightness K  < 19 in Taurus, so the candidates inside this limit all seem very promising as potential analogs to 2M1207b. The background star density is higher in the direction of Upper Sco since it lies in the direction of the galactic center, so even with our smaller sample size, we still expect three chance alignments with ρ < 2'' and K <19. However, star count models are not well constrained at K ≳ 14 since the counts are dominated by low-mass halo stars, a population that is not as well studied as the brighter disk stars. As a result, the contamination rate for faint stars could be uncertain by a factor of at least two to three. It would be prudent to measure common proper motion for any candidate companion fainter than K ∼ 14, especially since any genuine companion would have M  < 10–20 M_Jup, making it an extremely compelling discovery. The predicted background star contamination rate rises quickly for separations ≳3'', matching the many candidates we have discovered, so we provisionally adopt this separation as an outer limit at which it is worthwhile to test candidates for possible association. Future searches at wider separations should use seeing-limited data from publicly available surveys like UKIDSS, which will observe the majority of these targets in four NIR filters.

Finally, multiplicity surveys in young clusters and associations must also consider chance alignments between two unbound member stars. These chance alignments are extremely difficult to distinguish from genuine binary systems since all association members are young, at similar distance, and comoving to very high precision. The only solution is to treat their probability in a statistical sense, which we have already done in our treatment of wider binary systems (Kraus & Hillenbrand 2009b). We found that the probability of a chance alignment between two young association members is negligible for separations <10'', so we will proceed under the assumption that any pair of young stars constitutes a genuine binary system.

In Table 5, we list the mass ratios and component masses that we infer for the candidate companions brighter than K = 14 and closer than 3'' from their primary star, which we henceforth consider to be bound binary companions. For targets that were observed at multiple epochs, we list the independent estimates from each epoch. These properties were derived using the methods we describe in Section 3.3.

Binary population statistics are traditionally presented in terms of histograms of binary frequency versus separation or mass ratio, where the data is presented only for a range where the survey is complete. The analytic form of the preferred model is then fit to these histograms in order to infer the population properties. This approach has the virtue of simplicity, but estimating the probability density function (PDF) for the model's scale parameters is often difficult, especially if there are covariances between parameters. This method also is manifestly inadequate for handling heterogeneous data sets. If different stars have different detection limits, such as from being observed with different methods or under different atmospheric conditions, then simple histograms can be constructed only by appealing to completeness corrections that are themselves poorly constrained.

A better solution for working with heterogeneous data is to adopt a Bayesian approach, where the scale parameters of the model are assigned a prior PDF and that PDF is modified by each observation. This method exploits Bayes' theorem:

$$\begin{equation} P(\theta |O) \propto P(O|\theta) P(\theta), \end{equation} \tag{ 1 }$$

where θ represents the "model" (a set of scale parameters describing the functional form), O represents the observation, P(θ|O) is the posterior PDF for the model (as a function of its parameters) given the data, P(O|θ) is the probability of obtaining an observation as a function of the model parameters, and P(θ) is the prior PDF for the model (again, as a function of its parameters). In cases with multiple observations (such as a survey of many targets), the posterior function for one observation is then used as the prior function for the next observation.

Allen (2007, hereafter A07) developed the relevant techniques for applying Bayesian statistics to VLM multiplicity, and we have largely adopted his approach in this work and a parallel survey of solar-type stars in Taurus (Kraus et al. 2011). We specifically describe the binary population in terms of a binary fraction F, a power-law mass ratio distribution with exponent γ, and a log-normal separation distribution with mean log (s) and standard deviation σ_log (s). We have adopted the same Poisson likelihood function as A07, but we will use a moderately different set of prior distributions. As for A07, we use constant priors for γ and log (s). However, A07 seem to state that scale invariance requires the optimal unbiased prior for σ_log (s) to be proportional to 1/σ_log (s), whereas it actually should be proportional to 1/σ_s, or constant in σ_log (s).

The difference for F is more subtle. A07 used a prior proportional to 1/F, which was suggested to be suitable for a Poissonian variable. As we discussed in Kraus et al. (2011), the Poisson case is sometimes appropriate if there are many high-order multiples, though our newest analysis actually uses the Jeffries prior that is proportional to $1/\sqrt{F}$ (A. L. Kraus et al., in preparation). However, binary companions exclude other binary companions in similar orbits (spanning at least a decade of semimajor axis), so they are not genuinely Poissonian. As we show below, ultracool binaries typically only have small projected separations, perhaps not spanning much more than a decade of semimajor axis, so the exclusionary effect should be quite significant. We are attempting to use completely uninformed priors for this analysis, but a more sophisticated prior that declines more quickly than $1/\sqrt{F}$ should be considered for future analyses.

Another significant difference between our analysis and that of A07 is that we directly model the projected separation distribution, whereas A07 used the semimajor axis as a model parameter and then extrapolated a projected separation distribution using an assumed eccentricity distribution and randomly distributed inclinations and phase angles. The eccentricity distribution for binary systems is still largely unconstrained for the separations and mass ranges that we consider, but the most likely eccentricity distributions yield a separation distribution that is directly proportional to the semimajor axis distribution. As such, our results can be directly related to theoretical semimajor axis distributions once the eccentricity distribution is predicted by theory or measured by future surveys. We have also omitted the volume-completeness correction used by A07 to compensate for the overluminosity of similar-brightness binaries. The discovery surveys for most of our sample members were spatially limited, not flux- or volume-limited, so binary systems were equally likely to be detected. The high-resolution imaging techniques used in past surveys were themselves flux-limited, but we chose our LGSAO sample in part to compensate for this limit, so it should not significantly influence our results to invoke detections and detection limits from those past surveys where needed.

Our specific implementation of Bayesian analysis follows that of A07, defining two-dimensional (2D) functions of projected separation log (s) and mass ratio q that denote the number of observations sensitive to each set of log (s) and q (the "window function") and the corresponding number of companions with that set of parameters. We iterated our calculation over all mass ratios from 0 to 1 in steps of 0.01 and over all values of log (s) between 0.5 and 3.6 dex in steps of 0.1 dex; this range of log (s) was chosen to encompass the full range of projected (physics) separations for which imaging observations are sensitive, from 3 AU (due to the minimum resolution limit of the most sensitive surveys; Kraus et al. 2005, 2006) to 4500 AU (the maximum projected separation for which we have eliminated background star interlopers for most targets; Kraus & Hillenbrand 2007a, 2009b).

Our own sample comprises the vast majority of available measurements in the VLM regime (M ≲ 0.15 M_☉), but almost all of the measurements for higher-mass stars must be adopted from previous surveys. As we describe in Appendix B and list in Table 7, we have specifically adopted the detections and detection limits for all stars with M ≲ 0.5 M_☉ for previous surveys of Taurus-Auriga (Ghez et al. 1993; Simon et al. 1995; Kraus et al. 2006; Konopacky et al. 2007), Upper Sco (Köhler et al. 2000; Kraus et al. 2005; Biller et al. 2011), and Cha-I (Ahmic et al. 2007; Lafrenière et al. 2008). In each case, we converted the measured angular separations (in mas) and flux ratios (in Δm) for their detections and detection limits into physical quantities (separations in AU and mass ratios) using the same methods that we applied to our own sample. This ensures a uniform set of inputs for our analysis, whereas each survey's inferred system properties were derived using different combinations of association distances/ages and pre-main-sequence stellar evolutionary models.

Past studies have indicated that multiple star formation might be a mass-dependent process, so we have divided our sample into several bins. The stellar/substellar boundary represents a natural breaking point since it corresponds roughly to the M/L boundary for field objects, allowing for natural comparison to field samples. As we have indicated in our studies of solar-type multiplicity and wide multiplicity, systems with primary masses of ≳0.5 M_☉ tend to have fundamentally different binary parameters, featuring a log-flat separation distribution and many very wide systems, so we have adopted this limit as the maximum for consideration in our sample. The desire for similar number statistics per bin therefore dictated four mass bins of similar width: <0.07 M_☉, 0.07–0.15 M_☉, 0.15–0.30 M_☉, and 0.30–0.50 M_☉. A more rigorous treatment might incorporate a mass dependence directly into our fit parameters (F, γ, log (s), and σ_log (s)), but in the absence of any theoretical guidance on the functional form of this mass dependence, we will defer such analysis to a future study.

Our Bayesian analysis yields a PDF for all possible "models" that is defined across four dimensions, so we cannot present the full results in a 2D medium. However, any two parameters for which the covariance is small can be presented separately without discarding information. This independence allows us to instead present the results as a series of lower-dimensional surfaces, where PDF is integrated across the uncorrelated parameters in order to flatten its dimensionality. As we describe in the next section, our results can be described with a manageable number of 2D or 1D surfaces.

Finally, it is impossible to define a true PDF with only null detections, so we cannot use this analysis for the lowest-mass bin (M < 0.07 M_☉) since it includes no resolved binary systems. Our choice to use generic conjugate priors does not disallow arbitrarily extreme values, such as small mean separations or steep power laws. If we do not have enough constraints (i.e., detected binaries) to force the PDF to zero at all extrema of the scale parameters, then the integrated probability will diverge and render the PDF unnormalizable. Since we cannot estimate a well-defined probability for any particular set of parameters being "correct," we instead must settle for a weaker result: the probability that a given set of parameters would have yielded our null detection. This measurement is equally valid for ruling out parameter space, but does not carry any explicitly affirmative value; regions where the model is less improbable are not necessarily regions where the model is probable. The act of "flattening" the PDF to visualizable 1D or 2D figures is also not defined for this type of constraint, so in figures where the PDF for higher-mass bins is flattened, we will instead show a cross-section through the lowest-mass PDF where we adopt the field T dwarf parameters suggested by Burgasser et al. (2006): log (s) ∼ 0.6, σ_log (s) ∼ 0.3, and γ ∼ 4.2.

The four-dimensional posterior PDF for our Bayesian analysis can be flattened to present six 2D probability surfaces and four 1D probability curves, but to convey the useful conclusions, we only need surfaces for covarying parameters and curves for non-covarying parameters. For our results, we will present two probability surfaces (F versus log (s) and F versus σ_log (s)) and one probability curve (γ).

In general, we found little covariance between γ and any other parameter, which is largely a result of the shallow mass–luminosity relation for young stars and brown dwarfs; at a given separation, most observations are either unable to detect any companions or sensitive to companions with almost all mass ratios. There is significant covariance between the binary frequency F and the two parameters in the separation distribution, log (s) and σ_log (s). This degeneracy results from the inner working angle for most of the input data sets (∼5–10 AU) being of similar order as the mean separation, since our measurements are consistent with a range of binary frequencies as long as an appropriate fraction of the companions are "hidden" inside the detection limit with a smaller mean separation and correspondingly wider standard deviation. As we discuss further in Section 5.3, a comparison to similar constraints from simulated stellar populations (e.g., Bate 2012) also can provide a direct test of star formation models and indicate future directions for their improvement.

5.2.1. Frequencies and Separation Distributions for Stellar Multiplicity

In Figure 6, we show the PDFs for our three stellar-mass bins as projected onto the F-log (s) and F-σ_log (s) planes. In the two lower-mass bins (0.30–0.15 M_☉ and 0.15–0.07 M_☉), there is a significant degeneracy between the overall binary frequency and the mean separation, where a smaller mean separation is paired with a higher frequency. This degeneracy is unavoidable when fitting a normal distribution whose mean is near or outside the fitting region; the first derivative (i.e., the slope) of the distribution across the fitting region yields the standard deviation with little ambiguity, but distinguishing between the total amplitude of the curve and the distance (in standard deviations) to the mean requires measurements of both the number of measurements and the second derivative (i.e., the change in slope) across the fitting region. Measuring each successive derivative requires either more S/N or a wider fit regime. The highest-mass bin also shows some degeneracy with separation, but not to the same extent since its mean separation is outside of the typical inner working angles for many of the input surveys.

Zoom In Zoom Out Reset image size

The F-log (s) locus in the highest-mass bin (0.3–0.5 M_☉) is clearly distinct from the loci of the two lower bins (0.30–0.15 M_☉ and 0.15–0.07 M_☉), as its 90% confidence region does not overlap with the same regions for the other bins. This indicates that the mean separation and/or the binary frequency are significantly higher for 0.3–0.5 M_☉ stars. The strong degeneracies seen for the lower-mass bins make it difficult to draw any strong conclusions, but it appears that the frequency and/or mean separation for binary systems declines from the 0.15–0.30 M_sun regime to the 0.07–0.15 M_☉ regime. As we discuss in Section 5.2.2, this decline seems to continue in the substellar bin, though our constraint on its relative magnitude in each parameter is even weaker since we only have a null detection for that mass bin. The F-σ_log (s) loci are not as easy to interpret; the middle locus is biased to a larger standard deviation, but this might be the result of having four probable binary systems with separations of >1000 AU; if some of these pairs of stars are actually chance alignments of unrelated association members, then removing them would reduce the standard deviation and the mean separation by a significant amount.

The frequency–separation degeneracy must be addressed before we can draw any stronger conclusions. The most direct solution would be to increase the number statistics in our existing program, yielding a better estimate of the high-order derivatives in the separation distribution. However, this endeavor would be very observationally expensive; simulations show that even doubling our sample would not decrease the length of the degenerate locus, only its width. A less direct solution would be to expand the range of separations over which the distribution is constrained, either by observing at higher resolution (sampling more of the core separation distribution) or by searching for spectroscopic binaries (constraining the other wing of the separation distribution). We suggest that an RV survey would be significantly cheaper since it can exploit the multiplexing of wide-field multi-object spectrographs, plus the separation–frequency degeneracy that results from an RV survey's outer working angle should be perpendicular to the degeneracy from an imaging surveys' inner working angle. Such surveys are currently being pursued for nearby young populations like the sigma Ori cluster (Maxted et al. 2008) and the ONC (Tobin et al. 2009), and their results could be modeled with similar Bayesian techniques in order to produce constraints analogous to those shown in Figure 6. However, such a modeling effort is beyond the scope of the current work.

Finally, although Bayesian inference allows us to estimate the most general limits on the binary population, the degeneracies in those limits make it difficult to straightforwardly grasp the differences in our subsamples. We address this by forward-modeling from our four-dimensional PDF back into the range of separations and mass ratios where our observations could detect companions around most of our targets. The net result of this extrapolation is to implement a minor correction for incompleteness, but rather than adopting one assumed form for the underlying distributions, we implicitly integrated the correction over all possible distributions, weighted by the probability for each distribution. To this end, we have integrated over the entire four-dimensional PDF of each mass subsample to extrapolate the binary frequency at separations of 8–5000 AU and spanning all mass ratios of 0 < q < 1. We find that in order of declining mass, our three subsamples (0.5–0.3 M_☉, 0.30–0.15 M_☉, and 0.15–0.07 M_☉) have binary frequencies of 50⁺¹⁰ _{− 9}%, 29⁺⁷ _{− 6}%, and 21⁺⁷ _{− 6}% in this range of parameter space.

5.2.2. Limits on Substellar Multiplicity

As discussed above, we cannot directly constrain the parameters of the substellar binary population because we did not discover any such binaries in our sample. However, we can estimate the probability of a null detection as a function of the four parameters in our model, ruling out a large portion of parameter space. In Figure 7, we show our null detection probability surfaces for the substellar-mass bin in the F-log (s) and F-σ_log (s) planes. We can't integrate over the unplotted dimensions of our PDF since the integral diverges, so we instead show cross-sections for the most likely values as inferred by Burgasser et al. (2006): log (s) ∼ 0.6, σ_log (s) ∼ 0.3 dex, and γ ∼ 4.2. We chose these parameters because the T dwarf sample studied by Burgasser et al. more closely matches our mass range than the full sample of MLT dwarfs studied by A07.

Zoom In Zoom Out Reset image size

We find that for the given values of σ_log (s) and γ, we cannot rule out any mean separations ≲1 AU at >50% confidence. However, we can rule out combinations of increasing mean separation and decreasing binary frequency; if the mean separation is 2 AU, then the binary frequency is <11% at 50% confidence and <38% at 90% confidence. If the mean separation is 4 AU, which is the maximum value consistent with the results of Burgasser et al., then the corresponding frequency limits are <4% and <11%, respectively. Conversely, if the total binary frequency is ∼20%, then the 50% and 90% confidence limits on the mean separation are 1.6 AU and 2.8 AU. The corresponding probability surface for log (s) ∼ 0.6 and γ ∼ 4.2 is strongly concentrated at low frequencies since this mean separation is very close to the inner working angle of our LGSAO survey, and therefore at least half of all companions should have been detectable.

In summary, all of these limits for the substellar regime are extremely discrepant with respect to the confidence intervals for the higher-mass subsamples, which indicates that the mass-dependent tightening of binary systems continues into the substellar regime. There are no well defined and observationally supported models for how low-mass binary systems form, so it is difficult to infer the underlying justification for the continued decline of system separations and/or frequencies into the substellar regime. However, the trend for declining separations and frequencies in the field appears to be established at very early ages and must result directly from the formation process.

5.2.3. Mass Ratio Distributions

In Figure 8, we show the PDFs for our three stellar-mass bins as projected onto the γ-axis; we do not show any results for the substellar-mass bin because our null detection does not yield a useful constraint on its mass ratio distribution. Unlike for Figures 6 and 7, we decided to project the PDF onto an axis instead of a plane in order to display our constraints on the mass ratio distribution. There is no significant covariance between our constraints on γ and those for other parameters, so this choice simplifies our presentation.

Zoom In Zoom Out Reset image size

The 0.3–0.5 M_☉ subsample has a best-fit slope of γ = +0.18^+0.33 _{− 0.30}, a value which is consistent with the linearly flat mass ratio distribution found for higher-mass stars in young associations (Kraus et al. 2008, 2011). By contrast, the 0.15–0.30 M_☉ subsample has a steeper slope of γ = +1.02^+0.59 _{− 0.52}, a value intermediate between the flat slope of higher-mass stars and the typically very steep power laws (γ ∼ 2–4) seen for late-M stars and L/T brown dwarfs in the field.

Finally, the 0.07–0.15 M_☉ subsample has a best-fit slope that is similar to the 0.15–0.30 M_☉ bin, albeit with a very wide confidence interval, yielding γ = 0.96^+0.70 _{− 0.59}. This result at first appears to contradict the overall trend for a steeping mass ratio distribution with declining mass that is seen in the field and would seem to lead from the 0.15–0.30 M_☉ subsample. However, closer inspection of the sample suggests a possible solution. Of the seven binary systems in the 0.07–0.15 M_☉ bin with projected separations of ≲20 AU, all have mass ratios of ≳0.5. By contrast, of the five binary systems with projected separations of ≳25 AU (GG Tau B, CFHT-Tau-17, and CFHT-Tau-18 in Taurus, Hn 13 in Cha-I, and RX J1558.1-2405b in Upper Sco), all but Hn 13 have a mass ratio of ≲0.5. The corresponding limits on the mass ratio power-law exponent are γ = 4.0^+1.9 _{− 1.6} for the closer subset and γ = −0.3^+0.7 _{− 0.5} for the wider subset.

This trend is particularly intriguing because the five wider systems seem to approach or exceed the mass-maximum separation limit observed for field systems by Burgasser et al. (2003), who observed that for VLM binary systems in the field, there is a mass-dependent upper envelope to binary system separations, a_max ∼ 1400 M²_tot. The five systems have typical total masses of ∼0.15 M_☉, corresponding to maximum "allowed" separations of ∼30 AU. As such, they appear to be unusual in both their separation and their mass ratio.

It is tempting to consider whether the markedly different mass function for wider VLM binary systems is a result of a different formation history. For example, wider binary systems most likely form earlier in the collapse of the progenitor molecular core. At these earlier stages, there is still more material left in the circumstellar envelope that might preferentially accrete onto the more massive binary component, driving the mass ratio further from unity. By contrast, close binary systems most likely form in the final stage of collapse, after much of the circumstellar envelope has been accreted into the central mass and little would remain for preferential accretion. Thus, if fragmentation tends to yield similar-mass components, then the epoch of fragmentation would dictate how far the mass ratio could evolve from unity. Since lower-mass binary systems also tend to have smaller separations, this would naturally lead to the trend for lower-mass binaries to have mass ratio distributions which are increasingly peaked at unity.

This model does not explain why these systems do not have analogs in the field, so we must appeal to a separate trend to justify their absence. Most of our sample is drawn from environments that are much less dense than typical star clusters; indeed, all of our targets are in loose associations that are unbound and should disperse within the next 10–50 Myr. By contrast, most stars form in denser clusters (Lada & Lada 2003) that are much more dynamically active and will ionize loosely bound binary systems. For example, the separation distribution for solar-mass binary systems is truncated at separations of ∼300 AU in young clusters like the ONC (Köhler et al. 2006) and at ∼100 AU in older clusters like Praesepe (Patience et al. 2002). VLM binary systems in Praesepe with equivalent binding energy would have separations a factor of ∼3 lower (∼30 AU). Therefore, these systems might have counterparts in denser clusters, but those counterparts could be disrupted into their component singleton stars before reaching the field.

Finally, we must consider a more prosaic explanation as well. All of our targets are located at distances of ∼140 pc, so we cannot resolve binary systems closer than ∼5–10 AU. It is possible that our "unusually wide" binary systems are actually hierarchical multiples, where one component of the wide pair appears fainter (and thus less massive) because it is actually a close double comprised of two stars that each contain approximately half the mass of the primary. This would yield a total mass ratio close to unity in the wide pair, plus the higher total mass would allow for a correspondingly wider separation without violating the a_max–M_tot relation. Large field surveys are starting to uncover a significant number of the very rare systems that appear at first to violate this relation (Caballero 2007; Artigau et al. 2007; Radigan et al. 2009; Dhital et al. 2010), but follow-up high-resolution imaging has shown that at least some of them are hierarchical triples or even quadruples (Law et al. 2010).

Smoothed particle hydrodynamic (SPH) simulations are now capable of producing synthetic populations of young stars with statistically robust sizes (e.g., Bate 2009, 2012; Offner et al. 2010), including samples of synthetic binaries that are equivalent in size to our observed samples of real systems. These synthetic binary populations provide a natural comparison sample for the results of our survey. The only synthetic population large enough to investigate mass-dependent trends was recently reported by Bate (2012), who simulated a cluster of 183 stars and brown dwarfs, including 42 multiple systems. Their simulation only ran for 2 × 10⁵ yr (due to computational expense), and had to end before the cluster dispersed all of its gas or became unbound, so the synthetic population likely was still evolving at the conclusion. However, it still provides an illustrative comparison sample. To produce a direct comparison, we have subjected their binary sample to the same Bayesian formalism as our observed sample, with the caveats that we set no detection limits (since all binaries are known) and we treated the semimajor axis as equivalent to the projected separation (since they are statistically equal to within nearly unity; Dupuy & Liu 2011).

In Figures 6 and 8, we plot the confidence intervals corresponding to the posterior PDF as marginalized to the (F, μ_log (ρ)) and (F, σ_log (ρ)) planes and the γ line. These comparisons show that simulations reproduce a binary population that is broadly similar to that seen in observations; binaries have approximately the correct frequency, semimajor axes, and companion masses. This consistency is especially intriguing because the simulated environment is much denser than the regions we have observed, and hence dynamical interactions shape the synthetic population to a much greater degree. A more detailed comparison shows that not all features are consistent, though.

The simulations reproduce the most distinctive feature of real binary populations: a binary frequency which declines with decreasing primary mass. However, the separation distribution does not show the corresponding trend toward smaller semimajor axes, with a mean separation of ∼20–40 AU across the entire mass range. Furthermore, the standard deviation of the separation distribution only declines in the least massive bin, where it denotes a paucity of binaries at both >100 AU (in agreement with observations) and <5 AU (in sharp contrast with observations). The mass ratio distribution also does not match the expected trend for lower-mass systems to have mass ratios near unity. The lowest-mass bin shows a very strong bias toward such systems, but the highest-mass bin also shows a similar bias.

These discrepancies might be a result of not following the binary formation process to its conclusion. Accretion of additional circumstellar material should further modify the masses of the binary companions, and their semimajor axes also could be modified if the specific angular momentum of the accreted material is much higher or lower than for the binary components. Furthermore, dynamical interactions should disrupt or harden binary systems, and both effects would preferentially bias lower-mass binaries to have smaller semimajor axes than higher-mass binaries. The use of sink particles (with radius r = 0.5 AU) also could affect processes occurring on AU scales, inhibiting the simulated production of tight binary systems.

If longer simulations cannot improve the level of agreement, then it might suggest that other changes are required. The simulations now incorporate radiative feedback (Bate 2009, 2012), and other simulations have suggested that magnetic fields are not likely to be significant on small scales (Price & Bate 2009). One possible avenue is to simulate different initial conditions. The simulations of Bate (2012) currently begin with an isothermal sphere of uniform density, which does not match the configurations seen for pre-stellar environments like the Pipe Nebula (e.g., Lada et al. 2008). Another possible variable to change is the turbulent power spectrum, which might affect the properties of binary systems (e.g., Goodwin et al. 2004).