Issue 
A&A
Volume 599, March 2017



Article Number  A14  
Number of page(s)  33  
Section  Galactic structure, stellar clusters and populations  
DOI  https://doi.org/10.1051/00046361/201629398  
Published online  20 February 2017 
Multiplicity and clustering in Taurus starforming region
I. Unexpected ultrawide pairs of highorder multiplicity in Taurus
^{1} Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France
email: isabelle.joncour@univgrenoblealpes.fr; joncouri@gaia.astro.umd.edu
^{2} Department of Astronomy, University of Maryland, College Park, MD 20742, USA
^{3} Astronomy Department, University of California, Berkeley, CA 947203411, USA
Received: 26 July 2016
Accepted: 24 November 2016
Aims. This work analyses the spatial distribution of stars in Taurus with a specific focus on multiple stars and wide pairs in order to derive new constraints on star formation and early dynamical evolution scenarios.
Methods. We collected the multiplicity data of stars in Taurus to build an uptodate stellar/multiplicity catalog. We first present a general study of nearestneighbor statistics on spatial random distribution, comparing its analytical distribution and moments to those obtained from Monte Carlo samplings. We introduce the onepoint correlation Ψ function to complement the pair correlation function and define the spatial regimes departing from randomness in Taurus. We then perform a set of statistical studies to characterize the binary regime that prevails in Taurus.
Results. The Ψ function in Taurus has a scalefree trend with a similar exponent as the correlation function at small scale. It extends almost 3 decades up to ~60 kAU showing a potential extended wide binary regime. This was hidden in the correlation function due to the clustering pattern blending. Distinguishing two stellar populations, single stars versus multiple systems (separation ≤1 kAU), within Class II/III stars observed at high angular resolution, we highlight a major spatial neighborhood difference between the two populations using nearestneighbor statistics. The multiple systems are three times more likely to have a distant companion within 10 kAU when compared to single stars. We show that this is due to the presence of most probable physical ultrawide pairs (UWPs, defined as such from their mutual nearest neighbor property), that are themselves generally composed of multiple systems containing up to five stars altogether. More generally, our work highlights; 1) a new large population of candidate UWPs in Taurus within the range 1–60 kAU in Taurus and 2) the major local structural role they play up to 60 kAU. There are three different types of UWPs; either composed of two tight and comparatively massive stars (MM), by one single and one multiple (SM), or by two distant lowmass singles (SS) stars. These UWPs are biased towards high multiplicity and higherstellarmass components at shorter separations. The multiplicity fraction per ultrawide pair with separation less than 10 kAU may be as high as 83.5 ± 19.6%.
Conclusions. We suggest that these young premain sequence UWPs may be pristine imprints of their spatial configuration at birth resulting from a cascade fragmentation scenario of the natal molecular core. They could be the older counterparts, at least for those separated by less than 10 kAU, to the ≤0.5 Myr prestellar cores/Class 0 multiple objects observed at radio/millimeter wavelengths.
Key words: open clusters and associations: individual: Taurus / stars: variables: T Tauri, Herbig Ae/Be / binaries: general / binaries: visual / methods: data analysis / stars: statistics
© ESO, 2017
1. Introduction
The formation of stars results fundamentally from the runaway unbalance between overwhelming inwards selfgravity and outwards pressuregradient forces (thermal, turbulence, magnetic, radiation, etc.) inside a perturbed gas cloud. Star formation has been known to be a complex process since at least the first model proposed almost three decades ago of a single isolated lowmass star formed out of one dense gas core (Shu et al. 1987). It also involves radiative processes (heating and cooling), generation and decay of turbulence and magnetic fields, chemical evolution of molecules, as well as dust and feedback mechanisms of newly formed stars in their natal environment. A realistic picture must be even more complicated since the stars are rarely born in isolation (Lada & Lada 2003) and, moreover, they are generally not single but part of multiple systems (Mathieu 1994; Duchêne & Kraus 2013). Overall, the whole story of star formation spans over ten orders of magnitude lengthscale range. Giant gas clouds, having typical sizes from 10 pc to 100 pc (Dobbs et al. 2014), fragment, cool, and collapse to form dynamic clusters of young multiple objets whose sizes are approximately a solar radius. Since newly born stellar objects are closely linked to their gaseous progenitors (Hartmann 2002), their spatial distribution and their multiplicity properties offer important tests on star formation models.
Giant molecular clouds appear to be anything but uniform. They are highly substructured into a hierarchy of gas clumps and into highly intertwined filamentary networks, as strikingly revealed by the recent Herschel Gould Belt Survey (André et al. 2014). In turn these filaments shelter single or chains of molecular cores (Tafalla & Hacar 2015). Whether such heterogeneous patterns remain as the relics of the star formation processes in young stellar populations or are totally and quickly erased by subsequent dynamical evolution is an important question.
To seek such imprints of star formation, the giant Taurus mole cular cloud is an ideal nursery. Its proximity (137 pc, Torres et al. 2007), the youth (1–10 Myr) of its stellar population (Bertout et al. 2007; Andrews et al. 2013), and its low stellar density bring a complete census of approximately 350 young lowmass stars down to 0.02 M_{⊙} (Luhman et al. 2010; Rebull et al. 2010). As there is a small spatial dispersion of (proto)stars within the gaseous filament (Hartmann 2002; Covey et al. 2006), these stars are most probably born where they are now observed.
The current paradigm of star formation attributes a seemingly fundamental role to dense molecular cores thought to be the very cradles of stars (WardThompson et al. 2007). The longstanding scenario proposed so far to form one lowmass star from one collapsing single core (Shu et al. 1987) predicts that the dense cores set the final mass of young stars. This is supported by the correspondence of the shapes of the core mass function (CMF) and the initial mass function (IMF) of stars, although the peak of the IMF is shifted by a factor 2–3 towards lower mass. This shift may be due to a partial 30–40% efficiency gas to star conversion (Alves et al. 2007; Lada et al. 2008; Marsh et al. 2016; Könyves et al. 2015), although this estimate is higher than the starformation efficiency (1–10%) estimated for the whole cloud. Alternatively, such a shift between the core mass function and the initial mass function may also result from the multiplicity of stars born within one individual core that has undergone multiple fragmentation (Goodwin et al. 2008). This may lead to a coretostars efficiency as high as 100% (Holman et al.2013; see for a review Offner et al.2014).
The multiplicity appears to be a key feature in starforming regions. The companion frequency of young stars in Taurus is generally twice that of field stars (Duchêne & Kraus 2013), and the formation of coeval and wide binaries up to 5 kAU is a common outcome of the star formation process (Kraus & Hillenbrand 2009a). Furthermore, since multiplicity appears to decrease with age (see for a review Duchêne & Kraus2013) due to dynamical evolution, observing multiplicity at this young age means that the multiple systemcore picture of star formation is the rule rather than the exception. In this picture, the maximal size of multiple systems cannot exceed the size of their progenitor cores, typically approximately 20 kAU, that is, 0.1 pc (WardThompson et al. 2007). If the dense cores are de facto the smallest bricks to build a young stellar population, we then expect that spacing between stars as well as binary semimajor axes be related to parental molecular core size and spacing, provided that dynamical evolution has not erased early imprints. Furthermore, if the formation of multiple systems and single stars share a common scenario, we do not a priori expect any difference in spatial distribution between the two.
Several spatial studies have been performed in Taurus. Some were based on the stellar twopoint correlation function and the related mean surface density of companions to examine the scalefree behavior of clustering modes (Gomez et al. 1993; Larson 1995; Simon 1997; Gladwin et al. 1999; Hartmann 2002; Kraus & Hillenbrand 2008). Other studies based on firstnearestneighbor separation (1NNS) statistics analysis aimed either at comparing the spatial distribution of stars in Taurus to a random distribution (Gomez et al. 1993) or at studying their distribution as a function of their mass and spectral energy distribution (SED) (Luhman 2006; Luhman et al. 2010).
In this work, after a description of the stellar data catalog we have built to perform our studies (Sect. 2), we apply these two statistical tools (1NNS distribution and twopoint correlation function) to analyse the spatial distribution of stars in Taurus with a focus on multiples with respect to single stars, while introducing the onepoint correlation function Ψ (Sect. 3). We then assess the statistical properties of 1NNS couples while defining ultrawide pairs (UWP) as mutual nearest neighbors in the range of 1–100 kAU and we show the crucial role they play in spatial clustering features (Sect. 4). We make a synthesis of results and open a discussion (Sect. 5) before a summary to conclude the paper (Sect. 6).
2. Data
2.1. Input catalog
Although a recent catalog of Taurus members has been released including newly detected midinfrared Widefield Infrared Survey Explorer (WISE) sources (Esplin et al. 2014), we adopted the catalog containing 352 Taurus members that offers a full census of members down to 0.02 M_{⊙} (Luhman et al. 2010; Rebull et al. 2010), which we supplemented with stellar multiplicity data. The newly identified WISE members are mainly due to a wider coverage of the Taurus region compared to Spitzer’s, similar to the new candidates recently identified in Taurus using the Sloan Digital Sky Survey (Luhman et al. 2017). These new members belong primarily to the dispersed stellar population and have typically not been observed at high angular resolution (HAR), hence multiplicity information is not complete for these sources. Therefore, setting this population aside should not significantly affect our conclusions surrounding local pairing statistical properties. Moreover, the Luhman et al. (2010) catalog provides a full and coherent analysis of SEDs as obtained from the Spitzer Infrared Array Camera and Multiband Imaging Photometer, leading to a system of four Classes of SED (Class 0, I, II, and III objects) that we will use in our statistical studies. This classification is widely used to assess the evolutionary state of the material surrounding the star (Lada & Wilking 1984; Lada 1987; Andre et al. 1993; Greene et al. 1994). Further studies have revealed that pure SED analyses can be misleading (Carney et al. 2016) for individual objects. However, although uncertainties remain on the individual classification, there is a clear statistical sequence of ages as a function of SED classification when the population is taken as an ensemble.
Taurus is known to be the archetype of lowmass star formation in loose groups (Jones & Herbig 1979; Gomez et al. 1993; Kirk & Myers 2011). The star mass estimates were taken from the latter work exploiting the same catalog of stars from Luhman et al. (2010). We note that the stellar masses have been estimated from their spectral type while assuming a constant age of 1–2 Myr for all stars of the complex and using an ad hoc isochrone at 1–2 Myr composed of a sequence of several theoretical evolutionary tracks depending on the mass range (see Kirk & Myers2011, for details). Since a recent careful analysis based on optical spectroscopic study of young stars in starforming regions, including Taurus, has shown that spectral types measured over the past decade are mostly consistent with the new evaluation (Herczeg & Hillenbrand 2014), the main source of uncertainties surrounding mass estimation comes from the hypothesis of a unique age for all the members (1–2 Myr).
The mass was estimated for each star using the effective temperature converted from its spectral type and an adhoc isochrone at 1–2 Myr composed of a sequence of several theoretical evolutionary tracks depending on the mass range (see Kirk & Myers2011, for details). For a given spectral type, assuming a younger age for a star than reality will tend to underestimate its mass. But, even if the absolute value of the mass is subject to great uncertainties (30 to 50%), yet we do not expect a systematic differential bias in mass estimation. Besides, there is no indication of a mixture of components in the age distribution that would suggest the presence of distinct episodic generation of stars, nor any evidence of a systemic dynamical evolution that would divide the stars into distinct agebased populations. New generation of stars and dynamical evolution appear to be continuous processes that globally lead to a smooth dispersion of age. For instance, multiple systems and single stars have a same mean age and display the same age dispersion (White & Ghez 2001). Thus, assuming a single age for all the stars will bias the absolute value of each derived mass, but the relative difference between them, which is used in our analysis, should be much less affected.
Fig. 1 Spatial distribution of stars within the Taurus region superimposed on nearinfrared extinction (unit mag, gray scale) from Juvela & Montillaud (2016). Dashed polygon: restricted spatial window W_{in} in Taurus complex; blue filled dots: random sampling. Colorcoding for Taurus stars: red plus marks: clustered stars [1NNS ≤ 0.1°]; black plus marks: “inhibited” regime stars [0.1°< 1NNS ≤ 0.55°]; green plus marks: isolated stars [1NNS beyond 0.55°] (see Sect. 3). 
2.2. Multiplicity
To construct a census of stellar systems in Taurus, we began by compiling a list of all known multiple members of the region from almost three decades of increasingly higherangularresolution observations (see references in Table C.1). We also assigned a flag to each star when it has been observed at HAR.
In order to separately study the populations of single stars and multiple systems and to distinguish clustering from multiplicity, we grouped together all stars that are separated by less than 7′′ (0.92 kAU ~ 1 kAU) to form a single entity. Throughout this paper we refer to them simply as “multiple systems”. This threshold is, to a certain extent, arbitrary. However, two reasons motivate this choice. First of all, it is the lower threshold that defines wide binaries (separation >1 kAU) and secondly (the main reason), this threshold is close to the 6″ beam of the SpitzerMIPS data. Thus, we ensure that the census of neighbor stars beyond 1 kAU is completely independent of the Class of the Taurus members.
We assigned the total number of stars within a 1 kAU radius n_{∗} (for a single star n_{∗} = 1, for a binary n_{∗} = 2, ...) to each entry of the catalog, and we note that this parameter is only reliable for Class II and III type stars that have been observed at HAR, using techniques such as speckle interferometry imaging and adaptive optics with spatial resolution between 0.05″ to 2″. These techniques probe the immediate neighborhood of a star, at least beyond ~10 AU at the distance of Taurus (i.e., 5−20 AU depending on the technique). Very close binaries (i.e., ≲10 AU) most probably originate from a distinct scenario from other binaries, as they are probably not formed by fragmentation in situ (Bate et al. 2003). We thus consider spectroscopic binaries and the closest visual binaries (below ~20 AU) as single stars to ensure the consistency and completeness of our visual binary analysis, since no exhaustive spectroscopic binaries survey has been performed till now.
The resulting catalog (see Table C.1) contains 338 sources. Of these, 250 have been observed at HAR and amongst the latter 94 are multiple systems, resulting in a mean multiplicity fraction (MF) of 40%, which could truly be as high as 54%, if all the remaining stars that have not been observed at HAR are multiple systems. A similar multiplicity fraction is obtained when considering only Class II/III stars (i.e., MF = 40% ± 10%). The mean companion frequency for Class II/III stars observed at HAR is 48%. We note that a chisquared test shows that the proportion of Class I, II, and III in the star sample observed at HAR is statistically the same when considering either single stars or multiple stars. The disk fraction (the number of Class II stars over the number of Class II and III stars) is also statistically the same in the two populations. Using Class as a proxy for the age, we obtain, as already stated by White & Ghez (2001) that uses theoretical stellar evolution models, that the two populations share the same age and the same range of age dispersion.
2.3. Samples
The sources in Taurus, as was first pointed out by Gomez et al. (1993), are not randomly distributed on the sky (see Fig. 1). Some sets of stars appear grouped at a global scale in two structures located North (in L1507 molecular cloud) and SouthWest (in L1543 molecular cloud) with respect to the three main filaments, while some of them are tightly subgrouped at small scales within the main central filaments. Yet, other stars appear more spatially dispersed between the main structures. One of the goals of this work is to more precisely quantify their spatial distribution based on nearestneighbor analysis, now that the census of star population in Taurus has more than doubled since the seminal work of Gomez et al. (1993).
From our catalog we extracted three distinct samples (defined in Table 1) that we use to test different hypotheses throughout this work, notably for multiplicity testing.
Subsamples of stars in Taurus.
3. Local spatial analysis
In this section, we perform the spatial analysis of stars in the entire Taurus region and we assess the difference in the spatial distribution between multiple systems and single stars based on their respective first nearest neighbor separation (1NNS) statistics. Nearestneighbor statistics have been used for decades in Taurus to; (1) demonstrate departures from random uniform distributions (Gomez et al. 1993); (2) show that the spatial distribution of brown dwarfs and stars are undistinguishable (Luhman 2004); and (3) study the spatial distribution of stars as a function of their Class (Luhman et al. 2010). We start with some preliminary work on 1NNS statistics to advocate the use of the theoretical 1NNS distribution derived for a random distribution in an infinite medium as a reliable proxy for a random distribution enclosed in a finite and irregularly shaped window provided that the window is large and wellpopulated. We then define the onepoint correlation function Ψ, which complements the twopoint correlation function, to study the binary regime. This function allows us to quantify the departure of a population from a random uniform one and to identify distinct spatial regimes (clustering, inhibition, and dispersion).
Most of our work was performed using the free R software environment for statistical computing and graphics (R Core Team 2015) using a set of packages such as Spatstat, Hmisc, fields, FITSio, calibrate, xtable, astrolab, deming, mixtools, and magicaxis (Baddeley & Turner 2005; Harrell et al. 2015; Nychka et al. 2015; Harris 2013; Graffelman 2013; Dahl 2014; Chakraborty et al. 2014; Therneau 2014; Benaglia et al. 2009; Robotham 2015).
3.1. 1NNS study in Taurus
Moments of the 1NNS distribution in Taurus.
This subsection aims at analyzing the spatial distribution in Taurus based on a 1NNS analysis. The 1NNS of a star is, by definition, the distance to its nearest neighbor star. Throughout this paper, we deal with the projected 1NNS, computed for each star as the minimum value of projected angular distance between each pair of stars onto the celestial sphere (see Appendix A).
After estimating the mean surface density of stars, we first compare the Taurus 1NNS distribution in the three main filaments to that associated to the entire region. We then compare the theoretical random 1NNS distribution obtained for an infinite random spatial process to the 1NNS distribution derived from Monte Carlo samplings in a finite irregularly shaped window, and finally we compare Taurus 1NNS distribution to the 1NNS random distribution.
3.1.1. Mean stellar surface density in Taurus
Fig. 2 Empirical cumulative distribution of the 1NNS distribution for the whole region (black) and for the central part within W_{in} region (blue). The blue vertical line represents the clustering length, r_{c} = 0.1°. 
To derive an estimate of the mean surface density ρ_{w} of the spatial process in Taurus, we focus on its central part as the least heterogeneous area and define the spatial polygonal window W_{in} that encloses the central three main filaments (see Fig. 1), to get: (1)where A_{w} = 48.5deg^{2} (i.e. 289.57 pc^{2}) is the projected area of the window W_{in} that encloses N_{w} = 252 stars. This window defines our star sample S_{2} (see Table 1).
In the following paragraph, we estimate the uncertainty of the mean projected stellar surface density of Taurus. A conservative estimate of the uncertainty on the surface density may be obtained when choosing for the window, the convex hull of the stars enclosed within the polygonal window W_{in}. The convex hull of a set of points in the Euclidean space is defined as the smallest convex set of points that contains the whole set of points.The area of the convex hull, A ~ 33deg^{2}, being less than the polygonal window, the surface density estimate mechanically increases by a factor of 1.6 (N_{w}/A ~ 8 deg^{2}). Conversely, taking the convex hull area of the whole Taurus region (~202 deg^{2}) containing 338 objects lowers the mean surface density to ~2 deg^{2}. Combined, these extreme cases yield a probable range for the projected stellar density of Taurus of (2)
3.1.2. 1NNS distribution across Taurus
The 1NNS median value evaluated within the full region of Taurus (see Table 2) does not differ from the 1NNS median value evaluated within the three main filaments, and their 1NNS cumulative distributions do not differ either for length scale below r_{c} ~ 0.1° (see Fig. 2). This indicates that the spatial properties of the spatial distribution of stars are similar across the whole region of Taurus below this angular scale. Above the 0.1° scale, the two cumulative distributions diverge due to the presence of highly dispersed stars between the three main filaments and other main stellar concentrations in Taurus. As a result of these, the 1NNS distribution of the whole region has a more populated largescale tail. Nonetheless, up to the third quantile, the 1NNS quartiles are similar. As a matter of fact, due to its highly asymmetrical shape, the Taurus 1NNS distribution is better evaluated using quantiles rather than the mean and standard deviation that are more conveniently suited for symmetrical functions.
3.1.3. 1NNS distribution: theoretical random vs. Monte Carlo samplings
In this subsection, we compare the theoretical 1NNS distribution for a spatially random process in an infinite media to the 1NNS distribution obtained from Monte Carlo samplings in a finite, irregularlyshaped window to establish that the former can be used as a reliable proxy for the latter.
The random theoretical 1NNS distribution w_{R}(log r) derived for an infinite media (see Appendix B for the full development) is given by: (3)where ρ is the mean intensity of the random process (the average surface density of stars) and r is the 1NNS. We want to compare this theoretical distribution to the 1NNS histogram obtained from 10 000 random Monte Carlo samplings within the finite and irregularlyshaped window W_{in} (see Fig. 3). Their respective moments are similar (see Table 2), with theoretical moments computed from the following expressions: where , r_{0.5}, σ^{2}, σ^{2}, and γ are the mean, median, variance, standard deviation, and skewness factor (see Appendix B for their derivation). The moments from random Monte Carlo samplings were obtained from: (5)where Λ is the ordered sequence of the 1NNS from the smallest to the highest value, the standard errors of the moments and the standard error of the median are taken respectively from Ahn & Fessler (2003) and Kendall & Stuart (1977).
Fig. 3 Probability distribution function (PDF) of the 1NNS distributions observed in Taurus (solid blue histogram: within the limited window W_{in}, that is, the three main filaments; dashed blue histogram: the entire Taurus region) and predicted for a random distribution (black histogram: for a Monte Carlo sampling within the window W_{in}; red curve: unlimited theoretical random 1NNS distribution, see Eq. (3)). Both random populations are drawn using the same mean surface density ρ_{w} = 5deg^{2}. Increasing the mean density by a factor of two leads to a shift towards the left whose amplitude is represented by the green arrow. The Monte Carlo sampling in a finite window and the unlimited theoretical random 1NNS do not significantly differ. The 1NNS distributions computed either in the whole Taurus or within W_{in} do not significantly differ either. However, the Taurus and random 1NNS distributions are statistically highly inconsistent, with a pvalue of ~10^{16} for the KS test. The spatial distribution in Taurus is clearly far from being random. 
We thus conclude that finite size and edge window effects do not significantly alter the distribution of 1NNS from that obtained in an infinite medium at the level of our observational uncertainties. To prove this statement quantitatively, we use a bootstrapping technique, comparing 10 000 realizations of a) a sample of N ~ N_{w} = 252 1NNS obtained from standard rejection sampling technique using analytical random 1NNS distribution (Eq. (3)) to b) the 1NNS associated to random Monte Carlo samplings of 252 stars enclosed within W_{in}. The twosample KolmogorovSmirnov (KS) and AndersonDarling (AD) statistical tests indicate that the two 1NNS distributions cannot be statistically distinguished with a mean pvalue and standard deviation respectively of 0.36 ± 0.28 for the KS test and 0.34 ± 0.27 for the AD test.
We have therefore shown that the random theoretical 1NNS probability density function (Eq. (3)) may be used as a reliable proxy to describe a 1NNS random distribution even in the case of a random population enclosed in a finite and irregularly shaped window, such as W_{in}. We mention, as an asset, that using random theoretical 1NNS distribution proxy avoids large Monte Carlo sampling computations required to populate highly improbable bins of very small spacings associated with tightly clustered regions.
3.1.4. Spatial distribution in Taurus
The comparison between the Taurus and random 1NNS distributions clearly reveals an excess of short spacings in Taurus (see Fig. 3), for which the 1NNS distribution has smaller central values (median and mean) and a wider dispersion (see Table 2). The 2sample KS and AD tests, give a pvalue less than ≤2.2 × 10^{16} and ≤10^{40}, respectively, that the two samples are mutually consistent. Thus, we can firmly reject the hypothesis that the spatial distribution of stars within Taurus is random.
This result is not surprising, since the same conclusion was already reached by Gomez et al. (1993) based on a sample of 139 stars. But the expansion of the stellar census by nearly 200 new stars in Taurus strengthens this conclusion. The median value of the 1NNS distribution we derived is r_{0.5} ~ 0.067° (~0.16 pc, i.e. 33 kAU); almost half the value obtained by these authors. This is what we expect when the spatial distribution of additional stars follows a similar spatial pattern from the original population, since from Eq. (4d), the median value is reduced by a factor of ~.
This also underlines the modest usefulness of the absolute value of the median nearestneighbor distance to determine clustering property when the stellar census is incomplete. In the following section, we introduce an alternative criterion based on the onepoint correlation function to quantitatively assess local departures from pure randomness and to determine characteristic clustering scales.
3.2. One and twopoint correlation statistics
In order to identify a criterion that quantifies the departure of a spatial distribution of stars from randomness, we introduce the onepoint correlation function based on the 1NNS distribution. This function is a new tool to assess and quantify the binary and spatial clustering regimes of stars and aims at complementing the twopoint correlation statistics.
3.2.1. The twopoint and pair correlation functions
The twopoint correlation function (TPCF) ξ_{V}(r) is defined for a 3D spatial process as a measure of the probability to get any excess of stellar pairs with a separation distance vector r above random expectation: , where ρ_{V} is the mean volumic density of the uniform random process and dP is the infinitesimal joint probability to find a pair of stars respectively centered in the volume elements dV_{1} and dV_{2} and separated by r (Peebles 1980). Similarly, for an isotropic and stationary spatial process, the standard projected angular twopoint correlation function ξ(r) is then defined as the probability to find a pair of stars each separated by a projected angular distance r above random expectation: dP = ρ^{2}(1 + ξ(r))dΩ_{1}dΩ_{2}, where ρ is the density of stars per steradian and dP is the infinitesimal joint probability to find a pair of stars each within a unit solid angle dΩ_{1}, dΩ_{2} and separated by a projected angular distance r.
Besides the twopoint correlation function, the second order statistics may also be evaluated from related entities:

