© The Authors 2023

1. Introduction

Cluster formation is important for the study of many key questions in modern astrophysics. First, it is central to the star formation process itself (e.g., McKee & Ostriker 2007). It is commonly accepted that most stars in galaxies (70%−90%) form in groups, clusters, or hierarchical systems, and spend some time gravitationally bound with their siblings when still embedded in their progenitor molecular cloud (Lada & Lada 2003). In the Milky Way (MW) this evidence comes from the global clustered structure of the disk and spiral arms (e.g., Kounkel & Covey 2019) and from the similar star formation rates observed in embedded clusters (Lada & Lada 2003) and in the field (Miller & Scalo 1979). On a larger scale, this indication is supported by the good agreement between the mass density in stellar clusters and the average co-moving stellar density at the peak of the universal star formation density at redshift ∼2 (Madau & Dickinson 2014).

Secondly, cluster formation has many implications for the early interplay between stellar and gas dynamics, the possible formation of gravitational wave sources (Di Carlo et al. 2019; Banerjee 2021), and the dynamical properties of young star clusters (e.g., McMillan et al. 2007; Ballone et al. 2020, 2021; Livernois et al. 2021; Tiongco et al. 2022). Finally, cluster formation is fundamental for our understanding of the assembly process of galaxies in a cosmological context, as major star-forming episodes in galaxies are typically accompanied by significant star cluster production (Forbes et al. 2018) and the main properties of these systems are thus strictly connected with those of their hosts (Brodie & Strader 2006; Dalessandro et al. 2012).

However, the possible presence of unifying principles governing the formation of stellar clusters and whether they form through a monolithic event or as a result of a hierarchical process of star formation in which stars are formed across a continuous distribution of gas densities is still a matter of debate (Lada et al. 1984; Moeckel & Bonnell 2009; Kruijssen 2012; Banerjee & Kroupa 2014, 2015; Treviño-Morales et al. 2019; Kuhn et al. 2019; Pang et al. 2022). In addition, the observed presence in the most massive clusters of the multiple stellar populations characterized by specific abundance patterns in a number of light elements (see, e.g., Bastian & Lardo 2018; Gratton et al. 2019 for recent reviews) has raised new questions on the physical mechanisms at the basis of cluster formation, and their dependence on the environment and the formation epoch (Krumholz et al. 2019). Despite the tremendous observational and theoretical efforts in the recent years (Allison et al. 2010; Parker 2014; Adamo et al. 2015), our understanding of how and where star clusters form is still in its infancy.

Numerous recent observational studies of stellar clusters have significantly enriched our knowledge about the properties of these systems. The unprecedented kinematical mapping of the Galaxy and its stellar components, as secured by the Gaia mission (Gaia Collaboration 2021), is revolutionizing the field (Cantat-Gaudin & Anders 2020; Cantat-Gaudin et al. 2020; Castro-Ginard et al. 2022; Cantat-Gaudin 2022) enabling detailed studies of nearby star-forming regions and their use as ideal laboratories to shed new light on our understanding of cluster formation and early evolution (e.g., Beccari et al. 2018; Meingast et al. 2019; Lim et al. 2020). Observations have revealed significant structural and kinematical complexity, such as significant deviations from spherical symmetry, the presence of extended tails in both young (e.g., Meingast et al. 2021; Jerabkova et al. 2021; Pang et al. 2021) and old star clusters (e.g., Grillmair 2019; Bonaca et al. 2020), as well as evidence of internal rotation and/or radially anisotropic velocity distributions (e.g., Hénault-Brunet et al. 2012; Ferraro et al. 2018; Kamann et al. 2018; Vasiliev & Baumgardt 2021; Dalessandro et al. 2021a).

It has also been shown that young stellar systems often have a complex clumpy structure characterized by the presence of several stellar subsystems (e.g., Kuhn et al. 2019, 2020; Getman et al. 2019; Lim et al. 2020; Dalessandro et al. 2021b; Zeidler et al. 2021). While some of these systems are likely to dissolve, some may evolve into massive and long-lived clusters. Interestingly in this context, Dalessandro et al. (2021b) have found that the well-known clusters h and χ Persei are just components of an association of clusters embedded in a wide stellar halo of similar age. This structure, named LISCA I (where LISCA stands for Lively Infancy of Star Clusters and Associations), has provided the first detailed observational picture of an ongoing massive cluster hierarchical assembly. This is the first time that such a formation mechanism has been identified in the MW and it has important implications on our understanding of the environmental conditions (both locally and in the distant Universe) necessary to form massive stellar clusters. For many years, hierarchical cluster formation was invoked as the preferred dynamical route to form rotating and high-ellipticity star clusters (e.g., de Oliveira et al. 1998). More recently, it has been used to interpret the properties of particularly massive and dynamically complex clusters (e.g., Lee et al. 1999; Brüns & Kroupa 2011), and as an avenue to form clusters with multiple populations with different light elements abundances (e.g., Gavagnin et al. 2016; Hong et al. 2017). However, while hierarchical cluster formation is believed to work efficiently in high-density starburst galaxies (e.g., Bastian et al. 2011; Chandar et al. 2011), we are still missing an adequate understanding of its effectiveness in lower density environments, like the MW and the Magellanic Clouds.

Dynamical simulations (Bonnell et al. 2003; Ballone et al. 2020; Livernois et al. 2021) show that within a hierarchical assembly framework, the fragmentation of a molecular cloud may lead first to the formation of tens of small clumps (∼100 M_⊙). Then the surviving clumps merge to form a few more massive (10³ − 10⁴ M_⊙) and larger clusters. Finally, one or two clusters survive this hierarchical merger process and will eventually evolve into a single massive cluster (e.g., Fujii & Portegies Zwart 2016; Livernois et al. 2021).

As a part of a larger project aimed at constraining the occurrence of the hierarchical assembly process within local disk-like galaxies, and at testing whether it is able to form long-lived systems surviving the initial and turbulent few tens of million years of existence (e.g., Moeckel & Bate 2010; Gieles et al. 2006), we present a detailed photometric and kinematic study, mainly based on Gaia DR3 data (Gaia Collaboration 2023), of a region in the Galactic Perseus Arm including the clusters NGC 663 and NGC 654, that appears to be analogous to LISCA I. We also present a comparison of the observational results with the dynamical properties emerging in N-body models following the violent relaxation phase of a stellar system and its subsequent evolution.

The paper is structured as follows. The adopted data set is presented in Sect. 2. Sections 3 and 4 describe the physical properties of the area under study, its structure and kinematics, respectively. Section 5 presents the physical properties of star clusters belonging to the system. In Sect. 6 we discuss the total system’s mass. A comparison with a set of N-body simulations is described in Sect. 7. The main conclusions are drawn in Sect. 8.

2. Catalogs and preliminary analyses

2.1. The catalogs

From the Gaia Archive¹ we retrieved DR3 data for sources distributed within a large area on the sky (5° in radius) arbitrarily centered on the position of NGC 654 and having a five-parameter astrometric solution (i.e., sources with sky position, proper motion, and parallax measurements) and G < 19.5 mag. This catalog comprised 4.5 million sources.

We supplemented this data set with high-resolution optical and near-infrared spectra obtained with the HARPS-N (Cosentino et al. 2014) and GIANO-B (Oliva et al. 2012; Tozzi et al. 2016) spectrographs at the TNG as part of the SPA - Stellar Population Astrophysics: the detailed, age-resolved chemistry of the Milky Way disk Large Program (Program ID A37TAC13, PI: L. Origlia). Line-of-sight (LOS) velocities were obtained for all the observed stars, while detailed chemical abundances for the subsample of red supergiants were computed by Fanelli et al. (2022).

2.2. Clustering analysis

As any coherent stellar structure in the considered area should appear as an overdensity in the multi-dimensional phase-space of positions and velocities, we performed a clustering analysis on the whole catalog by means of the Hierarchical Density-Based Spatial Clustering of Application with Noise (HDBSCAN) algorithm (McInnes et al. 2017). For each star we used as inputs the Galactic coordinates, parallax, and proper motion components (ℓ,b, ϖ, μ_α*, μ_δ), and we set the HDBSCAN parameters as min_cluster_size = 40 and min_samples = 30. The parameter min_cluster_size sets a lower limit to the number of objects an overdensity should have to be identified as a cluster (hence we could not identify clusters with fewer than 40 members), while min_samples represents the number of sources used in determining the nearest neighbor distance for each source. Hence, increasing min_samples increases the mutual reachability distance among sources, and only the densest areas survive as clusters². Furthermore, HDBSCAN assigns a cluster membership probability to each star based on its distance from the neighboring stars. The closer the star is to the other cluster members, the higher the membership probability, and vice versa.

