Low mass young stars in the Milky Way unveiled by DBSCAN and Gaia EDR3. Mapping the star forming regions within 1.5 Kpc

It is by now well known that stars originate from collapse of cold molecular clouds, and mainly form in over-dense structures and clusters usually designated as star forming regions (SFRs). During the very early phases, young stellar objects (YSOs) can be identified in the near, mid far infrared (IR) and radio wavelengths because of the presence of the optically thick infalling envelope or circumstellar disk around the central star. In the subsequent pre-main sequence stage phase, they become visible also in the optical bands. But, when the final dispersal of the disk material occurs and non-accreting transition disks form, YSOs can no longer be identified in IR or radio surveys (Ercolano et al., 2021) and a complete census is possible only in the optical bands.

While a clean identification of YSOs is very hard using only optical photometry, an efficient way to systematically single out SFRs is by the identification of kinematical stellar groups having a common space motion. With an unprecedented astrometric precision and sky coverage, Gaia data offer the possibility to recognise the SFRs as common proper motion groups, at least within the Gaia observational limits.

Data from the Gaia mission are revolutionising our capability to map the youngest stellar populations of the Milky Way in the optical bands, which is one of the main core science goals for an overall understanding of the Galactic components. The youngest stellar component is crucial to better characterize the Galactic thin disk, and its spiral arms and to understand its origin.

The characterisation of individual SFRs and their dynamics are also fundamental to understand the local formation, evolution and dispersion of star clusters, as well as the star formation history and the Initial Mass Function. Finally, statistical studies of YSOs during the early years of their formation, when the proto-planetary discs are evolving and planets form, are crucial to shed light on planet formation theory.

With more than 1.3 billion stars with precise proper motions and astrometric (positions and parallaxes) and photometric measurements, Gaia DR2 data allowed several studies aimed to identify clustered populations of the Milky Way. Some of these studies have been dedicated to SFRs, associations and moving groups. Zari et al. (2018) presented an analysis of the clustered and diffuse young populations within 500 pc, using a combination of photometric and astrometric criteria. Analogously, Kerr et al. (2021) studied the solar neighbourhood by applying the HDBSCAN clustering algorithm (McInnes et al., 2017). They found 27 young groups, associations and significant substructures, associated to known clusters and SFRs, and released a catalogue including $\sim 3\times 10^{4}$ Gaia DR2 YSOs within 333 pc.

Cantat-Gaudin et al. (2018) started from a list of known clusters to assign them unsupervised membership and parameters. Other studies have been dedicated to systematically find open clusters in the Galaxy. Castro-Ginard et al. (2018) used the DBSCAN algorithm (Ester et al., 1996) to select a list of candidate open clusters (OC) which they then refined to identify real OCs with a well defined main sequence (MS). Other papers have been recently published to both discover new open clusters and derive their parameters (e.g. Cantat-Gaudin & Anders, 2020; Cantat-Gaudin et al., 2020; Castro-Ginard et al., 2020; Liu & Pang, 2019).

A recent attempt to find Galactic Plane (GP) clustered populations, including SFRs, has been made by Kounkel & Covey (2019) and Kounkel et al. (2020), again using Gaia DR2 data and the HDBSCAN unsupervised algorithm in 5D space ($l$, $b$, $\pi$, $\mu_{\alpha}*$ , $\mu_{\delta}$). In these works, the first limited to 1 Kpc and the second to 3 Kpc, they found clustered populations but also associations, moving groups and string-like structures, parallel to the GP, spanning hundreds of parsec in length. Clusters aged between 10 Myr and 1 Gyr have been found, with an onion-like approach, i.e. using the entire catalogue with different cutoffs in parallax and progressively merging the different catalogues.

A different approach has been adopted by Bica et al. (2019) who used infrared (IR) data from 2MASS, WISE, VVV, Spitzer, and Herschel surveys to compile a catalog of 10 978 Galactic star clusters, and associations, including 4 234 embedded clusters.

With the advent of Gaia Early Data Release 3 (EDR3), based on 34 months of observations¹¹1Gaia DR2 data were based on 22 months of observations, available photometric and astrometric measurements improved significantly. In particular, photometric improvements have been made in the calibration models, in the different photometric systems and in the treatment of the BP and RP local background flux (Riello et al., 2021).