the main surface density of companions (MSDC) above randomexpectation: Σ(r) = ρξ(r);

the pair correlation function g(r) as a measure of the probability to have any pair of stars separated by r: g(r) = 1 + ξ(r).
When the spatial distribution is random, all these functions are constant (ξ(r) = 0,g(r) = 1,Σ(r) = 0). Although different estimators may be used to compute the pair correlation function (Landy & Szalay 1993; Kerscher et al. 2000), we use the simple and natural one defined as: (6)where and are the histogram of separations between stars within Taurus and in a randomly distributed catalogue, respectively.
Fig. 4 Onepoint correlation function Ψ (Eq. (7), blue symbols) and estimated pair correlation function g (gray symbols) using the sample of stars within W_{in}. The pair correlation function has been estimated using the estimator given in Eq. (6) while performing 10 000 random Monte Carlo samplings within W_{in}. The 1σ vertical errorbars of Ψ and g are estimated from Poisson statistics. The solid horizontal blue line represents both the constant Ψ = 1 and ĝ = 1 functions, that is, the onepoint and pair correlation function expected for a spatial random distribution. Vertical black solid lines delimit three spatial regimes, clustering, inhibition, and dispersion, derived from the crossing of the Taurus Ψ function with the Ψ = 1 horizontal line associated to random spatial distribution. The vertical dashed purple line indicates the 7′′ spacing threshold used to group stars into a unique multiple system. The dotteddashed green line is the Ψ ∝ r^{1.5} function resulting from a linear regression fit taking into account vertical errorbars over ten 1NNS logarithm bins (Δlog r = 0.2). The Ψ function may reveal an extended binary regime from 1.6 kAU to 100 kAU in Taurus, which has remained unidentified beyond approximately 5 kAU in previous works that used the pair correlation function alone, since the latter reveals a power law break around that separation. 
Since the pair correlation function estimate is mainly sensitive to the shape and extension of the window, no straightforward theoretical expression can be derived for on general grounds, even in the random spatial distribution case. We then compute the pair correlation function using Eq. (6) while performing 100 000 random Monte Carlo samplings within W_{in} (see Fig. 4).
The slope and the shape of the resulting Taurus pair correlation fonction are comparable to those obtained in previous works studying the twopoint correlation function or the main surface density of companions (Gomez et al. 1993; Larson 1995; Simon 1997; Gladwin et al. 1999; Hartmann 2002; Kraus & Hillenbrand 2008). Together with our own, these studies show that the second order statistics exhibit at least two distinct powerlaws: a steep binary/multiple regime at small scale and a smoother clustering regime beyond a breaking point. Larson (1995) suggested that; (1) the observed break at 0.04 pc (8.25 kAU) can be the signature of the transition from a binary to a clustering regime; (2) it could correspond to the Jeans length for the gas in cool dense molecular cores; (3) the scalefree binary regime may reflect the scalefree fragmentation of collapsing clumps; and (4) the selfsimilar clustering regime may reflect the hierarchical, fractal, and spatial distribution of the gas from which the stars formed.
Later on, Kraus & Hillenbrand (2008) found an extended transition zone between the binary regime and the global clustering of stars in groups rather than a sharp break between the two regimes. This transition zone was found to begin at 0.02 pc (4.25 kAU) and end at 0.2 pc (42.5 kAU). They proposed that part of this transition zone, which is almost flat from 0.11 pc (23.52 kAU) to 0.2 pc (42.5 kAU), is due to the random motion of stars that may smooth out primordial stellar association and lead to a quasiconstant surface density of pairs.
However, some caution has to be taken when interpreting the twopoint correlation function beyond the binary regime, especially when the spatial distribution of stars deviates from isotropy. Indeed the standard form of the twopoint correlation function is obtained with the hypothesis of a homogeneous and isotropic distribution. However, the Taurus complex is not isotropic, as evidenced by its elongated gaseous structures in which stars are preferentially located. The break in the twopoint correlation function between the binary and the clustering regimes may thus be partly due to the filamentary geometrical anisotropy, rather than being only a signature of a hierarchical/fractal gas structure as proposed by Larson (1995).
Indeed, the observed break in the twopoint correlation function could also be linked to the observed width of filaments, which was found to be relatively universal at approximately 0.1 pc (André et al. 2014) even though Ysard et al. (2013) found that filament width may vary by up to a factor of approximately four. This spread of filament widths may also contribute to the transition elbow zone that is present between the binary and clustering regimes. Thus, the elbow in the correlation function that extends up to 0.25 pc may be related to geometrical spatial stellar structures (see Joncour et al., in prep.).
3.2.2. The onepoint correlation function
Similarly to the pair correlation function, we define the onepoint correlation function, Ψ, as the probability of having a closest star located at r from any chosen star at random, which in turn defines an equivalent local stellar density function σ: dP = σ(r)dV, with σ(r) = ρΨ(r). The Ψ function, which describes the spatial location of stars, gives first order variation trends of the spatial process relative to a random distribution. Instead, the pair correlation function is associated to secondorder characteristics describing the spatial colocation of stars. One estimator of the Ψ function may be defined from the ratio of the 1NNS Taurus distribution w_{Tau}(log r) to the 1NNS random distribution w_{R}(log r) obtained either by Monte Carlo simulations or from theoretical random probability density function (Eq. (3)): (7)We note that this function is usually called the nearest neighbor ratio when it is evaluated as the average over the whole range of 1NNS values. Here we evaluate it bin per bin of spacing. For a random spatial distribution, the onepoint correlation function reduces to unity (Ψ(r) = 1).
The onepoint correlation Ψ function can be fitted from 1.6 kAU to ~100 kAU, by a power law using Pearson’s chisquared linear regression taking into account vertical errorbars (see Fig. 4): (8)with at the 95% confidence level. This allows us to derive a fractal dimension D_{f} associated to the binary regime in Taurus. Since the random distribution in 2D increases as w_{R} ∝ r^{2} for (see Appendix B), the 1NNS distribution within Taurus w_{Tau} follows a shallow, rising powerlaw function with the following exponent: (9)We note, that this scalefree behavior of the Ψ function with an exponent lower than the projected 2D euclidian space dimension is equivalent to a uniform distribution in a 2D projected fractal space with the fractal dimension D_{f} = 0.5 (Sakhr & Nieminen 2006).
3.2.3. Departure from randomness
The one point correlation function Ψ may be used to assess departure from a spatially random distribution. From its definition, regimes where Ψ is below unity indicate a deficit of stars (inhibition) with respect to random, whereas where it is above unity reveals an excess of stars, either a clustering trend at small spacings (aggregation) or a dispersion trend at large spacings. The Ψ function evaluated for Taurus thus reveals three spatial regimes (see Fig. 4). The stellar clustering regime covers the range of spacings below ~0.1°, where Ψ(r) > 1. This defines the upper clustering length threshold r_{c}: (10)The clustering regime is followed by an “inhibition zone” where Ψ(r) < 1 associated with a typical length scale between ~0.1° (0.2 pc) and ~0.3° (0.7 pc). Finally, Taurus faces an excess of highly dispersed stars for typical length scale above ~0.3°. The “clustered stars” are tightly packed together in specific locations while the “inhibited stars” are more widely spread between the zones of clustered stars (see Fig. 1).
3.2.4. Ultrawide binary regime in Taurus?
The comparison between the pair correlation and the Ψ functions yields new insight into an extended binary regime. Both of them have the same scalefree trend at least in the four logarithm bins of separation (see Fig. 4), from 1.5 to 5 kAU. While the correlation function reaches an elbow from that point on, breaking from the binary regime to reach the clustering regime, the Ψ function extends the same powerlaw trend up to ~100 kAU. Thus, although it remains hidden in the correlation function because of mixing between binarity and clustering, the binary regime in Taurus appears to extend to much larger scales than previously realized.
3.3. Local neighborhood of multiples versus single stars
Fig. 5 First nearest neighbor separation fraction distribution of multiple systems (solid blue histogram) versus single stars (dashed black histogram) for Class II and III objects observed at high angular resolution (HAR) in Taurus region. Each bin represents the number density of stars per unit logarithmic interval in projected separation, that is, the number density fraction of 1NNS per (Δlog r = 0.2) interval. There is a marked excess of 1NNS in the range 1−10 kAU for the multiple systems with respect to single stars. 
In this subsection, we study the spatial distribution of multiple systems and single stars based on the comparison of their 1NNS statistics.
The S_{1} sample in the whole Taurus contains 96 stars flagged as visual multiple systems (in the range of 10 AU to 1 kAU) and 242 a priori single stars. To limit biases, we define the S_{3} sample composed exclusively of Class II and III stars that have been observed at HAR. From that subsample, we separately construct the 1NNS histograms of the 89 multiple systems and the 140 single stars (see Fig. 5) that appear very different one from each other. Both the KS and AD statistical tests indicate that we can reject the theory that they come from the same parent population with a confidence level of 99.8%. The 1NNS fraction of multiple systems having a projected separation less than 10 kAU is noticeably higher than for the single stars in the same range (see Fig. 5).
The excess fraction of 1NNS at close spacings is remarkably high: nearly 40% of the 1NNS of all multiple systems are located within the first five bins, that is, between 1 and 10 kAU. This is more than 2.5 times the fraction found in the same range of separations for the single stars (15%). It is also twice as high as the frequency of visual companions per decade of projected separation found in the 10−2000 AU separation range (see Duchêne & Kraus2013, Fig. 5). Moreover, 63% of the companions of multiple systems in the five first bins are themselves multiple systems (25% of the total multiple systems), against 43% for the singles (6% of total single stars sample).
The local neighborhood of multiple stars is thus more populated than that of simple stars: a multiple system has almost three times more chance of having a companion within 10 kAU, and that companion in turn has a higher chance of being a close multiple system itself than a single star.
3.4. Beyond local neighborhood
To test whether this striking difference between the neighborhood of multiple systems and single stars extends beyond the first neighbor, we analyse the second nearest neighbor distributions, again using the biasfree S_{3} sample. Although there is still a slight excess of companions for multiple systems within the first bin up to 10 kAU, this difference does not appear statistically significant: the KS and AD tests both give a pvalue of 0.15. We further compare their respective kNNS distributions, up to k = 8. They too cannot be statistically distinguished (pvalue of 1 for both KS and AD tests). This confirms the statistical identity between the neighborhood of the two populations, multiples and singles, beyond their first neighbor.
Since the 2NNS distribution is not statistically different for single stars and multiple systems, we can use its median value (0.13°, 0.3 pc or ~60 kAU) as a typical length scale at which these two populations cannot be statistically distinguished. This value is close to the upper clustering threshold r_{c} = 0.24 pc derived in the second section.
Interestingly, we also note that the distribution of the second nearest neighbor of multiple systems cannot be statistically distinguished from the first nearest neighbor of single stars either (pvalue of 0.1224 for the KS test). This suggests that, although the neighborhood of multiple systems has been found to be more populated within 10 kAU, this trend fades away beyond 60 kAU. Based on the knearest neighbor analysis (k ≥ 2), we therefore define a scale of 60 kAU beyond which the stellar environment is statistically identical for single stars and multiple systems.
4. From multiples to ultrawide pairs
Fig. 6 Distribution of ultrawide pair fraction as a function of separation in the three main Taurus filaments (i.e., for stars enclosed within W_{in}, solid blue histogram) versus pair fraction of random mutual pairs (obtained from 10 000 Monte Carlo samplings, thick black histogram) and 1NNS fraction of nonmutual pairs (dashed blue histogram). The solid blue line represent a lognormal fit to the distribution of nonmutual pairs in Taurus. The distribution of mutual pairs of Taurus does not result from a random process (pvalue of 10^{16}), and is statistically different from nonmutual pairs distribution in Taurus (pvalue of 10^{14}). Conversely, the 1NNS distribution of nonmutual pairs in Taurus is reasonably well fitted (pvalue of 0.2) by a longnormal function, with a mean value μ = 4.74 ± 0.04 and a standard deviation σ = 0.46 ± 0.03. 
In this section, we further analyze the reasons for such a divergence between neighborhoods of multiple systems and single stars.
4.1. Mutual nearest neighbors as ultrawide pairs
From the list of all nearest neighbor pairings (see Table C.2), we selected the subset of all mutual nearest neighbors, that is, pairs in which the nearest neighbors are reciprocal. This property serves to probe the most “connected” pairs. It has been used, for instance, as part of the process to identify physical binary candidates in simulations (Parker et al. 2009; Kouwenhoven et al. 2010).
We find that over the 338 first nearest neighbor pairings of stars in Taurus, 45% of pairings arenonmutual pairs (within this type of couples, only one star is the nearest neighbor of the other, while the second one has a different nearest neighbor) whereas 55% of the whole pairings are mutual pairs.
The 1NNS distribution for mutual pairs is markedly different from that of nonmutual pairs; the KS test returns a pvalue of 10^{14} when assuming, as a null hypothesis, that these two populations are drawn from a same parent population (see Fig. 6). The distribution for nonmutual pairs is reasonably well fitted by a lognormal function (pvalue = 0.2, showing that two distributions are statistically the same) with a mean μ = 4.74 ± 0.04 (i.e. 55 kAU) and a standard deviation σ = 0.46 ± 0.03 (~20 to 160 kAU range). This suggests that, unlike the mutual pairs, the distribution of 1NNS for nonmutual pairs in Taurus may be seen as a statistical realization of a multiplicative process, which is converted to an additive process in the log domain, for which the central limit theorem may be applied. Although it is beyond the scope of this paper to establish the physical origin of this multiplicative process, we can speculate that relative random motion and gravitational encounters between physically unrelated stars may contribute to randomization of the spatial location of individual stars. This is not the case for gravitational binaries or common proper motion pairs, however, as they are kept together at the same mutual nearest distance over time.
The mutual pairs in Taurus are very unlikely due to random mutual pairing (see Fig. 6). A KS test comparing the distribution of 1NNS separation of mutual pairs with the W_{in} window to that of the mutual first nearest pairs from an ensemble of random Monte Carlo Poisson samplings within the same window returns a pvalue lower than 10^{16}.
Moreover, it is worth noting that the mutual pairs we identified in the separation range 1–5 kAU, are already known to be true physical binaries, that is, gravitationally bound binaries (Kraus & Hillenbrand 2009b). This further validates our selection of mutual nearest neighbors to identify physically related pairs in that range. We thus expect that a large majority of these mutual pairs to be at least common propermotion pairs, and plausibly even gravitationally bound. From now on, we refer to these systems as ultrawide pairs (UWPs).
The presentday distribution of separations for these UWPs should reflect its counterpart at birth, provided that dynamical encounters are rare enough to keep it unchanged (see discussion Sect. 5). These mutual pairs covering separations from 1 to ~60 kAU are the only type of pairs (versus nonmutual pairs) found in the shortseparation end of the 1NNS distribution. Approximately 60% (resp. 90%) of 1NNS have separations below 20 kAU (resp. 60 kAU) with a relatively longuniform distribution. Thus, the continuous powerlaw trend at small and intermediate scales observed in the Ψ function is mainly due to the existence of these mutual pairs, which define the binary regime. They generate a scalefree clustering pattern and constitute a major structural spatial imprint in Taurus.
4.2. An Öpik law for the UWPs separation
Fig. 7 Ultrawide pair fraction F within the whole Taurus region (sample S_{1}) as a function of separation, plotted using Δδ_{P} = 0.25 dex bins. The solid blue line is a power law fit to the data (, with log α = −1.21 and β = 0.14) between 3 and ~60 kAU (vertical black dashed line). The gray area represents the 95% confidence band. In terms of the power law parameters, the corresponding parameter 95% confidence intervals are log α = [ −1.78,−0.65 ], β = [ 0,0.28 ]. 
In this subsection we discuss the distribution of UWP separation based on the data that we derived (see Table C.2).
The UWPs fraction F in Taurus follows a slightly increasing power law, (see Fig. 7), at least up to an upper cutoff of ~60 kAU (i.e., 0.27 pc). The robustness of this simple fit is confirmed by performing a power law fit to the empirical cumulative function that explains 97.73% of observed variations. Beyond this cutoff, there is a sharp decrease in the two last bins of separations. The first two bins overlap the study conducted by Kraus & Hillenbrand (2009b), who found that the separation of binaries below 4.2 kAU is loguniform in this range (F(log δ) = const. or, equivalently ), i.e., following Opik’s law.
4.3. Class pairing in UWPs
We first focus on the 37 mutual pairs with separation less than 10 kAU (from sample S_{1}), a regime where the distinction between multiple systems and single stars strongly suggests that a physical process is at play. Within this subpopulation, Class I and Class II are overrepresented (and, accordingly, Class III underrepresented) relative to the overall Taurus population (see Table 3). A chisquared test confirms that this difference is statistically significant (pvalue of 0.01). In other words, the stellar population within UWPs separated by less than 10 kAU is less evolved, i.e., younger, than the rest of the population in the Taurus molecular cloud.
Class fraction in Taurus.
Fig. 8 Distributions of SED Class pairings expected from a random pairing process based on the global proportions of Class I, II, and III stars in Taurus (left) and observed among Taurus UWPs (right). The SED Class pairing within UWPs deviates significantly from random pairing, suggesting coevality within these systems. 
In the following, we show that the class pairing found in these UWPs deviates from random pairing. The probability P_{ij} of getting an unordered pair (i,j) without replacement is given by: (11)where N is the the total number of stars in the sample under consideration, { i,j } = I,II,III is the SED Class of stars, P_{i} is the fraction of Class i type star, N_{i} is the number of stars of Class i, and γ = 1 and for i = j whereas γ = 2 and for (i ≠ j).
We compute the random Class pairing probability for each type of pairings from Eq. (11) and the number of stars of each Class N_{i} derived from the first line in Table 3. These probabilities are quite different from those observed among UWPs in Taurus (see Table 4 and Fig. 8). The main difference can be summarized as follows: there are more pairs of the same Class, that is, (I, I), (II, II) and (III, III) pairs, than expected from random pairing. Correspondingly, there are fewer pairs from mixed Classes, (I, II) and (II, III) pairs, and there are no (I, III) pairs at all despite the expectation that there should be ~10% of such pairs based on random pairing. These differences are found to be statistically significant at the 99.74% level based on a chisquared test.
Applied to the larger sample of UWPs covering the whole range of separation (sample P3), we obtain the same result as above with an even higher significance level (p = 10^{4}). These results indicate that the UWPs are most likely coeval.
Class pairings fraction in UWPs compared to random Class pairing.
4.4. High multiplicity fraction within UWPs
Considering that the stars are part of UWPs, there are almost equal numbers of single stars (50) and multiple systems (44). The resulting multiplicity fraction of 47% is ~15% higher than that of the sample of Class II and III stars observed at HAR that are not in UWPs, a result that is statistically significant (pvalue of 0.007 for a chisquared test). In other words, multiple systems are more often within UWPs than not. More interestingly, the multiplicity fraction within UWPs depends on the separation range: it declines from a maximum of ~70% at short separations (at approximately 5 kAU) to the average multiplicity fraction of 47% at large separation of approximately 6 × 10 kAU, above which the multiplicity fraction remains constant (see Fig. 9). Thus, with a confidence level higher than 95%, the multiplicity fraction within UWPs is higher when their separation is shorter.
Fig. 9 Multiplicity fraction of Class II and III stars observed at HAR and members of UWPs as a function of their separation. For a given value of δ, the black (resp. blue) symbol represents the multiplicity fraction computed over all separations smaller (resp. greater) than δ. Errorbars are standard errors computed from the standard formula: where n_{i} is the sample size, and N is the total population size of Class II and III stars observed at HAR within UWPs independently of their separation. The dashed black line represents the mean multiplicity fraction for stars within UWPs whereas the gray line marks the corresponding mean multiplicity fraction for stars that are not in UWPs. The red dotted lines define the interval within which the multiplicity fraction within UWPs is significantly (pvalue of 0.05 or less based on a chisquared test) higher than the mean value. 
4.5. Multiplicity pairing within UWPs
In this subsection, we look at the multiplicity nature (single or multiple) of individual components composing the UWPs. There are 47 UWPs in sample S_{3} with 12 of them being composed of two multiple systems (MM pairs), 15 of them of two single stars (SS pairs) and the remaining 20 pairs of one multiple system and one single star (SM pairs). The probabilities associated to random multiplicity pairings are computed from Eq. (11) and give 0.28, 0.22, and 0.50 to get a SS, SM or MM pair, respectively. Using a Pearson’s χ^{2} test, we cannot rule out a random multiplicity pairing (pvalue 0.59) of the observed UWPs over the whole range of separation. But, the observed probability is not uniform along with the separation. For instance, considering only UWPs whose separation is less than 10 kAU, more than double the expected MM pairs are observed (48% of the UWPs are of MM pair type, 19% of SS and 33% of SM). The same χ^{2} test indicates that a random multiplicity pairing can be rejected with a high significance level (pvalue 0.018). Thus the higher multiplicity fraction for tighter UWPs, as shown in the previous subsection, stems primarily from the fact that multiple systems tend to be paired with other multiple systems in that range of separation (≤10 kAU).
A comparison of the pair fraction distribution between MM, SM, and SS UWPs as a function of their separation reveals indeed that nearly 90% of MM pairs are concentrated below 10 kAU, the SM pairs are distributed essentially loguniformly over the whole range, and SS pairs are preferentially separated far apart (median separation at 35 kAU, see Fig. 10). As a consequence, MM systems are the main contributor to the observed “bump” in the 1NNS distribution of multiple systems outlined in the previous section. Of the 30 multiple stellar systems having their 1NNS less than 10 kAU, 20 of them are within multiple/multiple UWP (65%), 7 members (25%) are within single/multiple UWPs and 3 (10%) are nonmutual pairs. Therefore, the spatial neighborhood difference between multiple systems versus single stars comes down to a higher fraction of “tighter” MM UWPs with respect to the SS UWPs.
Fig. 10 SS (dashed blue), SM (dashed red) and MM (solid blue) pair fraction as a function of their separation. MM pairs are tighter (most are found below 10 kAU) whereas SS pairs are far apart (with a peak at approximately 50 kAU); on the other hand, SM pairs are relatively uniformly distributed. 
4.6. UWP multiplicity
In this subsection, we study the properties of UWPs, as a function of their multiplicity n_{∗}, defined as the sum of the total number of stars within each pair (2 ≤ n_{∗} ≤ 5).
For this analysis, we continue focusing on sample S_{3} (half the population of UWPs), which has the most reliable multiplicity data, since each of the components has been observed at HAR. Within it, we investigate the distribution of the UWPs as a function of their degree of multiplicity. The number of UWPs seems to be rather uniform up to a multiplicity of four (see Table 5), and exhibits a sharp decrease for a multiplicity of five. We note that this behavior differs from the geometric progression for premain sequence closer binaries in Taurus N_{n∗} ∝ b^{− n∗} = 2.6^{− n∗} (Duchêne & Kraus 2013).
Number and type of high order multiplicity in UWPs.
Following Correia et al. (2006), we define the multiplicity fraction per ultrawide binary (MF_{uw}) as the number of higher multiplicity systems over the total number of UWP systems of the sample with projected component separations typically higher than 1 kAU. Over the whole range of separation we obtain: (12)We similarly define the companion frequency per ultrawide binary (CF_{uw}) as: (13)We use the Poisson error to estimate the error made on MF_{uw} (resp. CF_{uw}), that is, we use the square root of the number of higher order multiple systems (resp. the total number of companions in higher order systems 2 × T + 3 × Q + ...) divided by the total number of systems (B, T, ...). We also estimate the multiplicity fraction and companion frequency per ultrawide binary for UWPs separated by less than 60 kAU (resp. 10 kAU): MF_{60} = 65.79 ± 13.16% and CF_{60} = 1.79 ± 0.22 (resp. MF_{10} = 80.95 ± 19.63% and CF_{10} = 2.33 ± 0.33).
The values we obtained are far above the values of multiplicity fraction (26.8 ± 8.1%) and companion frequency (0.68 ± 0.13) per wide binary found in young T Tauri wide binaries in the range 14–17 AU up to 1.7–2.3 kAU (Correia et al. 2006), highlighting the ubiquity of highorder ultrawide binaries in Taurus found in our sample.
Fig. 11 Mean separation of UWPs (sample S_{3}) and standard errorbars (±) as a function of their multiplicity n_{∗} = 2,3,... The separation is consistent with a geometric progression δ_{P} ∝ a^{− n∗} with a ~ 1.8 (shown as the blue line). 
Based on the previous subsection, and since tighter pairs are mostly composed by MM pairs, we expect them to have a higher multiplicity. Indeed, the mean projected separation δ_{P} of UWPs is negatively correlated with the degree of multiplicity (see Fig. 11). We obtain the best geometric progression fit as: (14)where is the mean separation of UWPs per bin of degree of multiplicity n_{∗}. Thus, UWPs composed of two single stars are wider on average than UWPs composed of a single star and a binary, which in turn are wider on average than UWPs composed of two multiple systems. In summary, tighter UWPs are biased towards higherorder multiplicity.
4.7. Ultrawide pairs and mass function
The mass function of single stars and primaries of multiple systems within UWPs depends on their type (see right Fig. 12). A KS test indicates that the mass distribution for SM pairs cannot be distinguished from that of either the SS or the MM pairs. However, a similar KS test shows that the difference in mass distribution between SS and MM pairs is significant with a confidence level of 97.8%.
Specifically, compared to SS pairs, MM pairs exhibit:

