© ESO, 2017

1. Introduction

OB stars are not distributed randomly in the sky, but cluster in loose, unbound groups, which are usually referred to as OB associations (Blaauw 1964). In the solar vicinity, OB associations are located near star-forming regions (Bally 2008), hence they are prime sites for large scale studies of star formation processes and of the effects of early-type stars on the interstellar medium.

At the end of the last century, the data of the Hipparcos satellite (ESA 1997) allowed the characterization of the stellar content and kinematic properties of nearby OB associations, deeply changing our knowledge and understanding of the solar vicinity and the entire Gould’s Belt (de Zeeuw et al. 1999). The canonical methods used for OB association member identification rely on the fact that stars belonging to the same OB association share the same mean velocity (plus a small random velocity dispersion). The common space velocity is perceived as a motion of the members toward a convergent point in the sky (for more details see, e.g., de Bruijne 1999; Hoogerwerf & Aguilar 1999). Unfortunately, the motion of the Orion OB association is directed primarily radially away from the Sun. For this reason, the methods of membership determination with the Hipparcos proper motions did not perform well in Orion.

The Orion star-forming region is the nearest (d ~ 400 pc) giant molecular cloud complex and it is a site of active star formation, including high mass stars. All stages of star formation can be found here, from deeply embedded protoclusters to fully exposed OB associations (e.g., Brown et al. 1994; Bally 2008; Briceno 2008; Muench et al. 2008; Da Rio et al. 2014; Getman et al. 2014). The different modes of star formation occurring here (isolated, distributed, and clustered) allow us to study the effect of the environment on star formation processes in great detail. Moreover, the Orion region is an excellent nearby example of the effects that young, massive stars have on the surrounding interstellar medium. The Orion-Eridanus superbubble is an expanding structure, probably driven by the combined effects of ionizing UV radiation, stellar winds, and supernova explosions from the OB association (Ochsendorf et al. 2015; Schlafly et al. 2015).

The Orion OB association consists of several groups of varying ages that are partially superimposed along our line of sight (Bally 2008) and extend over an area of ~ 30° × 25° (corresponding to roughly 200 pc × 170 pc). Blaauw (1964) divided the Orion OB association into four subgroups. Orion OB1a is located northwest of the Belt stars and has an age of about 8 to 12 Myr (Brown et al. 1994). Orion OB1b contains the Belt stars and has an age estimate ranging from 1.7 to 8 Myr (Brown et al. 1994; Bally 2008). Orion OB1c (Bally 2008, estimated age from 2 to 6 Myr) includes the Sword stars and is located directly in front of the Orion Nebula, M43, and NGC 1977. Hence, it is very hard to separate the stellar populations of OB1c and OB1d, the latter corresponding to the Orion Nebula Cluster (ONC; see, e.g., Da Rio et al. 2014). It is not clear whether the entire region is a single continuous star-forming event, where Ori OB1c is the more evolved stellar population emerging from the cloud in which group 1d still resides, or whether 1c and 1d represent two different star formation events (see, e.g., Muench et al. 2008). In subsequent studies, many more subgroups have been identified, such as 25 Ori (Briceño et al. 2007), σ Ori (Walter et al. 2008), and λ Ori (Mathieu 2008). Even though the σ Ori and 25 Ori subgroups are located in the direction of the Orion OB1a and OB1b subgroups, the former subgroups have different kinematic properties with respect to the traditional association members (Briceño et al. 2007; Jeffries et al. 2006); the λ Ori group (Mathieu 2008) formation could have been triggered by the expansion of the bubble created by Orion OB1a. Its age and distance from the center of OB1a are also similar to those of OB1c. More recently, Alves & Bouy (2012) and Bouy et al. (2014) reported the discovery of a young population of stars in the foreground of the ONC, which was however questioned by Da Rio et al. (2016), Fang et al. (2017), and Kounkel et al. (2017a). Finally, Kubiak et al. (2017) identified a rich and young population surrounding ϵ Ori.

In this study, we use the first Gaia data release (Gaia Collaboration 2016a,b), hereafter Gaia DR1, to explore the three-dimensional arrangement and age ordering of the many stellar groups between the Sun and the Orion molecular clouds; the overall goal is to construct a new classification and characterization of the young, non-embedded stellar population in the region. Our approach is based on the parallaxes provided for stars brighter than G ~ 12 mag in the Tycho-Gaia Astrometric Solution (TGAS; Michalik et al. 2015; Lindegren et al. 2016) subset of the Gaia DR1 catalog, and on the combination of Gaia DR1 and 2MASS photometry. These data are briefly described in Sect. 2. We find evidence for the presence of a young (age <20 Myr) population, loosely clustered around the following known groups: 25 Ori, ϵ Ori, and σ Ori, and NGC 1980 and the ONC. We derive distances to these subgroups and (relative) ages in Sect. 3. In Sect. 4 we use the Pan-STARRS1 photometric catalog (Chambers et al. 2016) to confirm our age ranking. Our results, which we discuss in Sect. 5 and summarize in Sect. 6, are the first step in using Gaia data to unveil the complex star formation history of Orion and give a general overview of the episodes and the duration of the star formation processes in the entire region.

2. Data

