© ESO 2019

1. Introduction

The tendency of O and B type stars to loosely cluster in the sky was recognised at the beginning of the 20th century by the pioneering studies summarised in Blaauw (1964). At the end of the last century, the data of the HIPPARCOS satellite allowed de Zeeuw et al. (1999), de Bruijne (1999), Hoogerwerf & Aguilar (1999), and many others to characterise the stellar content and the kinematic properties of nearby OB associations. OB associations have been long considered as expanding remnants of your star clusters as their members are widely dispersed over the sky (Brown et al. 1999; Lada & Lada 2003). The classical explanation for this is that star clusters are formed embedded within molecular clouds, where the gravitational potential of both the stars and the gas holds them together. When feedback disperses the gas left over from star formation, the cluster becomes supervirial and will expand and disperse, thus being visible for a short time as an OB association. While many observations support this model (Lada & Lada 2003, and references therein), it has been difficult to test whether OB associations are indeed expanding. Wright et al. (2016) and Wright & Mamajek (2018) studied the kinematics of the Cygnus OB2 and Scorpius-Centaurus associations respectively, and concluded that they were not formed by the disruption of individual star clusters. Wright & Mamajek (2018) further concluded that Sco-Cen was likely born highly sub-structured, with multiple small-scale star formation events contributing to the overall OB association, and not as a single, monolithic burst of clustered star formation. These conclusions can be related to the fact that the distribution of young stars within their parental molecular clouds is fractal, hierarchical, and follows the filamentary structures of the dense gas, both spatially (Gutermuth et al. 2008) and kinematically (Hacar et al. 2016). Clusters then form where filaments overlap (Myers 2009; Schneider et al. 2012; Hacar et al. 2016, 2017): their formation might be due to higher column densities or to the merging of filaments that have already formed stars. OB associations would therefore constitute the final stage of this star formation mechanism. They still keep memory of the parental gas sub-structure where they originated as they slowly disperse in the field.

At a distance of ∼380 pc (Zari et al. 2017), the Orion star forming region is the nearest site of active high-mass star formation. It is a benchmark for studying all stages and modes of star formation (Brown et al. 1994; Jeffries et al. 2006; Bally 2008; Briceno 2008; Muench et al. 2008; Da Rio et al. 2014, 2016; Getman et al. 2014; Hacar et al. 2016; Kubiak et al. 2017; Fang et al. 2017; Kounkel et al. 2017), in addition to the effect of star formation processes on the surrounding interstellar medium (Ochsendorf et al. 2015; Schlafly et al. 2015; Soler et al. 2018). Zari et al. (2017) used Gaia DR1 (Gaia Collaboration 2016a,b) to study the density distribution of the young, non-embedded stellar population in the sky, and obtained a first picture of the star formation history of the Orion region in terms of the various star formation episodes, their duration, and their effects on the surrounding interstellar medium. Even though proper motions where available for the Tycho-Gaia astrometric solution (TGAS, TGAS, Michalik et al. 2015) sub-set of Gaia DR1, they were not accurate enough to perform a precise kinematic analysis. Proper motions in Orion are indeed small as stars move on average radially away from the Sun. Furthermore, to derive the ages of the stellar populations, a single distance value was considered (d ∼ 380 pc) as parallax uncertainties were too large to resolve the spatial configuration of the groups that were identified. By combining the data of the second release of the Gaia satellite (hereafter Gaia DR2, Gaia Collaboration 2018) and APOGEE-2, Kounkel et al. (2018) study the entire Orion complex, providing a classification of the stellar population in five groups, and an analysis of their ages and kinematics. Kos et al. (2018) use Gaia DR2 parameters supplemented with radial velocities from the GALAH and APOGEE surveys to perform a clustering analysis towards the 25 Ori cluster region. They find that one cluster is significantly older (21 ± 2 Myr) compared to the rest of the region. Großschedl et al. (2018) investigate the 3D shape and orientation of the Orion A molecular cloud by analysing the distances of mid-infrared selected young stellar objects, and find that the cloud is elongated and oriented towards the galactic plane, and presents two different components one dense and star forming and one ∼75 pc long, more diffuse and star-formation quiet.

In this work, we use Gaia DR2 to study the three dimensional (3D) structure of the Orion OB association, we model the kinematics of the sub-groups that constitute it and we give estimates of their ages, to obtain a complete picture of the star formation history of the region and to put it in the broader context of the theories of sequential and triggered star formation. In Sect. 2 we present the data and describe how we select the young stellar population in Orion. In Sect. 3 we study its 3D configuration in Cartesian galactic co-ordinates, and we isolate young groups by making use of the DBSCAN clustering algorithm. In Sect. 4 we perform the kinematic analysis by using a maximum likelihood approach. In Sect. 5 we derive ages and extinctions of all the groups resulting from the analysis of Sect. 4. In Sect. 6 we discuss our findings. The conclusions of this work are summarised in Sect. 7.

2. Data

Following Zari et al. (2017), we select the sources with co-ordinates

(1)

and we restrict our sample to the sources with 1.5 <  ϖ <  5.0 mas. Since the Orion population moves mostly radially away from the Sun, we consider only stars with small proper motions:

(2)

We derive distances by inverting parallaxes, d = 1000/ϖ pc thus we restrict our sample to sources with ϖ/σ_ϖ >  5.0, following the recommendations in Bailer-Jones (2015). The effect of this cut is to exclude sources at faint magnitudes (G >  20 mag), but it does not introduce significant biases in the determination of distances to the clusters or the study of their 3D configuration.

2.1. Obtaining a “clean” sample

We apply the following cuts on the photometric and astrometric quality, based on Lindegren et al. (2018) complemented by the information contained on the Gaia known issues page¹. We select all the sources with RUWE < 1.4, following the slides by Lindegren et al². The renormalized unit weight error (RUWE) is defined as:

(3)

where: χ² is the astrometric goodness-of-fit in the AL direction (astrometric_chi2_al); N is the number of good observations AL (astrometric_n_good_obs_al); u₀(C, G) is an empirical normalization factor, which is a function of C = G_BP − G_RP and G. This cut seeks to remove sources with spurious parallaxes or proper motions.

We use the flux excess ratio:

(4)

