Press Release
Free Access
Issue
A&A
Volume 645, January 2021
Article Number A84
Number of page(s) 23
Section Galactic structure, stellar clusters and populations
DOI https://doi.org/10.1051/0004-6361/202038610
Published online 18 January 2021

© ESO 2021

1. Introduction

Star clusters are fundamental building blocks of galaxies. To this day, they are oftentimes thought to be the birthplace of most stars (e.g., Lada & Lada 2003; Parker & Goodwin 2007). Over the past decade, this view has been challenged on multiple occasions (e.g., Elmegreen 2008; Bressert et al. 2010; Kruijssen 2012; cf. Parker & Meyer 2012), where the Gaia mission (Gaia Collaboration 2016) and particularly the second data release (Gaia DR2; Gaia Collaboration 2018) have recently initiated a renaissance in our understanding of the predominant conditions for star formation (e.g., Wright & Mamajek 2018; Cantat-Gaudin et al. 2019a; Kuhn et al. 2019; Wright et al. 2019; Ward et al. 2020). These insights have been promoted mostly by investigating the large-scale structure and kinematics of massive OB associations with ages up to a few million years. However, exploring the full population of more evolved stellar structures becomes increasingly difficult as stellar systems quickly disperse and blend in with the Galactic field (Lada & Lada 2003; Pellerin et al. 2007). As a consequence, the dissociation of young stellar structures presents a challenging obstacle for our understanding of star formation demographics, and in particular, the dynamical evolution of stellar systems.

The dissolution process of open star clusters set inside a Galactic tidal field has been a popular subject in the astrophysical literature in the past decades. The long-term evolution of star clusters has been shown to be dependent on a great variety of variables that include, for instance, galaxy properties, the position in the host galaxy, orbital characteristics, nonaxisymmetric perturbations, tidal stripping, or radial migration (e.g., Spitzer 1958; Krumholz et al. 2019). On timescales of several hundred million years, dense star clusters dissolve by gradually releasing stars into tidal tails (e.g., Ernst et al. 2011; Renaud et al. 2011), features that have only recently been confirmed from an observational point of view (Meingast & Alves 2019; Röser & Schilbach 2019).

At earlier stages, theory suggests that a crucial phase is the transition of a cluster from an embedded to a gas-free state, usually referred to as residual gas expulsion (e.g., Baumgardt & Kroupa 2007; cf. Dale et al. 2015). Briefly summarized, star clusters are not expected to be in virial equilibrium after gas expulsion (Ernst et al. 2015). As a consequence, clusters experience a phase of violent relaxation, so that individual young stars can become gravitationally unbound, and start to expand. Theory establishes the critical factors for the degree of expansion and the disruption of young clusters to be the star formation efficiency (SFE, e.g., Lada et al. 1984; Shukirgaliyev et al. 2017), the timescales of gas expulsion (Dinnbier & Kroupa 2020), the dynamical properties of the cluster at the onset of gas expulsion (Goodwin 2009), primordial mass segregation (Brinkmann et al. 2017), and subsequent violent relaxation (Shukirgaliyev et al. 2018). The relative importance of these factors continues to be openly debated and remains unresolved, mostly because we lack unbiased observational constraints.

Understanding the early evolution of stellar systems is further complicated by the highly substructured nature of star formation. Several studies of young star-forming regions do indeed portray an irregular picture with substantial substructure (e.g., Cartwright & Whitworth 2004; Gutermuth et al. 2008; Kuhn et al. 2014). For instance, the nearest star cluster containing massive stars, the Orion nebula cluster, is in fact part of a much larger, highly elongated, substructured, and actively star-forming molecular cloud complex spanning about 90 pc (Großschedl et al. 2018; Hacar et al. 2018). The Orion cloud complex itself is embedded in the recently discovered Radcliffe Wave, an undulating 3 kpc long string of star-forming regions (Alves et al. 2020). Particularly in dense environments, Kruijssen et al. (2012) suggested that the cluster dissolution process could be dominated by tidal interactions in the greater birth environment of the natal molecular cloud complex.

In particular, a recent discovery highlights the complex interplay between the dynamical evolution of stellar systems and the physical conditions at birth. Meingast et al. (2019, hereinafter referred to as Paper I) reported the discovery of the large stellar stream Meingast 1 (see also Ratzenböck et al. 2020, hereinafter referred to as Paper IV) with an age of only about 120 Myr (Curtis et al. 2019). The stream seems to be more massive than but coeval with the Pleiades star cluster, but appears in stark contrast morphologically: Meingast 1 lacks a prominent cluster core, which means that if it formed as a monolithic cluster, it did not retain a centrally peaked structure. Instead, the population stretches across several hundred parsec and was only recently discovered because homogeneous radial velocity data across the sky became available. As a consequence, Meingast 1 is a prime example that highlights how different physical conditions at early stages can lead to vastly different dynamical evolution scenarios (see also Röser & Schilbach 2020).

Collecting observational evidence for the dissolution process of star clusters has been severely hampered by the inherent difficulties connected to the identification of stripped cluster members (e.g., Boss 1908; Artyukhina & Kholopov 1964; Eggen et al. 1987; Dalessandro et al. 2015; Yeh et al. 2019). The reason for this is that stars located beyond the tidal radius of a system are immersed in an overwhelmingly large number of unrelated field stars. Tracing population members out to several tidal radii and even beyond entails recovering individual sources orders of magnitude below the typical field star density. As a result, observers are confronted with the infamous problem of finding a needle in a haystack. Stripped cluster members cannot be identified with the classical approach of locating overdensities with respect to the background population. Instead, kinematic and chemical information must be used to identify stellar siblings reliably (Kamdar et al. 2019). While elemental abundances are only now becoming available in moderately sufficient quantities, the Gaia mission has ushered kinematic measurements into the age of big data. Nevertheless, identifying stellar populations across large regions of the sky remains challenging because Gaia primarily provides proper motion measurements, whereas 3D velocities are available only for a relatively small subset. In particular, projection effects and significant measurement errors in geometric distances still pose delicate challenges even for sophisticated machine-learning applications. In this paper we address these issues and introduce new methods for mitigating common limitations present in popular member-identification and cluster-analysis procedures.

The work presented here focuses on the analysis of ten nearby young open star clusters1: α Per (Melotte 20), Blanco 1, IC 2602, IC 2391, Messier 39, NGC 2451A, NGC 2516, NGC 2547, Platais 9, and the Pleiades. We selected these particular clusters because they are believed to be well known and demonstrate the capabilities of our methods extraordinarily well. Moreover, these clusters span a suitable range in age from about 30 to 300 Myr (see Table 1) so that their dissolution process may be analyzed, and they are located in the solar neighborhood. Figure 1 shows postage stamps for each cluster from the Digitized Sky Survey2 scaled to an edge length of 10 pc at the distance of the cluster core. These clusters all share a prominent appearance on the sky, where some of them have been known for centuries or even millennia. Despite their long-lived history in astrophysical research, we present unambiguous evidence that these apparent clusters comprise only the proverbial tip of the iceberg and are in fact each surrounded by a vast stellar halo, which we refer to as their corona3.

thumbnail Fig. 1.

Postage stamps for our ten target clusters, extracted from the Digitized Sky Survey; each stamp has a side length of 10 pc at the distance of the cluster core, and is displayed so that Galactic north is up, and east is left. From top left to bottom right, the clusters are sorted by their peak volume density.

Table 1.

Cluster parameters obtained from the literature.

2. Data

We based our analysis exclusively on Gaia DR2. To mitigate the effect of inaccurate measurements and outliers, we used the same quality cuts as employed in Paper IV. Specifically, the constraints we used are ϖ/σϖ > 10, fBP/σfBP > 10, fRP/σfRP > 10, astrometric_sigma5d_max < 0.5 mas, visibility_periods_used > 6, and the renormalized unit weight error, RUWE < 1.4 (Lindegren 2018). Here, ϖ, fBP, and fRP denote the parallax and the measured flux in the BP and RP bands, respectively, while σϖ, σfBP, and σfRP refer to the standard error in the corresponding parameters. To obtain a thorough census of the population of each cluster, we further limited parallaxes to values greater than 1 mas, corresponding to a geometric distance limit of about 1 kpc. These cuts together result in a remarkably clean if conservative sample comprising a total of 31 682 087 sources. All further steps in the data analysis, the results, and their analysis are based on this set of sources.

We preferred to calculate distances by inverting the individual parallax measurements as opposed to methods that take various distance biases into account. The reason for this choice was mainly our rigorous constraint of a 10σ significance of the parallax measurements. Furthermore, methods such as described by Bailer-Jones (2015), for example, impose a characteristic length-scale depending on the global galactic field distribution that is not directly applicable to cluster populations. Moreover, we did not apply a systematic parallax offset either as detailed in Lindegren et al. (2018), for example, because the quoted value is a global average and can be significantly different on local scales.

In the main text of this paper, we refer to a series of parameters and reference frames. Most notably, XYZ coordinates and their time differentials UVW refer to Cartesian coordinates in a Galactic reference frame centered on the Sun, where the velocity components are not corrected for the solar motion. In this frame, X and its differential U point toward the Galactic center, Y and V point in the direction of Galactic rotation, and Z and W are positive toward the Galactic north pole. In some cases, we prefer Galactocentric cylindrical coordinates with vR, vϕ, and vz that denote the radial, azimuthal, and vertical velocity components, respectively. On local scales, theses cylindrical velocities have a nearly linear relation to their Cartesian counterparts, but better represent structures on scales where galactic dynamics become non-negligible (e.g., curvature of orbits). The details for the definitions of the reference frames and other parameters are largely adopted from Paper I. The main differences are the switch to a new set of variables for the galactocentric frame as defined in the astropy v4.0 frame defaults and the change of the sign for the galactocentric radial velocity (vR). A comprehensive list of the definitions, including references and a short description, is given in Appendix B.

3. Methods

3.1. Cluster membership in the literature

The past few decades have seen several applications of various methods for identifying and characterizing stellar populations from astrometric data. Recent examples of techniques applied to Gaia astrometry are UPMASK (Krone-Martins & Moitinho 2014) and applications rooted in the popular clustering algorithms DBSCAN (Ester et al. 1996) and HDBSCAN (Campello et al. 2013; McInnes et al. 2017). These algorithms have been used in extensive studies of stellar populations in the Galactic disk, for example, by Cantat-Gaudin et al. (2018), Castro-Ginard et al. (2019), and Kounkel & Covey (2019), respectively.

These methods identify overdensities in the five-dimensional parameter space spanned by on-sky source coordinates (α,  δ), proper motions (μα*, μδ), and parallaxes (ϖ). To some extent, photometric information is also used in the identification of distinct structures and membership determination. The applications of these algorithms have led to the identification of several hundred new open clusters in our local kiloparsec (e.g., Castro-Ginard et al. 2019, 2020), and to the classification of numerous allegedly comoving string-like groups, see Kounkel & Covey (2019).

Ideally, the full six-dimensional parameters space (3D positions and 3D velocities) should be used to identify comoving stellar populations. However, this information can only be computed for the relatively small subset of sources for which radial velocity measurements are available in addition to positional and proper motion data. Conversely, identifying physically relevant structures in the five-dimensional parameter space, given by the 3D position and two tangential velocities is significantly more challenging because of projection effects across large regions (≳ 10°) of the sky. In particular, connecting sources in proper motions (or tangential velocities) poses a delicate challenge because the proper motion signature of a comoving population (i.e., a group of stars sharing the same space motion) largely depends on our perspective. This problem is illustrated in Fig. 2, which shows the expected tangential velocities of stars that are comoving with the cluster IC 2391 across the entire sky in a Galactic coordinate frame. The top and bottom panels in the figure show tangential velocities in the direction of Galactic longitude (vt, l) and along Galactic latitude (vt, b), respectively. Both panels show a complex pattern of tangential motions with fluctuating amplitudes and spatially inhomogeneous variations (i.e., the spatial derivatives of these patterns are also variable). Even at distances of only a few degrees from the cluster center, the measured tangential velocities can be highly variable, making comoving populations exceptionally hard to trace in tangential velocities with unsupervised clustering algorithms.

