First Phase Space Portrait of a Hierarchical Stellar Structure in the Milky Way

It is widely accepted that most (70%–90%) stars in galaxies form in groups, clusters, or hierarchies, and spend some time gravitationally bound with their siblings while still embedded in their progenitor molecular cloud (Lada & Lada 2003; Portegies Zwart et al. 2010). Indeed, young stars tend to be grouped together into a hierarchy of scales where smaller and denser subregions (with a typical dimension of a few parsecs) are inside larger and looser ones (up to hundreds of parsecs; Elmegreen & Hunter 2010). The majority of such systems will be disrupted in their first few million years of existence, due to mechanisms possibly involving gas loss driven by stellar feedback (Moeckel & Bate 2010) or encounters with giant molecular clouds (Gieles et al. 2006). Nonetheless, a fraction of proto-clusters will survive the embedded phase and remain bound over longer timescales. The process of clustered star formation has major implications for many fundamental astrophysical areas of research including the star formation process itself (McKee & Ostriker 2007), the early interplay between stellar and gas dynamics and the consequences of gas expulsion for the cluster disruption (Baumgardt & Kroupa 2007), the possible formation of gravitational wave sources (Di Carlo et al. 2019), and the dynamical properties of young star clusters (McMillan et al. 2007). Cluster formation has also key implications for our understanding of the assembly process of galaxies in a cosmological context. Indeed, 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 therefore strictly connected with those of their hosts (Brodie & Strader 2006; Dalessandro et al. 2012). Indeed, massive star clusters may play a role in the formation of galactic substructures, and their partial or total dissolution contributes to the overall mass budget and stellar population properties of the hosts.

However, in order to efficiently exploit stellar clusters as tracers of galaxy and large-scale structure formation, it is essential to understand the physical processes setting their initial properties such as mass, structure, and chemical composition, and how they evolve across the cosmic time. In the last decade, many observational and theoretical studies have greatly enriched our understanding of the formation and early evolution of star clusters (e.g., Goodwin & Whitworth 2004; Allison et al. 2010; Parker et al. 2014; Adamo et al. 2015, 2020). However, questions concerning the possible presence of unifying principles governing the formation of different stellar systems are still unanswered (e.g., Bonnell et al. 2003; Banerjee & Kroupa 2014, 2015, 2017; Longmore et al. 2014; Kuhn et al. 2019, 2020). In addition, the presence in old globular clusters of multiple stellar populations characterized by specific abundance patterns in a number of light elements (see, for example, Bastian & Lardo 2018 and Gratton et al. 2019 for recent reviews) has indeed raised new key questions about the physical mechanisms at the basis of clusters' formation, and their dependence on the environment and the formation epoch (Krumholz et al. 2019). As a matter of fact, in spite of tremendous observational and theoretical efforts, our understanding of the star clusters' formation history and the underlying physical processes is still in its infancy.

The unprecedented precision and sensitivity in measuring stellar parallaxes and proper motions secured by the Gaia mission (Gaia Collaboration et al. 2018, 2020) now provide key tools to study in detail nearby star-forming regions and to use them as ideal laboratories to shed new light on our understanding of cluster formation and early evolution. Indeed, several papers have been published recently on the subject (e.g., Beccari et al. 2018; Zari et al. 2018; Kuhn et al. 2019; Meingast et al. 2019; Lim et al. 2020).

In this context, in the present study, we focus on the neighborhood of the well-known young double clusters h and χ Persei (NGC 884 and NGC 869, respectively). This area is located in the Galactic Perseus Arm, which includes massive star-forming regions W3–W4–W5, with two giant H ii regions (W4 and W5), a massive molecular ridge with active formation (W3), and several embedded star clusters and/or associations (Carpenter et al. 2000). The paper is organized as follows. The adopted data set is presented in Section 2; in Sections 3 and 4 the main spectro-photometric and kinematic properties of the area are described, respectively. A comparison with a set of N-body simulations following the violent relaxation phase and its subsequent evolution is described in Section 5. The main conclusions are discussed in Section 6.