where I_X is the photometric flux in band X, to exclude sources with possible issues in the BP and RP photometry, affecting in particular faint sources in crowded areas. We apply Eq. (C.2) in Lindegren et al. (2018), which we report here for clarity:

(5)

Evans et al. (2018) and Arenou et al. (2018) mention that Gaia DR2 photometry is affected by some systematic errors. Evans et al. (2018) and Maíz Apellániz & Weiler (2018) propose corrections to mitigate these effects. We apply these corrections and we report them here for clarity:

• 2 ≤ G ≤ 6 mag: G_corr = −0.047344 + 1.16405G − 0.046799G² + 0.0035015G³

• 2 ≤ G ≤ 4 mag: G_BP, corr = G_BP − 2.0384 + 0.95282G − 0.11018G²

• 2 ≤ G ≤ 3.5 mag:

• 6 ≤ G ≤ 16 mag: G_corr = G − 0.0032(G − 6.0)

• G >  16 mag: G_corr = G − 0.032

In the rest of the paper we use the corrected G, G_BP, and G_RP magnitudes without using the subscript “corr”.

2.2. Selecting the young stellar population

Figure 1 (left) shows the M_G vs. G_BP − G_RP colour-magnitude diagram of the “clean” sample obtained in Sect. 2.1. Although faint, the pre-main sequence and the upper main sequence, indicating the presence of the young population in the region, are visible, and can be used to guide the selection of the young stellar populations towards Orion.

Fig. 1. Observed colour-magnitude diagram (left) and sky distribution (right) of the sources selected in Sect. 2. The solid black lines in the left panel isochrones defined in Eq. (6), which are used to select the young stellar population in Orion.

To select young stars, we use the PARSEC isochrones (Bressan et al. 2012; Tang et al. 2014; Chen et al. 2014) with A_V = 0.3 mag and age τ = 10 Myr to define the following region in the M_G vs. G_BP − G_RP colour-magnitude diagram (solid black lines in Fig. 1):

(6)

where the sub-script Iso indicates the isochrone values for G_BP, G_RP, and M_G. We choose A_V = 0.3 mag following Zari et al. (2017). The distribution in the sky of the sources selected in this fashion is shown in Fig. 1 (right). The regions in which we divide the field are also indicated, together with the sub-groups in which the Orion OB1 association is classically split: Orion OB1a, OB1b, OB1c, and OB1d. The same groups identified in Zari et al. (2017) and Kounkel et al. (2018) are visible, which confirms the correctness of the selection.

In Sect. 4 we focus on the kinematics of the Orion population. To complement the Gaia DR2 radial velocities we cross-matched our sources with the APOGEE DR14 catalogue (Abolfathi et al. 2018). The APOGEE synthetic heliocentric velocities (SYNTHVHELIO_AVG, an average of the individual measured RVs using spectra cross-correlations with single best-match synthetic spectrum) were used.

3. 3D distribution and identification of clusters

We first studied the three-dimensional (3D) distribution of sources using a similar approach as in Zari et al. (2018). In summary, we:

1. compute galactic Cartesian co-ordinates for all the sources, x_g, y_g, z_g;

2. define a volume, V = (800, 800, 350), centred in the Sun, and we divide it in 3 × 3 × 3 pc cubes;

3. compute the number of sources in each cube;

4. compute the source density D(x_g, y_g, z_g) by smoothing the distribution with a Gaussian filter, with width w = 2 pc;

5. normalise the density distribution from 0 to 1 by applying the sigmoidal logistic function:

   (7)

   with L = 2, κ = 4 pc, and D₀ = 0.

Figure 2 shows the density distribution of sources f(D) on the galactic plane for different values of z_g. Different density enhancements are visible, corresponding to well known-clusters. The first and second panel show stars in the Orion A molecular cloud. The Orion Nebula Cluster (ONC) corresponds to the most prominent density enhancement. The third panel is particularly interesting because it clearly shows the presence of a foreground population to the ONC, confirming the conclusions by Bouy et al. (2014). Some clusters corresponding to the Belt region also become visible, although the bulk of the population is located between Z = −116 pc and Z = −101 pc. The last three panels mainly show the λ Ori cluster. At Z = −92 pc the northern elongation of the 25 Ori group is visible. The density distribution looks elongated towards the line of sight: this is an effect of the parallax errors. The parallax error distribution is peaked at σ_ϖ = 0.046 mas, but presents a long tail towards larger values (the 84th percentile is 0.11 mas).

Fig. 2. Density distribution of the sources in Orion for different Z values. The orange stars indicate the positions of (from top to bottom in each panel): λ Ori, 25 Ori, ϵ Ori, and the Orion Nebula cluster. The white solid contours represent the 0.4, 0.5, 0.6 and 0.8 density levels (the density is normalised to have values from 0 to 1). The Sun is located at (X, Y) = (0, 0).

To isolate the members of each cluster, we considered only the sources within the density level f(D) = 0.5 of the 3D map shown in Fig. 2. This value is arbitrary and aims at selecting the densest regions of the maps. The clusters were then separated by using the DBSCAN algorithm³. As described for instance by Price-Jones & Bovy (2019), DBSCAN is a density-based clustering algorithm that views clusters as areas of high density separated by areas of low density in space, without requiring any prior assumption on the number of groups present. There are two parameters to the algorithm, min_samples and eps, which define the density of the clusters. Higher min_samples or lower eps values indicate higher densities necessary to form a cluster.

Clusters in Orion have different sizes and numbers of members, and therefore different densities: for this reason we applied the clustering algorithm twice. The first time we used min_samples=50 and eps=7 pc to isolate the main structures, shown in Figs. 3 and 4 (top), obtaining five groups. The group that encompasses 25 Ori, the Belt region and the Orion A foreground can be visibly divided in sub-groups. Thus we apply DBSCAN only to this group with different parameters: we find that min_samples=30 and eps=5 pc are the best values to separate all the sub-clusters (see Figs. 3 and 4, bottom).

Fig. 3. Distribution of the groups identified with the first (top) and second (bottom) iteration of DBSCAN in the planes (X, Y), (X, Z), and (Y, Z).

Fig. 4. Sky distribution of the groups identified with the first (top) and second (bottom) iteration of DBSCAN. The colours correspond to those in Fig. 3.

This method has the drawback of excluding stars that might be related to the star formation events in Orion, but are more dispersed than the rest of the population in 3D space (but could still be compact in proper motion space). This is further discussed in Sect. 6.