a deficit of mass below 0.1 M_{⊙}, with a frequency of only 5% (against 20%);

a deficit of stars with mass in the 0.1–0.6 M_{⊙} range, with a frequency of 25% for the MM pairs (against 50%);

an excess of stars with mass between 0.6 and 0.8 M_{⊙}, with a frequency of 45% (against 25%);

a predominance of stars more massive than 0.8 M_{⊙} (20% of the most massive stars in Taurus are in MM pairs but only 5% are in SS pairs).
Thus, 70% of stars within SS pairs have a mass of less than 0.6 M_{⊙} and ~65% of stars within MM pairs are more massive than 0.6 M_{⊙}.
Fig. 12 Boxplot for the distribution of the (primary) star mass within UWP (left), the sum of the masses of the two (primary) stars within each UWP (center), and the separation (right) of UWPs. The SS, SM and MM types of systems are shown separately to illustrate that MM pairs are tighter and more massive than SS pairs. 
As a result of these trends (see Fig. 12), the median primary mass of MM pairs is twice as high as the median primary mass in SS pairs (0.7 M_{⊙} compared to 0.3 M_{⊙}). Furthermore, they have a much smaller dispersion as indicated by the interquartile range for MM pairs (0.5–0.8 M_{⊙} versus 0.1–0.7 M_{⊙} for SS pairs, see left panel in Fig. 12). The median primary mass of MM pairs is intriguingly close to the unusual peak at 0.8 M_{⊙} in Taurus reported by Briceño et al. (2002). This peak was tentatively explained by a fragmentation/ejection hydrodynamical model (Goodwin et al.2004) of molecular cores that have unusual properties (extended envelope, a narrow range of core mass function with a peak near 5 M_{⊙}, and a low level of turbulence). The results of simulations show that 50% of the lowmass objects that form within the cores are ejected from their cradles to produce an extended population of lowmass stars and browndwarfs, whereas the remaining stars in multiple systems accrete the gas reach a mass of up to 0.8–1 M_{⊙}. This type of dynamical ejection model was also advocated by Reipurth & Clarke (2001), Boss (2001), Bate et al. (2003), Kroupa & Bouvier (2003a), Kroupa et al. (2003). However, Briceño et al. (2002) and Luhman (2006) argued that the absence of a widely dispersed browndwarf population is strong enough evidence to reject this type of model. Therefore up to now, there is no clear consensus to explain the peak at 0.8 M_{⊙} in the mass function of Taurus.
Fig. 13 Primary mass versus mass ratio within MM (filled dark blue circle), SS (filled light red circle), and SM (open light blue circle) UWPs. Intermediate (resp. smallest) size of marks: UWPs separated by less than 35 kAU (resp. less than 10 kA), the largest size devoted to wider UWPs beyond 35 kAU. Dashed lines: brown dwarfs (BD)/planet (0.012 M_{⊙}) and BD/star (0.075 M_{⊙}) limits. Dotted line: 0.8 M_{⊙} star limit. The superpositions are due to the discretized nature of mass determination models from the spectral type and evolutionary tracks. 
Summing the two (primary) masses within each UWP amplifies the contrast between the lower mass threshold needed to produce MM pairs with respect to SS pairs (see the middle panel in Fig. 12). The median sum of primary masses within MM pairs (1.5 M_{⊙}) is twice as high as that within SS pairs (0.7 M_{⊙}), but the width of the interquartile ranges are comparable (1.2–1.7 M_{⊙} and 0.4–0.8 M_{⊙} for MM and SS pairs, respectively). The difference in mass range between SS and MM pairs would be even more pronounced if the mass of inner companions within MM pairs were also included in the sum; unfortunately, the current multiplicity data are insufficient to evaluate this in a consistent manner across the entire sample. In summary, the median of the stellar mass in MM pairs is twice as high as that within SS pairs, and they are also ten times less divergent; the masses of stars within SS pairs spread a larger range unlike the primary masses in MM pairs. The SM pairs have intermediate properties between these two classes.
4.8. Mass ratio in UWPs
The study of the mass ratio of the secondary mass over the primary mass in the UWPs as a function of the primary mass shows several interesting features (see Fig. 13). First the secondary mass covers the whole range of mass, from equal mass down to brown dwarf type mass (0.075 M_{⊙}) suggesting the same continuous formation scenario from star to brown dwarf. Secondly, there is a clumping trend of MM UWPs towards the upper part of the diagram. Further, we note that there are no UWPs composed of two brown dwarfs.
Fig. 14 Correlation between multiplicity and primary (large size mark) & secondary (small size mark) mass of stars in MM (filled dark blue triangle), SS (filled light red triangle), and SM (light blue open triangle) UWPs. 
Some other features appear to be spurious. The apparent concentration of secondary brown dwarfs pairing with 0.6–0.8 M_{⊙} stars showing up in Fig. 13 is in fact due to the preeminence of these 0.6–0.8 M_{⊙} stars in the mass function. By performing 10 000 Monte Carlo mass pairing samplings amongst the mass sample of UWP stars, we obtain indeed a pvalue of 0.15 compatible with the hypothesis that this pattern is due to chance pairing. The next intriguing pattern that appears at first glance is the empty top right “wedge”, where no secondary stars heavier than 0.8 M_{⊙} lie in for primary stars heavier than 0.8 M_{⊙}. By performing a similar Monte Carlo sampling, we obtain a pvalue of 0.005, which seems to indicate that we can apparently reject the hypothesis that this occurs by chance. This result, however, is not based on a large sample size (median expectation: four pairs within that wedge), and the stellar mass determination from theoretical evolutionary tracks, as well as the spectral type classification, artificially stretches the width of this empty wedge. At this stage, we cannot therefore ensure that it is a reliable pattern.
The increase of the multiplicity goes along with the primary mass of the UWPs (see Fig. 14). This trend was also noted for closer binaries (Kraus & Hillenbrand 2007). The mass of the secondary star does not correlate with the multiplicity to the same extent.
5. Discussion
In this section, we discuss the possible nature of the UWPs identified in Taurus.
5.1. UWPs as physical pairs
The findings surrounding UWPs obtained in the previous section amount to an array of circumstantial evidence suggesting that these pairs are probably physical pairs:

