© ESO 2019

1. Introduction

Young star clusters (YSC) with resolved populations are ideal objects for studying and constraining star and star cluster formation. Most of the theoretical and observational works devoted to the study of YSCs seem to agree that their (pre-)stellar population, still embedded in the pristine molecular cloud, is mostly coeval with an intrinsic age spread of ≈0.5 − 1 Myr and a physical size of ≈1 pc (e.g Zinnecker et al. 1993; Lada & Lada 2003; Pfalzner 2011; Marks & Kroupa 2012; Getman et al. 2018a).

However, many authors report on the presence of an age spread up to several Myr in YSCs (Palla & Stahler 2000; Cargile & James 2010; Cignoni et al. 2010; Reggiani et al. 2011; Bell et al. 2013; Balog et al. 2016; Getman et al. 2018b). Whether the measured age spreads are real and, hence, related to the physics at play during the early stages of the formation of the stars in clusters (Larson 2003; Pflamm-Altenburg & Kroupa 2007, 2009; Klassen et al. 2017) or related to inaccurate evaluation of the impact of observational biases related to differential extinction, stellar variability and complex physical processes like episodic accretion (Baraffe et al. 2009; Da Rio et al. 2010; Jeffries et al. 2011) – is still an open debate.

Recently, Beccari et al. (2017, hereafter B17) reported on the detection of three well separated sequences of pre-main sequence (PMS) stars in an optical color-magnitude diagram (CMD) of an area of ∼1.5 deg radius centered on the Orion Nebula Cluster (ONC). The stars belonging to the three sequences, while being all centered around the Trapezium area, show a different spatial distribution with the apparently oldest (and most populous) population being more spatially spread around the center of the ONC with respect to the youngest one which shows an increasing concentration toward the center. The authors discussed the possibility that differential extinction and/or a population of unresolved binaries might be at the origin of such an observational feature in their CMD. The effect of differential reddening is safely excluded by the fact that the three populations show an identical distribution of visual extinction. While binaries are certainly present as unresolved sources, B17 demonstrated that in order to reproduce the separation in the color distribution of the populations, they require a mass-ratio distribution heavily skewed toward equal mass binaries. While such a population of binaries could not be safely excluded using the data in hand, B17 stressed that it would still imply an exotic population of binaries never observed before in a cluster. In addition, if unresolved binaries or triples were the reason for the discrete sequences in the CMD, then this would contradict that the binary sequence in this assumption is the most concentrated, because binaries are more easily broken up in denser environments (Kroupa 1995). Supported by recent results from infrared spectroscopy from Da Rio et al. (2016), they concluded that a real difference in age is most likely at the origin of the presence of three populations in the CMD. The populations of PMS stars in the ONC were thus formed in bursts of star formation separated in age by ≈1 Myr. We note that the projection effect (i.e., the emergence of apparent unresolved binaries in the inner region due to projection) is larger in the inner region of the ONC and this could conceivably fake a more centrally concentrated younger population. The stellar number density at the center is 4.7 × 10⁴ stars pc⁻³ (McCaughrean & Stauffer 1994) which implies an average projected separation between two stars at the densest place of about 6000 au. The resolution of our survey (OmegaCAM and Gaia DR2) is a couple of 100 au such that the occurrence of an unresolved line-of-sight (non-physical) binary is sufficiently small to be ignored.

Kroupa et al. (2018) proposed a theoretical explanation in support of multiple-bursts of star formation in the ONC, which has been further investigated by Wang et al. (2019). In summary, (1) the formation of a first generation, fueled at the hub of inflowing molecular filaments, is truncated by the feedback of massive stars formed as a part of the first generation; (2) these massive stars are stellar-dynamically ejected and as a consequence the inflow of gas still present in surrounding filaments is restored; (3) the next generation of stars forms. This process can be repeated until either the ionizing stars are not ejected or the filamentary gas reservoir is exhausted or destroyed. The model well predicts the spatial distribution of the three populations as described in B17.

In this paper we use the Gaia DR2 (Gaia Collaboration 2018) in combination with OmegaCAM photometry to safely isolate the ONC stellar population from fore- or background objects and to confirm or discard the detection of multiple sequences in the optical CMD. The impact of unresolved binaries will be studied in detail. If the presence of multiple populations in the ONC is confirmed, the case of the ONC may represent a fruitful testing site for theoretical models of star cluster formation including the potentially important process of its stellar-dynamical modulation. In addition we take advantage of the Gaia spatial resolving power and analyze apparent multiple system candidates using isolated (that is, with negligible fore- or background contamination) ONC star members only. Such a study is possible for the first time thanks to Gaia.

The manuscript is structured as follows: In Sect. 2 we introduce our data sets, then we apply selection criteria to obtain ONC members (Sect. 3). Using this sample we recover the sequences found by B17 and discuss their properties in Sect. 4. The study of the binary stars allowed by Gaia is described in Sect. 5, followed by a discussion and the conclusions in Sect. 6.

2. Data sets

2.1. The Gaia data

We used the python Astroquery package (Ginsburg et al. 2013) to retrieve the Gaia DR2 data (Gaia Collaboration 2018) from the Gaia science archive¹. We downloaded all the objects detected by Gaia that are within a radius of three degrees on the sky from the ONC (RA ≈ 83.75 deg., Dec ≈ −5.48 deg., 1 deg on the sky corresponds to ≈7 pc at the ONC distance of 400 pc) without any additional filtering. The data sample contains 278 444 targets, of which 12% have parallaxes measured with a relative error of 10% or smaller, while 15% do not have a measured parallax. The data sample is dominated by faint stars, with only 43% stars being brighter than 19 mag in the GaiaG-band.

We note that Gaia DR2 photometry in BP and RP filters is not suitable to study the dense and still in-gas embedded regions as the ONC. We follow the Gaia DR2 quality data check in the Appendix A and Sect. 3.3 to study Gaia DR2 photometric capabilities on our data set and conclude that for our analysis it is necessary to use additional photometric data.

2.2. The OmegaCAM catalog

We acquired a new set of deep multi-exposures with OmegaCAM, attached to the 2.6-m VST telescope in Paranal, with the aim of increasing the photometric accuracy toward the faintest region of the optical CMD shown in B17. We surveyed a 3 × 3 deg area around the center of the ONC using the r and i filters, that is, the same bands as adopted in B17. In particular we acquired 10 × 25 s exposures in both filters for each pointing under the proposal 098.C-0850(A), PI Beccari. As in B17, the whole data reduction process (from removal of detector signatures to the extraction and calibration of stellar magnitudes) was carried out at the Cambridge Astronomical Survey Unit (CASU). We downloaded the astrometrically and photometrically calibrated single band catalogs from the VST archive at CASU². We used a large number of stars in common with the AAVSO Photometric All-Sky Survey (APASS) to correct for possible residual offsets in r and i magnitudes between each single exposures. We finally require that a single object should be detected in at least seven out of ten images both in the r and i band in order to be accepted as a real star. Since the error on individual stellar sources is very small (see Fig. 2), the average of the magnitudes measured in each individual frame was adopted as stellar magnitude while the standard error was used as the associated photometric error. The final catalog includes 93846 objects homogeneously sampled in r and i bands down to r ≈ 21 − 22 mag on an 3 × 3 deg area around the cluster’s center.

2.3. The initial catalog

