Discovery and astrophysical properties of Galactic open clusters in dense stellar fields using Gaia DR2

Dividing the region in 4 sectors of same area diminishes the amount of stars in the vector point diagram (VPD) and colour-magnitude diagram (CMD), making it easier to find overdensities, as it increases the relative number of cluster stars over the field stars (i.e. the contrast) in both parameters spaces. Fig  2 shows the spatial filter applied to one of the analysed fields by dividing the full region in 4 smaller slices. The field is centred in the coordinates RA=257.632^∘ and DEC=-36.3409^∘ and contains $\sim$ 2 million sources. Each slice contains $\sim$4.0$\times 10^{5}$ sources (red), $\sim$7.4.$\times 10^{5}$ sources (blue), $\sim$2.7$\times 10^{5}$ sources (green) and $\sim$6.8$\times 10^{5}$ sources (purple).

For each sector, we build a $G$ versus $(G_{BP}-G_{RP})$ CMD and applied colour and magnitude filters. The first filter is a simple magnitude cut that keeps in our sample sources brighter than G = 18 mag (similar to the filter used in Cantat-Gaudin et al.  2018a). This procedure produces smaller working subsamples and makes the following procedures computationally faster.

The second filter divides the remaining sample ($G<18\,$mag) in two sub-samples by a colour value. This colour value is a free parameter and will depend in the local reddening as we visually divide the denser portions of the CMD in half. Usually, the first sub-sample keeps stars with $G_{BP}-G_{RP}$ smaller than a limiting value between 1.5 and 2.5 mag. The second sub-sample will contain stars with $G_{BP}-G_{RP}$ larger than this limit. This first filter is important to separate highly reddened stars, thus enhancing the cluster main sequence. A lower colour cut value is suitable to probe younger clusters with low reddening, while higher values are appropriate to probe older and/or more reddened ones.

Using region 3 of Fig 2 as an example, the colour and magnitude filters reduced the number of sources from $\sim$2.7$\times 10^{5}$ to $\sim$2.3$\times 10^{4}$ in the bluer side of the cut, and to $\sim$1.1$\times 10^{4}$ in the red side of the cut (see Fig 3).

With the subsamples settled, we analyze the stars distribution on sky charts and VPDs, visually searching for overdensities in both spaces. If any overdensity is detected in the VPD, we apply a box-shaped mask of size 1 mas yr^-1 on it and analyze the spatial and parallax distribution of stars inside it. If this overdensity represents the typical proper motion of stars in a cluster, we expect a gaussian distribution of the parallax values and a spatial bound between the stars (usually with a dense core surrounded by a sparse halo), both with few outliers. This procedure is shown in Fig.  4. The VPD (top-left) contains $\sim$2.3 $\times 10^{4}$ sources and after applying the proper motion filter (top-right), this number is reduced to $\sim$ 200 stars, as showed by black dots (bottom-left) and histogram frequencies (bottom right).

The majority of our new clusters were found by detecting a VPD overdensity, however there were some few cases where the colour cut revealed an overdensity more clearly in the sky chart, only becoming evident in the VPD after application of a spatial mask, as shown in the example below.

If an overdensity is detected on the sky charts, we apply an extraction mask of variable size, that depends on the visual concentration of the cluster. Usually, we adopt values between $\sim$5′ $\times$ 5′ and $\sim$15′ $\times$ 15′  around its visual centre and analyze the subsample proper motion and parallax distributions. If the overdensity detected on the sky chart represents a real cluster, the subsample spatially restricted must form a clump on the VPD space and a gaussian distribution of the parallax value, both with few outliers.

The Fig. 5 exemplifies such procedure, where we analyse the field 32 (Tab.4) comprising $\sim$1.5 million sources. The upper left and upper right panels show the data, before and after a colour cut defined by $G_{BP}-G_{RP}>2.5$, respectively. After the sectioning and colour filtering, as explained before, the sample of $\sim$2.9$\times 10^{5}$ sources reduces to $\sim$7.7$\times 10^{3}$ sources. The middle left panel shows a VPD of the colour filtered sample. Even though, in this case, it is not possible to note a clear overdensity in the VPD, it can be seen an obvious concentration in the filtered sky chart (middle right panel). After applying a spatial mask around an overdensity at edge of the region and analysing this subsample in a VPD (bottom left panel) and their parallax values (bottom right panel), we note the presence of a new cluster candidate UFMG 54.