nonrandom 1NNS distribution of mutual pairs (versuslognormal 1NNS distribution for nonmutual pairs);

approximate loguniform separation distribution (Öpik’s law distribution) for mutual pairs, extending up to ~60 kAU (i.e., 0.27 pc) at least, observed for binaries at closer separations (Kraus & Hillenbrand 2009b);

nonrandom Class pairings within mutual pairs suggesting coevality;

UWPs in the range of separation <5 kAU are already known to be true physical binaries (Kraus & Hillenbrand 2009b).
To be gravitationally bound, the internal energy of UWPs must be negative, that is, 1 / 2μΔv^{2}−GM_{1}M_{2}/r < 0, with μ being the reduced mass (μ = M_{1}M_{2}/ (M_{1} + M_{2}), M_{1} and M_{2} being the mass of the stars), G the gravitational constant, r and Δv the separation and the relative velocity between the two stars, respectively. Thus, the condition to meet for a pair to be bound can be rewritten as follows: (15)This relative velocity between the stars is low compared to the observed already low velocity dispersion in the Taurus stellar complex; 2–4 km s^{1} in the whole region (Frink et al. 1997; Rivera et al. 2015) and as low as 1–2 km s^{1} in stellar subgroups (Jones & Herbig 1979). However, as we will see in the following subsection, the low stellar density of Taurus and the low stellar velocity dispersion prevent efficient exchanges of energy during random encounters between wide pairs and other stars, such that the characteristic disruption time τ_{d} of a few ~10^{3} Myr is much longer than the age of the Taurus premain sequence population, that is, 1–10 Myr old (Kenyon & Hartmann 1995; Bertout et al. 2007; Andrews et al. 2013). Therefore, we consider in the following that UWPs are physically linked. This hypothesis will be further tested by results of the Gaia survey that will provide 3D spatial data and 2D/3D velocity measurements (proper motions of stars and radial velocity for part of them).
5.2. A short review on wide binaries: models and observations
From the theoretical side, different models and numerical simulations have been performed to produce multiple systems (see Reipurth et al. 2014, for a complete review on recent numerical simulation improvements). One notable result from simulations of star formation from collapsing gas is that they do not generate wide pairs at birth beyond 1 kAU because of the strong dynamical interactions that occur in simulations of the collapse of turbulent cores (Bate 2012).
Nonetheless, wide pairs are present among field stars, albeit in relatively small numbers: approximately 10% of all solartype field stars have a companion with separations larger than 1 kAU (Raghavan et al. 2010; Lépine & Bongiorno 2007; Longhitano & Binggeli 2010) and maximum separations as high as 20 kAU (Chanamé & Gould 2004; Shaya & Olling 2011; Dhital et al. 2015).
The lowbutnotnegligible fraction of visual wide pairs in the field could not solely be the result of random gravitational capture: the binary creation rate in the field is estimated to be of the order 4 × 10^{21} pc^{3} Gyr^{1} (Kouwenhoven et al. 2010) showing that random creation of bound pairs in the field, as well as in the starforming regions, is an extremely rare event.
An alternate mechanism must be at play to produce such pairs either at birth or as the result of subsequent dynamical evolution. But wide pairs with separations beyond 0.1 pc appear to be particularly unlikely pristine products of star formation process because their separations exceed the typical size of a collapsing cloud core, precluding a scenario based exclusively on core fragmentation.
Consequently, two scenarios based on dynamical evolution have been proposed to explain their origin. One possibility is that an initially compact triple system could dynamically unfold, with the tighter pair shrinking while the third component being ejected onto an extremely long orbit (Reipurth & Mikkola 2012). Alternatively, the gradual dissolution of a star cluster could leave behind wide binaries that are barely bound if two stars happen to be very close in the 6D phase space (Kouwenhoven et al. 2010; Moeckel & Bate 2010; Moeckel & Clarke 2011). For most star clusters, the dissolution time, when this “pairing” could occur, is of the order of a few Gyr (Lamers & Gieles 2006).
The model of Kouwenhoven et al. (2010) also predicts that these wide pairs are highorder multiples as a reflection of the high multiplicity of individual stars at birth. This matches nicely with our findings on UWPs in Taurus, however their mechanism also induces random pairing with respect to mass and Class because the binding pairs are generated by chance depending on the proximity of the objects once the gas of the cluster is removed. This is not what is observed within UWPs in Taurus, since we have shown in the previous section that the Classpairing is not random and that the distributions of mass is not the same in SS and MM pairs.
Another inadequacy of both the “pairing during cluster dissolution” theory and the “unfolding of triple systems” theory is that they occur on timescales that are orders of magnitude longer than the age of the Taurus population (0.1–0.5 Myr for Class I stars and a few Myr for Class II and III stars).
The large fraction of UWPs in Taurus is reminiscent of the surprisingly large number of UWPs (10–50 kAU) recently identified in slightly older moving groups (10–100 Myr) such as in the β Pic moving group (AlonsoFloriano et al.2015). A preference has also been identified for highorder multiplicity among these UWPs, as well as in AB Doradus, TW Hydrae and TucanaeHorologium moving groups (Torres et al. 2006; Elliott et al. 2016). Elliott & Bayo (2016) suggest that the binary population from close (0.1 AU) to very wide systems (100 kAU) in β Pic Moving Group can be accounted for by the internal dynamical interaction of triple systems as proposed by Reipurth & Mikkola (2012). But, as an alternative, these very wide pairs may be the survival part of the wide pairs population formed at an earlier time.
Some of these wide pairs may even survive longer since solartype stars within highorder multiplicity UWPs with periods that peak around 100 yr (separations about 36.5 kAU) have also been observed in the field (Tokovinin 2014) as well as for latespectraltype M1–M5 dwarf UWPs (Law et al. 2010). In the latter case, a very high value was found for the multiple fraction per ultrawide pair, that is, % for systems with separations >4 kAU, very close to what we obtain for the highordermultiple fraction derived in Taurus UWPs with separation less than 10 kAU.
On the other hand very wide pairs composed of extremely young (~0.1 Myr) Class 0 objects are relatively common (Looney et al. 2000; Tobin et al. 2010, 2016; van Kempen et al. 2012; Chen et al. 2013; Lee et al. 2015; Pineda et al. 2015). We thus conclude that the UWPs in Taurus must have been produced by a mechanism that acts extremely early on and, indeed, that these systems may well be the pristine outcome of star formation itself, as is suggested moreover by the coevality trend.
We contend that the population of UWPs that we identified in Taurus supports the idea that very wide binary systems are almost universally produced in starforming regions, even though the physical mechanism leading to their formation remains under debate.
5.3. A pristine origin for these UWPs
To argue for pristine UWPs, one must ensure that they can survive dynamical interactions. The destruction of multiple stars occurs either due to the intrinsic instability of the nonhierarchical multiple stellar system itself or due to the decay of a few Nbody systems driven by hard or soft encounters with other stars of that region. For the former mechanism, the timescale is very short, that is, a few crossing times, or approximately a few 10 kyr (Reipurth 2000), far less than the Taurus age (few Myr).
Disruption and dissolution of wide binary stars in the Galaxy by gravitational encounters with other stars, molecular clouds, and tidal forces due to gravitational potential of the Galaxy is a longstanding study (Chandrasekhar 1944; Heggie 1975, 1977; Retterer & King 1982; Bahcall et al. 1985; Jiang & Tremaine 2010). It was estimated that in the solar neighborhood, the halflife of a binary composed of two solartype stars separated by 31 kAU is 10 Gyr (Jiang & Tremaine 2010).
We expect that in a lowstellardensity region such as Taurus (i.e., 1–10 pc^{2}), dynamical interactions between stellar systems are rare and soft. We thus expect minor dynamical destruction of the wide pairs in Taurus, such that most probably wide binaries in Taurus reveal their pristine configuration. The following arguments on characteristic times associated with the destruction of binaries due to dynamical encounters will give some support to these expectations.
Binney & Tremaine (2008) have estimated the effects of dynamical encounters of binaries with random stars, allowing us to derive two destruction timescales for wide binaries. First, we consider catastrophic single encounters, characterized by a timescale (16)Secondly, we evaluate the destruction timescale in the diffusive regime: (17)In both cases, a is the semimajor axis, M_{tot} is the total mass of the pair, ρ_{c} is the stellar volumic density, M_{c} is the mass of the random encounter star, σ_{rel} is the stellar dispersion velocity, and k_{diff} = 0.085(b_{min}/a)^{2} is a diffusive parameter with b_{min} ~ a being the impact parameter cut off.
While the semimajor axis a cannot be estimated directly for these long period binaries, it is statistically correlated with the projected separation δ, as δ ≈ a (Leinert et al. 1993; Tokovinin & Lépine 2012). The closeness of the two values depends on the chosen eccentricity distribution type, either flat (f(e) = const.), as observed in solartype stars from the Hipparcos catalog by Raghavan et al. (2010) for wide binaries, or thermalized (f(e) = 2e), as first theoretically proposed by Jeans (1919), generalized to broader types of distribution by Ambartsumian (1937) and observed in field stars by Duquennoy & Mayor (1991).
We evaluate the volume density through ρ_{c} = ρ_{w}/ Δ ~ 1 / 25 ~ 0.04 pc^{3}, with ρ_{w} being the mean projected surface density (Eq. (1)) and Δ ~ 20 pc the depth of Taurus (Torres et al. 2007).
Both destruction times are much larger than the age of the young Taurus starforming region. Therefore, in such an environment, it is extremely unlikely that wide binaries, even with separation as large as 100 kAU, can be destroyed by dynamical encounters. Numerical simulations have indeed confirmed that loose associations such as Taurus provide an environment that is inefficient at binary disruption even at large separation (Kroupa & Bouvier 2003b; Marks & Kroupa 2012). Even if they use an initial distribution of binary separation up to 10 kAU, we expect that a same result would be obtained with an extended binary separation up to 100 kAU.
We may thus feel confident that the very wide pairs populations in Taurus may trace the initial spatial distribution at birth without being destroyed by dynamical encounters up to now.
Fig. 15 Schematics for a prolate core with a length d in the main direction and a width of d/n_{c} giving birth to two multiple systems whose primaries are separated by δ. 
5.4. Core properties for pristine UWPs
In the context of a pristine configuration for the wide pairs, we look for the initial conditions that could lead to the formation of such pairs in a core fragmentation model (see Fig. 15).
We thus make the hypotheses that; (1) these UWPs are pristine remnants of the starformation process; and (2) they form within a single core with a coretostar efficiency factor α_{C} ~ 40% (Könyves et al. 2015). The total mass of parental clump can be evaluated as M_{tot}(1 + α_{B}) /α_{C}, where M_{tot} is the sum of (primary) masses and α_{B}M_{tot} is the mass of inner companion stars (α_{B} = 0 for SS pairs). We also introduce a geometrical shape factor n_{c} to quantify the propension of clumps to be elongated and prolate rather than roundish or oblate (n_{c} ~ 2–4, Myers et al. 1991; Ryden 1996; Lomax et al. 2013). From those assumptions, we get the mean particle density ρ of an initial core as: (18)where μ = 2.33 is the mean molecular weight and m_{H} is the atomic hydrogen mass.
We proceed to compute the median density estimate of initial cores that would be required to produce each type of pair ρ_{SS},ρ_{SM},ρ_{MM}: (19a)(19b)(19c)We then question whether those kinds of cores would be instable or stable against any perturbation. As the Jeans instability length is given by: (20)where c_{s} = 0.2 km s^{1} is the sound speed at T = 10 K, we derive the following Jeans length estimates for the three types of pairs: (21a)(21b)(21c)where δ_{SS},δ_{SM},δ_{MM} are the median separations of the three types of UWPs.
The Jeans length is therefore slightly larger than the characteristic length of the SS fictitious clumps (Eq. (21a)), while it is smaller for the MM and SM ones (Eqs. (21b) and (21c)). In other words, the latter two appear to be instable against gravitational (Jeans) instability, while the former appears to be marginally stable.
Therefore, it is plausible that SM and MM pairs formed within the same core, while the SS pairs may have formed through an alternative scenario. We may still consider a pristine origin for them. The two stars may be born in two distinct but proximate cores. Or alternatively, the two stars may be born within the same lowdensity core that becomes supercritical due to the external potential of the Galaxy (BallesterosParedes et al. 2009b,a) or due to the tidal potential of the parent molecular cloud (Horton et al. 2001), both scenarios triggering enhanced fragmentation.
5.5. Towards a fragmentation cascade scenario
We note a negative correlation between the multiplicity within the UWPs and their separation (see Fig. 11): the multiplicity decreases with wider separation. The fragmentation of dense molecular cores into UWPs with at least one highmass (≳0.5M_{⊙}) component seems to have an impact on the probability that either one component, and most probably both, will further fragment. The extent of this fragmentation cascade would depend on the separation of the first two fragments. It suggest a competitive scenario between collapse and an efficient fragmentation process of the gas core/clump to create multiple systems rather than one single object. The densest and most massive molecular cores in Taurus would produce high multiplicity hierarchical MM UWPs. Based on the work of Kraus et al. (2011), this fragmentation cascade scenario may be inhibited at lower spatial scales, since they found that there is no relation between a 1–5 kAUwide binary and innermost highorder multiplicity distribution.
The probable primordial nature of UWPs suggests a scenario in which the turbulent fragmentation of a molecular filament forms overdensities (Pineda et al. 2015) that undergo a fragmentation cascade, producing wide pairs that are themselves closer multiple systems. These systems may be stable enough to survive 1–30 Myr in a lowdensity environment. This link between multiple systems and wide pairs has also been outlined by millimeter observations that show that fragmentation of clumps can proceed through separate envelopes or inside a common envelope (Looney et al. 2000).
5.6. Distributed and clustered modes of formation clues in Taurus
Taurus is considered as the archetypical starforming region that gives rise to isolated prestellar cores, with a typical density of 10^{5} cm^{3} and a size of 0.1 pc. The consensus model for star formation therefore predicts that this cloud will produce a distributed star population, in contrast with the clustered starformation mode that is associated with denser regions such as ρ Ophiuchus, with typical core densities of 10^{7} cm^{3} and sizes of 0.02–0.03 pc (WardThompson et al. 2007). In this context, it is worth noting the significant (more than two orders of magnitude) difference between the typical densities derived for the hypothesized cores that give birth to the SM pairs and MM pairs. These results suggest that, despite its overall low core densities, the Taurus molecular cloud may harbor the two modes of star formation, isolated and clustered.
6. Conclusion
In the Taurus starforming complex, we have identified an extended binary regime up to ~60 kAU using the one point correlation function, Ψ, which further allowed us to distinguish three spatial regimes. The clustering regime has an upper clustering threshold of 0.1° (0.24 pc, i.e., 50 kAU, at the distance of Taurus). The Ψ function appears to be scalefree (ψ = r^{1.53}) over three decades and extends the usual binary regime, as defined by the twopoint correlation function, to a very wide pair regime. We note that the value of the upper clustering length coincides with the local maximum of Taurus NH^{2+} molecular cores correlation function reported by Hacar et al. (2013). As the latter indicates a typical separation between the cores, we expect to see a clustering imprint of stars starting at this value, as is observed. This length is more than twice the universal typical width (0.1 pc) of filamentary structures highlighted by Herschel in Taurus and more generally in Gould Belt starforming regions (André et al. 2014). The clustering regime is followed by a regime of stellar inhibition between 0.1° and 0.5° (0.24 pc up to 1.2 pc), and ends beyond this with a third regime associated to highly isolated stars.
We have highlighted a major structural pattern in the spatial distribution of stars in Taurus based on their multiplicity status. Distinguishing single stars from multiple systems, defined as having at least one physical companion within 1 kAU, we highlight a major and unexpected difference in their spatial distribution based on 1NNS study: the stellar neighborhood in the range of 1−10 kAU of multiple stars is more “crowded” than the single star neighborhood: 40% of multiple systems have a companion within that range compared to 15% for single stars. The probability that a multiple system has a wide companion is threefold that of a single star. We have shown that this excess within 10 kAU is due to the very wide pairs candidates that are moreover generally composed themselves of two multiple systems.
We have identified a potential very wide pairs population in Taurus. Our work shows that 55% of our stellar catalog in Taurus are UWPs in the range of 1−60 kAU. These UWPs are composed preferentially of stars of the same Class, departing from random Class pairing, and, thus, pointing towards coeval pairs. The pair fraction of UWPs follows an almost flat Öpik function (r^{0.14}), extending what has been found for the young binaries at lower range. This coevality clue and Öpik law behavior coupled with the fact that their 1NNS distribution differs from random mutual pairs suggests that these may be true physical pairs; indeed, UWPs with separation less than 5 kAU have been found to be wide binaries (Kraus & Hillenbrand 2009b).
Amongst the UWP population in Taurus, we distinguish three different types of UWPs, depending on whether they are composed either of two multiple systems (MM pairs), a single star and a multiple system (SM pairs), or two single stars (SS pairs). Their properties differ in terms of (primary) mass and separation range. The MM pairs are composed of more massive stars (median mass: 0.65 M_{⊙}) and tighter (median value 3.5 kAU, 75% of them having separation below 9 kAU). The median mass of stars within SS pairs is 0.25 M_{⊙} and their median separation is 35 kAU (75% of them have a separation above 12 kAU). SM pairs have intermediate properties (median mass of primaries of 0.5 M_{⊙} and a uniform separation distribution covering the whole range 1−60 kAU). Multiplicity in these UWPs is a decreasing function of stellar primary mass. We also show that the multiplicity within UWPs increases as the separation of (primary) stars decreases, showing an increased highordermultiple fraction for the tightest UWP targets in Taurus. Class
These differences suggest a different scenario for the formation of UWPs in Taurus. The “massive” MM pairs most probably reflect the imprints of the starformation process within the same molecular core. Based on estimates of their hypothesized natal core properties, Taurus may generate both isolated and clustered star formation. The “lowmass” SS pairs could point to a different formation scenario or may result from a dynamical process, such as lowmass star ejection from multiple systems.
Obtaining parallaxes and proper motions for all Taurus members would help tremendously in identifying which of the UWPs are physical binaries, common proper motion pairs and those which may form from a dynamical process. This will ultimately be achieved thanks to the Gaia survey that will determine the astrometric quantities (angular position, proper motion, and parallax) of stars with an accuracy of 20 μ and at a brightness of 15 mag, and determine the radial velocity of stars brighter than 16 mag.
Unfortunately, the first Gaia release (Gaia Collaboration 2016) does not provide enough data on the whole sample of stars in Taurus (only ~5%). We expect that future releases of Gaia, perhaps as early as the second release, will allow us to firmly establish the physical status of the UWPs, in turn providing a means to discriminate between different scenarios for their formation, and more generally will provide invaluable information on the structure and dynamics of Taurus as a young association (Moraux 2016).
In summary, we identified a new category of ultrawide pairs (UWPs) in Taurus outlining a high order of multiplicity for tighter ones (separation less than 10 kAU). We suggest that part of these UWPs may be pristine imprints of their spatial configuration at birth and put forward a cascade fragmentation scenario for their formation. Furthermore, UWPs may constitute a potential link between the wide pairs recently identified in very young Class 0 type objects (≲0.5 Myr) and the somewhat older but still young moving groups (~20–30 Myr).
Acknowledgments
The authors would like to thank Lee Mundy for helpful discussions and Marc Pound for reading the first draft of this article. We also want to sincerely thank the anonymous referee who helped us to improve both the scientific and language aspects of this paper. This work was funded by the French national research agency through ANR 2010 JCJC 05011 DESC (PI E. Moraux).
References
 Ahn, S., & Fessler, A. 2003, Standard Errors of Mean, Variance, and Standard Deviation Estimators, Tech. Rep., Technical Report, Ann Arbor, MI, USA [Google Scholar]
 AlonsoFloriano, F. J., Caballero, J. A., CortésContreras, M., Solano, E., & Montes, D. 2015, A&A, 583, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Alves, J., Lombardi, M., & Lada, C. J. 2007, A&A, 462, L17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ambartsumian, V. A. 1937, Astron. Zh., 14, 207 [Google Scholar]
 Andre, P., WardThompson, D., & Barsony, M. 1993, ApJ, 406, 122 [NASA ADS] [CrossRef] [Google Scholar]
 André, P., Di Francesco, J., WardThompson, D., et al. 2014, Protostars and Planets VI (Tuscon: University of Arizona Press), 27 [Google Scholar]
 Andrews, S. M., Rosenfeld, K. A., Kraus, A. L., & Wilner, D. J. 2013, ApJ, 771, 129 [NASA ADS] [CrossRef] [Google Scholar]
 Baddeley, A., & Turner, R. 2005, J. Statist. Software, 12, 1 [CrossRef] [Google Scholar]
 Bahcall, J. N., Hut, P., & Tremaine, S. 1985, ApJ, 290, 15 [NASA ADS] [CrossRef] [Google Scholar]
 BallesterosParedes, J., Gómez, G. C., Loinard, L., Torres, R. M., & Pichardo, B. 2009a, MNRAS, 395, L81 [NASA ADS] [CrossRef] [Google Scholar]
 BallesterosParedes, J., Gómez, G. C., Pichardo, B., & VázquezSemadeni, E. 2009b, MNRAS, 393, 1563 [NASA ADS] [CrossRef] [Google Scholar]
 Bate, M. R. 2012, MNRAS, 419, 3115 [NASA ADS] [CrossRef] [Google Scholar]
 Bate, M. R., Bonnell, I. A., & Bromm, V. 2003, MNRAS, 339, 577 [NASA ADS] [CrossRef] [Google Scholar]
 Benaglia, T., Chauveau, D., Hunter, D. R., & Young, D. 2009, J. Stat. Software, 32, 1 [CrossRef] [Google Scholar]
 Bertout, C., Siess, L., & Cabrit, S. 2007, A&A, 473, L21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Binney, J., & Tremaine, S. 2008, in Galactic Dynamics, second edition (Princeton University Press) [Google Scholar]
 Boss, A. P. 2001, ApJ, 551, L167 [NASA ADS] [CrossRef] [Google Scholar]
 Briceño, C., Luhman, K. L., Hartmann, L., Stauffer, J. R., & Kirkpatrick, J. D. 2002, ApJ, 580, 317 [NASA ADS] [CrossRef] [Google Scholar]
 Burrows, C. J., Stapelfeldt, K. R., Watson, A. M., et al. 1996, ApJ, 473, 437 [NASA ADS] [CrossRef] [Google Scholar]
 Carney, M. T., Yıldız, U. A., Mottram, J. C., et al. 2016, A&A, 586, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chakraborty, A., & Ge, J. 2004, AJ, 127, 2898 [NASA ADS] [CrossRef] [Google Scholar]
 Chakraborty, A., Feigelson, E. D., & Babu, G. J. 2014, Astrolabe: Astronomy Users Library for R, r package version 0.1 [Google Scholar]
 Chanamé, J., & Gould, A. 2004, ApJ, 601, 289 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1943, Rev. Mod. Phys., 15, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1944, ApJ, 99, 54 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, X., Arce, H. G., Zhang, Q., et al. 2013, ApJ, 768, 110 [NASA ADS] [CrossRef] [Google Scholar]
 Connelley, M. S., Reipurth, B., & Tokunaga, A. T. 2008, AJ, 135, 2496 [NASA ADS] [CrossRef] [Google Scholar]
 Connelley, M. S., Reipurth, B., & Tokunaga, A. T. 2009, AJ, 138, 1193 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, S., Zinnecker, H., Ratzka, T., & Sterzik, M. F. 2006, A&A, 459, 909 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Covey, K. R., Greene, T. P., Doppmann, G. W., & Lada, C. J. 2006, AJ, 131, 512 [NASA ADS] [CrossRef] [Google Scholar]
 Daemgen, S., Bonavita, M., Jayawardhana, R., Lafrenière, D., & Janson, M. 2015, ApJ, 799, 155 [NASA ADS] [CrossRef] [Google Scholar]
 Dahl, D. B. 2014, xtable: Export tables to LaTeX or HTML, r package version 1.74 [Google Scholar]
 Dhital, S., West, A. A., Stassun, K. G., Schluns, K. J., & Massey, A. P. 2015, AJ, 150, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Di Folco, E., Dutrey, A., Le Bouquin, J.B., et al. 2014, A&A, 565, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dobbs, C. L., Krumholz, M. R., BallesterosParedes, J., et al. 2014, Protostars and Planets VI (Tuscon: University of Arizona Press), 3 [Google Scholar]
 Duchêne, G. 1999, A&A, 341, 547 [NASA ADS] [Google Scholar]
 Duchêne, G., & Kraus, A. 2013, ARA&A, 51, 269 [Google Scholar]
 Duchêne, G., Ghez, A. M., McCabe, C., & Weinberger, A. J. 2003, ApJ, 592, 288 [NASA ADS] [CrossRef] [Google Scholar]
 Duchêne, G., Bontemps, S., Bouvier, J., et al. 2007, A&A, 476, 229 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Duquennoy, A., & Mayor, M. 1991, A&A, 248, 485 [NASA ADS] [Google Scholar]
 Elliott, P., & Bayo, A. 2016, MNRAS, 459, 4499 [NASA ADS] [CrossRef] [Google Scholar]
 Elliott, P., Bayo, A., Melo, C. H. F., et al. 2016, A&A, 590, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Esplin, T. L., Luhman, K. L., & Mamajek, E. E. 2014, ApJ, 784, 126 [NASA ADS] [CrossRef] [Google Scholar]
 Frink, S., Röser, S., Neuhäuser, R., & Sterzik, M. F. 1997, A&A, 325, 613 [NASA ADS] [Google Scholar]
 Gaia Collaboration (Brown, A. G. A., et al.) 2016, A&A, 595, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ghez, A. M., Neugebauer, G., & Matthews, K. 1993, AJ, 106, 2005 [NASA ADS] [CrossRef] [Google Scholar]
 Ghez, A. M., White, R. J., & Simon, M. 1997, ApJ, 490, 353 [NASA ADS] [CrossRef] [Google Scholar]
 Gladwin, P. P., Kitsionas, S., Boffin, H. M. J., & Whitworth, A. P. 1999, MNRAS, 302, 305 [NASA ADS] [CrossRef] [Google Scholar]
 Gomez, M., Hartmann, L., Kenyon, S. J., & Hewett, R. 1993, AJ, 105, 1927 [NASA ADS] [CrossRef] [Google Scholar]
 Goodwin, S. P., Whitworth, A. P., & WardThompson, D. 2004, A&A, 419, 543 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Goodwin, S. P., Nutter, D., Kroupa, P., WardThompson, D., & Whitworth, A. P. 2008, A&A, 477, 823 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Graffelman, J. 2013, calibrate: Calibration of Scatterplot and Biplot Axes, r package version 1.7.2 [Google Scholar]
 Greene, T. P., Wilking, B. A., Andre, P., Young, E. T., & Lada, C. J. 1994, ApJ, 434, 614 [NASA ADS] [CrossRef] [Google Scholar]
 Hacar, A., Tafalla, M., Kauffmann, J., & Kovács, A. 2013, A&A, 554, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Harrell, F., E., Dupont, C., & others. 2015, Hmisc: Harrell Miscellaneous, r package version 3.160 [Google Scholar]
 Harris, A. 2013, FITSio: FITS (Flexible Image Transport System) utilities, r package version 2.00 [Google Scholar]
 Hartmann, L. 2002, ApJ, 578, 914 [NASA ADS] [CrossRef] [Google Scholar]
 Heggie, D. C. 1975, MNRAS, 173, 729 [NASA ADS] [CrossRef] [Google Scholar]
 Heggie, D. C. 1977, Rev. Mexicana Astron. Astrofis., 3, 169 [NASA ADS] [Google Scholar]
 Herczeg, G. J., & Hillenbrand, L. A. 2014, ApJ, 786, 97 [Google Scholar]
 Holman, K., Walch, S. K., Goodwin, S. P., & Whitworth, A. P. 2013, MNRAS, 432, 3534 [NASA ADS] [CrossRef] [Google Scholar]
 Horton, A. J., Bate, M. R., & Bonnell, I. A. 2001, MNRAS, 321, 585 [NASA ADS] [CrossRef] [Google Scholar]
 Ireland, M. J., & Kraus, A. L. 2008, ApJ, 678, L59 [NASA ADS] [CrossRef] [Google Scholar]
 Itoh, Y., Hayashi, M., Tamura, M., et al. 2005, ApJ, 620, 984 [NASA ADS] [CrossRef] [Google Scholar]
 Jeans, J. H. 1919, MNRAS, 79, 408 [NASA ADS] [CrossRef] [Google Scholar]
 Jiang, Y.F., & Tremaine, S. 2010, MNRAS, 401, 977 [NASA ADS] [CrossRef] [Google Scholar]
 Jones, B. F., & Herbig, G. H. 1979, AJ, 84, 1872 [NASA ADS] [CrossRef] [Google Scholar]
 Juvela, M., & Montillaud, J. 2016, A&A, 585, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kendall, M., & Stuart, A. 1977, The advanced theory of statistics, 4th edn., Vol. 1: Distribution theory (New York, NY: Macmillan) [Google Scholar]
 Kenyon, S. J., & Hartmann, L. 1995, ApJS, 101, 117 [NASA ADS] [CrossRef] [Google Scholar]
 Kerscher, M., Szapudi, I., & Szalay, A. S. 2000, ApJ, 535, L13 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Kirk, H., & Myers, P. C. 2011, ApJ, 727, 64 [NASA ADS] [CrossRef] [Google Scholar]
 Kohler, R., & Leinert, C. 1998, A&A, 331, 977 [NASA ADS] [Google Scholar]
 Konopacky, Q. M., Ghez, A. M., Rice, E. L., & Duchêne, G. 2007, ApJ, 663, 394 [NASA ADS] [CrossRef] [Google Scholar]
 Könyves, V., André, P., Men’shchikov, A., et al. 2015, A&A, 584, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Koresko, C. D. 2000, ApJ, 531, L147 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Kouwenhoven, M. B. N., Goodwin, S. P., Parker, R. J., et al. 2010, MNRAS, 404, 1835 [NASA ADS] [Google Scholar]
 Kraus, A. L., & Hillenbrand, L. A. 2007, ApJ, 662, 413 [NASA ADS] [CrossRef] [Google Scholar]
 Kraus, A. L., & Hillenbrand, L. A. 2008, ApJ, 686, L111 [NASA ADS] [CrossRef] [Google Scholar]
 Kraus, A. L., & Hillenbrand, L. A. 2009a, ApJ, 704, 531 [Google Scholar]
 Kraus, A. L., & Hillenbrand, L. A. 2009b, ApJ, 703, 1511 [NASA ADS] [CrossRef] [Google Scholar]
 Kraus, A. L., & Hillenbrand, L. A. 2012, ApJ, 757, 141 [NASA ADS] [CrossRef] [Google Scholar]
 Kraus, A. L., White, R. J., & Hillenbrand, L. A. 2006, ApJ, 649, 306 [NASA ADS] [CrossRef] [Google Scholar]
 Kraus, A. L., Ireland, M. J., Martinache, F., & Hillenbrand, L. A. 2011, ApJ, 731, 8 [NASA ADS] [CrossRef] [Google Scholar]
 Kroupa, P., & Bouvier, J. 2003a, MNRAS, 346, 369 [NASA ADS] [CrossRef] [Google Scholar]
 Kroupa, P., & Bouvier, J. 2003b, MNRAS, 346, 343 [NASA ADS] [CrossRef] [Google Scholar]
 Kroupa, P., Bouvier, J., Duchêne, G., & Moraux, E. 2003, MNRAS, 346, 354 [NASA ADS] [CrossRef] [Google Scholar]
 Lada, C. J. 1987, in Star Forming Regions, eds. M. Peimbert, & J. Jugaku, IAU Symp., 115, 1 [Google Scholar]
 Lada, C. J., & Lada, E. A. 2003, ARA&A, 41, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Lada, C. J., & Wilking, B. A. 1984, ApJ, 287, 610 [NASA ADS] [CrossRef] [Google Scholar]
 Lada, C. J., Muench, A. A., Rathborne, J., Alves, J. F., & Lombardi, M. 2008, ApJ, 672, 410 [NASA ADS] [CrossRef] [Google Scholar]
 Lamers, H. J. G. L. M., & Gieles, M. 2006, A&A, 455, L17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Landy, S. D., & Szalay, A. S. 1993, ApJ, 412, 64 [NASA ADS] [CrossRef] [Google Scholar]
 Larson, R. B. 1995, MNRAS, 272, 213 [NASA ADS] [Google Scholar]
 Law, N. M., Dhital, S., Kraus, A., Stassun, K. G., & West, A. A. 2010, ApJ, 720, 1727 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, K. I., Dunham, M. M., Myers, P. C., et al. 2015, ApJ, 814, 114 [NASA ADS] [CrossRef] [Google Scholar]
 Leinert, C., Zinnecker, H., Weitzel, N., et al. 1993, A&A, 278, 129 [NASA ADS] [Google Scholar]
 Leinert, C., Richichi, A., & Haas, M. 1997, A&A, 318, 472 [NASA ADS] [Google Scholar]
 Lépine, S., & Bongiorno, B. 2007, AJ, 133, 889 [NASA ADS] [CrossRef] [Google Scholar]
 Lomax, O., Whitworth, A. P., & Cartwright, A. 2013, MNRAS, 436, 2680 [NASA ADS] [CrossRef] [Google Scholar]
 Longhitano, M., & Binggeli, B. 2010, A&A, 509, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Looney, L. W., Mundy, L. G., & Welch, W. J. 2000, ApJ, 529, 477 [NASA ADS] [CrossRef] [Google Scholar]
 Luhman, K. L. 2004, ApJ, 617, 1216 [NASA ADS] [CrossRef] [Google Scholar]
 Luhman, K. L. 2006, ApJ, 645, 676 [NASA ADS] [CrossRef] [Google Scholar]
 Luhman, K. L., Allen, P. R., Espaillat, C., Hartmann, L., & Calvet, N. 2010, ApJS, 186, 111 [NASA ADS] [CrossRef] [Google Scholar]
 Luhman, K. L., Mamajek, E. E., Shukla, S. J., & Loutrel, N. P. 2017, AJ, 153, 46 [NASA ADS] [CrossRef] [Google Scholar]
 Marks, M., & Kroupa, P. 2012, A&A, 543, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Marsh, K. A., Kirk, J. M., André, P., et al. 2016, MNRAS, 459, 342 [NASA ADS] [CrossRef] [Google Scholar]
 Mathieu, R. D. 1994, ARA&A, 32, 465 [NASA ADS] [CrossRef] [Google Scholar]
 Moeckel, N., & Bate, M. R. 2010, MNRAS, 404, 721 [NASA ADS] [CrossRef] [Google Scholar]
 Moeckel, N., & Clarke, C. J. 2011, MNRAS, 415, 1179 [NASA ADS] [CrossRef] [Google Scholar]
 Monnier, J. D., Tannirkulam, A., Tuthill, P. G., et al. 2008, ApJ, 681, L97 [NASA ADS] [CrossRef] [Google Scholar]
 Moraux, E. 2016, EAS Pub. Ser., 8081, 73 [Google Scholar]
 Myers, P. C., Fuller, G. A., Goodman, A. A., & Benson, P. J. 1991, ApJ, 376, 561 [NASA ADS] [CrossRef] [Google Scholar]
 Nychka, D., Furrer, R., & Sain, S. 2015, Rfields: Tools for Spatial Data, r package version 8.21 [Google Scholar]
 Offner, S. S. R., Clark, P. C., Hennebelle, P., et al. 2014, Protostars and Planets VI (Tuscon: University of Arizona Press), 53 [Google Scholar]
 Padgett, D. L., Brandner, W., Stapelfeldt, K. R., et al. 1999, AJ, 117, 1490 [NASA ADS] [CrossRef] [Google Scholar]
 Parker, R. J., Goodwin, S. P., Kroupa, P., & Kouwenhoven, M. B. N. 2009, MNRAS, 397, 1577 [NASA ADS] [CrossRef] [Google Scholar]
 Peebles, P. J. E. 1980, in Ninth Texas Symposium on Relativistic Astrophysics, eds. J. Ehlers, J. J. Perry, & M. Walker, Annals of the New York Academy of Sciences, 336, 161 [Google Scholar]
 Pineda, J. E., Offner, S. S. R., Parker, R. J., et al. 2015, Nature, 518, 213 [NASA ADS] [CrossRef] [Google Scholar]
 R Core Team 2015, R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria [Google Scholar]
 Raghavan, D., McAlister, H. A., Henry, T. J., et al. 2010, ApJS, 190, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rebull, L. M., Padgett, D. L., McCabe, C.E., et al. 2010, ApJS, 186, 259 [NASA ADS] [CrossRef] [Google Scholar]
 Reipurth, B. 2000, AJ, 120, 3177 [CrossRef] [Google Scholar]
 Reipurth, B., & Clarke, C. 2001, AJ, 122, 432 [NASA ADS] [CrossRef] [Google Scholar]
 Reipurth, B., & Mikkola, S. 2012, Nature, 492, 221 [Google Scholar]
 Reipurth, B., & Zinnecker, H. 1993, A&A, 278, 81 [NASA ADS] [Google Scholar]
 Reipurth, B., Clarke, C. J., Boss, A. P., et al. 2014, Protostars and Planets VI (Tuscon: University of Arizona Press), 267 [Google Scholar]
 Retterer, J. M., & King, I. R. 1982, ApJ, 254, 214 [NASA ADS] [CrossRef] [Google Scholar]
 Rivera, J. L., Loinard, L., Dzib, S. A., et al. 2015, ApJ, 807, 119 [NASA ADS] [CrossRef] [Google Scholar]
 Robotham, A. 2015, magicaxis: Pretty Scientific Plotting with MinorTick and log MinorTick Support, r package version 1.9.4 [Google Scholar]
 Ryden, B. S. 1996, ApJ, 471, 822 [NASA ADS] [CrossRef] [Google Scholar]
 Sakhr, J., & Nieminen, J. M. 2006, Phys. Rev. E, 73, 036201 [NASA ADS] [CrossRef] [Google Scholar]
 Sartoretti, P., Brown, R. A., Latham, D. W., & Torres, G. 1998, A&A, 334, 592 [NASA ADS] [Google Scholar]
 Shaya, E. J., & Olling, R. P. 2011, ApJS, 192, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Shu, F. H., Adams, F. C., & Lizano, S. 1987, ARA&A, 25, 23 [Google Scholar]
 Simon, M. 1997, ApJ, 482, L81 [NASA ADS] [CrossRef] [Google Scholar]
 Simon, M., Ghez, A. M., Leinert, C., et al. 1995, ApJ, 443, 625 [NASA ADS] [CrossRef] [Google Scholar]
 Simon, M., Beck, T. L., Greene, T. P., et al. 1999, AJ, 117, 1594 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, K. W., Balega, Y. Y., Duschl, W. J., et al. 2005, A&A, 431, 307 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tafalla, M., & Hacar, A. 2015, A&A, 574, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Therneau, T. 2014, deming: Deming, ThielSen and PassingBablock Regression, r package version 1.01 [Google Scholar]
 Tobin, J. J., Hartmann, L., Looney, L. W., & Chiang, H.F. 2010, ApJ, 712, 1010 [NASA ADS] [CrossRef] [Google Scholar]
 Tobin, J. J., Looney, L. W., Li, Z.Y., et al. 2016, ApJ, 818, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Todorov, K., Luhman, K. L., & McLeod, K. K. 2010, ApJ, 714, L84 [NASA ADS] [CrossRef] [Google Scholar]
 Todorov, K. O., Luhman, K. L., Konopacky, Q. M., et al. 2014, ApJ, 788, 40 [NASA ADS] [CrossRef] [Google Scholar]
 Tokovinin, A. 2014, AJ, 147, 87 [NASA ADS] [CrossRef] [Google Scholar]
 Tokovinin, A., & Lépine, S. 2012, AJ, 144, 102 [NASA ADS] [CrossRef] [Google Scholar]
 Torres, C. A. O., Quast, G. R., da Silva, L., et al. 2006, A&A, 460, 695 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Torres, R. M., Loinard, L., Mioduszewski, A. J., & Rodríguez, L. F. 2007, ApJ, 671, 1813 [NASA ADS] [CrossRef] [Google Scholar]
 van Kempen, T. A., Longmore, S. N., Johnstone, D., Pillai, T., & Fuente, A. 2012, ApJ, 751, 137 [NASA ADS] [CrossRef] [Google Scholar]
 Vincenty, T. 1975, Survey Review, 22, 88 [CrossRef] [Google Scholar]
 WardThompson, D., André, P., Crutcher, R., et al. 2007, Protostars and Planets V (Tuscon: University of Arizona Press), 33 [Google Scholar]
 White, R. J., & Ghez, A. M. 2001, ApJ, 556, 265 [Google Scholar]
 Ysard, N., Abergel, A., Ristorcelli, I., et al. 2013, A&A, 559, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: Spherical geometry