4. Kinematics

In this section we study the kinematics of the groups selected in the previous section. We used an iterative maximum likelihood approach to determine (a) the average motion of the groups, (b) their velocity dispersion, and (c) (where possible) the presence of a linear expansion term. We used the method proposed by Lindegren et al. (2000) and applied in Reino et al. (2018) and Bravi et al. (2018), adding however a term to take into account a potential expansion of the cluster from its centre. The method is summarised in Sect. 4.1, tested in Appendix A, and the results are presented in Sect. 4.2. Here we use ICRS co-ordinates, which we differentiate from galactic co-ordinates by adding the subscript “I” when needed.

4.1. Method

Our method extends the maximum-likelihood method developed by Lindegren et al. (2000, L00) by adding measured radial velocities (see Reino et al. 2018) and by including a linear expansion term in the cluster velocity model. Following L00, we assume that the members of a cluster share the same three-dimensional space motion with a small isotropic dispersion term. Reino et al. (2018) extended L00’s method by:

– adding measured radial velocity, whenever available, as a fourth observable, besides trigonometric parallax and proper motion;

– making a transition from the χ² statistic used in L00, and denoted g, to a p value or 1 − CDF(g,DOF) as a goodness-of-fit statistic;

– using a mixed three- and four-dimensional likelihood function so that both stars with and without known radial velocity can be treated simultaneously.

Following L00, we include a linear expansion term in the cluster velocity model by writing the expected space velocity of a single star at position b_i as:

(8)

where b₀ is an arbitrary reference position, namely the point where the local velocity u(b) assumes the status of “centroid” velocity v₀. The co-ordinates of b₀ are therefore fixed in advance. The matrix T is simply a diagonal matrix of the form:

An expanding cluster will have κ >  0, from which an expansion age, τ = 1/(γκ) [Myr] can be derived (γ is a conversion factor of 1.0227 pc Myr⁻¹ km⁻¹ s, see for example Wright & Mamajek 2018).

The method is applied to the members of the clusters identified in Sect. 3. These clusters still contain “outliers”, that is real non-members, or members which have (slightly) discrepant astrometry (or radial velocities) as a result of unrecognised multiplicity, them escaping from the cluster, etc. Such outliers can be found, after maximising the likelihood function, by computing the p value (associated with a particular g value) for each star in the solution (Eq. (19) in L00). The largest outlier is removed from the sample and a new maximum likelihood solution is determined, until all g values are acceptably small (g_i ≤ g_lim or p_i ≥ p_lim). The stopping criterion is the same as in Reino et al. (2018), and is associated to a significance level p_lim = 0.0027. As noted in Reino et al. (2018), if one stops too early, real outliers will be left and the best-fit velocity dispersion will remain too high. On the contrary, one can keep on iterating and removing outliers until just two stars with very similar three-dimensional motions are left, severely underestimating the velocity dispersion. Astrometric data only can not distinguish between expansion or contraction of a cluster from a change in v₀ (see L00). Therefore when the fraction of measured radial velocities is lower than the 20% we do not estimate the expansion coefficient κ (implicitly assuming κ = 0). The threshold is conservative for certain groups, but the derived parameters are robust for all the groups.

4.2. Results

The results of the kinematic modelling code are given in Table 1. Being quite isolated with respect to the rest of the population, the λ Ori group (group A) is easy to identify and separate from the others, therefore the results do not require any specific clarification. This is not the case for the groups with 199 ° < l <  216°. We comment on the results for these groups by dividing them in three “regions” according to their sky distribution: the 25 Ori region, the Belt region, and the Orion A region.

4.2.1. 25 Ori

We define the 25 Ori region as:

(9)

which corresponds to the groups B₀ and B₆ identified by DBSCAN. The proper motions of the sources in the region (black dots in Fig. 5, left) separate in two clumps. This was shown also by Kos et al. (2018), who however apply a different classification scheme to separate the clusters in the region. The separation is also visible when considering the proper motion diagram of group B₀ (blue dots in Fig. 5, left). The number of sources is lower because the DBSCAN algorithm favours the high density groups (so when the density drops under a certain level the stars are considered as “noise stars” and not classified as members of any cluster).

Fig. 5. Left: proper motion diagram of the stars in the 25 Ori region defined in the text (black dots), of the stars belonging to group B0 (blue dots), and of the stars belonging to group B0, b (light blue empty squares) and B6 (pink crosses). Right: sky distribution in galactic co-ordinates of group B0, B0, b, and B6. The colours and symbols are the same as on the left.

We considered the sources selected by DBSCAN, and we isolated the second group (B_0, b, light blue squares in Fig. 5, left) by applying the following cuts in proper motion space:

(10)

We applied separately the kinematic modelling code to the two groups. The results are reported in Table 1. We also run the kinematic modelling code considering all the sources in the region, after separating the two groups using the same criteria of Eq. (10). The estimated parameters are consistent. The sky distribution of the sources of group B₀ and B_0, b is shown in Fig. 5 (right panel). While group B₀’s distribution shows a clump towards 25 Ori, and the northern elongation reported for instance by Lombardi et al. (2017) and Briceño et al. (2019), group B_0, b’s sources are scattered in the field and do not show any clear concentration. Together with the findings by Kos et al. (2018) in terms of ages (see also Sect. 5), this points to the conclusion that group B_0, b is slowly dispersing in the galactic field. Here we are limiting our samples to the 25 Ori region, but in principle members of the B_0, b group could be found spread over a larger area of the sky (and 3D space).

Group B₆ consists only of 30 members, none of which has a measured radial velocity, therefore we decided not apply the kinematic modelling code. The parallax distribution suggests that B₆ is closer to the Sun than group B₀, while the proper motion distribution does not show any difference with respect to group B₀. We suspect that group B₆ coincides with a small over-density of sources within group B₀, which gets classified as a separate group because of a local density drop. We ran the kinematic modelling code for groups B₀ and B₆ together: the estimated parameters are consistent with those found for group B₀ only, which supports our hypothesis.

4.2.2. Belt

Many of the clusters identified by DBSCAN (B₂, B₃, B₄, B₅, B₇ and E) are located in the Sky towards the Belt region. Figure 6 shows the proper motion diagram for the Belt region defined as

(11)