The analysis presented in this study is based on the content of Gaia DR1 (Gaia Collaboration 2016a; van Leeuwen et al. 2017), complemented with the photometric data from the 2MASS catalog (Skrutskie et al. 2006) and the Pan-STARRS1 photometric catalog (Chambers et al. 2016). Figure 1 shows the field selected for this study (1)We chose this field by slightly enlarging the region considered in de Zeeuw et al. (1999). We performed the cross-match using the Gaia archive (Marrese et al. 2017). The query is reported in Appendix B. In the cross-match with 2MASS, we included only the sources with photometry flag “ph_qual = AAA” and we requested the angular distance of the cross-matched sources to be < 1′′. We decided to exclude from our analysis the sources that are either young stars inside the cloud or background galaxies. We performed this filtering with a (J−K) versus (H−K_s) color-color diagram, where extincted sources are easily identified along the reddening band. Following Alves & Bouy (2012), we required that (2)The first condition is taken as the border between non-extincted and extincted sources. The second is meant to reject faint sources to make the selection more robust against photometric errors. The third condition excludes sources with dubious infrared colors (either bluer or redder than main sequence stars). The total number of Gaia sources in the field is N = 9 926 756. The number of stars resulting from the cross-match with 2MASS is N = 5 059 068, which further decreases to only N = 1 450 911 after applying the photometric selection. Figure 2 shows a schematic representation of the field. The stellar groups relevant for this study are indicated as black empty circles and red stars. The coordinates of the stars and clusters shown are reported in Table 1. H_α emission (Finkbeiner 2003) is shown with blue contours, while dust structures (Planck Collaboration XI 2014) are plotted in black.

3. Orion in Gaia DR1

In this section we identify and characterize the stellar population toward Orion. At first, we focus on the TGAS subsample and, after making a preliminary selection based on proper motions, we study the source distribution in parallax intervals. We notice the presence of an interesting concentration of sources toward the center of the field, peaking roughly at parallax ϖ = 2.65 mas (Sect. 3.1). The sources belonging to this concentration also create a sequence in the color magnitude diagrams made combining Gaia DR1 and 2MASS photometry (Sect. 3.2). These findings prompt us to look at the entire Gaia DR1. In the same color magnitude diagrams, we notice the presence of a young sequence, well visible between G = 14 mag and G = 18 mag, which we interpret as the faint counterpart of the TGAS sequence. We make a preliminary selection of the sources belonging to the sequence, and we study their distribution in the sky, finding that they correspond to the TGAS concentrations (Sect. 3.3). We refine our selection and finally we determine the ages of the groups we identify (Sect. 3.4).

3.1. Distances: the Tycho-Gaia subsample

Parallaxes and proper motions are available only for a subsample of Gaia DR1, namely the TGAS (Michalik et al. 2015; Lindegren et al. 2016). We considered all the TGAS sources in the field. Since the motion of Orion OB1 is mostly directed radially away from the Sun, the observed proper motions are small. For this reason, a rough selection of the TGAS sources can be made requiring (3)where μ_{α ∗} and μ_δ are the proper motions in right ascension and declination. The selection above follows roughly de Zeeuw et al. (1999). Figure 3 shows the distribution in the sky of the sources selected with Eq. (3)as a function of their parallax ϖ, from small (ϖ = 0 mas) to large parallaxes up until ϖ = 5 mas (therefore until d = 200 pc). The outline of the Orion A and B clouds and of the λ Ori dust ring is visible (compare with Fig. 1) in the first panel, which shows sources further away than d = 500 pc. This gives us confidence that the sorting of sources in distance (through parallax) is correct. The second panel in Fig. 3 shows stars with parallax 2 <ϖ< 3.5 mas, which corresponds to a distance 285 <d< 500 pc. Some source overdensities toward the center of the field, (l,b) ~ (205°,−18°), are clearly visible, and they are not due to projection effects but are indicative of real clustering in three-dimensional space. We studied the distribution in the sky of the sources with parallaxes 2 <ϖ< 3.5 mas using a kernel density estimation (KDE). The KDE is a non-parametric way to estimate the probability density function of the distribution of the sources in the sky without any assumption on their distribution. Furthermore, it smooths the contribution of each data point over a local neighborhood and it should therefore deliver a more robust estimate of the structure of the data and its density function. We used a multivariate normal kernel with isotropic bandwidth =0.4°. This value was chosen empirically as a good compromise between over- and undersmoothing physical density enhancements among random density fluctuations. To avoid projection distortions, we used a metric where the distance between two points on a curved surface is determined by the haversine formula. The details of the procedure are described in Appendix C.

To assess the significance of the density enhancements we assume that the field stars are distributed uniformly in longitude, while the source density varies in latitude. We thus averaged the source density over longitude along fixed latitude bins and we estimated the variance in source density using the same binning. The significance of the density enhancements is (4)where D(l,b) is the density estimate obtained with the KDE, ⟨D(b)⟩ is the average density as a function of latitude, and Var (D(b)) is the variance per latitude. Figure 4 shows the source probability density function and the black contours represent the S = 3 levels. Figure 5 shows the KDE of the parallax distribution of all the sources with 2 < ϖ < 3.5 mas and of those within the S = 3 contour levels (solid orange and blue dashed line, respectively). We used a Gaussian kernel with bandwidth =0.1 mas, which is comparable to the average parallax error (~ 0.3 mas). The distribution of the sources within the S = 3 contour levels peaks at ϖ ~ 2.65 mas. This supports the notion that the stars within the density enhancements are concentrated in space. To confirm the significance of the difference between the parallax distribution of the two samples, we performed N = 1000 realizations of the parallax density distribution (of both samples) by randomly sampling the single stellar parallaxes and then we computed the 5th and the 95th percentiles, which are shown as fine lines in 5. Finally, we noticed that the spread in the parallax distribution (~0.5 mas) is larger than the typical parallax error, therefore we can hypothesize that it is due to an actual distance spread of ~150 pc and not only to the dispersion induced by the errors.

