Structure, kinematics, and ages of the young stellar populations in the Orion region

TL;DR

This study uses Gaia DR2 data to analyze the 3D structure, kinematics, and ages of stellar populations in Orion, revealing complex, multi-episode star formation history with multiple sub-structures and ages.

Contribution

It provides a detailed 3D map and kinematic analysis of Orion's stellar groups, highlighting multiple star formation events and challenging simple sequential models.

Findings

Presence of an old (~15 Myr) population near 25 Ori

Detection of 12-15 Myr groups in the Belt region

Identification of a ~10 Myr population in front of Orion A cloud

Abstract

We present a study of the three dimensional structure, kinematics, and age distribution of the Orion OB association, based on the second data release of the Gaia satellite (Gaia DR2). Our goal is to obtain a complete picture of the star formation history of the Orion complex and to relate our findings to theories of sequential and triggered star formation. We select the Orion population with simple photometric criteria, and we construct a three dimensional map in galactic Cartesian coordinates to study the physical arrangement of the stellar clusters in the Orion region. The map shows structures that extend for roughly $150 \, \mathrm{pc}$ along the line of sight, divided in multiple sub-clusters. We separate different groups by using the density based clustering algorithm DBSCAN. We study the kinematic properties of all the groups found by DBSCAN first by inspecting their proper motion distribution, and then by applying a kinematic modelling code based on an iterative maximum likelihood approach, which we use to derive their mean velocity, velocity dispersion and isotropic expansion. By using an isochrone fitting procedure we provide ages and extinction values for all the groups. We confirm the presence of an old population ($\sim 15$ Myr) towards the 25 Ori region, and we find that groups with ages of $12-15 \, \mathrm{Myr}$ are present also towards the Belt region. We notice the presence of a population of $\sim 10$ Myr also in front of the Orion A molecular cloud. Our findings suggest that star formation in Orion does not follow a simple sequential scenario, but instead consists of multiple events, which caused kinematic and physical sub-structure. To fully explain the detailed sequence of events, specific simulations and further radial velocity data are needed.

Tables

Table 1. Table 1: Results of the kinematic modelling for the groups identified by DBSCAN. The first column indicates the region in the sky where the groups are located.

	#	$N$	$N_{RV}$	$v_{x,I}\left[\mathrm{km}{s}^{-1}\right]$	$v_{y,I}\left[\mathrm{km}{s}^{-1}\right]$	$v_{z,I}\left[\mathrm{km}{s}^{-1}\right]$	${\sigma }_v\left[\mathrm{km}{s}^{-1}\right]$	$\kappa \left[\mathrm{km}{s}^{-1}{\mathrm{pc}}^{-1}\right]$	$\varpi \left[\mathrm{mas}\right]$
$\lambda$ Ori	A	296	81	0.75 $\pm$ 0.05	27.4 $\pm$ 0.1	0.5 $\pm$ 0.05	0.73 $\pm$ 0.02	0.122 $\pm$ 0.007	${2.48}_{0.06}^{0.07}$
Orion A	C	1059	489	0.94 $\pm$ 0.06	27.1 $\pm$ 0.09	-2.68 $\pm$ 0.06	1.63 $\pm$ 0.03	0.07 $\pm$ 0.006	${2.5}_{0.08}^{0.06}$
Orion A	D	69	50	1.1 $\pm$ 0.1	21.8 $\pm$ 0.2	-3.7 $\pm$ 0.1	0.98 $\pm$ 0.07	-0.02 $\pm$ 0.02	${2.37}_{0.03}^{0.02}$
Belt	E	150	10	4.6 $\pm$ 0.3	32.4 $\pm$ 2.7	-1 $\pm$ 0.1	1.28 $\pm$ 0.06	-	${2.31}_{0.04}^{0.03}$
25 Ori	$B_0$	710	73	0.41 $\pm$ 0.04	19.2 $\pm$ 0.2	0.045 $\pm$ 0.03	0.74 $\pm$ 0.02	-	${2.84}_{0.08}^{0.08}$
25 Ori	$B_{0b}$	54	4	4.8 $\pm$ 0.4	25. $\pm$ 2.8	2.1 $\pm$ 0.1	0.38 $\pm$ 0.04	-	${2.78}_{0.06}^{0.06}$
Orion A	$B_1$	265	73	0.8 $\pm$ 0.1	26.8 $\pm$ 0.2	-3.0 $\pm$ 0.1	1.55 $\pm$ 0.05	0.03 $\pm$ 0.02	${2.74}_{0.07}^{0.04}$
Belt	$B_2$	174	48	0.17 $\pm$ 0.14	27.84 $\pm$ 0.3	-2.6 $\pm$ 0.1	1.6 $\pm$ 0.07	0.08 $\pm$ 0.03	${2.48}_{0.05}^{0.04}$
Belt	$B_3$	290	48	-0.74 $\pm$ 0.03	22.5 $\pm$ 0.2	-2.69 $\pm$ 0.05	0.79 $\pm$ 0.03	-	${2.78}_{0.08}^{0.07}$
Belt	$B_4$	46	10	-0.2 $\pm$ 0.3	26.5 $\pm$ 0.6	-3.1 $\pm$ 0.2	1.6 $\pm$ 0.1	0.07 $\pm$ 0.06	${2.62}_{0.02}^{0.02}$
Belt	$B_5$	248	12	0.9 $\pm$ 0.1	20.7 $\pm$ 0.9	-2.44 $\pm$ 0.06	0.7 $\pm$ 0.03	-	${2.79}_{0.06}^{0.09}$
Belt	$B_{5,b}$	622	48	1.2 $\pm$ 0.1	24.8 $\pm$ 0.2	-1.42 $\pm$ 0.04	0.82 $\pm$ 0.02	-	${2.75}_{0.36}^{0.22}$
25 Ori	$B_6$	40	0	-	-	-	-	-	-
Belt	$B_7$	30	0	-	-	-	-	-	-
Belt	$B_{7,b}$	441	63	4.9 $\pm$ 0.05	26.7 $\pm$ 0.2	-1.55 $\pm$ 0.05	0.95 $\pm$ 0.03	0.024 $\pm$ 0.003	${2.36}_{0.23}^{0.2}$
Belt	$B_8$	245	18	5.8 $\pm$ 0.08	28.5 $\pm$ 0.4	1.5 $\pm$ 0.06	0.9 $\pm$ 0.03	0.05 $\pm$ 0.05	${2.34}_{0.18}^{0.2}$
Orion A	F	116	17	5.8 $\pm$ 0.1	21.2 $\pm$ 0.5	0.3 $\pm$ 0.12	1. $\pm$ 0.06	-	${2.27}_{0.48}^{0.2}$