Fig. 6. Proper motion diagram for all the stars in the Belt region (grey dots). Each panel corresponds to one of the groups identified by DBSCAN (the colours are the same as in Figs. 3 and 4) except for group B7, which is indicated by black dots for representation purposes.

Proper motions in the Belt region present a high degree of sub-structure, indicating that the Belt hosts groups with different kinematic properties.

– Groups B₂ and B₄ are mostly located towards the σ Ori cluster (see Fig. 7) and ζ Ori. Group B₃’s members are spread towards ϵ Ori and δ Ori. The parameters estimated by the kinematic modelling code suggest that B₂ and B₄ have compatible v_y, I values, which are significantly different from those of group B₃. This is consistent with what is found by Jeffries et al. (2006), who already notice the presence of two kinematics components towards the cluster. The kinematic properties of group B₃ are similar to those of groups D, B₀ (not located in the Belt region, see Fig. 4), and B₅. We notice that group B₂’s velocity dispersion is large (∼1.6 km s⁻¹) compared for instance to that of group B₃ (0.41 ± 0.02 km s⁻¹). The proper motion distribution shows indeed some substructures, which cause the large value of the velocity dispersion. As mentioned above, the presence of kinematic substructure may indicate the co-existence of groups with different kinematics in the same area. An inspection of group B₂’s 3D configuration (see Fig. 3, in particular the X − Y projection) shows that the source distribution is not uniform, and seems to be divided into (at least two) elongated structures.

Fig. 7. Distribution in the sky of the stars in the Belt region. Each panel corresponds to one group identified by DBSCAN. The colour-coding is the same as in Fig. 6. Orange stars correspond to, from left to right: σ Ori, ζ Ori, ϵ Ori, η Ori, δ Ori.

– Group B₅ is located below the Belt, towards η Ori, and shares similar kinematics with group B₃, although they seem to be well separated in space (see Figs. 3 and 4). The proper motion distribution shows two clumps, similar to what is observed towards 25 Ori. We separate the smaller clump, which we refer to as B_5, b by using simple cuts in proper motion space:

(12)

In contrast to what we have done for group B_0, b, here we apply the conditions of Eq. (12) to all the sources in the Belt region, and not just those within the f(D) = 0.5 level of the 3D density map. This is the reason why the number of sources is higher than for group E (see Table 1). This choice is motivated by the fact that the visual inspection of the proper motion diagram suggests that the clump is more extended and the number of sources is larger than what found by DBSCAN. Further, the number of sources of the smaller clump is too small to retrieve the kinematic parameters accurately. The parameters estimated by the kinematic modelling code (see Table 1) show that group B₅ and group B_5, b have different kinematic properties, while having similar parallaxes. Comparing Figs. 6 and 10 one can notice that the region defined in Eq. (12) also includes sources classified as members of group B₂. The sky distribution of sources belonging to group B_5, b (see Fig. 8) shows indeed some sources clustering around σ Ori. Most of the sources however are located in the same region as group B₅, although they are spread throughout the entire longitudinal extent of the Belt region. This seems to suggest that group B_5, b is more extended than the Belt region, especially to lower galactic latitudes and longitudes. Similar conclusions can be drawn after studying the 3D distribution of group B_5, b (Fig. 8): some sources clump in the same area as group B₂ and B₄ (σ Ori), while others are located closer to group B₅. This explains why DBSCAN does not separate successfully groups B₅ and B_5, b: their members show different kinematics but are mixed in space.

Fig. 8. Top panel: sky distribution of the sources belonging to group B5, b (brown empty squares), group B5 (brown dots), and all the sources in the Belt region defined in the text (grey dots). The orange stars mark the position of σ Ori, ζ Ori, ϵ Ori, η Ori, and δ Ori. Bottom panels: distribution in 3D galactic co-ordinates of group B5, b (brown squares), and of all the sources belonging to the Belt region (grey dots).

– Group E is the most distant group in the entire Orion region (see Table 1 and Fig. 3). Since not many radial velocity measurements are available, the kinematic properties are determined with less accuracy than for the other groups, especially in the y_I direction. While v_y, I is comparable with those of group A, C, B₁, B₂, B₄ (and B_7, b and B₈, see below), the v_x, I component is different from the other groups. As for group B₅, the proper motions seem to be divided in two clumps, one of which does not correspond to any other DBSCAN groups. We select group B₈ by applying the following conditions:

(13)

Similarly as for group B_5, b, and with the same motivations, we consider again all the sources in the Belt region. The estimated kinematic parameters are reported in Table 1. The source distribution in the sky and in 3D Cartesian space is shown in Fig. 12, compared to that of group E. The sources are loosely distributed in the entire Belt region, although they seem to clump next to group E.

– DBSCAN identifies only 30 sources belonging to group B₇, none of them with a measured radial velocity, therefore the kinematic modelling code does not succeed in determining reliable parameters. Similarly to what was found for group B_5, b and B₈, when considering all the stars in the Belt area, we notice that many more sources clump in the same proper motion region that are excluded when we apply the condition f(D) > 0.5 or that are classified as “noise” stars by DBSCAN. We therefore select group B_7, b according to the following equations (see Fig. 11):

(14)

The number of sources is now much larger (see Table 1), and the parameters can be accurately determined. Figure 10 shows the source distribution in the sky and in Cartesian galactic co-ordinates. We notice that the sources are distributed in the sky towards the reflection nebulae M 78 and NGC 2071, where two groups of young stars are present and towards the centre of the Belt.

– Figure 9 shows the dust distribution towards the Belt region, where a bubble is visible (see Ochsendorf et al. 2014, 2015). Some of the groups we identified might be responsible for the origin of the Belt bubble. In particular groups E and B₈ are located in the sky within the dust structure shown in Fig. 9, at different distances. Group B₈ is slightly more diffuse than the bubble, but the central over-density is still located within the bubble boundaries. The stellar winds and the supernova explosions coming from these groups might be responsible for the creation of the bubble itself.

Fig. 9. Planck data and groups E (purple dots) and B8 (black dots). The orange star represents σ Ori.

4.2.3. Orion A