Figure 6 shows the median parallax over bins of 1° × 1° for the sources within the S = 3 levels. The stars associated with 25 Ori have slightly larger parallaxes than those in the direction toward the ONC, which implies smaller distances from the Sun. We computed the median parallaxes in 2° × 2° boxes centered in 25 Ori, ϵ Ori, and the ONC. We obtained

• 25 Ori: (d ~ 355 pc);

• ϵ Ori: (d ~ 362 pc);

• ONC: (d ~ 413),

where the quoted errors correspond to the 16th and 84th percentiles.

These values are consistent with the photometric distances determined by Brown et al. (1994): i.e., 380 ± 90 pc for Ori1a; 360 ± 70 pc for Ori OB1b; and 400 ± 90pc for OB1c. Using the Hipparcos parallaxes de Zeeuw et al. (1999) reported the mean distances to be 336 ± 16 pc for Ori OB1a, 473 ± 33 pc for Ori OB1b, and 506 ± 37pc for Ori OB1c. Distances to the ONC have been determined by, among others, Stassun et al. (2004), Hirota et al. (2007), Jeffries (2007), Menten et al. (2007), Sandstrom et al. (2007), Kim et al. (2008) and Kraus et al. (2009). These distance estimates range from to 437 ± 19 pc. The latest distance estimate was obtained by Kounkel et al. (2017b), who found a distance of 388 ± 5 pc using radio VLBA observations of young stellar objects. Thus, the TGAS distances are in agreement with the estimates above.

3.2. Color magnitude diagrams

We combined Gaia and 2MASS photometry to make color-magnitude diagrams of the sources within the S = 3 levels defined in Fig. 4. These sources define a sequence at the bright end of the color-magnitude diagram (black big dots in Fig. 7, left). The spread of the sequence does not significantly change using apparent or absolute magnitudes. This prompted us to look further at the entire field, using the entire Gaia DR1 catalog to find evidence of the faint counterpart of the concentration reported in Sect. 3.1. Figure 7 (left) shows a G versus G−J color magnitude diagram of the central region of the field, with coordinates Figure 7 (right) shows the same color magnitude diagram after unsharp masking. A dense, red sequence is visible between G = 14 mag and G = 18 mag. This kind of sequence (also reported, for example, by Alves & Bouy 2012) indicates the presence of a population of young stars. Indeed, the locus of the sequence is situated above the main sequence at the distance of Orion. Several basic characteristics can be inferred from the diagram:

• 1.

  The density of the sequence suggests that the population is rich.

• 2.

  The sequence appears not to be significantly affected by reddening, indicating that the sources are in front of or at the edges of the clouds.

• 3.

  The dispersion of the sequence is ~ 0.5 mag. This can be due to multiple reasons, such as the presence of unresolved binaries, presence of groups of varying ages or distances, or of field contaminants.

Since our field is large, the number of contaminants is high. Therefore, we decided to eliminate the bulk of the field stars by requiring the following conditions to hold (orange line in Fig. 7 left): (5)

3.3. Source distribution

We chose to study the distribution in the sky of the sources selected with Eq. (5)repeating the procedure explained in Sect. 3.1. We analyzed the source density using a multivariate normal kernel with isotropic bandwidth =0.3° and haversine metric. Figure 8 shows the normalized probability density function of the source distribution on the sky. The dashed contours represent the S = 3 levels of the TGAS density map. The density enhancements toward the center of the field are in the same direction as the groups shown in Fig. 2 and reported in Table 1. The density peak in (l,b) ~ (206°,−12.5°) is associated with the old open cluster NGC 2112 (age ~1.8 Gyr and distance ~ 940 pc; see, e.g., Carraro et al. 2008, and references therein).

Figure 9 shows D(l,b)−⟨D(b)⟩ (same notation as in Sect. 3.2) and the contours represent the S = 1 (gray) and S = 2 (black) significance levels. A certain degree of contamination is present, however, the groups are clearly separated from the field stars. Aside from the structures already highlighted in the TGAS map of Fig. 4, the following features are visible in the KDE of Fig. 9:

• The density enhancements toward λ Ori include not only the central cluster (Collinder 69; ~ (195°,−12°), but also some structures probably related to Barnard 30 (~ 192°,−11.5°) and LDN 1588 (~ 194.5°,−15.8°). Some small overdensities are located on the Hα bubble to the left of LDN 1588 and they do not correspond to any previously known group.

• The shape of 25 Ori is elongated and presents a northern and a southern extension, which are also present in the TGAS KDE of Fig. 4.

• South of ϵ Ori, a significant overdensity is present, possibly related to the Orion X group discovered by Bouy & Alves (2015).

• Around the center of the Orion dust ring (~ 214°,−13°) discovered by Schlafly et al. (2015) a number of density enhancements are present. These overdensities are also visible in the TGAS map of Fig. 4, but here they are more evident.

For the following analysis steps, we selected all the sources related to the most significant density enhancements, that is, those within the S = 2 contour levels shown in Fig. 9.

3.4. Age estimates

To determine the age(s) of the population(s) we identified, we performed a Bayesian isochrone fit using a method similar to that described in Jørgensen & Lindegren (2005) and, more recently, in Valls-Gabaud (2014). These authors used Bayesian theory to derive stellar ages based on a comparison of observed data with theoretical isochrones. Age (t) is one free parameter of the problem, but not the only free parameter: the initial stellar mass (m) and chemical composition (Z) are also considered as model parameters. We simplified the problem assuming a fixed value for Z. Using the same notation as Jørgensen & Lindegren (2005), the posterior probability f(t,m) for the age and mass is given by (6)where f₀(t,m) is the prior probability density and L the likelihood function. Integrating with respect to m gives the posterior probability function of the age of the star, f(t). We assume independent Gaussian errors on all the observed quantities with standard errors σ_i. The likelihood function is then written as with where n is the number of observed quantities, and q^obs and q(t,m) are the vectors of observed and modeled quantities. Following Jørgensen & Lindegren (2005), we write the prior as where ψ(t) is the prior on the star formation history and ξ(m) is the prior on the initial mass function. We assume a flat prior on the star formation history and a power law for the initial mass function (IMF) with a = 2.7. We chose a power law following Jørgensen & Lindegren (2005). We also tested other IMFs and find that the final results are not strongly dependent on the chosen IMF. We adopted the maximum of f(t) as our best estimate of the stellar age. We computed the confidence interval following the procedure explained in detail in Jørgensen & Lindegren (2005). It might happen that the maximum of f(t) coincides exactly with one of the extreme ages considered. In this case, only an upper or a lower bound to the age can be set and we call our age estimate ill defined. On the other case, if the maximum of f(t) falls within the age range considered, we call our age estimate well defined.

To perform the fit we compared the observed G magnitude and G−J color to those predicted by the PARSEC (PAdova and TRieste Stellar Evolution Code; Bressan et al. 2012; Chen et al. 2014; Tang et al. 2014) library of stellar evolutionary tracks. We used isochronal tracks from log (age/yr) = 6.0 (1 Myr) to log (age/yr) = 8.5 (200 Myr) with a step of log (age/yr) = 0.01. We chose the range above because we are mainly interested in young (age <20 Myr) sources. As mentioned above, we fixed the metallicity to Z = 0.02, following Brown et al. (1994). The isochronal tracks have an extinction correction of A_V = 0.25 mag. The correction was derived computing the average extinction toward the stars in Brown et al. (1994). We decided to fix the extinction to a single value mainly to keep the problem simple. Besides, we excluded most of the extincted sources when we applied the criteria of Eq. (2).

We applied the fitting procedure to all the stars resulting from the selection procedure in Sect. 3.3, fixing the parallax to the mean value derived in the Sect. 3.1, that is, ϖ = 2.65 mas. This choice is motivated primarily by the fact that with the current data quality is not possible to disentangle the spatial structure of the region precisely. We experimented with more sophisticated choices for the parallax values, however, even if these options lead to different single age estimates, they do not change the general conclusions of the analysis. In particular, the age ranking of the groups does not change.

Figure10 shows the color magnitude diagram of the sources with estimated age younger than 20 Myr. The gray crosses are the sources whose age is ill-defined, the black dots represent the sources with well-defined ages. Noteworthy, the sources with ill-defined age consist mainly of galactic contaminants, which we could then remove from our sample.

Figure 11 shows the density (obtained with a Gaussian kernel, with bandwidth =0.3°) of the source sky distribution as a function of their age, t. The densities are normalized to their individual maximum, so that their color scale is the same. The coordinates of the density enhancements change with time. This means that the groups we identified have different relative ages as follows:

• σOri. The peak associated wth σ Ori ((l,b) = (207,−17.5) deg) is in the first panel (1 <t< 3 Myr), and some residuals are present also in the second (3 <t< 5 Myr) and fourth (7 <t< 9 Myr) panels. Hernández et al. (2007), Sherry et al. (2008), and Zapatero Osorio et al. (2002) all estimate an age of 2−4 Myr, which is compatible with what we find. Instead, Bell et al. (2013) puts the cluster at 6 Myr.

• 25 Ori. The 25 Ori group ((l,b) = (20.1,−18.3) deg) appears in the third panel (5 <t< 7 Myr), peaks in the sixth panel (9 <t< 11 Myr), and then fades away. Briceño et al. (2007) found that the age of 25 Ori is ~ 7−10 Myr. Our age estimate is slightly older, but still fits the picture of 25 Ori being the oldest group in the region.

• Belt population. The population toward ϵ Ori ((l,b) ~ (205.2,−17.2) deg) becomes prominent for t> 9 Myr. Here, Kubiak et al. (2017) estimated the age to be older than ~ 5 Myr without any other constraint.

• ONC, NGC 1980, NGC 1981, and NGC 1977. The overdensities associated with NGC 1980, NGC 1981, NGC 1977, and the ONC ( centered in (l,b) ~ (209,−19.5) deg) are very prominent until the eighth panel of Fig. 11. In this last case it is difficult to disentangle exactly which group is younger, especially because the underlying data point distribution is smoothed by the kernel. The density enhancement in the first panel (1 <t< 3 Myr) is most likely related to the ONC and L1641 (Reggiani et al. 2011; Da Rio et al. 2014, 2016). The density enhancement associated with NGC 1977 peaks in the same age ranges (7 <t< 9 Myr) as that associated with NGC 1980, which however remains visible until later ages (15 <t< 20 Myr) and fades away only for t> 20 Myr. Finally, the density enhancement associated with NGC 1981 does not clearly stand out in any panel, excluding perhaps the density enhancements with age 11 <t< 13 Myr and 13 <t< 15 Myr. An interesting feature of the maps is the fact that the shape and position of the density enhancements related to NGC 1980 change with time. In particular, for early ages only one peak is present, while from ~7 Myr two peaks are visible. This is a further confirmation that the density enhancements in the first three age panels include L1641 and the ONC, which are indeed younger than the other groups. Bouy et al. (2014) derived an age ~ 5−10 Myr for NGC 1980 and NGC 1981.

The last panel shows the stars with estimated ages >20 Myr. The source distribution is uniform. These are field stars with estimated ages ranging from 20 to 200 Myr.

Our fitting procedure did not take into account the presence of unresolved binaries among our data. Since the sample includes pre-main sequence stars, the binary population could be mistaken for a younger population at the same distance. For example, the binary counterpart of a population with age t ~ 12 Myr falls in the same locus of the G−J versus G color magnitude diagram as a population with age t ~ 7 Myr. This means that the fit could mistake the unresolved binaries for a younger population, therefore the interpretation of Fig. 11 requires some care. Another caveat is related to the definition of the Gaia G band in the PARSEC libraries. Indeed, the nominal Gaia G passband (Jordi et al. 2010) implemented in the PARSEC libraries is different from the actual passband (cfr. Carrasco et al. 2016). This affects the values of G and G-J predicted by the PARSEC libraries and therefore our absolute age estimates, but does not influence the age ordering. The same can be said for the extinction. Choosing a different (constant) extinction value shifts the isochronal tracks and therefore the estimated age is different, but does not modify the age ranking. In conclusion, the age ranking we obtain is robust, and, even with all the aforementioned cautions, Fig. 11 shows the potential of producing age maps for the Orion region.

4. Orion in Pan-STARRS1

To confirm the age ordering we obtain with Gaia DR1, we applied the analysis described in Sect. 3 to the recently published Pan-STARRS1 photometric catalog (Chambers et al. 2016; Magnier et al. 2016).

Pan-STARRS1 has carried out a set of distinct synoptic imaging sky surveys including the 3π Steradian Survey and the Medium Deep Survey in five bands (grizy). The mean 5σ point source limiting sensitivities in the stacked 3π Steradian Survey in grizy are (23.3, 23.2, 23.1, 22.3, 21.4) magnitudes, respectively. For stars fainter than r ~ 12 mag, Pan-STARRS1 and Gaia DR1 photometric accuracies are comparable. Stars brighter than r ~ 12 mag have large photometric errors in the PanSTARRS filters, therefore we decided to exclude such stars from our sample. We considered the same field defined in Eq. (1)and we performed a cross-match of the sources with Gaia DR1 and 2MASS, using a cross-match radius of 1′′. We did not account for proper motions, since the mean epoch of the Pan-STARRS1 observations goes from 2008 to 2014 for the cross-matched stars and therefore the cross-match radius is larger than the distance covered in the sky by any star moving with an average proper motion of a few mas yr^-1. We obtained N = 88 607 cross-matched sources and we analyzed this sample with the same procedure explained in Sect. 3. Briefly, we first excluded the bulk of the field stars making a cut in the r−i versus r color magnitude diagram, (7)Then we performed the same JHK photometric selection as in Eq. (5), and we studied the on-sky distribution of the sources. We find some density enhancements, corresponding to those already investigated only with the Gaia DR1. We then smoothed the data point distribution in Galactic coordinates using a Gaussian kernel with bandwidth 0.3°. We selected all the sources within the S = 2 density levels and we estimated the single stellar ages with the same Bayesian fitting procedure described above. In this case, however, we did not use the Gaia and 2MASS photometry, but the r and i Pan-STARRS1 bands.

Figure 12 shows the on-sky distribution of the sources with similar ages. The age intervals used are the same as in Fig. 11. The density enhancements corresponding to known groups are visible. Moreover, by comparing Figs. 11 and 12, one can immediately notice that the same groups appear in the same age intervals except for the ϵ Ori group, which appears slightly older than with Gaia DR1 photometry. Indeed the ϵ Ori density enhancement peaks in 15 <t< 20 Myr with PanSTARRS photometry, while it is spread between 11 <t< 20 Myr with Gaia DR1. Another interesting feature of the Pan-STARRS1 age maps are the density enhancements below ϵ Ori. These structures appear prominently in the oldest age panels and might be related to the Orion X population (Bouy & Alves 2015).

These results strengthen our confidence in the age estimates obtained with Gaia photometry, in particular regarding age ordering.

5. Discussion

The present analysis confirms the presence of a large and diffuse young population toward Orion, whose average distance is d ~ 380 pc. The ages determined in Sect. 3.4 show that the groups are young (age <20 Myr) and not coeval. The age ranking determined using Gaia and 2MASS photometry (Fig. 8) is consistent with that determined using Pan-STARRS1 (Fig. 12).

Figures 9, 11, and 12 show some important features, which can potentially give new insights into our understanding of the Orion region.

The Orion dust ring. As already mentioned in Sect. 3.3, a number of overdensities are present toward the Orion dust ring discovered by Schlafly et al. (2015). The age analysis is not conclusive since many overdensities are not within S = 2. Unfortunately, there are no proper motions and/or parallaxes available for these sources (nor in Gaia DR1 nor in other surveys), and their distribution in the color magnitude diagram is not very informative. Additional clues about their origin will be hopefully provided by Gaia DR2.

The Orion Blue-stream. Bouy & Alves (2015) studied the three-dimensional spatial density of OB stars in the solar neighborhood and found three large stream-like structures, one of which is located toward l ~ 200° in the Orion constellation (Orion X). Figure 13 shows the position of the candidate members of the Orion X group as blue stars. Even though the candidate member center looks slightly shifted with respect to the density enhancements shown in the map, it is difficult to argue that these stars are not related to the young population we analyzed in this study. Bouy & Alves (2015) reported that the parallax distribution of the Orion X sources goes from ϖ ~ 3mas to ϖ ~ 6 mas (150 <d< 300 pc), which indicates that Orion X is in the foreground of the Orion complex. Bouy & Alves (2015) also proposed that the newly discovered complex could be older than Orion OB1 and, therefore, constitute the front edge of a stream of star formation propagating further away from the Sun.

To test this scenario we proceeded as follows. First we complemented the bright end of TGAS with Hipparcos data, and then we selected the stars using the proper motion criterion of Eq. (3)and with 3 <ϖ< 7 mas. In this way we restricted our sample to the stars probably kinematically related to the Orion OB association, but on average closer to the Sun. The density of the distribution of these sources in the sky is shown in Fig. 13, together with the Orion X candidate members. We selected the sources within the S = 2 levels (with S defined in Sect. 3), and we used the Bayesian isochronal fitting procedure to estimate the age of this population. Out of the 48 Orion X candidate members listed in Bouy & Alves (2015), only 22 are included in TGAS (the others are probably too bright). To perform the isochronal fit, we could actually use the measured parallax instead of one single value. The age distribution for the foreground sources is shown in Fig. 14 (orange histogram). As a comparison, the age distribution of the sources within the density enhancements and with 2 <ϖ< 3.5 mas is also shown (blue histogram). On average, the foreground population looks older, which is consistent with the picture that Bouy & Alves (2015) proposed. There are however two caveats:

• the age distributions are broad;

• the parallax errors are large and dominate the age estimate.

With future Gaia data releases we will be able to further study the Orion X population and more precisely characterize it.

25 Ori. As pointed out in Sect. 3.3 the 25 Ori group presents a northern extension (~ 200°,−17°) visible in the TGAS, Gaia DR1, and Pan-STARRS1 density maps. The northern extension parallax is only slightly larger than that of the 25 Ori group, and the age analysis suggests that the groups are coeval. With a different approach, Lombardi et al. (2017) have found evidence of the same kind of structure (see their Fig. 15). Gaia DR2 will be fundamental in discerning the properties of this new substructure of the 25 Ori group.

Theλ Ori group. In Sect. 3.3 we pointed out some overdensities located on the Hα bubble surrounding λ Ori, which are not related to known groups (to our knowledge). We further investigated the stars belonging to these overdensities, however, there are no parallaxes nor proper motions available for these sources and it is difficult to draw firm conclusions from the photometry only (also combining Gaia DR1 and Pan-STARRS1). In this case as well, we have to conclude that hopefully Gaia DR2 will clarify whether these groups are real or not.

NGC 1980 and the ONC. One of the most interesting features of the maps of Figs. 11 and 12 is the prominent density enhancement toward NGC 1980, NGC 1977, and the ONC. The density enhancement is not only concentrated in one panel, but persists throughout until it disappears in the last panel. This can be explained in at least two ways:

• There are multiple populations at roughly the same distance, withdifferent ages.

• There is only one population with a single age, however its spread along the line of sight is so large that using only one parallax value for the fit is not accurate enough.

Both explanations have supporters. Alves & Bouy (2012) suggested that NGC 1980 is not directly related to the ONC, i.e., these two young associations are not the same population emerging from the parental cloud but are instead distinct overlapping populations. On the other hand, based on the fact that the kinematic properties of NGC 1980 are indistinguishable from those of the rest of the population at the same position in the sky, Da Rio et al. (2016) have argued that NGC 1980 simply represents the older tail of the age distribution around the ONC in the context of an extended star formation event. Using isochronal ages, Fang et al. (2017) have found that the foreground population has a median age of 1−2 Myr, which is similar to that of the other young stars in Orion A. Furthermore, they confirmed that the kinematics of the foreground population is similar to that of the molecular clouds and other young stars in the region. These authors therefore argue against the presence of a large foreground cluster in front of Orion A. Kounkel et al. (2017a) have estimated that the age of NGC 1980 is ~ 3 Myr, which is comparable with the study by Fang et al. (2017), however they were not able to confirm or disprove whether NGC 1980 is in the foreground on the ONC. Finally, Beccari et al. (2017) have discovered three well-separated pre-main sequences in the r−i versus r color magnitude diagram obtained with the data of the wide field optical camera OmegaCAM on the VLT Survey Telescope in a region around the ONC. These sequences can be explained as a population of unresolved binaries or as three populations with varying ages. The populations studied by Beccari et al. are unlikely to be related to NGC 1980, however, if confirmed, they would constitute an example of non-coeval populations in the same cluster. Figure 11 shows that the group corresponding to NGC 1980 is well defined not only at very young ages (1 <t< 3 Myr), but at least until t ~ 15Myr. We discuss below the influence that unresolved binaries have on our age determination (indeed our fit does not account for them); the main point is that unresolved binaries influence the youngest age intervals, not the oldest. This would point toward the actual existence of two populations: the first related to the ONC and the second to the Alves & Bouy (2012) foreground population.

In conclusion, the ages of the stellar populations toward Orion show a gradient, which goes from 25 Ori and ϵ Ori toward the ONC and the Orion A and B clouds. The age gradient is also associated with a parallax gradient. Indeed, the older population toward 25 Ori and ϵ Ori is also closer to the Sun than the younger population toward the ONC (see also Fig. 6). Gaia DR2 will provide distances to the individual stars of each different group, and we will therefore be able to obtain more precise ages for these populations as well.

To study whether or not the parallax gradient influences the age determination, we performed the same Bayesian isochrone fit changing the parallax of each star according to its position, following Fig. 6. We also performed the analysis including a uniform prior on the parallax distribution and then marginalizing over the parallax. In both cases, the estimated ages for the single groups have some small variations, however our conclusions do not significantly change.

To test how our result depends on the set of isochrones we chose, we performed the fit again, using the MESA Isochrones and Stellar Tracks (MIST; Dotter 2016; Choi et al. 2016). We fixed the metallicity to Z_⊙ and we applied the usual extinction correction of A_V = 0.25 mag. Whilst in this case the single ages are in general estimated to be younger than the estimates based on the PARSEC models (e.g., the 25 Ori group peaks between 9 <t< 11 Myr), age ordering does not change significantly.

Finally, we studied the distribution in the sky of the coeval sources fainter than 14 mag. In this magnitude range we can remove the sources that are most likely Galactic contaminants. We found again the same groups and the same age ordering.

As mentioned above and in Sect. 3.4, the unresolved binary sequence could stand out as a separate, seemingly younger population, which would add further complications to the age determination of the group. The Bayesian fitting procedure does not take into account the presence of unresolved binaries. The net effect of this is that the unresolved binaries population is mistaken for a younger population. For example, the difference in magnitudes between the 5−7 Myr and the 13−20 Myr isochrones corresponds almost exactly to the 0.75 mag separating the primary sequence from the unresolved binary sequence. This is a major cause of age spread and it could greatly affect our age estimates, thus it appears even clearer that great care needs to be used when analyzing them. On the other hand, however, binary should affect all populations in the same way. This further supports the robustness of our relative age estimates.

Another intriguing problem is related to the relation between the density enhancements, diffusely distributed massive stars, and gas distribution. Figure 15 shows the S = 3,6, and 9 contour levels of the overdensities on top of an extinction map obtained with Planck data (Planck Collaboration XI 2014) probing the dark clouds. The older group 25 Ori is located far away from the gas, while the younger groups of λ Ori, σ Ori, ϵ Ori, NGC 1977, and NGC 1980 closely follow the clouds. Orion A and B are behind the density enhancements. The three-dimensional structure of the region is still unclear, and current data accuracy is not yet good enough to draw definite conclusions, especially at the distance and direction of Orion. The data quality however will improve in future Gaia data releases, and likewise our understanding of the region. In particular, precise parallaxes, proper motions, and radial velocities will allow us to address directly the recent discovery that the Orion clouds might be part of an ancient dust ring (Schlafly et al. 2015), the blue streams scenario proposed by Bouy & Alves (2015), and the complex nested shell picture unveiled by Ochsendorf et al. (2015).

6. Conclusions

In this paper we made use of Gaia DR1 (Gaia Collaboration (2016a,b; van Leeuwen et al. 2017) to study the stellar populations toward Orion. Our results are as follows:

• Using TGAS (Michalik et al. 2015; Lindegrenet al. 2016) we found evi-dence for the presence of a young population, at a parallaxϖ ~ 2.65 mas (d ~ 377 pc), loosely distributed around the following known clusters: 25 Ori, ϵ Ori and σ Ori, and NGC 1980, and ONC. The stars belonging to these groupings define a sequence in all the color magnitude diagrams constructed by combining Gaia DR1 and 2MASS photometry.

• We considered the entire Gaia DR1, again realizing color magnitude diagrams combining Gaia and 2MASS photometry for the entire field. Well visible between G = 14 mag and G = 18 mag, we found the low mass counterpart of the sources isolated with TGAS.

• After a preliminary selection to exclude field stars, we studied the distribution in the sky of the sources belonging to this sequence using a KDE. We found density enhancements in the sky distribution comparable to those in the TGAS sample.

• We estimated the ages of the sources within the density enhancements, using a Bayesian isochrone fitting procedure described in detail in Jørgensen & Lindegren (2005). We assumed all the stars to be at the same parallax, ϖ = 2.65. We found that the groupings have different ages. In particular, there is an age gradient going from 25 Ori (13−15 Myr) to the ONC (1 Myr).

• To consolidate our findings, we repeated the fitting procedure using the sources in common with the Pan-STARRS1 (Chambers et al. 2016; Magnier et al. 2016) r and i filters, finding the same age ordering as with Gaia DR1.

• We studied the distribution in the sky of the groups we found. In particular:

  • 1.

    The 25 Ori cluster presents a northern extension, which wasalso reported by Lombardiet al. (2017).

  • 2.

    Some of the density enhancements toward the λ Ori complex are related to known clusters (Col 69, B30, and LDN 1588), but some other overdensities on the left of the ring are new. Unfortunately, it was not possible to investigate these further since we have neither precise proper motions nor parallaxes.

  • 3.

    Some overdensities are also present within the Orion dust ring discovered by Schlafly et al. (2015), and these might be related to the star formation process out of which the ring was formed. In this case as well, however, more data are needed to confirm our speculations.

  • 4.

    The Orion X candidate members (Bouy & Alves 2015) are related to some of the density enhancements shown in Fig. 11. We studied the sky and age distribution of the TGAS sources with proper motions as in Eq. (3)and parallax 3 < ϖ < 7 mas, and we found that the stars with 2 < ϖ < 3.5 mas are on average younger than those with 3 < ϖ < 7 mas.

• We discussed the implications of the age ranking weobtained. We found that the estimated ages toward theNGC 1980 cluster span a broad range ofvalues. This could be due to either the presence of twopopulations coming from two different episodes of starformation or to a large spread along the line of sight of thesame population. Some confusion might also arise fromthe presence of unresolved binaries, which are not mod-eled in the fit and usually stand out as a younger popula-tion. We related our findings to previous works by Bouyet al. (2014), DaRio et al. (2016), and Fanget al. (2017).

• Finally, we link the stellar groups to the gas and dust features in Orion, albeit in a qualitative and preliminary fashion. Future Gaia releases will allow to address these questions in unparalleled detail.

Acknowledgments

We are thankful to the anonymous referee for comments that greatly improved the manuscript. This project was developed in part at the 2016 NYC Gaia Sprint, hosted by the Center for Computational Astrophysics at the Simons Foundation in New York City, and at the 2017 Heidelberg Gaia Sprint, hosted by the Max-Planck-Institut für Astronomie, Heidelberg. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC; https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This publication has made use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. The Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, the Queen’s University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation Grant No. AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation. This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration 2013). This work has made extensive use of IPython (Pérez & Granger 2007), Matplotlib (Hunter 2007), astroML (Vanderplas et al. 2012), scikit-learn (Pedregosa et al. 2011), and TOPCAT (Taylor 2005, http://www.star.bris.ac.uk/~mbt/topcat/). This work would have not been possible without the countless hours put in by members of the open-source community all around the world. Finally, C.F.M. gratefully acknowledges an ESA Research Fellowship.

References

Appendix A: color magnitude and color–color diagrams

In this appendix we show the color-color and color magnitude diagrams constructed combining Gaia DR1 and 2MASS photometry. The sources in the first panel are those remaining after

applying the 2MASS photometry quality selection cut (“ph_qual = AAA”). The other panels show the cuts of Eq. (2). We did not apply exactly the same photometric criteria as in Alves & Bouy (2012) because there is probably a typo in their Eq. (1) that causes 0 sources to be selected. However, Fig. A.1 looks similar to their Fig. 4.

Appendix B: ADQL queries

We report here the queries used to a) select the sources in our field and b) perform the cross-match with 2MASS.