For the present work we used photometric and astrometric information obtained from Gaia Early Data Release 3 (EDR3); we refer the reader to Gaia Collaboration et al. (2018, 2020) for details about the survey and available data. These data were complemented with high-resolution optical and near-infrared spectra obtained, respectively, with the HARPS-N (Cosentino et al. 2014) and GIANO-B (Oliva et al. 2012; Tozzi et al. 2016) spectrographs mounted on the Italian Telescopio Nazionale Galileo (TNG) operated on the island of La Palma at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias by the Fundacíon Galileo Galilei of the National Institute for Astrophysics (INAF). These observations are part of the TNG Large Program titled SPA—Stellar Population Astrophysics: the detailed, age-resolved chemistry of the Milky Way (MW) disk (Program ID A37TAC13, PI: L. Origlia), aimed at measuring detailed chemical abundances and radial velocities of the luminous stellar populations of the MW thin disk and its associated star clusters (Origlia et al. 2019). Within SPA we collected HARPS-N spectra at R ∼ 115,000 in the 378–691 nm range and GIANO-B spectra at R ∼ 50,000 in the 950–2450 nm range for 37 stars among blue and red supergiants and young main sequence stars in the analyzed area. HARPS-N spectra are automatically processed by the instrument data reduction software pipeline. GIANO-B spectra have been reduced using the GOFIO data reduction software code (Rainer et al. 2018).

2.1. Area Characterization and Cluster Identification

We selected stars in the Gaia-EDR3 catalog, distributed within a circular projected area with radius of 3 degrees from h and χ Persei (Figure 1), located at 1σ from their mean distance ($d={2.353}_{-0.197}^{+0.178}$ kpc, i.e., those stars having parallaxes in the range 0.385 mas < μ < 0.465 mas) and moving with similar mean velocities as the two clusters (${\mu }_{\alpha }\cos \delta =-0.64$ mas yr⁻¹ and μ_δ = −1.17 mas yr⁻¹) within a very strict tolerance range of ±0.4 mas yr⁻¹, corresponding to about 4 km s⁻¹ at their distance.

Zoom In Zoom Out Reset image size