The DBSCAN groups associated with the Orion A molecular cloud are those labelled B₁, C, and D. Group B₁ and C nearly occupy the same position in the sky and share very similar kinematic properties (see Table 1), however they are at different distances, with group B₁ being closer to the Sun than group C. This poses interesting questions about their origin: the two groups might be identified separately by DBSCAN just because of a local under-density of sources. In this case, the Orion A cloud would be even more elongated along the line of sight than previously thought (Großschedl et al. 2018). The radial velocities of the embedded sources in the Orion A molecular cloud are tightly related to the motion of the molecular gas in the cloud (Hacar et al. 2016). So, if the foreground is moving as the stars in the cloud, and stars in the cloud are coupled to the gas, the foreground group might have originated from the same cloud complex. The proper motion diagram of the three groups is shown in Fig. 13. We define the Orion A region as:

(15)

The proper motions of all the sources (grey dots in Fig. 13) in the region show a clump in μ_α*, μ_δ ∼ ( − 2.,1) (see also left panel of Fig. 14). We select the sources with proper motions:

(16)

(black dots in Fig. 14) and we study their distribution in the sky and on the X − Y plane in galactic Cartesian co-ordinates. We label this group as group F. Figure 14 (centre) shows that the sources are loosely distributed in the Orion A region, and seem to cluster at (l, b)∼(209, −19). Figure 14 (right) show that the members of group F are loosely spread at larger distances than the sources associated with the Orion A molecular cloud. We run the kinematic modelling on group F and we find the parameters reported in Table 1. We compare the proper motions of group F with those of the other groups, and we notice that they are roughly the same as those of group B₈ (see Fig. 10). Nevertheless the results of the kinematic modelling for the two groups are quite dissimilar. This could be due to the fact that, for both groups, the number of stars with measured radial velocity is small, and therefore the 3D velocity is not well constrained. An inspection of the parallax distribution of group F also shows a number of sources with small parallax (ϖ <  1.9 mas), which are most likely field contaminants.

Fig. 10. Distribution in the sky (top) and in 3D space (bottom) of the stars belonging to group B7, b (black dots) compared to those in group E (purple dots) The grey dots represent all the sources in the Belt region. The orange stars are the same as defined in Fig. 7.

Fig. 11. Proper motion diagram of all the sources in the Belt region. The orange rectangles are those defined in Eqs. (12) and (13).

Fig. 12. Distribution in the sky (top) and in 3D space (bottom) of the stars belonging to group B8 (black crosses) compared to those in group E (purple dots). The orange stars are the same as defined in Fig. 7.

Fig. 13. Proper motion diagram of the sources towards the Orion A molecular cloud, for group B1 (orange dots in the left panel), group C (green dots in central panel), and group D (red dots in right panel). The grey dots represent the proper motions of all the sources in the Orion A region (see text).

Fig. 14. Proper motion (left), sky distribution (centre), and distribution in the X − Y plane in Cartesian galactic co-ordinates of all the sources in the Orion A region (grey dots) and of those selected through Eq. (16) (black dots).

5. Ages

We determine ages (τ) and extinctions (A_V) of the groups we identified by performing an isochrone fit based on a maximum likelihood approach similar to the methods described in Jørgensen & Lindegren (2005), Valls-Gabaud (2014), and Zari et al. (2017).

Assuming independent Gaussian errors on all the observed quantities we can write the likelihood for a single star to come from an isochrone with certain properties θ = (τ, A_V, Z, …), as:

(17)

with:

(18)

where m is the stellar mass, n is the number of observed quantities, and q^obs and q(θ, m) are the vectors of observed and modelled quantities. To take into account the fact that stars are not distributed uniformly along the isochrone, we weight the jth likelihood with a factor w defined as:

(19)

where n_redder is the number of stars with G_BP − G_RP colour larger than that of the jth star and n_bluer is the number of stars with G_BP − G_RP smaller than that of the jth star. This choice gives larger weights to blue, massive stars, to take into account that they are fewer than the low-mass members of the clusters.

The likelihood for N coeval stars is just defined as:

(20)

Since we are interested in determining the ages and the extinctions of the groups, we fix the metallicity to Z  =  Z_⊙  =  0.0158 and we integrate Eq. (13) on the mass, so that the probability density function as a function of age τ and extinction A_V is given by:

(21)

To perform the fit we compare the observed G magnitude and G_BP − G_RP colour 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, using the passbands by Maíz Apellániz & Weiler (2018). We used isochronal tracks from log(age/yr) = 6.0 (1 Myr) to log(age/yr) = 8.0 (100 Myr), with a step of log(age/yr) = 0.05, and from A_V = 0 mag to A_V = 2.5 mag with a step of 0.1 mag.

Our fitting procedure does not take into account the presence of unresolved binaries, the photometric variability of young stars, the presence of circumstellar material, or potential age spreads within single groups. These effects can bias our age estimates and this issue is further discussed in Sect. 6.

We compute the age τ and the A_V for the groups identified by DBSCAN, and for the groups we selected in Sect. 4. The results are reported in Table 2. Figures 15 and 16 show the log-likelihood log L = log L_combined(τ, A_V) we obtain for group B₀, and the M_G vs. G − G_RP (left) and M_G vs. G_BP − G_RP (right) colour-magnitude diagrams (the colour-magnitude diagrams for the other groups are shown in Appendix B). The orange solid line corresponds to the best-fitting isochrone. As mentioned above, we perform the fit using the G_BP − G_RP colour, and we show the colour-magnitude diagram in G − G_RP as a quality check. We adopt the maximum of L_combined(τ, A_V) as our best estimate of the stellar age, and we compute the confidence intervals by evaluating the 16th and the 84th percentiles after marginalizing over A_V. Figure 15 shows a correlation between age and extinction: at large extinction values the isochrones move towards redder colours, and soon they do not intersect the upper main sequence. However they still can fit the low pre-main sequence.

Fig. 15. log L for the cluster B0. We note the correlation between age and extinction.

Fig. 16. MG vs. G − GRP (left) and MG vs. GBP − GRP (right) the colour magnitudes for group B0. The symbol sizes represent the weights assigned to each star. The solid orange line represents the best fit isochrone.

6. Dicussion

In this section we summarise and comment the results obtained in the previous sections and we put them in the broader context of the models of sequential star formation and triggering.

Kinematics. By considering the v_y, I velocities, we notice that we can roughly divide them in two groups, the first one with v_y, I ∼ 20 km s⁻¹ and the second one with v_y, I ∼ 26 − 27 km s⁻¹. We observe a loose correlation between velocity and distances (the farthest objects are also the fastest), while there is no correlation between velocity and age or distance and age.

