AN HST SURVEY FOR 100–1000 au COMPANIONS AROUND YOUNG STELLAR OBJECTS IN THE ORION MOLECULAR CLOUDS: EVIDENCE FOR ENVIRONMENTALLY DEPENDENT MULTIPLICITY

The HST survey targeted 280 protostars in the Orion A and B clouds that were identified with 1.2–24 μm photometry from the Spitzer Orion Survey (Kryukova et al. 2012; Megeath et al. 2012). This sample of protostars was selected to be the entire sample of Spitzer-identified protostars with predicted 70 μm fluxes in excess of 42 mJy; this flux cutoff was used to ensure they would be detectable with Herschel/PACS as part of the Herschel Orion Protostar Survey, an open-time key project selected for the Herschel mission (Megeath PI). The HST survey input catalog is not a complete sample of all Spitzer-identified protostars satisfying the 70 μm flux cutoff for two reasons. First, the methodology for identifying protostars in the Spitzer data was still evolving at the time the input catalog was defined; hence, several protostars in Orion were missed. Second, the flux cutoff was based on an extrapolation from the Spitzer 3.6–24 μm photometry, and in many cases the measured flux was substantially different. Furthermore, there is a small degree of contamination from reddened pre-main-sequence stars and galaxies. Hence, only 248 out of the original 280 sources are still classified as protostars. The subsequent IRTF protostar survey targeted 66 Spitzer-identified objects with the m_3.6μm < 10.4 mag.

As part of program GO 11548, 87 orbits of observations were successfully executed with the NICMOS camera in 2008 August and September before the failure of the cryocooler. With the deployment of the WFC3 in 2009 June, the program was redesigned for the near-IR channel of that instrument. The two primary changes were the elimination of the F205W band observations due to the wavelength cutoff of WFC3 at 1.75 μm and the inclusion of multiple protostars in a single pointing with the larger FOV of WFC3. The remaining protostars in the target catalog were observed in 126 orbits obtained between 2009 August and 2010 December.

2.2.1. The NICMOS Observations

A total of 87 protostars were observed with the NIC2 camera, which has a 192 × 192 FOV and a pixel size of 0075. Due to the obscuration of the molecular clouds, some of the guide stars were of marginal brightness; consequently, the guiding failed for nine protostars. Each protostar was imaged in a single orbit split between the F160W and F205W bands. The NIC-SPIRAL-DITH pattern was used with four points and a 1'' spacing; this pattern moved the source to the four corners of a square with 1'' sides in the focal plane. The MULTIACCUM mode was used to maximize the dynamic range of the observations. The predefined SAMP-SEQ = STEP16 sequence was selected for the MULTIACCUM; this mode samples up the ramp starting with readings spaced by 0.303 s and then increases the interval up to a maximum of 16 s. A sequence of 25 steps were chosen for the F160W band and 16 steps for the F205W band. This resulted in total integration times of 1215.4 s and 767.6 s in the F160W and F205W bands, respectively.

The initial calibrated frames were created with the calnica command from the IRAF NICMOS data-reduction tools provided by the Space Telescope Science Institute (STScI); this reconstructed the MULTIACCUM data into units of DN/S and performed the flat fielding and dark subtraction on the data. The resulting frames showed a spatially varying residual background. This background was most obvious at 2.05 μm and is due in part to thermal emission from the telescope. To reduce this artifact, which appeared to change very little over time, all images that had no detections of a protostar or only faint detections were used to generate a residual "sky" frame that was not contaminated by stars or nebulosity. The residual "sky" was then subtracted from all of the data. This procedure was done separately for the F205W and F160W frames.

Even after the subtraction of the residual "sky" frame, the four individual quadrants of the NICMOS array showed varying offset levels. The offset levels were particularly pronounced when a very bright source was present in the image. To match the offset, every image was individually adjusted by adding a constant value to each of the four quadrants with a custom interactive data language (IDL) program. These values were varied until a discontinuity was no longer apparent upon visual inspection of the image.

A generic bad-pixel mask was downloaded from the STScI and then was supplemented by finding additional bad pixels through visual inspection. The final masks included several dead pixels, the gap between the four detector quadrants, and the circular region of the array blocked by the NICMOS coronagraph. Finally, the four dithered images were combined with the calnicb command.

While the relative astrometry within a NICMOS image is precisely calibrated, the absolute pointing can vary by as much as 15. In Section 3, we describe the detection of point sources in the NICMOS frames and the measurement of their positions and magnitudes. To correct the absolute pointing, we compared the positions of the point sources in the NICMOS fields to analogous sources in the Spitzer Orion Survey. For every NICMOS source for which there was a reliable Spitzer counterpart, we calculated the offsets between the NICMOS and the Spitzer coordinates from Megeath et al. (2012, Figure 1). The Spitzer positions, which have been refined through a comparison to the 2MASS point-source catalog, have positional errors <03.⁸ In some of the NICMOS images there were no point sources detected; the protostars in these images were either too faint or were surrounded by bright nebulosity. To facilitate the position refinement of these data, and to average out uncertainties in the Spitzer positions, we applied a global position correction to the images in bulk, as opposed to correcting each object individually. We used the median offsets between the NICMOS and IRAC positions; these were δR.A. = 0161, δdecl. = −0376 before 2008 August 13, and δR.A. = 0298, δdecl. = −0628 after that date. Although these offsets improved the pointing significantly (Figure 1), there were six outliers where the total offset was greater than 1'' or the position offset had a different sign as compared to the average. In these cases, individual offsets were applied so that the NICMOS coordinates corresponded to those of the Spitzer counterparts.

Zoom In Zoom Out Reset image size

2.2.2. The WFC3 Observations

A total of five nights of observations were obtained at the NASA IRTF with the NSFCAM2 camera, which houses a 2048 × 2048 detector with a FOV of 80'' × 80'' and a pixel scale of 004 per pixel. Our observations were obtained in the L' band, centered at 3.76 μm, which was selected as a compromise between maximizing the flux from the observed protostars, all of which exhibit rising SEDs with increasing wavelength, and maximizing the sensitivity of the observations, which declines with wavelength due to the rising thermal background of the Earth's atmosphere. The typical seeing was 075, which corresponds to a spatial resolution of ∼300 au at 420 pc. We imaged 36 fields in 2010 January and 24 fields in 2010 December, all of them containing at least one confirmed protostar. A five-point dithering pattern put the target at the center of the frame and in each of the quadrants. Depending on the brightness of the central source, the dither pattern was repeated three to five times, resulting in 15 or 25 dithered frames per targeted field. Each frame had an exposure time of 30 s using 1 s integrations and 30 coadds. In total, we obtained images of 98 Spitzer-identified YSOs.

The NSFCAM2 data were reduced using custom IDL routines. The flat fields were constructed using frames targeting different sources that were imaged close in time but had differences in their air masses of ∼0.4. The frames were organized into pairs, and the low air mass frame of the pair was subtracted from the high air mass frame. The resulting images were then normalized and stacked, and a sky frame was constructed by taking the median value at each pixel in the stack. A flat field was created for all five nights of the observations; these were consistent to 2% from night to night within a given observing run.

To subtract out the sky contribution, all the flat-fielded frames toward a given source were then grouped into sets of five consecutive frames; these correspond to the five spatially offset fields of our five-point dither pattern. A constant offset was added to each image so that the median value of the pixels in each of the five frames was equal to the median pixel value of all five frames. We then stacked the images and produced a sky frame by calculating the median of the five values at each pixel location. This approach provided enough redundancy to eliminate stars and compact nebulosity as well as artifacts such as cosmic rays, but ensured that the intensity of the sky did not change substantially during the sky measurement. Each of the five frames was then subtracted by the sky frame.

NSFCAM2 at the time of the observations utilized an engineering-grade detector that possessed many low-responsivity or dead pixels. A bad-pixel mask was created by removing the pixels where the flat field showed more than 10% deviation from unity. The sky-subtracted, flat-fielded, and masked frames were then mosaicked together using a custom IDL routine provided by R. Gutermuth. This routine determined the offsets between overlapping frames by comparing the positions of sources in common between the frames. The images were then interpolated onto a mosaic grid using a bilinear interpolation. This was done iteratively, with the mosaic being assembled one image at a time. During each iteration, the signal levels of the mosaic and the overlapping sky frame were compared, and a constant value was subtracted from the new image so that it had the same signal level as the mosaic.