Field selection:

 
select gaia.source_id, gaia.ra, gaia.dec, gaia.l, gaia.b, gaia.phot_g_mean_mag, gaia.pmra, 
gaia.pmdec, gaia.parallax, gaia.pmra_error, gaia.pmdec_error, gaia.parallax_error 
from gaiadr1.gaia_source as gaia 
where gaia.l>=190.0 and gaia.l<=220.0 and gaia.b>=-30.0 and gaia.b<=-5.0
         

Cross-match with 2MASS:

 
select gaia.source_id, gaia.l, gaia.b, gaia.phot_g_mean_mag, gaia.phot_g_mean_flux, 
gaia.phot_g_mean_flux_error, gaia.parallax,gaia.parallax_error, 
gaia.pmra,gaia.pmdec,gaia.pmra_error,gaia.pmdec_error, 
tmass.j_m, tmass.j_msigcom, tmass.h_m, tmass.h_msigcom, tmass.ks_m, tmass.ks_msigcom, 
tmass.ph_qual from gaiadr1.gaia_source as gaia 
inner join gaiadr1.tmass_best_neighbour as xmatch 
on gaia.source_id = xmatch.source_id 
inner join gaiadr1.tmass_original_valid as tmass 
on tmass.tmass_oid = xmatch.tmass_oid 
where gaia.l > 190.0 and gaia.l < 220.0 
and gaia.b < -5.0 and gaia.b > -30.0 and xmatch.angular_distance < 1.0 
         