In the kinematic modelling code we included isotropic expansion, however expansion could be an-isotropic, as observed for example by Cantat-Gaudin et al. (2019) and Wright & Mamajek (2018), although expansion due to residual gas expulsion is usually thought to be isotropic. The expansion ages determined by using the formula τ_exp = 1/(γκ) give a loose indication of the group ages, and confirm the age ordering obtained by the isochrone fitting procedure. The results of the simulations that we performed to test the kinematic modelling code (see Appendix A) showed that the expansion parameter κ always resulted to be under-estimated, thus providing over-estimated expansion ages. This is consistent with the expansion ages obtained for the DBSCAN groups.

As mentioned in Sect. 3, by using the DBSCAN algorithm we preferentially select clusters that are dense in 3D space, and tend to neglect more diffuse groups. This effect is mitigated by the visual inspection of the proper motion diagrams of the DBSCAN groups, which we use to select groups with common kinematic properties that DBSCAN fails to retrieve. Further, one of the goals of the kinematic modelling code is to exclude outliers from the DBSCAN groups. Outliers are stars that do not share the same kinematic properties as the other cluster members: this implies that also stars that should be considered cluster members, such as binaries, are excluded from the DBSCAN groups.

These considerations suggest that the groups that we analyse are not complete in terms of membership. The aim of this study is however to characterise the global properties of the stellar population in the Orion region. A more detailed analysis of the physical properties for which a complete membership list is important, such as the initial mass function, is left to future studies.

Ages. The results obtained by fitting isochrones to the colour-magnitude diagrams of the groups isolated in Sect. 4 confirm the existence of the old population towards the 25 Ori group found by Kos et al. (2018), which corresponds to our group B_0, b. Kos et al. (2018) derive an age of 20 Myr, while we obtain an age of 15 Myr. This could depend on the different extinction values used or by a slightly different membership list. We also found that, towards the Belt, group E, B_5, b, B_7, b, and B₈ are older than 10 Myr, and that some older sources are also found in the Orion A region (group F). The population in front of the Orion A cloud (group B₁) is around 10 Myr old. The age is similar to the estimated age for the group related to the Orion A cloud (group C). However, the colour-magnitude diagram of group C (see Appendix B) shows that, not unexpectedly, many sources are brighter than the 10 Myr isochrone, and therefore likely younger.

A substantial luminosity spread has been observed in the colour-magnitude diagram of the stellar population towards the ONC (see for example Jeffries et al. 2011; Da Rio et al. 2010). This spread represents the combined effect of a real age spread, possibly due to the presence of multiple populations (Jerabkova et al. 2019; Beccari et al. 2017), and of an apparent spread caused by other physical effects that scatter the measured luminosities, such as stellar variability and scattered light from circumstellar material. Age spreads are not included in our data modelling, therefore our age estimate for group C should be considered as an upper estimate for the age of the stellar population towards the Orion A molecular cloud, which also contains younger sources. The older population is more numerous than the younger ones, and therefore our age estimates are biased toward older ages. The age estimate for group C and for all the other groups is very precise (see Table 2). This is partly an artefact of using a single isochrone set, and ignoring differential extinction as well as the effects mentioned above. The presence of unresolved binaries in our data is also not taken into account, and could introduce biases towards younger age estimates, as unresolved binaries appear brighter than single stars. This could be the case for example for groups B2 and B5 (see Fig. B.1). For the other groups the single star sequence is usually more numerous than the unresolved binary sequence, thus the fit results are weighted towards the single star sequence.

In terms of age ranking, our age estimates agree with those found by Kounkel et al. (2018): their Fig. 13 indicates indeed the presence of a diffuse older population, which however they find to be around 10 Myr old. The difference in the maximum age they obtain is due to a number of differences in our fitting procedure: for example, they use A_V = 0 mag and a previous version of the Gaia DR2 filters. Our results contradict instead what was found by Briceño et al. (2019), who derive an age sequence that agrees with the long-standing picture of star formation starting in the 25 Ori region (also called Orion OB1a) and sequentially propagating towards the Belt region (Orion OB1b and 1c) and the Orion A molecular clouds (Orion OB1d).

Sequential star formation and triggering in Orion. The main result of this work is that the star formation history of the Orion region is complex and fragmentary. The Orion region is composed of many subgroups with different kinematic properties. Star formation started around 15 Myr ago (or 20 according to Kos et al. 2018), and still continues in the Orion A and B molecular clouds. The groups that we observe at the present time are sometimes spatially mixed (such as in the 25 Ori region) but their kinematics retain traces of their different origin. Figure 17 shows a schematic view of the Orion region, which summarises our results. The arrows represent the velocity vectors (in galactic Cartesian co-ordinates and corrected for the solar motion following Schönrich et al. 2010) of the groups we identified, and are colour-coded by the group ages. The grey contours represent the stellar density integrated in the Z (left), Y (centre), and X direction. The Sun is at (X, Y, Z) = (0, 0, 0) pc.

Fig. 17. Contours: density distribution of the selected sources projected in the Galactic plane (left), in the (X, Z) plane (centre), and in the (Y, Z) plane (the Sun is at (X, Y, Z) = 0, 0, 0). The arrows represent the velocities (in Galactic co-ordinates) estimated in Sect. 3 for all the groups, and are corrected for the solar motion. The colours of the arrows represent the log(age/yr) obtained by fitting the colour-magnitude diagrams of the clusters in Sect. 4.

Cantat-Gaudin et al. (2019) studied the Vela OB association, finding that a large fraction of the young stars in the region are not concentrated in clusters, but rather distributed in sparse structures, elongated along the Galactic plane. Krause et al. (2018) performed a multi-wavelength analysis of the Scorpius-Centaurus association, and suggested a refined scenario to explain the age sequence of the sub-groups that form the association. Similar to these studies, we find that the star formation history of Orion is not consistent with simple sequential star formation scenarios. Further, the traditional groups in which the Orion OB association is sub-divided are not monolithic episodes of star formation, but exhibit significant kinematic and physical sub-structure.

