Close companions around young stars

TL;DR

This study analyzes the multiplicity of young stars using APOGEE-2 spectra, revealing binary properties, distributions, and their relation to stellar environment and disk presence, providing new insights into early stellar evolution.

Contribution

It introduces an autonomous method to identify spectroscopic binaries in a large, homogeneous sample of young stars, expanding understanding of their multiplicity characteristics.

Findings

399 binaries identified out of 5007 sources

Mass ratio distribution is mostly uniform with an excess of twins

No strong trend of multiplicity fraction with cluster age

Abstract

Multiplicity is a fundamental property that is set early during stellar lifetimes, and it is a stringent probe of the physics of star formation. The distribution of close companions around young stars is still poorly constrained by observations. We present an analysis of stellar multiplicity derived from APOGEE-2 spectra obtained in targeted observations of nearby star-forming regions. This is the largest homogeneously observed sample of high-resolution spectra of young stars. We developed an autonomous method to identify double lined spectroscopic binaries (SB2s). Out of 5007 sources spanning the mass range of $\sim$0.05--1.5 \msun, we find 399 binaries, including both RV variables and SB2s. The mass ratio distribution of SB2s is consistent with a uniform for $q<0.95$ with an excess of twins with $q>0.95$. The period distribution is consistent with what has been observed in close binaries ($<10$ AU) in the evolved populations. Three systems are found to have $q\sim$0.1, with a companion located within the brown dwarf desert. There are not any strong trends in the multiplicity fraction (MF) as a function of cluster age from 1 to 100 Myr. There is a weak dependence on stellar density, with companions being most numerous at $\Sigma_*\sim30$ stars/pc$^{-2}$, and decreasing in more diffuse regions. Finally, disk-bearing sources are deficient in SB2s (but not RV variables) by a factor of $\sim$2; this deficit is recovered by the systems without disks. This may indicate a quick dispersal of disk material in short-period equal mass systems that is less effective in binaries with lower $q$.

Tables

Table 1. Table 1: Catalog of the sources observed by APOGEE-2 survey towards the young clusters and star-forming regions

2MASS	$\alpha$	$\delta$	Region	YSO	N	RV	$v_{rot}\sin i$	$T_{eff}$	$\log g$	Veiling	Binary
ID	deg.	deg.		class	epoch	km s-1	km s-1	K	dex		SB2?	${\chi }^2$	Slope	Lit. ${\chi }^2$	Lit. Slope	q
2M03192494+4859402	49.85391998	48.99450302	alphaper	3	4	-0.489 $\pm$ 0.536	20.38 $\pm$ 0.97		5.223 $\pm$ 0.067	0.289 $\pm$ 0.038	1	0.3	0.4	0.3	0.4
2M03193068+4903021	49.87783813	49.05060196	alphaper	3	4	-9.223 $\pm$ 0.176	0.076 $\pm$ 0.56		5.090 $\pm$ 0.054	0.290 $\pm$ 0.025	1	0.4	0.4	0.4	0.4
2M03193117+4941171	49.87989044	49.68811035	alphaper	3	3	0.588 $\pm$ 1.881	37.76 $\pm$ 3.00		5.161 $\pm$ 0.132	0.195 $\pm$ 0.068	1
2M03220975+4834024	50.54065704	48.56735229	alphaper	3	4	15.332 $\pm$ 0.214	14.82 $\pm$ 0.48		5.346 $\pm$ 0.057	0.403 $\pm$ 0.024	2	1641.8	64.5	1641.8	64.5	0.89 $\pm$ 0.02

Table 2. Table 2: Stellar properties measured from the IN-SYNC pipeline across the individual epochs

2MASS	MJD	RV	$v_{rot}\sin i$	$T_{eff}$	$\log g$	Veiling	SNR
ID		km s-1	km s-1	K	dex
2M03173879+4850489	57821	-0.77 $\pm$ 0.35	0.086 $\pm$ 0.896	3617 $\pm$ 72	4.97 $\pm$ 0.14	0.52 $\pm$ 0.06	41.27
2M03173879+4850489	58097	-0.00 $\pm$ 0.39	0.056 $\pm$ 0.621	3464 $\pm$ 63	4.92 $\pm$ 0.15	0.49 $\pm$ 0.07	37.86
2M03173879+4850489	58148	-1.11 $\pm$ 0.38	0.148 $\pm$ 1.802	3507 $\pm$ 50	5.09 $\pm$ 0.13	0.41 $\pm$ 0.06	35.81
2M03173879+4850489	58179	-1.29 $\pm$ 0.61	0.208 $\pm$ 2.319	3448 $\pm$ 75	5.03 $\pm$ 0.22	0.42 $\pm$ 0.10	24.25

Table 3. Table 3: Observed clusters and star-forming regions

Cluster	Ageaaa - Luhman et al. (2016); b - Luhman (2018); c - Venuti et al. (2018); d - Kounkel et al. (2018); e - Balachandran et al. (2011); f - Bouvier et al. (2018)	N∗	Nbin	$T_{eff}$,5-50-95%	$\log g$,5-50-95%	$\Delta$t	Nvisits	RV	${\mu }_{\alpha }$	${\mu }_{\delta }$	$\pi$
	Myr			(K)	(dex)	(days)		km s-1	mas yr-1	mas yr-1	mas
IC348	3 (a)	299	24	2750–3800–4990	3.9–4.5–4.8	538	4	10–25	0–10	-15–0	1.5–5
NGC 1333	1 (a)	132	7	2680–3980–5060	3.2–4.3–5.1	438	4	5–25	5–20	-20–0	2–6
Taurus	2 (b)	199	37	2550–3770–5220	3.4–4.4–5.0	357	3	10–30	-5–20	-35– -10	4.8–10
NGC 2264	3 (c)	265	28	3210–4190–5010	3.6–4.6–4.9	1892	6	10–35	-8–2	-10–2	0.8–2.5
ONC	2 (d)	1274	137	2990–3980–4690	3.5–4.4–4.8	1457	2	14–36	-5–5	-5–5	1.5–3.5
L1641	2 (d)	711	40	2990–4000–4780	3.8–4.5–4.8	1476	3	14–36	-5–5	-5–5	1.5–3.5
Orion B	1 (d)	308	23	3010–4170–5310	3.5–4.3–4.7	605	3	14–36	-5–5	-5–5	1.5–3.5
Orion CbbFollowing nomenclature from (Kounkel et al., 2018); sources without membership assignment were split on the basis of their RVs.	4 (d)	360	34	3060–4210–5540	3.9–4.6–5.0	628	3	14–25	-5–5	-5–5	1.5–3.5
Orion DbbFollowing nomenclature from (Kounkel et al., 2018); sources without membership assignment were split on the basis of their RVs.	6 (d)	473	19	3324–4280–5710	4.4–4.7–5.0	628	2	25–36	-5–5	-5–5	1.5–3.5
$\lambda$ Ori	4 (d)	379	33	3240–4130–5250	4.1–4.6–4.9	710	4	14–36	-5–5	-5–5	1.5–3.5
$\alpha$ Per	50 (e)	152	9	3020–3900–5700	4.7–5.1–5.5	358	4	-15–15	15–35	-35– -20	4–8
Pleiades	125 (f)	455	36	3060–4210–5770	4.8–5.2–5.5	2190	3	0–20	10–30	-55– -35	6–9

Table 4. Table 4: Deconvolved CCFs

2MASS	MJD	FWHM1	AMP1	$v_1$	Flag1
ID		km s-1		km s-1
2M03173879+4850489	57821	20.1$\pm$ 0.9	0.77$\pm$0.03	-2.11$\pm$ 0.36	4
2M03173879+4850489	58097	17.9$\pm$ 1.1	0.63$\pm$0.03	-0.16$\pm$ 0.46	4
2M03173879+4850489	58148	18.2$\pm$ 1.4	0.69$\pm$0.05	-2.24$\pm$ 0.60	4
2M03173879+4850489	58179	25.9$\pm$ 0.7	0.50$\pm$0.01	-4.35$\pm$ 0.31	4

Table 5. Table 5: Flag quality

Flag	Fraction of	Summary
	multiply deconvolved spectra
1	42.4%	Failed FWHM/amplitude test; could be noise
2	34.8%	Fail the symmetry test; could be falsely multiply deconvolved
3	11.4%	Inconclusive SB2/Spotted star pair
4	11.4%	The primary peaks in all stars, and the secondary peaks of bona fide SB2s

Table 6. Table 6: V m ​ a ​ x subscript 𝑉 𝑚 𝑎 𝑥 V_{max}  and S m ​ a ​ x subscript 𝑆 𝑚 𝑎 𝑥 S_{max}  cuts to identify close binaries

Period	$V_{max}$	$S_{max}$
(day)	km s-1	km s-1
2000	2.6	21
1000	4.6	26
500	7.7	33
200	13	46
100	20	62
50	27	78
20	40	110
10	60	140
5	90	160

Equations