During the first observing run, two different standard stars were observed: HD 40335 and SAO 112626. In the second run, we expanded to five standard stars: HD 40335, SAO 112626, GL105.5, HD 44612, and HD 1160. During both runs, one standard, HD 40335, was observed throughout the night at different air masses. The data were reduced in an identical manner to the science targets, but the photometry was measured for the individual frames instead of using the mosaicked data. We adopted the tabulated L'-band magnitudes from the UKIRT faint infrared standard catalog,⁹ with an exception of GL105.5, the photometry for which came from the UKIRT bright infrared standard catalog.¹⁰ We determined the zero-point magnitudes as a function of air mass for each night.

The images were registered using absolute coordinates from the Spitzer Orion Survey (Megeath et al. 2012). After identifying counterpart point sources in the Spitzer data, a linear transformation between the pixel values of the NSFCAM2 sources and the absolute coordinates were determined for each of the NSFCAM2 mosaics using least-squares fits. By adopting a linear transformation, we disregarded any nonlinear distortion in the IRTF data; this was necessary given the small number of sources in common. The adopted linear coordinate transformation is only accurate to half an arcsecond around the edges of the mosaic; this level of precision is adequate for the analysis in this paper.

Point sources were identified with PhotVis, a tool for IDL that facilitates both automated and interactive source finding and photometry (Gutermuth et al. 2008). The sources were found using an automated search with a 20σ detection limit. After the automated source finding, each field was visually examined, previously unidentified sources were added, and false detections of noise and galaxies were eliminated. We adopted a zero point of 24.51 mag. To estimate the uncertainty, we adopted the gain for the exposure time of 2496.17 s, as given by the HST documentation. We limited our catalog to point sources and excluded extended or nebulous sources.

Point-spread function (PSF) fitting photometry was then performed using the IDL implementation of DAOPHOT. We used a modified version of DAOPHOT designed to mask out saturated pixels in the images (Kryukova et al. 2012). This allowed us to remove the contribution of bright stars, find faint companions, and perform simultaneous PSF fitting of the companions using a second iteration of DAOPHOT. In addition, the PSF fitting recovered the fluxes for several sources that were saturated. Fitting the data to the WFC3 PSF also allowed us to distinguish contamination from slightly extended galaxies initially identified as point sources by the automatic detection routine.

The WFC3 photometry is consistent with the 2MASS H-band photometry (Figure 2). In the cases where the 2MASS and WFC3 photometry disagreed, the WFC3 magnitudes are fainter, with the differences highest near the detection limit of the 2MASS data. This is expected since the lower angular resolution 2MASS point-source photometry can be contaminated by extended emission, and the faintest 2MASS sources are most affected by contamination.

Zoom In Zoom Out Reset image size

In contrast to the WFC3 data, the much narrower NICMOS FOV typically contained only a few sources. This allowed us to identify point sources through visual inspection with the ATV package in IDL (Barth 2001). The photometry was measured at the positions identified by ATV using the procedure "aper" from the IDL astronomy library (Landsman 1993). We adopted an aperture of 6.5 pix (04875) with a sky annulus from 10 pix (075) to 20 pix (15). The zero-point magnitude and aperture corrections are from the NICMOS documentation. However, a comparison between the small number of sources that were observed with both WFC3 and NICMOS showed a constant offset between the F160W magnitudes. Thus, we corrected the NICMOS zero point by −0.276, which was the median difference in photometry between stars detected in both the WFC3 and NICMOS data sets with 12 < m_F160W < 15. We applied the same offset for both F160W and F205W NICMOS photometry. Comparisons of the uncorrected NICMOS photometry to the WFC3 and 2MASS photometry are shown in Figure 2.

We subtracted the PSFs of the primary sources to search for faint objects hidden by the wings of the PSFs. Unlike the WFC3 images, the NICMOS images do not contain enough point sources to construct a PSF directly from the image; for this reason we did not use DAOPHOT, which relies on a PSF generated from the data itself. Instead, we created a PSF model with TinyTim for both the F160W and F205W bands (Suchkov & Krist 1998). A custom IDL PSF fitting routine was used to remove the contribution of the point source from the data; this was used to search for close companions. For sources where a companion was detected, simultaneous fits were used to determine the magnitudes for both objects. We adopted the values from the prefitted aperture photometry for the uncertainties of those sources. For the five secondary sources that were not detected prior to fitting the adopted uncertainties, we assigned them the uncertainty of a source found by NICMOS with a similar magnitude.

The PhotVis routine in IDL was used to identify point sources in the NSFCAM2 mosaics and measure their photometry. We searched for sources manually because the automatic finding routine identified too many spurious sources, due to noise spikes and artifacts in the data, particularly toward the edges of the frames. PSF fitting to the NSFCAM2 data was not performed because of the variations of the PSF with both position on the detector and between observations, due to changes in the seeing. A 50 pixel (2'') aperture, with a sky annulus from 50 to 73 (292) pixels, was typically used; however, for close binaries, we reduced the aperture to less than half the separation of the sources. In these cases, we maintained the same size for the sky annulus. We then calculated an aperture correction for these photometry results using several isolated stars in the same images. The uncertainty was calculated from the variation in the aperture correction between different isolated stars in the image.

Individual point-source catalogs were generated for all of the three cameras; these were then merged with each other and with the catalog of all point sources in the Spitzer Orion Survey (Megeath et al. 2012). For the NICMOS observations, there were 78 fields with usable data. Due to the high levels of extinction through the Orion molecular clouds, particularly in regions of high gas column density that typically surround protostars, and the small FOV of NICMOS, six of the fields did not have detections. A total of 72 frames observed in F205W and a total of 63 in fields in F160W had at least one detected object. The NICMOS catalog contains F205W photometry for 123 point sources and F160W photometry for 98 point sources. The WFC3 catalog contains 4984 point sources, including 128 HOPS objects, of which 123 are protostars and five are young stars with disks (Furlan et al. 2016) In addition to this, the larger FOV encompassed six additional protostars and 192 young stars with disks that were not part of the HOPS catalog. The final point-source catalog from NSFCAM2 has 337 YSOs.

A total of 178 protostars have point-source photometry in the F205W or F160W bands. An additional 23 were detected with NSFCAM2; these were either too red to detect with HST or, in the case of five objects, were not observed with HST. The remaining HOPS protostars that are not part of our final catalog were either too faint to observe at our targeted wavelengths or were extended and did not contain a point source. Since every source was visually identified or found using an algorithm that identified sources with peaks 20σ over the surrounding noise level, we did not filter the data based on their uncertainties or magnitudes. For the WFC3 catalog, we find that 99% of the sources have m₁₆₀ ≤ 23.3 mag and σ₁₆₀ ≤ 0.14 mag. The 95% limit for the NICMOS sources are m₁₆₀ ≤ 21.3 mag and m₂₀₅ ≤ 20.8 mag (we use the 95% since there are only about 100 sources in the NICMOS data). For the IRTF L' band, 99% of the sources are brighter than 14.1 mag and have uncertainties ≤0.16 mag. The photometry of dusty YSOs and their companions extracted from the final point-source catalog are presented in Tables 1, 2. We note that several of the sources imaged by both NICMOS and WFC3 show significant differences in their m₁₆₀ mag; these will be studied in future papers to search for variability among the observed YSOs and their companions.

To search for an excess density of sources, we utilize the mean surface density of companions (Larson 1995), which is related to the two-point correlation function (Bate et al. 1998). In this approach, the density of objects is determined in concentric annuli centered on each Spitzer-identified dusty YSO within the NICMOS, WFC3, and NSFCAM2 fields. The density is then averaged over all the dusty YSOs using the equation

$$\begin{eqnarray}{{\rm{\Sigma }}}_{\mathrm{comp}}({r}_{i}) & = & \displaystyle \frac{1}{{N}_{\mathrm{YSO}}}\displaystyle \sum _{j=1}^{{N}_{\mathrm{YSO}}}\\ & & \times \;\displaystyle \frac{{N}_{\mathrm{comp}}(j,{r}_{i}-{\rm{\Delta }}i/2.\lt r\leqslant {r}_{i}+{\rm{\Delta }}i/2)}{\pi [{({r}_{i}+{\rm{\Delta }}i/2)}^{2}-{({r}_{i}-{\rm{\Delta }}i/2)}^{2}]},\end{eqnarray} \tag{ 1 }$$

where the Spitzer-identified YSOs at the center of the annuli (hereafter: primaries) are numbered j = 1 to N_YSO, Δi is the width of the ith annuli that extends from ${r}_{i}-{\rm{\Delta }}i/2$ to ${r}_{i}+{\rm{\Delta }}i/2$, and N_comp is the number of NICMOS, WFC3, or NSFCAM2 sources within the ith annulus centered on the jth Spitzer YSO.

