© ESO, 2017

1. Introduction

The formation of stars results fundamentally from the runaway unbalance between overwhelming inwards self-gravity and outwards pressure-gradient 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 low-mass 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 length-scale 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 low-mass 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 (Ward-Thompson et al. 2007). The long-standing scenario proposed so far to form one low-mass 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 star-formation 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 core-to-stars 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 star-forming 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 system-core 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 (Ward-Thompson 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 semi-major 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 two-point correlation function and the related mean surface density of companions to examine the scale-free 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 first-nearest-neighbor separation (1-NNS) 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 (1-NNS distribution and two-point correlation function) to analyse the spatial distribution of stars in Taurus with a focus on multiples with respect to single stars, while introducing the one-point correlation function Ψ (Sect. 3). We then assess the statistical properties of 1-NNS couples while defining ultra-wide 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 mid-infrared Wide-field 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 low-mass 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 star-forming 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 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). 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 age-based 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.

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 higher-angular-resolution 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 Spitzer-MIPS 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 chi-squared 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 South-West (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 nearest-neighbor 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.

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 (1-NNS) statistics. Nearest-neighbor 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 1-NNS statistics to advocate the use of the theoretical 1-NNS 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 well-populated. We then define the one-point correlation function Ψ, which complements the two-point 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. 1-NNS study in Taurus

This subsection aims at analyzing the spatial distribution in Taurus based on a 1-NNS analysis. The 1-NNS of a star is, by definition, the distance to its nearest neighbor star. Throughout this paper, we deal with the projected 1-NNS, 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 1-NNS distribution in the three main filaments to that associated to the entire region. We then compare the theoretical random 1-NNS distribution obtained for an infinite random spatial process to the 1-NNS distribution derived from Monte Carlo samplings in a finite irregularly shaped window, and finally we compare Taurus 1-NNS distribution to the 1-NNS random distribution.

3.1.1. Mean stellar surface density in Taurus

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: ${\mathit{\rho }}_{\mathrm{w}}\mathrm{=}{\mathit{N}}_{\mathrm{w}}\mathit{/}{\mathit{A}}_{\mathrm{w}}\mathrm{~}\mathrm{5}\hspace{0.17em}{\mathrm{deg}}^{-2}\mathrm{~}\mathrm{1}\hspace{0.17em}{\mathrm{pc}}^{-2}\mathit{,}$(1)where A_w = 48.5deg² (i.e. 289.57 pc²) is the projected area of the window W_in that encloses N_w = 252 stars. This window defines our star sample S₂ (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², 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²) 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 $\mathit{\rho }\mathrm{~}\mathrm{5}\mathrm{\pm }\mathrm{3}\hspace{0.17em}\mathrm{stars}\mathit{/}{\mathrm{deg}}^{-2}\mathit{.}$(2)

3.1.2. 1-NNS distribution across Taurus

The 1-NNS median value evaluated within the full region of Taurus (see Table 2) does not differ from the 1-NNS median value evaluated within the three main filaments, and their 1-NNS 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 1-NNS distribution of the whole region has a more populated large-scale tail. Nonetheless, up to the third quantile, the 1-NNS quartiles are similar. As a matter of fact, due to its highly asymmetrical shape, the Taurus 1-NNS distribution is better evaluated using quantiles rather than the mean and standard deviation that are more conveniently suited for symmetrical functions.

3.1.3. 1-NNS distribution: theoretical random vs. Monte Carlo samplings

In this subsection, we compare the theoretical 1-NNS distribution for a spatially random process in an infinite media to the 1-NNS distribution obtained from Monte Carlo samplings in a finite, irregularly-shaped window to establish that the former can be used as a reliable proxy for the latter.

The random theoretical 1-NNS distribution w_R(log r) derived for an infinite media (see Appendix B for the full development) is given by: ${\mathit{w}}_{\mathit{R}}\mathrm{\left(}\log \mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{2}\hspace{0.17em}\ln \mathrm{\left(}\mathrm{10}\mathrm{\right)}\hspace{0.17em}\mathit{\pi }\hspace{0.17em}\mathit{\rho }\hspace{0.17em}{\mathit{r}}^{\mathrm{2}}\hspace{0.17em}\exp \mathrm{(}\mathrm{-}\mathit{\pi \rho }{\mathit{r}}^{\mathrm{2}}\mathrm{)}\mathit{,}$(3)where ρ is the mean intensity of the random process (the average surface density of stars) and r is the 1-NNS. We want to compare this theoretical distribution to the 1-NNS histogram obtained from 10 000 random Monte Carlo samplings within the finite and irregularly-shaped window W_in (see Fig. 3). Their respective moments are similar (see Table 2), with theoretical moments computed from the following expressions: $\begin{array}{ccc}& & \\ & & \\ & & \\ & & \end{array}$where $\overline{\mathit{r}}$, r_0.5, σ², σ², 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: $\begin{array}{ccc}& & \overline{\mathit{r}}\mathrm{=}\frac{\mathrm{1}}{\mathit{N}\mathrm{-}\mathrm{1}}\sum _{\mathrm{1}}^{\mathit{N}}{\mathit{r}}_{\mathit{i}}\mathrm{\pm }{\left[{\mathit{N}}_{\mathit{i}}\mathrm{-}\mathrm{1}\right]}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\hspace{0.17em}{\mathit{\sigma }}_{\mathrm{Tau}}\mathit{,}\\ & & {\mathit{\sigma }}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{(}\mathit{N}\mathrm{-}\mathrm{1}\mathrm{)}}\sum _{\mathrm{1}}^{\mathit{N}}\mathrm{(}{\mathit{r}}_{\mathit{i}}\mathrm{-}\overline{\mathit{r}}{\mathrm{)}}^{\mathrm{2}}\mathrm{\pm }\sqrt{\mathrm{2}}{\left[\mathit{N}\mathrm{-}\mathrm{1}\right]}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{\sigma }}^{\mathrm{2}}\mathit{,}\\ & & \mathit{\sigma }\mathrm{=}\sqrt{\frac{\mathrm{1}}{\mathit{N}\mathrm{-}\mathrm{1}}\sum _{\mathrm{1}}^{\mathit{N}}\mathrm{(}{\mathit{r}}_{\mathit{i}}\mathrm{-}\overline{\mathit{r}}{\mathrm{)}}^{\mathrm{2}}}\mathrm{\pm }{\left[\mathrm{2}\hspace{0.17em}\mathrm{(}\mathit{N}\mathrm{-}\mathrm{1}\mathrm{)}\right]}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathit{\sigma ,}\\ & & \mathit{\gamma }\mathrm{=}\frac{\mathit{N}}{\mathrm{\left(}\mathrm{\right(}\mathit{N}\mathrm{-}\mathrm{1}\mathrm{\left)}\mathrm{\right(}\mathit{N}\mathrm{-}\mathrm{2}\mathrm{\left)}\mathrm{\right)}}\sum _{\mathrm{1}}^{\mathit{N}}\mathrm{(}{\mathit{r}}_{\mathit{i}}\mathrm{-}\overline{\mathit{r}}{\mathrm{)}}^{\mathrm{3}}\mathrm{\pm }{\mathrm{(}\mathit{N}\mathrm{-}{\mathrm{1}}^{\mathrm{)}}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathit{\gamma ,}\\ & & \end{array}$(5)where Λ is the ordered sequence of the 1-NNS 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).