As shown in Figure 1(a), in addition to h and χ Persei, other clusters, clumps, and elongated structures are clearly visible in the area after applying the described proper motion and parallax selections. Two of them are the known open clusters NGC 957 and Basel 10. To identify the additional cluster-like components, we used the clustering algorithm Density-Based Spatial Clustering of Applications with Noise (DBSCAN). DBSCAN is a density-based clustering algorithm that interprets clusters as regions of high density separated by areas of low density in space, without requiring any additional prior. Only two parameters need to be set in the algorithm, which are eps and min_samples. The parameter eps defines the radius of neighborhood around a point x, and min_samples represents the minimum number of neighbors within the eps radius. Higher min_samples or lower eps values indicate higher densities necessary to identify a subsystem as a cluster. Inspecting Figure 1, it is quite evident that clusters in the considered region have different sizes and numbers of members, and therefore different densities. For this reason, the clustering analysis needs to be run more than once. We first focused on the central region, where h and χ Persei are located. We applied DBSCAN on a subsample of stars having magnitude G < 16 and by imposing min_samples = 45 and eps = 0.04. This choice allowed us to isolate the main structures near the center of the map and to identify three main groups. Two of them are obviously h and χ Persei, the most dense and easy to identify, but a third component, very close to h Persei—that we will call hereafter N 1—has been also identified. Its location is highlighted in Figure 2. Then, we moved to the more external part, where less dense structures are present. We applied the DBSCAN algorithm only to such stars with the same magnitude cuts and using different parameters, min_samples = 41 and eps = 0.06, and we identified in this way four additional structures: the already known Basel 10 and NGC 957, and the two new additional systems that we name N 2 and N 3 (see Table 1 for a summary of the clusters' main properties). We note that the identification of the new clusters is quite solid against different settings of the DBSCAN search. In particular, we checked that results do not vary when using different eps values in the range 0.03–0.09.

Zoom In Zoom Out Reset image size

Table 1. Coordinates and Mean Velocities of the Clusters Identified in LISCA I

Cluster	R.A.(J2000)	Decl.(J2000)	${\mu }_{\alpha }\cos (\delta )$	μδ	LOS–V
	(hh:mm:ss)	(dd:mm:ss)	(mas yr−1)	(mas yr−1)	(km s−1)
χ Per	02:22:07.281	57:09:48.62	−0.621	−1.164	−45.3
h Per	02:18:58.429	57:07:43.01	−0.662	−1.170	−45.6
NGC 957	02:33:23.532	57:33:17.32	−0.293	−1.120	−44.8
Basel 10	02:19:23.301	58:17:25.58	−0.466	−0.881	−53.0
N 1	02:18:23.612	57:12:40.79	−0.682	−1.165	−47.1
N 2	02:29:52.466	58:08:19.30	−0.433	−0.923	−49.7
N 3	02:28:44.542	57:41:37.94	−0.437	−0.957	−46.2

Download table as:  ASCIITypeset image

The presence of these additional clusters is remarkably evident in the two-dimensional (2D) density map shown in Figure 1(b), which has been obtained by transforming the distribution of star positions into a smoothed surface density function through the use of a Gaussian kernel with width of $3^{\prime}$ (see, for example, Dalessandro et al. 2015). Interestingly, all star clusters in the complex show some degree of elongation, which is particularly relevant in the case of N 3. In addition, it is worth noting the presence of a significant overdensity bridging the regions between h and χ Persei and extending out to NGC 957 and N 2, which may be suggestive of ongoing tidal interactions or the remnants of primordial filamentary structures, which are now commonly found in the MW disk (Kounkel & Covey 2019). Figure 1(b) also shows the isodensity curves as obtained from 3σ to 50σ, where σ is the "background" density dispersion observed for stars located at about 3 degrees from h and χ Persei in the southeast quadrant. They clearly suggest the presence of a diffuse low-density stellar halo. The existence of a possible halo stellar distribution in the proximity of h and χ Persei (at distances < 0.5 ° from the two clusters) was already proposed by Currie et al. (2010). Thanks to the Gaia astrometric information, here we find that the halo actually extends with a rather irregular shape for at least 3 ° from h and χ Persei at a 3σ significance level.

Are these comoving clusters and diffuse halo (along with its substructures) part of a single system? Or are they just independent neighbors? The answers to these questions can have key implications for our understanding of their formation and evolution, and, more in general, on stellar cluster formation and early dynamical processes. To address them we will examine two main lines of evidence that will be detailed in Sections 3 and 4.

First, we derived stellar cluster and halo star chemical abundances by performing a detailed analysis of the high-resolution near-infrared spectra obtained with GIANO-B within the SPA-TNG Large Program. Detailed chemical abundances and abundance patterns for all of the most important metals are computed by using the MARCS model atmospheres (Gustafsson et al. 2008) and the spectral synthesis technique implemented in the Turbospectrum code (Alvarez & Plez 1998; Plez 2012) for the red supergiants in the sample. A detailed and comprehensive description of the analysis and results will be presented in a forthcoming paper (C. Fanelli et al. 2021, in preparation). For the purpose of this study, we report that all of the analyzed stars share the same metallicity within the uncertainties (which are typically of the order of 0.05 dex), with a mean value of [Fe/H] = − 0.29 ± 0.03 dex and about solar-scaled [α/Fe]. We compared these estimates with predictions from the Besançon Galaxy model (Robin et al. 2003) for stars in this region and selected in distance and proper motions as done before, and with results from surveys of the Galactic disk (Hayden et al. 2014; Mikolaitis et al. 2014) at the Galactocentric distance of h and χ Persei. The comparison distributions have a similar mean metallicity to that of stars analyzed here, but they are significantly broader, extending for about 1.5 dex in the range −1.0 < [Fe/H] < 0.3 dex and therefore suggesting that the lack of significant spread is not the expected outcome of the disk population properties in this area.

Second, we constrained cluster and halo ages by comparing the differential reddening corrected Gaia (G, B_p–R_p) color–magnitude diagrams (CMDs) with a suitable set of isochrones (Figure 3). Extinction values have been assigned to each star in the Gaia catalog by interpolation with the Schlegel et al. (1998) extinction maps and corrected to estimates by Schlafly & Finkbeiner (2011). An accurate age derivation is certainly a challenge for young and scarcely populated clusters (in particular for N 1, N 2, and N 3), as the turn-off region is significantly affected by low number counts (and relative fluctuations), stellar rotation, and binarity (Li et al. 2019), which can make the adoption of a direct isochrone fitting approach unreliable. However, we stress that here we are mainly interested in investigating the presence of significant age differences among clusters. Currie et al. (2010) found that h and χ Persei and their surrounding stars are coeval with an age t_age ∼ 14 Myr. We extended such an analysis by using a set of PARSEC isochrones (Bressan et al. 2012). We adopted the derived mean metallicity ([Fe/H] = −0.29) and distance (d = 2.353 kpc) and an average extinction value A_V = 1.65 mag. We found that all systems are compatible with being almost coeval with ages in the 15–25 Myr range represented by the three isochrones shown in each CMD of Figure 3. While uncertainties on the distance can surely impact the derived ages, we stress that the main source of uncertainty is the limited number of bright stars needed to anchor the best-fit isochrone. Interestingly, we find that also halo stars share the same age as clusters, in agreement with previous findings (Currie et al. 2010). Only Basel 10 is possibly slightly older (∼30 Myr).

Zoom In Zoom Out Reset image size

4.1. Kinematical Properties

The mean motion of each star cluster in the plane of the sky has been obtained by using Gaia-EDR3 proper motions of stars within confidence radii in the range of $16\buildrel{\,\prime}\over{.} 5\mbox{--}25^{\prime}$ from each clusters' center (Table 1) and by adopting Equation (2) from Gaia Collaboration et al. (2018). Figure 1(a) shows the derived cluster velocity vectors with respect to the entire system mean motion. All velocities were corrected for the effects of perspective expansion/contraction by using the procedures described in van de Ven et al. (2006), which, however, turn out to be negligible (of the order of a few one-hundredths of km s⁻¹) with respect to the typical stellar motions. By construction, all clusters' average motions reciprocally coincide within ∼3 km s⁻¹ in the sky component of the velocity vector. We then split the entire system in a 20 × 20 regular cell grid with each cell having a side of $16\buildrel{\,\prime}\over{.} 5$ and including only stars with G < 16. In each cell, we derived the average motions of stars with respect to the entire system mean motion. Also in this case, we properly accounted for the effects of perspective expansion/contraction. The resulting velocity map is shown in Figure 4(a), where an almost regular expansion pattern is clearly visible. Such a kinematical behavior might be due to a combination of effects associated with the dynamics of the violent relaxation phase and the dynamical response of the cluster to the mass loss from supernova explosions and gas expulsion, all of which can cause a significant cluster radial expansion. Indeed, a significant fraction of nearby young (t < 5 Myr) clusters and associations has been recently found to show a similar behavior based on Gaia DR2 data (Kuhn et al. 2019). We stress here that the expansion of a stellar system is a transient process not associated with a single evolutionary path and eventual fate. We further discuss this point in Section 4.2. We also note that, in the specific case of the analyzed region, this behavior can be in part linked to the so-called "Hubble Flow-like" large-scale motion that affects regions extending for several tens of degrees in the Perseus Arm and which has been observed (Román-Zúñiga et al. 2019) to be moving the Perseus Arm away from the W3–W4–W5 massive star formation zones with a velocity of 15 km s⁻¹ kpc⁻¹.

Zoom In Zoom Out Reset image size

The panels in Figure 4(b) show the smoothed density distributions of both the radial and tangential velocities as a function of the distance from the system barycenter. The radial component shows a quite narrow distribution slowly increasing as a function of the distance and reaching a peak of ∼2 km s⁻¹ at $R\sim 100^{\prime}$, thus confirming the expansion signal. The tangential velocity component shows a broader distribution, with a slowly declining trend as a function of the distance that can be indicative of a possible large-scale (of the order of some degrees) rotation pattern with an amplitude of 2–3 km s⁻¹ (corresponding to about 0.2–0.3 mas yr⁻¹).

To check for the line-of-sight velocity distribution, we complemented Gaia proper motions with LOS–Vs obtained with HARPS-N within the SPA-TNG Large Program for the 37 bright stars fulfilling the proper motions and parallax selection criteria described before (light blue symbols in Figure 1(a); see Section 3). Accurate (at better than 1 km s^-1) heliocentric radial velocities for the observed stars have been obtained by means of standard cross-correlation technique (Tonry & Davis1979) as implemented in the IRAF ¹⁰ (Tody 1986) task FXCOR and suitable synthetic templates computed with the SYNTHE code (Kurucz 1973). The median line-of-sight velocity of all the analyzed stars is −43.9 ± 3.5) km s^-1. The resulting LOS–V distribution of cluster and halo stars is shown in Figure 5. For each cluster, we obtained the mean velocity by averaging out the LOS–Vs of stars falling within $16\buildrel{\,\prime}\over{.} 5\mbox{--}25^{\prime}$ from the clusters' center. Our analysis shows that the clusters are moving at an average LOS–V_clusters ∼ −47.4 km s⁻¹ with σ ∼ 2.9 km s⁻¹. The stars (17) not directly associable with clusters are here tentatively assigned to the halo. The agreement between the line-of-sight motion of such stars and that of the cluster is satisfactory, as shown in Figure 5, with an average velocity LOS–V_halo ∼ −43.4 km s⁻¹ and σ ∼ 6.3 km s⁻¹. Figure 5 also shows a comparison with the LOS–V distribution obtained from the Besançon model for the same sample of stars used before. The observed LOS–V distribution is significantly narrower than the one expected for field disk stars within the same field of view, thus making it unlikely to be just a random extraction of disk star velocities in the area.