To determine N_comp within the annuli, we first identified a sample of primary sources by selecting the Spitzer-identified YSOs that appear as point sources in the HST data. This excludes sources that are too reddened to be detected in the F160W, F205W, and L'-band images and extended nebulous sources in which the central protostar is not detected as a point source. To find the number of sources as a function of radii from these primary sources, we use the point-source catalog extracted from the NICMOS, WFC3, and NSFCAM2 data.

In Figure 3, we plot the Σ_comp(r_i) for the protostars and pre-main-sequence star primaries separately. The pre-main-sequence primaries are only those that can be identified by their infrared excess and thus have disks. At larger separations, the plot shows a constant surface density of sources due to chance alignments with other YSOs and background stars. It also shows a clear spike in the density of sources within 1000 au of the primaries; we identify this spike as the signature of a population of companions. We note that there may also be companions at larger distances; however, the surface density of those companions is not high enough to rise above the surface density of the line-of-sight contamination. Sources at separations <80 au cannot be resolved; therefore, we are only able to identify companions between 80 and 1000 au.

Zoom In Zoom Out Reset image size

To further confirm this result, we also show the results of the Emark routine from the R package (Baddeley & Turner 2003; Schlather et al. 2004). Emark performs an analysis of a marked-point process where each dusty YSO identified in the Spitzer catalog is assigned a mark of 1 and the remaining sources a mark of 0. The routine determines the expectation value of the marks given the presence of an additional source at a specified projected separation. The increase of the expectation value at small separations indicates that sources with a neighboring source within 1000 au are more likely to be YSOs, once again demonstrating that there is an enhanced density of sources around the YSOs.

Between 80 and 1000 au, we find that there are 29 candidate binary and one candidate ternary systems around the protostellar primaries and 27 candidate binary and one candidate ternary systems around the pre-main-sequence primaries. The primaries and companions are displayed in Table 1 with their positions, separations, and HST or IRTF magnitudes. Three of the companions around protostars are at separations less than 100 au (HOPS 281, 268, and 24 at separations of 80, 96, and 99.9 au, respectively; see Table 1), and one of the binaries around pre-main-sequence stars has a separation less than 100 au (MGM 3376, separation 95 au). All of the detected YSOs that do not have a companion are listed in Table 2.

We present the numbers of primaries and candidate companions between 100 and 1000 au in Table 3; we will further motivate these limits in the following subsection. The numbers are determined for three different permutations of our combined IRTF, NICMOS, and WFC3 data set: the combined set from all three cameras, the NICMOS and WFC3 data, and the WFC3 data only. While the combined data set from all three cameras (hereafter the combined sample) gives us the most companions, the IRTF angular resolution is lower than that of the WFC3 and NICMOS data, motivating a combined look at the data from the two HST cameras alone (hereafter the HST sample). Finally, the WFC3 data set has the highest sensitivity and has the largest FOV, and hence it is worthwhile to look at the data from this camera alone (hereafter the WFC3 sample). All but one of the pre-main-sequence stars with disks are found only in the WFC3 sample; the remaining one is in the IRTF sample. The images of the candidate multiples obtained with the three cameras are shown in Figure 4 through 11.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3.  Number of Primaries, Companions, and Contaminants vs. Sample

Sample	Protostars	Cand.a	Cont.b	Comp.c	Pre-msd	Cand.a	Cont.b	Comp.c
Combined	201	28	152	20.8${}_{-2.8}^{+3.2}$	198	28	168	20.1${}_{-3.1}^{+2.9}$
HST	178	27	142	${20.3}_{-2.3}^{+2.7}$	197	28	168	${20.1}_{-3.1}^{+2.9}$
WFC3	129	21	109	15.9${}_{-1.9}^{+2.1}$	197	28	168	20.1${}_{-3.1}^{+2.9}$

Notes.

^aNumber of sources between projected separations of 100 and 1000 au. ^bNumber of sources between projected separations of 2000 and 5000 au. ^cNumber of companions calculated using r_inner = 100 au. ^dPre-main-sequence stars identified by the presence of an IR excess due to a dusty disk.

Download table as:  ASCIITypeset image

The detection of faint companions can be limited by confusion with the structure in the wings of the primary's PSF. To determine the detection limits for companions as a function of projected separation, we conducted a fake star analysis on the WFC3 data. We created a PSF using the IDL implementation of the DAOPHOT package utilizing a densely populated field in the WFC3 data. The PSF was then added to the vicinity of four isolated stars in different WFC3 fields. The PSFs were added at different projected separations, ranging from 45 to 1000 au with a 5 au increment. At each separation, the PSF was added to 12 positions. This process was repeated with fake stars that were 0–8 mag fainter than the four stars, in steps of 0.5 mag. The fake stars were then recovered with PhotVis using the same search parameters as for the WFC3 data. The fraction of sources identified for each combination of projected separations and Δm₁₆₀ was then determined, where Δm₁₆₀ is the fake star magnitude minus the primary magnitude. The results at each separation were averaged over concentric annuli with a uniform spacing in separation. Since the magnitudes of the faintest companions we can detect depend on the brightness of the primary, the detection fractions are also averaged over uniform intervals of Δm₁₆₀.

The curves of the projected separations versus Δm₁₆₀ determined from this analysis are shown in Figure 12. Overlaid are the separations and Δm₁₆₀ for the the primaries and all the candidate companions. In addition, double stars found around the remaining stars in the observed fields are shown: these may be visual binaries around pre-main-sequence stars without disks, or they may be optical doubles due to random coincidences between stars in the line of sight. In either case, they demonstrate that the discovery space delimited by the fake star test is consistent with the observed visual binaries and optical double stars.

Zoom In Zoom Out Reset image size

The 99% completeness limits for the WFC3 data range from Δm₁₆₀ ∼ 2.5 mag at 200 au, to 5 mag at 500 au, and to almost 7 mag at 1000 au. The resulting mass limits depend on the brightness of the primary, the amount of extinction, and the age of the source. Given that the last two cannot be ascertained directly from these data, we cannot uniquely determine the mass ratio probed. However, a Δm₁₆₀ range of 2.5 mag translates into a mass range of approximately a factor of 3–5, a range of 5 mag implies a mass range covering approximately a factor of 15, and 7 mag entails a mass range in excess of a factor of 50 (using the pre-main-sequence tracks of Baraffe et al. 1998).

To determine if the completeness values are comparable for NICMOS, we repeated the analysis on the NICMOS data for a smaller number of separations and mass ratios. Unlike the WFC3 data, the sources in the NICMOS data were identified by visual inspection. The analysis showed that the completeness of NICMOS was similar to that of the WFC3. A more detailed analysis will be done as part of a future paper of the separation and Δmag distributions. In contrast, the magnitude range and separations probed by the IRTF data are much more limited; these data were primarily obtained to provide color information for the companions. For this reason, we do not include the sources only imaged by the IRTF in the following analysis. We also limit our analysis to the candidate companions with projected separation of 100–1000 au. Since the densities of sources at larger radii are dominated by either field stars or other YSOs that happen to be found in projection near the targeted source, we cannot look for companions at wider separations with these data. For separations less than 100 au, we are only sensitive to mass ratios close to unity and thus highly incomplete.

The number of companions must be corrected for chance coincidences with stars in the line of sight. To do this, we estimate the surface density of the contaminants by counting the number of sources detected in the annulus between 2000 and 5000 au and summing over all the primaries. At these separations, the mean surface density of companions is flat and does not increase with decreasing separation, indicating that we no longer detect signifiant numbers of companions (Figure 3). The density is measured near enough to the primary that it provides an accurate assessment of the surface-density contaminants. We then calculate the number of companions from the equation

$$\begin{eqnarray}&&{N}_{\mathrm{comp}}={N}_{\mathrm{cand}}-{N}_{\mathrm{cont}}\left(\displaystyle \frac{{{\rm{\Omega }}}_{\mathrm{comp}}}{{{\rm{\Omega }}}_{\mathrm{cont}}}\right),\end{eqnarray} \tag{ 2 }$$

where N_cand is the number of candidate companions with projected separations from their primaries between 100 and 1000 au, and N_cont is the number of sources in the region between 2000 and 5000 au.