We thus conclude that finite size and edge window effects do not significantly alter the distribution of 1-NNS 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 1-NNS obtained from standard rejection sampling technique using analytical random 1-NNS distribution (Eq. (3)) to b) the 1-NNS associated to random Monte Carlo samplings of 252 stars enclosed within W_in. The two-sample Kolmogorov-Smirnov (KS) and Anderson-Darling (AD) statistical tests indicate that the two 1-NNS distributions cannot be statistically distinguished with a mean p-value 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 1-NNS probability density function (Eq. (3)) may be used as a reliable proxy to describe a 1-NNS 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 1-NNS 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 1-NNS distributions clearly reveals an excess of short spacings in Taurus (see Fig. 3), for which the 1-NNS distribution has smaller central values (median and mean) and a wider dispersion (see Table 2). The 2-sample KS and AD tests, give a p-value 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 1-NNS 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 ~$\sqrt{\mathrm{338}\mathit{/}\mathrm{139}}\mathrm{~}\mathrm{1.5}$.

This also underlines the modest usefulness of the absolute value of the median nearest-neighbor distance to determine clustering property when the stellar census is incomplete. In the following section, we introduce an alternative criterion based on the one-point correlation function to quantitatively assess local departures from pure randomness and to determine characteristic clustering scales.

3.2. One and two-point correlation statistics

In order to identify a criterion that quantifies the departure of a spatial distribution of stars from randomness, we introduce the one-point correlation function based on the 1-NNS distribution. This function is a new tool to assess and quantify the binary and spatial clustering regimes of stars and aims at complementing the two-point correlation statistics.

3.2.1. The two-point and pair correlation functions

The two-point 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: $\mathrm{d}\mathit{P}\mathrm{=}{\mathit{\rho }}_{\mathit{V}}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{\xi }}_{\mathit{V}}\mathrm{(}{r}\mathrm{\left)}\mathrm{\right)}\mathrm{d}{\mathit{V}}_{\mathrm{1}}\mathrm{d}{\mathit{V}}_{\mathrm{2}}$, 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₁ and dV₂ and separated by r (Peebles 1980). Similarly, for an isotropic and stationary spatial process, the standard projected angular two-point 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 = ρ²(1 + ξ(r))dΩ₁dΩ₂, 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Ω₁, dΩ₂ and separated by a projected angular distance r.