Zoom In Zoom Out Reset image size

4.2. Mass and Structural Properties

A detailed characterization of the individual clusters in the region will be the subject of a related paper (E. Dalessandro 2021, in preparation). Here, we provide further details on the diffuse stellar halo. To determine the projected density profile of the halo we used stars with G < 16 and selected in proper motions and parallaxes as described before. We also subtracted the contribution of each cluster by removing all stars located at a distance $r\lt 25^{\prime}$ from each cluster center. We divided the field of view in 16 concentric annuli centered on the system center and split in a number (from 2 to 4 depending on the azimuthal coverage) of subsectors. We then counted the number of stars lying within each subsector and divided it by the subsector area. The stellar density in each annulus was finally obtained as the average of the subsector densities, and the standard deviation among the subsector densities was adopted as the error. The resulting density profile is shown in Figure 6. The background-corrected profile is characterized by a typical profile with a high-density core and low-density halo that can be well fit by a single-mass King model (King 1966), thus suggesting that the diffuse stellar halo resembles a low-density cluster-like stellar structure. To determine the physical parameters of the halo we used single-mass, spherical, and isotropic King models (King 1966). This is a single-parameter family, where the shape of the density profile is uniquely determined by the dimensionless parameter W₀ (or concentration c). To determine the best-fit model we compared the background-subtracted surface density profile with the King model family and found the halo profile to be best fit by a low-density model (c = 0.84) with a large core and half-mass radii ${r}_{c}={1770.0}_{-194.7}^{+182.4}^{\prime\prime}$ and ${R}_{h}={3381.1}_{-304.5}^{+287.1}^{\prime\prime}$.