thumbnail Fig. 2.

Tangential velocity signature for sources comoving with the cluster IC 2391 in Galactic projection. The top panel displays velocities along Galactic longitude (vl), and the bottom panel shows velocities along Galactic latitude (vb). The position of the cluster core is marked with a black star. The tangential velocities depict a complex, highly variable pattern and illustrate the challenge of identifying comoving populations in large regions of the sky.

Because of the challenges outlined above, it is important to highlight the value of confirming and interpreting the physical significance of identified structures. Here, Gaia data directly facilitate two important measures: (a) photometry that is consistent with a single-age population in the form of a well-defined main and possibly giant sequence in observational Hertzsprung-Russel diagrams (HRD); (b) coherent 3D space motions verified by a velocity dispersion that is significantly smaller than measured for the Galactic field. In reference to point (b), star clusters in general, but also their associated tidal features in particular, are expected to feature velocity dispersions on the order of only a few km s−1 (Küpper et al. 2010a,b; Chumak et al. 2010). In addition, young stellar populations are also expected to feature very small velocity dispersions (σv3D ≲ 5 km s−1; e.g., Kuhn et al. 2019).

Velocity dispersions calculated directly from observational data are inflated by individual outliers and measurement errors (in particular from radial velocities). To remedy this problem, in all our calculations of velocity dispersions, including other published member catalogs, we first removed all sources outside a 3σ range in any individual velocity component and then deconvolved the resulting distribution with a single 3D Gaussian (Bovy et al. 2011). Moreover, we calculated the velocity dispersion from galactocentric cylindrical components to mitigate additional inflation effects resulting from curving orbits of large structures. The 3D velocity and σv3D can only be computed from the subset of sources in a given sample that includes radial velocity data.

Figure 3 displays the velocity dispersions (left panel) and apparent sizes (right panel) of all structures with a minimum of ten members identified in Cantat-Gaudin et al. (2018) and Kounkel & Covey (2019), respectively. While Cantat-Gaudin et al. (2018) specifically targeted known clusters for member identification, Kounkel & Covey (2019) deployed HDBSCAN in a largely unsupervised approach. Figure 3 shows that Cantat-Gaudin et al. determined groups with overall very small velocity dispersions, where about two-thirds of the structures feature values below 5 km s−1, but with relatively limited sizes. Conversely, Kounkel & Covey (2019) established physically larger structures, but also with significantly larger (deconvolved) velocity dispersions. About 90% of the identified structures feature kinematically hot values above 5 km s−1, which raises questions about the effectiveness of the authors’ unsupervised approach and challenges the physical relevance of many identified structures. Furthermore, this situation changes only marginally when young groups (< 100 Myr) are considered alone. The cumulative deconvolved velocity dispersion for these young groups is shown in Fig. 3. For this sample, the vast majority of structures also shows values that far exceed the expected velocity dispersions for young populations.

thumbnail Fig. 3.

Cumulative histograms for the deconvolved 3D velocity dispersions (left panel) and apparent sizes (right panel) of identified structures as published by Cantat-Gaudin et al. (2018; red) and Kounkel & Covey (2019; blue). Cantat-Gaudin et al. characterized members close to known clusters with relatively small velocity dispersions. Kounkel & Covey identified much larger structures, but more than 90% of the populations feature velocity dispersions greater than 5 km s−1. Young structures (< 100 Myr) in the sample from Kounkel & Covey (2019; dashed blue line) also feature mostly relatively large velocity dispersions.

The particular open star clusters in our sample have been studied frequently. Nevertheless, the above-outlined results and their inherent limitations call for a more thorough treatment of the member selection procedure. In light of the complex task of identifying connected structures in five dimensions, we introduce a new procedure in the next section that mitigates several biases in a physically motivated approach.

3.2. Selecting comoving populations with tangential velocities

Perspective effects on tangential velocities (see Fig. 2) need to be taken into account when comoving structures are identified in large regions of the sky. For this reason, we implemented a supervised method that separates spatial from kinematic dimensions and also corrects for highly variable tangential velocity signatures. The method is inspired by the popular convergent point technique (e.g., van Leeuwen 2009) and requires prior knowledge about the bulk 3D space motion of the target population. Our method computes for each source the difference between the actual measured tangential velocity from Gaia and the expected tangential velocity, assuming the bulk 3D velocity of a particular cluster population. In the following detailed description of the method, this difference is referred to as ) for the right ascension and declination components, respectively, and is expressed in km s−1. Any source that shares the exact same space velocity with the input vector (the cluster bulk velocity) therefore fulfills km s−1. Sources with marginally different 3D velocities feature slightly larger offsets to the origin in this transformed parameter space. As a result, our method facilitates efficient filtering of kinematically incompatible sources and can be easily followed-up with well-established clustering methods in the XYZ space.

As input, our supervised method requires establishing the bulk cluster velocity at a given 3D space coordinate. We used the cluster membership tables published by Cantat-Gaudin et al. (2018), where we retained only sources with membership probabilities greater than 0.8. We then computed mean positions and velocities after applying a 3σ clipping iteratively for five times. The mean velocities in galactocentric cylindrical coordinates for each cluster are listed in Table 1. Based on the calculated mean XYZ position, we sliced the input catalog (see Sect. 2) to retain only sources in a generously large box with edge lengths of 300 pc in X, 500 pc in Y, and 100 pc in Z, centered on the 3D cluster positions determined before. These limits were determined iteratively and were chosen as a compromise between not confining our selection by these limits and reducing contamination in the transformed tangential velocity space. Within our already stringent constraints on data quality, this procedure typically leaves several hundred thousand sources in the sliced box of each cluster.

The next paragraphs outline the method behind our membership determinations in detail. Figure 4 illustrates all steps in the selection procedure in reference to the α Per group and serves mainly as a guidance for our explanations. The left and central columns show the spatial arrangement of sources in Galactic Cartesian coordinates, and the column on the right displays kinematic information. The panels in the top row feature kernel density maps of the sliced box around the cluster and furthermore indicate the line of sight and the cluster’s Galactic orbit. The XYZ and maps were constructed with Epanechnikov kernels with sizes of 5 pc and 0.25 km s−1, respectively. These particular kernel sizes are not physically motivated, but instead were chosen for visualization purposes. Already in the 3D position data, α Per appears as an easily discernible overdensity at the center of the panels, together with other prominent clusters (in this case, the Pleiades and Hyades). Similarly, the transformed velocity space in the right-hand side panel also features a well-defined overdensity near the origin, depicting sources with a kinematic profile similar to that of α Per.

thumbnail Fig. 4.

Member selection process for the α Per cluster. The left and center columns display Galactic XYZ coordinates centered on the Sun, and the right-hand side column shows deprojected tangential velocities. The top row shows kernel density maps centered on α Per that distinctly highlight the cluster populations. The line of sight (dashed line) and the Galactic orbit of α Per (solid line) are also shown. The bottom panels present the member selection procedure, in which first a radial filter is applied in deprojected tangential velocities (red sources), followed by DBSCAN clustering in the XYZ parameter space to produce the final selection (blue dots).

In the second step, our procedure filters sources that are kinematically incompatible with the cluster bulk motion. To this end, we experimented with several methods, ranging from simple thresholding to other popular clustering methods such as DBSCAN and OPTICS (Ankerst et al. 1999). However, for several of our target clusters (and also associations and clusters outside of our sample), the distribution in the plane did not resemble such clearly separable overdensities as in the case of α Per in the upper right panel of Fig. 4. As a consequence, the algorithms performed poorly in many cases, and we instead chose to apply a hard limit of 1.5 km s−1 about the cluster bulk motion. This step drastically reduced the number of remaining sources by more than two orders of magnitude to a few thousand kinematically compatible sources. This reduced selection is displayed in the bottom panels of Fig. 4.

In a third and last step we extracted overdensities in the XYZ space from the kinematically filtered sample. To this end, we used DBSCAN with ϵ = 15 pc and minPts = 5. Furthermore, we only retained core points from the DBSCAN algorithm in order to keep the contamination at the population perimeters as low as possible. We experimented with several parameter combinations, and with the aid of the resulting observational HRDs, again chose a compromise between a too confined selection (only cluster core) and too many contaminating field stars. The resulting final selection for α Per is displayed in the bottom panels of Fig. 4. For most of our target clusters, this procedure resulted in a single, clearly separable population. Only in very few cases did our DBSCAN setup also label other (assumed to be unrelated) populations, which we chose to not investigate further here.

The choice of the free parameters for our method (filter radius in , DBSCAN ϵ, and minPts) is accompanied by a set of caveats. For instance, the hard clipping in the transformed tangential velocity space only represents a subset of the kinematically compatible population because members may very well lie beyond this radius. Moreover, in reference to the top right panel in Fig. 4, the overdensity in velocities is clearly not circular and seems to extend beyond the 1.5 km s−1 limit. The overdensity instead appears to take an elliptical shape, which is also confirmed in the final selection in the bottom right panel. Other disadvantages relate to the fixed parameter setup for DBSCAN when searching for overdensities in XYZ. For instance, all field stars that by chance share a similar velocity and are located in the same volume are selected as well. Moreover, measurement errors in geometric distances along the line of sight result in elongated structures, thereby distorting the XYZ space. This effect depends on the measurement errors in geometric distance and is more pronounced for more remote populations. Taken together, our parameter choice was optimized for the target clusters and should be carefully evaluated when this method is extended to other populations.

3.3. Deconvolution of the spatial distribution

The spatial arrangement of our selections for the ten target clusters appears distorted through significant measurement errors in distances along the line of sight. For example, for the most remote cluster population in our sample, NGC 2516, the cluster core is at a distance of about 410 pc. Even with our rigorous quality criteria, we expect standard errors in distance up to 40 pc. With typical cluster core sizes of only a few parsec, we consequently observe significant elongation effects in the cluster appearance. However, for reliable determinations of physical cluster sizes and cluster kinematics, for instance, this effect needs to be accounted for.

To infer the underlying true spatial arrangement of the populations, we employed a novel algorithm that is described in detail in Meingast (in prep.). Here, we provide a brief summary of its functionality. In a first step, the algorithm constructs covariance matrices for each source with respect to the measurement errors in XYZ positions. This is done by randomly sampling the observed parallax from the (symmetric) error distribution. Errors in on-sky position (α, δ) are ignored because they are negligibly small. The XYZ distribution, informed by the covariance matrices of the given spatial arrangement, is then deconvolved with a mixture of Gaussians, a process referred to as extreme deconvolution in the literature (XD; Bovy et al. 2011). The dimensionality of the mixture model (i.e., the number of Gaussian components) is determined with the Bayesian information criterion (Schwarz 1978). The number of components for our target clusters is typically between 3 and 5.

This step results in the description of the spatial configuration of a population in the form of a probability distribution composed of a mixture of 3D Gaussians. In a second step, the method infers the geometric distance of each source along the line of sight by using the mixture model as distance prior together with a Gaussian likelihood from the parallax ϖ measurement by Gaia. To obtain the deconvolved distances, we sampled the posterior with the Markov chain Monte Carlo ensemble sampler (MCMC) emcee (Foreman-Mackey et al. 2013, 2019).

The results of this procedure are visualized in Fig. 5. The figure shows the XY distribution of the central parts of NGC 2516, where the luminance value of the points is scaled proportionally to volume density. The left panel displays the positions when the Gaia parallaxes are directly inverted. In this view, the cluster appears highly elongated along the line of sight. The panel on the right-hand side displays the median posterior positions in the MCMC sampling. Our method clearly infers a good approximation of the true shape of the cluster core.

thumbnail Fig. 5.