We assume that N_cont gives the number of contaminating stars, either foreground or background stars or other Orion YSOs in the line of sight, which we scale by the ratios of the solid angles of the companion and contamination regions to estimate the number of contaminating objects in the companion region. If we ignore the effect of incompleteness, the ratio of solid angles becomes the area of the annulus extending from 100 to 1000 au for the companions over the area of the annulus extending from 2000 to 5000 au for the contaminants:

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{\Omega }}}_{\mathrm{comp}}}{{{\rm{\Omega }}}_{\mathrm{cont}}}=\left(\displaystyle \frac{{1000}^{2}-{100}^{2}}{{5000}^{2}-{2000}^{2}}\right)=0.047.\end{eqnarray} \tag{ 3 }$$

The number of contaminants in the 2000–5000 au annulus and the corrected number of companions are given in Table 3 for each of the three samples. The contamination is the total number of contaminants for all of the primaries. The number of contaminants is small because of the high extinction toward the YSOs: there are 0.04 companions per primary between 100 and 1000 au. Given the small number of companions, however, as many as 26% of the companions may be contaminants. Thus, the primary uncertainty in the number of companions is due to the contamination subtraction. In Appendix A , we discuss how we measure the uncertainty in the contamination subtraction using the Poisson statistics of the line-of-sight contamination.

Since the reduction of completeness near the primary can also limit the detection of contaminants, the number of contaminants at 100–1000 au can be overestimated. Due to confusion with the primary's PSF, fainter stars can be detected in the 2000–5000 au annulus, where we measure the surface density of contaminants, than can be detected within 1000 au of the primary. This will result in an oversubtraction in the number of contaminants, and therefore the numbers of companions in Table 3 are lower limits to the actual numbers of companions. The oversubtraction will be partially mitigated by the uniform spatial distribution of contaminants, which implies that most of the contamination will be at the largest separations that subtend the largest angle on the sky.

To correct for this bias in the WFC3 sample, we repeat our estimates of the number of companions and the CFs. We use the 25% and 99% completeness curves for the WFC3 data to determine the inner radii for calculating Ω_comp (Figure 12). The correction is calculated by determining Δm₁₆₀ between each object in the 2000–5000 au annulus and the targeted primary star in the center of the annulus and then adopting the criteria that contaminating objects in the line of sight can only be detected outside a limiting radius where the completeness is ≥99% for that value of Δm₁₆₀. We determine the ratio with the equation

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{\Omega }}}_{\mathrm{comp}}}{{{\rm{\Omega }}}_{\mathrm{cont}}}=\displaystyle \frac{1}{{N}_{\mathrm{cont}}}\displaystyle \sum _{i=1}^{{N}_{\mathrm{cont}}}\;\left(\displaystyle \frac{{1000}^{2}-{r}_{i}^{2}}{{5000}^{2}-{2000}^{2}}\right),\end{eqnarray} \tag{ 4 }$$

where r_i is set to r_99% for the ith object in the 2000–5000 au annulus, down to a minimum value of 100 au. The number of contaminants is reduced by a factor of 0.46 for protostars and by a factor of 0.42 for pre-main-sequence stars. This analysis was repeated using the Δm₁₆₀ versus radius curve for a 25% level of completeness by adopting r_25% for the inner radius (Figure 12). This reduces the number of contaminants by a factor of 0.6 for protostars and 0.56 for pre-main-sequence stars. By reducing the estimated number of line-of-sight contaminants, the number of companions to protostars increases from 15.9 to 18 for r_25% and to 18.6 for r_99%, and the number of companions to pre-main-sequence stars increases from 20.1 to 23.6 and 24.7 for r_25% and r_99%, respectively. This is an increase of 1.35σ for the protostars and 1.5σ for the pre-main-sequence stars. The number of candidate companions and the level of contamination for the different choices of r_inner are given for the WFC3 sample in Table 4. Again, the determination of the uncertainties is described in Appendix A .

Table 4.  WFC3 Sample: Number of Primaries and Companions vs. Surface Density

ΣYSO	rinner	Protostars	Cand.a	Cont.b	Comp.c	Pre-msd	Cand.a	Cont.b	Comp.c
Alle	100 au	129	21	109	15.9${}_{-1.9}^{+2.1}$	197	28	168	20.1${}_{-3.1}^{+2.9}$
Alle	r25%	129	21	65.1	17.9${}_{-1.9}^{+2.1}$	197	28	94.1	23.6${}_{-2.6}^{+2.4}$
Alle	r99%	129	21	50.3	18.6${}_{-1.6}^{+1.4}$	197	28	70.8	24.7${}_{-1.7}^{+1.3}$
>45 pc−2	100 au	56	12	51	9.6${}_{-1.6}^{+1.4}$	123	19	93	14.6${}_{-2.6}^{+2.4}$
<45 pc−2	100 au	73	9	58	6.3${}_{-1.3}^{+1.7}$	74	9	75	5.5${}_{-1.5}^{+1.5}$
>45 pc−2	r25%	56	12	31.5	10.5${}_{-1.5}^{+1.5}$	123	19	52.2	16.5${}_{-1.5}^{+1.5}$
<45 pc−2	r25%	73	9	33.6	7.4${}_{-1.4}^{+1.6}$	74	9	41.9	7.0${}_{-1.0}^{+1.0}$
>45 pc−2	r99%	56	12	26.6	10.7${}_{-0.7}^{+1.3}$	123	19	39.6	17.1${}_{-1.1}^{+1.9}$
<45 pc−2	r99%	73	9	23.7	7.9${}_{-0.9}^{+1.1}$	74	9	31.2	7.5${}_{-1.5}^{+1.5}$

Notes.

^aNumber of sources between projected separations of 100 and 1000 au. ^bNumber of sources between 2000 and 5000 au. ^cNumber of companions calculated using r_inner. ^dPre-main-sequence stars identified by the presence of an IR excess due to a dusty disk. ^eIncludes both high and low YSO surface density regions.

Download table as:  ASCIITypeset image

The incidence of multiplicity among a population of stars can be quantified by the companion fraction, or CF, which is defined by the ratio of the number of companions to the total number of single stars and systems:

$$\begin{eqnarray}&&\mathrm{CF}=\displaystyle \frac{B+2T+3Q}{S+B+T+Q},\end{eqnarray} \tag{ 5 }$$

where S is the number of single-star systems, B the number of binary systems, T the number of ternary systems, and Q the number of quaternary systems (e.g., Ghez et al. 1997). (This equation can be expanded to include larger systems, but there are only binary and ternary systems in our sample.) The CF can be calculated by finding the total number of candidate companions with separations between 100 and 1000 au and then subtracting the estimated number of line-of-sight contaminants. Alternatively, we can think of the CF as the integral of the mean surface density of companions minus the average density of line-of-sight contaminants, over an integration range of 100–1000 au (Figure 3). Multiplicity fraction (MF) is also frequently used in the literature; this is given by the fraction of stars hosting multiple systems and counts each system only once. Since the CF is derived from the surface density of point sources from which we have identified the Orion companions, it is the natural statistic to use in our analysis. The CF is then

$$\begin{eqnarray}&&\mathrm{CF}=\displaystyle \frac{{N}_{\mathrm{comp}}}{{N}_{\mathrm{prim}}}-\displaystyle \frac{{N}_{\mathrm{cont}}}{{N}_{\mathrm{prim}}}\left(\displaystyle \frac{{{\rm{\Omega }}}_{\mathrm{comp}}}{{{\rm{\Omega }}}_{\mathrm{cont}}}\right),\end{eqnarray} \tag{ 6 }$$

where N_prim is the number of primaries. The calculation of the uncertainties for the CF is described in Appendix B .

The CFs are given in Table 5 for an r_inner = 100 au. The CF of all the dusty YSOs is ${10.3}_{-1.0}^{+1.0}\%$ for the combined sample, ${10.8}_{-0.9}^{+1.0}\%$ for the HST sample, and ${11}_{-1.2}^{+1.2}\%$ for the WFC3 sample. Thus, the CF does not change significantly between the different samples. For the HST sample, the CFs are ${11.4}_{-1.3}^{+1.5}\%$ for protostars and ${10.2}_{-1.6}^{+1.5}\%$ for the pre-main-sequence stars. Although the CF for pre-main-sequence stars is slightly lower, the CFs are consistent within the uncertainties. The CFs for different choices of r_inner are given for the WFC3 sample in Table 5. The CF increases from ${12.3}_{-1.4}^{+1.7}\%$ to ${13.9}_{-1.5}^{+1.6}\%$ and ${14.4}_{-1.3}^{+1.1}\%$ for the protostars, using r_25% and r_99%, respectively (Table 6). For the pre-main-sequence stars, the CF increases from ${10.2}_{-1.6}^{+1.5}\%$ to ${12}_{-1.3}^{+1.2}\%$ and ${12.5}_{-0.8}^{+1.2}\%$ for r_25% and r_99%, respectively (Table 6). The CFs for the pre-main-sequence stars and protostars both increase and remain consistent, and the absence of a significant difference between the two evolutionary stages persists.