Zoom In Zoom Out Reset image size

We analyzed the mass function (MF) of the likely halo stars. Stellar masses were derived by interpolation as a function of the G-band magnitude with the 20 Myr isochrone used in Figure 3. Starting from the observed MF we derived the stellar mass of the entire system. We normalized both a Salpeter (1955) and a Kroupa (2001) theoretical initial mass function (IMF) to the high-mass portion (roughly corresponding to G < 12) of the observed MF, where the photometric completeness of Gaia is expected to be close to 100% in relatively low-density environments like those surrounding h and χ Persei. Then, the total mass is simply given by the integral of the normalized IMF in the mass range (0.1 < m < 8) M_⊙. We find that the total mass ranges from M_tot ∼ 6 × 10⁵ M_⊙, when a Salpeter IMF is adopted, to 1.5 × 10⁵ M_⊙ when a Kroupa IMF is used instead, making it compatible with the mass regime of the so-called Young Massive Clusters (YMCs) and with the peak of the present-day mass distribution of Galactic globular clusters. We find that the mass is equally split between clusters and the diffuse stellar halo. Stellar clusters have masses ranging from ∼ 2 × 10⁴ M_⊙ as for h and χ Persei (in quite good agreement with estimates from Currie et al. 2010) and ∼ 500 M_⊙ for N 3. By adopting the same heliocentric distance as before and a Galactocentric distance of 9.7 kpc, we find that the formal Jacobi radius of the system (which gives an approximate indication of the expected possible extension of the system) ranges from $\sim 140^{\prime}$ to $\sim 200^{\prime}$ (corresponding to 91–132 pc), depending on the adopted mass, thus nicely matching the observed properties of the analyzed structure.