In this paper, we use Gaia EDR3 data to systematically identify the low mass component of SFRs in the Galaxy, with ages approximately $<10$ Myr and within a distance limit of $\sim$1.5 Kpc, imposed by our data selection. We focus our analysis on very young clusters, by exploiting the significant progress achieved with Gaia EDR3 data. A full exploitation of the Gaia data and the results presented here would require further data, such as spectroscopic determination of individual stellar parameters, such as effective temperatures, gravities, and stellar luminosities as well as rotational and radial velocities, crucial to derive masses, ages and 3D space velocities. Even though the results presented here cannot be used at this stage to determine the IMF, Star Formation history and 3D kinematics of the SFRs, they can serve to trace the very young Galactic stellar component within 1.5-2 Kpc through a systematic method that homogeneously identifies the bulk population of the SFRs. Such results can be used both for statistical as well as for individual detailed analyses. The paper is organised as follows: we describe in Sect. 2 the requirements adopted to select the Gaia EDR3 data and in Sect. 3 the photometric selection applied to obtain the starting sample of the YSO candidates. In Sect.  4 we describe the method adopted to identify SFRs and stellar clusters, the criteria adopted to validate them and the age classification. Our results and the discussion are presented in Sect. 5 and  6, respectively; finally, our summary and conclusions are presented in Sect 7. In the Appendix A we show the effects of the reddening in the Gaia color-absolute magnitude diagrams, in the Appendix B we estimate the effect of multiplicity in the selection of the YSOs while in the Appendix C we describe the comparison of specific regions with the literature.

In this analysis, we use the Gaia EDR3 data (Gaia Collaboration et al., 2016, 2021) which provide precise astrometry and kinematics ($l$, $b$, $\pi$, $\mu_{\alpha}*$ , $\mu_{\delta}$ ) as well as excellent photometry in three broad bands ($G$, $G_{\rm BP}$, $G_{\rm RP}$). Since our analysis is focused on the Galactic midplane, where most of the YSOs are expected to be found, we selected sources within $|b|<$30^∘. We limit our selection to 7.5$<G\leq 20.5$. The limit $G=$7.5 has been chosen to discard objects with magnitudes derived from saturated CCD images, while $G=$20.5 is the limit to include most of the objects with magnitude $G$ uncertainties lower than 0.2 mag. This range includes the young low mass populations ($0.1\lesssim M/M_{\odot}\lesssim 1.5$) of the known SFRs within the distance set by the limiting magnitude. In addition, we considered only positive parallax values. This choice does not introduce any bias since we do not expect to investigate stars with very small parallaxes that could have negative values (Luri et al., 2018). Finally, we imposed a relative parallax error lower than 20%, to discard stars with a poorly constrained distance, and, to take into account the Gaia EDR3 systematics, we also considered the renormalized unit weight error (RUWE), (Lindegren et al., 2021b), expected to be $<1.4$ for sources where the single-star model provides a good fit to the astrometric observations.

To summarise, data of our interest were selected from the Astronomical Data Query Language (ADQL) interface of the ESA Gaia Archive²²2https://gea.esac.esa.int/archive/ using the following restrictions:

We included in the query also a photometric condition aimed to include the Pre-Main Sequence (PMS) region of the $M_{G}$ vs. $G-G_{\rm RP}$ color-absolute magnitude diagram (CAMD) where all very young stars (t$\lesssim$10 Myr) are expected to be found. We split our selection in two samples, namely bright and faint, according to the following criteria:

These limits are drawn as solid blue and green lines in Fig. 1. We note that in this work, for the reddening uncorrected absolute magnitudes, we adopted the definition $M_{G}=G+5\,\rm{Log}(\pi)-10$, based on the inverted Gaia EDR3 parallaxes, since, as shown in Piecka & Paunzen (2021), within $<$2 kpc, the inverse-parallax method gives results comparable to distances derived by the Bayesian approach (Bailer-Jones et al., 2021).

The minimum value $M_{G}=5$ was set to avoid the upper region of the color-absolute magnitude diagram, where the overlap of the Upper MS or PMS stars of the SFRs with giants, MS or Turn-off stars, is expected to be very high, especially if the reddening is not corrected. This implies a cut of the massive population of the SFRs but it does not represent an issue for our investigation since we are mainly interested in the rich low mass component of these populations.