Table 5.  CFs vs. Sample

Sample	Protostarsa	Pre-msa,b	Mergeda,c
Combined	10.4${}_{-1.4}^{+1.6}$	${10.1}_{-1.6}^{+1.5}$	${10.3}_{-1.0}^{+1.0}$
HST	11.4${}_{-1.3}^{+1.5}$	${10.2}_{-1.6}^{+1.5}$	${10.8}_{-0.9}^{+1.0}$
WFC3	12.3${}_{-1.4}^{+1.7}$	${10.2}_{-1.6}^{+1.5}$	${11.0}_{-1.2}^{+1.2}$

Notes.

^aCalculated using an r_inner = 100 au. ^bPre-main-sequence stars with disks. ^cMerged sample of dusty YSOs.

Download table as:  ASCIITypeset image

Table 6.  WFC3 Sample: CF vs. Surface Density

${{\rm{\Sigma }}}_{\mathrm{YSO}}$	rinner	Protostars	Pre-ms Stars.a	Merged.b
Allc	100 au	${12.3}_{-1.4}^{+1.7}\%$	${10.2}_{-1.6}^{+1.5}\%$	${11.0}_{-1.2}^{+1.2}\%$
Allc	r25%	${13.9}_{-1.5}^{+1.6}\%$	${12.0}_{-1.3}^{+1.2}\%$	${12.7}_{-0.8}^{+0.8}\%$
Allc	r99%	${14.4}_{-1.3}^{+1.1}\%$	${12.5}_{-0.8}^{+1.2}\%$	${13.3}_{-0.7}^{+0.8}\%$
>45 pc−2	100 au	${17.1}_{-2.8}^{+2.5}\%$	${11.9}_{-1.3}^{+1.9}\%$	${13.5}_{-1.2}^{+1.6}\%$
<45 pc−2	100 au	8.6${}_{-1.7}^{+2.4}\%$	7.4${}_{-2.0}^{+2.1}\%$	8.0${}_{-1.9}^{+1.5}\%$
>45 pc−2	r25%	${18.8}_{-2.7}^{+2.6}\%$	${13.4}_{-1.3}^{+1.2}\%$	${15.1}_{-1.1}^{+1.1}\%$
<45 pc−2	r25%	${10.2}_{-1.9}^{+2.2}\%$	9.5${}_{-1.4}^{+1.3}\%$	9.8${}_{-1.0}^{+1.1}\%$
>45 pc−2	r99%	${19.2}_{-1.3}^{+2.2}\%$	${13.9}_{-0.9}^{+0.7}\%$	${15.6}_{-1.0}^{+1.2}\%$
<45 pc−2	r99%	${10.8}_{-1.2}^{+1.5}\%$	${10.2}_{-2.1}^{+2.0}\%$	${10.5}_{-1.0}^{+1.1}\%$

Notes.

^aPre-main-sequence stars identified by the presence of an IR excess due to a dusty disk. ^bIncludes both protostars and pre-main-sequence stars with disks. ^cIncludes both high and low YSO surface density regions.

Download table as:  ASCIITypeset image

A comparison of the CF in regions of high and low stellar density is motivated by the observed variations in the protostellar luminosity function with stellar density discovered by Kryukova et al. (2012). Although the changes in the protostellar luminosity function may not be directly related to changes in the multiplicity, the variations in the luminosity function indicate that the low-mass star-formation process is altered by the environments found in regions of low and high stellar density. To compare the incidence of multiplicity in regions of low and high stellar density, we divide our sample of YSOs by adopting an analysis similar to that of Kryukova et al. (2012). Specifically, we used the Spitzer Orion Survey point-source catalog of Megeath et al. (2012) to determine the projected distance to the fourth-nearest Spitzer YSO, hereafter NN5, for each primary or single star in the HST survey.¹¹ The fourth-nearest neighbor was chosen since it is not biased by binary, ternary, or quaternary systems and yet still measures the local density of YSOs surrounding a particular object. When computing NN5, we do not include any of the sources that were identified as likely companions since none of them have been resolved in the Spitzer sample. We then use the NN5 to divide the YSOs into high- and low-density regions by whether or not their nearest-neighbor density, as given by $3/(\pi {r}_{\mathrm{NN}5}^{2})$, exceeds the threshold surface density (see Megeath et al. 2016 for the definition of the nearest-neighbor density adopted in this paper).

In Figure 13, we plot the CFs of the HST sample for the high and low surface density regions for threshold densities corresponding to NN5 = 25,000 au (65 pc⁻²), 30,000 au (45 pc⁻²), and 35,000 au (33 pc⁻²); these values bracket the NN5 where the number of objects in high- and low-density regions are equal. For each of these, we include error bars for the line-of-sight contamination. In this figure, we set r_inner = 100 au, which implies there may be an oversubtraction of the contamination. The subtraction of the line-of-sight contamination ensures that the variations in the CF are not due to a higher number of chance coincidences with other YSOs in regions of higher density. We find that, for the range of threshold NN5 values, the protostars and young stars with disks found in dense regions have a higher CF than the low-density regions. Note that the CFs for the protostars, pre-main-sequence stars with disks, and combined dusty YSOs in the high-density regions are consistent with each other within the uncertainties, as are the CFs of the protostars, disks, and dusty YSOs in the low-density regions.

Zoom In Zoom Out Reset image size

A more rigorous analysis can be done with the WFC3 data using different values of r_inner. In Table 6, we show the number of primaries and companions in the low- and high-density regions as a function of r_inner. The resulting ratios,

$$\begin{eqnarray}&&R=\displaystyle \frac{\mathrm{CF}({\rm{\Sigma }}\gt 45\;{{\rm{pc}}}^{-2})}{\mathrm{CF}({\rm{\Sigma }}\lt 45\;{{\rm{pc}}}^{-2})},\end{eqnarray} \tag{ 7 }$$

are calculated using a Bayesian parameter estimation in Appendix B and presented in Table 7. Since the ratio of CFs does not appear to be strongly dependent on the threshold NN5 value, we use the 45 pc⁻² density threshold only. The spatial distribution of YSOs in our sample based on this threshold is presented in Figure 14. As we adjust the value of r_inner from 100 au to r_99%, the median value of R decreases from 1.96 to 1.74 for protostars, from 1.5 to 1.28 for pre-main-sequence stars, and from 1.64 to 1.45 for the merged samples. Each of these three samples consistently shows an elevated CF in the dense environments, although the change in CF is highest for protostars.

Zoom In Zoom Out Reset image size

The primary uncertainty in the CF is Poisson fluctuations in the amount of contamination (Appendix A ); hence, we estimate the significance of this enhancement with respect to these fluctuations using three methods. First, our parameter estimation determines the probability distribution of R (Table 7, Appendix B ). These show that the probability for R ≤ 1 depends on both R_inner and the sample. For protostars, the probability decreases from 0.01 to 0.001 as we change r_inner from 100 au to r_99%. Similarly, for stars with disks, the probability decreases from 0.08 to 0.01 as we change r_inner. Finally, for the merged sample, the probability declines from 0.008 to 0.0004. Parameter estimation techniques, however, are not always appropriate for hypothesis testing. Thus, we also employ a Bayesian hypothesis test to determine the probability of R ≤ 1 (Appendix B ). We find for protostars that the probability of R ≤ 1 ranges from 0.05 to 0.007 as we change r_inner from 100 au to r_99%, for pre-main-sequence stars the probability ranges from 0.22 to 0.008, and for the merged sample, the probability ranges from 0.036 to 0.004. Hence, the significance of the result remains. Finally, we employ a frequentist hypothesis test to determine the probability of the null hypothesis, that the CF is the same in both environments (Appendix B ). The probability that R = 1 changes ranges from 0.008 to 0.0007 for protostars as we change r_inner from 100 au to r_99%, from 0.07 to 0.04 for disks, and from 0.009 to 0.0008 for the merged sample. We conclude that the variation in the CF between high- and low-density regions is significant and is highly unlikely to result from fluctuations in the contamination.