We emphasize that the observed radial variation of the halo is not expected for a generic sample of field stars. The radial profile we found therefore provides further evidence that the halo identified in our analysis is one of the dynamical components of the complex cluster structure.

In addition, by combining all of the structural and kinematic information collected so far we can calculate both the one-dimensional (σ_1D) and the virial velocity dispersions (${\sigma }_{\mathrm{vir}}^{2}={GM}/\eta {R}_{h}$, see, e.g., Equation (19) in Kuhn et al. 2019), where σ_vir represents the velocity dispersion of a virialized stellar system. Kuhn et al. (2019) used these two quantities for an approximate determination of the dynamical state of 28 young (1–5 Myr) nearby clusters. Following Kuhn et al. (2019), we have adopted η ∼ 10 as appropriate for a stellar system with a Plummer star density distribution. We find that σ_1D = 1.54 ± 0.06 km s⁻¹ and σ_vir ranges from 1.29 km s⁻¹ to 2.63 km s⁻¹ for R_h = 33811 and assuming a total mass in the 1.5–6 × 10⁵ M_⊙ range. While these values provide only an approximate description of the system's dynamical state, they suggest it should be gravitationally bound.

The purely observational results collected so far show that (a) all stars and clusters distributed close to h and χ Persei define a quite complex structure with evidence of possible mutual interactions in the form of intra-cluster overdensities and/or bridges, (b) they have the same metallicity, α-element abundance, and age within a small interval, and (c) the stellar density distribution of the surrounding diffuse stellar halo resembles that of a cluster-like stellar system. The combination of these pieces of evidence clearly suggests that stars distributed within a few degrees from h and χ Persei are part of a common, substructured, and recently formed or still forming stellar system that we have named LISCA I. ¹¹

Several theoretical studies have explored the early dynamical evolution of star clusters and tried to establish a connection between the initial structural and kinematical properties of molecular clouds to the formation of single and multiple stellar clusters (see, for example, Goodwin & Whitworth 2004; Allison et al. 2010; Banerjee & Kroupa 2014, 2015; Parker et al. 2014; Fujii & Portegies Zwart 2016; Arnold et al. 2017; Sills et al. 2018). In order to provide an example of the evolutionary history behind a system with structural and kinematical properties similar to those found in LISCA I, we show here the results of one of the realizations of a large suite of direct N-body simulations aimed at investigating the violent relaxation phase of a cluster, starting from initial conditions characterized by nonvanishing total angular momentum and a fractal density distribution, to heuristically represent a turbulent giant molecular cloud in a variety of kinematic and structural initial states.

While a comprehensive description and phase-space analysis of this suite of N-body simulations, along with a comparison with previous works, will be presented in a separate article (A. L. Varri et al. 2021, in preparation), here we provide a short summary of the key properties of the survey. The initial conditions for these simulations were generated by using McLuster (Kuepper et al. 2011), with both homogeneous and inhomogeneous mass distributions. The values of the fractal dimension explored in the survey range from 1.6 to 3.0, with D = 3.0 corresponding to a homogeneous stellar system. As for the stellar velocities, the systems have been initialized with a variety of "temperatures" of the stellar distribution, with a ratio of the kinetic energy due to random stellar motions to the potential energy of the system ranging from 0.1 ("cold") to 0.25 ("warm"). With these initial conditions, the system is typically initially out of virial equilibrium, undergoes the so-called violent relaxation process (Aarseth et al. 1988), and eventually settles into an equilibrium configuration in a few dynamical times. Such a process naturally produces a characteristic pressure anisotropy signature (van Albada 1982; Trenti et al. 2005; Vesperini et al. 2014) in the velocity space of the resulting configuration. In addition, we have also explored the effects of the presence of some primordial angular momentum (Gott 1973; Boily et al. 1999). Internal rotation was added to the initial conditions such that the ratio between the kinetic energy due to ordered motions and the potential energy could assume a range of values from 0.0 (nonrotating) to 0.5 (moderately rotating), and the initial rotation curve is always approximately solid-body. All N-body simulations have N = 65,536 equal-mass particles and the simulations have been performed with the direct summation code Starlab (Portegies Zwart et al. 2001). The initial mass distribution of the simulation depicted in Figure 7 is characterized by a value of the fractal dimension parameter D = 2.4. Such a definition corresponds to significant deviations from a uniform spatial distribution of stellar positions. The initial kinematic state is "cold" (random kinetic energy/potential energy = 0.1) and moderately rotating (ordered kinetic energy/potential energy = 0.5).