Table 2. Table 2: Age estimates for the groups identified in Section 3 and 4. The column log ⁡ ( age / yr ) age yr \mathrm{\log(age/yr)} ( τ 𝜏 \tau ) indicates the (log-)age estimated by the isochrone fitting procedure. The column τ e ​ x ​ p subscript 𝜏 𝑒 𝑥 𝑝 \tau_{exp} indicates the expansion ages determined by using the formula τ e ​ x ​ p = 1 / ( γ ​ κ ) subscript 𝜏 𝑒 𝑥 𝑝 1 𝛾 𝜅 \tau_{exp}=1/(\gamma\,\kappa) for the groups for which it is possible to determine the expansion parameter κ 𝜅 \kappa . The number of stars N 𝑁 N is different than in Table 3 because by applying the kinematic modelling we remove kinematic outliers from the groups.

#	$N$	$\log \left(\mathrm{age}/\mathrm{yr}\right)$	$\tau \left[\mathrm{Myr}\right]$	$A_V\left[\mathrm{mag}\right]$	${\tau }_{exp}\left[\mathrm{Myr}\right]$
A	274	${6.75}_{0.01}^{0.03}$	${5.6}_{0.1}^{0.4}$	0.4	8.0
C	943	${6.9}_{0.01}^{0.03}$	$8_{0.04}^{0.5}$	0.2	14.0
D	60	${6.85}_{0.02}^{0.03}$	$7_{0.2}^{0.6}$	1.3	-
E	139	${7.05}_{0.005}^{0.04}$	${11.2}_{0.1}^1$	0.5	-
$B_0$	622	${7.05}_{0.005}^{0.04}$	${11.2}_{0.1}^1$	0.2	-
$B_{0b}$	44	${7.15}_{0.004}^{0.1}$	${14}_{0.25}^3$	0.4	-
$B_1$	246	${7.0}_{0.01}^{0.03}$	${10}_{0.23}^{0.7}$	0.4	32.6
$B_2$	154	${6.6}_{0.01}^{0.03}$	$4_{0.1}^{0.3}$	0.3	12.2
$B_3$	221	${6.9}_{0.01}^{0.04}$	$8_{0.04}^{0.7}$	0.2	-
$B_4$	44	${6.6}_{0.01}^{0.03}$	$4_{0.1}^{0.3}$	0	14
$B_5$	234	${6.9}_{0.01}^{0.04}$	$8_{0.04}^{0.7}$	0.2	-
$B_{5,b}$	605	${7.05}_{0.005}^{0.03}$	${11.2}_{0.1}^1$	0.2	-
$B_{7,b}$	418	${7.05}_{0.005}^{0.04}$	${11.2}_{0.1}^1$	0.3	40
$B_8$	237	${7.15}_{0.004}^{0.04}$	${14}_{0.25}^{1.5}$	0.3	-
F	108	${7.05}_{0.005}^{0.03}$	${11.2}_{0.1}^1$	0.3	-