We used the C³ tool from Riccio et al. (2017) in order to identify the stars in common between the OmegaCAM and the Gaia catalog. C³ is a command-line open-source Python script that, among several other options, can cross-match two catalogs based on the sky positions of the sources. We used the C³best matching option and circular selection region with a radius of 50 arcsec that is used for sky partition. We found 84 022 targets in common between the Gaia and the OmegaCAM data. The majority of targets that are in the OmegaCAM catalog but not in the Gaia one are faint (r ⪆ 21) and blue objects (r − i ⪅ 1.2). We note that we also loose a fraction of objects in the most crowded regions, like the center of the ONC (see Appendix A). Hereafter we refer to the Gaia DR2 and OmegaCAM cross-matched catalog as the initial catalog.

3. The ONC members with Gaia

As described in the previous section, our initial catalog includes optical photometry from OmegaCAM and astrometric informations from Gaia DR2 for objects in the range of mag 12 ≲ r ≲ 20 within 3 × 3 deg around the center of the ONC. We now propose a roadmap which allows users to fully exploit the potential offered by the Gaia DR2 astrometric information in efficiently disentangling the stellar populations of a cluster from field objects in the CMD. Such a strategy is adopted here to specifically identify and study the PMS population(s) in the ONC, but we are confident that such criteria can be applied to perform similar studies of stellar populations in close-by young clusters.

3.1. Constraining the reliable ONC members: Parallax selection – C1

We first use the parallaxes provided by Gaia DR2 to reliably identify ONC members. We show as gray histogram in Fig. 1 the distribution of the parallaxes ϖ for all the stars belonging to our initial catalog. The stellar population belonging to the ONC is easily recognized through a well defined peak at ϖ_peak ≈ 2.5 mas corresponding to a distance d ≈ 400 pc, in full agreement with previous estimates of the ONC distance (e.g., Menten et al. 2007). This plot already witnesses the unprecedented opportunity offered by the Gaia DR2 catalog to safely identify and hence study the details of the stellar populations in clusters.

Fig. 1. Left panel: distribution of the parallax, ϖ, of all the objects in the initial catalog (gray histogram) and the bona fide ONC members selected with C1 (black histogram; see Sect. 3). The dashed line show the Gaussian fit to the distribution of parallaxes of the stars in the initial catalog. The vertical line indicates the ±2σ distance from the peak of the distribution at ≈2.52. The top horizontal axis shows corresponding values, d, in parsecs using d [pc] = 1/(ϖ [mas] × 10−3), valid for those objects that have small relative errors (Gaia Collaboration 2018). Right panel: distributions of the RA and Dec proper motions (gray and black histograms, respectively). To select the objects that have proper motions consistent with the bulk of the ONC we fit each distribution with a Gaussian function shown by dashed lines. The gray horizontal bars indicate the 1σ range of the respective fit. The top axis shows the values of the proper motions in km s−1, calculated assuming a distance of 400 pc.

In order to isolate the ONC members, we performed a fit to the parallax distribution (dashed line in Fig. 1) using a Gaussian function in combination with an exponential function. The latter is appropriate to describe the fore- or background sources. The fit shows a peak at

$$\begin{array}{c}\hfill {\varpi }_{\mathrm{ONC}}=2.52\pm 0.15\mathrm{mas},\end{array}$$

where the uncertainty is 1-σ of the Gaussian fit. Note that with σ we here indicate the uncertainty related to the Gaussian fit while elsewhere in the paper we indicate as σ_ϖ the error associated to the value of ϖ as given in the Gaia catalog. To ensure that we only keep the best candidates, we retain as bona fide members all the objects whose parallax falls in the range defined by ±2σ distance from the peak of the distribution (that is 2.22 and 2.82 mas – the vertical lines in Fig. 1). In addition, only the targets having a relative parallax error within 10% are considered as this ensures good quality of the astrometric solution (see also Gaia Collaboration 2018). The final selection criterion C1 is purely based on the parallaxes and is defined as:

$$\begin{array}{c}\hfill 2.22-3{\sigma }_{\varpi }\le \varpi \le 2.82+3{\sigma }_{\varpi }\mathrm{and}{\sigma }_{\varpi }/\varpi \le 0.1,\end{array}$$(1)

where σ_ϖ is the uncertainty on each single parallax given in Gaia DR2.

The black solid line in Fig. 1 shows the distribution of the parallaxes of the 4988 stars selected using the parallax filtering C1 (Eq. (1)). We include the stellar population in a ≈100 pc region from the ONC along the line of sight. Such a region is somehow larger than the physical extension of the ONC star forming region on the sky (e.g., B17 or our Fig. 8). We note, however, that having 10% uncertainty in parallax at the distance of the ONC results in ≈40 pc uncertainty on the distance inferred from the parallax.

By adopting the criterion C1 (Eq. (1)), we uniquely rely on the use of parallaxes, implying that targets with available optical photometry from OmegaCAM but not having precise enough cataloged parallaxes are not considered. Such objects are usually faint or are residing in crowded fields, like the ONC center. In addition multiple or close by objects for which a single star’s astrometric solution was not found would also lack cataloged astrometric parameters. We stress here that 12% of the stars in the initial photometric catalog do not have parallaxes. This fraction is reduced to 4% if we only consider objects that are in the CMD position occupied by the young stellar population of the ONC.

3.2. Constraining the reliable ONC members: Proper motion selection – C2

As previously mentioned, by adopting the C1 (Eq. (1)) criterion we isolate the stars in a region with a depth of ≈100 pc around the ONC. Since this area is much larger than the expected physical size of the ONC (Hillenbrand 1997), it is not unexpected that the use of the C1 (Eq. (1)) criterion alone does not fully filter out fore- or back-ground objects. Since the stars belonging to a given cluster do in general show similar values of proper motions (μ_α, δ), we can use the μ_α, δ available from Gaia DR2 to further constrain the ONC membership of the stars selected with the C1 (Eq. (1)) criterion. We show in the right panel of Fig. 1 the distributions of μ_α and μ_δ (gray and black histograms, respectively) of all the stars in the C1 (Eq. (1)) selected data set.

It is interesting to note that the distributions of μ_α, δ are not Gaussian. The histograms show a prominent main peak together with one or more sub-peaks populating the wings of the distributions. Such sub-peaks might be indicative of a system of sub-clusters or streams located in the vicinity of the main cluster and likely born out of the same molecular cloud. We will investigate this interesting possibility and the detailed nature of such sub-peaks in a forthcoming work, while here we only focus on the main population of the ONC star cluster. To further isolate the bona fide ONC members we hence fit a Gaussian function to the main peaks only. The two Gaussian functions are shown as dashed lines in the right panel of Fig. 1 with peaks at

$$\begin{array}{cc}\hfill {\mu }_{\alpha }^{\mathrm{ONC}}& =1.1\pm 0.6\mathrm{mas}/\mathrm{yr},\hfill \\ \hfill {\mu }_{\delta }^{\mathrm{ONC}}& =0.3\pm 0.5\mathrm{mas}/\mathrm{yr}.\hfill \end{array}$$(2)

These values are consistent with the ones derived by Kuhn et al. (2019) for the inner most 378 selected members, ${\mu }_{\alpha \star }^K=1.51\pm 0.11$ mas yr⁻¹, ${\mu }_{\delta }^K=0.5\pm 0.12$ mas yr⁻¹. The selection criterion nσC2 is defined in the PM space and allows us to include only the stars whose PM in RA and in Dec fall inside a nσ range as follows,

$$\begin{array}{cc}\hfill {\mu }_{\alpha }^{\mathrm{ONC}}-n{\sigma }_{\mathrm{F}{\mu }_{\alpha }}& \le {\mu }_{\alpha }\le +{\mu }_{\alpha }^{\mathrm{ONC}}+n{\sigma }_{\mathrm{F}{\mu }_{\alpha }},\hfill \\ \hfill {\mu }_{\delta }^{\mathrm{ONC}}-n{\sigma }_{\mathrm{F}{\mu }_{\delta }}& \le {\mu }_{\delta }\le +{\mu }_{\delta }^{\mathrm{ONC}}+n{\sigma }_{\mathrm{F}{\mu }_{\delta }},\hfill \end{array}$$(3)