We identified 131 clustered systems within the full 5° field of view. To exclude spurious detections and select only systems that can be classified as clusters to a high significance level, we followed the post-processing approach described by Hunt & Reffert (2021), which uses the nearest-neighbors distance as a proxy for the local density. Only structures with a median value of the nearest-neighbors distance smaller than that of field stars at a 3σ level according to a Mann-Whitney statistics (Mann & Whitney 1947) were flagged as true stellar clusters. Out of 131 putative clusters, 54 systems fulfilled these criteria and were retained for the subsequent analysis. Recent open clusters catalogs (e.g., Cantat-Gaudin & Anders 2020; Cantat-Gaudin et al. 2020; Castro-Ginard et al. 2022) list 45 clusters in the region with more than 40 members (which is our threshold for identification). Interestingly, we recovered all the known clusters but two (hence 11 unknown structures were identified by this study), namely UBC 186 and UPK 265. We verified that UPK 265 could have been recovered by slightly changing the input parameters we set for the clustering analysis; however, it would have been excluded by the preliminary parallax selection (according to its value reported by Cantat-Gaudin et al. 2020) described below. The case of UBC 186 is more interesting. A careful investigation of its members reveals significant overlap with NGC 581 (128 out of 131 of NGC 581 members are in common with UBC 186; Cantat-Gaudin et al. 2020). Our analysis was able to properly identify both NGC 581 and another nearby structure that was labeled as UBC 186 by Cantat-Gaudin et al. (2020). However, the latter was flagged as a false detection by the adopted post-processing routine. It is important to note here, however, that the following analysis and the results of this paper do not depend on the inclusion or exclusion of any specific substructure or cluster.

Starting from the sample of 54 structures, we performed a preliminary selection to identify clusters sharing 3D position and 2D velocity with NGC 654 (ϖ = 0.31 ± 0.05 mas, μ_α* = −1.1 ± 0.1 mas yr⁻¹ and μ_δ = −0.3 ± 0.1 mas yr⁻¹, obtained by Cantat-Gaudin et al. 2020, using Gaia DR2 data), retaining only those with distance D = 2.8 − 3.2 kpc and co-moving within about 5.5 km s⁻¹ (corresponding to 0.38 mas yr⁻¹ at 3 kpc), according to their median parallax and proper motion estimated from Gaia DR3 data. Nine clusters (comprising NGC 654 itself) were selected in this way. We note in passing that none of the 11 previously unknown structures fulfilled these criteria.

Finally, we determined physically motivated selections in parallax and proper motion with the aim of selecting all the sources in the field of view sharing 3D position and 2D velocity with the nine clusters. Specifically, we inferred the intrinsic cluster distributions in parallax and proper motion (by using only stars with membership probability higher than 90%) by means of a Gaussian mixture modeling technique (the Extreme Deconvolution³ package developed by Bovy et al. 2011), thereby properly accounting for errors and correlation between measurements. In Fig. 1 we show the distributions inferred in this way in the parallax (top panel) and proper motion components (bottom panel).

Fig. 1. Inferred distributions in parallax (top panel) and proper motion (bottom panel) from likely (> 90%) cluster members. Different stellar clusters are in different colors. Dashed gray lines in the bottom panel are iso-probability contours at the 1, 2, and 3σ levels. The red dashed circle and vertical lines represent the range in proper motion and parallax inside which stars were selected.

We thus retained all the sources in the Gaia catalog with proper motions and parallaxes compatible, within 3σ, to the cluster distribution. Selected stars share similar distances, ϖ  ∈  [0.285;  0.407] mas, which corresponds to D  ∈  [2.46;  3.51] kpc, and co-move within 0.536 mas yr⁻¹ (about 7.5 km s⁻¹). In Fig. 2 we show the 2D density map of the region along with iso-density contours. The iso-density curves highlight the presence of small-scale clumpy structures corresponding to the identified stellar clusters (labeled in blue) as well as a lower-density diffuse halo extending for at least 3° from NGC 663 and NGC 654 and co-moving with the clusters.

Fig. 2. Spatial distribution in Galactic coordinates of stars selected in proper motion and parallax. Star cluster names are shown in blue, while black lines are iso-density contours enclosing 11.8% (solid), 39.3% (dashed), and 67.5% (dash-dotted) of the normalized star density distribution. The black arrow shows the direction of the Galactic center (GC).

2.3. Completeness of the Gaia catalog

The estimate of the Gaia catalog completeness is certainly a challenge due to, for instance, a composite data reduction pipeline and a complex satellite scanning law. Moreover, it does not depend only on the telescope properties themselves, but also on the physical properties of the observed regions such as crowding and extinction. This issue was first tackled by Everall et al. (2021, for the DR2) and Everall & Boubert (2022, for the EDR3) who directly modeled Gaia’s reduction pipeline and scanning law. More recently, Cantat-Gaudin et al. (2023) adopted an empirical approach to estimate the photometric completeness in the G band by comparing the Gaia catalog with the Dark Energy Camera Plane Survey (Schlafly et al. 2018; Saydjari et al. 2023).

For the purpose of this study, we need to assess the probability that a source is included in the catalog with magnitude G measure and a five-parameter solution. The resulting joint probability is

$$\begin{array}{c}\hfill p(5\mathrm{params}·G)=p(G)\times p(5\mathrm{params}|G).\end{array}$$(1)

We retrieved the first term on the right-hand side from the completeness maps of Cantat-Gaudin et al. (2023) (see Fig. 3 for the 2D map of photometric completeness computed for G = 19.5 mag), whereas the p(5 params | G) term in Eq. (1) was estimated from the number count ratios between sources with the five-parameter solution (k) compared to the total number of sources (n) for a given sky patch and magnitude bin

$$\begin{array}{c}\hfill p(5\mathrm{params}|G)=\frac{k+1}{n+2},\end{array}$$(2)

Fig. 3. Two-dimensional map of the p(G) term in Eq. (1) for G = 19.5 mag, computed with the GaiaUnlimited package (Cantat-Gaudin et al. 2023). White points show the location of star clusters.

following the documentation of the GaiaUnlimited project⁴. We computed Eq. (2) for a regular spatial grid in a region of 10° around NGC 654 assuming a spatial bin size of Δδ = 0.2° (thus it follows Δα = Δδ cos δ_NGC654 ≃ 0.42° in order to obtain a square grid) and for magnitude bins of ΔG = 0.2 mag down to G ≤ 19.5 mag.

In Fig. 4 we show as examples the 2D p(5params | G) maps computed for G = 18.5 mag (left panel) and G = 19.5 mag (right panel).

Fig. 4. Catalog completeness in the five-parameter solution estimated from star count ratios (p(5params|G) term in Eq. (2)) for G = 18.5 mag (left panel) and G = 19.5 mag (right panel). White points give the positions of star clusters.

The comparison between Figs. 3 and 4 clearly shows that the selection of sources with the five-parameter solution has a dominant impact on the final completeness. Indeed, p(5params | G) significantly drops for G > 18.5 mag, while p(G) remains almost equal to 1 down to G = 19.5 mag. Hence, in the following we assumed p(G) = 1, and thus simplifying Eq. (1) and Eq. (2) into:

$$\begin{array}{c}\hfill p(5\mathrm{params}·G)\simeq \frac{k+1}{n+2}\forall G<19.5\mathrm{mag}.\end{array}$$(3)

In the subsequent analyses we correct stellar counts for incompleteness according to Eq. (3).

3. Physical properties of the observed area

In Sect. 2.2 we identified a region encompassing nine star clusters embedded in a low-density and diffuse stellar halo (see Fig. 2) lying within strict ranges in 2D velocity, position, and parallax by construction. In this section we characterize the physical properties of the area based on the Gaia photometry and the spectroscopic data.

3.1. Differential reddening

Available Galactic extinction maps (e.g., Schlegel et al. 1998 and recalculations from Schlafly & Finkbeiner 2011) report a quite significant and strongly variable (E(B − V)∼0.5 − 3.5 mag) extinction along the LOS for the region under investigation. Here we provide an independent estimate of the differential reddening based on a suitable color-color diagram and following the approach adopted by Dalessandro et al. (2018). In particular, combining the Gaia G band with the r, i, z photometric bands from the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS data release 2, Chambers et al. 2016), we constructed the (G − r) versus (i − z) color-color diagram. This diagram turned out to be the most suitable choice as the evolutionary sequences run almost orthogonally to the reddening vector in these colors.