The majority of the clusters studied in this work exhibit parallaxes values between $\sim$0.2 and 0.7 mas, where there is too much contamination by star fields. Therefore, diagrams of proper motions in right ascension or declination versus parallax were not very useful to find visual overdensities. The parallax is taken into account when we perform the membership assessment (Sect 4.2).

In both cases, after the candidate detection, we used the subsamples to build histograms of the quantities RA, DEC, $\mu_{\alpha}^{*}$, $\mu_{\delta}$ and parallax and computed their mode values. We use these values in the cluster analysis by restricting our sample inside a smaller parameter space. It is important to emphasize that none of the filters applied to find the clusters (colour filters, proper motion filters or spatial filters) are kept in the subsequent analysis steps. Instead, the modal quantities computed are used as first guess and then refined. For example, the clusters centres are refined by building radial profiles, the proper motion final values are calculated after assessing membership, etc.

With the centres and proper motion mean values computed in the previous analysis, we followed a similar procedure adopted in Ferreira et al. (2019) and selected stars inside 30′ from the objects centres and with proper motions inside a square box around its mean values. This proper motion filter removes the excess of field stars, increasing the contrast of the cluster population over the field and allowing the construction of the radial density profiles (RDPs). Our objects exhibit parallax smaller than 1.2 mas. In this case, a 1 mas yr^-1 box around the mean proper motion value is expected to encompass more than 3$\sigma$ of the distribution, assuming an average proper motion dispersion of 0.13 mas yr^-1, as found by (Cantat-Gaudin & Anders, 2020). Therefore, we adopted a 1 mas yr^-1 box. When the field contamination was not critical, we allowed a 2 mas yr^-1 proper motion box aiming at recovering more members, which was employed for seven clusters. On the other hand, for six critical cases, where the clusters candidates were less populous, we needed to adopt box size values of 0.8 mas yr^-1, to enhance the contrast against the field.

At this point we have defined our working sample for the radial profile procedure: best quality Gaia data within the square proper motion box and 30′ from the objects centres. We did not find any systematic effect introduced by this procedure that could bias our determination of the structural parameters. Afterwards, we built RDPs on which we estimated the size of each object, the local background density and refined the centres previously computed. We built the RDP by counting the number of stars in concentric rings around the initial centre and dividing it by the rings area. Independent RDPs were made by using different ring widths and then merged into a single, final RDP.

We refined the centre values by making small adjustments to the coordinates previously set, accepting the solution with maximum central density as the new centre values. We fitted a straight line to the stellar density for rings well beyond the objects core, typically between $\sim 15$′ and $\sim 25$′, to compute the sky background level ($\sigma_{bkg}$) and its uncertainty. The limiting radius ($r_{lim}$), defined as the radius where the stellar density reaches the sky level, resulted between 4′ and 20′  for all clusters. Fig. 7 shows examples of the RDPs and the derived sky density levels used to determine the limiting radius for some cluster candidates.

Finally, the clusters’ structural parameters were obtained by weighted fittings of the King (1962) analytical model over their RDP. We present some of the profiles fitted in Fig 7. The error of each RDP datum correspond to its Poisson uncertainty and was used as weight in the fitting process. The structural parameters: central density ($\sigma_{0}$), core radius ($r_{c}$) and tidal radius ($r_{t}$) are presented in the Table 2. The number $N$ of cluster members is indicated in the last column of this table. This table is also available electronically through Vizier¹¹1http://cdsarc.u-strasbg.fr/vizier/cat/J/MNRAS/vol/page.