where nσ_F μα and nσ_F μδ are the number of σ values of the Gaussian to the μ_α, δ distributions, respectively. The 1-σ values denoted as σ_F μα and σ_F μδ are shown in the right panel of Fig. 1 as horizontal light and dark gray bars. Their estimated values are written in Eq. (2). While the combined use of the C1 (Eq. (1)) and C2 (Eq. (3)) criteria is quite powerful in removing fore- or backgrounds objects, it certainly also removes genuine ONC members. Still we stress here that such objects show PM properties that deviate from the main ONC populations. Such objects are likely binary systems with separations ≈0.1″ − 0.4″ or stars sitting in regions affected by severe crowding where the performance of Gaia degrades and the astrometric solution can be spurious (Arenou et al. 2018; Lindegren et al. 2018). In the first case, the use of C1 (Eq. (1)), and even more so C2 (Eq. (3)), very nicely serves one of the main goals of this study to verify the presence of multiple populations in the ONC reported by B17 (see Sect. 3.3). In the second case, removing genuine ONC stars that are mostly sitting in the center of the ONC where most of the stars belonging to the youngest population are found, might delete statistically significant kinematic signatures. Hence, the rigorous selection criteria adopted in this work might potentially weaken or even remove any signature of multiple sequences in the CMD. On the other hand, we are confident that, by applying these filters, if any signature of multiple populations is still observed, this would be a solid observational support toward the presence of multiple and sequential bursts of star formation in the ONC.

3.3. Constraining the reliable ONC members: Gaia quality check – C3

There are several factors that can affect the accuracy of the astrometric and photometric parameters measured by the Gaia satellite (Lindegren et al. 2018; Evans et al. 2018; Arenou et al. 2018). Among others, stellar crowding in a high density environment, faintness of the stars and stellar multiplicity are of particular interest for our work. Such effects are at the origin of the fact that for a number of objects observed with OmegaCAM, parallaxes were not available or had too large uncertainty (C1; see above). The Gaia catalog offers several parameters that can be used to verify the quality of the measures associated to given objects (Arenou et al. 2018). In the following we describe a few parameters that are of particular interest for our work:

• Duplicated source criterion: The current version of the Gaia catalog cannot provide any reliable measurements for sources whose separation on the sky is lower than 0.4 − 0.5″ (Lindegren et al. 2018; Arenou et al. 2018). Such objects are flagged in the Gaia DR2 as duplicated sources (duplicated_source == True). Hence, it is immediately understandable that this parameter can be used to identify objects whose measurements are affected by high stellar density (crowding) or unresolved multiplicity. According to Arenou et al. (2018) the average separation among duplicated targets is 0.019″, corresponding to ≈8 au at the distance of the ONC (≈400 pc). In summary, when we apply the criterion C3 we remove all the stars which appear with the duplicated_source == True in the Gaia catalog.

• Astrometric: Arenou et al. (2018) defined the following criterion:

  $$\begin{array}{c}\hfill u\le 1.2max(1,exp(-0.2(G-19.5))),\end{array}$$(4)

  where $u=\sqrt{{\chi }^2/\nu }$, χ²= astrometric_chi2_al and ν= (astrometric_n_good_obs_al −5) suggesting that it can be used to check for spurious astrometric solutions. In the first panel of Fig. A.1, the u parameter is shown as a function of the apparent GaiaG magnitude. The red dashed line shows Eq. (4) for the equal sign. Following the study of Arenou et al. (2018) we will reject from our study all objects whose value of u falls above this line.

The Duplicate source criterion together with the Astrometric criterion are called the C3 criterion. We show in the left panel of Fig. 2 the CMD of all the stars included in our initial catalog. In the right panel of the same figure we instead show the CMD of the 852 bona fide ONC stars selected using the C1 (Eq. (1)) and 1σC2 (Eq. (3)) and C3 filtering criteria described above. The population of PMS stars belonging to the ONC is identified with unprecedented accuracy in our final catalog. A main population of PMS is clearly detected while a second parallel sequence is also clearly visible. We emphasize here that using the criteria C1+C2+C3 we already removed many unresolved binary systems.

Fig. 2. Left panel: CMD for the initial catalog. The red crosses show 3-σ color and magnitudes errors. Right panel: CMD of the parallax and proper motion selected ONC members (criteria C1 (Eq. (1)) and 1σC2 (Eq. (3)) and C3). The color lines are the best fitting isochrones from the Pisa stellar evolutionary model of the three populations.

As previously mentioned, the potential biases introduced by the filtering certainly affect the star counts in the very central part of the cluster because of the presence of stellar crowding and multiple systems. In the next section we will use the catalog of bona fide members to quantitatively investigate the presence of multiple sequences as reported in B17, and their relation to multiplicity. We stress here that the youngest (reddest) populations detected in B17 are also more concentrated toward the center of the cluster. Therefore the biases introduced by the C1 (Eq. (1)), C2 (Eq. (3)) and C3 criteria will certainly affect the statistical significance of the detection of the redder (and possibly younger) sub-populations.

4. Identification of multiple sequences in the CMD

In order to identify the presence of multiple and parallel sequences of PMS objects in the ONC, we applied the same approach as described in B17. First, we calculate the main ridge line of the bluest and most populated PMS population in the range of magnitudes 15.5 <  r <  17 (left panel of Fig. 3). We use the mean ridge line as reference in the (r, r − i) CMD and we then calculate the projected distance in r − i colors of each star from the reference line (see the central panel of Fig. 3). We show in the right panel of Fig. 3 the distribution of the perpendicular distances in color from the main ridge line, Δ(r − i). The black histogram shows the distribution of C1+1σC2+C3 selected data sample. To complement Fig. 3 that is showing projected distribution Δ(r − i) from the main ridge we show, on the request of the anonymous referee, a distribution in vertical distances from the main ridge line Δ(r), see Fig. 4.

Fig. 3. Left panel: CMD of the C1+C2−1σ+C3 selected data sample. The main ridge line is shown as a gray line. Middle panel: CMD rotated along the slope of the main ridge line. Right panel: histograms of the distances in color from the main ridge line. The black solid line corresponds to the data sample after applying all selection criteria (C1, C2 1-σ and C3). The vertical line at Δ(r − i)=0 shows the position of the mean ridge line which serves as a reference line. It also represent the location of the single stars belonging to the main population of the ONC. The two vertical lines Δ(r − i) > 0 indicate the expected positions of the equal mass binary and triple systems, respectively. The red dashed histogram shows the unresolved binary stars from Tobin et al. (2009). The solid-line red histogram shows the distribution of the unresolved binaries that passed the selection criteria (C1, C2 1-σ and C3).

Fig. 4. Complementary plot to Fig. 3 showing the distribution of magnitudes from the main ridge line instead of the perpendicular distance in color and magnitude that is on Fig. 3. We show the position of the main ridge line (at 0.0) as well as unresolved equal-mass binaries and equal-mass triples stars, using vertical lines. The positions of the three peaks relative to the main ridge line and equal-mass binaries or triples is quantitatively comparable to the Fig. 3.

Clearly the distribution in color of the ONC members is described by at least two peaks with a gap in the middle. This feature fully confirms the observation in B17. As already noted by B17, the number of stars belonging to the third population (i.e., the reddest peak) is quite low. We would like to stress here once more that the use of the C1 (Eq. (1)), C2 (Eq. (3)) and C3 filtering preferentially removes stars located in crowded regions. We hence expect that our selection criteria mostly lower the number of stars in the reddest or youngest population that, according to B17 is the one that is most concentrated toward the ONC center.