Spatial configuration of the central regions of NGC 2516 before (left) and after (right) deconvolving source distances. The original distribution clearly shows a substantial elongation effect along the line of sight (dotted gray line), while the extreme deconvolution (XD) infers the expected spherical arrangement of the cluster core.

3.4. Contamination

In this method, the tangential velocity filter effectively removes most field star contaminants. However, all sources that pass this filter and are clustered with the subsequent application of DBSCAN will become part of the final member selection. In this way, some stars that are part of the Galactic field or other populations which, by chance, share similar tangential velocities, will inevitably also contaminate our population selection. In the following description, we refer to these contaminating sources as background population.

To estimate the number of contaminating background sources, we exploited the specific workflow of our method. To this end, we compared the finally selected cluster population with all remaining sources after the initial constraint in tangential velocities (i.e., we compared the red and blue sources in Fig. 4). We then obtained a contamination estimate by comparing stellar volume densities of our extracted cluster populations ρcl to the median volume density of all remaining background sources ⟨ρbg⟩ in the corresponding Gaia subset. Specifically, we parameterized the volume density at the location of each source via the distance to its seventh nearest neighbor and chose the median background density as comparison to avoid edge effects. With these two parameters at hand, we obtained fractional contamination estimates for each source (fc, i) via

(1)

Following this definition, the parameter fc can be viewed as a measure of the contrast between cluster and background stars. For example, fc takes a value of 0.5 for a cluster member located at a position with equal cluster and background volume densities. On the other hand, fc is close to 0 for sources that are surrounded mostly by their population siblings, while rogue cluster members will have values close to 1.

Figure 6 displays the cumulative distribution of fc for each population. The figure reveals that the large majority of the determined cluster members show exceptionally low values of fractional contamination. Typically, the majority of sources is located at positions where members outnumber background stars by a factor of 100 or more (fc <  10−2). The fractional contamination increases toward the outskirts of a population, whereas the lowest values for fc are found at the center of the most compact clusters (Blanco 1, Pleiades, and NGC 2516). Blanco 1 is characterized overall by very low values of fc, which can be attributed to its location about 200 pc below the Galactic plane, where the field star population is already thinned out. This is also illustrated in Fig. A.2, which shows the source distribution for each cluster in the XY plane, color-coded with individual contamination estimates.

thumbnail Fig. 6.

Cumulative distributions of the source contamination fractions fc for each cluster. The contamination fraction fc denotes the volume density contrast for each source with respect to the surrounding background population.

In addition, the parameter fc also enables an estimate of the total contamination by summing over all fractional contamination measures for each source. The measured total contamination by background stars amounts to well below 10% for all populations, with a maximum of 8% for IC 2391, values as low as 1% for the Pleiades and Blanco 1, and a median of 3% across all ten target clusters. The estimated number of total contaminating sources is listed in Table 2 for each population, and fc values for each source are available at the CDS (Table A.1).

Table 2.

Computed cluster parameters and statistics.

4. Results

This section contains a presentation of our results from an observational point of view. All results presented in this section are based on our application of the methods outlined in Sect. 3 to the ten target clusters listed in Table 1. This procedure includes rigorous filters that are used to retain only sources that are both kinematically compatible with the cluster bulk motion and are spatially connected to the cluster cores. Moreover, we also minimized the effect of structure elongation due to large measurement errors in the geometric distances along the line of sight. We describe subjects including the observational HRDs of the clusters (Sect. 4.1), their spatial distribution (Sect. 4.2), present-day mass functions (Sect. 4.3), radial profiles (Sect. 4.4), and cluster kinematics (Sect. 4.5). Each individual topic is discussed only broadly, and our results should serve as a presentation of the observables.

4.1. Observational Hertzsprung-Russel diagrams

At the start of our investigation, we created observational HRDs for the member selection of each target cluster. In these diagrams, the main sequence of a single-age population appears as a narrow distribution, while a random selection would result in more broadly scattered absolute magnitudes and colors. We therefore used the HRDs mainly as a diagnostic tool to validate our selection and to show that the extracted source populations are likely coeval.

Figure 7 displays the HRDs of the ten target clusters, where each cluster is displayed in its respective color in addition to all sources in the initial 300  ×  500  ×  100 pc search box of each group. The diagrams show that the populations appear as exceptionally narrow main sequences. Likewise, the sequences appear in stark contrast to their respective background population. This observation is further confirmed by performing a two-dimensional version of the Kolmogorov-Smirnov test (Peacock 1983) between the determined cluster members and all sources in their respective search boxes in the color-absolute magnitude space. These tests deliver p-values that are equivalent to 0 with respect to the floating-point precision. The same test on an equal-sized random draw of the background sample typically results in p >  0.5. Another indicator for a clean selection procedure are the frequently well-visible stellar multiple sequences, even down to the lowest stellar masses. Without further consideration, we take note that especially NGC 2516 shows a remarkably distinct stellar multiple sequence. Following these findings and also including the low contamination rate as detailed in Sect. 3.4, we conclude that our procedure indeed extracts coeval populations.

thumbnail Fig. 7.

Observational HRDs for the ten target clusters. Members of each cluster are displayed as colored points. The gray distribution in the background comprises all Gaia sources in the initial search box of each population. Our selection results in remarkably clean and narrow main sequences that frequently also show well-visible binary sequences. This clear separation from the broad background distribution acts as validation of the successful application of our member selection procedure.

For Messier 39, α Per, and NGC 2516, we also find giant stars in their respective HRDs. The location of the sources on the giant branch for Messier 39 and α Per do not match an expected evolutionary track for these clusters, which likely makes them contaminants. The three sources on the giant branch in the NGC 2516 selection are HD 64320, HD 65662, and CD-60 1965. All of them are classified as giants (Houk & Cowley 1975; Keenan & McNeil 1989; Sowell 1987) and match an isochrone for the cluster age as listed in Table 1. Moreover, some selections also contain white dwarfs, with particular objects of interest in the very young clusters NGC 2451A and IC 2391. If these sources are part of the population, their progenitors must have had a mass similar to early B-type stars to already have produced white dwarfs.

We furthermore find several sources located between the main sequence and the white dwarf locus. These are particularly well visible at the faint end of the selection for the Pleiades. A closer examination reveals that these sources are almost exclusively characterized by an above-average BP/RP flux excess, indicating issues with the photometry. When the blue band is replaced with the G band in the HRDs (i.e., mG − mRP are plotted on the abscissa), this problem is mostly solved, which further indicates problems with the BP photometry for the few affected sources. We note here that our selection method does not include any photometric criteria, which means that these few scattered sources are likely not contaminants.

4.2. Spatial distribution

Figure 8 shows the spatial arrangement of the extracted populations in a Galactic Cartesian coordinate frame centered on the Sun4. For the purpose of highlighting the cluster cores, the lightness values of the colors represent volume density, with higher densities corresponding to lower lightness values. Evidently, the cluster populations extend far beyond their long-known central regions and are surrounded by vast stellar halos we term coronae, reaching sizes up to several hundred parsec in diameter. As a result, the cluster cores, together with their massive coronae, form large, coeval, and comoving extended cluster populations, each encompassing up to hundreds of thousands of cubic parsec of space. Moreover, the groups appear to be spatially intertwined, where NGC 2451A, IC 2391, and IC 2602, and α Per and the Pleiades share a considerable volume in which the populations mix.

thumbnail Fig. 8.

Spatial distribution of our selection for the ten clusters in heliocentric Galactic coordinates. The color scale is proportional to volume density, accentuating the dense central regions of the populations. In this view, the Sun is located at the origin, and the directions toward the Galactic center (positive X) and of the generic Galactic rotation (positive Y) are indicated with arrows. Cluster cores are surrounded by enormous stellar coronae with sizes up to several hundred parsec. In addition, we observe a characteristic pattern where extended cluster populations form elongated shapes whose leading arm is tilted toward the inner Galaxy. An interactive 3D version of this figure is also available online.

A closer inspection of the XY plane in Fig. 8 shows that the cluster populations altogether form highly elongated shapes, where the parts closer to the Galactic center form leading arms in the direction of Galactic rotation. In contrast, the sections located farther away from the Galactic center appear to lag behind. We also observe a variety of position angles of the coronae with respect to their orientation in the Galactic plane. While NGC 2451A, for example, is visibly more tilted toward the Galactic center, most of the other groups appear to have a similar (but not exactly the same) orientation, being more elongated in the direction of Galactic rotation. We note that the alignment of the coronae is clearly not correlated with the line of sight from our vantage point, an artificial effect that often originates in problems with respect to significantly enhanced errors in parallax or poor method deployment.

For most populations, we observe a correlation between vertical tilt (best visible in the top right panel in Fig. 8) and vertical velocity (listed in Table 1). In the groups with positive vertical velocity, the leading parts of Platais 9, α Per, IC 2391, and NGC 2516 are also oriented toward the Galactic north pole. Messier 39, NGC2451 A, NGC 2547, and the Pleiades have negative vertical velocity, and their leading parts are tilted toward the Galactic south pole. IC 2602 has the highest vertical velocity, but does not show a similar trend, and interestingly, the situation is reversed for Blanco 1 (negative vertical velocity, leading part tilted toward the Galactic north pole). It therefore appears that the groups are preferentially oriented toward their direction of motion, but this does not seem to be a universal occurrence.

Following the deconvolution procedure, we also redetermined the 3D position of the population density maxima (i.e., the cluster cores) from the Gaussian mixture model. These new distances are listed in Table 2. With reference to Fig. 8 we find that the cluster cores seem to be preferentially located toward the center of the populations. Only NGC 2516 seems to be an exception: the trailing arm of the corona appears less prominent than the leading part. We attribute this to observational limits and our constraints on data properties, however, that significantly thin out source densities at distances beyond ∼400 pc, which is the distance to the cluster core of NGC 2516.

4.3. Mass functions

In this section, we derive total mass estimates for each population in order to lay the foundation for subsequently establishing physical properties. Determining the total mass requires an extrapolation for low-mass sources because they are beyond a distance-dependent observational completeness limit.

We first obtained individual masses for each observed source by interpolating their absolute magnitudes and colors on a given isochrone. To this end, we used PARSEC isochrones (Bressan et al. 2012) with the Weiler (2018) revised Gaia DR2 passbands and adopted ages and (assumed to be constant) line-of-sight extinctions from Bossini et al. (2019). Figure 9 shows the resulting measured mass functions for each target population as colored histograms. As expected, sensitivity limits cause a truncation at the low-mass end that is particularly well visible for NGC 2516, the most remote population in our sample. Different observed maxima at the high-mass end are also visible, where some populations show a well-sampled distribution up to several solar masses, while others (e.g., NGC 2516) are apparently truncated. In contrast to low-mass stars that can drop out of our selection because of sensitivity limits, the high-mass range can be affected by both incomplete catalog coverage (too bright sources) and stellar evolution processes (supernovae in the past).

thumbnail Fig. 9.

Present-day mass functions (colored histograms) and the best-fitting truncated IMFs (solid black histogram) for each population. The completeness limits on the low-mass end appear different for each cluster because of the distance-dependent sensitivity limits. The completeness at the high-mass end is more difficult to evaluate because individual stellar evolutionary histories are governed by the independent ages of the groups. If the mass distribution follows the displayed truncated standard IMF, we estimate the missing mass fraction to be in a range from about 15% for IC 2602 to almost 45% for NGC 2516.