The uncertainties in the number of companions and CFs in Tables 3–6 are the uncertainties from the contamination correction: no additional sources of uncertainty are included. In the literature, the number of companions typically includes an additional uncertainty typically given by $\sqrt{{{PN}}_{\mathrm{prim}}}=\sqrt{{N}_{\mathrm{comp}}}$, or more correctly by $\sqrt{(1-P){{PN}}_{\mathrm{prim}}}$ (Burgasser et al. 2003). This additional uncertainty comes from modeling the presence of a companion as the result of a Bernoulli trial. In this model, the number of companions is described by the binomial distribution with the number of trials equal to the number of primaries, the number of successful trials equal to the number of companions, and P being the probability of having a companion. Modeling the presence of a companion as a Bernoulli trial has two possible physical interpretations. First, the formation of a companion can be approximated as a random process, which has a probability P of creating a companion between 100 and 1000 au. Alternatively, if the observed primaries are drawn from a much larger population of primaries, in which a fraction P have companions between 100 and 1000 au, then the choice of a primary with a companion could once again be treated as a Bernoulli trial. In both cases, the uncertainty in the number of companions, ignoring the uncertainty from contamination subtraction, is given by $\sqrt{(1-P){{PN}}_{\mathrm{prim}}}$ (Burgasser et al. 2003).

The modeling of the presence of a binary as a Bernoulli trial may not always be appropriate. In the first interpretation, the presence of a companion is assumed to be the result of a random process with probability P; this may not be true if the presence of a companion is strongly influenced by external, environmental factors. In the second interpretation, the observed sample is assumed to be drawn from a much larger population. This may be relevant for binary surveys of field stars; however, since we observed most of the protostars in the Orion molecular clouds outside the ONC, the assumption that we are drawing primaries randomly from a much larger population is not appropriate for our analysis of Orion protostars. It is, however, appropriate for our sample of pre-main-sequence stars that is drawn from a much larger population.

For these reasons, we have not included binomial statistics in our analysis of the CF and the ratio of the CFs in different environments. This motivates the question whether the difference in the CF would be significant if binomial statistics were included. To address this question, we define the companion star probability, or CP, as the probability of a star having a companion. We assume the presence of a companion can be treated as a random process similar to a Bernoulli trial, but with the additional outcomes of a ternary and quaternary system. The CP can be written as

$$\begin{eqnarray}&&\mathrm{CP}={P}_{\mathrm{binary}}+2{P}_{\mathrm{ternary}}+3{P}_{\mathrm{quatenary}},\end{eqnarray} \tag{ 8 }$$

where P_binary is the probability of a binary system, P_ternary is the probability of a ternary system, and P_quaternary is the probability of a quaternary system. (Additional terms can be added for larger systems.) The average number of companions is given by

$$\begin{eqnarray}&&{N}_{\mathrm{comp}}=\mathrm{CP}\times {N}_{\mathrm{prim}}.\end{eqnarray} \tag{ 9 }$$

If we adopt a particular value for the CP, the resulting CFs follow a binomial distribution with a probability equal to the CP.

To examine whether the CP varies between the high- and low-density samples, we repeat the analysis applied to the CF for the CP (Appendix C ). First, we use parameter estimation to estimate the value of R_CP by substituting CP for CF in Equation (7). In Table 7, we find that the median value of R_CP is within 4% of the value given by the R determined with the CFs. However, the significance of the result is lower: the probability that R_CP ≤ 1 is 0.11–0.10 for protostars, 0.22–0.26 for disks, and 0.10–0.11 for the merged sample.

If we apply a Bayesian hypothesis test to the CP, we must consider the probability that CP_high > CP_low, that CP_high = CP_low, and that CP_high < CP_low (in contrast, CF_high cannot equal CF_low since the number of contaminating stars can only be an integer and the resulting CFs of high- and low-density regions cannot be exactly the same; see Appendix B ). Since the R_CP = 1 hypothesis has half the number of parameters (i.e., a single value of CP versus independent values of CP for the high- and low-density regions), the resulting Ockham factor favors the R_CP = 1 hypothesis: the probability that R_CP = 1 is 0.48–0.58 for protostars, 0.60–0.63 for disks, and 0.58–0.61 for the merged sample. These assume that the prior probability for CP is uniformly distributed between 0 and 1. Finally, we can use our frequentist hypothesis test (Appendix C ) to find the probability that R_CP = 1. This probability is 0.10–0.12 for protostars, 0.17–0.20 for disks, and 0.07–0.09 for the merged sample. Thus, despite having one of the largest surveys for multiplicity to date, a larger sample is needed. We note that, if we used the combined HST sample and r_inner = 100 au, the frequentist hypothesis test for the merged sample gives a probability as low as 0.03 that R_CP = 1, suggesting that the differences in the CP become more significant as we increase the sample size. A future investigation will employ the complete HST data, taking into account the different sensitivity levels, to reexamine variations in the CP and its significance.

The analysis of the mean frequency of companions for pre-main-sequence stars was pioneered by Larson (1995) using existing catalogs of the Taurus–Auriga region. Over spatial scales of 5 au to 5 pc, Larson found that the mean surface density of companions could be fit by two power laws that intersected at a separation of 0.04 pc. The inner power law fit the data for projected separations less than 0.04 pc and had an exponent of −2.15. The sources included in this part of the power law were considered to be part of multiple systems. This power law is in approximate agreement with the flat distribution of orbital periods per logarithmic interval of period for solar-type stars in the solar neighborhood (Duquennoy & Mayor 1991; Bate et al. 1998). In contrast, the power law for separations greater than 0.04 pc was much shallower, with an exponent of −0.62. Larson (1995) proposed that this power law reflected the fractal spatial distribution of young stars in the Taurus cloud as a whole. In this interpretation, sources within 0.04 pc were true companions that were formed within one Jeans length of each other, while objects at larger separations were young stars in the Taurus–Auriga region that were not gravitationally bound companions. Using the HST/NICMOS observation of the Lynds 1688 region of the Ophiuchus molecular cloud, Allen et al. (2002) found a similar dual power law structure for the mean surface density of companions in that cloud. Simon (1997) compiled data sets on three star-forming regions and fit dual power laws to the mean surface densities of the Taurus, Ophiuchus, and Orion clouds. The power laws were similar to those found by Larson (1995), but the break point between the power laws varied from 12,000 au in Taurus to 400 au in Orion. Simon (1997) proposed that the separation was dependent on the Jeans mass and the stellar density. Bate et al. (1998) clarified the dependence on stellar density by showing that the break point occurs at the point at which the surface density of companions exceeds that of young stars in the line of sight. This break point occurs at shorter and shorter separations as the surface density of line-of-sight stars increases. Because we are comparing column densities, the break point depends not only on the volume density of young stars but also on the depth of the observed star-forming region and the corresponding depth of the column over which the surface density is measured. This length is much larger than the orbital separations in multiple-star systems, and the break point can occur at a separation much smaller than the maximum orbital separations of gravitationally bound systems.

To compare our data to these previous studies, we merge our HST with the Spitzer Orion Survey of Megeath et al. (2012) to look at the mean surface density of companions between 100 au and 5 pc. For separations ≤1000 au, these sources come from the HST data. For separations ≥10,000 au, these sources are the dusty YSOs from the Spitzer survey. In the in-between region of overlap, we eliminated sources that were likely faint background objects and instituted a F160W magnitude cutoff of 20 mag. The mean surface density of companions from this combined data set is shown in Figure 17. As found previously, the mean surface density of companions shows a steep power law at small separations and a much shallower power law at large separations. However, there are several differences. First, we find a plateau in the mean surface density that extends from 1000 to 20,000 au (i.e., from 25 to 50''). Since the extent of the plateau corresponds to approximately half the size of the WFC3 field (135'' × 127''), we find it likely that the plateau is an observational artifact. This artifact arises from the inclusion of HST sources for which there is not sufficient Spitzer photometry to distinguish between a dusty YSO and a background source. Since these sources could be faint background stars and galaxies, many of the objects in the plateau are likely to be background objects detected by WFC3. Second, we find that the break point between the inner power law and the plateau is at 1000 au, or 0.005 pc, much less than the 0.04 pc found by Larson (1995). As explained by Bate et al. (1998), this point is set by where the surface density of point sources in the line of sight equals the surface density of companions. In the case of the Orion data, this occurs when the surface density of companions equals the density of sources in the plateau. Hence, the smaller distance in Orion is not due to a physical difference but to a higher level of line-of-sight contamination.

Zoom In Zoom Out Reset image size