In order to further reduce the fraction of contaminants we used also the condition $G-G_{\rm RP}>0.58$, that is the minimum expected unreddened color for low mass (M$\lesssim$1.2M_⊙) PMS (age $\leq$ 10 Myr) stars.

Our photometric selection as well as the subsequent analysis are based on the $G-G_{\rm RP}$ colors. This choice allows us to avoid the use of the $G_{\rm BP}$ magnitudes that for $G\gtrapprox$20 are strongly affected by the application of the minimum flux threshold, which overestimates the mean BP flux. This issue also affects the RP flux, but with a considerably lower effect in $G_{\rm RP}$ than in $G_{\rm BP}$ (Riello et al., 2021).

Once the data have been retrieved by the ESA Gaia Archive, parallax values were corrected by the zero point bias reported in Lindegren et al. (2021a), using the Python code available to the community³³3https://gitlab.com/icc-ub/public/gaiadr3\_zeropoint, that is a function of source magnitude, colour, and celestial position.

In addition, we performed a further data filtering by considering only objects with error in the $G-G_{\rm RP}$ smaller than 0.14 mag. Standard errors in the magnitudes were computed by using the propagations of the flux errors with the formulas:

where $FG$, $FG_{\rm BP}$ and $FG_{\rm RP}$ are the mean fluxes in the $G$, $BP$ and $RP$ bands, respectively, and $\sigma(G_{0})=0.0027553202$, $\sigma(G_{BP_{0}})=0.0027901700$ and $\sigma(G_{RP_{0}})=0.0037793818$, are the Gaia EDR3 zero point uncertainties⁴⁴4See https://www.cosmos.esa.int/web/gaia/edr3-passbandsforfurtherdetails.

In this section, we describe and discuss how we performed the final photometric selection of the sample used as input for the subsequent clustering analysis, that is based on the astrometric and kinematic Gaia EDR3 parameters, as will be described in Sec. 4.

By considering the typical complexity of the environment of young stars and the dependence of the reddening law from the stellar effective temperature due to the large spectral range covered by the Gaia bands (Anders et al., 2019), we do not attempt to correct colors and magnitudes for reddening and absorption but we use their observed values. This is certainly one of the main source of contamination by older stars that will be overcome as will be discussed later.

Our goal is to start from a complete sample, including all potential YSOs with ages $<10$ Myr, at least in the photometric range set as described in Sect. 2. In particular, we selected the objects with $M_{G}$ falling on the red side of the solar metallicity 10 Myr isochrone computed using the PISA models (Dell’Omodarme et al., 2012; Randich et al., 2018; Tognelli et al., 2018, 2020), in the $M_{G}$ vs. $G-G_{\rm RP}$ diagram shown in Fig. 1. To check if the selected photometric limit is compliant with our requirements, we compared it with the reddening uncorrected CAMD of some SFRs and young clusters for which membership has been recently derived by Jackson et al. (2022), based on the 3D kinematics of the spectroscopic targets. We find that the adopted 10 Myr isocrone delimits the PMS region of clusters, such as NGC 2264, Lambda Ori, Lambda Ori B35 and Rho Ophiuchi, that are in the age range ($t<10$ Myr) of our main interest. However, also members of $\sim$20 Myr old clusters, such as Gamma Velorum, fall completely in the selected photometric region, while members of $\sim$50 Myr old clusters, such as NGC 2451b, fall partially in the selected photometric region at $M_{G}\gtrsim 9$. Going to clusters with ages $t>50$ Myr the overlapping region occurs at fainter magnitudes.

Since the adopted isochrone is limited to 0.1 M_⊙, corresponding to $M_{G}$=10.7, the photometric limit at fainter magnitudes has been extrapolated using a linear extrapolation. To check the position of such extrapolation, we compared it with the empirical sequence by Pecaut & Mamajek (2013) for which mean stellar colors and effective temperatures are given down to M and L spectral types and that can be used as upper limit to the region we are interested on. Our photometric limit approaches such sequence and crosses it at $M_{G}\sim 13$. This ensures us to set an inclusive photometric selection close to the MS at the lowest mass tail. In fact, even though this implies the inclusion of stars older than 10 Myr, on the other hand, it avoides a bias against the selection of very young stars.