To facilitate an extrapolation to total present-day cluster masses, we compared the measured distributions to Kroupa (2001) initial mass functions (IMF). Specifically, we decided to sample IMFs truncated to the highest measured stellar mass because in our selection we can hardly distinguish between mass loss due to stellar evolution and an observational bias in the Gaia database. The limit on the low-mass end was set to 0.03 M, thus also accounting for potentially missing mass below the hydrogen-burning limit. We then obtained the best-fit mass function by minimizing the total mass difference between a given (truncated) IMF and the measured total masses. Because of completeness limits in the Gaia data, only the mass range between 0.5 and 2 M was used for this comparison. The lower limit of 0.5 M corresponds to an M1 dwarf with an absolute G magnitude of about 8.8 mag (Pecaut & Mamajek 2013; Evans et al. 2018). For a distance modulus of μ = 8.5 mag (i.e., a distance of 500 pc), this corresponds to an apparent G magnitude of 17.3 mag for the lower mass limit, which is well within the Gaia completeness threshold (Arenou et al. 2018) and even leaves an ample margin for interstellar extinction. The resulting best-fit mass functions are displayed in Fig. 9. The measured total mass Mtot and the total mass of the best-fit mass function are listed in Table 2. The missing mass fraction due to incompleteness at the low-mass end is between 15% (IC 2602) and 45% (NGC 2516), with a mean of about 30% across all ten clusters.

4.4. Radial profiles

Our completeness-corrected mass estimates, together with the deconvolved spatial distribution, enabled a determination of physically relevant geometric parameters, specifically, tidal and half-mass radii (rt and rh). For our purpose, we use the terms of the Jacobi and tidal radii interchangeably (cf. Binney & Tremaine 2008) and follow the relation outlined by King (1962) and Ernst et al. (2010) with

(2)

where G is the gravitational constant, M is the enclosed total mass, Ω is the circular frequency of a near-circular orbit at the distance of the population to the Galactic center, and β = κ/Ω, with κ depicting the epicyclic frequency. Here, we calculated the circular and epicyclic frequencies for each population with the MWPotential2014 potential as defined in galpy (Bovy 2015).

We continued to calculate cumulative enclosed total masses as a function of distance from the density maximum of each population. To account for sources beyond the sensitivity limit of Gaia, the measured enclosed total masses were multiplied by the individual incompleteness factors as determined from the mass functions in Sect. 4.3. At the same time and using the above equation, we also calculated tidal radii for the enclosed mass at a given radial distance to the cluster center. This procedure provided tidal radii as a function of distance from the population centers (or conversely, as a function of enclosed mass). The intersection of the (completeness-corrected) radial mass profile and the tidal radius profile then yielded the actual tidal radius for each system. The values for the determined tidal radii, along with the system half-mass radii, are reported in Table 2. Figure 10 illustrates this procedure. The measured tidal radius (i.e., the intersection of the profiles) is marked. The statistical errors for the tidal radii are most likely dominated by the mass incompleteness factor. Varying this factor by 20% results in changes of rt at about a 5% level, or about 0.5 pc.

thumbnail Fig. 10.

Radial mass profiles and tidal radii for each cluster population. The colored lines depict the enclosed mass fraction as a function of distance from the density maximum of each population. The solid black lines trace the tidal radius for a given enclosed total mass content. The intersection of these profiles identifies the tidal radius. The mass profiles themselves depict a variety of shapes, with some clusters showing clear kinks near the tidal radius. For the majority of the clusters, the mass profiles also reveal that the bulk of their stellar mass content is located beyond the tidal radii in the stellar coronae.

We also verified the determined tidal radii by comparing our selection to the surrounding Galactic field distribution. Specifically, the tidal radius should correspond to the radius, where the volume density of a population equals the density of the surrounding field. Similar to our approach to estimating the contamination fraction (Sect. 3.4), we parameterized the volume density with the seventh nearest neighbor of each source. While for the contamination we only considered sources that survive the tangential velocity constraint, for this comparison we retained all Gaia sources in the vicinity of each group. Figure 11 shows the ratio of the population volume density (ρcluster) to the volume density in the field (ρfield) as a function of distance from the center of a group, normalized to the calculated value for the tidal radius. For all groups and within about 10%, the determined tidal radii correspond to the radii where the cluster volume density drops below the density in the field. In addition, we observe that the volume density contrast of our selection compared to field stars ranges across about five orders of magnitude from well above 100 down to 10−3 for individual sources.

thumbnail Fig. 11.

Population volume densities normalized to the surrounding field as a function of distance from the cluster center, parameterized by the tidal radii. The observed volume density vanishes into the field at the tidal radius. Our member selection recovers sources across five orders of magnitude in volume density, down to three orders of magnitude below the Galactic field at the outermost rim of the coronae.

Inspecting Fig. 10 more closely, we find subtle differences in the radial mass profiles. Several groups, and in particular, α Per, NGC 2547, Blanco 1, and IC 2391, show a pronounced kink in their profiles with very steep gradients toward the center. Other groups appear to have a much flatter profile with a less pronounced core structure and a smoother transition between the central regions and the coronae beyond the tidal radius. In reference to Fig. 11, the coronae themselves appear to follow a smooth distribution out to a few tidal radii before hitting the noise limit of our method with the current Gaia data release. As an exception to the general appearance of the profiles, even the central parts of Platais 9 appear to lie only marginally above mean field density, and the radial profile is significantly flatter as well.

Figure 10 also shows that for several clusters, the enclosed mass fraction is well below 50% at the tidal radius. Comparing the measured tidal and half-mass radii in Table 2, we find that the half-mass radii of only four populations (Blanco 1, NGC 2547, NGC 2516, and the Pleiades) are smaller than their tidal radii. This finding indicates that the bulk of the stellar mass content is located beyond the perimeter where the gravitational potential of the clusters dominates. The coronae are therefore largely gravitationally unbound with respect to the cluster cores. Figure 12 visualizes the statistics on the mass fractions inside and outside the tidal radius, and the corresponding fractions are listed in Table 2. Here, the statistics for NGC 2516 may be biased because we likely only recover parts of the trailing arm of its corona because of sensitivity limits and quality constraints on the data. Moreover, NGC 2547 has very similar half-mass and tidal radii, which leaves only Blanco 1 and the Pleiades with a clear excess of mass content inside their respective tidal radii (for Blanco 1, see also Zhang et al. 2020).

thumbnail Fig. 12.

Total mass content for each population. The statistics are displayed separately for the mass content inside and outside the tidal radius in light- and dark-shaded bars, respectively. For most clusters, the bulk of the mass content is located in their respective coronae.

4.5. Cluster kinematics

Our extensive membership determination also enables a more detailed look into the kinematic profile of each group. In particular, with the subset of sources that includes radial velocities in Gaia, we can determine 3D velocity vectors and subsequently characterize the expansion (or contraction) rate along different spatial axes. When the cluster samples are reduced to sources with all necessary measurements for calculating space motions, between about 10% and 20% of the full member list are retained, with a minimum of 30 sources in Messier 39 and a maximum of 242 sources for the Pleiades. Based on these subsets, we first constructed linear expansion profiles where we compared each spatial dimension (XYZ) with its corresponding velocities (UVW). A positive correlation between these components indicates expansion along the given spatial dimension, while a negative correlation indicates contracting motions.

The fitting procedure requires particular attention when linear relations between spatial and velocity components are constructed: the distribution of sources in the position-velocity planes (e.g., XU) does not follow a linear correlation where the dispersion along the axes can be described with a normal distribution. Especially outliers with small measurement errors can affect a weighted fitting procedure, which leads to heavily biased slope determinations. Similarly, introducing another parameter to account for the scatter (as done by, e.g., Wright et al. 2019) does not yield reliable results either because this additional scaling factor is different for all sources. For this reason, we decided to iteratively clip outliers prior to fitting the data by removing (a) sources with low signal-to-noise ratios (e.g., |U/σU|< 3) and (b) all sources outside a 3σ range about the median in each velocity component. During the subsequent fitting, we then ignored measurement errors and determined the significance of the slopes by bootstrapping with 500 iterations. The fits themselves were obtained with the astropy linear least-squares fitting engine.

Table 2 lists the resulting slopes kxu, kyv, and kzw and their errors for the XU, YV, and ZW planes, respectively. In addition, Fig. A.5 shows the fits in all position-velocity combinations for each population, with the data values displayed relative to the group medians on both the abscissa and ordinate. The figure shows a clear correlation between the error magnitude and the dominating measurement parameter for a given velocity dimension. For instance, if the line of sight of a cluster aligns with the direction of motion in Galactic rotation, the V component is mostly determined by radial velocity measurements. This is the case for Messier 39, where we find many V measurements below our significance threshold. Similarly, Blanco 1 is located toward the direction of the Galactic south pole. As a consequence, the measured radial velocities mostly affect the vertical velocity W, which again leads to many statistically insignificant measurements in this component.

The determined slopes reveal several interesting findings. We first note that the statistical errors of the fits are relatively large. Because of the frequently large scatter in the position-velocity diagrams, this is an expected effect and emphasizes the reliability of our bootstrapped fitting errors. Nevertheless, all populations, except for NGC 2547 and the Pleiades, show slopes different from 0 in at least one dimension at a 2σ significance level, six out of the ten clusters at a 3σ level. Moreover, most statistically relevant slopes are positive, revealing largely expanding configurations. This is particularly the case for the XU plane, which means that all groups except for the Pleiades and to some degree also NGC 2547 appear to expand along this dimension. For Messier 39 and IC 2602, we furthermore observe negative slopes in the ZW plane, indicating contracting motions along the Z-axis. For Messier 39, this measurement comes close to a 2σ significance, while it exceeds the 3σ level for IC 2602. This finding is of particular interest because both clusters show expanding motions along other spatial axes at the same time, revealing the highly dynamic states of the populations. Interestingly, we do not measure any significant expanding or contracting motions for the Pleiades: of our ten target populations, the fits for the Pleiades sources altogether result in values that are largely incompatible with expanding or contracting motions at all.

Because the cluster coronae stretch across hundreds of parsec, it is likely that the tidal field of the Galaxy, and specifically, the differential rotation pattern, affect the corona shapes and their dynamics. To investigate the nature of the expansion signature in reference to Galactic rotation, we further explored the relation between source kinematics and their position in the Galaxy. For this comparison, we chose to use angular velocities about the Galactic center (denoted as ω) because of the large size of the cluster populations. In this particular case, angular velocities have a distinct advantage over space motions because sources with similar space velocities in the direction of Galactic rotation, but different galactocentric radii, will nonetheless drift apart over time as a result of different angular velocities.

For an effective comparison of our measurements, we computed the differential angular velocity gradient (dω/dR) of the clusters about the Galactic center and compared the results to the Milky Way rotation curve. These gradients were computed through linear fits of the angular velocities as a function of distance to the Galactic center (Fig. A.4). Figure 13 displays the results of this procedure. We show the measured gradients in angular velocities for each cluster as a function of the Galactocentric radius. The indicated error in the ordinate is the statistical (1σ) error of the fit, while the error bar on the abscissa denotes the physical extent of each group. All populations, except for IC 2391, show an anticorrelation, in the sense of decreasing angular velocity with increasing Galactocentric radius. In the same figure, we show the gradient resulting from the differential rotation pattern of the Galaxy as derived from the rotation curve of the MWPotential2014 potential (Bovy 2015). In this comparison, all clusters, except for IC 2391, clearly appear to have an angular velocity gradient significantly different from 0 that is either compatible with or larger than the intrinsic differential velocity field of the Galaxy. Interestingly, the clusters with values that exceed the differential rotation of the Galaxy itself are the oldest in our sample: Messier 39, Platais 9, and NGC 2516.

thumbnail Fig. 13.

Measured angular velocity gradients compared to the Milky Way rotation curve. Most populations show a significant gradient along their galactocentric radial axes that is largely compatible with or larger than the intrinsic velocity field of the Galaxy.

We furthermore explored the internal velocity distribution of each cluster to investigate the isotropy of the velocity field and to identify preferred directions of motion. Figure 14 displays a position-angle histogram in the UV plane relative to the bulk motion of the cluster core. Each wedge is a histogram bin with a size of 15°. All clusters, with the exception of Messier 39 because of the small number of sources, feature a well-visible preferential axis of motion. For instance, Blanco 1 and Pleiades members mainly move along the X-axis, while stars in Platais 9 and NGC 2516 move mostly along the Y-axis. At the same time, the velocity distribution frequently appears to be asymmetric. To test the symmetry of the velocity field, we performed a Wilcoxon signed-rank test (Wilcoxon 1945) for each cluster and on each UVW velocity component. For all clusters, except for Blanco 1 and IC 2391, we find statistically significant (p <  0.05) evidence for asymmetry in at least one velocity component. Furthermore, a comparison to Fig. 8 shows that there is no correlation between the dominant internal direction of motion and the spatial orientation of a cluster.