Zoom In Zoom Out Reset image size

A representative example of the early dynamical evolution explored in such N-body simulations is illustrated in Figure 7 and can be summarized as follows: (a) at early stages the distribution of stars preserves the initial inhomogeneities and tens of small clumps ( ∼100 M_⊙) can be identified (Figure 7(a)). Several of these clumps are then destroyed because of dynamical interactions, and some of their stars contribute to the surrounding field and to the formation of a gravitationally bound halo (Figures 7(b), (c)); (b) at a later time, the surviving clumps merge to form a few more massive (10³–10⁴ M_⊙) and larger clusters (Figure 7(d), (e)); (c) finally, one or two clusters survive this hierarchical merger process and will eventually evolve into a single massive cluster (Figure 7(f)).

We emphasize that the snapshots shown here are not meant to provide a dynamical model tailored to the evolution of the LISCA I system, but rather illustrate the typical dynamical evolution of a hierarchical cluster assembly process in the presence of a realistic degree of kinematic complexity. Such complexity, which is often neglected in star cluster formation studies, has a key role in determining the properties of the violent relaxation phase, the longevity of the population of subclumps, and the structural and kinematic properties of the resulting end-product massive cluster. The six snapshots in Figure 7 show the system after about 1.4, 3.0, 4.6, 5.4, 6.7, and 12.6 freefall timescales (t_ff). Assuming an initial mass of M_i ∼ 6 × 10⁵ M_⊙ and a 3D half-mass–radius r_h ∼ 40 pc, the snapshot in Figure 7(b) would represent the structure of the system at the approximate age of the LISCA I system (20–25 Myr). Therefore, the observational properties of h and χ Persei and their surrounding structures fit well within an early stage of the hierarchical cluster assembly process (Figures 7(b), (c)), which may potentially lead the system to the formation of a stellar system of some 10⁵ M_⊙ within a short timescale (∼100 Myr; Figure 7(f)).

Also, the kinematical patterns corresponding to the early snapshots of Figure 7 (see Figure 8) nicely resemble those observed in LISCA I. In particular, we note that also the expansion pattern observed in LISCA I is compatible with the results from our suite of simulations and shown in the lower panel of Figure 8 to compare with Figure 4(b). As discussed in Section 4.1, although there might be different physical mechanisms responsible for the expansion of the system, the (qualitative) agreement with the simulations suggests that it is compatible with the dynamical evolution of a system in the early phases of a hierarchical cluster assembly.

Zoom In Zoom Out Reset image size

Finally, we point out that although our simulations explore the evolution of systems starting with a clumpy mass distribution similar to one that could arise from a hierarchical star formation process, structural properties like those observed in the current dynamical phase of LISCA I could also result from the growth of density fluctuations in a system emerging with a more homogeneous initial mass distribution from a monolithic formation event (see, e.g., Aarseth et al. 1988).

The results presented in this paper suggest that we might have caught in the act, for the first time in our "courtyard," an ongoing hierarchical cluster assembly process for which the structural and kinematical properties of the involved components are described in detail.

It is hard to make strong predictions about the final by-product of LISCA I and the fraction of mass that will be actually retained during its evolution, as it may depend on many factors including the nature and duration of the observed expansion. However, the results obtained here suggest that the formation of small stellar structures and their subsequent growth driven by dynamical interactions might have strongly contributed to shape the observed properties of LISCA I, thus possibly representing a viable process to form massive and long-lived stellar systems also in relatively low-density environments. This mechanism was long believed to work only in massive starburst galaxies, such as M83 (Bastian et al. 2011) and M51 (Chandar et al. 2011), but to be inefficient in local, lower density disk galaxies such as the MW or M31 (Krumholz et al. 2012).