We derived differential extinction star by star by minimizing differences along the reddening vector with respect to a reference system. As a reference, we chose the median color-color distribution of likely main sequence stars (G > 12 mag or G_BP − G_RP < 0.5 mag) belonging to the cluster NGC 581. NGC 581 stars are distributed on average at bluer colors than other stars in the field, thus suggesting they are located in a region with relatively small extinction (color excess for these stars was derived from Schlegel et al. 1998; Schlafly & Finkbeiner 2011: E(B − V)_NGC581 ≃ 0.54 mag). Afterward, for each star, we obtained the median colors of the closest 50 likely main sequence (G > 12 mag or G_BP − G_RP < 0.5 mag) neighbor stars, and we determined the distance of this median value to the reference point along the reddening vector (using coefficients from Cardelli et al. 1989). The extinction value corresponding to the derived distance is then assigned to the specific star. To all the sources that do not fulfill the criteria of being likely main sequence stars and to those that do not have a counterpart in the PanSTARRS catalog, we assign the median reddening of the closest 50 neighbors.

As a representative example, we show in Fig. 5 the observed CMD of NGC 663 members as well as those obtained using extinction values from Schlafly & Finkbeiner (2011) (in red) and obtained in this work (in blue). The differential reddening corrections derived in this work nicely squeeze the sequence in the CMD compared to the observed one, and thus confirm the robustness of our estimates. On the contrary, those from Schlafly & Finkbeiner (2011) significantly spread the sequence and move stars at nonphysical colors, reaching (G_BP − G_RP)₀ = −2 mag, thus suggesting that the adopted values for the E(B − V) variations are likely overestimated.

Fig. 5. Color-magnitude diagram for NGC 663 members. In black the observed photometry, while in red and blue are shown the distributions of stars after correcting for differential reddening. For the former Galactic extinction maps by Schlegel et al. (1998), Schlafly & Finkbeiner (2011) are used, and for the latter the reddening corrections estimated in this work.

Finally, in Fig. 6 we show the resulting reddening map, as derived in this work, in which each star in the catalog is color-coded according to the inferred extinction. We note here that while differential reddening might play a role in shaping specific features of the iso-density contours shown in Fig. 2 (i.e., missing sources due to locally higher extinction artificially produces underdense regions), it is unlikely that it impacts the overall observed density gradient across the field of view as low-extinction regions, such as b ≲ −3° (see Fig. 6), still result as underdense compared to the diffuse halo. In Fig. 7 we show the observed (left panel) and differential reddening corrected (right panel) CMDs of the full catalog for comparison.

Fig. 6. Two-dimensional reddening map in Galactic coordinates. The extinction was computed star by star from color-color diagrams (see text for further details). White areas correspond to regions devoid of stars.

Fig. 7. Observed (left panel) and differential reddening corrected (right panel) CMD for the full catalog. In blue we show stars with LOS velocity measurements from Gaia DR3 (big squares for stars considered in the LOS analysis), whereas in red are stars targeted by high-resolution spectroscopy. The subsample with chemical abundances is circled in red.

3.2. Cluster and halo ages

Determining the ages of young (< 100 Myr) sparsely populated star clusters is certainly a challenge. At these ages, the color-magnitude distribution of turn-off stars is strongly affected by stellar rotation (Li et al. 2019). Moreover, the low number of stars and short evolutionary timescales of massive stars (> 10 M_⊙) might hamper a detailed estimate of the bright and blue main sequence termination, which in turn would bias the age inferred by standard methods such as isochrone fitting. Notwithstanding these limitations, here we attempt to tackle this issue by adopting a specific approach that is only marginally sensitive to stellar rotation and minimizes the impact of low-number statistics. In particular, we used a set of synthetic simple stellar populations obtained from the PARSEC database (Bressan et al. 2012) with [Fe/H] ≃ (−0.30 ± 0.01) dex (corresponding to the mean metallicity of the area Fanelli et al. 2022) and sampling the age range 1–100 Myr with a regular step of 1 Myr. We compared them with the observed cumulative luminosity functions (CLFs) in the G band, after correcting for differential reddening⁵ and completeness, as described in Sect. 2.3.

In order to account for number fluctuations, we randomly picked several times (≃100) a (virtually) independent sample of N stars from each synthetic population. The number of extracted stars N was set to be the number of objects in the synthetic population with G₀ > 6 mag after applying a normalization to the luminosity function in the range 12.5 < G₀ < 16 mag (at least 2 mag fainter than the main sequence termination for populations younger than 100 Myr at a distance of about 3 kpc). For each extraction, we then constructed the CLF, and we determined the median CLF of all the extractions (as well as its corresponding 68% credible region). The median CLF obtained for different ages was then compared to the observed one (see, e.g., Fig. 8) and the best fit was defined as the one that minimizes the χ² statistics. Comparison with the CLF was carried out up to G = 18 mag (corresponding to about G₀ ≃ 16 mag), below which the catalog’s completeness drops (see Fig. 4).

Fig. 8. Cumulative luminosity function for the stellar cluster NGC 654 (gray histogram), along with the normalized histogram of a 32 Myr old synthetic population computed with the PARSEC model (blue histogram). The error bars show the standard deviation of the model’s count fluctuations due to several extractions (see text for further details). Bins are 0.5 mag wide.

Furthermore, when comparing the synthetic CLF to the observed one, we looked for supergiant stars in the range G₀ < 8 mag and (G_BP − G_RP)₀ > 0.2 mag. If present, these stars provide strong constraints on the age, and hence we limited the analysis only to those ages that are able to explain the presence of evolved stars at the observed magnitudes. This allowed us to inform the fitting procedure about the likely young age of the system, even in absence of bright, blue main sequence stars. We note that with this procedure we assigned a narrower uniform prior to the cluster’s age. If red supergiants were not present we did not apply any selection on the age.

Color-magnitude diagrams for each cluster are shown in Fig. 9 along with best-fit isochrones. The nine clusters have ages in a narrow range of 14 − 44 Myr. The only exception is Berkeley 6, for which we derived an age of ${95}_{-15}^{+4}$ Myr. We also found a slight mismatch between isochrones and data, especially visible in the clusters Riddle 4 and NGC 654. This discrepancy likely results from local underestimations of the differential reddening, which in turn biases the age inference toward older ages.

Fig. 9. Color-magnitude diagrams of cluster members (black dots) corrected for differential reddening. For each cluster, the best-fit isochrone is shown in blue and the median age along with the 68% credible interval is also reported. The shaded areas show the region inside which stars are flagged as supergiants.

Nevertheless, typical errors in age estimates are about 10 − 15 Myr (Fig. 9). They account only for uncertainties arising from the fitting procedure, although errors in the differential reddening, distance, and incorrect membership assignment might also be important sources of uncertainties. However, we note that we are mainly interested in constraining relative ages rather than absolute ones.

A comparison with ages from the literature shows qualitatively overall agreement as they are in the range 15 − 38 Myr for the clusters under study (Cantat-Gaudin et al. 2020). The only exception is Berkeley 6 for which Cantat-Gaudin et al. (2020) report an age of about 200 Myr, consistent with the system being older than other clusters.

Finally, the same analysis was carried out for the stellar halo (i.e., all the stars that did not belong to any cluster according to the membership probabilities assigned by the clustering algorithm), finding that its age ($\sim {16}_{-1}^{+1}$ Myr) is consistent with the ages of the clusters embedded within it.

3.3. Line-of-sight velocity distribution and iron content

We investigated the LOS velocity and the metallicity distributions in the region using the TNG-GIARPS spectroscopic catalog presented in Sect. 2.1 (shown in red in Fig. 7), and we compared them with those expected for the surrounding Galactic field obtained from the Besançon Milky Way model (Robin et al. 2003) after applying the same parallax and proper motion selections. We computed LOS velocities for 24 stars, five of which are cluster members, while the remaining 19 belong to the halo. Among them, chemical abundances are available for the seven red supergiants (double red circles in Fig. 7), one of which belongs to NGC 581.