thumbnail Fig. 14.

UV position angle histogram relative to the bulk motion of each cluster. Each wedge covers an angle of 15°, and its size and face color are proportional to the number of sources in the given direction (lighter colors and longer wedges denote more sources). Most clusters show a preferential symmetric axis of motion, but with an overall asymmetric internal velocity field.

We also attempted to determine the isotropy in radial or tangential velocity components. Specifically, the anisotropy parameter β, such as used by Baumgardt & Kroupa (2007) or in a different form by Vesperini et al. (2014), is a popular parameterization in star cluster simulations to constrain formation and evolution physics. However, because of several biases in our measurements, we did not arrive at reliable values for β. Specifically, low number statistics, individual outliers, and large differences in measurement errors for individual velocity components resulted in values and errors for β that did not allow any meaningful interpretation. For instance, using the definition of Baumgardt & Kroupa (2007), our measured β would have been compatible with both radial and tangential anisotropy within the statistical errors. On the other hand, the definition of β by Vesperini et al. (2014) includes measured velocity dispersions in each 3D velocity component. These measurements, however, are heavily affected by our perspective, where one component can be largely dominated by the measured heliocentric radial velocities, for instance, while the other components are mostly determined by the much better constrained proper motions. Conversely, determining velocity anisotropy from proper motions alone, as done by, for example, Wright et al. (2019), is not possible either because of the large on-sky extent of the coronae.

5. Discussion

5.1. Current cluster morphology

In Sect. 4.2 we showed that the compact cluster cores are surrounded by large comoving and co-eval coronae (see Fig. 8). To understand the observed morphology in general terms, we draw parallels to tidal structures of more evolved clusters and additionally take the observed kinematic profile into account. A particular characteristic pattern of the extended structures is their similar orientation in the XY plane. The leading arm of all clusters with prominent coronae is always oriented toward the Galactic center. In reference to tidal structures in older star clusters, stars on outbound Galactic orbits form a trailing arm because the angular velocity decreases about the Galactic center. Conversely, stars on inbound galactic orbits overtake the cluster, thus forming a leading arm (Kroupa 2008; Meingast & Alves 2019). Moreover, the measured angular velocity profiles agree well with this analogy, where we indeed find a matching gradient along the galactocentric radial axes (Figs. 13 and A.4). Consequently, stars that escape the gravitational potential of a cluster form highly elongated structures whose leading parts are tilted toward the Galactic center primarily as a consequence of Galactic differential rotation. Compared to this effect, the internal velocity field of each cluster can have a variety of preferred directions and appears to be negligible with respect to the shape of the formed coronae (Fig. 14).

Our results also illustrate that the characteristic tilt of the groups in the XY plane is not universal: especially NGC 2451A appears to be distinct from the other populations. Because we find populations younger than NGC 2451A, we conclude that the larger tilt angle of this population is not a consequence of age. Furthermore, based on our measurements, we also conclude that NGC 2451A does not show unusual values regarding cluster kinematics, such as Galactic orbit, galactocentric radial velocity, or velocity dispersion. This means that differential Galactic rotation cannot be the only determining factor for the different observed tilts and elongations. The current shape of NGC 2451A may therefore be a consequence of birth conditions, including size and orientation of or velocity gradients within the birth molecular cloud, for instance. Moreover, the angular velocity gradients of some clusters are far higher than the value attributable to the rotation curve. As a consequence, some clusters have intrinsic velocity fields that amplify this signal so that their expansion rates exceed the inherent rate enforced by the Milky Way.

In general, we do not see a direct correlation between the age of a population and its current spatial arrangement: even in the youngest clusters, we observe a variety of shapes and sizes. For instance, IC 2391 and NGC 2451A already appear as heavily elongated structures. In contrast, the only marginally younger cluster NGC 2547, with an age of about 30 Myr the youngest cluster in our sample, features a less pronounced corona. Taking this argument further, the 300 Myr cluster Messier 39 is surrounded by a smaller coronal structure than some of its counterparts in the sample that are about five to ten times younger. This finding is not surprising because environmental and physical conditions at birth and during the earliest stages of the cluster evolution are important.

We also see significant overlap between the populations, where particularly IC 2602, IC 2391, and NGC 2451A form an intertwined network. Because we already observe such a significant overlap with only ten analyzed clusters and conservative selection criteria, we can only speculate about the complex 3D picture that all populations in this volume must portray. Our findings also offer a new perspective on our earlier results published in this paper series, where Fürnkranz et al. (2019) also reported evidence for temporal overlapping stellar populations. Because of the rigid assumption of impact parameters of only 20 pc of the authors, the conclusion on the encounter timescale of about 250 Myr may be severely overestimated, and cohabitation events may instead occur much more often, perhaps even most of the time.

5.2. Origin of the coronae

The current appearance of the cluster superstructures in both spatial and kinematic aspects appears to be largely affected by differential rotation in the Galaxy. The emergence of the coronae in the first place, however, can be attributed to different mechanisms. In one scenario, the stars that we observe today in the coronae were initially part of a single, much more compact cluster and only became unbound after formation. In this case, the stars in the coronae gradually left the cluster core regions to steadily contribute to the growth of the tidal structure. In another less well-studied scenario, the initial configuration could have resembled a more complex configuration, where star formation occurred in a larger hub and the currently observed cluster core was already decoupled from other star-forming parts in the birth molecular cloud. If the velocity gradients in such a hub were small enough, the entire structure could have stayed together for several dozen million years and eventually have formed the large coronae that we observe today.

In the scenario of initially compact clusters, important factors are the initial configuration, and according to theory, the phase of residual gas expulsion and subsequent violent relaxation. In this case, a testable prediction refers to the spatial configuration of the systems, as well as to the observed velocity structure. Specifically, if the observed coronae are indeed a consequence of the dynamical evolution of initially compact clusters, the (remnant) cluster cores would indeed be found preferentially centered on the superstructures, and the surrounding coronae would appear to be largely symmetric (ignoring perturbations since their formation). Similarly, the velocity field should appear largely symmetric about the cluster core, indicating that stars leave the cluster along a preferred spatial axis. While we do find that the cluster core is preferentially located at the center of the coronae (see Fig. 8), we find mostly asymmetric velocity distributions about the cluster cores (Fig. 14). However, these measurements could be biased by our workflow. In our method, we searched for member stars that are kinematically compatible specifically with the cluster cores. As a consequence, we largely enforce symmetry in tangential velocity space, which is certainly mirrored in 3D velocities to some degree. Because of this potential bias with respect to symmetry, our measurements of clearly asymmetric velocity fields gain even more emphasis. At the same time, however, it is not clear if this bias could be systematically carried over to the spatial domain because the spatial domain on its own should not be affected by this problem as our application of DBSCAN in XYZ coordinates has no prior information on cluster locations or their shapes.

Prominent examples of detailed numerical investigations in reference to the long-term evolution of star clusters can be found in Baumgardt & Kroupa (2007), the paper series by Bekdaulet Shukirgaliyev (Shukirgaliyev et al. 2017, 2018, 2019), and Dinnbier & Kroupa (2020). Morphologically, the models in Shukirgaliyev et al. (2018) and Dinnbier & Kroupa (2020) are indeed reminiscent of the coronae presented in this paper. However, depending on the initial configuration (e.g., SFE, the ratio of rh to rt, or the dynamical state before gas expulsion), the models allow for a wide range of observed conditions at different cluster ages. As a consequence of our lack of knowledge about the initial conditions, as well as of the large variety of possible outcomes (and under the assumption that this scenario is indeed true), we are currently not in a position to present unambiguous evidence that this scenario did indeed occur.

A distinctive example can be made of the Pleiades, which we observe to be very compact; more than 80% of its mass lies within the tidal radius of the cluster. With our methods we are able to recover sources down to volume number densities of about 10−3pc−3 (see Fig. A.3). Based on this already excellent recovery threshold, the current observed cluster shape is roughly compatible with multiple setups as reported by Dinnbier & Kroupa (2020). Depending on the gas-expulsion timescale and the compactness of the initial cluster, models with SFE  = 1/3, 2/3, and even without any primordial gas (SFE  = 1) could resemble our determined cluster shape for the Pleiades. For this particular case, an optimization of our member identification to further improve the recovery rate seems to be necessary to draw meaningful conclusions. An improvement by another order of magnitude down to volume densities 10−4pc−3 could reveal more details and facilitate a better comparison to the Dinnbier & Kroupa (2020) models. To reach such limits, however, it will likely not be enough to improve the precision of the Gaia measurements. A significant boost in the recovery rate would be achieved by including additional features in the parameter space, such as a reliably measured third velocity component, or photometric information and population ages. In particular, 4MOST (de Jong et al. 2012), which will be installed at the VISTA telescope in the near future, will likely provide radial velocities for a significant fraction of cluster members.

In another formation scenario, the observed superstructures could be remnants of initially larger configurations. In such a setting, the observed cluster cores did not form in isolation, but instead were part of a larger complex with multiple simultaneous star-forming events. As a consequence, the observed structures could have formed in an already extended arrangement that is only amplified by the differential Galactic rotation. From an observational point of view, recent discoveries of young stellar structures on scales of hundreds of parsec provide an interesting reference point. For instance, Bouy & Alves (2015) explored large stream-like configurations of OB stars that seem to connect many very young clusters and OB associations in the solar neighborhood. Embedded in these so-called blue streams, Beccari et al. (2020) identified a 35 Myr old relic stellar filament that seems to bridge the gaps between several known and newly identified clusters (e.g., NGC 2547, NGC 2451B, Collinder 132, and Collinder 140). Because of the extent and age of the structure, the authors argued that it is most likely a remnant of a star-forming filament and not a consequence of tidal forces (for another example, see Jerabkova et al. 2019; Tian 2020). This observational evidence highlights that the observed coronae could be a consequence of the filamentary, highly substructured nature of star-forming regions. In reference to Beccari et al. (2020), we note that NGC 2547 is part of their identified structure (see also Cantat-Gaudin et al. 2019a,b). Our analysis of this cluster, however, does not produce unambiguous evidence that the cluster is embedded in a larger structure. NGC 2547 does appear kinematically well distinct from the larger complex that includes NGC 2451B, Collinder 132, and Collinder 140: Comparing the mean 3D velocity of NGC 2547 to NGC 2451B, we find similar vϕ and vZ components, but a difference in vR of about 5 km s−1.

From a theoretical point of view, Dinnbier & Kroupa (2020) also considered the scenario of star formation in a larger complex, but in a largely idealized setting. Their results indicate that in such a case, the stars that formed distributed throughout the molecular cloud would indeed form an elongated structure with sizes up to several hundred parsec in diameter. However, the bulk of these stars would be located below our volume density recovery threshold. Moreover, in their experiment, they focused on the Pleiades, and the above-described conditions only apply to very specific circumstances and furthermore are only presented for a population age of 125 Myr. At earlier stages in the cluster evolution, it is very likely that some imprint of the initial configuration is better preserved and may very well be visible within our constraints.

6. Conclusions and summary