We run the queries using the ESA Gaia archive (https://archives.esac.esa.int/gaia). On the archive, we suggest that the user creates a personal account. This indeed allows the user to save queries and store data (up to 1 GB).

Appendix C: Kernel density estimation on the sphere

The referee pointed out that the kernel density estimation carried out on flat projections of the Orion sky field will suffer from area distortions, and suggested the use of the von Mises-Fisher (vMF) kernel, which is intended for analyses on the unit sphere. This kernel is given by the following equation for a two-dimensional unit sphere (i.e., for three-dimensional unit vectors x): (C.1)where the unit vector m represents the mean direction for the kernel and is the modified Bessel function of the first kind and order 1/2, and m^Tx indicates the inner product of m and x. See Banerjee et al. (2005) for details. The parameter κ ≥ 0 indicates the concentration of the kernel around the direction indicated by m. The normalization constant can be rewritten by considering that (C.2)(Eq. (10.39.1) in Olver et al. 2010), which leads to:(C.3)The exponent in the kernel contains the inner product of m and x and this can also be written as (C.4)where (α,δ) represent the ICRS coordinates of the points on the sky, and ρ is the angle between m and x. On the unit sphere this angle also represents the distance along a great circle between the points m and x, also known as the “haversine distance”. The value of ρ can also be calculated using the haversine function (hav) given by (C.5)The formula for ρ then becomes (C.6)\newpage

This can be verified by writing out both sides of the equation in terms of (1−cosθ)/2.

To continue, we note that the half width at half maximum of the vMF kernel expressed in terms of ρ (ρ_HWHM) is given by (C.7)where Eq. (C.7) follows from (C.8)as the maximum of f(x | m,κ) occurs when cosρ = 1.

Equivalently, for a given ρ_HWHM the corresponding value of κ is(C.9)In our kernel density estimates of source distributions on the sky the kernel sizes are of order 1 degree (0.017 radians) or less. This is already in the regime where to good accuracy cosρ ≈ 1−ρ²/ 2. At the same time the value of κ becomes very large (~ 4550; see Eq. (C.9) for ρ_HWHM = 0.017), such that sinhκ → exp(κ)/2. Hence the vMF kernel becomes approximately (C.10)This is in fact a two-dimensional normal distribution with standard deviations along the two principal axes, where in the small angle regime one can write ρ² = (Δαcosδ)² + Δδ², with Δα = α−α_m and Δδ = δ−δ_m. This shows that in our case (with kernel sizes of a degree or less), the vMF kernel can be approximated as a two-dimensional Gaussian in terms of the haversine distance.

Our implementation of the kernel density estimate is in Python and makes use of the sklearn.neighbors. KernelDensity module in the scikit-learn package by specifying that the “haversine” metric should be used during the fitting stage of the density estimate (using the parameters kernel=“Gaussian” and metric=“haversine”).

All Tables

All Figures