We note that for the photometric selection, the minimum and maximum $M_{G}$ associated to each observed star have been computed by considering the 1 $\sigma$ parallax uncertainties, that are dominating with respect to the magnitude uncertainties. The photometric selection with respect to the reference isochrone has been performed by considering the compatibility of $M_{G}$ magnitudes, with respect to their minimum and maximum values, i.e. they were selected if either their minimum or maximum value lie inside the selection region. At the end of this selection we remain with a catalogue of 18 057 300  Gaia EDR3 entries.

Performing a photometric selection as inclusive as possible, as we have done, implies the introduction of a significant contamination by old field or open cluster stars, mainly due to the uncorrected reddening , binarity or to the overlapping photometric region, at the low mass range, where the sensitivity of the $G-G_{\rm RP}$ colors in distinguishing PMS or MS stars becomes very low.

However, the contamination by field stars does not represent a strong issue for our clustering analysis, since they are not expected to share similar astrometric and kinematic properties. In addition, since we aim at investigating the low mass component of the SFRs, that is also the most dominant ($\gtrsim 80$% Lada, 2006), the statistical contrast with respect to field contaminants is expected to be favorable to detect them.

A more tricky effect of our inclusive photometric selection is that also clusters older than 10 Myr can partially fall in the selected region and be recognized as candidate clusters in the subsequent analysis. As shown in Fig. 1, at faint absolute magnitudes ($M_{G}>9$), the model-computed isochrones are not very sensitive to stellar ages, and tend to overlap, especially in the $M_{G}$ vs. $G-G_{\rm RP}$ diagram. In addition, spectral synthesis of M dwarf stars suffers from the accuracy of the adopted atmosphere models and/or from incomplete molecular data. The model-predicted colors of very low mass stars are therefore uncertain. A further complication is the observed discrepancy between radii and colors of low mass stars, likely due to the distorting effects of magnetic activity and starspots on the structure of active stars (Somers et al., 2020; Franciosini et al., 2021). All these effects cause a spread of the low mass MS and can bring magnitudes and colors of $\sim 100$ Myr old stars to fall in the region selected by us as compatible to stars with $t<10$ Myr. For all these reasons, as discussed for example in Jeffries et al. (2017), the ages judged from ”standard” isochrones are almost certainly underestimated due to a systematic bias.

At faint magnitudes, the fraction of old cluster members falling in the adopted photometric region decreases with cluster ages. Hence clusters of about 20-30 Myrs will be almost completely included in our selected sample, while at the age of 100-500 Myr only the low mass tail will be included. However, because of the adopted photometric limit, the low mass tails will be included only for relatively close clusters ($d<$500 pc).

As already mentioned before, a partial contamination by old cluster members in our photometric sample can occur also for bright stars ($M_{G}<9-10$) if their reddening or a binary status gives them observed magnitudes and colors compatible with the selected photometric region. As shown in the Appendix A, the effects of using colors and magnitudes uncorrected for reddening are expected to be more severe for reddened stars with spectral types earlier than G, in comparison with later spectral types, in the sense that the selected sample is expected to be contaminated mainly by these objects, falling in the brightest part of the photometric region adopted in this work. The implications of this contingency will be discussed in the following sections.

Finally, we also considered the possible effects due to the metallicity on the selection by considering 10 Myr isochrones for metallicity lower or higher than solar. The comparison shows that while YSOs with over solar (Z=0.020, [Fe/H]=0.2) or sub solar (Z=0.005, [Fe/H]=-0.45) metallicities would fall in the selected photometric region, very metal poor YSOs (Z=0.001, [Fe/H]=-1.10) would remain outside. However, as recently found by Spina et al. (2017) at Galactocentric radii from $\sim$6.5 kpc to 8.70 kpc, young open clusters and SFRs have close-to-solar or slightly sub-solar metallicities and therefore we conclude that no SFRs are expected to be missed for metallicity effects with our photometric assumptions.

Based on the adopted photometric selection, our data set encompasses substantially all YSOs with age $t\lesssim 10$ Myr and observed $M_{G}>5$, including the most reddened ($A_{V}<3-4$) that can be detected with Gaia. Even YSOs with accretion (e.g. Gullbring et al., 1998) or seen in scattered light (Bonito et al., 2013) or flares in M-type stars (e.g. Mitra-Kraev et al., 2005) are expected to be included in our sample. In fact, these phenomena affect the $G_{\rm BP}-G$ or the $G_{\rm BP}-G_{\rm RP}$ colors, causing the stellar colors becoming bluer than their photospheric colors, while, on the contrary, their effect on the $G-G_{\rm RP}$ colors goes in the same direction of the reddening, causing these latter colors to become redder.