We have analyzed ten nearby, prominent young open star clusters with Gaia DR2: α Per, Blanco 1, IC 2602, IC 2391, Messier 39, NGC 2451A, NGC 2516, NGC 2547, Platais 9, and the Pleiades (Fig. 1). We started our investigation by demonstrating the challenges of reliably identifying comoving stellar populations in proper motion space in large regions of the sky (Figs. 2 and 3). Modern machine-learning tools offer outstanding and easily deployed possibilities for analyzing the wealth of data coming from Gaia. However, consistency checks and reliable contamination estimates must be included to allow for physically meaningful interpretations. We introduced a novel cluster member identification tool that bypasses the intricacies coming from projection effects in proper motion space (Fig. 4). Following the member identification, we deconvolved the spatial distribution with a mixture of Gaussians to mitigate the effect of measurement errors in geometric distances (Fig. 5). This method will be made publicly available in a follow-up paper.

The main results of this work can be summarized as follows:

  1. We found that except for the Pleiades, all target clusters are surrounded by a vast stellar halo that we refer to as the corona. These coronae, together with the cluster cores, form extended cluster populations with sizes ≳100 pc (Figs. 8 and A.1). The coronae are comoving with the cluster cores and appear coeval in the observational HRDs.

  2. The spatial appearance of the coronae is reminiscent of tidal structures found in older clusters whose leading arm is oriented toward the inner Galaxy and that have an approximately symmetric trailing arm. This similarity is further confirmed by comparing the angular velocities of the cluster members about the Galactic center to the rotation curve of the Milky Way. All clusters, except for IC 2391, show angular velocity gradients across the coronae that are compatible with or exceed the intrinsic differential rotation pattern of the Galaxy (Fig. 13).

  3. The extracted cluster populations present a remarkably clean main sequence with stellar multiple sequences for most of the clusters and have 3D velocity dispersions typical of young populations (median 1.4 km s−1; Figs. 6 and 7). We determine a contamination rate in the extracted populations of a few percent (median 3%).

  4. When we compared our selection to the Galactic stellar field, our method recovered cluster members down to volume densities between two and three orders of magnitude below the background stellar field (Fig. 11). This corresponds to about 10−3pc−3 in volume number density (Fig. A.3).

  5. We constructed mass functions of all extended structures (core and corona) (Fig. 9). The total determined system masses enabled a robust calculation of the tidal and half-mass radii of the clusters (Fig. 10). The radial cluster profiles reveal that the bulk of the stellar mass for many of the clusters is located beyond their tidal radius (Fig 12).

  6. A concise analysis of the velocity field of each cluster reveals a highly dynamic picture. All clusters, with the exception of NGC 2547 and the Pleiades, show statistically significant expansion along at least one spatial axis. For Messier 39 and IC 2602, we find evidence for simultaneous expansion along their Galactic X (approximately galactocentric radial) axis and contraction along their Galactic vertical axis. The fits are presented in Fig. A.5.

  7. The velocity field of the extended populations is aligned along a spatial axis that is unique to each cluster. The general direction of motion with respect to the cluster cores appears to be asymmetric for all clusters, except for Blanco 1 and IC 2391 (Fig. 14).

The clusters we presented are relatively young and span ages across an order of magnitude from about 30 to 300 Myr. We argue that it is most plausible that our sample of cluster cores and coronae is the outcome of a combination of different effects. The factors that contribute to their current appearance likely include (a) the complex morphology and kinematics of star-forming molecular clouds, (b) the ability of young stars to preserve their natal kinematic properties, (c) the consequences of residual gas expulsion and subsequent violent relaxation, and (d) differential rotation in the Galactic disk.

For the youngest populations (NGC 2451A, IC 2602, IC 2391, and NGC 2547), stars that formed throughout the entire parental molecular cloud complex would still comprise large parts of the coronae in the form of relic structures. For the more evolved clusters with ages of about or above 100 Myr (Platais 9, Messier 39, α Per, Blanco 1, NGC 2516, and the Pleiades), this natal imprint could be beyond our current detection abilities, and the observed extended structures would be largely made up of stripped cluster stars.

Finally, our application of a new method to ten prominent nearby open star clusters revealed an unexpected and exciting new perspective on their spatial and dynamical structure. Our findings expose more questions than answers and open several avenues for follow-up research to extend our knowledge on the evolution of stellar structures in the Galaxy.


1

Throughout this paper, the term cluster will be used interchangeably with other expressions (e.g., system, structure) to refer to an entire comoving and coeval population, regardless of the gravitational boundedness.

3

The term corona has been used several times throughout the past decades in reference to star clusters, loosely referring to the less dense outer regions of open clusters (Artyukhina & Kholopov 1964). Here, we use the term corona to generally describe the stellar population of star clusters beyond their tidal radius.

4

An interactive 3D visualization of the spatial arrangement as visible in Fig. 8 of the populations is available online, also including a comparison to the catalogs from Cantat-Gaudin et al. (2018) and Kounkel & Covey (2019). In addition, Fig. A.1 shows the on-sky distribution of the populations.