Table 3. Table 3: Results of the tests of the kinematic modelling for cluster A.

	$v_{x,I}\left[\mathrm{km}{s}^{-1}\right]$	$v_{y,I}\left[\mathrm{km}{s}^{-1}\right]$	$v_{z,I}\left[\mathrm{km}{s}^{-1}\right]$	${\sigma }_v\left[\mathrm{km}{s}^{-1}\right]$	$\kappa \left[\mathrm{km}{s}^{-1}{\mathrm{pc}}^{-1}\right]$
Initial values	-5.0	45.0	6.0	0.3	0.1
A.1	-5.9	45.6	5.57	0.3	0.1
	-5.88 $\pm$ 0.03	45.57 $\pm$ 0.05	5.56 $\pm$ 0.027	0.32 $\pm$ 0.01	0.1 $\pm$ 0.01
A.2	-5.9	45.6	5.57	1.0	0.1
	-5.96 $\pm$ 0.08	45.5 $\pm$ 0.1	5.54 $\pm$ 0.08	1.01 $\pm$ 0.03	0.06 $\pm$ 0.02
A.3
20/200 missing radial velocities	-5.873 $\pm$ 0.03	45.6 $\pm$ 0.06	5.586 $\pm$ 0.03	0.31 $\pm$ 0.01	0.1 $\pm$ 0.01
100/200 missing radial velocities	-5.91 $\pm$ 0.035	45.55 $\pm$ 0.07	5.564 $\pm$ 0.03	0.3 $\pm$ 0.01	0.1 $\pm$ 0.1
190/200 missing radial velocities	-6.0 $\pm$ 0.4	45.0 $\pm$ 1.0	5.5 $\pm$ 0.3	0.3 $\pm$ 0.01	0.1 $\pm$ 0.03

Table 4. Table 4: Results of the tests of the kinematic modelling for cluster B.

	$v_{x,I}\left[\mathrm{km}{s}^{-1}\right]$	$v_{y,I}\left[\mathrm{km}{s}^{-1}\right]$	$v_{z,I}\left[\mathrm{km}{s}^{-1}\right]$	${\sigma }_v\left[\mathrm{km}{s}^{-1}\right]$	$\kappa \left[\mathrm{km}{s}^{-1}{\mathrm{pc}}^{-1}\right]$
Initial values	0.0	20.	0.0	0.3	0.1
B.1	0.66	19.73	0.53	0.3	0.1
	0.65 $\pm$ 0.03	19.83 $\pm$ 0.063	0.51 $\pm$ 0.026	0.36 $\pm$ 0.013	0.07 $\pm$ 0.01
B.2	0.66	19.73	0.53	1.0	0.1
	0.65 $\pm$ 0.07	19.8 $\pm$ 0.1	0.56 $\pm$ 0.07	1.0 $\pm$ 0.03	0.05 $\pm$ 0.01
B.3
20/200 missing radial velocities	0.69 $\pm$ 0.03	19.80 $\pm$ 0.06	0.57 $\pm$ 0.02	0.33 $\pm$ 0.01	0.05 $\pm$ 0.006
100/200 missing radial velocities	0.67 $\pm$ 0.03	19.87 $\pm$ 0.08	0.52 $\pm$ 0.03	0.38 $\pm$ 0.01	0.05 $\pm$ 0.007
190/200 missing radial velocities	0.67 $\pm$ 0.043	19.89 $\pm$ 0.232	0.56 $\pm$ 0.025	0.32 $\pm$ 0.01	0.09 $\pm$ 0.008

Table 5. Table 5: Results of the tests of the kinematic modelling for cluster B, with missing radial velocities and outliers.

	$v_{x,I}\left[\mathrm{km}{s}^{-1}\right]$	$v_{y,I}\left[\mathrm{km}{s}^{-1}\right]$	$v_{z,I}\left[\mathrm{km}{s}^{-1}\right]$	${\sigma }_v\left[\mathrm{km}{s}^{-1}\right]$	$\kappa \left[\mathrm{km}{s}^{-1}{\mathrm{pc}}^{-1}\right]$
No missing radial velocities and 20 outliers
Initial parameters	0.0	20.	0.0	0.3	0.1
Real values	0.66	19.73	0.53	0.3	0.1
First iteration	0.62 $\pm$ 6.	19.88 $\pm$ 6.	0.71 $\pm$ 6	88.09 $\pm$ 2.4	-0.17 $\pm$ 0.4
After exclusion procedure	0.63 $\pm$ 0.03	19.75 $\pm$ 0.06	0.52 $\pm$ 0.03	0.33 $\pm$ 0.01	0.05 $\pm$ 0.006
20 missing radial velocities and 20 outliers
Initial parameters	0.0	20.	0.0	0.3	0.1
Real values	0.66	19.73	0.53	0.3	0.1
First iteration	0.64 $\pm$ 4.85	19.88 $\pm$ 5.2	0.6 $\pm$ 4.85	72. $\pm$ 2.	0.4 $\pm$ 0.3
After exclusion procedure	0.68 $\pm$ 0.03	19.83 $\pm$ 0.06	0.55 $\pm$ 0.03	0.37 $\pm$ 0.01	0.07 $\pm$ 0.007