We also note that, while Gaia DR3 (Gaia Collaboration 2023) provides LOS velocity for 1164 selected stars, their color-magnitude distribution (shown in Fig. 7) suggests they are mostly field interlopers. Nevertheless, some bright stars that are likely members of the system have LOS measurements from Gaia. We therefore selected those stars with G₀ < 8 mag (removing objects belonging to the older disk population) and with rv_expected_sig_to_noise> 5 and rv_renormalised_gof< 2 (thus selecting sources with reliable LOS velocity; see Katz et al. 2019). Out of the 1164 stars, only 4 fulfilled these criteria and were thus included in the catalog (shown as larger blue squares in Fig. 7).

In Fig. 10 we show the distributions in LOS velocity (left panel, constructed with 24 stars from the high-resolution spectroscopic catalog plus 4 stars from Gaia DR3) and in [Fe/H] abundance (right panel, for the seven red supergiants), superimposed to the distributions of the surrounding Galactic field. The intrinsic widths of the observed distributions were inferred by means of a maximum likelihood approach (accounting for individual errors on measurements), and in Fig. 10 we report their median values along with the 68% credible intervals. In particular, we obtained a LOS velocity dispersion ${\sigma }_{\mathrm{LOS}}=8_{-1}^{+1}$ km s⁻¹ (to be compared with σ_LOS, MW = 20 km s⁻¹) and a mean velocity $⟨v_{\mathrm{LOS}}⟩=-{41}_{-2}^{+2}$ km s⁻¹, whereas for the metallicity we obtained a dispersion ${\sigma }_{[\mathrm{Fe}/\mathrm{H}]}=0.{008}_{-0.005}^{+0.009}$ dex (opposed to σ_[Fe/H] = 0.2 dex for the Galactic field) and a mean metallicity $⟨[\mathrm{Fe}/\mathrm{H}]⟩=-0.{30}_{-0.01}^{+0.01}$ dex.

Fig. 10. LOS velocity (left panel) and iron-over-hydrogen abundance (right panel) distributions for members of the selected structure (gray histogram) and for a Milky Way model (blue histogram). The intrinsic dispersions are also shown in the top left corners. All the distributions are scaled by a constant factor for visualization purposes only.

Interestingly, the observed distributions of the stars in the region are significantly narrower than those expected for a randomly selected group of co-moving Galactic stars, thus strengthening the evidence of kinematic coherence and suggesting a significant chemical homogeneity. In addition, the consistency in both LOS velocity and chemical content between cluster and halo stars further validates the assumption of a physical and coherent structure embedding the star clusters.

The literature data on the bulk clusters’ LOS velocity support the kinematic coherence of all clusters but one. The LOS velocity reported for Berkeley 6 (v_LOS ≃ −89 ± 52 km s⁻¹, Tarricq et al. 2021) is significantly lower than the system’s bulk velocity. However, we note that this value is still compatible within the huge uncertainty as only two stars were used for its estimate. Moreover, Spina et al. (2021) measured the iron content for one member of Berkeley 6 to be around [Fe/H] ∼ − 0.179 dex, significantly higher than [Fe/H] ≃ − 0.3 dex, although we note again that chemical abundances were derived for only one cluster member whose membership probability is < 30% (Spina et al. 2021). Therefore, better constraints on the stellar membership, age, and 3D velocity are needed before drawing any conclusions about the role of Berkeley 6 in the system.

4. Structural and kinematic properties of the diffuse stellar halo

4.1. Density distribution

We constructed the number density profile of the diffuse stellar halo. We took as the system’s center the center of mass of stars with M ≥ 2 M_⊙. First, celestial coordinates were converted into local Cartesian ones, assuming the centroid of the system as an initial guess for the system’s center. After that, the center of mass was computed in Cartesian coordinates and it was converted back into celestial coordinates, obtaining (α_CM; δ_CM) = (26.4559;  61.7865) degree.

We binned stars radially with respect to this center and we set the width of each radial annulus to contain 2500 sources each. Radial shells were then split into four angular sectors where the density was computed simply as the ratio of the number of stars to the sector’s area. The final shell density and error were the mean and standard deviation of the four measurements, respectively. Finally, we also accounted for Poissonian error in each bin by summing in quadrature to the standard deviation a term $1/\sqrt{N_{\mathrm{shell}}}$, with N_shell being the number of stars within the shell.

In Fig. 11 we show the number density profile for sources out to 8° from the system’s center of mass and with G ≤ 18 mag. When studying the density distribution, we temporarily extended the catalog up to 8° from the system’s center in order to assess the background density, while the selection in G was a good compromise between the catalog’s completeness and statistics. Interestingly, the observed density resembled a cluster-like profile over about a factor of 10 in density. At about R ≳ 6°, the density profile flattens, and we estimated the background density as the weighted mean of bins at distances larger than 6° from the adopted center, obtaining Σ_{⋆,background} ≃ 1.5 × 10⁻⁵ stars arcsec⁻², which was then subtracted from the observed profile.

Fig. 11. Stellar number density profile of an 8° region around the system’s center of mass and considering stars brighter than G = 18 mag. The observed profile is shown in gray, while the intrinsic profile (after background subtraction) is shown in black. The Plummer model (in red) and the King model (in blue) are also shown, along with the corresponding 68% credible regions constructed from the posterior samples.

Finally, we fitted the density distribution within 5° using the King (King 1962) and Plummer (Plummer 1911) models, which are typically adopted to reproduce stellar cluster density profiles. All the free parameters were constrained assuming a χ² likelihood and uniform priors (in logarithm), and exploring the parameter space with a Markov chain Monte Carlo (MCMC) technique using the Python package emcee⁶ (Foreman-Mackey et al. 2013).

In Fig. 11 we therefore show the density profile (before the background subtraction in gray, and after in black) along with the best-fit models and the associated errors. Both models provided a nice description of the data.

4.2. Kinematic properties

We investigated the kinematic properties of the stellar halo by further selecting stars fulfilling the following astrometric quality selection criteria (Lindegren et al. 2021): ruwe ≤1.4, astrometric_gof_al ≤ 1 and astrometric_excess_noise ≤ 1 mas (if astrometric_excess_noise_sig >  2), thus excluding those sources for which the standard five-parameter solution does not provide a reliable fit of the observed data.

First, we accounted for perspective effects induced by the system’s bulk motion (van Leeuwen 2009) on the μ_α* and μ_δ components; we thus corrected the velocities for each star assuming a bulk average motion of $(\overline{{\mu }_{\alpha \ast }};\overline{{\mu }_{\delta }})=(-1.14;-0.33)$ mas yr⁻¹ (estimated using Gaia data for sources with reliable astrometry) and the mean LOS velocity obtained from the spectroscopic catalog supplemented with Gaia DR3 data (see Sect. 3.3). Owing to the large area of the sky covered by our data, the magnitude of the perspective correction resulted as non-negligible, reaching up to 0.2 mas yr⁻¹ (about 2.8 km s⁻¹) at 5° from the center, and hence caution must be taken when interpreting results at such large angular scales.

Looking at the distribution of velocities on the plane of the sky may offer a first glimpse of the dynamic state of the system. We thus performed a centroidal Voronoi tessellation (Cappellari & Copin 2003), exploiting the density profile shown in Fig. 11, such that each bin contains about the same number of stars.

Radial velocities⁷ in each bin were inferred by means of an MCMC exploration assuming as likelihood (see, e.g., Pryor & Meylan 1993; Raso et al. 2020)

$$\begin{array}{cc}& ln\mathcal{L}=-\frac{1}{2}{\displaystyle \sum _i}[\frac{{(v_{\mathrm{R},\mathrm{i}}-\overline{v_{\mathrm{R}}})}^2}{{\sigma }_{\mathrm{R}}^2+e_{\mathrm{R},\mathrm{i}}^2}+ln({\sigma }_{\mathrm{R}}^2+e_{\mathrm{R},\mathrm{i}}^2)\hfill \\ \hfill & +\frac{{(v_{\mathrm{T},\mathrm{i}}-\overline{v_{\mathrm{T}}})}^2}{{\sigma }_{\mathrm{T}}^2+e_{\mathrm{T},\mathrm{i}}^2}+ln({\sigma }_{\mathrm{T}}^2+e_{\mathrm{T},\mathrm{i}}^2)],\hfill \end{array}$$(4)