Acknowledgments

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 research made use of Astropy, (http://www.astropy.org) a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018). This research has made use of “Aladin sky atlas” developed at CDS, Strasbourg Observatory, France (Bonnarel et al. 2000). We also acknowledge the various Python packages that were used in the data analysis of this work, including NumPy (van der Walt et al. 2011), SciPy (Jones et al. 2001), scikit-learn (Pedregosa et al. 2011), Matplotlib (Hunter 2007), and Plotly (Plotly Technologies Inc. 2015). This research has made use of the SIMBAD database operated at CDS, Strasbourg, France (Wenger et al. 2000). This research has made use of the VizieR catalogue access tool, CDS, Strasbourg, France (Ochsenbein et al. 2000; DOI: 10.26093/cds/vizier).

References

  1. Alves, J., Zucker, C., Goodman, A. A., et al. 2020, Nature, 578, 237 [Google Scholar]
  2. Ankerst, M., Breunig, M. M., Kriegel, H. P., & Sander, J. 1999, Proceedings of the 1999 ACM SIGMOD International Conference on Management of Data, SIGMOD ’99 (New York, NY, USA: Association for Computing Machinery), 49 [Google Scholar]
  3. Arenou, F., Luri, X., Babusiaux, C., et al. 2018, A&A, 616, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Artyukhina, N. M., & Kholopov, P. N. 1964, Sov. Astron., 7, 840 [Google Scholar]
  5. Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123 [Google Scholar]
  7. Bailer-Jones, C. A. L. 2015, PASP, 127, 994 [NASA ADS] [CrossRef] [Google Scholar]
  8. Barrado y Navascués, D., Stauffer, J. R., & Patten, B. M. 1999, ApJ, 522, L53 [NASA ADS] [CrossRef] [Google Scholar]
  9. Baumgardt, H., & Kroupa, P. 2007, MNRAS, 380, 1589 [NASA ADS] [CrossRef] [Google Scholar]
  10. Beccari, G., Boffin, H. M. J., & Jerabkova, T. 2020, MNRAS, 491, 2205 [CrossRef] [Google Scholar]
  11. Bennett, M., & Bovy, J. 2019, MNRAS, 482, 1417 [NASA ADS] [CrossRef] [Google Scholar]
  12. Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton University Press) [Google Scholar]
  13. Bonnarel, F., Fernique, P., Bienaymé, O., et al. 2000, A&AS, 143, 33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Boss, L. J. 1908, AJ, 26, 31 [NASA ADS] [CrossRef] [Google Scholar]
  15. Bossini, D., Vallenari, A., Bragaglia, A., et al. 2019, A&A, 623, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Bouy, H., & Alves, J. 2015, A&A, 584, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Bovy, J. 2015, ApJS, 216, 29 [Google Scholar]
  18. Bovy, J., Hogg, D. W., & Roweis, S. T. 2011, Ann. Appl. Stat., 5, 1657 [NASA ADS] [CrossRef] [Google Scholar]
  19. Bressert, E., Bastian, N., Gutermuth, R., et al. 2010, MNRAS, 409, L54 [NASA ADS] [Google Scholar]
  20. Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127 [Google Scholar]
  21. Brinkmann, N., Banerjee, S., Motwani, B., & Kroupa, P. 2017, A&A, 600, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Campello, R. J. G. B., Moulavi, D., & Sander, J. 2013, in Advances in Knowledge Discovery and Data Mining, eds. J. Pei, V. S. Tseng, L. Cao, H. Motoda, & G. Xu (Berlin, Heidelberg: Springer, Berlin Heidelberg), 160 [Google Scholar]
  23. Cantat-Gaudin, T., Jordi, C., Vallenari, A., et al. 2018, A&A, 618, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Cantat-Gaudin, T., Jordi, C., Wright, N. J., et al. 2019a, A&A, 626, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Cantat-Gaudin, T., Mapelli, M., Balaguer-Núñez, L., et al. 2019b, A&A, 621, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Cartwright, A., & Whitworth, A. P. 2004, MNRAS, 348, 589 [NASA ADS] [CrossRef] [Google Scholar]
  27. Castro-Ginard, A., Jordi, C., Luri, X., Cantat-Gaudin, T., & Balaguer-Núñez, L. 2019, A&A, 627, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Castro-Ginard, A., Jordi, C., Luri, X., et al. 2020, A&A, 635, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Chumak, Y. O., Platais, I., McLaughlin, D. E., Rastorguev, A. S., & Chumak, O. V. 2010, MNRAS, 402, 1841 [NASA ADS] [CrossRef] [Google Scholar]
  30. Cummings, J. D., & Kalirai, J. S. 2018, ApJ, 156, 165 [NASA ADS] [CrossRef] [Google Scholar]
  31. Curtis, J. L., Agüeros, M. A., Mamajek, E. E., Wright, J. T., & Cummings, J. D. 2019, AJ, 158, 77 [NASA ADS] [CrossRef] [Google Scholar]
  32. Dahm, S. E. 2015, ApJ, 813, 108 [NASA ADS] [CrossRef] [Google Scholar]
  33. Dale, J. E., Ercolano, B., & Bonnell, I. A. 2015, MNRAS, 451, 987 [NASA ADS] [CrossRef] [Google Scholar]
  34. Dalessandro, E., Miocchi, P., Carraro, G., Jílková, L., & Moitinho, A. 2015, MNRAS, 449, 1811 [NASA ADS] [CrossRef] [Google Scholar]
  35. de Jong, R. S., Bellido-Tirado, O., Chiappini, C., et al. 2012, in Ground-based and Airborne Instrumentation for Astronomy IV, SPIE Conf. Ser., 8446, 84460T [CrossRef] [Google Scholar]
  36. De Silva, G. M., D’Orazi, V., Melo, C., et al. 2013, MNRAS, 431, 1005 [NASA ADS] [CrossRef] [Google Scholar]
  37. Dinnbier, F., & Kroupa, P. 2020, A&A, 640, A85 [EDP Sciences] [Google Scholar]
  38. Dobbie, P. D., Lodieu, N., & Sharp, R. G. 2010, MNRAS, 409, 1002 [NASA ADS] [CrossRef] [Google Scholar]
  39. Drimmel, R., & Poggio, E. 2018, Res. Notes Am. Astron. Soc., 2, 210 [NASA ADS] [CrossRef] [Google Scholar]
  40. Eggen, O. J. 1987, in NATO Advanced Science Institutes (ASI) Series C, eds. G. Gilmore, & B. Carswell, NATO ASI Ser. C, 207, 211 [Google Scholar]
  41. Elmegreen, B. G. 2008, ApJ, 672, 1006 [NASA ADS] [CrossRef] [Google Scholar]
  42. Ernst, A., Just, A., Berczik, P., & Petrov, M. I. 2010, A&A, 524, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. Ernst, A., Just, A., Berczik, P., & Olczak, C. 2011, A&A, 536, A64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Ernst, A., Berczik, P., Just, A., & Noel, T. 2015, Astron. Nachr., 336, 577 [NASA ADS] [CrossRef] [Google Scholar]
  45. Ester, M., Kriegel, H. P., Sander, J., & Xu, X. 1996, Proc. of 2nd International Conference on Knowledge Discovery, 226 [Google Scholar]
  46. Evans, D. W., Riello, M., De Angeli, F., et al. 2018, A&A, 616, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
  48. Foreman-Mackey, D., Farr, W., Sinha, M., et al. 2019, J. Open Source Software, 4, 1864 [NASA ADS] [CrossRef] [Google Scholar]
  49. Fürnkranz, V., Meingast, S., & Alves, J. 2019, A&A, 624, L11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  50. Gaia Collaboration (Prusti, T., et al.) 2016, A&A, 595, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Goodwin, S. P. 2009, Ap&SS, 324, 259 [NASA ADS] [CrossRef] [Google Scholar]
  53. Gravity Collaboration (Abuter, R., et al.) 2018, 615, L15 [Google Scholar]
  54. Großschedl, J. E., Alves, J., Meingast, S., et al. 2018, A&A, 619, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Gutermuth, R. A., Myers, P. C., Megeath, S. T., et al. 2008, ApJ, 674, 336 [NASA ADS] [CrossRef] [Google Scholar]
  56. Hacar, A., Tafalla, M., Forbrich, J., et al. 2018, A&A, 610, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Houk, N., & Cowley, A. P. 1975, University of Michigan Catalogue of Two-dimensional Spectral Types for the HD Stars, I (Department of Astronomy: University of Michigan) [Google Scholar]
  58. Hünsch, M., Weidner, C., & Schmitt, J. H. M. M. 2003, A&A, 402, 571 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  59. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
  60. Jeffries, R. D., & Oliveira, J. M. 2005, MNRAS, 358, 13 [NASA ADS] [CrossRef] [Google Scholar]
  61. Jeffries, R. D., James, D. J., & Thurston, M. R. 1998, MNRAS, 300, 550 [NASA ADS] [CrossRef] [Google Scholar]
  62. Jerabkova, T., Boffin, H. M. J., Beccari, G., & Anderson, R. I. 2019, MNRAS, 489, 4418 [CrossRef] [Google Scholar]
  63. Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open source scientific tools for Python [Google Scholar]
  64. Juarez, A. J., Cargile, P. A., James, D. J., & Stassun, K. G. 2014, ApJ, 795, 143 [NASA ADS] [CrossRef] [Google Scholar]
  65. Kamdar, H., Conroy, C., Ting, Y.-S., et al. 2019, ApJ, 884, L42 [CrossRef] [Google Scholar]
  66. Keenan, P. C., & McNeil, R. C. 1989, ApJS, 71, 245 [NASA ADS] [CrossRef] [Google Scholar]
  67. Kharchenko, N. V., Piskunov, A. E., Röser, S., Schilbach, E., & Scholz, R. D. 2005, A&A, 438, 1163 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
  68. Kharchenko, N. V., Piskunov, A. E., Roeser, S., Schilbach, E., & Scholz, R. D. 2013, VizieR Online Data Catalog, J/A+A/558/A53 [Google Scholar]
  69. King, I. 1962, AJ, 67, 471 [NASA ADS] [CrossRef] [Google Scholar]
  70. Kounkel, M., & Covey, K. 2019, AJ, 158, 122 [Google Scholar]
  71. Krone-Martins, A., & Moitinho, A. 2014, A&A, 561, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. Kroupa, P. 2001, MNRAS, 322, 231 [NASA ADS] [CrossRef] [Google Scholar]
  73. Kroupa, P. 2008, Initial Conditions for Star Clusters (Dordrecht: Springer, Netherlands), 760, 181 [Google Scholar]
  74. Kruijssen, J. M. D. 2012, MNRAS, 426, 3008 [NASA ADS] [CrossRef] [Google Scholar]
  75. Kruijssen, J. M. D., Maschberger, T., Moeckel, N., et al. 2012, MNRAS, 419, 841 [NASA ADS] [CrossRef] [Google Scholar]
  76. Krumholz, M. R., McKee, C. F., & Bland -Hawthorn, J. 2019, ARA&A, 57, 227 [NASA ADS] [CrossRef] [Google Scholar]
  77. Kuhn, M. A., Feigelson, E. D., Getman, K. V., et al. 2014, ApJ, 787, 107 [Google Scholar]
  78. Kuhn, M. A., Hillenbrand, L. A., Sills, A., Feigelson, E. D., & Getman, K. V. 2019, ApJ, 870, 32 [Google Scholar]
  79. Küpper, A. H. W., Kroupa, P., Baumgardt, H., & Heggie, D. C. 2010a, MNRAS, 407, 2241 [NASA ADS] [CrossRef] [Google Scholar]
  80. Küpper, A. H. W., Kroupa, P., Baumgardt, H., & Heggie, D. C. 2010b, MNRAS, 401, 105 [NASA ADS] [CrossRef] [Google Scholar]
  81. Lada, C. J., & Lada, E. A. 2003, ARA&A, 41, 57 [NASA ADS] [CrossRef] [Google Scholar]
  82. Lada, C. J., Margulis, M., & Dearborn, D. 1984, ApJ, 285, 141 [NASA ADS] [CrossRef] [Google Scholar]
  83. Lindegren, L. 2018, GAIA-C3-TN-LU-LL-124 [Google Scholar]
  84. Lindegren, L., Hernández, J., Bombrun, A., et al. 2018, A&A, 616, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  85. Lyra, W., Moitinho, A., van der Bliek, N. S., & Alves, J. 2006, A&A, 453, 101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  86. McInnes, L., Healy, J., & Astels, S. 2017, J. Open Source Software, 2, 205 [CrossRef] [Google Scholar]
  87. Meingast, S., & Alves, J. 2019, A&A, 621, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  88. Meingast, S., Alves, J., & Fürnkranz, V. 2019, A&A, 622, L13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  89. Mermilliod, J. C. 1981, A&A, 97, 235 [NASA ADS] [Google Scholar]
  90. Meynet, G., Mermilliod, J. C., & Maeder, A. 1993, A&AS, 98, 477 [Google Scholar]
  91. Navascues, D. B., Stauffer, J. R., & Jayawardhana, R. 2004, ApJ, 614, 386 [NASA ADS] [CrossRef] [Google Scholar]
  92. Naylor, T., & Jeffries, R. D. 2006, MNRAS, 373, 1251 [NASA ADS] [CrossRef] [Google Scholar]
  93. Ochsenbein, F., Bauer, P., & Marcout, J. 2000, A&AS, 143, 23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  94. Parker, R. J., & Goodwin, S. P. 2007, MNRAS, 380, 1271 [NASA ADS] [CrossRef] [Google Scholar]
  95. Parker, R. J., & Meyer, M. R. 2012, MNRAS, 427, 637 [Google Scholar]
  96. Peacock, J. A. 1983, MNRAS, 202, 615 [NASA ADS] [Google Scholar]
  97. Pecaut, M. J., & Mamajek, E. E. 2013, ApJS, 208, 9 [Google Scholar]
  98. Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, J. Mach. Learn. Res., 12, 2825 [Google Scholar]
  99. Pellerin, A., Meyer, M., Harris, J., & Calzetti, D. 2007, ApJ, 658, L87 [NASA ADS] [CrossRef] [Google Scholar]
  100. Platais, I., Kozhurina-Platais, V., & van Leeuwen, F. 1998, AJ, 116, 2423 [NASA ADS] [CrossRef] [Google Scholar]
  101. Platais, I., Kozhurina-Platais, V., Barnes, S., et al. 2001, AJ, 122, 1486 [NASA ADS] [CrossRef] [Google Scholar]
  102. Plotly Technologies Inc. 2015, Collaborative Data Science [Google Scholar]
  103. Randich, S., Tognelli, E., Jackson, R., et al. 2018, A&A, 612, A99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  104. Ratzenböck, S., Meingast, S., Alves, J., Möller, T., & Bomze, I. 2020, A&A, 639, A64 [CrossRef] [EDP Sciences] [Google Scholar]
  105. Reid, M. J., & Brunthaler, A. 2004, ApJ, 616, 872 [NASA ADS] [CrossRef] [Google Scholar]
  106. Renaud, F., Gieles, M., & Boily, C. M. 2011, MNRAS, 418, 759 [NASA ADS] [CrossRef] [Google Scholar]
  107. Röser, S., & Schilbach, E. 2019, A&A, 627, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  108. Röser, S., & Schilbach, E. 2020, A&A, 638, A9 [CrossRef] [EDP Sciences] [Google Scholar]
  109. Schönrich, R., Binney, J., & Dehnen, W. 2010, MNRAS, 403, 1829 [NASA ADS] [CrossRef] [Google Scholar]
  110. Schwarz, G. 1978, Ann. Stat., 6, 461 [Google Scholar]
  111. Shukirgaliyev, B., Parmentier, G., Berczik, P., & Just, A. 2017, A&A, 605, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  112. Shukirgaliyev, B., Parmentier, G., Just, A., & Berczik, P. 2018, ApJ, 863, 171 [NASA ADS] [CrossRef] [Google Scholar]
  113. Shukirgaliyev, B., Parmentier, G., Berczik, P., & Just, A. 2019, MNRAS, 486, 1045 [CrossRef] [Google Scholar]
  114. Silaj, J., & Landstreet, J. D. 2014, A&A, 566, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  115. Sowell, J. R. 1987, ApJS, 64, 241 [NASA ADS] [CrossRef] [Google Scholar]
  116. Spitzer, L., Jr 1958, ApJ, 127, 17 [NASA ADS] [CrossRef] [Google Scholar]
  117. Stauffer, J. R., Hartmann, L. W., Prosser, C. F., et al. 1997, ApJ, 479, 776 [NASA ADS] [CrossRef] [Google Scholar]
  118. Stauffer, J. R., Schultz, G., & Kirkpatrick, J. D. 1998, ApJ, 499, L199 [NASA ADS] [CrossRef] [Google Scholar]
  119. Tadross, A., Werner, P., Osman, A., & Marie, M. 2002, New Ast., 7, 553 [NASA ADS] [CrossRef] [Google Scholar]
  120. Tian, H.-J. 2020, ApJ, 904, 196 [CrossRef] [Google Scholar]
  121. Vande Putte, D., Garnier, T. P., Ferreras, I., Mignani, R. P., & Cropper, M. 2010, MNRAS, 407, 2109 [NASA ADS] [CrossRef] [Google Scholar]
  122. van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Comput. Sci. Eng., 13, 22 [Google Scholar]
  123. van Leeuwen, F. 2009, A&A, 497, 209 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  124. Vesperini, E., Varri, A. L., McMillan, S. L. W., & Zepf, S. E. 2014, MNRAS, 443, L79 [NASA ADS] [CrossRef] [Google Scholar]
  125. Ward, J. L., Kruijssen, J. M. D., & Rix, H.-W. 2020, MNRAS, 495, 663 [CrossRef] [Google Scholar]
  126. Weiler, M. 2018, A&A, 617, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  127. Wenger, M., Ochsenbein, F., Egret, D., et al. 2000, A&AS, 143, 9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  128. Wilcoxon, F. 1945, Biom. Bull., 1, 80 [CrossRef] [Google Scholar]
  129. Wright, N. J., & Mamajek, E. E. 2018, MNRAS, 476, 381 [Google Scholar]
  130. Wright, N. J., Jeffries, R. D., Jackson, R. J., et al. 2019, MNRAS, 486, 2477 [CrossRef] [Google Scholar]
  131. Yen, S. X., Reffert, S., Schilbach, E., et al. 2018, A&A, 615, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  132. Yeh, F. C., Carraro, G., Montalto, M., & Seleznev, A. F. 2019, AJ, 157, 115 [NASA ADS] [CrossRef] [Google Scholar]
  133. Zhang, J., Zhao, J., Oswalt, T. D., et al. 2019, ApJ, 887, 84 [CrossRef] [Google Scholar]
  134. Zhang, Y., Tang, S.-Y., Chen, W. P., Pang, X., & Liu, J. Z. 2020, ApJ, 889, 99 [CrossRef] [Google Scholar]

Appendix A: Supplementary plots and tables