We do not observe any clear age gradient nor any clear evidence of triggering in the kinematic properties of the groups (such as those predicted for instance by Hartmann et al. 2001). As Cantat-Gaudin et al. (2019) suggest, the difference in velocity that are observed might be the result of galactic shear, or the consequence of a velocity pattern already imprinted in the filaments belonging to the parent molecular cloud these young populations formed from. The disposition in space of the clusters might reflect the structure of their parental molecular clouds: however this should be confirmed by specific simulations of the star formation process in the Orion region.

7. Conclusions

In this work we studied the 3D structure, the kinematics, and the age ordering of the young stellar groups of the Orion star forming region, making use of Gaia DR2.

• We selected young sources by applying simple cuts in the M_G vs. G_BP − G_RP colour-magnitude diagram, and we studied their density distribution in 3D galactic co-ordinates.

• We constructed 3D density maps normalised between 0 and 1, and we selected only the sources above a threshold of 0.5. We then applied the DBSCAN clustering algorithm to identify groups in 3D space and we analysed their properties in terms of ages and kinematics.

• We inspected the proper motions of all the groups. Certain single groups in 3D space showed sub-structures in their proper motion distribution. In this case we further sub-divided the groups, making simple cuts based on the proper motion distribution. We retrieved average motions, velocity dispersion, and isotropic expansion for all the groups identified by applying a kinematic modelling code.

• We compared the 3D velocities of all the groups, and we found evidence of kinematic sub-structures.

• We computed ages and extinctions for all the groups by using a 2D maximum likelihood approach. We found that star formation in Orion started around 15 Myr ago in two groups, one towards the Belt region, and one towards the 25 Ori region.

• We did not find any clear age gradient, or any evidence of sequential star formation propagating from the 25 Ori region towards the Belt region and the Orion A and B molecular gas.

In conclusion, the picture of the Orion that we obtained from this study is that of a highly sub-structured ensemble of young stars with different ages, with several kinematic groups, mixed in 3D space and overlapping in the sky. These results do not agree well with sequential star formation models, and would require designated specific simulations to be fully explained.

The limited number of radial velocities available for most of the groups, as well as their large uncertainties, does not allow to characterise fully the internal kinematics of the clusters, or establish the presence of an-isotropic expansion. Future, ground based spectroscopic surveys could provide precise radial velocities for a large sample of sources, which, combined with the next Gaia releases, will allow to better probe the internal kinematics of young clusters and OB associations.

Acknowledgments

We thank the referee for their comments, which improved the manuscript. 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 project was developed in part at the 2018 NYC Gaia Sprint, hosted by the Center for Computational Astrophysics at the Simons Foundation in New York City. This work has made extensive use of Matplotlib (Hunter 2007), 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.

References

Appendix A: Testing the kinematic modeling code with simulated clusters

We generate a sample of N = 200 stars which mimics the kinematics properties of young clusters and we test our code by changing (a) the position of the sample (in particular its distance to the Sun), (b) the velocity dispersion, and (c) the expansion coefficient (κ) value. In particular we are interested in the ability of the code to retrieve the correct value for κ, especially when not all the radial velocities of the cluster members are provided.

A.1. Simulation set up

The simulate star positions are drawn from Gaussian distributions with σ = 2 pc. The velocity of each simulated star is drawn following the same assumption as in L00, that is from a Gaussian distribution centred in v₀ with a small velocity dispersion σ. We include expansion following Eq. (9), chosing κ = 0.1 km s⁻¹ pc⁻¹.

We obtain the observed quantities (positions, parallax, proper motions, and radial velocities)⁴ by adding typical Gaia errors in the Orion region drawn from Gaussian distribution with widths 0.1 mas, 0.1 mas yr⁻¹, and 3 km s⁻¹ respectively.

A.2. Simple tests

We simulate two clusters at different distances and with different velocities (see Tables 1 and 2, respectively): cluster A is similar in terms of kinematics v_0, I = ( − 5.0,  45.0,  6.0) km s⁻¹ and distance (x_0, I, y_0, I, z_0, I) = (17.89,  42.14,  13.16) pc to the Hyades cluster; cluster B is instead resembling the 25 Ori cluster: (x_0, I, y_0, I, z_0, I) = (52.96,  343.97,  10.21) pc and v_0, I = (0.0,  20.0,  0.0) km s⁻¹. We run the simulations in five different scenarios for both the simulated clusters:

1. σ_v = 0.3 km s⁻¹ and κ = 0.1 km s⁻¹ pc⁻¹.

2. σ_v = 1.0 km s⁻¹ and κ = 0.1 km s⁻¹ pc⁻¹;

3. σ_v = 0.3 km s⁻¹, κ = 0.1 km s⁻¹ pc⁻¹, and a fraction f = [10%,50%,95%] of stars without measured radial velocities.

The average velocities are always retrieved quite correctly in both cases; σ and κ are retrieved correctly for cluster A, however we notice that for cluster B the value of κ is usually underestimated, while σ is usually slightly over-estimated. When the number of observed radial velocities is too low, the expansion parameter can not be retrieved as it can not be separated from v₀ from astrometric data only. In the cases when this happens, we do not give any estimate for the expansion term κ. When there are no radial velocities available the velocity is very poorly constrained, especially for cluster B: in this case we do not give estimates for the velocities. When 10% or 50% of the measured radial velocities are missing, the errors on the estimated parameters are of the same order of magnitude as in the other cases were all the kinematic data are available. However, not unexpectedly, when only 5% of the radial velocities is available, the error on the v_y parameter is roughly one order of magnitude larger than in the other cases.

A.3. Realistic tests

In the real case it is likely that the clusters selected with the DBSCAN algorithm have both stars without measured radial velocities and kinematic outliers. We therefore further tested our code for cluster in two cases (see Table A.1). In the first one we include 20 kinematic outliers in our simulated clusters: the kinematic outliers have a broader spatial distribution than the simulated cluster members (σ = 5 pc), and their velocities are drawn from a Gaussian distribution with mean 20 km s⁻¹ in x_I, y_I, z_I, and dispersion σ_v = 10 km s⁻¹. In the second one we include 20 kinematic outliers and we remove the 10% of measured radial velocities. In both cases, after the exclusion procedure the parameters are retrieved correctly. We notice that also in this case the expansion coefficient κ is under-estimated (roughly by a factor of 2), while σ_v is slightly over-estimated.

A.4. Initial conditions