Most of our cluster candidates are projected against dense stellar fields. For this reason, in order to derive membership likelihoods to the stars, we applied a routine described in Angelo et al. (2019) that evaluates statistically the overdensity of cluster stars in comparison to those in a nearby field in the 3D astrometric space ($\mu_{\alpha}^{*}$, $\mu_{\delta}$, $\varpi$). We call the region centred in the cluster’s equatorial coordinates and restricted by its $r_{lim}$ as the cluster region. The comparison field is restricted by a ring-like region with inner radius $r_{lim}$+3′. The outer radius is such that the field area is equal to 3 times the cluster region.

In order to improve the decontamination algorithm performance and minimize the presence of outliers (i.e., small groups of field stars with non-zero memberships), we build a box-like filter in the astrometric space centred in the peak values of $\mu_{\alpha}^{*}$, $\mu_{\delta}$ and $\varpi$ computed in the Sec 2. We restricted the sample in box with width of 3 mas yr^-1 for $\mu_{\alpha}^{*}$ and $\mu_{\delta}$ and 1 mas for $\varpi$.

The 3D astrometric space is divided in cells with widths of 1.0 mas yr^-1, 1.0 mas yr^-1 and 0.1 mas in $\mu_{\alpha}^{*}$,$\mu_{\delta}$ and $\varpi$, respectively. These cell dimensions are large enough to accommodate a significant number of stars, but small enough to detect local fluctuations in the stellar density through the whole data domain.

A membership likelihood ($l_{\textrm{star}}$) is computed for cluster stars within each cell by using normalized multivariate Gaussian distributions, which formulation (see eq. 1, Angelo et al., 2019) incorporates the uncertainties and correlations between the astrometric parameters. The same procedure is performed for stars in the comparison field. The degree of similarity between the two groups (cluster and field stars) within a given 3D cell is established by using an entropy-like function $S$ (see eq. 3, Angelo et al., 2019). Stars within cells for which $S_{\textrm{cluster}}<S_{\textrm{field}}$ are considered possible members.

For cluster stars within those cells, we additionally computed an exponential factor ($L_{\textrm{star}}$) which evaluates the overdensity of cluster stars in the 3D space relatively to the whole grid configuration (see eq. 4, Angelo et al., 2019). This procedure was adopted to ensure that only significant local overdensities, statistically distinguishable from the distribution of field stars, would receive apreciable membership likelihoods.

Finally, cell sizes are increased and decreased by 1/3 of their original sizes in each dimension and the procedure is repeated. The final likelihood for a given star corresponds to the median of the set of $L_{\textrm{star}}$ values taken through the whole grid configurations. Fig. 8 exemplifies the application of such procedure to the cluster UFMG6.

To build the final members sample, we adopted a membership probability cut by keeping only stars with membership probabilities higher than 70$\%$. This is important to clean the sample of less probable members. It makes evident the presence of the clusters stars loci in the CMDs and the existence of a clump in the 3D astrometric space ($\mu_{\alpha}^{*}$,$\mu_{\delta}$ and $\varpi$), revealing the cluster nature.

We employed solar metallicity PARSEC-COLIBRI models (Marigo et al., 2017) to perform isochrone fittings on the decontaminated samples (i.e., stars with likelihoods > 70%) to determine age, distance and reddening. We adopted a reddening law (Cardelli et al., 1989; O’Donnell, 1994) to convert $E(G_{BP}-G_{RP})$ to $E(B-V)$. Using a set of several solar metallicity isochrones covering a range of ages, we visually inspected the match between the model and the cluster stars loci, ensuring a good fitting in specific evolutionary regions such as the main sequence, the turnoff, and the giant clump. To associate the uncertainties, we changed the reddening and distance modulus $(m-M)$ simultaneously searching for the maximum acceptable deviations from the best values fitted.

The resulting values of the colour excess, distance module and age for all the clusters are present in the Table 3. The full table containing the astrophysical parameters and member lists derived for the clusters stars in this work are available as electronic tables through Vizier. Fig. 9 shows the best-fitting isochrones to the field-cleaned CMDs of some of the clusters investigated.