In this part of the appendix, we present a range of plots and tables that mainly provide supplementary information to the main article. Table A.1 contains the first ten sources of our selection for each cluster (the full table is available at the CDS). All raw Gaia data can be obtained by cross-matching either the Gaia DR2 source ID or the set of equatorial coordinates. In addition, we also list a set of variables specific to our analysis. These include the distance d that was used for the member identification (raw), as well as the finally adopted distance after extreme deconvolution (xd). We also list the Galactic Cartesian velocities UVW, and for each source, the contamination fraction fc and the name of the associated cluster. Table A.2 lists a set of unique age determinations for the ten clusters discussed in this manuscript. This list is not complete, but should indicate the large variance with respect to age determinations for the populations. Our adopted age in this manuscript can be found in Table 1.

Table A.1.

First ten sources sorted by source ID.

Table A.2.

Compilation of cluster ages from the literature.

Figure A.1 displays the selected sources in Galactic coordinates in a stereographic projection. Individual sources are displayed as colored circles. The sizes of the data points are proportional to the measured distance, where smaller circles represent more remote sources. This scaling is calculated individually for each panel and cannot be compared between the clusters. Figure A.2 displays for each cluster the source distribution in the Galactic Cartesian XY plane. The color scale for each point shows the calculated contamination fraction for each source. This measure increases with distance from the density maxima of the clusters. Figure A.3 displays the measured volume number density as a function of distance to the cluster center. The measured volume density for each source was parameterized with the distance to the seventh nearest neighbor d7NN so that (including the source itself, there are eight sources in the given volume spanned by the distance to the seventh nearest neighbor). For all populations, we can recover sources down to about 10−3 pc−3. Figures A.5 and A.4 provide the results of our linear fitting procedures for the XYZUVW linear expansion model and angular velocity gradient about the Galactic center, respectively. In both diagrams, all sources included in the fit are displayed as colored points, whereas sources that did not pass the applied 3σ significance threshold are marked as crosses.

thumbnail Fig. A.1.

On-sky view of clusters in Galactic coordinates in a stereographic projection. The luminance value of the symbol colors is proportional to volume density, with higher densities corresponding to darker shades. The symbol sizes are proportional to distance, where the size shrinks for increasing distances. The symbol sizes are scaled for each panel individually and are not comparable between clusters.

thumbnail Fig. A.2.

Source distribution in the XY plane (centered on the Sun) for each cluster, color-coded with the contamination fraction fc. We observe a strong correlation between fc and the distance from the cluster centers, resulting from much lower volume number densities in the coronae.

thumbnail Fig. A.3.

Volume number density profiles as a function of radius from the density maxima of the clusters. The volume density was parameterized with the seventh nearest neighbor of each source. The solid colored lines represent a smoothed running median that is plotted in addition to the values for individual sources. For reference, the tidal radius of each cluster is marked with a vertical gray dashed line.

thumbnail Fig. A.4.

Galactocentric angular velocity ω as a function of galactocentric radius for each cluster. The solid gray lines show 500 independent bootstrapped fitting results. The dashed lines are the adopted linear fits that were determined as the median of all realizations. The crosses are individual measurements below 3σ significance and were not included in the fit.

thumbnail Fig. A.5.

Linear expansion model fits. Each row represents an individual cluster. The columns from left to right show the relation and linear fits for the parameter combinations XU, YV, and ZW, respectively. All values are shown relative to the density maxima or bulk motions of the clusters. The crosses are individual measurements below 3σ significance and were not included in the fit.

Appendix B: Coordinate system definitions

The following list summarizes our parameter setup for the calculation of galactocentric cylindrical velocities:

  • αGC = 266.4051°; δGC = −28.936175°; right ascension and declination of the Galactic center (Reid & Brunthaler 2004).

  • R0 = 8122 pc; distance from the vantage point to the Galactic center (Abuter 2018).

  • ϕ0 = 180°; position angle of the Sun in galactocentric cylindrical coordinates.

  • Z0 = 20.8 pc; height of the Sun above the Galactic plane toward the Galactic north pole (Bennett & Bovy 2019).

  • vR, ⊙ = 12.9 km s−1, vϕ, ⊙ = 245.6 km s1, vZ, ⊙ = 7.78 km s−1; solar velocity in galactocentric cylindrical coordinates (Drimmel & Poggio 2018).

  • (U, V, W)solar = (11.1, 12.24, 7.25) km s−1; barycentric velocity of the Sun relative to the local standard of rest (Schönrich et al. 2010).

  • vR; galactocentric radial velocity component. Positive when moving away from the Galactic center.

  • vϕ; galactocentric azimuthal velocity component. Positive in the direction of Galactic rotation.

  • vZ; galactocentric vertical velocity component. Positive toward the Galactic north pole.

Supplementary Material

Interactive 3D version of Fig. 8 (Access here)

All Tables

Table 1.

Cluster parameters obtained from the literature.

Table 2.

Computed cluster parameters and statistics.

Table A.1.

First ten sources sorted by source ID.

Table A.2.

Compilation of cluster ages from the literature.

All Figures

thumbnail Fig. 1.

Postage stamps for our ten target clusters, extracted from the Digitized Sky Survey; each stamp has a side length of 10 pc at the distance of the cluster core, and is displayed so that Galactic north is up, and east is left. From top left to bottom right, the clusters are sorted by their peak volume density.

In the text
thumbnail Fig. 2.

Tangential velocity signature for sources comoving with the cluster IC 2391 in Galactic projection. The top panel displays velocities along Galactic longitude (vl), and the bottom panel shows velocities along Galactic latitude (vb). The position of the cluster core is marked with a black star. The tangential velocities depict a complex, highly variable pattern and illustrate the challenge of identifying comoving populations in large regions of the sky.

In the text
thumbnail Fig. 3.

Cumulative histograms for the deconvolved 3D velocity dispersions (left panel) and apparent sizes (right panel) of identified structures as published by Cantat-Gaudin et al. (2018; red) and Kounkel & Covey (2019; blue). Cantat-Gaudin et al. characterized members close to known clusters with relatively small velocity dispersions. Kounkel & Covey identified much larger structures, but more than 90% of the populations feature velocity dispersions greater than 5 km s−1. Young structures (< 100 Myr) in the sample from Kounkel & Covey (2019; dashed blue line) also feature mostly relatively large velocity dispersions.

In the text
thumbnail Fig. 4.

Member selection process for the α Per cluster. The left and center columns display Galactic XYZ coordinates centered on the Sun, and the right-hand side column shows deprojected tangential velocities. The top row shows kernel density maps centered on α Per that distinctly highlight the cluster populations. The line of sight (dashed line) and the Galactic orbit of α Per (solid line) are also shown. The bottom panels present the member selection procedure, in which first a radial filter is applied in deprojected tangential velocities (red sources), followed by DBSCAN clustering in the XYZ parameter space to produce the final selection (blue dots).

In the text
thumbnail Fig. 5.

Spatial configuration of the central regions of NGC 2516 before (left) and after (right) deconvolving source distances. The original distribution clearly shows a substantial elongation effect along the line of sight (dotted gray line), while the extreme deconvolution (XD) infers the expected spherical arrangement of the cluster core.

In the text
thumbnail Fig. 6.

Cumulative distributions of the source contamination fractions fc for each cluster. The contamination fraction fc denotes the volume density contrast for each source with respect to the surrounding background population.

In the text
thumbnail Fig. 7.

Observational HRDs for the ten target clusters. Members of each cluster are displayed as colored points. The gray distribution in the background comprises all Gaia sources in the initial search box of each population. Our selection results in remarkably clean and narrow main sequences that frequently also show well-visible binary sequences. This clear separation from the broad background distribution acts as validation of the successful application of our member selection procedure.

In the text
thumbnail Fig. 8.

Spatial distribution of our selection for the ten clusters in heliocentric Galactic coordinates. The color scale is proportional to volume density, accentuating the dense central regions of the populations. In this view, the Sun is located at the origin, and the directions toward the Galactic center (positive X) and of the generic Galactic rotation (positive Y) are indicated with arrows. Cluster cores are surrounded by enormous stellar coronae with sizes up to several hundred parsec. In addition, we observe a characteristic pattern where extended cluster populations form elongated shapes whose leading arm is tilted toward the inner Galaxy. An interactive 3D version of this figure is also available online.

In the text
thumbnail Fig. 9.

Present-day mass functions (colored histograms) and the best-fitting truncated IMFs (solid black histogram) for each population. The completeness limits on the low-mass end appear different for each cluster because of the distance-dependent sensitivity limits. The completeness at the high-mass end is more difficult to evaluate because individual stellar evolutionary histories are governed by the independent ages of the groups. If the mass distribution follows the displayed truncated standard IMF, we estimate the missing mass fraction to be in a range from about 15% for IC 2602 to almost 45% for NGC 2516.

In the text
thumbnail Fig. 10.

Radial mass profiles and tidal radii for each cluster population. The colored lines depict the enclosed mass fraction as a function of distance from the density maximum of each population. The solid black lines trace the tidal radius for a given enclosed total mass content. The intersection of these profiles identifies the tidal radius. The mass profiles themselves depict a variety of shapes, with some clusters showing clear kinks near the tidal radius. For the majority of the clusters, the mass profiles also reveal that the bulk of their stellar mass content is located beyond the tidal radii in the stellar coronae.

In the text
thumbnail Fig. 11.

Population volume densities normalized to the surrounding field as a function of distance from the cluster center, parameterized by the tidal radii. The observed volume density vanishes into the field at the tidal radius. Our member selection recovers sources across five orders of magnitude in volume density, down to three orders of magnitude below the Galactic field at the outermost rim of the coronae.

In the text
thumbnail Fig. 12.

Total mass content for each population. The statistics are displayed separately for the mass content inside and outside the tidal radius in light- and dark-shaded bars, respectively. For most clusters, the bulk of the mass content is located in their respective coronae.

In the text
thumbnail Fig. 13.

Measured angular velocity gradients compared to the Milky Way rotation curve. Most populations show a significant gradient along their galactocentric radial axes that is largely compatible with or larger than the intrinsic velocity field of the Galaxy.

In the text
thumbnail Fig. 14.

UV position angle histogram relative to the bulk motion of each cluster. Each wedge covers an angle of 15°, and its size and face color are proportional to the number of sources in the given direction (lighter colors and longer wedges denote more sources). Most clusters show a preferential symmetric axis of motion, but with an overall asymmetric internal velocity field.

In the text
thumbnail Fig. A.1.

On-sky view of clusters in Galactic coordinates in a stereographic projection. The luminance value of the symbol colors is proportional to volume density, with higher densities corresponding to darker shades. The symbol sizes are proportional to distance, where the size shrinks for increasing distances. The symbol sizes are scaled for each panel individually and are not comparable between clusters.

In the text
thumbnail Fig. A.2.

Source distribution in the XY plane (centered on the Sun) for each cluster, color-coded with the contamination fraction fc. We observe a strong correlation between fc and the distance from the cluster centers, resulting from much lower volume number densities in the coronae.

In the text
thumbnail Fig. A.3.

Volume number density profiles as a function of radius from the density maxima of the clusters. The volume density was parameterized with the seventh nearest neighbor of each source. The solid colored lines represent a smoothed running median that is plotted in addition to the values for individual sources. For reference, the tidal radius of each cluster is marked with a vertical gray dashed line.

In the text
thumbnail Fig. A.4.

Galactocentric angular velocity ω as a function of galactocentric radius for each cluster. The solid gray lines show 500 independent bootstrapped fitting results. The dashed lines are the adopted linear fits that were determined as the median of all realizations. The crosses are individual measurements below 3σ significance and were not included in the fit.

In the text
thumbnail Fig. A.5.

Linear expansion model fits. Each row represents an individual cluster. The columns from left to right show the relation and linear fits for the parameter combinations XU, YV, and ZW, respectively. All values are shown relative to the density maxima or bulk motions of the clusters. The crosses are individual measurements below 3σ significance and were not included in the fit.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.