We used three Gaussian functions (gray line in the figure) to fit the color distribution and in particular the three peaks. We then use these functions to select members belonging to each population. Given the fact that the second and third sub-populations are overall populated by a very low number of stars (54 and 15, respectively) we decided to proceed further with our analysis by dividing the ONC populations into an “old” and “young” population (respectively blue and orange bar in Fig. 3 and afterwards). Hence the stars belonging to the second and third peak are grouped and studied together.

In Fig. 5 we show the distribution on the sky of the two sub-populations (blue and orange points). The two-dimensional Kolmogorov–Smirnov test (Peacock 1983; Fasano & Franceschini 1987; Press et al. 2007) done on the measured distribution suggests that the two populations are not extracted from the same parental population at ≈2.3σ level of significance. We would like to stress here that the data set shown in the figure suffers from severe incompleteness, especially in the central region, because of the Gaia DR2 loss of sensitivity in regions affected by stellar crowding. As extensively discussed by B17, the use of the OmegaCAM data-set alone (i.e., without applying the filtering used in this work) indicate the presence of a prominent concentration of all three populations toward he same center.

Fig. 5. The RA – Dec distribution of the first and second sequence identified in Fig. 3. To compare on-the-sky distributions of the “old” and “young” sequences we use two-dimensional Kolmogorov–Smirnov test (Peacock 1983; Fasano & Franceschini 1987; Press et al. 2007) resulting in p = 0.01 (≈2.3σ).

Gaia offers us the possibility to study, for the first time, the parallax distribution and the proper motions of the two sub-populations. We plot in Fig. 6 the histogram of the parallax (left panel) and proper motions (right panel) distributions of the blue and orange populations. The 2D KS test indicates that the proper motions of the blue and orange sub-populations are not extracted from the same parental population at only ≈1.3σ. For the parallax distributions, which are 1D, the Anderson-Darling statistical test³ indicates that they are not extracted from the same parental distribution at ≈2.3σ (see Alves & Bouy 2012, for suggested line-of-sight complex distribution of stellar populations). Since the members of the ONC were selected such that they have the same proper motion values (see Sect. 2.2) the potential difference in proper motions for the different populations might have been erased. However, it is the use of proper motions that offers the most reliable way to constrain the membership in the ONC with Gaia DR2 as the parallax uncertainties are still too large. The 2.3σ difference is thus suggestive but not conclusive evidence of a real difference in proper motions and needs to be investigated further with future releases of the Gaia catalog.

Fig. 6. Left panel: distribution of the parallax, ϖ, of the identified main-ridge sequence (blue), second +third sequence (orange). The top horizontal axis shows corresponding values, d, in parsec using d [pc] = 1/(ϖ [mas] × 10−3), valid for those objects that have small relative errors (Gaia Collaboration 2018). Right panel: RA and Dec proper motions with the histograms of projected distributions on the sides of the plot.

We used the Anderson-Darling test to explore the possibility that the comparison of the distributions of proper motions are not affected by different distributions in magnitudes among the different stellar populations. The test indicates that the distributions of magnitudes of the bluest population and the two reddest ones combined are extracted from parental populations which cannot be distinguished significantly (p = 0.4 for the magnitudes not to being drawn from the same parental distribution).

To summarize, the currently available data indicate that the ONC objects are gradually more centrally concentrated (proportionally to Δ(r − i)) and suggest potential differences in the parallax distributions between the sequences. Under the assumption that the bluer population is older with respect to the orange one (see B17) the current analysis indicate that the older population of the ONC might be closer to us with respect to the younger one. Both of these findings seems to support the hypothesis that a genuine difference in age is at the origin of the distribution in the CMD of the two sub-populations.

4.1. Multiple sequences in the ONC with binary star analysis

Here we would like to study the presence of unresolved multiple systems among the catalog used to obtain the histograms Fig. 3. Ideally, the main ridge line of the main population (corresponding to Δ(r − i)=0) represents the position of single stars. In reality single stars will have certain scatter due to intrinsic photometric uncertainty and possible age spreads.

Unresolved binary system will have a total magnitude in the range of the primary component up to being brighter by 0.75 mag in the case of equal mass binaries. Thus the upper theoretical limit for Δ(r − i) of unresolved binaries is constructed by projecting the theoretical location on the CMD of the equal mass binary sequence along the main ridge line, as done for the observed stars. The same exercise can be done for the equal mass triple systems. The two gray vertical lines in the right panel of Fig. 3 represent the theoretical maximum values in Δ(r − i) of any unresolved binary and triple systems. But this is not realistically the case because the second and third peaks are too defined, since an age dispersion would soften the peaks if they were due to unresolved binaries and triples. That is, in order to explain the outlying systems through a small age spread, one would need an even more unnatural mass-ratio distribution than already excluded by B17. While it is possible that the objects falling in between the reference line at Δ(r − i)=0 and the one indicating the position of multiples are indeed unresolved binaries, one should invoke a peculiar mass ratio distribution in order to reproduce the observed color distribution and the presence of a gap (see B17 for explicit calculations). The separation of the main ridge line and the second peak remain the same as in B17 and thus we do not repeat the calculations here. We would need to do exactly the same as in B17 with those same data, since the data used here are subject to, among other biases as stated above, the central incompleteness that certainly affects the shape of the distributions in Δ(r − i).

Tobin et al. (2009) published a catalog of 135 spectroscopic binaries (SBs) in the ONC. All the SBs are recovered in our initial catalog. We show as a gray dashed histogram in Fig. 3 the distribution in Δ(r − i) of the unresolved binary stars from Tobin et al. (2009) falling in the magnitude range 15.5 <  r <  17 (64 out of 135). Among the 64 unresolved binaries, 30 populate the Δ(r − i) range coincident with the bluest and oldest population of the ONC, while 34 are in the reddest peak. This sample allow us to test how efficient is our filtering methods in identifying unresolved multiple systems.

In fact we applied to the catalog of 64 unresolved binaries from Tobin et al. (2009) and retrieved in our photometry the C1, C2 and C3 selection criteria. The distribution in the Δ(r − i) plane of the catalog of filtered binaries that survive to our selection criteria is shown as solid red line in the right panel of Fig. 3. Only 16 stars out of the initial 64 (25%) passed our filtering. In particular, nine stars among the 30 objects populating the bluest (30%) peak remains undetected, while only seven among the reddest 34 systems passed our selection criteria (20%). This test demonstrates that our filtering method largely removes the unresolved binary stars, especially those that are in the position of the second peak (that is q ⪆ 0.5). So it indeed shows that our filtering method is able to remove binaries that are present in the second sequence. To confirm this we performed the following experiment: we define as the general population, the sample that has been parallax selected with a 3-sigma cut in proper motion, that is, we keep only the likely members of the regions (noting that such a sample is incomplete). For this general population the above-applied filtering keeps 60% of the stars belonging to the main peak (population one) and 40% of the stars belonging to the second and third peak – proving (a) that the filtered number of targets is smaller for the general population data sample than for the binary star catalog by Tobin et al. (2009), in agreement with the statement that the applied filtering is able to remove binary stars, and (b) that the removed fraction of stars in the range where we expect unresolved binary or triple stars is higher, in line with the previous point (a).

We note that the filtering removes a substantial fraction of multiple systems. In particular, Gaia cannot find astrometric solutions typically for binaries. The exact fraction of binaries filtered out by which of the criteria is however a non-trivial problem and cannot be quantified here as this would depend on detailed modeling of the space craft and the stellar population at the distance of the ONC, see for example Michalik et al. (2014).