Besides the two-point 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: $\mathit{ĝ}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{{\mathit{N}}_{\mathrm{Tau}}^{\mathrm{Pairs}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{{\mathit{N}}_{\mathrm{R}}^{\mathrm{Pairs}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\mathit{,}$(6)where ${\mathit{N}}_{\mathrm{Tau}}^{\mathrm{Pairs}}\mathrm{\left(}\mathit{r}\mathrm{\right)}$ and ${\mathit{N}}_{\mathrm{R}}^{\mathrm{Pairs}}\mathrm{\left(}\mathit{r}\mathrm{\right)}$ are the histogram of separations between stars within Taurus and in a randomly distributed catalogue, respectively.

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 ${\mathit{N}}_{\mathrm{R}}^{\mathrm{Pairs}}\mathrm{\left(}\mathit{r}\mathrm{\right)}$ 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 two-point 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 power-laws: 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 scale-free binary regime may reflect the scale-free fragmentation of collapsing clumps; and (4) the self-similar 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 quasi-constant surface density of pairs.

However, some caution has to be taken when interpreting the two-point correlation function beyond the binary regime, especially when the spatial distribution of stars deviates from isotropy. Indeed the standard form of the two-point 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 two-point 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 two-point 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 one-point correlation function

Similarly to the pair correlation function, we define the one-point 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 second-order characteristics describing the spatial co-location of stars. One estimator of the Ψ function may be defined from the ratio of the 1-NNS Taurus distribution w_Tau(log r) to the 1-NNS random distribution w_R(log r) obtained either by Monte Carlo simulations or from theoretical random probability density function (Eq. (3)): $\mathrm{\Psi }\mathrm{\left(}\log \mathit{r}\mathrm{\right)}\mathrm{=}\frac{{\mathit{w}}_{\mathrm{Tau}}\mathrm{\left(}\log \mathit{r}\mathrm{\right)}}{{\mathit{w}}_{\mathit{R}}\mathrm{\left(}\log \mathit{r}\mathrm{\right)}}\mathrm{·}$(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 1-NNS values. Here we evaluate it bin per bin of spacing. For a random spatial distribution, the one-point correlation function reduces to unity (Ψ(r) = 1).

The one-point correlation Ψ function can be fitted from 1.6 kAU to ~100 kAU, by a power law using Pearson’s chi-squared linear regression taking into account vertical errorbars (see Fig. 4): $\mathrm{\Psi }\mathrm{=}\frac{{\mathit{w}}_{\mathrm{Tau}}}{{\mathit{w}}_{\mathit{R}}}\mathrm{\propto }{\mathit{r}}^{\mathrm{-}\mathit{\alpha }}\mathit{,}$(8)with $\mathit{\alpha }\mathrm{=}\mathrm{1.}{\mathrm{5}}_{\mathrm{0.3}}^{\mathrm{+}\mathrm{0.2}}$ 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² for $\mathit{r}\mathrm{\ll }\mathrm{1}\mathit{/}\sqrt{\mathit{\pi \rho }}$ (see Appendix B), the 1-NNS distribution within Taurus w_Tau follows a shallow, rising power-law function with the following exponent: ${\mathit{D}}_{\mathrm{f}}\mathrm{~}\mathrm{2}\mathrm{-}\mathit{\alpha }\mathrm{=}\mathrm{0.5.}$(9)We note, that this scale-free 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: ${\mathit{r}}_{\mathrm{c}}\mathrm{=}{\mathrm{0.1}}^{\mathrm{◦}}\mathrm{~}\mathrm{0.24}\hspace{0.17em}\mathrm{pc}\mathrm{~}\mathrm{50}\hspace{0.17em}\mathrm{kAU}\mathit{.}$(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. Ultra-wide 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 scale-free 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 power-law 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

In this subsection, we study the spatial distribution of multiple systems and single stars based on the comparison of their 1-NNS statistics.

The S₁ 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₃ sample composed exclusively of Class II and III stars that have been observed at HAR. From that sub-sample, we separately construct the 1-NNS 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 1-NNS 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 1-NNS at close spacings is remarkably high: nearly 40% of the 1-NNS 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 bias-free S₃ 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 p-value of 0.15. We further compare their respective k-NNS distributions, up to k = 8. They too cannot be statistically distinguished (p-value 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 2-NNS 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 (p-value 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 k-nearest 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 ultra-wide pairs

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 ultra-wide 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 arenon-mutual 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 1-NNS distribution for mutual pairs is markedly different from that of non-mutual pairs; the KS test returns a p-value 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 non-mutual pairs is reasonably well fitted by a log-normal function (p-value = 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 1-NNS for non-mutual 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 1-NNS 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 p-value 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 proper-motion pairs, and plausibly even gravitationally bound. From now on, we refer to these systems as ultra-wide pairs (UWPs).

The present-day 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 non-mutual pairs) found in the short-separation end of the 1-NNS distribution. Approximately 60% (resp. 90%) of 1-NNS have separations below 20 kAU (resp. 60 kAU) with a relatively long-uniform distribution. Thus, the continuous power-law 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 scale-free clustering pattern and constitute a major structural spatial imprint in Taurus.

4.2. An Öpik law for the UWPs separation

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, $\mathit{F}\mathrm{\propto }{\mathit{\delta }}_{\mathrm{P}}^{\mathrm{0.14}}$ (see Fig. 7), at least up to an upper cut-off 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 cut-off, 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 log-uniform in this range (F(log δ) = const. or, equivalently ${\mathit{F}}_{\mathit{O}}\mathrm{\left(}\mathit{\delta }\mathrm{\right)}\mathrm{\propto }\frac{\mathrm{1}}{\mathit{\delta }}$), 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₁), a regime where the distinction between multiple systems and single stars strongly suggests that a physical process is at play. Within this sub-population, Class I and Class II are over-represented (and, accordingly, Class III under-represented) relative to the overall Taurus population (see Table 3). A chi-squared test confirms that this difference is statistically significant (p-value 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.

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: ${\mathit{P}}_{\mathit{ij}}\mathrm{=}\mathit{\gamma }\mathrm{·}{\mathit{P}}_{\mathit{i}}\mathrm{·}{\mathit{P}}_{\mathit{j}}\mathrm{=}\mathit{\gamma }\mathrm{·}\mathrm{(}{\mathit{N}}_{\mathit{i}}\mathit{/}{\mathit{N}}^{\mathrm{)}}\mathrm{·}\mathrm{(}{\mathit{N}}_{\mathit{j}}^{\mathrm{\prime }}\mathit{/}\mathrm{(}\mathit{N}\mathrm{-}\mathrm{1}{\mathrm{\left)}}^{\mathrm{\right)}}\mathit{,}$(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 ${\mathit{N}}_{\mathit{j}}^{\mathrm{\prime }}\mathrm{=}{\mathit{N}}_{\mathit{i}}\mathrm{-}\mathrm{1}$ for i = j whereas γ = 2 and ${\mathit{N}}_{\mathit{j}}^{\mathrm{\prime }}\mathrm{=}{\mathit{N}}_{\mathit{j}}$ 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 chi-squared 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.

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 (p-value of 0.007 for a chi-squared 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.

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₃ 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 χ² test, we cannot rule out a random multiplicity pairing (p-value 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 χ² test indicates that a random multiplicity pairing can be rejected with a high significance level (p-value 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 log-uniformly 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 1-NNS distribution of multiple systems outlined in the previous section. Of the 30 multiple stellar systems having their 1-NNS 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 non-mutual 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.

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₃ (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 pre-main sequence closer binaries in Taurus N_n∗ ∝ b^{− n_∗} = 2.6^{− n_∗} (Duchêne & Kraus 2013).

Following Correia et al. (2006), we define the multiplicity fraction per ultra-wide 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: $\mathit{M}{\mathit{F}}_{\mathrm{uw}}\mathrm{=}\frac{\mathit{T}\mathrm{+}\mathit{Q}\mathrm{+}\mathit{...}}{\mathit{B}\mathrm{+}\mathit{T}\mathrm{+}\mathit{Q}\mathrm{+}\mathit{...}}\mathrm{=}\mathrm{65.6}\mathrm{\pm }\mathrm{11.85}\mathrm{\%}\mathit{.}$(12)We similarly define the companion frequency per ultra-wide binary (CF_uw) as: $\mathit{C}{\mathit{F}}_{\mathrm{uw}}\mathrm{=}\frac{\mathrm{2}\mathrm{\times }\mathit{T}\mathrm{+}\mathrm{3}\mathrm{\times }\mathit{Q}\mathrm{+}\mathit{...}}{\mathit{B}\mathrm{+}\mathit{T}\mathrm{+}\mathit{Q}\mathrm{+}\mathit{...}}\mathrm{=}\mathrm{1.77}\mathrm{\pm }\mathrm{0.19.}$(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 ultra-wide binary for UWPs separated by less than 60 kAU (resp. 10 kAU): MF₆₀ = 65.79 ± 13.16% and CF₆₀ = 1.79 ± 0.22 (resp. MF₁₀ = 80.95 ± 19.63% and CF₁₀ = 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 high-order ultra-wide binaries in Taurus found in our sample.

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: $\overline{{\mathit{\delta }}_{\mathrm{P}}\mathrm{\left(}{\mathit{n}}_{\mathrm{\ast }}\mathrm{\right)}}\mathrm{\propto }{\mathit{a}}^{\mathrm{-}{\mathit{n}}_{\mathrm{\ast }}}\mathrm{~}{\mathrm{1.8}}^{\mathrm{-}{\mathit{n}}_{\mathrm{\ast }}}\mathit{,}$(14)where $\overline{{\mathit{\delta }}_{\mathrm{P}}\mathrm{\left(}{\mathit{N}}_{\mathrm{\ast }}\mathrm{\right)}}$ 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 higher-order multiplicity.

4.7. Ultra-wide 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_⊙.

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 low-mass objects that form within the cores are ejected from their cradles to produce an extended population of low-mass stars and brown-dwarfs, 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 brown-dwarf 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.

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.

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 p-value 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 p-value 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:

• non-random 1-NNS distribution of mutual pairs (versuslognormal 1-NNS distribution for non-mutual pairs);

• approximate log-uniform 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);

• non-random 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²−GM₁M₂/r < 0, with μ being the reduced mass (μ = M₁M₂/ (M₁ + M₂), M₁ and M₂ 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: $\begin{array}{ccc}\mathrm{\Delta }\mathit{v}& \mathit{<}& {\left(\mathrm{2}\mathit{G}\mathrm{(}{\mathit{M}}_{\mathrm{1}}\mathrm{+}{\mathit{M}}_{\mathrm{2}}\mathrm{)}\mathit{/}\mathit{r}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\\ & \mathit{<}& \end{array}$(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³ Myr is much longer than the age of the Taurus pre-main 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 solar-type 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 low-but-not-negligible 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 star-forming 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 high-order 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 Class-pairing 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 (Alonso-Floriano et al.2015). A preference has also been identified for high-order multiplicity among these UWPs, as well as in AB Doradus, TW Hydrae and Tucanae-Horologium 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 solar-type stars within high-order 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 late-spectral-type M1–M5 dwarf UWPs (Law et al. 2010). In the latter case, a very high value was found for the multiple fraction per ultra-wide pair, that is, $\mathrm{7}{\mathrm{7}}_{-22}^{\mathrm{+}\mathrm{9}}$% for systems with separations >4 kAU, very close to what we obtain for the high-order-multiple 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 star-forming 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 non-hierarchical 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 long-standing 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 half-life of a binary composed of two solar-type stars separated by 31 kAU is 10 Gyr (Jiang & Tremaine 2010).

We expect that in a low-stellar-density 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 $\begin{array}{ccc}{\mathit{\tau }}_{\mathrm{D}}& \mathrm{=}& \sqrt{\mathrm{3}\mathit{/}\mathrm{14}}\hspace{0.17em}{\mathit{\pi }}^{-1}\hspace{0.17em}{\mathit{M}}_{\mathrm{tot}}^{\mathrm{1}\mathit{/}\mathrm{2}}\hspace{0.17em}\mathrm{(}{\mathit{G}}^{\mathrm{1}\mathit{/}\mathrm{2}}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{c}}\hspace{0.17em}{\mathit{M}}_{\mathrm{c}}\hspace{0.17em}{\mathit{a}}^{\mathrm{3}\mathit{/}\mathrm{2}}{\mathrm{)}}^{-1}\mathrm{;}\\ {\mathit{\tau }}_{\mathrm{D}}& \mathrm{~}& \end{array}$(16)Secondly, we evaluate the destruction timescale in the diffusive regime: $\begin{array}{ccc}{\mathit{\tau }}_{\mathrm{d}}& \mathrm{~}& \mathrm{0.085}\hspace{0.17em}\mathrm{\left(}{\mathit{\sigma }}_{\mathrm{rel}}{\mathit{M}}_{\mathrm{tot}}{\mathit{b}}_{\min }^{\mathrm{2}}\mathrm{\right)}\mathrm{(}\mathit{G}{\mathit{M}}_{\mathrm{c}}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{c}}{\mathit{a}}^{\mathrm{3}}{\mathrm{)}}^{-1}\mathrm{;}\\ {\mathit{\tau }}_{\mathrm{d}}& \mathrm{~}& \end{array}$(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)² is a diffusive parameter with b_min ~ a being the impact parameter cut off.

While the semi-major 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 solar-type 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 star-forming 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.

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 star-formation process; and (2) they form within a single core with a core-to-star 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 α_BM_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: $\mathit{\rho }\mathrm{=}{\mathit{M}}_{\mathrm{tot}}\mathit{/}\mathrm{(}\mathit{\mu }\hspace{0.17em}{\mathit{m}}_{\mathrm{H}}\mathrm{)}\hspace{0.17em}{\mathrm{(}\mathrm{4}\mathit{/}\mathrm{3}\hspace{0.17em}\mathit{\pi }\hspace{0.17em}{\mathit{n}}_{\mathit{C}}^{\mathrm{2}}{\mathit{\delta }}^{\mathrm{3}}\mathrm{)}}^{-1}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{\alpha }}_{\mathrm{B}}\mathrm{)}\mathit{/}{\mathit{\alpha }}_{\mathit{C}}\mathit{,}$(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: $\begin{array}{c}\\ {\mathit{\rho }}_{\mathrm{SS}}\mathrm{~}{\mathrm{10}}^{\mathrm{4}}\hspace{0.17em}\mathrm{c}{\mathrm{m}}^{-3}\hspace{0.17em}\left[\frac{{\mathit{M}}_{\mathrm{tot}}}{\mathrm{0.72}\hspace{0.17em}{\mathit{M}}_{\mathrm{\odot }}}\right]{\mathrm{[}{\frac{\mathrm{0.17}\hspace{0.17em}\mathrm{pc}}{\mathit{\delta }}}^{\mathrm{]}}}^{\mathrm{3}}\frac{\mathrm{[}{\mathit{n}}_{\mathit{C}}\mathit{/}\mathrm{4}{\mathrm{]}}^{\mathrm{2}}}{\mathrm{[}{\mathit{\alpha }}_{\mathit{C}}\mathit{/}\mathrm{0.5}\mathrm{]}}\mathit{,}\end{array}$(19a)$\begin{array}{c}\\ {\mathit{\rho }}_{\mathrm{SM}}\mathrm{~}{\mathrm{10}}^{\mathrm{5}}\hspace{0.17em}\mathrm{c}{\mathrm{m}}^{-3}\hspace{0.17em}\left[\frac{{\mathit{M}}_{\mathrm{tot}}}{\mathrm{1.24}\hspace{0.17em}{\mathit{M}}_{\mathrm{\odot }}}\right]{\mathrm{[}{\frac{\mathrm{0.12}\hspace{0.17em}\mathrm{pc}}{\mathit{\delta }}}^{\mathrm{]}}}^{\mathrm{3}}\frac{\mathrm{[}{\mathit{n}}_{\mathrm{C}}\mathit{/}\mathrm{4}{\mathrm{]}}^{\mathrm{2}}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{[}{\mathit{\alpha }}_{\mathrm{B}}\mathit{/}\mathrm{0.2}\mathrm{]}\mathrm{)}}{\mathrm{[}{\mathit{\alpha }}_{\mathrm{C}}\mathit{/}\mathrm{0.5}\mathrm{]}}\mathit{,}\end{array}$(19b)$\begin{array}{c}\\ {\mathit{\rho }}_{\mathrm{MM}}\mathrm{~}\mathrm{3}\mathrm{\times }{\mathrm{10}}^{\mathrm{7}}\hspace{0.17em}\mathrm{c}{\mathrm{m}}^{-3}\hspace{0.17em}\left[\frac{{\mathit{M}}_{\mathrm{tot}}}{\mathrm{1.47}\hspace{0.17em}{\mathit{M}}_{\mathrm{\odot }}}\right]{\mathrm{[}{\frac{\mathrm{0.017}\hspace{0.17em}\mathrm{pc}}{\mathit{\delta }}}^{\mathrm{]}}}^{\mathrm{3}}\frac{\mathrm{[}{\mathit{n}}_{\mathrm{C}}\mathit{/}\mathrm{4}{\mathrm{]}}^{\mathrm{2}}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{[}{\mathit{\alpha }}_{\mathrm{B}}\mathit{/}\mathrm{0.2}\mathrm{]}}{\mathrm{[}{\mathit{\alpha }}_{\mathrm{C}}\mathit{/}\mathrm{0.5}\mathrm{]}}\mathrm{·}\end{array}$(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: $\mathit{\lambda }\mathrm{=}{\mathit{c}}_{\mathrm{s}}{\left(\frac{\mathit{\pi }}{\mathit{G\rho }}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}$(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: ${\mathit{\lambda }}_{\mathrm{SS}}\mathrm{=}\mathrm{1.2}\left[\frac{{\mathit{\delta }}_{\mathrm{SS}}}{\mathrm{0.17}}\right]\hspace{0.17em}\mathrm{pc}\mathit{,}$(21a)${\mathit{\lambda }}_{\mathrm{SM}}\mathrm{=}\mathrm{0.7}\hspace{0.17em}\left[\frac{{\mathit{\delta }}_{\mathrm{SM}}}{\mathrm{0.12}}\right]\mathrm{pc}\mathit{,}$(21b)${\mathit{\lambda }}_{\mathrm{MM}}\mathrm{=}\mathrm{0.2}\hspace{0.17em}\left[\frac{{\mathit{\delta }}_{\mathrm{MM}}}{\mathrm{0.017}}\right]\mathrm{pc}\mathit{,}$(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 low-density core that becomes supercritical due to the external potential of the Galaxy (Ballesteros-Paredes 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 high-mass (≳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 kAU-wide binary and innermost high-order multiplicity distribution.

The probable primordial nature of UWPs suggests a scenario in which the turbulent fragmentation of a molecular filament forms over-densities (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 low-density 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 star-forming region that gives rise to isolated prestellar cores, with a typical density of 10⁵ 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 star-formation mode that is associated with denser regions such as ρ Ophiuchus, with typical core densities of 10⁷ cm^-3 and sizes of 0.02–0.03 pc (Ward-Thompson 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 star-forming 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 scale-free (ψ = r^-1.53) over three decades and extends the usual binary regime, as defined by the two-point 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²⁺ 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 star-forming 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 1-NNS 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 1-NNS 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 high-order-multiple 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 star-formation 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 “low-mass” SS pairs could point to a different formation scenario or may result from a dynamical process, such as low-mass 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 ultra-wide 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 0501-1 DESC (PI E. Moraux).

References

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: $\mathrm{\Delta }\mathit{\sigma }\mathrm{=}\arccos (\sin {\mathit{\alpha }}_{\mathrm{1}}\sin {\mathit{\alpha }}_{\mathrm{2}}\mathrm{+}\cos {\mathit{\alpha }}_{\mathrm{1}}\cos {\mathit{\alpha }}_{\mathrm{2}}\cos \mathrm{\Delta }\mathit{\delta })\mathit{,}$(A.1)where φ_i,λ_i (i = 1,2) are, respectively, the declination and right ascension of the star i, and Δδ = | δ₂−δ₁ | their absolute difference in declination. Since this above equation is ill-conditioned for small Δσ, it is better to use the following Vincenty formula (Vincenty 1975) applied to a sphere $\begin{array}{ccc}& & \mathrm{\Delta }\mathit{\sigma }\mathrm{=}\\ & & \arctan \left(\frac{\sqrt{{\left(\cos {\mathit{\alpha }}_{\mathrm{2}}\sin \mathrm{\Delta }\mathit{\delta }\right)}^{\mathrm{2}}\mathrm{+}{\left(\cos {\mathit{\alpha }}_{\mathrm{1}}\sin {\mathit{\alpha }}_{\mathrm{2}}\mathrm{-}\sin {\mathit{\phi }}_{\mathrm{1}}\cos {\mathit{\alpha }}_{\mathrm{2}}\cos \mathrm{\Delta }\mathit{\delta }\right)}^{\mathrm{2}}}}{\sin {\mathit{\alpha }}_{\mathrm{1}}\sin {\mathit{\alpha }}_{\mathrm{2}}\mathrm{+}\cos {\mathit{\alpha }}_{\mathrm{1}}\cos {\mathit{\alpha }}_{\mathrm{2}}\cos \mathrm{\Delta }\mathit{\delta }}\right)\mathrm{·}\end{array}$(A.2)

Appendix B: Analytics for k-Nearest Neighborhood statistics

Appendix B.1: Theoretical probability density function (PDF) and cumulative distribution of 1-nearest neighbor separation (1-NNS) for spatial random distribution

In this section, we derive the theoretical distribution of 1-NNS 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: $\mathit{w}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{d}\mathit{r}\mathrm{=}\left(\mathrm{1}\mathrm{-}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{r}}\mathit{w}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{d}\mathit{r}\right)\mathrm{\times }\mathrm{2}\mathit{\pi r}\mathrm{d}\mathit{r\rho ,}$(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: $\frac{\mathrm{d}}{\mathrm{d}\mathit{r}}\left(\frac{\mathit{w}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathrm{2}\mathit{\pi r\rho }}\right)\mathrm{=}\mathrm{-}\mathrm{2}\mathit{\pi r\rho }\frac{\mathit{w}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathrm{2}\mathit{\pi r\rho }}\mathrm{·}$(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: $\mathit{w}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{2}\mathit{\pi \rho r}\exp \mathrm{(}\mathrm{-}\mathit{\pi \rho }{\mathit{r}}^{\mathrm{2}}\mathrm{)}\mathit{.}$(B.3)The PDF integral is unity (i.e., we have ${{}^{\mathrm{\int }}}_{\mathrm{0}}^{\mathrm{+}\mathrm{\infty }}\mathit{w}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{d}\mathit{r}\mathrm{=}{{}^{\mathrm{[}}\mathrm{-}\exp \mathrm{(}\mathrm{-}\mathit{\pi \rho }{\mathit{r}}^{\mathrm{2}}{\mathrm{\left)}}^{\mathrm{\right]}}}_{\mathrm{0}}^{\mathrm{+}\mathrm{\infty }}\mathrm{=}\mathrm{1}$).

Furthermore, from the 1-NNS PDF (Eq. (B.3)), we get the related cumulative distribution function W(r) as: $\begin{array}{ccc}\mathit{W}\mathrm{\left(}\mathit{r}\mathrm{\right)}& \mathrm{=}& \int \begin{array}{c}\mathit{r}\\ \mathrm{0}\end{array}\mathrm{2}\mathit{\pi \rho r}\exp \mathrm{(}\mathrm{-}\mathit{\pi \rho }{\mathit{r}}^{\mathrm{2}}\mathrm{)}\mathrm{d}\mathit{r}\\ & \mathrm{=}& \end{array}$(B.4)that allows us to define the characteristic value r_q associated to the quantile q from: $\mathit{q}\mathrm{=}\mathit{W}\mathrm{\left(}{\mathit{r}}_{\mathit{q}}\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{-}\exp \mathrm{(}\mathrm{-}\mathit{\pi \rho }{\mathit{r}}_{\mathit{q}}^{\mathrm{2}}\mathrm{)}$, to get: ${\mathit{r}}_{\mathit{q}}\mathrm{=}\sqrt{\ln \mathrm{(}\mathrm{1}\mathit{/}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{q}\mathrm{\left)}\mathrm{\right)}\mathit{/}\mathrm{\left(}\mathit{\pi \rho }\mathrm{\right)}}\mathit{.}$(B.5)

Appendix B.2: Moments of first nearest neighbor distribution

Appendix B.2.1: Mean

The theoretical mean $\overline{{\mathit{r}}_{\mathit{t}}}$ of the first nearest neighbor distance in an infinite medium is computed from the PDF as: $\begin{array}{ccc}\overline{{\mathit{r}}_{\mathit{t}}}& \mathrm{=}& E\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{+}\mathrm{\infty }}\mathit{r}\hspace{0.17em}\mathit{w}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{d}\mathit{r}\\ & \mathrm{=}& \end{array}$(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 $\overline{{\mathit{r}}_{\mathit{t}}}$ in a random process of intensity ρ in an infinite medium, i.e., $\overline{{\mathit{r}}_{\mathit{t}}}$ is strictly equal to half of the inverse square root of the density: $\overline{{\mathit{r}}_{\mathit{t}}}\mathrm{\simeq }\frac{\mathrm{1}}{\mathrm{2}\sqrt{\mathit{\rho }}}\mathrm{·}$(B.7)

Appendix B.2.2: Variance

From the PDF w(r), we now compute the theoretical variance ${\mathit{\sigma }}_{\mathit{t}}^{\mathrm{2}}$ (σ_t being the standard deviation) of the first nearest neighbor distance for a random Poisson process, $\begin{array}{ccc}{\mathit{\sigma }}_{\mathit{t}}^{\mathrm{2}}& \mathrm{=}& \int \begin{array}{c}\mathrm{+}\mathrm{\infty }\\ \mathrm{0}\end{array}\mathrm{(}\mathit{r}\mathrm{-}\mathit{r̅}{\mathrm{)}}^{\mathrm{2}}\mathit{w}\mathrm{(}\mathit{r}\mathrm{)}\hspace{0.17em}\mathrm{d}\mathit{r}\\ & \mathrm{=}& \mathrm{2}\hspace{0.17em}\mathit{\pi }\hspace{0.17em}\mathit{\rho }\hspace{0.17em}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{+}\mathrm{\infty }}\mathit{r}\hspace{0.17em}{\left(\mathit{r}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}\sqrt{\mathit{\rho }}}\right)}^{\mathrm{2}}\exp \mathrm{(}\mathrm{-}\mathit{\pi \rho }{\mathit{r}}^{\mathrm{2}}\mathrm{)}\hspace{0.17em}\mathrm{d}\mathit{r}\mathit{.}\end{array}$(B.8)The integration gives: $\begin{array}{ccc}{\mathit{\sigma }}_{\mathit{t}}^{\mathrm{2}}& \mathrm{=}& \frac{-1}{\mathrm{4}\mathit{\pi \rho }}[\exp \mathrm{(}\mathrm{-}\mathit{\pi \rho }{\mathit{r}}^{\mathrm{2}}\mathrm{)}(\mathrm{2}\mathit{\pi }\exp \mathrm{\left(}\mathit{\pi \rho }\hspace{0.17em}{\mathit{r}}^{\mathrm{2}}\mathrm{\right)}\mathrm{erf}\mathrm{\left(}\sqrt{\mathit{\pi \rho }}\hspace{0.17em}\mathit{r}\mathrm{\right)})\\ & & \mathrm{+}\mathit{\pi }\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{2}\sqrt{\mathit{a}}\mathit{r}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{4}{]}_{\mathrm{0}}^{\mathrm{+}\mathrm{\infty }}\mathit{.}\end{array}$(B.9)We then compute series expansions of the integral at x = 0: $\begin{array}{ccc}\underset{\mathit{r}\mathrm{\to }\mathrm{0}}{\lim }{\mathit{\sigma }}_{\mathit{t}}^{\mathrm{2}}& \mathrm{=}& \mathrm{-}\frac{\mathrm{4}\mathrm{+}\mathit{\pi }}{\mathrm{4}\mathit{\pi \rho }}\mathrm{+}\frac{\mathit{\pi }{\mathit{r}}^{\mathrm{2}}}{\mathrm{4}}\mathrm{-}\frac{\mathrm{2}}{\mathrm{3}}\mathit{\pi }\sqrt{\mathit{\rho }}{\mathit{r}}^{\mathrm{3}}\\ & & \mathrm{-}\frac{\mathrm{1}}{\mathrm{8}}\mathrm{(}\mathit{\pi }\mathrm{-}\mathrm{4}\mathrm{)}\mathit{\pi \rho }{\mathit{r}}^{\mathrm{4}}\mathrm{+}\mathit{O}\mathrm{\left(}{\mathit{r}}^{\mathrm{5}}\mathrm{\right)}\mathit{,}\end{array}$(B.10)and then at x = + ∞: $\begin{array}{ccc}\underset{\mathit{r}\mathrm{\to }\mathrm{+}\mathrm{\infty }}{\lim }{\mathit{\sigma }}_{\mathit{t}}^{\mathrm{2}}& \mathrm{=}& \mathrm{-}\mathrm{1}\mathit{/}\mathrm{\left(}\mathrm{2}\mathit{\rho }\mathrm{\right)}(\exp \mathrm{(}\mathrm{-}\mathit{\pi }{\mathit{r}}^{\mathrm{2}}\mathrm{)}[\mathrm{2}\mathit{\rho }{\mathit{r}}^{\mathrm{2}}\mathrm{-}\mathrm{2}\sqrt{\mathit{\rho }}\\ & & \mathrm{+}\mathrm{(}\mathrm{1}\mathit{/}\mathrm{2}\mathrm{+}\mathrm{2}\mathit{/}\mathit{\pi }\mathrm{)}\mathrm{-}\mathrm{(}\sqrt{\mathit{\rho }}\hspace{0.17em}\mathit{\pi r}{\mathrm{)}}^{-1}\\ & & \mathrm{+}\mathrm{(}\mathrm{2}{\mathit{\rho }}^{\mathrm{2}\mathit{/}\mathrm{3}}{\mathit{\pi }}^{\mathrm{2}}\mathit{r}\mathrm{3}{\mathrm{)}}^{-1}\mathrm{+}\mathit{O}\mathrm{(}{\mathit{r}}^{-4}\mathrm{)}]\mathrm{+}\mathrm{1})\mathit{,}\end{array}$(B.11)to finally get: $\begin{array}{ccc}{\mathit{\sigma }}_{\mathit{t}}^{\mathrm{2}}& \mathrm{=}& \mathrm{-}\mathrm{(}\mathrm{2}\mathit{\rho }{\mathrm{)}}^{-1}\mathrm{+}\mathrm{(}\mathrm{4}\mathrm{+}\mathit{\pi }\mathrm{\left)}\mathrm{\right(}\mathrm{4}\mathit{\pi \rho }{\mathrm{)}}^{-1}\\ & \mathrm{=}& \end{array}$(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: $\begin{array}{ccc}{\mathit{\gamma }}_{\mathit{t}}& \mathrm{=}& \mathrm{(}E\mathrm{\left(}{\mathit{r}}^{\mathrm{3}}\mathrm{\right)}\mathrm{-}\mathrm{3}\overline{{\mathit{r}}_{\mathit{t}}}E\mathrm{\left(}{\mathit{r}}^{\mathrm{2}}\mathrm{\right)}\mathrm{+}\mathrm{2}\overline{{\mathit{r}}_{\mathit{t}}}{\mathrm{3}}^{\mathrm{)}}\mathit{/}\mathrm{(}{\mathit{\sigma }}_{\mathit{t}}{\mathrm{)}}^{\mathrm{3}}\\ & \mathrm{=}& \mathrm{(}\mathrm{3}\mathit{/}\mathrm{\left(}\mathrm{4}\mathit{\pi }{\mathit{\rho }}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{\right)}\mathrm{-}\mathrm{3}\overline{{\mathit{r}}_{\mathit{t}}}\mathit{/}\mathrm{\left(}\mathit{\pi \rho }\mathrm{\right)}\mathrm{+}\mathrm{2}\overline{{\mathit{r}}_{\mathit{t}}}{\mathrm{3}}^{\mathrm{)}}\mathit{/}\mathrm{(}{\mathit{\sigma }}_{\mathit{t}}{\mathrm{)}}^{\mathrm{3}}\\ & \mathrm{=}& \mathrm{2}\sqrt{\mathit{\pi }}\mathrm{(}\mathit{\pi }\mathrm{-}\mathrm{3}\mathrm{)}\mathit{/}\mathrm{(}\mathrm{4}\mathrm{-}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}\\ & \mathrm{\approx }& \mathrm{0.63.}\end{array}$(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(1-NNS) in a random process. Let y = Log(r), where r is the 1-NNS variable, and thus we have: $\mathit{w}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{d}\mathit{r}\mathrm{=}\mathit{w}\mathrm{\left(}\mathit{y}\mathrm{\right)}\mathrm{d}\mathit{y,}$(B.14)where w(r) is the PDF of the 1-NNS variable (Eq. (B.3)) and w(y) is the PDF of the logarithm variable. We get the following moments for the random k1-NNS theoretical distributions: $\mathit{w}\mathrm{\left(}\mathit{y}\mathrm{\right)}\mathrm{=}\mathit{w}\mathrm{\left(}\mathit{r}\mathrm{\right)}\left|\frac{\mathit{\partial r}}{\mathit{\partial y}}\right|\mathrm{=}\mathit{w}\mathrm{\left(}\mathit{r}\mathrm{\right)}\frac{\mathrm{d}\mathit{r}}{\mathrm{d}\mathit{y}}\mathrm{·}$(B.15)From this relation, we get the PDF of the first neighbor distance logarithm for a random process as: $\mathit{w}\mathrm{\left(}\mathit{y}\mathrm{\right)}\mathrm{=}\mathrm{2}\hspace{0.17em}\ln \mathrm{\left(}\mathrm{10}\mathrm{\right)}\hspace{0.17em}\mathit{\pi }\hspace{0.17em}\mathit{\rho }\hspace{0.17em}{\mathit{r}}^{\mathrm{2}}\hspace{0.17em}\exp \mathrm{(}\mathrm{-}\mathit{\pi \rho }{\mathit{r}}^{\mathrm{2}}\mathrm{)}\mathit{,}$(B.16)where it reaches its maximun value (mode) at ${\mathit{r}}_{\mathrm{mod}}\mathrm{=}\sqrt{\mathit{\pi \rho }}$ . From Taylor expansion for x ≪ 1, we get: $\mathit{w}\mathrm{\left(}\mathit{y}\mathrm{\right)}\mathrm{~}\mathrm{2}\hspace{0.17em}\ln \mathrm{\left(}\mathrm{10}\mathrm{\right)}\hspace{0.17em}\mathit{\pi }\hspace{0.17em}\mathit{\rho }\hspace{0.17em}{\mathit{r}}^{\mathrm{2}}\hspace{0.17em}\mathrm{+}\mathit{O}\mathrm{\left(}{\mathit{x}}^{\mathrm{4}}\mathrm{\right)}\mathit{,}$(B.17)as we expect in 2D for a random distribution, the first nearest neighbor 1-NNS distribution function increases as a squared power law for $\mathit{r}\mathrm{\ll }\mathrm{1}\mathit{/}\sqrt{\mathit{\pi \rho }}$. From Eq. (B.16) above, we derive its logarithmic expression: $\log \mathit{w}\mathrm{\left(}\mathit{y}\mathrm{\right)}\mathrm{=}\log \mathrm{(}\mathrm{2}\hspace{0.17em}\ln \mathrm{(}\mathrm{10}\mathrm{)}\hspace{0.17em}\mathit{\pi }\hspace{0.17em}\mathit{\rho }\hspace{0.17em}\mathrm{)}\mathrm{+}\mathrm{2}\hspace{0.17em}\mathit{r}\hspace{0.17em}\mathrm{-}\mathrm{(}\mathit{\pi \rho }\mathit{/}\ln \mathrm{10}\hspace{0.17em}\mathrm{)}{\mathit{r}}^{\mathrm{2}}\mathit{.}$(B.18)

In a Poisson point process in a d-dimensional space with intensity ρ, the distance r between a point and its kth neighbor is distributed according to the generalized Gamma distribution: $\begin{array}{c}\\ {\mathrm{𝒫}}_{\mathit{k}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{\approx }\frac{\mathit{d}}{\mathit{r}}\hspace{0.17em}\frac{{\left[\mathit{\rho }{\mathit{V}}_{\mathrm{d}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\right]}^{\mathit{k}}}{\mathrm{\Gamma }\mathrm{\left(}\mathit{k}\mathrm{\right)}}\mathrm{·}\exp \left(\mathrm{-}\mathit{\rho }{\mathit{V}}_{\mathrm{d}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\right)\mathit{.}\\ \end{array}$(B.19)where V_d(r) is the volume ball or radius r in d-space. $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{V}}_{\mathrm{d}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{c}}_{\mathrm{d}}{\mathit{r}}^{\mathit{d}}\end{array}\mathit{,}& & \end{array}$(B.20)and c_d is the unit volume ball in d-space given by: ${\mathit{c}}_{\mathrm{d}}\mathrm{=}{\mathit{\pi }}^{\mathit{d}\mathit{/}\mathrm{2}}\mathit{/}\mathrm{\Gamma }\mathrm{(}\mathit{d}\mathit{/}\mathrm{2}\mathrm{+}\mathrm{1}\mathrm{)}\mathit{,}$(B.21)and Γ is the Gamma function. So, in a planar (2D) distribution, we get: ${\mathrm{𝒫}}_{\mathit{k}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}\mathrm{(}\mathit{\pi \rho }{\mathrm{)}}^{\mathit{k}}}{\mathrm{\Gamma }\mathrm{\left(}\mathit{k}\mathrm{\right)}}\mathrm{·}{\mathit{r}}^{\mathrm{2}\mathit{k}\mathrm{-}\mathrm{1}}\mathrm{·}\exp \mathrm{(}\mathrm{-}\mathit{\rho \pi }{\mathit{r}}^{\mathrm{2}}\mathrm{)}\mathit{.}$(B.22)The cumulative distribution is then: ${\mathrm{𝒲}}_{\mathit{k}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{r}}{\mathrm{𝒫}}_{\mathit{k}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{d}\mathit{r}\mathit{.}$(B.23)So for the second nearest neighbor cumulative distribution we get: ${\mathrm{𝒲}}_{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{r}}{\mathrm{𝒫}}_{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{d}\mathit{r}\mathrm{=}\mathrm{1}\mathrm{-}\mathrm{(}\mathit{\pi \rho }{\mathit{r}}^{\mathrm{2}}\mathrm{+}\mathrm{1}\mathrm{)}\hspace{0.17em}\exp \mathrm{(}\mathrm{-}\mathit{\pi }{\mathit{r}}^{\mathrm{2}}\mathrm{)}\mathit{,}$(B.24)and for the third nearest neighbor cumulative distribution, we have: ${\mathrm{𝒲}}_{\mathrm{3}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{r}}{\mathrm{𝒫}}_{\mathrm{3}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{d}\mathit{r}\mathrm{=}\mathrm{1}\mathit{/}\mathrm{2}\mathrm{\left[}\exp \mathrm{\right(}\mathrm{-}\mathit{\pi }{\mathit{r}}^{\mathrm{2}}\mathrm{)}\mathrm{(}\mathrm{-}\mathit{\pi \rho }{\mathit{r}}^{\mathrm{2}}\mathrm{(}\mathit{\pi \rho }{\mathit{r}}^{\mathrm{2}}\mathrm{+}\mathrm{2}\mathrm{)}\mathrm{-}{\mathrm{2}}^{\mathrm{)}}\mathrm{+}{\mathrm{2}}^{\mathrm{]}}\mathit{.}$(B.25)

Appendix C: Catalogs

This section gathers the catalogs described respectively in Sects. 2 and 4.1 of the paper.

All Tables

All Figures