As a reminder, this section gives the formulae to obtain the angular distance Δσ between two stars on the celestial sphere, given by the spherical law of cosines: (A.1)where φ_{i},λ_{i} (i = 1,2) are, respectively, the declination and right ascension of the star i, and Δδ =  δ_{2}−δ_{1}  their absolute difference in declination. Since this above equation is illconditioned for small Δσ, it is better to use the following Vincenty formula (Vincenty 1975) applied to a sphere (A.2)
Appendix B: Analytics for kNearest Neighborhood statistics
Appendix B.1: Theoretical probability density function (PDF) and cumulative distribution of 1nearest neighbor separation (1NNS) for spatial random distribution
In this section, we derive the theoretical distribution of 1NNS for a random point process in an infinite medium, following the work done in 3D by Chandrasekhar (1943). The probability w(r)dr that the first neighbor of a star is located at a distance between r and r + dr in 2D, is given by the product of the probability that there is no stars in a disk of radius r multiplied by the probability of having stars in a shell area between r and r + dr: (B.1)where ρ is the mean intensity of the process in an infinite medium, that is, the average number of stars per unit of surface. Equation (B.1) may be written as: (B.2)Once integrated, from the equation above, we get the probability density function (PDF) of the first neighbor distance for a random Poissonnian process: (B.3)The PDF integral is unity (i.e., we have ).
Furthermore, from the 1NNS PDF (Eq. (B.3)), we get the related cumulative distribution function W(r) as: (B.4)that allows us to define the characteristic value r_{q} associated to the quantile q from: , to get: (B.5)
Appendix B.2: Moments of first nearest neighbor distribution
Appendix B.2.1: Mean
The theoretical mean of the first nearest neighbor distance in an infinite medium is computed from the PDF as: (B.6)where =E(r) is the expected value of r, and erf(x) is the error function. Taking the value at 0 and infinity, we end up with a simple analytical result of first nearest neighbor distance theoretical mean in a random process of intensity ρ in an infinite medium, i.e., is strictly equal to half of the inverse square root of the density: (B.7)
Appendix B.2.2: Variance
From the PDF w(r), we now compute the theoretical variance (σ_{t} being the standard deviation) of the first nearest neighbor distance for a random Poisson process, (B.8)The integration gives: (B.9)We then compute series expansions of the integral at x = 0: (B.10)and then at x = + ∞: (B.11)to finally get: (B.12)
Appendix B.2.3: Skewness
The skewness of a distribution quantifies its asymmetry, and is defined from the third moment of the distribution. When it is positive, the peak of the distribution is shifted to the left, that is, there is a tail towards the right. When it is negative, it is the other way round. The skewness may be defined from the moments of the distribution: (B.13)We note that the skewness of the first nearest neighbor distribution of a random spatial distribution is constant, that is, independent of the density/intensity of the random process.
Appendix B.3: Log PDF
Based on transformation rules, we then derive the theoretical PDF for the Log(1NNS) in a random process. Let y = Log(r), where r is the 1NNS variable, and thus we have: (B.14)where w(r) is the PDF of the 1NNS variable (Eq. (B.3)) and w(y) is the PDF of the logarithm variable. We get the following moments for the random k1NNS theoretical distributions: (B.15)From this relation, we get the PDF of the first neighbor distance logarithm for a random process as: (B.16)where it reaches its maximun value (mode) at . From Taylor expansion for x ≪ 1, we get: (B.17)as we expect in 2D for a random distribution, the first nearest neighbor 1NNS distribution function increases as a squared power law for . From Eq. (B.16) above, we derive its logarithmic expression: (B.18)
In a Poisson point process in a ddimensional space with intensity ρ, the distance r between a point and its kth neighbor is distributed according to the generalized Gamma distribution: (B.19)where V_{d}(r) is the volume ball or radius r in dspace. (B.20)and c_{d} is the unit volume ball in dspace given by: (B.21)and Γ is the Gamma function. So, in a planar (2D) distribution, we get: (B.22)The cumulative distribution is then: (B.23)So for the second nearest neighbor cumulative distribution we get: (B.24)and for the third nearest neighbor cumulative distribution, we have: (B.25)
Appendix C: Catalogs
This section gathers the catalogs described respectively in Sects. 2 and 4.1 of the paper.
Taurus stars catalog.
Catalog of the first nearest neighbor couples of stars.
All Tables
All Figures
Fig. 1 Spatial distribution of stars within the Taurus region superimposed on nearinfrared extinction (unit mag, gray scale) from Juvela & Montillaud (2016). Dashed polygon: restricted spatial window W_{in} in Taurus complex; blue filled dots: random sampling. Colorcoding for Taurus stars: red plus marks: clustered stars [1NNS ≤ 0.1°]; black plus marks: “inhibited” regime stars [0.1°< 1NNS ≤ 0.55°]; green plus marks: isolated stars [1NNS beyond 0.55°] (see Sect. 3). 