Equations

$$190^{\circ }<l<220^{\circ },-30^{\circ }<b<-5^{\circ },$$

$$|{\mu }_{\alpha \ast }|<10\mathrm{mas}\mathrm{y}{\mathrm{r}}^{-1},|{\mu }_{\delta }|<10\mathrm{mas}\mathrm{y}{\mathrm{r}}^{-1}.$$

$$\mathrm{RUWE}=\sqrt{{\chi }^2/(N-5)}/u_0(C,G)$$

$$E=(I_{\mathrm{BP}}+I_{\mathrm{RP}})/I_G,$$

$$1.0+0.015(G_{\mathrm{BP}}-G_{\mathrm{RP}}{)}^2<E<1.3+0.06(G_{\mathrm{BP}}-G_{\mathrm{RP}}{)}^2.$$

$$G_{BP}-G_{\mathrm{RP}}-0.2$$

$$G_{BP}-G_{\mathrm{RP}}+0.5$$

$$f(D)=\frac{L}{1+e^{-\kappa (D-D_0)}}-1,$$

$${\mathbf{u}}_i={\mathbf{v}}_0+\mathbf{T}({\mathbf{b}}_i-{\mathbf{b}}_0),$$

$$\mathbf{T}=\left[\begin{array}{ccc}\kappa & 0& 0\\ 0& \kappa & 0\\ 0& 0& \kappa \end{array}\right]$$

$$199^{\circ }<l<203^{\circ }-20^{\circ }<b<-15^{\circ },$$

$${\mu }_{\alpha \ast }<0\mathrm{mas}\mathrm{y}{\mathrm{r}}^{-1}{\mu }_{\delta }>-1\mathrm{mas}\mathrm{y}{\mathrm{r}}^{-1}.$$

$$203^{\circ }<l<207.5^{\circ }-21^{\circ }<b<-13^{\circ }.$$

$$0.3\mathrm{mas}\mathrm{y}{\mathrm{r}}^{-1}$$

$$-0.8\mathrm{mas}\mathrm{y}{\mathrm{r}}^{-1}$$

$$-2.2\mathrm{mas}\mathrm{y}{\mathrm{r}}^{-1}$$

$$0.4\mathrm{mas}\mathrm{y}{\mathrm{r}}^{-1}$$

$$-2.2\mathrm{mas}\mathrm{y}{\mathrm{r}}^{-1}$$

$$-2\mathrm{mas}\mathrm{y}{\mathrm{r}}^{-1}$$

$$207.5^{\circ }<l<216^{\circ },-22^{\circ }<b<-17^{\circ }.$$

$$-2.5\mathrm{mas}\mathrm{y}{\mathrm{r}}^{-1}$$

$$0\mathrm{mas}\mathrm{y}{\mathrm{r}}^{-1}$$

$$L(\mathbf{\theta },m)=\prod _{i=1}^n\left(\frac{1}{(2\pi {)}^{1/2}{\sigma }_i}\right)\times \exp \left(-{\chi }^2/2\right),$$

$${\chi }^2=\sum _{i=1}^n{\left(\frac{q_i^{\mathrm{obs}}-q_i(\mathbf{\theta },m)}{{\sigma }_i}\right)}^2,$$

$$w=\sqrt{\frac{n_{redderj}}{n_{bluerj}+1}},$$

$$L_{combined}(\mathbf{\theta },m)=\prod _{j=1}^N{L}_j(\mathbf{\theta },m{)}^{w_j}$$

$$L_{combined}(\tau ,A_V)=\prod _{j=1}^N\int L_j(\tau ,A_V,m)dm$$

Full text

1 Introduction