where v_X, i and e_X, i with X ∈ {R, T} are the radial (R) and tangential (T) components of the velocity and error for the ith star, respectively. Furthermore, we assumed uniform priors in the logarithms of the velocity dispersion (σ_R and σ_T) and uniform priors in the mean velocities ($\overline{v_{\mathrm{R}}}$ and $\overline{v_{\mathrm{T}}}$). In Fig. 12 we show the mean velocity vectors for each tile, and we color-coded stars in each Voronoi bin according to the mean radial velocity inferred in the bin. The directions of the arrows clearly show a contraction of the external regions (reaching speeds up to ∼4 − 5 km s⁻¹), also confirmed by the color distribution of stars in Fig. 12. Interestingly, in the central regions (R < 1 − 2 deg), a mild expansion on the order of ≃1 km s⁻¹ is observed, mainly visible in the purplish bins.

Fig. 12. Two-dimensional distribution of stars in local Cartesian coordinates (x; y). Stars are color-coded according to the inferred mean radial component of the velocity in their Voronoi bin, while black arrows show the mean velocity vector in each bin. Radial velocity is defined as positive if pointing away from the center, thus $\overline{v_{\mathrm{R}}}>0$ means expansion and $\overline{v_{\mathrm{R}}}<0$ means contraction. The velocity scale is shown in the bottom left corner.

The same pattern emerges when computing the mean radial velocity $\overline{v_{\mathrm{R}}}$ in spherical shells (see Eq. (4)), as shown in Fig. 13. In particular, the innermost 2° show a flat slightly positive profile ($\overline{v_{\mathrm{R}}}/{\sigma }_{\mathrm{R}}>0$) indicating central expansion, although the radial motion is highly dominated by random motion ($\overline{v_{\mathrm{R}}}/{\sigma }_{\mathrm{R}}<0.1$). However, moving toward larger radii, the contraction ($\overline{v_{\mathrm{R}}}/{\sigma }_{\mathrm{R}}<0$) becomes increasingly more prominent and the radial motion more ordered.

Fig. 13. Radial profile of the ratio of the radial mean velocity to the radial velocity dispersion computed in spherical shells. Black dots show the median values, while quoted errors are the 16th and 84th percentiles of the distributions. The dashed horizontal line shows the zero expansion–contraction level.

The robustness of the contraction pattern (Figs. 12 and 13) was tested against the assumption of a particular bulk LOS velocity when accounting for perspective expansion (since we had only a few measurements). We considered the worst-case scenario in which stars followed the LOS velocity distribution expected for the MW field stars in the region ($\overline{v_{\mathrm{LOS}}}\simeq -35$ km s⁻¹ and σ_LOS, MW ≃ 20 km s⁻¹, see Fig. 10).

The transverse motions were then corrected by randomly assigning to each star a LOS velocity extracted from the MW-like distribution. We iterated this procedure 100 times finding that the velocity pattern observed was weakly affected by our mean bulk motion assumption and the contraction showed in Fig. 13 was always recovered.

4.3. Mass dependence analysis

The presence of a mass spectrum has a non-negligible role in the dynamics of both young and old stellar systems. Old stellar clusters are indeed known to naturally develop mass segregation due to two-body interactions that cause significant kinetic energy exchange among stars and cause massive stars to sink toward the cluster’s center (Binney & Tremaine 2008). However, evidence of mass segregation has also been found in younger Galactic clusters (e.g., Hillenbrand & Hartmann 1998; Gouliermis et al. 2004; Stolte et al. 2006; Evans & Oh 2022), thus possibly implying a connection with the early stages of cluster formation (McMillan et al. 2007; Allison et al. 2009; Livernois et al. 2021). In addition, numerical simulations showed that during the violent relaxation phase, young stellar systems can start developing a dependence of the kinematical properties (rotation, velocity dispersion) on the stellar mass (see, e.g., Livernois et al. 2021). The investigation of possible mass-dependent dynamical properties is therefore crucial to shed further light on the dynamics of young stellar systems.

We estimated stellar masses using theoretical M–G₀ relation for zero age main sequence stars. This relation was obtained from the PARSEC models for a population of 14 Myr in age (the youngest age we estimate) with [Fe/H]= − 0.3 dex. Stellar masses were therefore derived by interpolation of this relation, and in Fig. 14 we show the spatial distribution of stars color-coded by their mass.

Fig. 14. Spatial distribution of stars color-coded by mass. Masses were obtained by mass-absolute magnitude relation (see text for further details).

Figure 14 indicates that massive stars (lighter symbols) are more centrally concentrated than lower mass ones (darker symbols). We thus quantitatively investigated this feature by looking at the cumulative profiles in different mass bins, after correcting for completeness (see Fig. 4) by assigning to each star a weight = 1/completeness. We limited the analysis to stars more massive than 1.1 M_⊙ to avoid very low completeness values. Figure 15 shows that the massive stars exhibited a more centrally concentrated spatial distribution than the lower mass ones. In addition, a Kolmogorov-Smirnov test confirmed the statistical significance of this result (p < 0.003) for every combination of mass bins. In Appendix A we show the completeness distributions for each mass bin.

Fig. 15. Cumulative radial profiles constructed in four different mass bins from stars with masses between 1.1–2 M⊙ (in red) up to stars more massive than 10 M⊙ (in dark blue). Stellar counts were corrected for incompleteness star by star.

Finally, we looked for evidence of a dependence of the kinematic properties on the stellar masses. In Fig. 16 we show the velocity dispersion profiles constructed for different mass bins, specifically for stars with mass M ≤ 2 M_⊙ (red symbols), 2 M_⊙  < M ≤ 4 M_⊙ (light blue symbols) and M > 4 M_⊙ (dark blue symbols).

Fig. 16. Velocity dispersion profiles in the radial (left panel) and tangential (right panel) components, in three different mass bins depicted with different colors: red for the least massive bin and dark blue for the most massive.

Our analysis reveals mild evidence of a dependence of the velocity dispersion on the stellar mass within about 2° from the center in both the radial and tangential components: more massive stars show a slightly smaller velocity dispersion. At radii R > 2 − 3 deg the velocity dispersion profiles become indistinguishable and no signs of equipartition are found. The trend is consistent with that expected for a stellar system that has started to evolve toward energy equipartition during its early evolutionary phases and is in general agreement with the trends found in the simulations presented in Livernois et al. (2021). We discuss this point further in Sect. 7.

Finally, we observed an increase in the velocity dispersion moving away from the center in every mass bin (variations up to ∼0.6 km s⁻¹, see Fig. 16). Several effects might be at play in driving such a pattern, for instance deviation from spherical symmetry (see, e.g., the color-coded map in Fig. 12, whose expansion pattern was clearly nonspherical), a nonconstant mean velocity within bins (as noted by Da Rio et al. 2017), and possibly the presence of residual field interlopers.

We also observed a clear dependence of the radial velocity component on the stellar mass. In particular, stars more massive than about ≳2 M_⊙ exhibit a higher positive mean radial velocity (up to about 1 km s⁻¹) than lower mass stars within the innermost 3° (see Fig. 17), while at larger radii massive stars contract toward the center ($\overline{v_{\mathrm{R}}}<0$) with a similar slope to lower mass stars, albeit with a higher normalization.

Fig. 17. Mean radial velocity profiles for stars in three different mass bins: M ≤ 2 M⊙ (red symbols), 2 < M ≤ 4 M⊙ (light blue symbols), and M > 4 M⊙ (dark blue symbols). Positive values of $\overline{v_{\mathrm{R}}}$ should be interpreted as expansion and negative values as contraction.

We thus further investigated this feature in the 2D plane using the Voronoi tessellation, giving particular attention to non-spherically symmetric features and putative links with the spatial distribution of stars. In Fig. 18 we show 2D maps of mean radial velocity for stars with mass above (left panel) and below (right panel) 2 M_⊙. Expansion ($\overline{v_{\mathrm{R}}}>0$) is depicted in purple, contraction ($\overline{v_{\mathrm{R}}}<0$) is depicted in blue, and black lines highlight iso-probability contours of the density distribution of the respective populations.

Fig. 18. Two-dimensional maps of mean radial velocity for high-mass stars (M ≥ 2 M⊙, left panel) and low-mass stars (right panel). Stars belonging to the same Voronoi bin are color-coded according to the mean radial velocity, while black lines are iso-density contours of the respective populations enclosing about 12% (solid line), 39% (dashed line), and 67% (dash-dotted line) of the underlying density distribution.

Interestingly, a clear connection between expansion and density comes up when looking at massive stars (left panel of Fig. 18), suggesting that higher-density regions expand faster than lower ones, whereas no clear connection is found for the low-mass population (right panel). In addition, we found consistent features to what is observed in Fig. 17: massive stars expand faster than lower mass ones, which in turn had higher contraction speeds in the outskirts.