To test whether the initial conditions of the minimisation have an impact on the estimated parameters, we performed 100 runs with initial guesses for the mean cluster velocity components, the velocity dispersion, and the expansion term κ drawn randomly from a Gaussian distribution centred on the mean parameters, with dispersion equal to the 20% of their real values. Reino et al. (2018) performed similar tests on the Hyades cluster (which as said above is kinematically similar to our cluster A), finding essentially no dependence from the estimated parameters from the initial conditions. Thus, we repeat these tests only on our simulated cluster B.

B.1: σ_v = 0.3 km s⁻¹ and κ = 0.1 km s⁻¹ pc⁻¹. We find that in general the minimisation results do not strongly depend on the initial parameters, however if the velocity dispersion σ_v is over-estimated and (or) the velocity in the x_I component is under- or over-estimated then the velocity in the y_I component is also under- or over-estimated.

B.2: σ_v = 1.0 km s⁻¹ and κ = 0.1 km s⁻¹ pc⁻¹. We find that the minimisation results do not depend on the initial parameters in any case. This is reassuring, as the values for σ_v in the clusters considered here are larger than 0.3 km s⁻¹. In the cases with σ_v = 1.0 km s⁻¹ and missing radial velocities (for 20, 100, and 190 stars respectively), the estimated parameters are retrieved correctly for any choice of initial conditions, except for the expansion parameter κ, that is underestimated. If outliers are present, the parameters are retrieved correctly after the exclusion procedure.

Appendix B: Colour magnitude diagrams

Figure B.1 shows the colour magnitude diagram for the groups that we identified in Sect. 4.

Fig. B.1. MG vs. GBP − GRP colour magnitude diagram for the groups selected in Sect. 4. The blue solid lines correspond to the best fitting isochrones, derived in Sect. 5.

All Tables

All Figures

	Fig. 1. Observed colour-magnitude diagram (left) and sky distribution (right) of the sources selected in Sect. 2. The solid black lines in the left panel isochrones defined in Eq. (6), which are used to select the young stellar population in Orion.
In the text

	Fig. 2. Density distribution of the sources in Orion for different Z values. The orange stars indicate the positions of (from top to bottom in each panel): λ Ori, 25 Ori, ϵ Ori, and the Orion Nebula cluster. The white solid contours represent the 0.4, 0.5, 0.6 and 0.8 density levels (the density is normalised to have values from 0 to 1). The Sun is located at (X, Y) = (0, 0).
In the text

	Fig. 3. Distribution of the groups identified with the first (top) and second (bottom) iteration of DBSCAN in the planes (X, Y), (X, Z), and (Y, Z).
In the text

	Fig. 4. Sky distribution of the groups identified with the first (top) and second (bottom) iteration of DBSCAN. The colours correspond to those in Fig. 3.
In the text

	Fig. 5. Left: proper motion diagram of the stars in the 25 Ori region defined in the text (black dots), of the stars belonging to group B0 (blue dots), and of the stars belonging to group B0, b (light blue empty squares) and B6 (pink crosses). Right: sky distribution in galactic co-ordinates of group B0, B0, b, and B6. The colours and symbols are the same as on the left.
In the text

	Fig. 6. Proper motion diagram for all the stars in the Belt region (grey dots). Each panel corresponds to one of the groups identified by DBSCAN (the colours are the same as in Figs. 3 and 4) except for group B7, which is indicated by black dots for representation purposes.
In the text

	Fig. 7. Distribution in the sky of the stars in the Belt region. Each panel corresponds to one group identified by DBSCAN. The colour-coding is the same as in Fig. 6. Orange stars correspond to, from left to right: σ Ori, ζ Ori, ϵ Ori, η Ori, δ Ori.
In the text

Fig. 8. Top panel: sky distribution of the sources belonging to group B5, b (brown empty squares), group B5 (brown dots), and all the sources in the Belt region defined in the text (grey dots). The orange stars mark the position of σ Ori, ζ Ori, ϵ Ori, η Ori, and δ Ori. Bottom panels: distribution in 3D galactic co-ordinates of group B5, b (brown squares), and of all the sources belonging to the Belt region (grey dots).

In the text

	Fig. 9. Planck data and groups E (purple dots) and B8 (black dots). The orange star represents σ Ori.
In the text

	Fig. 10. Distribution in the sky (top) and in 3D space (bottom) of the stars belonging to group B7, b (black dots) compared to those in group E (purple dots) The grey dots represent all the sources in the Belt region. The orange stars are the same as defined in Fig. 7.
In the text

	Fig. 11. Proper motion diagram of all the sources in the Belt region. The orange rectangles are those defined in Eqs. (12) and (13).
In the text

	Fig. 12. Distribution in the sky (top) and in 3D space (bottom) of the stars belonging to group B8 (black crosses) compared to those in group E (purple dots). The orange stars are the same as defined in Fig. 7.
In the text

	Fig. 13. Proper motion diagram of the sources towards the Orion A molecular cloud, for group B1 (orange dots in the left panel), group C (green dots in central panel), and group D (red dots in right panel). The grey dots represent the proper motions of all the sources in the Orion A region (see text).
In the text

	Fig. 14. Proper motion (left), sky distribution (centre), and distribution in the X − Y plane in Cartesian galactic co-ordinates of all the sources in the Orion A region (grey dots) and of those selected through Eq. (16) (black dots).
In the text

	Fig. 15. log L for the cluster B0. We note the correlation between age and extinction.
In the text

	Fig. 16. MG vs. G − GRP (left) and MG vs. GBP − GRP (right) the colour magnitudes for group B0. The symbol sizes represent the weights assigned to each star. The solid orange line represents the best fit isochrone.
In the text

Fig. 17. Contours: density distribution of the selected sources projected in the Galactic plane (left), in the (X, Z) plane (centre), and in the (Y, Z) plane (the Sun is at (X, Y, Z) = 0, 0, 0). The arrows represent the velocities (in Galactic co-ordinates) estimated in Sect. 3 for all the groups, and are corrected for the solar motion. The colours of the arrows represent the log(age/yr) obtained by fitting the colour-magnitude diagrams of the clusters in Sect. 4.

In the text

	Fig. B.1. MG vs. GBP − GRP colour magnitude diagram for the groups selected in Sect. 4. The blue solid lines correspond to the best fitting isochrones, derived in Sect. 5.
In the text