If LISCA I will be able to continue its evolution as schematically shown in Figure 7, it can produce a YMC or, alternatively, given its extension (R_h ∼ 35 pc), it could also be identified as the progenitor of a so-called diffuse star cluster (also known as "faint fuzzy"; Larsen & Brodie 2000) that is characterized by a mass comparable to that of globular clusters while being significantly more extended. These systems have a strongly debated origin (Fellhauer & Kroupa 2002); they are observed in external galaxies (Huxor et al. 2005), but have remained elusive in the MW so far. In any case, these results may give a boost to the possibility of studying massive cluster formation by resolving individual stars in our neighborhood, thus constraining the physical mechanisms with a level of detail that cannot be achieved for distant systems. Interestingly in this respect, the properties of the interstellar medium in the Galactic disk are similar to those in the disks of other nearby galaxies (Longmore et al. 2013). As a consequence, our understanding of star formation and feedback in close young clusters and associations can probe star formation across much of the local universe.

Recent results obtained by Gaia (Cantat-Gaudin et al. 2018; Castro-Ginard et al. 2018) have shown that the Galactic disk structure is composed of large numbers of stellar groups, hierarchies, and young clusters, with most of them distributed along filamentary, string-like structures with similar properties to those observed in the Perseus Arm and in LISCA I (see, for example, Kuhn et al. 2019). This creates the possibility to check whether the hierarchical assembly process framed here can be efficient within disk-like galaxies such as the MW. Our study would therefore introduce a possible new picture in which the different appearance and properties of young stellar aggregates ranging from clumps and filaments as observed in the Galaxy spiral arms (e.g., Kounkel & Covey 2019; Kuhn et al. 2019), to double/multiple clusters as found in the Magellanic Cloud disks (e.g., Dieball et al. 2002; Mucciarelli et al. 2012; Dalessandro et al. 2018) and finally to massive globular clusters, can be interpreted as the morphological evidence of different evolutionary phases of the same hierarchical assembly process. This has important implications for our understanding of the environmental conditions (both locally and in the distant universe) necessary to form massive stellar clusters and, as a consequence, on the evolution of stellar cluster properties over cosmic time. In this respect, it is extremely interesting to note that the hierarchical build-up of massive stellar clusters may play an important role during the formation of the multiple stellar populations observed in old globular clusters and could be one of the key ingredients necessary to explain the variety of stellar population properties (Bekki 2017; Bailin 2018; Howard et al. 2019). Finally, the results obtained for LISCA I and the N-body simulations also show that the kinematic state of the initial molecular cloud plays a key role in prolonging the survival timescales of sparse and low-mass clumps and associations, thus mitigating the discrepancy between their observed and expected numbers in the Galaxy and providing further support to a primordial origin of the "kinematic complexity" currently emerging in old globular clusters (Kamann et al. 2018; Lanzoni et al. 2018).

The authors thank the anonymous referee for the careful reading of the paper and the useful comments and suggestions. E.D. and L.O. acknowledge financial support from the project Light-on-Dark granted by MIUR through PRIN2017-2017K7REXT contract and from the Mainstream project SC3K—Star clusters in the inner 3 kpc (1.05.01.86.21) granted by INAF. A.L.V. acknowledges support from a UKRI Future Leaders Fellowship (MR/S018859/1). M.B. acknowledges the financial support to this research by INAF, through the Mainstream Grant 1.05.01.86.22 assigned to the project Chemo-dynamics of globular clusters: the Gaia revolution. S.S. gratefully acknowledges financial support from the European Research Council (ERC-CoG-646928, Multi-Pop). M.F. acknowledges the financial support of the Istituto Nazionale di Astrofisica (INAF), Osservatorio Astronomico di Roma, and Agenzia Spaziale Italiana (ASI) under contract to INAF: ASI 2014-049-R.0 dedicated to SSDC. 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.

Software: DBSCAN, FXCOR (Tonry & Davis 1979), GOFIO (Rainer et al. 2018), IRAF (Tody 1986), MARCS (Gustafsson et al. 2008), McLuster (Kuepper et al. 2011), Starlab (Portegies Zwart et al. 2001), SYNTHE (Kurucz 1973), Turbospectrum (Alvarez & Plez 1998; Plez 2012).