The tendency of O and B type stars to loosely cluster in the sky was recognized at the beginning of the 20th century by the pioneering studies summarized 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. Due to the fact that their members are widely dispersed over the sky, OB associations have been long considered as expanding remnants of young star clusters (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 and 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, as, while they are slowly dispersing in the field, they still keep memory of the parental gas sub-structure where they originated.

At a distance of $\sim 380\,\mathrm{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 (e.g., Brown et al., 1994; Jeffries et al., 2006; Bally, 2008; Briceno, 2008; Muench et al., 2008; Da Rio et al., 2014; Getman et al., 2014; Da Rio et al., 2016; 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 et al., 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 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\sim 380\,\mathrm{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 et al., 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, and find that one cluster is significantly older ($21\pm 2\,\mathrm{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 $\sim 75\,\mathrm{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 Section 2 we present the data and describe how we select the young stellar population in Orion. In Section 3 we study its 3D configuration in Cartesian galactic coordinates, and we isolate young groups by making use of the DBSCAN clustering algorithm. In Section 4 we perform the kinematic analysis by using a maximum likelihood approach. In Section 5 we derive ages and extinctions of all the groups resulting from the analysis of Section 4. In Section 6 we discuss our findings. The conclusions of this work are summarised in Section 7.

2 Data

Following Zari et al. (2017), we select the sources with coordinates

[TABLE]

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

[TABLE]

We derive distances by inverting parallaxes, $d=1000/\varpi\,\mathrm{pc}$ thus we restrict our sample to sources with $\varpi/\sigma_{\varpi}>5.0$, following the recommendations in Bailer-Jones (2015). The effect of this cut is to exclude sources at faint magnitudes ($G>20\,\mathrm{mag}$), but it does not introduce significant biases in e.g., 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 (https://www.cosmos.esa.int/web/gaia/dr2-known-issues):

• •

The Renormalized Unit Weight Error (RUWE) is defined as:

[TABLE]

where:

• –

$\chi^{2}$ 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_{0}(C,G)$ is an empirical normalization factor, which is a function of $C=G_{\mathrm{BP}}-G_{\mathrm{RP}}$ and $G$.

and we select all the sources with RUWE ¡ 1.4, following the slides by Lindegren et al. (see https://www.cosmos.esa.int/web/gaia/dr2-known-issues). This cut seeks to remove sources with spurious parallaxes or proper motions.

• •

We use the flux excess ratio:

[TABLE]

where $I_{\mathrm{X}}$ is the photometric flux in band $X$, to exclude sources with possible issues in the $\mathrm{BP}$ and $\mathrm{\mathrm{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:

[TABLE]

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\leq G\leq 6$ mag:

$G_{\mathrm{corr}}=-0.047344+1.16405G-0.046799G^{2}+0.0035015G^{3}$

• •

$2\leq G\leq 4$ mag:

$G_{\mathrm{BP,corr}}=G_{\mathrm{BP}}-2.0384+0.95282G-0.11018G^{2}$

• •

$2\leq G\leq 3.5$ mag:

$G_{\mathrm{\mathrm{RP},corr}}=G_{\mathrm{\mathrm{RP}}}-13.946+14.239G_{\mathrm{RP}}-4.23G_{\mathrm{RP}}^{2}+0.4532G_{\mathrm{RP}}^{3}$

• •

$6\leq G\leq 16$ mag:

$G_{\mathrm{corr}}=G-0.0032(G-6.0)$

• •

$G>16$ mag:

$G_{\mathrm{corr}}=G-0.032$

In the rest of the paper we use the corrected $G$, $G_{\mathrm{BP}}$, and $G_{\mathrm{RP}}$ magnitudes without using the subscript ”corr”.

2.2 Selecting the young stellar popoulation

Figure 1 (left) shows the $M_{G}$ vs. $G_{BP}-G_{\mathrm{RP}}$ colour-magnitude diagram of the ’clean’ sample obtained in Section 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.

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\,\mathrm{mag}$ and age $\tau=10\,\mathrm{Myr}$ to define the following region in the $M_{G}$ vs. $G_{BP}-G_{\mathrm{RP}}$ colour-magnitude diagram (solid black lines in Fig. 1):

[TABLE]

We choose $A_{V}=0.3\,\mathrm{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 Section 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 study the three-dimensional (3D) distribution of sources using a similar approach as in Zari et al. (2018). In summary, we:

compute galactic Cartesian coordinates for all the sources, $x_{g},y_{g},z_{g}$; 2. 2.

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

compute the number of sources in each cube; 4. 4.

compute the source density $D(x_{g},y_{g},z_{g})$ by smoothing the distribution with a Gaussian filter, with width $w=2\,\mathrm{pc}$; 5. 5.

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

[TABLE]

with $L=2$, $\kappa=4\,\mathrm{pc}$, and $D_{0}=0$.

Fig. 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\,\mathrm{pc}$ and $Z=-101\,\mathrm{pc}$. The last three panels mainly show the $\lambda$ Ori cluster. At $Z=-92\,\mathrm{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 $\sigma_{\varpi}=0.046\,\mathrm{mas}$, but presents a long tail towards larger values (the $84th$ percentile is $0.11\,\mathrm{mas}$).

To isolate the members of each cluster, we first consider 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 are then separated by using the DBSCAN algorithm 111We use the scikit-learn implementation of the algorithm (Pedregosa et al., 2011). As described e.g. 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 need to apply the clustering algorithm twice. The first time we use 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 paramenters: 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).

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 Section 6.

4 Kinematics

In this section we study the kinematics of the groups selected in the previous section. We use 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 use 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 Section 4.1, tested in Appendix A, and the results are presented in Section 4.2. Note that here we use ICRS coordinates, which we differentiate from galactic coordinates 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 $\chi^{2}$ statistic used in L00, and denoted $g$, to a $p$ value or $1-\mathrm{CDF}($g$,\mathrm{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 $\bm{b}_{i}$ as:

[TABLE]

where $\bm{b}_{0}$ is an arbitrary reference position, namely the point where the local velocity $\bm{u}(\bm{b})$ assumes the status of ’centroid’ velocity $\bm{v}_{0}$. The coordinates of $\bm{b}_{0}$ are therefore fixed in advance. The matrix $\bm{T}$ is simply a diagonal matrix of the form:

[TABLE]

An expanding cluster will have $\kappa>0$, from which an expansion age, $\tau=1/(\gamma\kappa)\,\mathrm{[Myr]}$ can be derived ($\gamma$ is a conversion factor of 1.0227 $\mathrm{pc\,Myr^{-1}\,km^{-1}\,s,}$ see e.g. Wright & Mamajek, 2018).

The method is applied to the members of the clusters identified in Section 3. These clusters still contain ’outliers’, i.e. real non-members, or members which have (slightly) discrepant astrometry (and/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}\leq g_{lim}$ or $p_{i}\geq p_{lim}$). The stopping criterion is the same as in Reino et al. (2018), i.e. that 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 $\bm{v}_{0}$ (see L00). Therefore when the fraction of measured radial velocities is lower than the 20% we do not estimate the expansion coefficient $\kappa$ (implicitly assuming $\kappa=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 give in Table 1. Being quite isolated with respect to the rest of the population, the $\lambda$ 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^{\circ}<l<216^{\circ}$. 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:

[TABLE]

which corresponds to the groups $B_{0}$ and $B_{6}$ 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_{0}$ (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).

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:

[TABLE]

We applied separately the kinematic modelling code to the two groups. The results are reported in Table 1. Note that 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_{0}$ and $B_{0,b}$ is shown in Fig. 5 (right panel). While group $B_{0}$’s distribution shows a clump towards 25 Ori, and the Northern elongation reported by e.g. 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 Section 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_{6}$ 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_{6}$ is closer to the Sun than group $B_{0}$, while the proper motion distribution does not show any difference with respect to group $B_{0}$. We suspect that group $B_{6}$ coincides with a small over-density of sources within group $B_{0}$, which gets classified as a separate group because of a local density drop. We ran the kinematic modelling code for groups $B_{0}$ and $B_{6}$ together: the estimated parameters are consistent with those found for group $B_{0}$ only, which supports our hypothesis.

4.2.2 Belt

Many of the clusters identified by DBSCAN ($B_{2},B_{3},B_{4},B_{5},B_{7}$ and E) are located in the Sky towards the Belt region. Fig. 6 shows the proper motion diagram for the Belt region defined as

[TABLE]

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_{2}$ and $B_{4}$ are mostly located towards the $\sigma$ Ori cluster (see Fig. 7) and $\zeta$ Ori. Group $B_{3}$’s members are spread towards $\epsilon$ Ori and $\delta$ Ori. The parameters estimated by the kinematic modelling code suggest that $B_{2}$ and $B_{4}$ have compatible $v_{y,I}$ values, which are significantly different from those of group $B_{3}$. 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_{3}$ are similar to those of groups D, $B_{0}$ (not located in the Belt region, see Fig. 4), and $B_{5}$. We notice that group $B_{2}$’s velocity dispersion is large ($\sim 1.6\,\mathrm{km\,s^{-1}}$) compared for instance to that of group $B_{3}$ ($0.41\pm 0.02\,\mathrm{km\,s^{-1}}$). 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_{2}$’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.

• •

Group $B_{5}$ is located below the Belt, towards $\eta$ Ori, and shares similar kinematics with group $B_{3}$, although they seem to be well separated in space (see Fig. 3 and 4). The proper motion distribution shows two clumps, similar to what is observed towards 25 Ori. We separate the the smaller clump, which we refer to as $B_{5,b}$ by using simple cuts in proper motion space:

[TABLE]

In contrast to what we have done for group $B_{0,b}$, here we apply the conditions of Eq. • ‣ 4.2.2 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_{5}$ and group $B_{5,b}$ have different kinematic properties, while having similar parallaxes. Comparing Fig. 6 and Fig. 10 one can notice that the region defined in Eq. • ‣ 4.2.2 also includes sources classified as members of group $B_{2}$. The sky distribution of sources belonging to group $B_{5,b}$ (see Fig. 8) shows indeed some sources clustering around $\sigma$ Ori. Most of the sources however are located in the same region as group $B_{5}$, 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_{2}$ and $B_{4}$ ($\sigma$ Ori), while others are located closer to group $B_{5}$. This explains why DBSCAN does not separate successfully groups $B_{5}$ and $B_{5,b}$: their members show different kinematics but are mixed in space.

• •

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_{1}$, $B_{2}$, $B_{4}$ (and $B_{7,b}$ and $B_{8}$, see below), the $v_{x,I}$ component is different from the other groups. As for group $B_{5}$, 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_{8}$ by applying the following conditions:

[TABLE]

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_{7}$, 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_{8}$, 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_{7b}$ according to the following equations (see Fig. ):

[TABLE]

The number of sources is now much larger (see Table 1), and the parameters can be accurately determined. Fig. 10 shows the source distribution in the sky and in Cartesian galactic coordinates. We notice that the sources are distributed in the sky towards the reflection nebulae M78 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_{8}$ are located in the sky within the dust structure shown in Fig. 9, at different distances. Group $B_{8}$ 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.

4.2.3 Orion A

The DBSCAN groups associated with the Orion A molecular cloud are those labelled $B_{1}$, C, and D. Group $B_{1}$ 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_{1}$ 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:

[TABLE]

The proper motions of all the sources (grey dots in Fig. 13) in the region show a clump in $\mu_{\alpha*},\mu_{\delta}\sim(-2.,1)$ (see also left panel of Fig. 14). We select the sources with proper motions:

[TABLE]

(black dots in Fig. 14) and we study their distribution in the sky and on the $X-Y$ plane in galactic Cartesian coordinates. We label this group as group F. Fig. 14 (centre) shows that the sources are loosely distributed in the Orion A region, and seem to cluster at $(l,b)\sim(209,-19)$. Fig. 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_{8}$ (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 ($\varpi<1.9\,\mathrm{mas}$), which are most likely field contaminants.

5 Ages

We determine ages ($\tau$) 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 $\bm{\theta}=(\tau,A_{V},Z,...)$, as:

[TABLE]

with:

[TABLE]

where $m$ is the stellar mass, $n$ is the number of observed quantities, and $\bm{q}^{\mathrm{obs}}$ and $\bm{q}(\bm{\theta},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 $j$th likelihood with a factor $w$ defined as:

[TABLE]

where $n_{redder}$ is the number of stars with $G_{BP}-G_{\mathrm{RP}}$ colour larger than that of the $j$th star and $n_{bluer}$ is the number of stars with $G_{BP}-G_{\mathrm{RP}}$ smaller than that of the $j$th 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:

[TABLE]

Since we are interested in determining the ages and the extinctions of the groups, we fix the metallicity to $Z=Z_{\odot}=0.0158$ and we integrate Eq. 13 on the mass, so that the probability density function as a function of age $\tau$ and extinction $A_{V}$ is given by:

[TABLE]

To perform the fit we compare the observed $G$ magnitude and $G_{BP}-G_{\mathrm{RP}}$ 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, using the passbands by Maíz Apellániz & Weiler (2018). We used isochronal tracks from $\mathrm{\log(age/yr)=6.0}$ (1 Myr) to $\mathrm{\log(age/yr)=8.0}$ (100 Myr), with a step of $\mathrm{\log(age/yr)=0.05}$ , and from $A_{V}=0\,\mathrm{mag}$ to $A_{V}=2.5\,\mathrm{mag}$ with a step of $0.1\,\mathrm{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 Section 6.2

5.1 Results

We compute the age $\tau$ and the $A_{V}$ for the groups identified by DBSCAN, and for the groups we selected in Section 4. The results are reported in Table 2. Figures 15 and 16 show the log-likelihood $\log L=\log L_{combined}(\tau,A_{V})$ we obtain for group $B_{0}$, and the $M_{G}$ vs. $G-G_{\mathrm{RP}}$ (left) and $M_{G}$ vs. $G_{BP}-G_{\mathrm{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_{\mathrm{RP}}$ colour, and we show the colour-magnitude diagram in $G-G_{\mathrm{RP}}$ as a quality check. We adopt the maximum of $L_{combined}(\tau,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.

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.

6.1 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}\sim 20\,\mathrm{km\,s^{-1}}$ and the second one with $v_{y,I}\sim 26-27\,\mathrm{km\,s^{-1}}$. 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. (2018) 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 $\tau_{\mathrm{exp}}=1/(\gamma\kappa)$ 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 $\kappa$ 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 Section 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.

6.2 Ages

The results obtained by fitting isochrones to the colour-magnitude diagrams of the groups isolated in Section 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_{8}$ 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_{1}$) 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 e.g. 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). 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\,\mathrm{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).

6.3 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 coordinates 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)\,\mathrm{pc}$.

Cantat-Gaudin et al. (2018) 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. (2018) 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 study 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 select young sources by applying simple cuts in the $M_{G}$ vs. $G_{BP}-G_{\mathrm{RP}}$ colour-magnitude diagram, and we study their density distribution in 3D galactic coordinates.

• •

We normalise our 3D density map between 0 and 1, and we select only the sources above a threshold of 0.5. We then apply the DBSCAN clustering algorithm to identify groups in 3D space and we analyse their properties in terms of ages and kinematics.

• •

We first inspect the proper motions of all the groups. In some cases we find that single groups in 3D space show sub-structures in their proper motion distribution. In this case we further sub-divide the groups, making simple cuts based on the proper motion distribution. We then apply a kinematic modelling code that we use to retrieve average motions, velocity dispersion, and isotropic expansion for all the groups identified.

• •

By comparing the 3D velocities of all the groups, we find evidence of kinematic sub-structures.

• •

We compute ages and extinctions for all the groups by using a 2D maximum likelihood approach.We find 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 do 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 obtain 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.

Acknowledgements.

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.

Appendix A Testing the kinematic modelling 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 ($\kappa$) value. In particular we are interested in the ability of the code to retrieve the correct value for $\kappa$, especially when not all the radial velocities of the cluster members are provided.

A.1 Simulation set up

The simulated star positions are drawn from Gaussian distributions with $\sigma=2\,\mathrm{pc}$. The velocity of each simulated star is drawn following the same assumption as in L00, i.e. from a Gaussian distribution centred in $\bm{v_{0}}$ with a small velocity dispersion $\sigma$. We include expansion following Eq. 9, chosing $\kappa=0.1\,\mathrm{km\,s^{-1}\,pc^{-1}}$.

We obtain the observed quantities (positions, parallax, proper motions, and radial velocities)222To do the transformation we make use of the pygaia routine phaseSpaceToAstrometry. by adding typical Gaia errors in the Orion region drawn from Gaussian distribution with widths $0.1\,\mathrm{mas}$, $0.1\,\mathrm{mas\,yr^{-1}}$, and $3\,\mathrm{km\,s^{-1}}$ 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 $\bm{v_{0,I}}=(-5.0,45.0,6.0)\,\mathrm{km\,s^{-1}}$ and distance $(x_{0,I},y_{0,I},z_{0,I})=(17.89,42.14,13.16)\,\mathrm{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)\,\mathrm{pc}$ and $\bm{v_{0,I}}=(0.0,20.0,0.0)\,\mathrm{km\,s^{-1}}$ . We run the simulations in five different scenarios for both the simulated clusters:

$\sigma_{v}=0.3\,\mathrm{km\,s^{-1}}$ and $\kappa=0.1\,\mathrm{km\,s^{-1}\,pc^{-1}}$. 2. 2.

$\sigma_{v}=1.0\,\mathrm{km\,s^{-1}}$ and $\kappa=0.1\,\mathrm{km\,s^{-1}\,pc^{-1}}$; 3. 3.

$\sigma_{v}=0.3\,\mathrm{km\,s^{-1}}$, $\kappa=0.1\mathrm{km\,s^{-1}\,pc^{-1}}$, and a fraction $f=[10\%,50\%,95\%]$ of stars without measured radial velocities.

The average velocities are always retrieved quite correctly in both cases; $\sigma$ and $\kappa$ are retrieved correctly for cluster A, however we notice that for cluster B the value of $\kappa$ is usually underestimated, while $\sigma$ 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 $\bm{v_{0}}$ from astrometric data only. In the cases when this happens, we do not give any estimate for the expansion term $\kappa$. 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 3). 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 ($\sigma=5\,\mathrm{pc}$), and their velocities are drawn from a Gaussian distribution with mean $20\,\mathrm{km\,s^{-1}}$ in $x_{I},y_{I},z_{I}$, and dispersion $\sigma_{v}=10\,\mathrm{km\,s^{-1}}$. 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 $\kappa$ is under-estimated (roughly by a factor of 2), while $\sigma_{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 $\kappa$ 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: $\sigma_{v}=0.3\,\mathrm{km\,s^{-1}}$ and $\kappa=0.1\,\mathrm{km\,s^{-1}\,pc^{-1}}$. We find that in general the minimisation results do not strongly depend on the initial parameters, however if the velocity dispersion $\sigma_{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: $\sigma_{v}=1.\,\mathrm{km\,s^{-1}}$ and $\kappa=0.1\,\mathrm{km\,s^{-1}\,pc^{-1}}$. We find that the minimisation results do not depend on the initial parameters in any case. This is reassuring, as the values for $\sigma_{v}$ in the clusters considered here are larger than $0.3\,\mathrm{km\,s^{-1}}$. In the cases with $\sigma_{v}=1.\,\mathrm{km\,s^{-1}}$ 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 $\kappa$, that is underestimated. If outliers are present, the parameters are retrieved correctly after the exclusion procedure.

Appendix B Color magnitude diagrams

Fig. 18 shows the colour magnitude diagram for the groups that we identified in Section 4.