A final difference is that the inner power law is less steep; the slopes of the inner power laws from the HST Orion survey are −1.43 ± 0.79 and −1.49 ± 0.44 for the protostars and pre-main-sequence stars, respectively, compared to −2.15 for Taurus (Larson 1995). This may be due to incompleteness at short separations; however, the uncertainties are large, and the difference is of a marginal statistical significance. Finally, the outer power law is steeper than those found in other regions. Part of this may result from the plateau, which has forced us to fit the outer power law over a different range of separations than that used by Larson (1995) or Simon (1997). Alternatively, the difference may be due to structural differences in the distribution of young stars in Orion compared to Taurus or Ophiuchus. Nevertheless, our result is qualitatively consistent with the interpretations of previous authors, reinforcing our claim that the mean surface density of companions within 1000 au is tracing companions and therefore can be used to characterize multiplicity in the Orion molecular clouds.

We find no evidence for a change in the CF between the protostellar and pre-main-sequence systems. Our CF is also consistent with that found by the HST/ACS survey of pre-main-sequence stars in the ONC: between 100 and 660 au, both Reipurth et al. (2007) and this study find a CF of 0.09. Furthermore, the CF for stars between 100 and 1000 au for solar-type (G) stars is 0.13 (Duquennoy & Mayor 1991).¹² The fact that the CF for the Orion protostars and pre-main-sequence stars is similar to that for solar-type field stars may result from the observed YSOs in the Orion sample having masses lower than solar-type stars, which would be expected for a standard IMF. Since main-sequence stars at 0.5 M_⊙ have CFs lower than that for G stars (Duchêne & Kraus 2013), the dusty YSO CFs may be larger than the field star CF for the appropriate mass range, and the agreement between the dusty YSO CFs and that of solar-type stars may be spurious. However, given the uncertainty in the masses of the Orion YSOs, our observed CFs are consistent with those of field stars, and we find no clear evidence for evolution.

As mentioned in the introduction, previous studies provided evidence for a large decrease in multiplicity from the Class 0 protostellar phase to the main sequence. This change is observed at the >50 au separations that can be resolved in submillimeter and IR observations. In an analysis of Submillimeter Array (SMA) data toward 33 Class 0 protostars, Chen et al. (2013) find an MF and CF of 0.64 ± 0.08 and 0.91 ± 0.05, respectively, in various star-forming regions, with a median resolution of 25 or ∼600 au for a typical distance of 240 pc. For a separation range of 100–1000 au, the CF is 0.24 ± 0.09; since their survey is incomplete at these separations, this CF should be considered a lower limit. Near-IR surveys have found the CF for more evolved Class I protostars. Duchêne et al. (2004) determine the CF between 110 and 1400 au to be 0.23 ± 0.09 for the Taurus Class I protostars and 0.29 ± 0.07 for similar protostars in Ophiuchus. Subsequent higher angular resolution observations gave CFs of 0.32 ± 0.06 between 45 and 1400 au in Taurus and 0.47 ± 0.08 between 14 and 1400 in Ophiuchus for a combined sample of 189 protostars (Duchêne et al. 2007). If we limit the separation range to 100–1000 au, the CFs are 0.14 ± 0.08 in Taurus and 0.23 ± 0.10 in Ophiuchus. Connelley et al. (2008b) find a CF across various star-forming regions to be 0.43 between 100 and 4500 au; this reduces to ∼0.24 if we limit the separations to 100–1000 au. While some of these values are higher than those we find in Orion over a comparable separation range, given the different completeness limits for various surveys and the large uncertainties due to the relatively small sample sizes, it is difficult to draw a meaningful comparison between them.

There are several potential reasons why the Orion data do not exhibit a changing CF between protostars, pre-main-sequence stars, and the field. First, the masses of the protostars, pre-main-sequence stars, and their companions have not been determined; hence, we cannot directly compare our Orion CF to the main-sequence CF, which is a strong function of mass (Lada 2006; Duchêne & Kraus 2013). Furthermore, when comparing protostars and pre-main-sequence stars in Orion, the extinction from protostellar envelopes limits the detection of faint, very low mass companions. The faint magnitudes of the protostars and the companions in Figure 15 clearly demonstrate the effect of extinction on protostellar companions, although this is compensated for, at least in part, by the high sensitivity of the HST at these wavelengths. This extinction could reduce the CF of the protostars relative to that for pre-main-sequence stars, thereby masking a drop in the CF. Third, by selecting protostars that appear as a point source at 1.60 μm, our sample is biased toward objects at the later stages of protostellar evolution. Connelley et al. (2008b) find a decrease in the CF with spectral index, with flat-spectrum sources showing a distinctly lower CF. Thus, if the evolution of multiplicity occurs during the protostellar phase (e.g., Reipurth 2000; Reipurth et al. 2010), the lack of evolution apparent in our Orion sample could result from a sample biased toward protostars around which this evolution has already occurred. Fourth, the primary variation in the CF may be a drop at separations less than 100 au or greater than 1000 au (Connelley et al. 2008b); there is currently no clear evidence in any survey of star-forming regions of a drop in the CF(s) or MFs between 100 and 1000 au. Finally, as we will discuss in the next section, differences in the CF may also have an environmental component. King et al. (2012) find that pre-main-sequence stars in Taurus have a higher CF than other regions. Thus, studies of the evolution of the CF between protostars and field stars must take into account such environmental variations.

Given the potential biases listed above, our data do not rule out that there is a decrease in the CF with evolution. Upcoming observations with the VLA and ALMA, as well as observations at wavelengths longer than 1.6 μm with the James Webb Space Telescope (JWST) or ground-based adaptive optics, may help resolve this issue by detecting companions around protostars at even earlier stages of their development.

Our survey of Orion shows a dependence of the incidence of 100–1000 au multiplicity on environment as defined by the local YSO density. In Section 4.5, we found that the CF increases by a factor of ∼1.5 between environments of low stellar density and high stellar density. This increase occurs for both the protostar and the pre-main-sequence star samples. This is inconsistent with previous surveys, which either found no differences between high and low stellar density regions or found that high stellar density regions had systematically lower CFs. Duchêne et al. (2004) compared the multiplicity of Class I and flat-spectrum protostars in the Taurus Molecular Cloud, which is characterized by low stellar densities, to that of protostars in the Ophiuchus cloud, which hosts dense clustered star formation. Over a separation range of 110–1400 au, they found no significant difference in the CFs or MFs. Comparisons of the CFs and MFs of pre-main-sequence stars have shown a different picture. King et al. (2012) find that the MF of the Taurus region is higher than in the ONC, but that the Chameleon I region, IC348 cluster, and Ophiuchus cluster are within ∼1σ of Orion. In their review, Duchêne & Kraus (2013) highlight this difference by showing that the CFs per decade of separation are 0.20 and 0.08 in low-density associations and dense clusters, respectively. Interestingly, a survey of protostars in the Ophiuchus, Taurus, Serpens, and the Orion clouds by Connelley et al. (2008b) showed that the CF in Orion was different from that for the other clouds. Specifically, they found that Orion has a much higher CF at separations between 100 and 200 au than the CF for the combined sample of protostars in the other clouds. The lack of consistency in these results motivates surveys within a single complex where biases due to different distances and sensitivities can be minimized.

As discussed in the introduction, the observed variations between low- and high-density regions have been primarily explained as the effect of the cluster environment on the evolution of pre-main-sequence systems after formation (King et al. 2012; Marks & Kroupa 2012). There are several reasons why variations in the CFs we observed in Orion are unlikely to be the result of interactions between stars in a dense, clustered environment. First, the spatial distribution of the YSOs found in high and low spatial density regions within the Orion cloud do not show that high-density sources are concentrated in large clusters; instead, we find that regions of high YSO density are located throughout the Orion complex and are found both in large clusters and in small groups (Figure 14). Second, if the variations are due to the interaction with neighboring stars during flybys, we would expect a lower CF between 100 and 1000 au in dense regions (Marks & Kroupa 2012; Parker & Meyer 2014). Such a decrease is apparent in the HST/ACS observations of the ONC that show a deficiency of wide systems in the dense inner region of Orion (Reipurth et al. 2007). This is the exact opposite of what we find. Finally, if the changes of the CF resulted from interactions between YSOs in dense environments, the effect should be stronger for the older pre-main-sequence stars. This is not seen: disks exhibit a smaller change in the CF between high- and low-density environments. For these reasons, we argue that the increased CF in dense regions is due to the enhanced formation of companions between 100 and 1000 au. These companions must be formed by the time the central protostar is detectable by 1.60 μm imaging.