Cantat-Gaudin et al. (2020) have obtained astrophysical parameters for 1867 open clusters previously identified with Gaia DR2 astrometry, including many discovered open clusters present in LP2019 and CG2020. Their reported parameters (private communication) for the 29 clusters in common with our sample tend to agree with our values for close (d < 2 kpc), relatively low reddening (E(B-V) < 0.9) and relatively old (log(t) > 8.8) or young (log(t) < 8.0) clusters. There are a few discrepancies which we attribute to the very reddened, and scarce nature of some clusters, often lacking a sizeable population of giant stars needed to properly constrain their astrophysical parameters.

The dispersion in proper motion is correlated with parallax for coherent stellar systems. Such a relationship is presented in Fig. 10 for the Gaia DR2 confirmed OCs (gray circles) according with Cantat-Gaudin & Anders (2020), our sample of recently discovered OCs (open red circles) and the newly discovered OCs (filled cyan circles). All OCs studied here have CMDs and structural properties compatible with real physical systems. This is corroborated by the fact that they lie in the diagram region where the bulk of confirmed OCs is.

Fig. 11 displays the distribution of the derived parameters for the 25 newly discovered OCs (filled cyan circles) compared with the recently discovered ones (open red circles). Panel (a) shows the reddening as a function of the heliocentric distance; panel (b) is the distribution of logarithmic age with the reddening; panel (c) presents the core radius versus the logarithmic age; and panel (d) shows the average stellar density of members within the tidal radius as a function of the concentration parameter. The properties of these two OC groups are similar, corresponding to low concentrated OCs, the majority of them located within 3 kpc from the Sun, with ages ranging from 30 Myr to 3.2 Gyr and reddening limited to $E(B-V)=2.5$.

Most of the OCs analysed in this work are located in the IV quadrant of the Galactic disc at low latitudes, therefore they are immersed in dense fields characteristic of the inner Milky Way. Fig. 12 shows the position of our sample projected onto the Galactic disc (panel a) and perpendicularly to it (panel b). The OCs are depicted in two age groups, below and above 400 Myr, indicated by open blue circles and filled red circles, respectively. The OCs tidal radii are represented by the symbol size. A schematic drawing of the spiral arms (Vallée, 2017), the four quadrant limits, the Galactic centre, the Sun’s location (at 8.34 kpc; Reid et al., 2014) and the central bar are shown in Fig. 12(a), while the OCs’ vertical distribution and scale-height of 60 pc (gray band; Bonatto et al., 2006) are indicated in Fig. 12(b).

The challenging task of separating fiducial members from field stars in a high density region of the Galactic disc was only possible due to the high precision of Gaia DR2 astrometric data and to the developed method optimised for cluster detection. In this method, as explained in Sect. 2, we employed a sequence of steps in the search of stars within delimited fields whose astrometric and photometric coherence are compatible with a star cluster. It is worth noticing that specialized automated searches with different algorithms using Gaia DR2 (Castro-Ginard et al., 2018, 2019; Cantat-Gaudin et al., 2018b, 2019;LP2020;Cantat-Gaudin & Anders, 2020;CG2020) are useful and necessary tools to detect such objects. Nevertheless, improvements are still needed for the detection of low concentration systems immersed in crowded fields. More akin to our technique, Sim2019 performed visual inspection of the Gaia DR2 astrometric and photometric data, resulting in the discovery of 207 new OCs. The 25 new entries reported in the present analysis amounts to $\sim$ 45 percent of our sample. Based on this, (we speculate that) a non-negligible fraction of OCs may have been missed by automatic detection algorithms.

Most of the studied OCs in the present work are inside the solar circle and have ages between 30 Myr and 3.2 Gyr (average and dispersion of $\log[t(yr)]=8.49\pm 0.52$), as can be inferred from Fig. 13, which presents the vertical distribution as a function of the Galactocentric distance for the Kharchenko et al. (2013)’s OC catalogue matched to Gaia DR2 confirmed OCs (Cantat-Gaudin & Anders, 2020). In this Figure, our sample is indicated by filled circles, with colours attributed according to the colourbar.