4.2. PMS isochrones

Stellar models have been computed using the most recent version of the Pisa stellar evolutionary code (see e.g., Tognelli et al. 2011) described in detail in previous papers (see Randich et al. 2018; Tognelli et al. 2018, and references therein). We simply recall the main aspects relevant for the present analysis. The reference set of models adopts a solar calibrated⁴ mixing length parameter, namely α_ML = 2.0. However, for the comparison we also used a much lower value, namely α_ML = 1.00 (low convection efficiency models) which produces cooler models with inflated radii (see e.g., Tognelli et al. 2012, 2018). Evolutionary tracks of mass ≥1.20 M_⊙ have been computed adopting a core convective overshooting parameter β_ov = 0.15, following the recent calibration of the Pisa models by means of the TZ Fornacis eclipsing binary (Valle et al. 2017).

From the models (in the mass range 0.1−20 M_⊙) we obtained isochrones in the age range [0.1, 100] Myr with a variable spacing in age (down to δt = 0.1 Myr), allowing for a very good age resolution. For the reference set of models we used [Fe/H] = +0.0, corresponding to an initial helium abundance Y = 0.274 and metallicity Z = 0.013 given the Asplund et al. (2009) solar heavy-element mixture.

We have converted the theoretical results into the OmegaCAM r and i absolute AB magnitudes using bolometric corrections we computed using a formalism similar to that described in Girardi et al. (2002), employing the MARCS 2008 synthetic spectra library (Gustafsson et al. 2008), which is available for T_eff ∈ [2500,  8000] K, log g ∈ [ − 0.5,  5.5], and [Fe/H] ∈[−2.0, +0.5] completed with the Castelli & Kurucz (2003) models for T_eff >  8000 K.

Next we compared the evolutionary models with the data in order to identify the best fitting isochrones for the observed PMS populations. This was performed by using the same Bayesian maximum likelihood technique described in Randich et al. (2018). Briefly, the method computes the distance of each star in the observational plane from a given isochrone extracted from our database. The square of the distance is used to define the likelihood of a single star. Then, the total likelihood is defined as the product of all the single star likelihoods.

The parameters to be recovered in this work are the age t and the reddening E(B − V). Thus, along with the grid in age, we also built up a fine grid in E(B − V). To derive the extinction in a given band (i.e., r and i) we adopted the Cardelli et al. (1989) extinction law by re-computing the extinction coefficients A_r/A_V and A_i/A_V for the OmegaCAM filters employed. In the following we have performed the parameters recovery under two assumptions: (1) the spread is caused by unresolved binary stars and (2) the spread is actually due to multiple populations.

In the first case, we derived the age and reddening running the recovery over the whole dataset, allowing for the presence of a binary sequence. We have imposed that the binary sequence has a fixed secondary-to-primary mass ratio q = 0.8 and 1. However, to achieve such a large separation as visible in the CMD, a value of q = 1 is preferred.

In the second case, that is, in case of multiple populations, we adopted another strategy. First, we selected the members belonging to each population (see Sect. 4): the selection identified three populations. We imposed the constraint that all the three populations contained in our sample are affected by the same reddening, in agreement with B17.

Such a condition has been implemented in our recovery code in four steps. First, we obtained the total likelihood – which depends on the age and E(B − V) – for each population. Then, we marginalized over the age to obtain a likelihood which is a function of only E(B − V). This step produces the distribution of the possible E(B − V) for each population. We multiply the E(B − V) distributions obtained at the previous step (for each population), to obtain the E(B − V) distribution for the whole dataset under the assumption that the best E(B − V) is the same. Then, as the best E(B − V) of the whole dataset we chose the maximum of such a distribution. As a final step, we run again the recovery for each population, but using this value of E(B − V) as a prior (Gaussian prior), to derive only the age of each population.

In both cases, to obtain the confidence interval on t and E(B − V), we adopted a Monte Carlo simulation, that is, we perturbed independently each datum within its uncertainty range, to create a given number N_pert = 100 of representations of the ONC population. For each representation we derived the age and reddening using the method discussed above, to obtain a sample of ages and reddening. Then, we defined the best value as the mid of the ordered sample (in t and E(B − V)) and the upper and lower extreme of confidence interval as the 84th and 16th percentile of the distribution, respectively.

The best fitting parameters are,

${\mathrm{age}}_{\mathrm{pop}1}=4.5_{-1.2}^{+1.5}$ Myr,

${\mathrm{age}}_{\mathrm{pop}2}=2.1_{-0.4}^{+0.5}$ Myr and

${\mathrm{age}}_{\mathrm{pop}3}=1.4_{-0.2}^{+0.3}$ Myr,

with E(B − V) = 0.052 ± 0.020 begin consistent with E(B − V) from B17. The corresponding χ² is χ² = 3.6. The best fitting isochrones are shown in Fig. 2 (right panel).

4.3. The apparently old “scattered” objects in the CMD

After applying selection criteria C1 (Eq. (1)), 1-σC2 (Eq. (3)) and C3 we expect to be left with mainly ONC members which are PMS stars. However in right panel of Fig. 2 there is clearly a noticeable group of stars that do not appear to be PMS stars, that is, the objects are bluer than the majority of stars. Manara et al. (2013) studied two such objects in detail using broadband, intermediate resolution VLT X-shooter spectra combined with an accurate method to determine the stellar parameters and the related age of the targets. They show that the two selected stars are actually as young as the bulk of the ONC stars. They conclude that, if only photometry is used as an age estimator then especially for bands sensitive to the presence of accretion, like the here used r and i, several accreting objects may appear scattered from the bulk PMS population on the CMD. That is accretors are most likely the explanation for the presence of the blue-wards scattered points in the CMD. We will investigate these bluer targets in detail in a follow-up work (Beccari et al., in prep.).

5. Apparent binary and triple systems

In the previous section we investigated the impact of unresolved binaries on the detection of multiple populations among the PMS of the ONC. We are now interested in exploiting the potential of the Gaia DR2 catalog in combination with the OmegaCAM photometry to investigate resolved multiple systems. The seeing limited OmegaCAM photometry aims at resolving binaries having separations ⪆0.7″. While the resolution of Gaia DR2 is 0.4″, Ziegler et al. (2018) show that Gaia is only able to reliably recover all binaries down to 1.0″ at magnitude contrast as large as 6 mag.

One method to investigate resolved apparent binary candidates is to construct the so-called Elbow plot (e.g., Larson 1995; Gladwin et al. 1999). The method consists in calculating the distribution of the number density of observed targets (Σ) as a function of on-the-sky separation (θ). As shown in Gladwin et al. (1999) the presence of an elbow in the density distribution described above is caused by the presence of resolved binaries. In short, for each star in our catalog we divide the surrounding area of sky into a set of annuli, by drawing circles of radius θ_i centered on the star, with θ_i = 2θ_i − 1 and i ≥ 1. θ₀ is chosen to be well below the smallest separation in the sample. Next we count the number of companion stars in each annulus and we calculate the Σ(θ)_i as follows,

$$\begin{array}{c}\hfill \mathrm{\Sigma }{(\theta )}_i=\frac{N_i}{\pi N_{\mathrm{tot}}({\theta }_i^2-{\theta }_{i-1}^2)},\end{array}$$(5)

where N_i is the number of stars in the ith bin, N_tot is the total number of stars and θ = (θ_i + θ_i − 1)/2 is the angular separation of measured sources.