We stress, however, that the constraint $M_{G}<5$, adopted to strongly reduce the contamination due to reddened Turn-off or MS stars, makes the selected photometric sample incomplete for the massive stellar component of the SFRs. A further expected missing stellar component is that of binary systems of the clusters, due to the restriction of the Gaia data to RUWE$<$1.4 (see Appendix B). In addition, since available data do not allow us obtaining reliable corrections for the reddening affecting colors and magnitudes of the selected YSOs, accurate stellar parameters such as individual stellar ages and masses will not be derived in the subsequent analysis. However, even though the results we aim to achieve are not suitable for investigations based on complete young populations or accurate stellar parameters, they are expected to trace the dominant component of the SFRs, i.e. their low mass population, and will be crucial to have an overall systematic view of the Galactic SFRs located within 1-2 Kpc from the Sun, as well as for detailed individual or statistical investigations of these YSOs.

In the previous sections, we described how overdensities in the 5D parameter space ($l$, $b$, $\pi$, $\mu_{\alpha}*$, $\mu_{\delta}$) have been identified, starting from a photometrically selected sample, that covers the expected PMS region of YSOs with ages $t<10$ Myr.

Since no attempt has been made to correct for interstellar reddening, the starting sample was contaminated also from older reddened stars. Another possible reason for the contamination of older stars derives from the adopted strategy to select the starting sample in the plane $M_{G}$ vs. $G-G_{\rm RP}$, where the sensitivity to stellar ages of the available isochrones is quite low for the low mass population. In fact, for faint and very low-mass stars, isochrones get closer and closer for ages larger than about 50-100 Myr and, consequently, it is difficult to separate young populations from the older ones. As a result, the DBSCAN clustering algorithm, adopted to resolve spatially concentrated and/or co-moving stellar populations located at the same distance, can select also clusters older than 10 Myr.

A pattern match procedure has been adopted to disentangle SFRs and young clusters from older and photometrically unphysical clusters. We found 354  SFRs with ages $t\lesssim 10$ Myr and 322  young clusters with ages approximatively in the range 10-100 Myr. We discuss here these validated findings in the context of the GP structure within $\sim$1.5 kpc from the Sun.

The maps of the young stellar clusters recognised by the DBSCAN clustering algorithm, most of them already known in the literature, have been shown in the previous sections, and specific spatial and kinematic details have been presented for some of them.

To identify clusters extended on scales larger than the 5^∘$\times$5^∘  boxes used in the analysis, we merged adjacent clusters with consistent proper motions and distances. This procedure has been applied to identify extended SFRs as a whole, as in the case of the Orion Complex or Sco OB2 UCL, with r₅₀ equal to $\simeq 17$^∘and $\simeq 15$^∘, respectively, that are among the most extended structures resolved in this work. In several cases, it identifies clusters that encompass multiple populations, as in the case of NGC 2264, that has been identified as a unique structure including also the close cluster IC 446 and other YSOs in the surrounding region. A more in depth analysis of the two clusters shows that their proper motions can be distinguished into slightly different sub-populations and thus our overall procedure to define clusters tends to include multiple sub-populations sharing similar properties, likely associated to the progenitor molecular cloud.

The question of the cluster and subcluster identification is a very complex issue that can be dealt at different spatial precision levels, required for a given analysis, as done, for example, for the MYStiX project, in Feigelson (2018), where a parametric statistical regression approach, providing hierarchical ellipsoid structures, has been adopted. The evidence of a wide range of central surface densities found in the MYStIX maps is in agreement with the different spatial morphology of the SFRs identified in this work.

Figure 10 shows the spatial distribution of the young stellar clusters found in this work, in three different bins of distance, [100, 600] pc, [600, 2000] pc and [100, 2000] pc. The young clusters are drawn by distinguishing them in the age bins $t<10$ Myr, 10 Myr$<t<100$ Myr and $t<100$ Myr. Note that clusters with 10 Myr$<t<100$ Myr have been found only in the solar neighbourhood ($<600$ pc) and thus are shown only in the [100, 600] pc distance range.