One possible reason why YSOs in high stellar density regions have a higher incidence of multiplicity is that they are systematically more massive. Observations of field stars show a strong increase of CF with mass (Duchêne & Kraus 2013). In the Orion clouds, Kryukova et al. (2012) found the luminosity function for protostars in regions of high stellar density was biased to higher luminosities than the luminosity function for regions of low stellar density. To assess the effect this may have on our analysis, we compare the observed magnitudes of the systems (which are usually dominated by the primaries; Figure 15) in high and low stellar density regions of the Orion clouds. We do this in two ways: by comparing the 24 μm magnitudes of the combined sample, and by comparing the dereddened J-band magnitudes of the pre-main-sequence stars.

We choose the 24 μm magnitude since both pre-main-sequence stars with disks and protostars show strong emission at these wavelengths. In comparison, observations of pre-main-sequence stars are generally not available at longer wavelengths, while shorter wavelength observations can be affected by the extinction and scattering in protostellar envelopes. Furthermore, this wavelength is usually not strongly affected by extinction. Thus, although there is not a simple conversion between mass and the 24 μm magnitude, it provides the means to compare the two samples to look for systematic differences in their observable properties. As shown in Figure 18, we do not see a systematic difference in the m₂₄ magnitude distributions of the high- and low-density regions. We determine the m₂₄ distributions of the low and high stellar density samples using the three different critical separations in Section 4.5; Kolmogorov–Smirnov (KS) tests comparing the two distributions yield probabilities of 0.46, 0.60, and 0.22 that the two are drawn from the same parent distribution for critical NN5 thresholds of 25,000 au, 30,000 au, and 35,000 au, respectively.

Zoom In Zoom Out Reset image size

A second approach is to compare the dereddened J-band magnitudes of the pre-main-sequence stars. The J-band magnitudes of pre-main-sequence stars are dominated by their photospheric emission and are strongly dependent on the masses of the stars. We use the 2MASS near-IR colors to deredden the sources following the approach of Gutermuth et al. (2008) and using the extinction law of Flaherty et al. (2007) as modified for Orion by Megeath et al. (2012). This gives us the photospheric J-band magnitudes of the pre-main-sequence stars. This magnitude depends on both the mass and age of the sources as they descend on pre-main-sequence tracks. We expect that different distributions of stellar mass would be reflected in different dereddened m_J distributions, except in the unlikely case that differences in the luminosity due to different mass distributions are canceled out by different distributions of ages. Again, as shown in Figure 19, the distributions are very similar. We compare the dereddened J-band magnitude distributions for sources in high and low stellar density regions using the three different NN5 thresholds in Section 4.5. In this case, KS tests give probabilities of 0.997, 0.96, and 0.61 that they are drawn from the same parent distribution for thresholds of 25,000 au, 30,000 au, and 35,000 au, respectively. We conclude that the environmental variations in the CFs are not likely to be the result of systematically different masses in our two samples.

Zoom In Zoom Out Reset image size

For these reasons, we find it likely that the variations in the CF are primordial. It is not understood, however, why YSOs would form more companions in regions of high stellar density. Part of the reason may have to do with the correlation between the gas column density and the stellar surface density (Gutermuth et al. 2011; Lada et al. 2013; Lombardi et al. 2013). Regions with a higher density of YSOs also tend to have a higher column density of molecular gas. This may also imply a higher volume density of gas leading to shorter Jeans lengths, smaller cores, and more prompt or turbulent fragmentation on short physical scales that lead to multiple systems (see also Tohline 2002; Connelley et al. 2008b; Offner et al. 2010). Alternatively, protostars in regions of high gas density may be subject to higher infall rates; this may lead to the formation of a larger number of companions in the outer, gravitationally unstable regions of protostellar disks fed by infall (Kratter et al. 2010a, 2010b). Although disks tend to have radii of 200 au and less (Mann et al. 2014), interactions in the disks may eject these objects into orbits extending beyond the disks, thereby populating the separations probed in this survey (Basu & Vorobyov 2012). However, neither of these suggestions can explain the high CF measured in Taurus, which shows stellar surface densities and gas column densities similar to the low stellar density regions in Orion. This suggests that other physical factors may come into play. A detailed comparison of the physical properties of the molecular gas in Taurus with those in Orion is needed to identify what factors can lead to a high CF in Taurus and a low CF in the low stellar density regions in Orion. Finally, we have not ruled out random fluctuations as being the cause for the difference. Although the CP varies in a manner similar to the CF, the null hypothesis of a constant CP for all protostars has not been excluded. If the presence of a companion around a protostar is a random process with a probability equal to CP, we cannot rule out that the observed environmental variations in the CP for both protostars and disks are due to random fluctuations at better than the 4%–10% level, depending on the test applied (see Section 4.6). In this interpretation, the enhanced formation of companions in dense regions would be fortuitous and not the result of environmental conditions. Larger samples are needed to reliably distinguish between random fluctuations and deterministic, environmentally driven variations in the CF.

The CF between 100 and 1000 au may also be influenced by complex dynamics. Simulations of nonhierarchical triples predict rapid declines in the CF between 100 and 1000 au during the first 100,000 years of the protostellar phase as one of the companions is ejected or launched into a wider orbit (Reipurth 2000; Reipurth et al. 2010). This has several implications for our interpretation. First, such chaotic processes may effectively be modeled as a random process and therefore provide a mechanism for the random fluctuation interpretation. Alternatively, changes in the initial separation due to environment can affect the subsequent evolution of nonhierarchical ternary or higher order systems, resulting in different CFs. For example, simulations show that the timescale for the dynamical evolution of nonhierarchical systems decreases when the initial mean separation of the system is reduced (Reipurth et al. 2010). Finally, for our sample of more evolved, near-IR protostars and pre-main-sequence stars, the rapid decrease in the CF may have already occurred, as implied by the similar CFs that we observe for the protostars and pre-main-sequence stars. Hence, if this rapid evolution of systems is present in Orion, then the CF during early phases of protostellar evolution may have been higher.

The separations probed by our survey are expected to be strongly effected by dynamical effects in clusters (Marks & Kroupa 2012; Parker & Meyer 2014). How might variations in the primordial CF supported by this paper affect the interpretation of previous work on the evolution of multiplicity in the ONC and other clusters? As mentioned in the previous section, the CF for separations of 100–660 au for the combined high and low density is 0.085, identical to the value found in the ONC by Reipurth et al. (2007). However, we now find that the CF is higher in dense regions. If we take the CF between 100 and 660 au, it is 0.098 for dense regions and 0.071 for low stellar density regions using a threshold NN5 of 30,000 au and the HST sample. Thus, if the initial CF for the star formation in the ONC was similar to what we find in dense regions, there has been a decrease in the CF for the ONC. This supports in part the conclusion of Marks & Kroupa (2012) and Reipurth et al. (2007) that the CF at 100–1000 au decreases due to interactions between stars in the cluster environment. However, we find that the initial multiplicity of the young stars is not invariant, as assumed by Marks & Kroupa (2012). Thus, these results imply a more complicated path to the field CF, where the distribution of environments affects both the initial, primordial CF created by the star-formation process and the subsequent evolution of the CF in clustered environments.

Future and upcoming studies promise many tests of the possible mechanisms for an environmentally varying CF. Ongoing measurement of the masses of the primaries and the companions in both high- and low-density environments can test the formation mechanism. If variations in the CF result from the enhanced formation of companions in protostellar disks located in regions of high gas and star density, then the mass function should be weighted to lower mass objects (Kratter et al. 2010a, 2010b; Basu & Vorobyov 2012). If, instead, the companions are produced by prompt and turbulent fragmentation, then the companion mass function should be similar to the IMF of the primary stars. Better statistics are required, particularly to test whether random fluctuations can explain the observed variations in the CF. An ongoing HST/WFC3 snapshot survey of protostars within 500 pc is now being executed that will increase the sample of protostars observed by HST by more than 30%. Furthermore, approved surveys of Orion Class 0 sources with the VLA and all Orion protostars with ALMA will allow us to examine multiplicity for very young Class 0 sources. Younger companions may also be probed with the JWST, which promises higher angular resolution at longer wavelengths. Expanding the range of separations may provide new tests of environmentally dependent multiplicity and its origin. The VLA, ALMA, and Spitzer can identify wide binaries with separations greater than 1000 au and probe environmental differences at these separation ranges (e.g., Chen et al. 2013). The JWST, VLA, and ALMA, as well as AO observations on 8+ meter telescopes, also promise data on companions at separations within 100 au, providing the ability to search for environmentally dependent multiplicity at shorter separations.