In the text 
Fig. 2 Empirical cumulative distribution of the 1NNS distribution for the whole region (black) and for the central part within W_{in} region (blue). The blue vertical line represents the clustering length, r_{c} = 0.1°. 

In the text 
Fig. 3 Probability distribution function (PDF) of the 1NNS distributions observed in Taurus (solid blue histogram: within the limited window W_{in}, that is, the three main filaments; dashed blue histogram: the entire Taurus region) and predicted for a random distribution (black histogram: for a Monte Carlo sampling within the window W_{in}; red curve: unlimited theoretical random 1NNS distribution, see Eq. (3)). Both random populations are drawn using the same mean surface density ρ_{w} = 5deg^{2}. Increasing the mean density by a factor of two leads to a shift towards the left whose amplitude is represented by the green arrow. The Monte Carlo sampling in a finite window and the unlimited theoretical random 1NNS do not significantly differ. The 1NNS distributions computed either in the whole Taurus or within W_{in} do not significantly differ either. However, the Taurus and random 1NNS distributions are statistically highly inconsistent, with a pvalue of ~10^{16} for the KS test. The spatial distribution in Taurus is clearly far from being random. 

In the text 
Fig. 4 Onepoint correlation function Ψ (Eq. (7), blue symbols) and estimated pair correlation function g (gray symbols) using the sample of stars within W_{in}. The pair correlation function has been estimated using the estimator given in Eq. (6) while performing 10 000 random Monte Carlo samplings within W_{in}. The 1σ vertical errorbars of Ψ and g are estimated from Poisson statistics. The solid horizontal blue line represents both the constant Ψ = 1 and ĝ = 1 functions, that is, the onepoint and pair correlation function expected for a spatial random distribution. Vertical black solid lines delimit three spatial regimes, clustering, inhibition, and dispersion, derived from the crossing of the Taurus Ψ function with the Ψ = 1 horizontal line associated to random spatial distribution. The vertical dashed purple line indicates the 7′′ spacing threshold used to group stars into a unique multiple system. The dotteddashed green line is the Ψ ∝ r^{1.5} function resulting from a linear regression fit taking into account vertical errorbars over ten 1NNS logarithm bins (Δlog r = 0.2). The Ψ function may reveal an extended binary regime from 1.6 kAU to 100 kAU in Taurus, which has remained unidentified beyond approximately 5 kAU in previous works that used the pair correlation function alone, since the latter reveals a power law break around that separation. 