The spectro-photometric (i.e., age and metallicity), structural (i.e., density profiles and mass segregation), and kinematic (i.e., contraction and equipartition) evidence collected so far suggests that the stars selected within a few degrees of NGC 654 are not just a group of co-moving stars; rather, the nine identified clusters and the extended low-density halo surrounding them are part of a common, substructured, and still-forming massive stellar system. Following the definition introduced in Dalessandro et al. (2021b), we named it LISCA II.

5. Structure and kinematic of the embedded star clusters

In this section we present the structural and kinematic properties of the nine star clusters composing LISCA II.

First, we constructed the number density radial profiles for all the clusters (following the same approach described in Sect. 4) with respect to their center of mass, and we fitted them with the Plummer (Plummer 1911) and King (King 1962) models. Figure 19 shows the observed profiles and the models, normalized to the clusters’ central densities and observed projected half-mass radii (R_hm) enabling a quantitative comparison among different clusters. The quantity R_hm is defined as the projected radius that encloses half of the total cluster’s mass directly obtained from the radial distribution of member stars.

Fig. 19. Number density profiles (black symbols) for the nine star clusters. Also shown are the model predictions for the Plummer (red) and King (blue) models. Shaded areas represent the 68% credible regions. All the profiles, both observed and models, were scaled to the central density predicted by the Plummer model Σ0,  plummer and to the observed projected half-mass radius Rhm.

Every star cluster exhibits a cluster-like profile that, within the errors, is equally well modeled by both the Plummer and King models, with the only exception of NGC 581 (see Fig. 19), which exhibits a sharper truncation not captured by the Plummer model.

In addition, we present the kinematic properties of star clusters: Fig. 20 shows the inferred 1D velocity dispersion, which is defined as

$$\begin{array}{c}\hfill {\sigma }_{1\mathrm{D}}\equiv \sqrt{\frac{{\sigma }_{\mathrm{R}}^2+{\sigma }_{\mathrm{T}}^2}{2}},\end{array}$$(5)

Fig. 20. One-dimensional velocity dispersion profiles for seven out of nine star clusters. Velocity dispersion was normalized to the same quantity computed for all the cluster’s members ($\overline{{\sigma }_{1\mathrm{D}}}$, i.e., without binning the stars), while projected distances from the cluster center of mass were normalized to the observed half-mass radius.

with σ_R and σ_T being the projected radial and tangential components of the velocity dispersion inferred assuming the likelihood in Eq. (4) and sampling the parameter space with an MCMC technique.

For two clusters, namely Riddle 4 and Berkeley 6, we could not compute reliable kinematic profiles since few stars (51 and 84 respectively) fulfilled the astrometric quality selections presented in Sect. 4.

All the other clusters exhibit rather flat dispersion profiles that can be hardly explained by equilibrium models, for instance those adopted to fit the stellar density distributions. Several effects might be at play in producing the observed flat dispersion profiles, such as contaminants from the stellar halo and dynamical heating due to tidal interactions with other clusters and substructures possibly taking place during the system’s early evolution. We discuss dynamical heating in more detail in Sect. 7.

To estimate the total mass of each cluster we followed the same approach described in Dalessandro et al. (2021b), and we normalized a Kroupa (2001) and a Salpeter (1955) initial mass function (IMF) in the range M > 4 M_⊙, such that the number of member stars matched the number predicted by direct integration of the IMF

$$\begin{array}{c}\hfill N_{\star \mathrm{observed}}={\displaystyle {\int }_{m_{\min }}^{m_{\max }}}\mathrm{IMF}(m)\mathrm{d}m,\end{array}$$(6)

with m_min = 4 M_⊙ and m_max = 11 − 14 M_⊙ for every cluster but Berkeley 6, whose maximum stellar mass is approximately 7 M_⊙ according to its inferred age. Once the normalization was obtained, we computed the total visible mass by integrating the IMF in the range M_min − M_max = 0.09 − 14.73 M_⊙

$$\begin{array}{c}\hfill M_{\mathrm{tot}}={\displaystyle {\int }_{M_{\min }}^{M_{\max }}\mathrm{IMF}(m)m\mathrm{d}m,}\end{array}$$(7)

respectively the minimum and maximum stellar mass for a 14 Myr-old simple stellar population with [Fe/H]= − 0.3 dex (Bressan et al. 2012). The derived cluster masses are in the range ≃0.5 − 5.6 × 10³ M_⊙ (≃0.7 − 8.8 × 10³ M_⊙) according to a Kroupa (Salpeter) IMF.

The main kinematic and structural properties of the clusters are summarized in Table 1. Specifically, for each cluster we included the coordinates of the center of mass (for stars more massive than 2 M_⊙), the mean parallax and proper motion (of the distributions shown in Fig. 1) along with errors, the Plummer scale length and the King core radius (obtained from the models shown in Fig. 19 and converted to parsec using the mean parallax also reported in Table 1), the total masses inferred by assuming either a Kroupa or a Salpeter IMF, and finally the total number of member stars for each cluster.

6. Total system mass

The total mass of LISCA II was estimated by using Eq. (7) and roughly the same approach as for the single clusters. However, in this case, we had to assume a radial extension of the system within which to integrate the stellar masses. To this end we used the Jacobi radius (R_J; see, e.g., Binney & Tremaine 2008) as a first-order physically plausible radial extension of the system. The radius R_J is simply defined as

$$\begin{array}{c}\hfill R_{\mathrm{J}}=R_0{\left(\frac{M_{\mathrm{LISCAII}}}{3M_{\mathrm{MW}}(<R_0)}\right)}^{1/3},\end{array}$$(8)

with R₀ ≃ 10.2 kpc and $M_{\mathrm{MW}}(<R_0)=10.{98}_{-0.10}^{+0.12}\times {10}^{10}$ M_⊙ being the galactocentric distance of the system and the MW mass enclosed within that radius, respectively (Cautun et al. 2020), and M_LISCAII is the total system mass. Starting from Eq. (8), we estimated both R_J and the mass of the system by using an iterative procedure until final convergence was reached. Depending on the assumed IMF we obtain

$$\begin{array}{cc}& R_{\mathrm{J},\mathrm{Kroupa}}=1.07{}^{°}\simeq 55\mathrm{pc},\hfill \\ \hfill & R_{\mathrm{J},\mathrm{Salpeter}}=1.35{}^{°}\simeq 70\mathrm{pc},\hfill \end{array}$$(9)

with the corresponding enclosed (< 2R_J) masses being

$$\begin{array}{cc}& M_{\mathrm{LISCAII},\mathrm{Kroupa}}=6.4\times {10}^4{M}_{\odot },\hfill \\ \hfill & M_{\mathrm{LISCAII},\mathrm{Salpeter}}=1.2\times {10}^5{M}_{\odot }.\hfill \end{array}$$(10)

It is important to note here that in this case R_J is only meant to provide a general indication of the spatial scale related to the strength of the tidal field at the location of LISCA II. In the complex case of a clumpy system far from a spherical configuration in dynamical equilibrium like that of LISCA II, the detailed implications of the effects of the tidal field truncation for this system would require a tailored set of simulations. We also note that the kinematical properties revealed by our analysis (Sect. 4) indicate that in this system all the stars, including those beyond the present estimate of R_J, are strongly contracting toward the center of the system, which is just the opposite of what is expected for stars escaping the system.

7. Comparison with N-body simulations of early cluster evolution

In this section we present an analysis of the dynamical properties of one of the N-body models studied in Livernois et al. (2021). The simulations of Livernois et al. (2021) explored the early evolution and violent relaxation phases of young rotating star clusters and followed their evolution from a hierarchical structure to a final monolithic equilibrium configuration.

We note that our goal here is not to build a detailed model of the LISCA II system, but rather to gather further insight into our observational analysis. Hence, here we just provide a theoretical example of the general dynamical properties expected in a hierarchical stellar cluster undergoing its early evolutionary phases. The model we analyzed for this paper is the F025 model of Livernois et al. (2021). The model is fully described in Livernois et al. (2021), but we summarize its main features for the purposes of this study. The model starts with 10⁵ stars with mass range 0.08 − 100 M_⊙, initially following a fractal distribution with a fractal dimension equal to 2.6; the system is initially dynamically cold and undergoes the collapse and subsequent structural oscillations typical of the violent relaxation phase. As a reference timescale for the presentation of our results, we adopt the system free-fall timescale, t_ff, which corresponds approximately to the timescale needed for the system to reach its maximum contraction during its initial collapse. Assuming an initial mass of 5  ×  10⁴ − 10⁵ M_⊙ and an initial radius of ∼50 − 60 pc, the t_ff would roughly correspond to 20 − 35 Myr. To capture the model at multiple evolutionary stages, we focus our attention on the snapshots at the following values of t/t_ff: 0.9 (denoted TFF09), 1.1 (TFF11), 1.3 (TFF13), 1.5 (TFF15), and 1.7 (TFF17).