We show in Fig. 7 the Elbow plot as calculated using the stars sampled with the Gaia DR2 catalog alone (black dashed line) and Gaia DR2 cross-matched with the OmegaCAM catalog (black solid line). We applied the C1 and 2σC2 criterion to both catalogs. The C1 and 2σC2 criterion reduces significantly the contamination by objects not belonging to the ONC, but on the other hand keeps the data sample as complete as possible. The slope of the elbow is −1.70 ± 0.05 which suggests incomplete data sample toward smaller separations as slope values for complete samples in binary-devoted studies are ≈ − 2 (Gladwin et al. 1999). This is consistent with the recovery completeness found by Ziegler et al. (2018) showing that Gaia DR2 is complete for multiple systems separated by at least 1 arcsec. Prosser et al. (1994) used the Hubble Space Telescope to probe the densest central region of the ONC and Simon (1997) then created the elbow plot. The data are plotted as gray crosses in Fig. 7. The difference between the results of Simon (1997) and our study is most likely caused by the different field of view, as Simon (1997) looked at the inner-most part of the ONC, as well as by the fact that using Gaia, we are able to successfully remove fore- or background contaminators.

Fig. 7. Elbow plot showing the surface number density as a function of separation on the sky. The black dashed line shows the original Gaia data sample, the black solid line shows the cross-matched, that is, initial, catalog after the C1 and 2-σC2 are applied. We fit the elbow part of the diagram obtaining a slope value −1.7 ± 0.05. The gray data points and the DM slope are from Simon (1997) having a smaller field of view by using the HST. The gray line is the Duquennoy & Mayor (1991) distribution of binary stars of G-type primaries that has been reformulated to the number density Σ by Gladwin et al. (1999) (their Eq. (5)). The vertical lines (from left to right) are the resolution limit of Gaia DR2, the resolution limit of OmegaCAM, respectively, the first red line shows (1 × 103 au) as the definition of wide binaries, and the last line shows the start of the Elbow feature (3 × 103 au). The color histograms are constructed for central circular regions (see Fig. 8) to see how the potential of detection of wide binaries changes with the chosen spatial region.

From the Elbow plot we can see that the overabundance of stars, suggesting the presence of multiple systems, starts at ≈3000 au. Interestingly, Scally et al. (1999), based on a proper motion study, suggested that there should be no wide binaries with separations larger than 1000 au. Therefore we decided to investigate the indications of the presence of wide binaries (binaries with apparent separations larger than 1000 au) in our data in more detail.

To investigate the effect of crowding on the Elbow plot we divide the sampled region in four circular regions. Three of the regions are concentric and centered on the ONC center, where the stellar crowding is more severe and hence the completeness expected to be lower. A fourth circle has been centered at 1 deg distance toward the south with respect to the ONC. Such region is located well in the outskirt of the cluster, where the stellar distribution is quite homogeneous and the stellar density very low. We show the location of the circles in the sky in Fig. 8 while in Fig. 7 we show with the same color code the Elbow plot only for objects within each circle. Clearly, the detection of wide binaries in the central region is challenged by most likely the high stellar density. The more we expand toward the external regions, the detection becomes more and more significant (see the blue line).

Fig. 8. On-the-sky distribution of the wide (on-the-sky component separation between 1000 and 3000 au) binaries. Only one component of the binary star candidates is plotted. The red circle shows the region, 15 arcmin = 0.25 deg, investigated by Scally et al. (1999) who found only four candidates. The orange (radius of 0.35 deg) circle shows the region we selected to test how the Elbow plot depends on the investigated spatial scale. The green circle is the (1 deg) large testing region and the blue circle (radius 0.35 deg) is outside where the stellar density is lower.

Next, we use parallax and proper motion cuts to minimize potential fore- or background contaminants (C1 and C2 but assuming 3σ tolerance as here we aim to remain as complete as possible). There are 86 candidates having apparent separation smaller than 3000 au as defined by the Elbow plot, and larger than 1000 au. Such pairs with separations between 1000 and 3000 au are our wide binaries candidates.

As additional criterion for two stars being a binary candidate we implement the relative proper motion criterion introduced by Scally et al. (1999). That is, we check the relative proper motions of all the candidates and require their value to be 0 within the 3σ_μ proper motion uncertainty (the average uncertainty, σ_μ, is larger than the escape velocity from a system with 1000 au separation assuming Solar masses of both components). After this we have 60 wide binary candidates from which ten are in the inner 15 arcmin region used by Scally et al. (1999). Scally et al. (1999) found only four candidates in the very same region. All the binary star candidates fulfilling the separation and the pm criteria are listed in Table A.1, including also binary candidates with separations smaller than 1000 au.

In Fig. 9 we construct the normalized binary fraction distributions defined as the number of binaries over the total number of objects, and compare it with the previous study of Duchêne et al. (1825). We extend the separation range of the Duchêne et al. (1825) study and find excess of wide binaries in comparison to Scally et al. (1999). We compare the central 15 arcmin around the ONC, that is, the same area studied by Scally et al. (1999), with the full region studied here and find comparable binary fraction.

Fig. 9. Binary fraction, defined as the number of binaries over the total number of objects, normalized by the bin width. The black points are data for the ONC from Duchêne et al. (1825). The black histogram shows fb for the spatially not constrained region, that is Gaia DR2 data with C1 and C2 selections. The gray dashed histogram shows the data only for the central part that is the same as in the study by Scally et al. (1999). The last bin showing wide binaries in the range 1000−3000 au represents newly discovered potential binaries that seemed to be absent (Scally et al. 1999) (see Fig. 8). The limit 3000 au is motivated by the beginning of the elbow in Fig. 7.

The main reason why we are able to detect more of these wide apparent binary star candidates than Scally et al. (1999) is likely due to the fact that Gaia DR2 allows us, for the first time, to obtain a clean sample of the ONC region only. The contamination of other sources is otherwise too high. We note that due to the crowding at the very center – in the ONC – the detection of the wide pairs are more likely due to chance projections. Our findings are sharpened up by the proper motion study, nevertheless obtaining the 3D velocity information would provide an additional test.

5.1. Apparent triple systems

In the Gaia data sample, that is Gaia only catalog with C1 parallax cut and 2-σC2 proper motion cut, we identified five apparent triple systems and only one of these has component relative proper motions consistent with zero. To identify triple system candidates we at first identified the closest potential binary companion candidate and quantified the center of light of the candidate binary. We then find its closest neighbor. If the closest neighbor has a projected separation smaller than 3000 au (the maximal separation captured by our data) then the three neighbors are classified as an apparent triple star candidate. The identified candidates are summarized in Table 1.

6. Summary and conclusions

The ONC star forming region has been studied in detail using the data from Gaia DR2 in combination with deep OmegaCAM photometric data. We summarize our main findings bellow: Firstly, we used parallaxes and proper motions to separate the ONC stellar population from background and foreground contamination using Gaia DR2. Cross-matching with the OmegaCAM photometric catalog provides high-quality photometry in r and i filters. This fact turns to be critical in studying the properties of the stellar populations in the ONC as Gaia DR2 photometry does not provide reliable magnitudes in crowded regions (see Appendix A). Using these data we are able to construct the CMD and detect two well separated PMSs with a suggestion for the existence of a third one, despite Gaia DR2 data being incomplete in the central high density regions where the members of the second and third sequences are mostly located.

To explore the effect of unresolved binary stars we applied filtering as thoroughly as possible. The adopted filtering criteria are able to remove basically all unresolved binary stars contributing to the second sequence. In order to probe the capability of our filtering criteria to isolate unresolved binaries, we identified in our catalog 64 unresolved binaries from Tobin et al. (2009). We applied our filtering selection to these stars and demonstrate that we are ab;e to safely identify and filter-out 80% of unresolved binaries.