$$\sum _{k=1}^N{\mathrm{AMP}}_k\exp \left[-4\ln 2(x-v_k{)}^2/{\mathrm{FWHM}}_k^2\right]$$

$$e_{max}(P)=1-{\left(P/k\right)}^{-2/3}$$

$$\begin{gathered} MF(<10AU)=16.5\times (308.9-119.6H+ \\ 17.46H^2-1.131H^3+0.02735H^4)\times \\ (7.634-2.965T_{eff}\times 10^{-3}+3.189T_{eff}^2\times 10^{-7}) \end{gathered}$$

Full text

binaries: spectroscopic, stars: pre-main sequence

††software: IN-SYNC pipeline (Cottaar et al., 2014), GaussPy (Lindner et al., 2015), PHOENIX spectral library (Husser et al., 2013), PARSEC-COLIBRI isochrones (Marigo et al., 2017), TOPCAT (Taylor, 2005)

1 Introduction

Approximately half of the solar-type stars in the solar neighborhood are found in multiple stellar systems (Duquennoy & Mayor, 1991; Raghavan et al., 2010). The period distribution of companions in these systems is log-normal, and the multiplicity fraction (MF) of close companions with separations $<$10 AU is $\sim$20% (Moe & Di Stefano, 2017). The total MF is mass dependent, however, decreasing to $\sim$30% for M dwarfs (e.g,. Ward-Duong et al., 2015), and increasing to almost 100% for OB stars (Sana et al., 2014). In addition to this dependence on primary mass, the MF has been found to depend on metallicity, in that metal poor solar-type stars are more likely to have a close companion (Moe et al., 2018; Badenes et al., 2018). A recent review of multiplicity has been conducted by Duchêne & Kraus (2013).

Binary systems form early in a star’s lifetime; companions are commonly observed around the youngest protostars (e.g., Chen & Arce, 2010; Tobin et al., 2016). The primordial multiplicity distribution may be altered through dynamical interactions between the young stars, such as through Kozai-Lidov interactions (e.g., Moe & Kratter, 2018), but by the time these stars disperse from their birth environments, their multiplicity properties are expected to resemble, on average, those found in the field population.

A number of studies of multiplicity in nearby star-forming regions have been conducted, typically through high resolution imaging. The current diffraction limited observations, however, struggle to resolve companions closer than 10 AU even in the closer ($<$200 pc) regions. However, a number of differences have been observed between the MF measured in these star-forming regions and the field. Most notably, the Taurus Molecular clouds are known for having a very high MF, in excess of 60–70% (e.g., Kraus et al., 2011), significantly higher than $\sim$40% found in the field (Raghavan et al., 2010). A similar excess has also been found recently in the Orion Nebular Cluster (ONC), where companions with separations of 10–100 AU are twice as numerous as the rate that is found in the field (Duchêne et al., 2018). At larger separations, however, the MF in the ONC drops to the levels consistent with the measurements in the field, possibly due to earlier dynamical interactions with other cluster members (Reipurth et al., 2007). In the larger Orion Molecular Clouds, however, outside of the ONC, the companions with separations of 100–1000 AU have been found to be $\sim$2 times more common in the densely populated stellar regions than in the more diffuse regions (Kounkel et al., 2016b), suggesting a density dependent mechanism for the fragmentation of protostellar envelopes.

In order to analyze the distribution of companions at the separations closer than 10 AU, spectroscopic studies have been conducted to search for radial velocity (RV) variable systems (e.g., Mathieu, 1992; Melo, 2003; Maxted et al., 2008; Tobin et al., 2009, 2015; Kuruwita et al., 2018). Most of these studies, however, have focused on individual regions, and comparing their results is difficult due to differences in temporal coverage. MF measured by these studies was broadly consistent with what is observed in the field, but small sample sizes make a more robust analysis difficult. A direct comparison has been made between the Chameleon I and Taurus (Nguyen et al., 2012) regions, and between the ONC and NGC 2264 (Kounkel et al., 2016a), which have found MFs that are largely consistent. On the other hand, Jaehnig et al. (2017) analyzed the APOGEE IN-SYNC observations of Orion, Perseus, and the Pleiades, finding an evidence for a possible evolution of the MF with age.

Fernandez et al. (2017) conducted a first large study of the double lined spectroscopic binaries (SB2s) in the young clusters in the IN-SYNC data. They identified 104 potential systems, doubling the number of known young SB2s, although their approach was largely reliant on the visual examination of the data, and thus it was difficult to perform any robust statistics on the sample.

The presence of multiples does affect the early stages of stellar evolution because of the effect binaries have on protoplanetary disks. Companions at separations closer than 200 AU are more likely to form from disks rather than core fragmentation (Tobin et al., 2018; Moe et al., 2018). Disks around closely-separated binaries tend to be less luminous at sub-mm wavelengths (e.g., Jensen et al., 1996; Harris et al., 2012; Akeson et al., 2019) than around single stars. Disks in multiple systems also have shorter lifetimes than disks around single stars (e.g., Cieza et al., 2009; Kraus et al., 2011), indicating that companions dynamically disrupt both the circumstellar and circumbinary disks. However, a number of close circumbinary planets have been discovered (e.g., Doyle et al., 2011; Orosz et al., 2012), indicating some circumbinary disks must survive long enough to form planets.

In this paper, we present a study of the MF as a function of age, mass, separation, environment, and the presence of a protoplanetary disk. We base this study on the analysis of the APOGEE observations of nearby star forming regions and young clusters. We search this dataset for RV variables111The term SB1 is not used to avoid the confusion, because some of the sources identified as RV variable are also found to be SB2s and SB2s with separations up to 10 AU, and report here the largest uniformly observed sample of close companions around young stars. In Section 2 we discuss the APOGEE observations, data processing, and the sample construction. In Section 3 we identify RV variables, and in Section 4 we discuss the methods for an automated identification of SB2s. In Section 5 we construct a sample of synthetic binaries from which we infer the completeness of our search methods as a function of each system’s orbital parameters. In Section 6 our measured and synthetic samples to test for differences between the multiple properties of young stars and those in the field. Finally, in Section 7 we discuss the implications of the observations on the formation and evolution of young close multiples.

2 Sample definition

2.1 APOGEE Observations

High resolution near-infrared spectra of several nearby star-forming regions and young clusters have been obtained by the Apache Point Observatory Galactic Evolution Experiment (APOGEE) spectrograph, which is mounted on the 2.5 m Sloan Digital Sky Survey (SDSS) telescope (Gunn et al., 2006; Blanton et al., 2017). The APOGEE spectrographs can observe up to 300 sources simultaneously across a 1.5∘ field of view, and the instrument covers the spectral range of 1.51–1.7 $\mu$m with the typical R$\sim$22,500 (although it may vary between the individual fibers, Wilson et al., 2010; Majewski et al., 2017). The typical limiting magnitude for the observations analyzed here is H$\sim$13 mag.

The analysis in this work focuses on YSOs in the Orion Complex (Da Rio et al., 2016, 2017; Kounkel et al., 2018), NGC 1333 (Foster et al., 2015), IC 348 (Cottaar et al., 2015), Taurus Molecular Clouds, NGC 2264, $\alpha$ Per, and Pleiades. The observations originally began as part of the SDSS-III IN-SYNC survey with APOGEE, and later transitioned into the SDSS-IV Young Clusters Survey with APOGEE-2 with an expanded list of regions, increased footprint on the sky, and additional epochs covering a larger baseline for some of the existing targets. Sources were selected for targeting using the presence of infrared excess, optical variability, and when feasible, identification as a YSO in the literature (Cottle et al., 2018). The goals of the APOGEE-2 survey was to create a homogeneously selected sample that is representative of the underlying population, although it may not be complete. Sources cannot be observed simultaneously if they are within 72” of each other, otherwise a fiber collision would occur. In densely concentrated regions, this required a strategy of using multiple plates to cover the same region in the sky but containing different sets of sources. When it was impossible to fill the targets with just the candidate YSOs, the plate included other field sources.

In total, 53,452 spectra were obtained for 14,823 stars across $\sim$150 sq. deg. field of view over 6.5 years, from September 2011 to March 2018 (Table 1). A third of the data was made public as part of SDSS DR14 (Abolfathi et al., 2018); the remaining spectra will be made released in DR16.

2.2 Stellar Parameters

The data were processed and stellar parameters (RV, $T_{eff}$, $\log g$, $v_{rot}\sin i$, and veiling - continuum excess in the spectrum due to accretion) were extracted from all the individual spectra using the pipeline developed as part of the IN-SYNC survey (Cottaar et al., 2014), and an uncertainty weighted average was computed for each of the parameters measured for each star over all visits. These parameters are separate from those computed by the APOGEE’s primary ASPCAP pipeline (García Pérez et al., 2016), as ASPCAP is primarily calibrated for the red giants, and is not as successful in extracting accurate stellar parameters for dwarfs and pre-main sequence stars. RVs are largely consistent between both the IN-SYNC and the ASPCAP pipelines, however. There are some systematics in the $T_{eff}$values measured with the IN-SYNC pipeline, most notable of which occurs at $\sim$3500 K, with this temperature being strongly disfavored, resulting in a gap in the overal distribution, which is apparent in Figure 1. A full description of the data processing methods applied to these spectra is presented in Kounkel et al. (2018). All of the measurements of the stellar properties for all the individual epochs for all sources across the footprint of the survey are included in Table 2.