7.1. Mass segregation and bulk internal motion

We start our analysis with the study of mass segregation in the N-body model and the LISCA II system. Here we focus on an analysis specifically aimed at detecting mass segregation on a local scale, which might provide further insight into the dynamics of systems characterized by the presence of clumps and substructures like those studied here. To quantify the level of mass segregation, we first calculate the local surface number density for each star using the distance of the sixth nearest star of any mass (Casertano & Hut 1985); the median local surface density of different mass bins, normalized by the median local surface density of the entire cluster, is plotted against the median stellar mass of each bin in Fig. 21. All snapshots show clear evidence of a local surface density increasing with the stellar mass, implying the presence of local-scale mass segregation, where massive stars are migrating toward the centers of the subclusters they are members of. A similar trend is also present in LISCA II, which shows significant local-scale mass segregation, complementing the global-scale mass segregation found in Fig. 15. This trend appears to evolve with time in the model snapshots; the two latest snapshots have the strongest trend of local-scale mass segregation.

Fig. 21. Local number density of stars as a function of stellar mass (see Sect. 7.1). The local density is normalized by the median local surface density of the entire cluster. The increase in local density with stellar mass indicates mass segregation on local scales across all clusters shown.

Focusing on the bulk internal motion of the system, we analyzed the radial velocity profile of the simulation data as a comparison to Fig. 13. In Fig. 22, we plot the radial velocity profile normalized by the velocity dispersion of the cluster, including only stars with m > 2 M_⊙. As the cluster evolves, we see the different regions of the cluster transition between expansion and contraction. The TFF09 snapshot is characterized by a trend similar to that found in LISCA II: an expansion of the inner regions and a strong contraction in the outer regions.

Fig. 22. Projected radial velocity profiles for each snapshot of the N-body simulation and the LISCA II system. The projected radius is normalized to the projected half-mass radius of the system in each snapshot, and the radial velocity is normalized to the radial velocity dispersion of all stars in the system.

7.2. Dynamics of subclusters

For insight into the possible dynamical evolution of the subclusters within LISCA II, we analyzed the dynamical properties of selected clumps from our N-body models. For each snapshot in the N-body data, we selected a few clumps in spherical 3D regions and determined their centers as the location of the maximum local density. These clumps are in dynamically active environments, and no clustering metrics were found to be appropriate across all snapshots for clump identification.

We start by showing in Fig. 23 the surface density profile (top panel) and LOS velocity dispersion profile (bottom panel) for the clumps selected in the TFF09 snapshot. We fit a Plummer model based on the surface density profile and enclosed mass within one projected half-mass radius, and overplot the best-fit model lines. All clumps show radial variation of the surface density profiles following the general shape of the Plummer model; on the other hand, the velocity dispersion profiles are flatter and more elevated than is expected from the best-fit Plummer models. This is the manifestation of the tidal heating in the cluster environment and is similar to what is seen in Fig. 20.

Fig. 23. Surface density profiles (top) and 1D velocity dispersion profiles (bottom) for the TFF09 snapshot. The best-fit Plummer model for each clump is also shown (dashed lines). The best-fit Plummer model is determined by fitting the surface density profile within one projected half-mass radius. The surface density and 1D velocity dispersion are both normalized to the corresponding values in the innermost radial bin of each clump.

We repeated the above analysis for all of our snapshots, and plot in Fig. 24 the average ratio of the surface density (top panel) and LOS velocity dispersion (bottom panel) to the best-fit Plummer models across the three biggest clumps in each snapshot. The surface density fits well out to approximately one projected half-mass radius, outside of which the clumps generally have a higher density than the best-fit Plummer model; these deviations from the best-fit Plummer model can be attributed to the high density of the surrounding environment and the perturbations due to interactions in the cluster environment. The effects of these perturbations are clearly visible in the LOS velocity dispersion profiles of all the snapshots analyzed. All the velocity dispersion profiles deviate from the profiles expected from the best-fit Plummer models across all radii with a dependence on the projected radius that varies in different snapshots.

Fig. 24. Radial profiles of the ratio of the clump radial profiles to the best-fit Plummer models for surface density (top) and 1D velocity dispersion (bottom). This ratio is averaged over the three largest clumps within each snapshot.

Finally, we summarized the relation between structural and kinematic perturbations by plotting, in Fig. 25, the ratio of density and LOS velocity dispersion to the corresponding values of these quantities from the best-fit Plummer models at R_h for the three biggest clumps in each snapshot. This plot clearly shows that the fingerprints of the highly active environment, where each clump is undergoing rapid interactions, mergers, and fragmentations, are more evident in the kinematic properties, as illustrated by the larger deviations of the velocity dispersion from the expected equilibrium values.

Fig. 25. One-dimensional velocity dispersion vs. projected surface density for the three largest clumps in each snapshot (see legend), each calculated at one projected half-mass radius. The 1D velocity dispersion and the projected surface density are both normalized by the value expected from the Plummer model best fitting the surface density profile. Point sizes are proportional to the total mass of each clump.

8. Discussion and conclusions