The distribution of SFRs ($t<10$ Myr) within 600 pc is dominated by the presence of big young structures crossing the GP such as the Orion and Perseus Complexes, Gamma Velorum (Pozzo 1), Lac OB1, under the GP, BH 23 (corresponding to Theia 80 in Kounkel & Covey, 2019) and RSG 8 close to the GP, Serpens, Alessi 62, Collinder 359 and Rho Ophiuchi, over the GP. The clusters with ages 10 Myr$<t<100$ Myr in the same distance range appear definitely more diffuse. Apart from the well known Sco-Cen association covering $\sim$60^∘  in longitudes, we detected as unique complex the likewise huge association in the Vela-Puppis region, including Trumpler 10, $\gamma$ Velorum, NGC 2457, NGC 2451B, as well as the associations around NGC 2232, Roslund 5 and Alessi 19.

Their positions appear to be connected to the clusters with $t<10$ Myr since they follow a spatial pattern crossing or very close to that of the SFRs. This suggests that they likely belong to a common star formation process encompassing at least two generations of YSOs, with the first generation including extended populations of dissolving young clusters and associations.

The large Sco-Cen association is connected to the Vela and Orion Complexes, confirming what already found by Bouy & Alves (2015) with Hipparcos data. These three regions are described there as three large-scale stream-like structures.

Going towards larger distances ($d\gtrsim 600$ pc), the SFRs show a more regular pattern, approximately parallel to the GP. The most prominent SFRs are ASCC 32 and Cep OB3b in the Cepheus, respectively under and over the GP at distance of $\sim$800-900 pc. Among the most distant SFRs with more than 300 members and distance $\gtrsim 1400$ pc, we detected NGC 2244, NGC 6530, NGC 6531 , NGC 2362 and FSR 0442.

The overall distribution of YSOs in SFRs with $d\lesssim 600$ pc traces a complex 3D pattern in the solar neighbourhood. In particular, in the Z vs. X edge-on Galactic projection (see bottom/left panel in Fig. 4 and top/left panel in Fig. 10) we find evidence of a projected inclined structure, mainly traced by the Orion, Vela OB2 and Rho Ophiuchi star forming complexes in the third and fourth Galactic quadrants and, by the Serpens, Lacerta OB1 and Perseus in the first and second Galactic quadrants. However, the SFRs falling in the Cepheus region do not follow this pattern. A global view of these structures and their spatial correlation with the surrounding nebular emission, suggests a pattern consistent with the results found in Molinari et al. (2010), where massive proto-clusters and entire clusters of YSOs in active SFRs are associated to clouds that collapse into filaments.

As already found in Zari et al. (2018), current data reveal a very complex 3D structure that cannot be simply described with the Gould Belt, i.e. the giant flat structure, inclined by $\sim 20$^∘  with respect to the Galactic Plane, first pointed out by Gould (1879). This insight was already suggested by Guillout (2001), who presented the first detection of the Gould Belt late type star population, and who proposed the alternative scenario of a Gould disk.

A more detailed representation of the young Galactic component in the Solar Neighborhood has been recently proposed by Alves et al. (2020), who determined the 3D distribution of all local cloud complexes by deriving accurate distances to about 380 lines of sight. They suggested that such 3D distribution could be described by a damped sinusoidal wave, that they call Radcliffe Wave, with an amplitude of $\sim 160$ pc and a period of $\sim 2$ Kpc. It crosses Orion, around a minimum, Cepheus (crest), North America and Cygnus X. This structure is separated and distinct from a second structure, indicated as ”split”, crossing Sco-Cen, Aquila and Serpens clouds. They propose that the Gould Belt is a projection effect of two linear cloud complexes. The spatial distribution of YSOs associated to SFRs that has been identified in our work shows much more complex and diffuse structures but, the two elongated linear structures suggested by Alves et al. (2020) approximately cross the borderline of the two separated structures visible in the X, Y map of Fig. 4, delimited by the SFRs indicated by Alves et al. (2020). This brings us to confirm that the local young Galactic component is very complex. While our data are broadly consistent with the Alves et al. (2020) findings, further investigations, including a more detailed analysis of the kinematics of the structures, based on the 3D space coordinates ($X,Y,Z$) and velocities ($U,V,W$) (e.g. de Zeeuw et al., 1999) are required to confirm the scenario and to find additional insights about their origin.