In the text 
Fig. 5 First nearest neighbor separation fraction distribution of multiple systems (solid blue histogram) versus single stars (dashed black histogram) for Class II and III objects observed at high angular resolution (HAR) in Taurus region. Each bin represents the number density of stars per unit logarithmic interval in projected separation, that is, the number density fraction of 1NNS per (Δlog r = 0.2) interval. There is a marked excess of 1NNS in the range 1−10 kAU for the multiple systems with respect to single stars. 

In the text 
Fig. 6 Distribution of ultrawide pair fraction as a function of separation in the three main Taurus filaments (i.e., for stars enclosed within W_{in}, solid blue histogram) versus pair fraction of random mutual pairs (obtained from 10 000 Monte Carlo samplings, thick black histogram) and 1NNS fraction of nonmutual pairs (dashed blue histogram). The solid blue line represent a lognormal fit to the distribution of nonmutual pairs in Taurus. The distribution of mutual pairs of Taurus does not result from a random process (pvalue of 10^{16}), and is statistically different from nonmutual pairs distribution in Taurus (pvalue of 10^{14}). Conversely, the 1NNS distribution of nonmutual pairs in Taurus is reasonably well fitted (pvalue of 0.2) by a longnormal function, with a mean value μ = 4.74 ± 0.04 and a standard deviation σ = 0.46 ± 0.03. 