To aid our identification of multiple systems, we supplement the stellar parameters extracted with the IN-SYNC pipeline with an analysis of the a cross correlation analysis. A cross correlation function (CCF) was computed for each spectrum by correlating the normalized spectra with normalized template from the PHOENIX spectral library (Husser et al., 2013) that best matched in $T_{eff}$ and $\log g$ to the source. The template was broadened by a 10 km s*-1* rotation kernel to achieve a smooth profile as 10 km s*-1* is a typical $v_{rot}\sin i$ for the pre-main sequence stars in the sample. The sky lines and the known bad pixels were masked. And, as the spectra observed by APOGEE are split into three separate chips with a chip gap spanning $\sim$50 Å, the CCF was computed for all the chips separately and added together. Correlation was done across 401 velocity channels, separated by 1 km s*-1*, centered at the velocity measured by the IN-SYNC pipeline at the given epoch rounded to the nearest integer km s*-1*.

Not all of the spectra obtained in the fields covered by the program were YSOs. If excess fibers were available in a given field, they were usually assigned to field red giants, which comprise the main targets of the APOGEE-2 survey. Additionally, the preliminary analysis of the yield of the young cluster targeting methods quantified a contamination rate of 5–10% (Cottle et al., 2018) due to inclusion of nearby dwarfs or distant giants.

2.3 Membership identification

To curate the sample limited only to the members of the corresponding star forming regions and young clusters, several cuts were applied. The catalog was cross-matched with Gaia DR2 (Gaia Collaboration et al., 2018; Lindegren et al., 2018), requiring agreements in positions of better than 1”. For sources that had Gaia astrometric solutions (95% of the total sample), we required the source to have weighted average RV or a combination of parallax and proper motions to be consistent with the typical cluster properties to within 15–25 km s*-1*  and parallax to within several mas, depending on the cluster, to fully encompass the Gaussian distribution of positions and velocities within a given region (Table 3).

Because Gaia DR2 does not have a prescription for the astrometric binaries (which on average would have longer periods than the spectroscopic binaries), the astrometric solutions for some of them may fall outside of these accepted thresholds. High extinction that is frequently associated with YSOs may further degrade the astrometry, even for single stars. To avoid unduly biasing our sample in the presence of these effects we did not impose any quality cut on the Gaia measurements, neither in terms of the uncertainties, nor in terms of the excess astrometric noise, or any other alternative metric. On the other hand, spectroscopic binaries that are not observed with the full phase coverage may appear to have an average RV that is different from the true average by up to several dozen km s*-1*. By requiring a consistency in either (and not necessarily both) of those properties, we attempt to retain all the possible YSO multiples in the sample. Sources that were not part of Gaia DR2 (5% of the sample) are most likely too reddened to have an optical detection. They are included without any RV constraint imposed onto these objects.

While these loose cuts may allow for an inclusion of some field stars as part of the cluster membership, this is statistically less detrimental to the analysis in this work than the exclusion of possible multiples. Any field stars will have a field probability of being a binary with the field orbital parameters, and their number in the sample is highly unlikely to be large enough to significantly alter the measured distribution in each cluster. On the other hand, the sample of spectroscopic binary YSOs is significantly smaller, so that any accidental exclusion of some of them would have a much stronger effect.

To provide the most inclusive sample of cluster members, we also examined the location of the sources in the $T_{eff}$/$\log g$ and $[G_{B}-G_{R}]/M_{G}$ diagrams were examined (Figure 1,2). Low mass YSOs are very distinct from more evolved stars in these diagrams. Even if sources failed both the distance and the RV cuts, they were included in our analysis if they fall into the appropriate portion of these diagrams. The reason why it was done is because, in some cases, e.g., hierarchical triple systems, it is possible for multiplicity to affect both RV and astrometric measurements.

Another constraint on the curated sample was to impose a T_{eff}$$<6000 K limit. Sources hotter than 6000 K have spectra dominated by hydrogen lines, resulting in a broad CCF and more uncertain RV measurements, making them poorly suited for a uniform analysis of spectroscopic multiplicity. Moreover, it becomes significantly more difficult to distinguish hotter YSOs from the more evolved field stars from their position in the HR diagram. Finally, we required SNR of the spectra to be $>20$, and that CCF could be decomposed into at least a single Gaussian component (Section 4).

In total, the curated sample contains 5007 stars/19127 individual spectra. Characteristic properties of stars in each cluster are given in Table 3: the number of stars in each individual region, the 5th and 95th percentile of $T_{eff}$ and $\log g$  baseline spanned by APOGEE, and a typical number of epochs per star. The cumulative distribution of the parameters is shown in the Figure 3.

2.4 Disk classification

To investigate the relationship between multiplicity and protoplanetary disks, the evolutionary classification on the state of the protoplanetary disk was obtained from the works of Megeath et al. (2012) for Orion A and B, Hernández et al. (2007) for $\sigma$ Ori, Suárez et al. (2017) for 25 Ori, Hernández et al. (2010) for $\lambda$ Ori, Rebull et al. (2010) for Taurus, Rapson et al. (2014) for NGC 2264, and Luhman et al. (2016) for NGC 1333 and IC348. These catalogs were further supplemented by the WISE classification from Marton et al. (2016), and the remaining sources were assumed to be Class IIIs.

In the sample, there are 1882 disk bearing Class II sources, and 3125 are diskless Class III sources. There does not appear to be any biases in targeting sources of either evolutionary class; the Class II to Class III ratio among the sources in the curated sample is representative of what is typical in the corresponding clusters. It should be noted that there are a few Class I YSOs in the sample, but not enough to do any robust statistics on them separately. As they also have protoplanetary disks, they are grouped together with Class II systems.

In general, spectra of Class II and Class III systems have similar signal to noise ratio, similar fluxes, and similar shape of their CCF profiles. The only quantitative difference between them is the veiling, which has been measured by the pipeline. As mentioned in Kounkel et al. (2018), the measurement of veiling in the spectrum has a number of systematic effects as a function of $T_{eff}$, because the pipeline is subject to degeneracy between veiling and metallicity (Figure 4). In general, measurements of veiling $<$0.6 may not necessarily be reliable. However, with only a few exceptions, Class II systems do have systematically higher veiling than Class III systems.

3 RV variables

In this section we analyze the results of the IN-SYNC pipeline, identify and correct the systematic offsets in the measured RVs, and identify RV variables.

3.1 Identification of RV variables