The unprecedented quality of the Gaia DR3 data (Gaia Collaboration 2023, supplemented by high-resolution spectra (Fanelli et al. 2022) obtained as part of the SPA-TNG large program, allowed us to identify the LISCA II system in the Perseus complex.

The spectro-photometric, structural, and kinematical properties of this system are in generally good agreement with those theoretically expected from the early dynamical evolution of a massive molecular cloud that experienced violent relaxation and is now in the process of hierarchically assembling its stellar constituents and evolving toward a monolithic structure. In particular, the observed evidence of mass segregation on a local and global scale, the mass-dependent kinematic properties, and out-of-equilibrium internal kinematics of individual subclusters, as well as the observation of a dominant contraction pattern toward the system center mainly driven by the external regions of LISCA II and of a milder central expansion, are compatible with what is expected in the early evolutionary phases of stellar systems assembling as a coherent massive structure by N-body models. The properties of these hierarchical stellar systems are shaped by a combination of large-scale variations of the system’s potential and smaller scale interactions of individual subclusters and clumps.

Although more detailed models and additional data would be necessary to further explore the possible fate of this system, the evidence collected in this paper suggests that LISCA II is a good candidate to evolve into a young massive (10⁴ − 10⁵ M_⊙) star cluster on a timescale of ∼100 Myr, corresponding to a few free-fall times. These results make LISCA II the second structure, the first being LISCA I (Dalessandro et al. 2021b), ever found in the MW in the process of hierarchically assembling in a massive stellar cluster. LISCA II is located at only ∼6 deg from LISCA I, with which it shares similar chemical composition ([Fe/H] =  − 0.30 dex), age (t ∼ 20 Myr), and overall mass.

In conclusion, the present analysis has provided a comprehensive characterization of the process of cluster assembly with a level of detail that cannot be achieved in external galaxies or at high redshift (where the progenitors of the oldest clusters formed), thus showing that probing cluster formation in local environments can help shed light on the physical processes involved in massive cluster formation and their role in determining the cluster’s dynamical properties. Moreover, we further showed that hierarchical cluster assembly is a viable process, even in low-density environments such as the MW (former observational evidence was mainly in high-density environments, e.g., Bastian et al. 2011; Chandar et al. 2011), and a statistical assessment of its effectiveness on Galactic scales is the subject of an ongoing study. It is interesting to note in this respect that the possible observed internal age spreads of the stellar populations belonging to LISCA I and II (∼10 Myr) nicely fit the observed trend (Parmentier et al. 2014) between the cluster formation environment stellar density and final cluster age internal variations, which possibly results from the different duration of the star formation processes and the link between their efficiency and the systems’ free-fall time. This further strengthens the idea that clusters with different present-day properties likely underwent similar formation processes.

Acknowledgments

The authors thank the anonymous referee for the careful reading of the paper and the useful comments and suggestions. A.D.C., E.D. and L.O. acknowledge financial support from the project Light-on-Dark granted by MIUR through PRIN2017-2017K7REXT contract. E.D. acknowledges support from the Indiana University Institute for Advanced Study through the Visiting Fellowship program. This work uses data from the European Space Agency (ESA) space mission Gaia. Gaia data are being processed by the Gaia Data Processing and Analysis Consortium (DPAC). Funding for the DPAC is provided by national institutions, in particular, the institutions participating in the Gaia Multi-Lateral Agreement (MLA).

References

Appendix A: Completeness distributions in different mass bins

In this Appendix we show the completeness distributions for the five-parameter solutions in the four different mass bins used to investigate the presence of mass segregation (see Figure 15). As can be seen, the higher the mass, the higher the completeness. The distribution width is due to both a finite mass range in the bin and differential reddening, which link the estimated mass to the observed G magnitude.

Fig. A.1. Completeness distribution in four mass bins. Colors are consistent with Figure 15 and the mass range is shown in the upper left corner.

All Tables

All Figures

	Fig. 1. Inferred distributions in parallax (top panel) and proper motion (bottom panel) from likely (> 90%) cluster members. Different stellar clusters are in different colors. Dashed gray lines in the bottom panel are iso-probability contours at the 1, 2, and 3σ levels. The red dashed circle and vertical lines represent the range in proper motion and parallax inside which stars were selected.
In the text

	Fig. 2. Spatial distribution in Galactic coordinates of stars selected in proper motion and parallax. Star cluster names are shown in blue, while black lines are iso-density contours enclosing 11.8% (solid), 39.3% (dashed), and 67.5% (dash-dotted) of the normalized star density distribution. The black arrow shows the direction of the Galactic center (GC).
In the text

	Fig. 3. Two-dimensional map of the p(G) term in Eq. (1) for G = 19.5 mag, computed with the GaiaUnlimited package (Cantat-Gaudin et al. 2023). White points show the location of star clusters.
In the text

	Fig. 4. Catalog completeness in the five-parameter solution estimated from star count ratios (p(5params|G) term in Eq. (2)) for G = 18.5 mag (left panel) and G = 19.5 mag (right panel). White points give the positions of star clusters.
In the text

	Fig. 5. Color-magnitude diagram for NGC 663 members. In black the observed photometry, while in red and blue are shown the distributions of stars after correcting for differential reddening. For the former Galactic extinction maps by Schlegel et al. (1998), Schlafly & Finkbeiner (2011) are used, and for the latter the reddening corrections estimated in this work.
In the text

	Fig. 6. Two-dimensional reddening map in Galactic coordinates. The extinction was computed star by star from color-color diagrams (see text for further details). White areas correspond to regions devoid of stars.
In the text

	Fig. 7. Observed (left panel) and differential reddening corrected (right panel) CMD for the full catalog. In blue we show stars with LOS velocity measurements from Gaia DR3 (big squares for stars considered in the LOS analysis), whereas in red are stars targeted by high-resolution spectroscopy. The subsample with chemical abundances is circled in red.
In the text

	Fig. 8. Cumulative luminosity function for the stellar cluster NGC 654 (gray histogram), along with the normalized histogram of a 32 Myr old synthetic population computed with the PARSEC model (blue histogram). The error bars show the standard deviation of the model’s count fluctuations due to several extractions (see text for further details). Bins are 0.5 mag wide.
In the text

	Fig. 9. Color-magnitude diagrams of cluster members (black dots) corrected for differential reddening. For each cluster, the best-fit isochrone is shown in blue and the median age along with the 68% credible interval is also reported. The shaded areas show the region inside which stars are flagged as supergiants.
In the text

	Fig. 10. LOS velocity (left panel) and iron-over-hydrogen abundance (right panel) distributions for members of the selected structure (gray histogram) and for a Milky Way model (blue histogram). The intrinsic dispersions are also shown in the top left corners. All the distributions are scaled by a constant factor for visualization purposes only.
In the text

Fig. 11. Stellar number density profile of an 8° region around the system’s center of mass and considering stars brighter than G = 18 mag. The observed profile is shown in gray, while the intrinsic profile (after background subtraction) is shown in black. The Plummer model (in red) and the King model (in blue) are also shown, along with the corresponding 68% credible regions constructed from the posterior samples.

In the text

Fig. 12. Two-dimensional distribution of stars in local Cartesian coordinates (x; y). Stars are color-coded according to the inferred mean radial component of the velocity in their Voronoi bin, while black arrows show the mean velocity vector in each bin. Radial velocity is defined as positive if pointing away from the center, thus $\overline{v_{\mathrm{R}}}>0$ means expansion and $\overline{v_{\mathrm{R}}}<0$ means contraction. The velocity scale is shown in the bottom left corner.

In the text

	Fig. 13. Radial profile of the ratio of the radial mean velocity to the radial velocity dispersion computed in spherical shells. Black dots show the median values, while quoted errors are the 16th and 84th percentiles of the distributions. The dashed horizontal line shows the zero expansion–contraction level.
In the text

	Fig. 14. Spatial distribution of stars color-coded by mass. Masses were obtained by mass-absolute magnitude relation (see text for further details).
In the text

	Fig. 15. Cumulative radial profiles constructed in four different mass bins from stars with masses between 1.1–2 M⊙ (in red) up to stars more massive than 10 M⊙ (in dark blue). Stellar counts were corrected for incompleteness star by star.
In the text

	Fig. 16. Velocity dispersion profiles in the radial (left panel) and tangential (right panel) components, in three different mass bins depicted with different colors: red for the least massive bin and dark blue for the most massive.
In the text

	Fig. 17. Mean radial velocity profiles for stars in three different mass bins: M ≤ 2 M⊙ (red symbols), 2 < M ≤ 4 M⊙ (light blue symbols), and M > 4 M⊙ (dark blue symbols). Positive values of $\overline{v_{\mathrm{R}}}$ should be interpreted as expansion and negative values as contraction.
In the text

Fig. 18. Two-dimensional maps of mean radial velocity for high-mass stars (M ≥ 2 M⊙, left panel) and low-mass stars (right panel). Stars belonging to the same Voronoi bin are color-coded according to the mean radial velocity, while black lines are iso-density contours of the respective populations enclosing about 12% (solid line), 39% (dashed line), and 67% (dash-dotted line) of the underlying density distribution.

In the text

	Fig. 19. Number density profiles (black symbols) for the nine star clusters. Also shown are the model predictions for the Plummer (red) and King (blue) models. Shaded areas represent the 68% credible regions. All the profiles, both observed and models, were scaled to the central density predicted by the Plummer model Σ0,  plummer and to the observed projected half-mass radius Rhm.
In the text

	Fig. 20. One-dimensional velocity dispersion profiles for seven out of nine star clusters. Velocity dispersion was normalized to the same quantity computed for all the cluster’s members ($\overline{{\sigma }_{1\mathrm{D}}}$, i.e., without binning the stars), while projected distances from the cluster center of mass were normalized to the observed half-mass radius.
In the text

	Fig. 21. Local number density of stars as a function of stellar mass (see Sect. 7.1). The local density is normalized by the median local surface density of the entire cluster. The increase in local density with stellar mass indicates mass segregation on local scales across all clusters shown.
In the text

	Fig. 22. Projected radial velocity profiles for each snapshot of the N-body simulation and the LISCA II system. The projected radius is normalized to the projected half-mass radius of the system in each snapshot, and the radial velocity is normalized to the radial velocity dispersion of all stars in the system.
In the text

Fig. 23. Surface density profiles (top) and 1D velocity dispersion profiles (bottom) for the TFF09 snapshot. The best-fit Plummer model for each clump is also shown (dashed lines). The best-fit Plummer model is determined by fitting the surface density profile within one projected half-mass radius. The surface density and 1D velocity dispersion are both normalized to the corresponding values in the innermost radial bin of each clump.

In the text

	Fig. 24. Radial profiles of the ratio of the clump radial profiles to the best-fit Plummer models for surface density (top) and 1D velocity dispersion (bottom). This ratio is averaged over the three largest clumps within each snapshot.
In the text

Fig. 25. One-dimensional velocity dispersion vs. projected surface density for the three largest clumps in each snapshot (see legend), each calculated at one projected half-mass radius. The 1D velocity dispersion and the projected surface density are both normalized by the value expected from the Plummer model best fitting the surface density profile. Point sizes are proportional to the total mass of each clump.

In the text

	Fig. A.1. Completeness distribution in four mass bins. Colors are consistent with Figure 15 and the mass range is shown in the upper left corner.
In the text