In the text 
Fig. 7 Ultrawide pair fraction F within the whole Taurus region (sample S_{1}) as a function of separation, plotted using Δδ_{P} = 0.25 dex bins. The solid blue line is a power law fit to the data (, with log α = −1.21 and β = 0.14) between 3 and ~60 kAU (vertical black dashed line). The gray area represents the 95% confidence band. In terms of the power law parameters, the corresponding parameter 95% confidence intervals are log α = [ −1.78,−0.65 ], β = [ 0,0.28 ]. 

In the text 
Fig. 8 Distributions of SED Class pairings expected from a random pairing process based on the global proportions of Class I, II, and III stars in Taurus (left) and observed among Taurus UWPs (right). The SED Class pairing within UWPs deviates significantly from random pairing, suggesting coevality within these systems. 

In the text 
Fig. 9 Multiplicity fraction of Class II and III stars observed at HAR and members of UWPs as a function of their separation. For a given value of δ, the black (resp. blue) symbol represents the multiplicity fraction computed over all separations smaller (resp. greater) than δ. Errorbars are standard errors computed from the standard formula: where n_{i} is the sample size, and N is the total population size of Class II and III stars observed at HAR within UWPs independently of their separation. The dashed black line represents the mean multiplicity fraction for stars within UWPs whereas the gray line marks the corresponding mean multiplicity fraction for stars that are not in UWPs. The red dotted lines define the interval within which the multiplicity fraction within UWPs is significantly (pvalue of 0.05 or less based on a chisquared test) higher than the mean value. 

In the text 
Fig. 10 SS (dashed blue), SM (dashed red) and MM (solid blue) pair fraction as a function of their separation. MM pairs are tighter (most are found below 10 kAU) whereas SS pairs are far apart (with a peak at approximately 50 kAU); on the other hand, SM pairs are relatively uniformly distributed. 

In the text 
Fig. 11 Mean separation of UWPs (sample S_{3}) and standard errorbars (±) as a function of their multiplicity n_{∗} = 2,3,... The separation is consistent with a geometric progression δ_{P} ∝ a^{− n∗} with a ~ 1.8 (shown as the blue line). 

In the text 
Fig. 12 Boxplot for the distribution of the (primary) star mass within UWP (left), the sum of the masses of the two (primary) stars within each UWP (center), and the separation (right) of UWPs. The SS, SM and MM types of systems are shown separately to illustrate that MM pairs are tighter and more massive than SS pairs. 

In the text 
Fig. 13 Primary mass versus mass ratio within MM (filled dark blue circle), SS (filled light red circle), and SM (open light blue circle) UWPs. Intermediate (resp. smallest) size of marks: UWPs separated by less than 35 kAU (resp. less than 10 kA), the largest size devoted to wider UWPs beyond 35 kAU. Dashed lines: brown dwarfs (BD)/planet (0.012 M_{⊙}) and BD/star (0.075 M_{⊙}) limits. Dotted line: 0.8 M_{⊙} star limit. The superpositions are due to the discretized nature of mass determination models from the spectral type and evolutionary tracks. 

In the text 
Fig. 14 Correlation between multiplicity and primary (large size mark) & secondary (small size mark) mass of stars in MM (filled dark blue triangle), SS (filled light red triangle), and SM (light blue open triangle) UWPs. 

In the text 
Fig. 15 Schematics for a prolate core with a length d in the main direction and a width of d/n_{c} giving birth to two multiple systems whose primaries are separated by δ. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.