To gain further insights on the star formation history of the SFRs, it is crucial to derive more accurate stellar ages. However, we do not attempt to derive stellar ages of the selected YSOs, for several reasons: first of all, the lack of a suitable photometric system. In fact, the large Gaia EDR3 G and RP photometric bands used for this work are not sensitive to the fundamental stellar parameters (effective temperatures, stellar ages, etc…), especially for low mass stars. However, future Gaia releases, overcoming issues related to the BP bands at faint magnitudes, could be crucial to this aim. Second, the lack of spectroscopic data needed to derive individual stellar reddening values, to appropriately place these YSOs on the HR diagram. Alternative ways, such as the use of 3D reddening maps (e.g. Bovy et al., 2016; Lallement et al., 2019), require careful analysis, since the integrated extinction tends to be underestimated in the molecular clouds, where SFRs are typically located. A detailed analysis is deferred to future works based on the combination of Gaia and spectroscopic data from available surveys, such as Gaia-ESO (Gilmore et al., 2012; Randich et al., 2013), LAMOST (Zhao et al., 2012) or APOGEE (Majewski et al., 2017), or future surveys such as WEAVE (Dalton et al., 2012) and 4MOST (Guiglion et al., 2019).

We used the machine learning unsupervised clustering algorithm DBSCAN to systematically identify all SFRs with ages $t\lesssim 10$ Myr, within $\sim$1.5 Kpc from the Sun. The density-based clustering algorithm has been applied to the Gaia EDR3 positions, parallaxes and proper motions of a photometrically-selected starting sample.

A pattern match procedure based on a template data set including typical clusters detected within the photometric sample has been used to distinguish very young clusters from the contaminant old clusters and from photometrically unphysical clusters. We provide here a catalogue with the main parameters (positions, spatial extent, median distance and number of members) of the 354  SFRs with ages $t\lesssim 10$ Myr. The parameters of the 322  young clusters with ages 10 Myr$\lesssim t\lesssim 100$ Myr are also given. We provide also the list of the 124 440  plus 65 863  YSOs found in the SFRs and the young clusters, respectively, including mainly late type K-M stars. A substantial number of YSOs have been recognized for the first time. Based on the comparison of our list of YSOs in the well known region Sco-Cen and NGC2264, we roughly estimate that, within our observational limits, the completeness of the census of cluster members obtained with our analysis is $\gtrsim$85%, at least in very rich and concentrated SFRs. For low density regions, such as the Taurus-Auriga association (see Appendix C), this completeness figure is expected to be around 50%. The mass function coverage of each cluster, strongly depends on the cluster distance, and is set by the observational limit.

Compact regular clusters, as well as SFRs in large complexes such as, for example, Taurus, Orion, Sco-OB2, Perseus, and Cygnus, have been identified with high efficiency, as estimated from the comparison with other available catalogues (see Appendix C ).

The overall distribution of these clusters in the Galaxy context shows that they are distributed along a very complex 3D pattern that seems to connect them at least within 500-600 pc. Outside of this distance, the clusters appear to be more regularly and closely distributed along the GP.

As far as we know, the catalogue of YSOs presented in this work is the sole all-sky catalogue based on the most recent Gaia EDR3 data, which benefit of major improvements with respect to Gaia DR2. This catalogue represents a step forward in the census of SFRs and can be used, for example, for further detailed interpretations of their spatial distribution in the context of the spiral arm model (Reid et al., 2019), since it covers a substantial region crossed by the Local arm and, marginally, some regions of the Perseus and Sagittarius-Carina arms (Poggio et al., 2021). Future and photometric deep surveys, such as the Rubin Legacy Survey of Space and Time (LSST) will allow to extend these limits.

We note that these results are not at this stage suitable for studies such as star formation history, cluster dynamics, based on the full space 3D velocity determination, or IMF, since the census of the SFRs is not complete and accurate masses and ages, as well as radial velocities can not be determined, until further data are available. Nevertheless, the dominant component of the SFRs has been detected and thus these results can be used as driving samples for the extraction of complete populations from Gaia data, by relaxing the stringent constraints adopted in this work.

Finally, the SFRs identified in this work are defined well enough to allow detailed studies of circumstellar disk evolution and direct imaging of young giant planets, based on multiband analisys of available or future additional observations (X-rays or IR or radio), targeting some of the individual clusters.