We confirm the previous result of B17 that a single stellar populations with binaries (or triples) can not alone describe the observed distribution of stars in the CMD. Our study supports the hypothesis that the PMS population of the ONC is better described as three stellar populations with different ages. We use PMS models calculated using the Pisa stellar evolutionary code to find the best fitting isochrones describing the population on the CMD. We estimate an age difference of ∼3 Myr between the youngest (∼1.4 Myr) and the oldest (∼4.5 Myr) stellar population. In line with B17 the youngest sequences appear to be more concentrated toward the ONC center with respect to the oldest population. We perform a series of statistical tests with the aim of identifying kinematic differences between the old and young stellar population (see Sect. 4). We find a mild indication that the individual sequences might have different parallax distributions with the older population slightly closer to us along the line of sight with respect to the youngest one. While this result is in principle in line with the youngest population being still more embedded in the molecular clouds, we note that the uncertainties in the Gaia DR2 data are still too big to allow us to perform such a detailed study with solid statistical significance. More precise data (or additional data – like radial velocities) will be necessary to confirm these results.

Interestingly, even adopting a rigorous filtering strategy via accurate parallaxes and proper motions selection criteria, the stars belonging to the ONC seem to occupy a ≈100 pc region in the line-of-sight direction and ≈20 pc on the sky. Hence this study well extend behind the canonical size of the ONC star cluster, which is only few a pc large (Hillenbrand 1997). While it is well know that the Orion Nebula is in fact a complex star forming region, hosting several episodes of star formations, we stress here that the precision of the Gaia DR2 parallaxes implies, for a 10% relative error, a 40 pc uncertainty at the distance of the ONC.

Secondly, we use our final catalog of bona fide ONC member to study the population of apparent resolved multiple objects in the region. We analyze apparent multiple system candidates using the elbow plot. The position of the elbow suggests that wide binaries that we are able to capture have a separation of up to ≈3000 au (≈7 arcsec). However, this depends on the investigated region and becomes less clear in the dense ONC center. Thus we identify all the targets having a separation smaller than 3000 au and identical proper motions within the measurement uncertainty. This allows us to detect in total 91 targets out of which 60 have a separation between 1000 au and 3000 au. Scally et al. (1999) suggested that wide binaries (semi-major axis larger than 1000 au corresponding to an orbital period longer than 28 000 yr for a system mass of 1.3 M_⊙) are largely absent from the ONC implying that ONC-like clusters cannot be at the origin of the majority of field stars. Here we have found evidence for a large population of wide binaries in the ONC with projected separation up to 3000 au, or an orbital period of ≈144 000 yr for a system mass of ≈1 M_⊙.

Using detailed stellar-dynamical models of the ONC (expanding, collapsing and dynamical equilibrium models) Kroupa et al. (1999) showed that, depending on the dynamical state of the ONC, wide binaries are expected to be present. Moreover, the way binary fraction varies radially can put new constraints on the dynamical state of the ONC. Evidence for the breakup of wide binaries after falling through the cluster has been found in the survey for optical binaries in the ONC by Reipurth et al. (2007). Even more realistic N-body models have been computed by Kroupa (2001) who assumed the ONC to be post-gas expulsion, that is, to have just emerged from the embedded phase. These authors also show that the ONC ought to have a population of wide binaries depending on how concentrated the ONC was prior to the onset of the removal of the residual gas (compare the upper and lower panels in their Fig. 10). Therefore, if there were no wide binaries, then the pre-gas-expulsion ONC would have corresponded to a model with a half-mass radius close to 0.2 pc, while the presence of wide binaries would require it to have been close to 0.45 pc. While firm conclusions cannot yet be reached on the dynamical state of the ONC, in particular also because the dynamical state of very young clusters with masses near to that of the ONC may be complex (Kroupa et al. 2018; Wang et al. 2019), this discussion shows the importance of surveying for binaries in the ONC and for correlating them with their ages, positions and motion vectors relative to the ONC.

Acknowledgments

TJ thanks Hans Zinnecker and Pavel Kroupa for interesting and fruitful discussions and useful comments to the manuscript. In addition TJ acknowledges discussions with Josefa Großschedl and Stefan Meingast, who found important formatting error in the Table A.1. We thank Frédéric Arenou for clarification of the photometric excess parameter in the Gaia DR2 catalog. PGPM, ET and SD acknowledge INFN “iniziativa specifica TAsP”. CFM acknowledges an ESO Fellowship. Based on observations collected at the European Southern Observatory under ESO program 098.C-0850(A).

References

Appendix A: Comparison of OmegaCAM photometry with Gaia DR2 photometry

In the following text we demonstrate that additional photometric data set is needed in dense region as ONC to be able to study the stellar populations. We follow Gaia DR2 release work by Arenou et al. (2018) and show that the photometric information in Gaia Gp and Rp filters gets affected by crowding and is not suitable for precise study. Arenou et al. (2018) introduced the photometric excess criterion,

$$\begin{array}{c}\hfill 1.0+0.015{(G_{\mathrm{BP}}-G_{\mathrm{RP}})}^2<E<1.3+0.06{(G_{\mathrm{BP}}-G_{\mathrm{RP}})}^2,\end{array}$$

where E is defined as the flux ratio in the three different Gaia passbands,

$$\begin{array}{c}\hfill E=\frac{I_{\mathrm{BP}}+I_{\mathrm{RP}}}{I_G}·\end{array}$$(A.1)

The Gaia DR2 presents the photometric set in broad band photometric filters G, G_BP and G_RP. The photometry measurements suffer from significant systematic effects for the faint sources (G ≥ 19), in crowded regions, nearby sources (such as multiple systems) included. The G band is measured using the astrometric CCDs with a window of typical size of 12 × 12 pixels (along-scan × across-scan). However it is 60 × 12 pixels for the BP and RP spectral bands. Thus for a single isolated point-source the G band flux has a comparable value with the sum of the integrated fluxes from the BP and RP bands as expected from the bands’ wavelength coverage. On the other hand, if the source is contaminated by a nearby source or it is extended so it is not fully covered by the G band window, then the G flux is not comparable with the sum of BP and RP fluxes (Gaia Collaboration 2016). The photometric excess criterion described above can hence be used to check the accuracy of the photometry in the Gaia DR2 catalog (Gaia Collaboration 2018).

We show in Fig. A.1 (third panel from the left) the distribution of the excess factor E as a function of G_BP − G_RP color for the entire catalog (gray points) and for the stars selected using the criterion C1 (black points). The red dashed lines indicate the selection limits on this plane as defined using the equation from Arenou et al. (2018) and reported above. According to the E limits defined by Arenou et al. (2018), only the objects falling in between the two lines should be retained. In the last panel on the right of Fig. A.1 we show the distribution of E as a function of RA Clearly, by filtering the data using such a parameter would completely exclude the central region of the ONC star forming region where the distribution of E clearly peaks. This is a direct consequence of the fact that the photometric magnitudes G_bp and G_rp of a given source available in the Gaia DR2 release can be affected by strong contamination by close-by (3.1 × 2.1 arcs², Evans et al. 2018) companions in regions affected by high stellar density. We hence decided to not use the E filtering in our analysis and to use the high quality optical photometry available through the OmegaCAM data.

Fig. A.1. First panel: astrometric u parameter as a function of G magnitude. The red dashed line is the u selection criterion from Arenou et al. (2018). Second panel: shows how this factor depends on RA. The u parameter is affected by crowding in the central parts of the ONC. Third panel: photometric excess, E, parameter as a function of the color in Gaia DR2 filters. The region between the dashed red lines shows the E selection criterion from Arenou et al. (2018). Fourth panel: E as a function of RA. demonstrating that E and therefore Gaia DR2 photometry in bp and rp filters are significantly affected by crowding. In all panels, the black points are the data from the initial catalog.