Among the 59 OCs in our sample, 27 are older than 500 Myr. Among the 25 new discoveries, 10 are older than 500 Myr. This suggests that nearly equal numbers of young and old clusters contribute to the inner disc uncompleteness. The disc flaring towards larger Galactocentric distances (see Fig. 13) is a well-known feature of the Milky Way morphology, detected not only with OCs (e.g. Bonatto et al., 2006), but also using field stars (e.g. Amôres et al., 2017) and H I (e.g. Kalberla & Dedes, 2008).

If surviving the early gas ejection phase, OCs change their structure with time in response to external dynamical effects combined with stellar evolution mass loss and internal gravitational interactions between the stars (Portegies Zwart et al., 2010; Renaud, 2018; Krumholz et al., 2019, and references therein). These processes lead the OCs to dissolution through evaporation with loss of stars that in some cases produce tidal tails (Bergond et al., 2001). External factors like the tidal Galactic field strength and OC encounters with the spiral arms (Martinez-Medina et al., 2016), molecular clouds (Spitzer, 1958; Gieles et al., 2006) and disc shocking (Ostriker et al., 1972; Martinez-Medina et al., 2017) are more relevant in the denser environment of the inner Galactic regions. As an OC revolves around the Galaxy in $\sim$ 200 Myr at the Galactocentric distance of Sun’s orbit, its content is depleted by the aforementioned effects, which are more effective the older the OC is. For OCs orbiting closer to the Galactic centre, the external tidal field is stronger and the disrupting gravitational encounters proceed in a shorter time-scale. On the other hand, in the outer disc, OCs can expand without losing stars, surviving longer. Therefore, old clusters tend to be found in the outer disc (van den Bergh & McClure, 1980; Lynga, 1982; Janes et al., 1988).

Our sample follows the general distribution of tidal radius with age expected for OCs orbiting in the inner Galactic disc. This can be seen in Fig. 14, where Gaia DR2 confirmed OCs (Cantat-Gaudin & Anders, 2020) were matched to Kharchenko et al. (2013)’s catalogue from which ages and tidal radii were extracted for 1140 OCs. The colourbar indicates the Galactocentric distance from Cantat-Gaudin & Anders (2020), calculated with a heliocentric distance of $8.34$ kpc (Reid et al., 2014). Our results for the OC sample characterized in the present analysis are plotted in Fig. 14 as filled circles.

Being inserted in the inner Galaxy under relevant Galactic tidal field stresses and higher probability of encounters with molecular clouds and spiral arms, our sample contains mostly moderate size OCs, with $r_{t}=7.9\pm 3.7$ pc, compared to the following overall averages and dispersion for the confirmed sample of OCs: $r_{t}=8.3\pm 4.9$ pc for $R_{G}\leq 9.5$ kpc (663 OCs) and $r_{t}=13.0\pm 7.9$ pc for $R_{G}>9.5$ kpc (477 OCs). The interplay of tidal radius, age and galactocentric distance reveals the connection between OC structure evolution and the Galactic environment where they live. On the observational side, few works were dedicated so far to this aspect of the OC population in the Milky Way (e.g. Nilakshi et al., 2002; Chen et al., 2004; Angelo et al., 2020).

Fig. 15 shows the distribution of derived parameters for our sample OCs compared with that for Gaia DR2 confirmed OCs. The literature structural parameters, age and reddening were extracted from Kharchenko et al. (2013)’s catalogue while the distance was obtained from Cantat-Gaudin & Anders (2020). From a qualitative analysis of the structural parameter distributions, as shown in Fig. 15(a,b,c), our sample contains OCs with the mode of the core radii larger than that for the literature sample. Given that the distribution of tidal radii is not significantly different of that for the literature confirmed OCs, our sample OCs tend to be less concentrated. The ages present similar distributions, while most of our sample OCs have larger reddenings and smaller distances, as can be seen in Fig. 15(d,e,f), reflecting their low Galactic latitudes, where dust is concentrated and the stellar field is crowded. With many of the recently discovered clusters being closer than $\sim 2$ kpc, it seems that their number is far from complete as discussed by Cantat-Gaudin et al. (2018b) and Cantat-Gaudin & Anders (2020).