Identification of RV variables was performed by applying two criteria to the IN-SYNC pipeline extracted RVs. The first criterion is a reduced $\chi^{2}>16$, which identifies sources that have a significant scatter in the typical velocity of a star. This is the same criterion that was used in Kounkel et al. (2016a), and it translates to a 4$\sigma$ detection of RV variability. The second criterion is derived from the shape of the linear fit to all the available RV measurements for each system as a function of time. We flag sources as RV variable if the slope of their fitted RVs is inconsistent with 0 by $>4\sigma$. This second criterion helps to identify systems with long periods and that have a larger number of observations; these systems may be missed just with a $\chi^{2}$ cut. For both of these criteria, we required a minimum of 3 epochs, as the determination from just two epochs is both incomplete and more likely to result in false positives (Section 5). However, even with these criteria, it is still possible that some of the systems identified may be false positives, both due to the random counting statistics as well as some effects outside of multiplicity that may affect RVs (such as starspots, discussed more in Section 4.2

As was noted in Kounkel et al. (2018), some of the RVs measured from lower SNR spectra by the IN-SYNC pipeline are not necessarily representative of the true velocity of the star. The spurious RVs are usually shifted by $>50$ km s*-1*, making the source appear to be strongly RV variable despite the stability of the remaining measurements. While SNR$>$20 cut eliminates most of these spurious RVs, a handful of these poor measurements do occur in spectra with SNR$>$20 and thus persist in our sample. To confirm the accuracy of all the measurements in the curated sample, we visually examined the CCFs of all the sources identified as RV variable, as well as all the sources that had pipeline extracted RVs inconsistent with the position of the primary Gaussian component by more than 5 km s*-1*. In total, 20 out of $>$19K observations in our sample ($\sim$0.1%) were affected by this issue, and their RVs were corrected to the RV determined from the CCF that we compute independently of the IN-SYNC pipeline. This check was not performed for the sources that have failed the membership criteria.

3.2 Systematics

Afterwards, a consistency check was performed to find an epoch specific systematic offset across the entire observing field. We calculated a difference between all the individual RV measurements (at a given epoch for a given field) and the weighted averages of the sources that have at least 3 epochs and have not been identified as RV variables. If an epoch had more than 50 such sources, the median offset was calculated for the entire epoch. There are 144 (out of 171) epochs for which this could be performed; 75% of them have a systematic offset $<$0.1 km s*-1*, and for 92% the offset is $<$0.2 km s*-1*, which is smaller than the typical $\sigma_{RV}$ as measured by the IN-SYNC pipeline (which is typically $\sim$0.2–0.3 km s*-1*). For only a few epochs the systematic offset is larger than the typical uncertainties of the observations. These offsets are removed from the RVs.

As was noted by Cottaar et al. (2014), RVs of the sources with low $T_{eff}$ appear to have a strong systematic offset. For each star forming region, we average the RVs of all the sources that are hotter than 3500 K (excluding those that are identified as multiples) to get the region typical RV. In Figure 5, we compare the RV offset of all sources relative to their cluster mean as a function of $T_{eff}$. Sources cooler than 3400 K appear to be systematically redshifted by $\Delta$RV=12.84-0.0038\times$$T_{eff}, which results in an offset as large as 4 km s*-1* at T_{eff}$$\sim2400 K. This effect is present regardless of the RV extraction method, as the RVs extracted by the IN-SYNC pipeline are generally consistent with RVs inferred from the CCFs, as well as with those reported in the APOGEE ASPCAP catalog. Most likely this is due to the synthetic spectra not carrying a sufficient precision in the energy levels for the lower $T_{eff}$ spectra.

3.3 Literature RVs

To improve the temporal baseline coverage of the data, we cross-match the sources in the curated sample with optical RV surveys, namely Dolan & Mathieu (2001), Fűrész et al. (2006, 2008), Briceño et al. (2007), Flaherty & Muzerolle (2008), González Hernández et al. (2008), Maxted et al. (2008), Sacco et al. (2008), Nguyen et al. (2012), Hernández et al. (2014), Kounkel et al. (2016a, 2017), and Jackson et al. (2018). However, this creates a very inhomogeneous dataset, as some of these surveys do not report individual RV measurements and instead average several measurements together. Many of them also do not list the dates of the individual observations, or even the related uncertainties. We use these data only for the purpose of identifying candidate RV variables in Lit $\chi^{2}$ and Lit Slope columns of Table 1, and do not rely on them for our further analysis of the MF or sample completeness in each region.

The RVs that are listed in Table 2 are corrected to remove all of the aforementioned systematics (i.e, effect of temperature, and epoch dependent offsets). In total, 123 sources can be identified as RV variables in the APOGEE curated sample on the basis of the $\chi^{2}$, 161 from the slope, and 174 from either method, out of 2774 sources in the curated sample that have at least 3 epochs. Including the archival data, there are 205 sources that can be identified as RV variables from the $\chi^{2}$ cut, 255 from the slope (only including those observations with a well-defined date), and 297 in total, out of 3352 sources with at least 3 epochs (Table 1).

4 SB2

In this section we develop and test a metric for an autonomous identification of SB2s and discuss some of the limiting cases that may confuse this metric.

4.1 Initial decomposition

Identification of SB2s was done on the shape of the CCF. To characterize the coherent signal in the CCF, and to relate it to the physical stellar properties, we also computed CCFs of the synthetic spectra of synthetic binary systems. The preliminary sample was constructed only for the purposes of calibrating and testing the pipeline, as well as exploring the initial effects on the parameter space, and it is not astrophysically significant beyond this; a more detailed synthetic sample with the full analysis of completeness is later discussed in Section 5.

The synthetic binaries were constructed using the PHOENIX spectral library (Husser et al., 2013), combining two templates with $T_{eff}$ and $\log g$ set by interpolating two randomly selected stellar masses onto a 5 Myr PARSEC-COLIBRI isochrone (Marigo et al., 2017). Both spectra were convolved with a velocity kernel corresponding to a random $v_{rot}\sin i$(same for both stars), offset from each other by a velocity that corresponded to random separations, modulated by a random factor that accounts for the orientation of the system and eccentricity. Then the spectra were added together, flux scaled to represent a random distance, and the combined spectrum was interpolated over the wavelength range of a typical APOGEE spectrum. Random flux noise that is typical of the empirical APOGEE spectra was applied. In total, we computed 90,000 synthetic spectra of binary stars, as well as 10,000 spectra of the control single stars.

CCFs were deconvolved into Gaussians using a Python package GaussPy (Lindner et al., 2015). The package relies on the Autonomous Gaussian Deconvolution which fits a profile described by a sum of $k$ Gaussians:

[TABLE]

where $\mathrm{AMP}_{k}$, $v_{k}$, and $\mathrm{FWHM}_{k}$ describe the amplitude, mean position, and width of each of the Gaussian components.

The initial estimates of all the parameters, as well as the number of components ($N$) that are being fitted are chosen by the procedure through derivative spectroscopy. The only parameter that is supplied to the procedure is $\alpha$ which filters the derivative of the function and controls the relative balance between real variance and noise. GaussPy may struggle to find a solution if the fitted function does not level out to the baseline of zero. Therefore, to improve GaussPy’s performance and to limit the number of peaks that would be fitted, from each CCF we subtracted the median, or 20% of the peak, whichever was largest. The CCF was also extended by 100 additional points set to zero along each end. GaussPy was trained on the CCFs produced by the synthetic spectra, optimizing the recovery of Gaussians with $v$ and FWHMs to accurately reflect the injected radial and rotational velocities. Based on this training exercise, we adopt $\log\alpha=$1.5 for the analysis, as it provided the optimal balance of filtering out the noise structure from the CCF and fitting the primary peaks.

There may be other techniques that could be used to decompose CCF into Gaussians. For example, Finite Mixture Models such as ‘mclust’ rely on well-established statistical methods based on maximum likelihood estimation (McLachlan & Peel, 2000; Fraley & Raftery, 2002). Performance of the GaussPy has not yet been tested against such techniques. Nonetheless, it was chosen due to the ease of implementation.

We report all of the GaussPy fitted parameters in the Table 4. After the deconvolution, three cuts (described below) were applied to filter the sample. The reliability of the identified components is encoded in numerical quality flags on a scale from 1 to 4, depending on the success of passing these filters, with the quality flag of 4 being the most reliable 5.

First, Gaussians with an amplitude $<0.15$ or $>3$, as well as a FWHM$<1$ or $>500$ km s*-1* were rejected. This cut removes the spurious CCF fits that do not correspond to realistic stellar properties. The components that failed these cuts were assigned a quality flag of 1. In the total sample, out of the initial sample of 3804 spectra that had been identified to have multiple components, only 2192 satisfied this cut (1525 of 2298 spectra in the curated sample).

A second cut removed sources whose Gaussian fits erroneously separated a single smooth peak of the CCF into multiple Gaussians. These spurious fits were identified by computing the CCFs’s antisymmetric noise ($\sigma_{a}$ in $R=h/\sqrt{2}\sigma_{a}$ from Tonry & Davis, 1979) within a 30 km s*-1* radius ($\sigma_{30}$) centered on the Gaussian with the largest amplitude. Figure 6a shows $\sigma_{30}$ as a function of the velocity separation $\Delta v$ between the multiple Gaussians recovered from a single spectrum, including both synthetic binaries and a sample of synthetic control single stars. We use the location of multiple components extracted from CCFs for control stars to identify areas of parameter space in Figure 6a where spurious extraction reside. Eliminating spurious extractions from empirical CCFs in this region as well, we discarded secondary components if the difference in the velocity of the multiple peaks $\Delta v$ is smaller than the FWHM of the primary peak and $\log_{10}(\Delta$RV$)<-0.32258\log_{10}(\sigma_{30})+1.07419$. While there are some CCFs that do appear asymmetric upon visual examination that would fail this cut (e.g. bottom right panel in Figure 7), this limited the false positive identification of multiple significant components extracted from synthetic single star CCFs to 0.65$\pm$0.08%. The secondary components that failed this cut were assigned the quality flag of 2. Within the footprint of young clusters, 728 sources/1325 visits were deconvolved into multiple components that pass both cuts, of which 573 sources/1022 visits were part of the curated sample.

4.2 Spotted star analysis

An asymmetric CCF that can be characterized by a multi-Gaussian fit may not necessarily indicate the presence of a binary. Young stars with convective envelopes do have a strong degree of magnetic activity, and significant portions of their surface are often covered by spots (e.g., Bouvier & Bertout, 1989; Scholz et al., 2005). A spot’s apparent velocity will depend both on its location, and the star’s rotation period. In many cases, the spot will have a radial velocity that differs significantly from the stellar systemic velocity, producing a line/CCF profile that resembles an SB2 (e.g. similar to the right panel of Figure 7).

In order to investigate the effect that spots would have on the CCF, we generated synthetic spectra of artificially spotted stars in a manner similar to the spectra of artificial binaries. Instead of broadening a synthetic spectrum by a kernel that would correspond to a given $v_{rot}\sin i$ as was done previously, a stellar disk was spatially split into 500$\times\log$ $v_{rot}\sin i$ evenly spaced regions. The spectrum from each region was then shifted to the apparent velocity due to the projected rotation velocity, and all spectra were coadded. For non-spotted stars, this simply approximates convolution with a broadening kernel. Then, a number of spots (up to 10) were generated, covering a total area up to 50% of the stellar disk. These spots are centered at randomly chosen regions, and all of the regions covered by a spot are set to the temperature of the spot ($T_{spot}$). $T_{spot}$ is chosen for each $T_{eff}$ according to the $T_{spot}$/$T_{eff}$ relation given by Berdyugina (2005). A total of 10,000 spotted spectra of single stars were generated, and their CCFs were analyzed using the same procedure as described above.

In this synthetic sample, 8.2$\pm$0.3% of the CCFs were deconvolved to have multiple components. While many of the CCFs produced for the synthetic spotted stars had a high degree of high frequency structure that can only be indicative of spots, there were also a significant number of clean CCFs with a slight asymmetry that, upon blind examination, cannot be distinguished from the CCFs of SB2s with low $\Delta v$. The recovery rate and CCF shape did not depend strongly on the number of spots injected, the total area covered by spots, or on the spot’s latitudinal or longitudinal positions (which was tracked only in a single spot systems). However, spots cannot produce multiple peaks in the CCF whose $\Delta v$ is significantly larger than the star’s $v_{rot}\sin i$. Therefore, a third cut was introduced, flagging secondary components with FWHM${}_{1}>0.8\Delta$RV km s*-1* or FWHM${}_{1}>0.25\Delta$RV+27.5 km s*-1* as likely spot signatures, this removed 97.8$\pm$0.5% of the synthetic spotted sample that remained from the previous cut. Potential components that trigger this likely spot cut were assigned a quality flag of 3, and the remaining binary candidates that pass all the previous tests were assigned a flag of 4. It should be noted that a number of bona fide SB2s likely fail this cut: in the spotless synthetic sample, 22% of the previously identified binaries are rejected. Further monitoring would be required to separate spotted stars and true binaries more conclusively (by detecting the orbital period, strong variability in the shape of the CCF, or through identifying spots/eclipses in the stellar light curve). This cut may not be necessary for the analysis of the stars on the main sequence. Nonetheless, we restricted our sample of SB2s to only those sources for which we can conclusively discard the possibility of the second peak originating due to star spots, i.e., the sources that had multiple components with the quality flag of 4 in at least one epoch.

4.3 Final sample and stellar properties of SB2s

We identify 229 sources/432 visits as SB2s in the full sample (14,823 sources), of which 141 sources/268 visits are in the curated sample of 5007 sources (Figure 7). Only 20% of SB2s can also be identified as RV variables (32% of the sources with at least 3 epochs). Comparing these results to the sample of 104 SB2s from Fernandez et al. (2017) who used the earlier APOGEE/IN-SYNC sample in their analysis, we identify 65 of them as SB2s with a flag of 4, 24 with a flag of 3, 1 with a flag of 2, and 6 with a flag of 1. Only 8 systems are not recovered at all, usually due to the secondary component having too low of an amplitude for GaussPy to recover.

We compare the RV extracted using GaussPy to those RVs that have been extracted by the IN-SYNC pipeline. If the ratio of the amplitudes of the two peaks is AMP2/AMP${}_{1}>0.6$ then typically the pipeline RV was positioned in between two components (Figure 8, left panel). This occurs regardless of the velocity separation between the peaks. Conversely, in the high flux ratio regime (AMP2/AMP${}_{1}<0.6$) the pipeline RV agrees with the strongest component.

The $v_{rot}\sin i$ values determined by the IN-SYNC pipeline can also be used to identify SB2s. The IN-SYNC $v_{rot}\sin i$ for SB2s are often inflated and unstable from one epoch to the next (Figure 8). By contrast, the FWHMs extracted from SB2 CCFs accurately reflect the actual $v_{rot}\sin i$ used to construct the synthetic sample (Section 5). In the data, there is also a good agreement between the measured $v_{rot}\sin i$ and FWHM in the sources not identified as SB2s. In the empirical spectra previously flagged as SB2s, however, the IN-SYNC pipeline $v_{rot}\sin i$ is typically larger by a factor of $\sim 2$ than that implied by the FWHM of the CCF peak, because the IN-SYNC pipeline is unable to resolve two components and conflates both of them into one. The IN-SYNC pipeline measures a $v_{rot}\sin i$ value for each epoch individually, and a weighted average is then computed for each system. For stars that have not been flagged as SB2s, individual measurements are usually consistent with each other, within the errors, but for the SB2s the scatter is typically $>$10 km s*-1* (Figure 8, right panel). This is the case for almost 80% of the identified SB2s; conversely, only 25% of the systems that do have \Delta$$v_{rot}\sin i$$>10 km s*-1* are not identified as SB2s by previous cuts. The reason for this strong variability is because SB2 systems tend to have short periods. Thus the components corresponding to each star in the binary would have different RV in different observations, and the width Gaussian envelope that would conflate both components would vary between the epochs.

It is difficult to determine the reliability of the IN-SYNC pipeline estimates for other stellar parameters such as $\log g$, $T_{eff}$, or veiling for the sources identified as SB2s, but the ranges of these parameters appear to be representative of what is typically found in the individual clusters. The sources identified as binaries (both RV variables, but especially SB2s) do appear to be brighter, however. The location of the binary sequence is ill-defined for pre-main sequence stars because stars of similar masses have a large range of stellar luminosities due to rapid evolution along pre-main sequence tracks. But, spectroscopic binaries do appear to be systematically ‘younger’ on average on the HR diagram compared to the sources not flagged for multiplicity because of the excess flux from the companions (Figure 2). Individual binaries, however, are not guaranteed to have higher fluxes compared to the rest of the population, as this is a function of the flux ratio between the secondary and the primary. Similarly, not all of the sources that lie on the binary sequence will be identified as spectroscopic binaries.

A few of the systems we identify as SB2s and RV variable have been previously resolve through high resolution imaging. Most notable of these systems is 2M04132722+2816247 (=V1096 Tau), for which both stars have been detected with radio interferometry (Galli et al., 2018). Other systems include 2M04184061+2819155 (=V892 Tau, Monnier et al., 2008) and 2M04352089+2254242 (=FF Tau, Kraus et al., 2008).

5 Completeness

To analyze the dependence of the recovered stars on the configuration of the system, we constructed a set of synthetic binary stars with stellar properties matched to those of the curated sample, and multiplicity properties consistent with the field population as measured by Moe & Di Stefano (2017). A total of 950 binary and 50 single systems were generated for each star in the curated sample, for a total of $\sim$5 million systems with $\sim$20 million synthetic spectra. Period distribution in the sample is log-normal, with $\mu_{\log P}=4.6$ and $\sigma_{\log P}=2.2$, where P is in days. This distribution is largely consistent with the more commonly used Raghavan et al. (2010) for the close binaries, with the primary difference occurring at wider systems due to the difference in accounting for tertiary companions. The minimum period was set at 2 days, which is close to the minimum period observed in YSOs. The main sequence stars do have companions with shorter separations, however, because young stars are still inflated compared to their older counterparts, they cannot yet sustain shorter orbits yet. Some of these companions will migrate closer in through the magnetic breaking and tides. The distribution of mass ratios is uniform above $q=0.08M_{\odot}/M_{1}$, with an excess of twins ($0.95<q<1.00$) found at close separations; the excess is characterized at 20% probability at the separations $a<0.1$ AU, and 0% at $a>200$ AU (on top of the underlying uniform distribution), decreasing linearly with respect to $\log a$ (Figure 9a). The eccentricity distribution is uniform from 0 to $e_{max}$, where $e_{max}$ is defined by

[TABLE]

and $k$ is the tidal circulization period, which was set to 2 days. For computational efficiency, $e_{max}$ was capped at 0.98 (Figure 9b).

The orbit of each synthetic binary was sampled at epochs consistent with the actual observations of the reference source to obtain the RVs of the simulated primary and companion, computed at each epoch separately. Appropriate PHOENIX synthetic spectra (Husser et al., 2013) were chosen for each synthetic system: for the primary, both $T_{eff}$ and $\log g$ were chosen as the best template match to the reference source; the parameters of the secondary were interpolated from the PARSEC-COLIBRI isochrones (Marigo et al., 2017) for a given estimate of age, primary mass, and mass ratio. Both spectra were doppler shifted to the corresponding RVs, broadened by the kernel corresponding to the IN-SYNC $v_{rot}\sin i$ of the reference source (if a reference source was identified to be an SB2, the FWHM of the CCF peak that corresponds to the primary component was used instead). The flux was then scaled from the stellar radius to the appropriate distance222In the case of Taurus there are two separate populations projected on top of each other (Luhman, 2018), although the distance and proper motions are the only criterions in which it is possible to reliably separate them, and thus in absence of Gaia astrometry for some of the sources for the sake of simplicity both of them were treated as a single region. Spectra of both components were co-added, interpolated over the relevant wavelenth range, and the pixel noise profile of the reference source was applied. The CCFs of the resulting spectra were computed, and deconvolved into components for SB2 identification using the same routines as applied to the actual APOGEE observations. It proved impractical to process $2\times 10^{7}$ synthetic spectra through the IN-SYNC pipeline to provide RV estimates that could be used, as the IN-SYNC pipeline takes several hours to process a single epoch spectrum. Instead, a single Gaussian was fit to the CCF, and the Gaussian centroid is used as an approximation of the single measured RV in the spectrum. A check was performed to confirm that the recovered properties are consistent with Figure 8. Then, the intrinsic properties of the recovered systems, both SB2s and RV variables, are examined.

The recovered $\Delta v$ of SB2s with a quality flag of 4 is well matched by the original input (Figure 10), with a false positive identification rate among the single stars of only 0.06%. The recovered FWHM of flag 4 sources is also typically representative of the injected $v_{rot}\sin i$. Our recovery procedure overestimates the RV separation of a substantial number of SB2s with quality flag 3, particularly among rapid rotators with the injected v_{rot}\sin i$$\sim60 km s*-1*. This occurs because such binary systems may produce asymmetries in the wings of the CCF even with the intrinsic $\Delta v$ of just a few km s*-1*. While the Gaussians in the best fit do reproduce the shape of the CCF, and the sources involved are real binaries (the false positive rate among single stars is 2.8%, assuming that none of the systems are spotted), neither the recovered RV separation nor the FWHM are representative of the true RVs of the components or the $v_{rot}\sin i$ in Flag=3 systems.

Next we examine the distribution of the orbital parameters of the recovered SB2s, defining completeness as a ratio between the number of systems identified as a binary and the total number of synthetic binary systems, marginalized over a given parameter. In terms of the overall probability of detection of a system as a spectroscopic binary, for SB2s, almost no detections occur at $\Delta v$ between the primary and the secondary (S_{max}$$\equiv\mid v_{1}-v_{2}\mid) below 20 km s*-1*. The completeness of the recovered sample rises rapidly between 20<$$S_{max}$$<40 km s*-1*, and then remains constant at S_{max}$$>40 km s*-1*. Among the RV variables, it is possible to detect the amplitude of RV variability between epochs (V_{max}$$\equiv max(RVjd1)-min(RVjd2)) down to 1 km s*-1*, but the recovered fraction reaches its maximum at V_{max}$$\sim10 km s*-1*  and remains constant at higher $V_{max}$. The highest probability of detection occurs at very short periods, with almost no binaries with separations beyond 10 AU identified (Figure 11). Highly eccentric systems increase the likelihood of an identification of SB2s, because only a single epoch is needed to catch such a system at the pericenter to detect the maximum RV separation. On the other hand, RV variables with high eccentricity are more difficult to identify, as multiple epochs with sufficiently high RV offsets are needed. The effect of eccentricity on recovery of both RV variables and SB2s is less pronounced for systems with periods $<$200 days (Figure 12). A high mass ratio increases the probability for detection of SB2s. Among the RV variables, however, detections are maximized at $q\sim 0.5$. At higher $q$, the companion makes a stronger contribution to the spectrum, even if the RV separation is too small for the system to be identified as an SB2. In these cases, the pipeline measured RV is influenced by both the primary and the secondary, and tends toward the stable barycentric RV of the system. Binaries whose components have low $v_{rot}\sin i$ are more reliably recovered, as their narrow spectral features produce similarly sharp CCF peaks, which allow more accurate RV measurements and a better separation of components even at low $\Delta v$. Multiple epochs do help with identification, but not significantly - sources with 10 epochs are only 1.5 times more likely to be caught in a favorable configuration for detection compared to a single epoch. RV variables with more than 3 epochs are not more easily identified via $\chi^{2}$, although there is a modest increase of RV variables identified via the slope.

It should be noted that because the synthetic sample was constructed using the real observations as a basis for the distributions of stellar properties, any unresolved SB2s in the data would have a somewhat inflated $v_{rot}\sin i$, resulting in a lower recovery rate for these sources, thus potentially underestimating the total MF in the model compared to the data. It is not expected, however, that this effect should be significant as the difference in completeness between the true and a slightly inflated $v_{rot}\sin i$ is small, and only a few sources will be affected.

As all the sources are individually represented in the model, we also compare the completeness of the overall sample as it is split among the individual regions. In general, the dependence of the completeness on the underlying stellar parameters is consistent among all of the regions, in shape if not in the absolute normalization. The only exception to this is the mass ratio, because it is driven not only by the typical masses of the stars in a given sample, but also by age (Figure 13). This does not strongly affect RV variables, but appears to have a strong effect on SB2s. In the synthetic spectra is easier to recover sources across all mass ratios in younger systems with lower $\log g$. On the other hand, as the stars shrink to their main sequence sizes, the sensitivity to low $q$ systems decrease significantly. This occurs because in the younger systems, the companion tends to have a lower surface gravity than the primary, thus the luminosity ratio between the secondary and the primary is higher than what it will be in a more evolved system, where the trend of $\log g$ is reversed.

Among the SB2s, the amplitude of the components in the CCF typically corresponds to the ratio of the effective temperatures as a function of the $T_{eff}$ of the primary (Figure 14). Instead of the ratio of the amplitudes, it is also possible to construct a similar relation using the integral of the individual Gaussian components, although with a greater degree of scatter. The trends observed in Figure 14 are not as apparent when the $T_{eff}$ ratio is replaced by $q$: the direction of dispersion with the relation to $T_{eff}$ of the primary reverses half-way through the range.

6 Results

In this section we analyze the trends in the data in relation to the model of field multiplicity.

6.1 Mass ratio

For SB2s that have multiple epochs and substantial $\Delta v$ between components, we measured mass ratios by fitting a slope to RVprimary vs RVsecondary (Wilson, 1941). We included in the fit all of the deconvolved components down to Flag=1, assuming that a robust Flag=4 detection was available in at least one epoch. When components were unresolved, both the primary and the secondary were assigned to the same central RV. Any apparent SB3s and SB4s (8 sources, most likely compact higher order multiple systems) were excluded, as the sample does not contain enough epochs to reliably fit their relative motions. If a component was apparent during visual examination but was not recovered by the pipeline, the peak of that component was fitted manually. Fits to determine $q$ could be performed for 78 systems, the resulting mass ratios are included in Table 1. In general, there was no confusion in tracking the companions. Primary/secondary misidentification resulted in fits that returned q$>1$ for 16 stars; the identification were consistent across all epochs, however, and we correct the error by reporting the inverse of the derived mass ratio (i.e., $q=1/q_{derived}$). The only confusion occurred in epochs in which companions were unresolved if the CCF had a low SNR component not associated with a real source.

Previously, Fernandez et al. (2017) compared the distribution of mass ratios among young SB2s detected by APOGEE with results from the other binary surveys. They found that equal mass binaries dominated their sample, with a tail stretching towards lower $q$. They found that the overall distribution is consistent with what has been observed by Pourbaix et al. (2004) for more evolved stars, but they interpreted the excess of twins as a detection bias, where high mass systems have a larger orbital velocity amplitude than lower mass ratio systems. Tokovinin (2000) and Moe & Di Stefano (2017), however, concluded that at short separations (within 200 AU), the overall distribution of $q$ is not entirely uniform but that there is a significant excess of companions with $0.95<q<1$ on top of the uniform distribution.

Examining the mass ratios measured from the data, we find a similar distribution to that observed by Fernandez et al. (2017), with a large excess of twins (Figure 15). Examining the completeness of the recovery of synthetic binaries as SB2s, we find that the probability of detection does rise somewhat for systems with high mass ratio, but it is not the most dominant factor: the completeness only changes from $\sim$5% of all binaries at $q\sim 0.4$ to 7% at $q\sim 0.9$. The reason for the apparent discontinuity at $q=0.95$ in the total distribution in Figure 11 is due to the fact that this excess of twins at short separations is included in the model population we adopt for our completeness calculations, and it is the small separations, not the near-equal mass ratios, that drives these twins to be twice as likely to be identified as SB2s. Without the intrinsic separation dependent excess of twins, such systems would not dominate the sample, i.e., the observational biases cannot account for their presence. Therefore, the interpretation used by Fernandez et al. (2017) regarding their presence does not apply, and the model of mass ratios from Moe & Di Stefano (2017) does accurately represent the data.

There are no strong differences in $q$ as a function of the evolutionary class, nor as a function of period (Section 6.5), nor is there a trend as a function of $T_{eff}$. However, curiously, there are 3 sources that have a peculiarly low $q\sim 0.1$; this is significantly lower than the minimum of 0.2 that is recovered in the synthetic observations. Two of these sources, 2M03281336+4856154 and 2M03301892+4943348, are located in $\alpha$ Per, the third, 2M03432619+2602308, is found in the Pleiades. The primaries of these systems are $\sim$0.2 M⊙objects, making their companions to be a part of the brown dwarf desert.

6.2 Normalization

The total MF for SB2s across all epochs is 2.8$\pm$0.2%, this MF rises only slightly to 3.1$\pm$0.3 in SB2s with at least 3 epochs. Among RV variables with at least 3 epochs, the MF is 6.3$\pm$0.5%. Because there is some overlap between the sources classified as both SB2s and RV variables, the combined multiplicity in the sample is 8.3$\pm$0.5%.

In order to compare the MF in all of the regions in a uniform manner, we use the recovered fraction of the synthetic binaries as a proxy for what is expected in each region given the underlying model of the distribution of companions. This allows us to take into the account any differences that originate from the different observing strategies between individual regions. We do not report the MFs measured in individual regions but rather compare if the measurement is over or under-abundant in comparison to the model.

For computational expediency, almost no single stars were generated as a part of the synthetic comparison sample. However, in order to calculate the MF in the model that is most representative of the data, it is necessary to incorporate single stars in some way. We account for the presence of the single stars by scaling down the MF from the synthetic sample. We originally considered a normalization that would maintain 16.5% of systems with separations less than 10 AU (with MF defined as $\frac{B+T}{S+B+T}$, where S is the number of single systems, B is the number of binaries, T is the number of triples with the given separation range, etc). This relation is interpolated from Moe & Di Stefano (2017), corrected for the typical masses ($\sim$0.6 M⊙) and metallicities (Z=0) of the stars in the primary sample. This scales the overall synthetic MF by a factor of 2.4.

However, the MF does have a strong dependence on the mass of the primary (e.g., Duchêne & Kraus, 2013; Raghavan et al., 2010; Ward-Duong et al., 2015), with the overall MF increasing by a factor of $\sim$1.7 from M to G dwarfs, with the effect most pronounced in wide binaries. Comparatively modest increase is observed in the MF of close systems, with only a factor of $\sim$1.2 increase from 0.6 to 1 M⊙, although it jumps by a factor of $\sim$1.9 from 1 to 2 M⊙ primaries (Moe & Di Stefano, 2017; Murphy et al., 2018). While most of the sample consists of K type objects (56% of sources have 3900<$$T_{eff}$$<4800 K, correlates to masses of $\sim$0.5–1 M⊙), overall the sample spans from close to the substellar boundary ($T_{eff}$=2300 K) up to early G type stars ($T_{eff}$= 6000 K, with masses as much as 2 M⊙ – in YSOs, by the time they reach the main sequence they will become hotter). The exact distribution of primary masses in the APOGEE sample may have an effect on the recovered MF. Furthermore, the APOGEE observations are limited in brightness down to $H\sim 13$ mag. Due to the Malmquist bias of the second kind, binaries may be overrepresented among the low mass sources in our sample because they were more likely to meet or exceed our targeting limit due to being brighter than their single counterparts. Due to the complex targeting strategy, it is difficult to reproduce this entirely through forward modeling, but close to the magnitude limit, it may artificially raise the MF by 10–100% in a given mass or flux range.

In Figure 16 we compare the MF as a function of both $T_{eff}$ and the $H$ band. In each of the bins (here, and in the subsequent subsections), we restrict the sample from the model to only those synthetic sources that were produced in reference to the sources that fall into a given bin. Both distributions show a decrease in MF from hotter to cooler stars by a factor of $\sim$2 from 4000 to 6000 K. Through performing linear regression, the rising slope of MF relative to the H band fluxes is inconsistent with being flat at 3.7$\sigma$ level with H$<$13, and 2.8$\sigma$ relative to $T_{eff}$ for T_{eff}$$>3800 K. This is largely consistent with what is expected from physical models of the dependence of multiplicity on primary mass. While each individual region has a different conversion from mass to $T_{eff}$ and apparent H mag due to different distances and ages, this trend holds across all of the regions. However, it is impossible to correct the increase in the MF in both H band fluxes or just $T_{eff}$ by using a normalization that relies on just one of these parameters, necessitating including both of them. This suggests that in addition to the astrophysical reasons for the increase, there might also be an effect from the increased SNR that helps to identify binaries that is not entirely replicated by the model.

The MF also increases for the faintest stars, with the strongest effect occurring at the observational magnitude limit, consistent with the Malmquist bias. This translates to an increase in MF for sources with T_{eff}$$<3500 K in the full sample, and while this is also observed in all of the individual regions covered by the survey, depending on the distance the location of the break would be different.

To correct for both of these effects, we apply a normalization scaling

[TABLE]

to the model. This was obtained by fitting the MF dependence on H band across the entire range, and then fitting the deficiencies in $T_{eff}$. It should be noted, however, that this normalization has very little effect on the results of the subsequent sections as compared to a simpler $MF(<10AU)=16.5\%$, as both mix high and low mass stars together and average out any strong $T_{eff}$ or flux dependence.

6.3 Overall distribution between the regions

Overall, the distribution of spectroscopic binaries averaged across all regions appears to be consistent with the model; this holds for both SB2s and RV variables (Figure 17). The same applies to all of the individual regions in the survey – while individual regions may show a slight non-statistically significant excess or deficit in one type of companions, those deviations average out in the combined sample. There are only a few exceptions:

• •

SB2s (but not RV variables) are deficient in NGC 1333 by a factor of 5.3 (4.3$\sigma$)

• •

RV variables (but not SB2s) are in excess in Taurus and L1641 by a factor of 1.6 (1.5–2$\sigma$)

• •

In Orion C, both SB2s and RV variables are in excess by a factor of 1.95 (1.9–2.2$\sigma$)

• •

the MF measured from RV variables is somewhat elevated (by as much as $\sim$50%) compared to what is measured from SB2s, although the trend is reversed in the older regions ($\alpha$ Per and Pleiades).

The excess of companions in Orion C may be explained by the membership assignment. Because this region and Orion D are projected on top of each other in the plane of the sky, there is some difficulty in distinguishing them apart. When possible, we used membership from Kounkel et al. (2018), otherwise a weighted average RV cut of 25 km s*-1*was used. While this cut is sufficient to reliably separate single stars in these populations, in the case of spectroscopic binaries this results in the systems that are currently being blueshifted being assigned to the Orion D, and those that are redshifted being assigned to the Orion C. Because Orion D has more sources in total by a factor of 1.3, the companions will be disproportionally represented in Orion C, resulting in an elevated MF. Averaging these two regions together brings them both into a closer agreement with the model.

6.4 YSO Class

We examine the effect of the multiplicity on the presence of the protoplanetary disk. Out of 1882 disk-bearing Class II sources, 59 are found to be RV variables, and 38 are SB2s. Out of 3125 diskless Class III sources, 116 are RV variables, and 103 are SB2s. Class IIs are deficient in SB2s compared to the model by the factor of 1.9 (4.2$\sigma$, Figure 17); this deficit is recovered in Class III, which show a slight excess over the model of a factor of 1.15 (1.3$\sigma$). On the other hand, RV variables have no effect on the protoplanetary disks: both Class II and Class IIIs show a similar excess of Class III SB2s. The physical implications of this are discussed in Section 7.

This does help to shed light on the origin of some of the minor differences in the distribution of companions in the total sample (Figure 17). This includes the deficit of SB2s in NGC 1333, as it is one of the youngest clusters in the study. However, this cannot account for the excess of companions in Taurus and L1641.

It should be noted that while veiling is systematically different between Class II and Class III systems, introducing it as an additional parameter in the construction of the synthetic binaries does not affect the recovery probability of close binaries.

The difference between RV variables and SB2s as a function of the evolutionary type does not appear to be an effect of RVs being artificially scattered by means other than a companion (e.g. star spots). To test this we also performed a comparison restricted only to the systems with low v_{rot}\sin i$$<20 km s*-1*. These slow rotators have a narrow profile of their CCF that does not show as much variability due to spots in the synthetic spectra. While spots may alter RVs, the resulting shift is smaller than the typical RV uncertainties in such systems. The difference between SB2s and RV Variables between Class II and Class III systems remained consistent when the sample was limited only to slow rotators.

Similarly, the difference cannot be explained by the systematic difference in $v_{rot}\sin i$ between Class II and Class III systems. Class II systems do tend to have $v_{rot}\sin i$ that is higher than that of Class III systems, but only by a few km s*-1*  not enough to cause significant discrepancy, and the individual $v_{rot}\sin i$measurements are propagated to the model, with the comparison relying only on those synthetic spectra that are referenced to the sources that are considered in the individual bin.

6.5 Close systems

We use the data to test the log-normal period distribution of the model population with $\mu_{\log P}=4.6$ and $\sigma_{\log P}=2.2$, and the minimum period of 2 days.

In order to limit the sample to just the closest systems, for each system identified as a binary in both the data and the synthetic sample, we measured $V_{max}$ for the systems identified as RV variables, and $S_{max}$ for SB2s. These properties are highly dependent on the inclination of the system, the cadence of the observations, and the probability of catching the system near the maximum velocity separation. These factors reduce the observed velocity amplitude, making $V_{max}$ and $S_{max}$ only lower limits to the intrinsic maximum RV separation of the system. However, sufficiently high amplitudes in either of these parameters will nonetheless identify systems with shorter periods.

Using the synthetic spectra, we measured the $V_{max}$ and $S_{max}$ cuts that rule out systems with periods longer a particular threshold with 66% certainty (Table 6). With this criterion, however, it is impossible to conclusively separate sources into individual period bins, only to give a probability that the period of the system is smaller than a given threshold. For the shortest period systems this introduces an increasingly more stringent cut, making it more likely to discard a system whose actual period is shorter than the desired upper limit. As such, for example, at 200 days it is possible to identify 72% of all such systems with these cuts, but only 32% of systems for periods within 10 days can be recovered, increasing a false negative fraction.

Examining the MF as a function of the cuts in $V_{max}$ and $S_{max}$ shows that the period distribution in the data is consistent with the model at nearly all separations. The exception to this are sources with periods less than 5 days, which show a weak 1$\sigma$ deficit in both RV variables and SB2s (Figure 18). This deficit is not entirely a product of the model’s 2 day minimum period, however, as a raising the minimum period would affect other bins with higher maximum period (as they also include systems with short periods), making them less compatible with the model. It is possible that there may be a physical explanation for the difference between the model and the data for these short period systems. For example, because YSOs are inflated compared to main sequence stars, YSOs should originally have wider separations, but some would evolve dynamically and harden to be closer together, lessening the deficit at short periods. It is also possible that there might be some differences in the manner RVs are measured in the synthetic sample could predominantly affect the most widely resolved systems. Additionally, because we probe only the likely maximum period and not the exact period distribution of each system, it is possible to use different statistical weighting for converting from velocity to period, which could push this sub-5 day deficit to either shorter or somewhat longer periods.

Separating the sample according to their evolutionary classification shows a similar distribution; while there might be an excess or a deficit that is seen in the total sample, it remains consistent at all separations with exception of the shortest periods.

6.6 Stellar Density

Because the sample spans a large variety of environments, we examine the effect of multiplicity on the local stellar density in the primordial population.

A measurement of stellar density at the position of each source in the sample is a non-trivial task because of the non-uniform membership list in each of the regions in question. Disparate distances and disparate ages result in different completeness limit in each cluster. No single method is capable of identifying all of the members of a population. Moreover, fundamentally, stellar density is a three-dimensional problem. While Gaia DR2 made incredible strides in making stellar distances accessible, how deeply positioned a given source might be in a given cluster still cannot be known as precisely as a relative position in the plane of the sky.

Nonetheless, we attempt to estimate stellar density for the sources in our sample. We begin with identifying the population corresponding to each young cluster and star forming region using the Gaia DR2 data, by using TopCat (Taylor, 2005) to make a rough selection in the position, parallax, proper motion, and color-magnitude space that includes each cluster’s sources in the curated catalog. We also made a cut perpendicular to the main sequence at the magnitude limit of the furthest cluster in the sample (i.e., NGC 2264), and discarded all sources lower than this limit in the closer regions. Similarly, we discard all high mass sources brighter than the sources that correspond to T_{eff}$$\sim6000 K in each cluster. Unfortunately, any sources that are too heavily extinct, or sources that have irregular kinematics despite being cluster members, cannot be counted in this sample.

We then measure a projected separation between the sources in the curated sample and the resulting catalog. This separation is converted into parsecs using the parallax measurement of the source in question. To estimate the local density, we record the projected separation to the fourth nearest neighbor (NN4), similarly to the approach used by Kounkel et al. (2016b). The relative distribution of densities is shown in figure 19. Again, however, it should be cautioned that this is just an estimate, made worse by any possible line of sight effects, such as may be the case in Orion C/D and in Taurus.

Looking at multiplicity as a function of NN4 (Figure 20), MF appears to have a maximum at NN4$\sim$0.2 pc, or $\Sigma_{*}\sim 30$ pc*-2* that is primarily driven by Class III SB2s, remaining mostly flat with a possible decline at the lowest densities of $\Sigma_{*}\sim 30$ pc*-2*

7 Discussion

The close binary fraction and properties of young stars provide an invaluable insight into the mechanisms behind formation and evolution of the multiples. Using APOGEE spectra of pre-MS stars in various star-forming environments and accounting for selection effects, we have found their close binary fraction is consistent with the field population with separations $<10$ AU. This result is consistent with models in which the majority of close binaries form during the embedded Class 0/I phase within the first $\sim$1–3 Myr (Moe & Kratter, 2018; Bate, 2012). The measured binary fractions inferred from both SB2s and RV variables are consistent with each other, validating our identification pipelines and corrections for incompleteness are robust. We also find the bias-corrected pre-MS spectroscopic binary fraction increases by a factor of $\sim$2 with increasing temperature and brightness up to $T_{eff}$= 6000 K, which is similar to the observed change in the field close binary fraction with respect to primary mass.

There might be a slight $\sim$30% deficit at the shortest of periods P$<$5 day. The observed deficit suggests a minority of very close binaries migrated during the main-sequence phase via tides (Fabrycky & Tremaine, 2007; Moe & Kratter, 2018).

The mass ratio distribution of young close binaries is also largely consistent with their field counterparts. Specifically, the mass ratios are uniformly distributed with a small excess of twins with $q>0.95$. Close pre-MS binaries also follow the short-end tail of a log-normal period distribution similar to that observed for low-mass and solar-type MS binaries in the field. There might be a slight deficit at the shortest of periods $P<$5 days.

The overall distribution is largely mirrored by the individual regions. There does not appear to be a dependence in MF on age, as it appears to be consistent between young clusters (e.g., ONC), as well as those that are more evolved (e.g. $\alpha$ Per, Pleiades). This contradicts what has been observed previously by Jaehnig et al. (2017), however, their trend has been primarily driven by a small number of RV variables in Pleiades.

Within each cluster younger than a few Myr, however, there appears to be a strong deficit of SB2s among sources with protoplanetary disks. The most extreme deficit is found in NGC 1333. This deficit of SB2s in disked systems is present across all separations. Close binaries are known to to affect protostellar disks – disks are half as common in visual binaries with projected separations of $<$40 AU than they are around single stars, and their frequency may further decrease at closer separations (Jensen et al., 1996; Cieza et al., 2009; Kraus et al., 2012; Harris et al., 2012; Rodriguez et al., 2015). Thus, the decreased MF in the Class II SB2s is consistent with the interpretation interpretation that close binaries either disrupt or accelerate the evolution of their natal disks. On the other hand, RV variables do not appear to affect disks, and it is consistent with Kuruwita et al. (2018) conclusions regarding the observations of Class II RV variables in Upper Scorpius and Upper Centaurus-Lupus.

It is not clear why there is such a stark difference between RV variables and SB2s regarding their effect on the protoplanetary disks. Observationally, the primary difference between these systems is their mass ratios. A possible interpretation is that more equal mass companions are the result of accreting a substantial fraction of the disk, thereby reducing the disk mass and accelerating the transition to the Class III phase.

Few regions, most notably, Taurus and L1641, do show a slight elevation in MF among RV variables (but not SB2s, thus affecting only the systems with lower $q$). It is not immediately apparent why this might be the case.

There may be some weak trends with stellar density, with multiples being most common at the projected $\Sigma_{*}\sim$30 pc*-2*, and declining at the lower end of the density distribution. While this trend may not be statistically significant, such observations pose important constraints in analyzing the degree to which the primordial binaries are disrupted in dense environments.

We acknowledge Adrian Price-Whelan for the discussion, and Claire Murray for the clarification of the functionality of GaussPy. M.K. and K.C. acknowledge support provided by the NSF through grant AST-1449476, and from the Research Corporation via a Time Domain Astrophysics Scialog award (#24217). M.M. and K.M.K. acknowledge financial support from NASA under grant ATP-170070. C.R-Z acknowledges support from program UNAM-DGAPA-PAPIIT IN 108117, Mexico. A.R-L acknowledges financial support provided in Chile by Comisión Nacional de Investigación Científica y Tecnológica (CONICYT) through the FONDECYT project 1170476 and by the QUIMAL project 130001. KPR acknowledges CONICYT PAI Concurso Nacional de Insercio ̵́n en la Academia, Convocatoria 2016 Folio PAI79160052. Support for JB is provided by the Ministry for the Economy, Development and Tourism, Programa Iniciativa Cientica Milenio grant IC120009, awarded to the Millennium Institute of Astrophysics (MAS). Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is www.sdss.org. SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.