To show the effect of crowding on the Gaia photometric filters directly we plot the comparison of the CMD for the Gaia photometry and the OmegaCAM one. For that we use the C1 (Eq. (1)) selected data sample. In Fig. A.2 we select the one degree region around the ONC for which we construct the CMD (black points). We do the same for only the inner region (central 0.2 deg) plotted using red points. We clearly see that the OmegaCAM photometry plots mainly PMS objects with no obvious difference between the black and red data point. Contrary to this, Gaia photometry in the bp and rp filters shows significant scatter and deviations from the PMS objects.

Fig. A.2. Top panels: CMD using Gaia photometric filters (left) and using OmegaCAM filters (right). The black points show the one-degree region from the ONC center and the red points the inner 0.2° region. Clearly, the Gaia photometry in bp and rp filters is affected by the central crowding. Bottom panel: RA and Dec plot showing the one-degree region and the 0.2° central selection (red points). In all panels, the black points are the C1 and 2-σC2 selected objects.

All Tables

All Figures

Fig. 1. Left panel: distribution of the parallax, ϖ, of all the objects in the initial catalog (gray histogram) and the bona fide ONC members selected with C1 (black histogram; see Sect. 3). The dashed line show the Gaussian fit to the distribution of parallaxes of the stars in the initial catalog. The vertical line indicates the ±2σ distance from the peak of the distribution at ≈2.52. The top horizontal axis shows corresponding values, d, in parsecs using d [pc] = 1/(ϖ [mas] × 10−3), valid for those objects that have small relative errors (Gaia Collaboration 2018). Right panel: distributions of the RA and Dec proper motions (gray and black histograms, respectively). To select the objects that have proper motions consistent with the bulk of the ONC we fit each distribution with a Gaussian function shown by dashed lines. The gray horizontal bars indicate the 1σ range of the respective fit. The top axis shows the values of the proper motions in km s−1, calculated assuming a distance of 400 pc.

In the text

	Fig. 2. Left panel: CMD for the initial catalog. The red crosses show 3-σ color and magnitudes errors. Right panel: CMD of the parallax and proper motion selected ONC members (criteria C1 (Eq. (1)) and 1σC2 (Eq. (3)) and C3). The color lines are the best fitting isochrones from the Pisa stellar evolutionary model of the three populations.
In the text

Fig. 3. Left panel: CMD of the C1+C2−1σ+C3 selected data sample. The main ridge line is shown as a gray line. Middle panel: CMD rotated along the slope of the main ridge line. Right panel: histograms of the distances in color from the main ridge line. The black solid line corresponds to the data sample after applying all selection criteria (C1, C2 1-σ and C3). The vertical line at Δ(r − i)=0 shows the position of the mean ridge line which serves as a reference line. It also represent the location of the single stars belonging to the main population of the ONC. The two vertical lines Δ(r − i) > 0 indicate the expected positions of the equal mass binary and triple systems, respectively. The red dashed histogram shows the unresolved binary stars from Tobin et al. (2009). The solid-line red histogram shows the distribution of the unresolved binaries that passed the selection criteria (C1, C2 1-σ and C3).

In the text

Fig. 4. Complementary plot to Fig. 3 showing the distribution of magnitudes from the main ridge line instead of the perpendicular distance in color and magnitude that is on Fig. 3. We show the position of the main ridge line (at 0.0) as well as unresolved equal-mass binaries and equal-mass triples stars, using vertical lines. The positions of the three peaks relative to the main ridge line and equal-mass binaries or triples is quantitatively comparable to the Fig. 3.

In the text

	Fig. 5. The RA – Dec distribution of the first and second sequence identified in Fig. 3. To compare on-the-sky distributions of the “old” and “young” sequences we use two-dimensional Kolmogorov–Smirnov test (Peacock 1983; Fasano & Franceschini 1987; Press et al. 2007) resulting in p = 0.01 (≈2.3σ).
In the text

Fig. 6. Left panel: distribution of the parallax, ϖ, of the identified main-ridge sequence (blue), second +third sequence (orange). The top horizontal axis shows corresponding values, d, in parsec using d [pc] = 1/(ϖ [mas] × 10−3), valid for those objects that have small relative errors (Gaia Collaboration 2018). Right panel: RA and Dec proper motions with the histograms of projected distributions on the sides of the plot.

In the text

Fig. 7. Elbow plot showing the surface number density as a function of separation on the sky. The black dashed line shows the original Gaia data sample, the black solid line shows the cross-matched, that is, initial, catalog after the C1 and 2-σC2 are applied. We fit the elbow part of the diagram obtaining a slope value −1.7 ± 0.05. The gray data points and the DM slope are from Simon (1997) having a smaller field of view by using the HST. The gray line is the Duquennoy & Mayor (1991) distribution of binary stars of G-type primaries that has been reformulated to the number density Σ by Gladwin et al. (1999) (their Eq. (5)). The vertical lines (from left to right) are the resolution limit of Gaia DR2, the resolution limit of OmegaCAM, respectively, the first red line shows (1 × 103 au) as the definition of wide binaries, and the last line shows the start of the Elbow feature (3 × 103 au). The color histograms are constructed for central circular regions (see Fig. 8) to see how the potential of detection of wide binaries changes with the chosen spatial region.

In the text

Fig. 8. On-the-sky distribution of the wide (on-the-sky component separation between 1000 and 3000 au) binaries. Only one component of the binary star candidates is plotted. The red circle shows the region, 15 arcmin = 0.25 deg, investigated by Scally et al. (1999) who found only four candidates. The orange (radius of 0.35 deg) circle shows the region we selected to test how the Elbow plot depends on the investigated spatial scale. The green circle is the (1 deg) large testing region and the blue circle (radius 0.35 deg) is outside where the stellar density is lower.

In the text

Fig. 9. Binary fraction, defined as the number of binaries over the total number of objects, normalized by the bin width. The black points are data for the ONC from Duchêne et al. (1825). The black histogram shows fb for the spatially not constrained region, that is Gaia DR2 data with C1 and C2 selections. The gray dashed histogram shows the data only for the central part that is the same as in the study by Scally et al. (1999). The last bin showing wide binaries in the range 1000−3000 au represents newly discovered potential binaries that seemed to be absent (Scally et al. 1999) (see Fig. 8). The limit 3000 au is motivated by the beginning of the elbow in Fig. 7.

In the text

Fig. A.1. First panel: astrometric u parameter as a function of G magnitude. The red dashed line is the u selection criterion from Arenou et al. (2018). Second panel: shows how this factor depends on RA. The u parameter is affected by crowding in the central parts of the ONC. Third panel: photometric excess, E, parameter as a function of the color in Gaia DR2 filters. The region between the dashed red lines shows the E selection criterion from Arenou et al. (2018). Fourth panel: E as a function of RA. demonstrating that E and therefore Gaia DR2 photometry in bp and rp filters are significantly affected by crowding. In all panels, the black points are the data from the initial catalog.

In the text

Fig. A.2. Top panels: CMD using Gaia photometric filters (left) and using OmegaCAM filters (right). The black points show the one-degree region from the ONC center and the red points the inner 0.2° region. Clearly, the Gaia photometry in bp and rp filters is affected by the central crowding. Bottom panel: RA and Dec plot showing the one-degree region and the 0.2° central selection (red points). In all panels, the black points are the C1 and 2-σC2 selected objects.

In the text