3D Morphology of Open Clusters in the Solar Neighborhood with Gaia EDR 3. II. Hierarchical Star Formation Revealed by Spatial and Kinematic Substructures

Stars are formed in molecular clouds where the environment is compacted with dense gas (Lada & Lada 2003). When the gravitational force inside molecular clouds overcomes the internal pressure after the cloud cools down, the gas starts to collapse and form stars. This contraction process can result in the formation of filaments that form a network feeding new materials into hubs at the intersections (Schneider et al. 2010; Krause et al. 2020). This mass transfer along the filaments is referred to as conveyor belt (Longmore et al. 2014; Krumholz et al. 2019). Therefore, the interstellar medium inside the molecular cloud is highly hierarchically structured. Evidence from infrared and submillimeter observations shows that the pre- and protostellar cores' positions are strongly correlated with the hierarchy inside the clouds (Elmegreen & Falgarone 1996; Testi et al. 2000; Gutermuth et al. 2009; Könyves et al. 2015; Rathborne et al. 2015; Treviño-Morales et al. 2019). This observational evidence indicates that star formation is taking place across a continuous density distribution, from the densest hub regions to broader and lower-density filamentary networks (Treviño-Morales et al. 2019). The aforementioned observational results indicate a need for improvement of the monolithic cluster formation mechanism (Lada et al. 1984). In the monolithic cluster formation scenario, star clusters were first born in the central densest region of molecular clouds. After a violent gas expulsion, stars in the dense clusters will expand and eventually disperse into the field, i.e., the infant mortality (Lada & Lada 2003).

To explain the star formation taking place in the hierarchical interstellar medium, Kruijssen (2012) developed the theoretical framework of hierarchical star formation. In this paradigm, only star formation in the high-density region can produce gravitationally bound star clusters. Stellar groups formed in low-density regions have filamentary substructures, and quickly become unbound after the residual gas is expelled. This model predicted that one-third of the field star population come from bound clusters. Recent Gaia observations tend to support this scenario and suggest that only about 8%–27% of the stars in the solar neighborhood have originated from bound clusters (Anders et al. 2021).

Young stellar groups are the ideal targets to test the hierarchical formation theory. Their youth allows us to trace the morphological and kinematic features inherited from their parent molecular cloud (Vázquez-Semadeni et al. 2017; Wright & Mamajek 2018; Ward et al. 2020; Krause et al. 2020). Previous studies concluded that stars of OB associations are born in the low-density regions of molecular clouds, and that their filamentary substructures were formed in situ (Wright et al. 2016; Wright & Mamajek 2018; Ward & Kruijssen 2018; Ward et al. 2020). The elongated filamentary substructures of stellar groups, the relic of filaments, are often associated with young (≲100 Myr) open clusters (Jerabkova et al. 2019; Cantat-Gaudin et al. 2019b; Kounkel & Covey 2019; Beccari et al. 2020; Pang et al. 2021b; Kerr et al. 2021; Wang et al. 2022) where significant kinematic substructures are also found (Wright et al. 2016; Wright & Mamajek 2018; Pang et al. 2021b). Young open clusters still have the morphological and kinematic substructures that support the nature of their unevolved dynamical status. Evidence of initial substructure in both morphology and kinematics is erased within a few global crossing times (Goodwin & Whitworth 2004).

According to the conveyor belt mechanism in the framework of hierarchical formation, not only gas, but also stars are transferred into the hub regions. The infalling flow brings smaller groups of stars along the filaments to the dense hub region (Vázquez-Semadeni et al. 2017; Treviño-Morales et al. 2019). These infalling subgroups will continue to merge when the dynamics of the intergroups reach subvirial (Goodwin & Whitworth 2004). The morphology of the merged subgroups is then temporarily preserved. After sufficient relaxation in the same potential well of the hub region, subgroups mix and lose their memory of the infalling dynamics from the filament. However, if the dynamical state of the intergroups transits from subvirial to supervirial during the infalling process, subgroups will just pass by each other. Numerical simulations have shown that star clusters form from the merging of substellar groups at a dynamically subvirial stage (McMillan et al. 2007; Allison et al. 2009; Moeckel & Bonnell 2009). Moreover, multiwavelength observations of clusters showing evidence of subgroup mergers (Kuhn et al. 2015) support these simulation results. Nonetheless, more observations of young clusters (<100 Myr) are required to fortify (strengthen) the hierarchical formation mechanism for star clusters.

In our first paper of this series in studying the 3D morphology of open clusters in the solar neighborhood (Pang et al. 2021a, hereafter Paper I), we found three out of 13 target clusters already show signs of elongated filamentary substructures, i.e., young clusters NGC 2232, NGC 2547, and NGC 2451B with an age of about 25–60 Myr. These filamentary substructures may be connected to a larger hierarchical network, e.g., NGC 2232 is embedded in a string structure of a few hundred parsec long (Kounkel & Covey 2019). Our follow-up study (Pang et al. 2021b) investigated five hierarchical clustering groups in the Vela OB2 region that showed a shell-like morphology related to kinematic substructures in the proper motion (PM) space. The significant expansion found in these hierarchical groups probably leads to the fate of dissolution. It is desirable to further search for more hierarchical groups and young clusters harboring spatial and kinematic substructures to understand the early formation process and evolution of star clusters. In this study, we investigate the morphology and kinematics of 85 open clusters (Table 1) located in the solar neighborhood primarily within 500 pc based on Gaia Early Data Release 3 (EDR 3) data (Gaia Collaboration et al. 2021). In our samples, 60 out of the 85 clusters are younger than 100 Myr. We follow the method developed in Paper I to quantify stellar groups' spatial distribution and kinematic distribution in either three- or two-dimensional space.

Table 1. Parameters of 85 Open Clusters

Cluster	R.A.	Decl.	Dist.	RV	${\mu }_{\alpha }\cos \delta$	μδ	Age	Mcl	rh	rt	N	fm	Flag
	(deg)		(pc)	(km s−1)	(mas yr−1)		(Myr)	(M⊙)	(pc)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)
Alessi 3	109.139147	−46.074168	277.3	1.22	−9.804	11.945	631	124.6	6.8	7.0	165	0.52	t
Alessi 5	160.927572	−61.087405	397.7	11.13	−15.354	2.535	52	244.0	4.4	8.7	279	0.89
Alessi 9	265.946512	−47.070299	207.6	−6.59	10.036	−9.109	265	61.7	4.5	5.5	79	0.62
Alessi 20	2.970330	58.532158	423.4	−5.74	7.904	−2.378	9	240.4	11.1	8.7	349	0.32	f1
Alessi 20 gp1 a	356.035534	56.078582	411.9	−10.68	8.130	−3.347	12	173.2	13.6	7.8	265	0.18	f2
Alessi 20 isl1 a	11.906507	50.262109	458.9	−3.21	7.391	−3.108	100	97.6	7.5	6.4	107	0.41
Alessi 24	260.842675	−62.670437	484.2	12.10	−0.441	−8.948	88	119.1	4.6	6.9	138	0.73
Alessi 62	284.025507	21.581477	618.9	12.77	0.256	−1.120	692	144.1	5.7	7.3	139	0.63
ASCC 16	81.160907	1.558096	346.4	20.95	1.299	−0.013	11	241.7	6.0	8.7	373	0.74
ASCC 19	82.011177	−1.975022	354.7	24.50	1.118	−1.173	9	198.5	7.8	8.2	323	0.53	f2
ASCC 58	153.712881	−55.126490	477.1	11.69	−13.353	2.703	52	254.9	9.9	8.9	328	0.44	f2
ASCC 105	295.444115	27.493132	557.8	−15.46	1.426	−1.602	74	67.5	5.9	5.7	67	0.48	f1
ASCC 127	346.610023	65.131092	374.4	−11.89	7.424	−1.785	15	183.0	11.9	7.9	260	0.29	f1
BH 99	159.413798	−59.062415	446.9	13.08	−14.470	1.089	81	556.7	8.3	11.5	730	0.63	f2
BH 164	222.290304	−66.446488	419.8	−0.03	−7.404	−10.689	65	194.5	5.3	8.1	243	0.72
Blanco 1 b	1.536331	−29.975190	236.4	6.21	18.708	2.606	100	342.9	6.7	10.2	703	0.64	t
Collinder 69	83.809113	9.809913	399.0	25.89	1.300	−2.113	14	401.8	8.1	10.3	620	0.63	f1
Collinder 135 b	109.553612	−37.079084	302.8	16.12	−10.101	6.165	40	253.2	9.7	8.9	377	0.41	f2
Collinder 140	110.944634	−32.063703	384.5	19.07	−8.073	4.769	50	179.6	10.2	7.9	241	0.42	f2
Collinder 350	267.003004	1.452913	367.8	−14.79	−4.933	−0.060	589	149.4	6.6	7.4	156	0.51	t
Coma Berenices b	186.034554	25.539697	85.6	0.04	−12.018	−8.940	700	101.6	4.7	6.8	158	0.59	t
Group X	217.766195	54.432658	99.6	−6.50	−16.051	−2.802	400	99.5	14.5	6.8	187	0.11	d
Gulliver 6	83.207952	−1.711134	413.3	30.27	0.128	−0.269	7	168.5	7.3	7.7	284	0.54	f1
Gulliver 21	106.889853	−25.471498	652.5	41.25	−1.875	4.226	275	83.3	4.9	6.1	74	0.68
Huluwa 1 b	122.390262	−47.017897	354.7	17.00	−6.374	9.360	12	724.0	14.6	12.6	1294	0.39	f1
Huluwa 2 b	121.163632	−48.884799	398.8	20.65	−5.525	8.226	11	467.4	15.3	10.9	743	0.26	f1
Huluwa 3 b	117.397635	−46.600642	394.7	21.76	−4.706	8.960	11	372.8	7.9	10.1	588	0.66
Huluwa 4 b	125.994565	−41.178685	341.8	19.83	−7.115	10.021	10	180.7	14.8	7.9	347	0.28	f1
Huluwa 5 b	126.933768	−34.930350	355.0	15.69	−7.007	10.779	8	60.8	4.7	5.5	102	0.65	f1
IC 348	56.093028	32.173637	316.5	16.99	4.415	−6.407	5	142.2	4.1	7.3	211	0.73
IC 2391 b	130.229784	−52.990392	151.3	15.01	−24.641	23.309	50	140.2	2.5	7.6	219	0.99
IC 2602 b	160.515429	−64.443674	151.4	16.57	−17.691	10.695	45	188.1	3.7	8.4	318	0.91
IC 4665 b	266.568565	5.608075	347.4	−13.19	−0.845	−8.543	36	158.5	6.0	7.9	197	0.72
IC 4756	279.617685	5.430164	473.7	−24.37	1.236	−4.991	955	508.0	8.2	11.2	497	0.70	t
LP 2371	82.056845	1.632299	367.0	29.36	−0.549	0.746	20	81.0	5.1	6.1	103	0.67
LP 2373	83.844838	−5.670678	386.9	21.01	1.621	−1.097	9	98.0	7.9	6.5	153	0.39	f1
LP 2373 gp1 a	81.170979	−2.515377	335.6	22.52	1.104	−0.277	11	187.4	12.0	8.0	341	0.23	f1
LP 2373 gp2	81.836831	1.892265	349.3	17.77	1.543	−0.545	9	543.0	11.6	11.4	884	0.47	f1
LP 2373 gp3	83.509434	−0.498267	349.1	20.10	1.671	−0.932	6	111.3	9.0	6.7	177	0.34	f1
LP 2373 gp4	84.260923	−2.053777	363.0	20.67	1.708	−1.310	6	296.2	9.9	9.3	490	0.48	f2
LP 2383	95.335972	−16.218225	364.3	24.81	−5.420	5.068	50	280.7	12.8	9.2	475	0.34	f2
LP 2388	127.066134	−47.931643	497.3	24.94	−5.923	6.884	200	149.4	6.8	7.4	199	0.53
LP 2428	43.296942	68.875038	436.1	−1.96	1.428	−7.976	200	111.0	6.0	6.7	142	0.58
LP 2429	84.460533	57.169259	479.5	−6.36	−3.222	−3.992	1150	148.0	6.1	7.4	179	0.61	t
LP 2439	103.540197	−5.896183	283.9	23.13	−7.336	−2.457	25	143.0	6.5	7.3	243	0.54	f2
LP 2441	279.297560	−14.225387	280.4	−21.37	−1.772	−9.437	75	187.8	8.9	8.0	226	0.46	f2
LP 2442	250.101032	−39.560185	175.8	1.45	−11.929	−21.288	15	318.7	6.2	9.6	662	0.70	f2
LP 2442 gp1 a	242.331352	−22.907537	138.9	−7.47	−11.731	−24.786	8	111.6	8.7	6.7	220	0.29	f2
LP 2442 gp2 a	243.556287	−23.095881	141.1	−6.73	−8.884	−24.562	8	151.8	7.9	7.5	363	0.47	f2
LP 2442 gp3 a	241.119806	−22.113086	143.4	−6.21	−12.199	−23.174	8	64.0	6.4	5.6	136	0.38	f2
LP 2442 gp4 a	244.532627	−25.119324	153.7	−4.34	−10.774	−21.905	8	113.0	10.2	6.8	284	0.26	f2
LP 2442 gp5 a	242.398266	−20.348396	153.9	−5.61	−9.609	−21.354	8	76.9	9.8	6.0	163	0.26	f2
Mamajek 4	276.232342	−51.301191	449.9	−28.90	4.413	−21.456	372	281.5	17.7	9.2	320	0.21	t
NGC 1901	79.555700	−68.181508	417.8	0.67	1.681	12.694	850	124.8	6.7	7.0	144	0.52
NGC 1977	83.847140	−4.973568	392.3	25.28	1.120	−0.547	3	107.7	4.7	6.7	202	0.76
NGC 1980	83.804049	−5.865419	383.6	23.50	1.119	0.344	6	753.2	7.6	12.7	1298	0.81
NGC 2232 b	96.721725	−4.630480	320.0	26.37	−4.626	−1.851	25	205.8	6.8	8.6	281	0.57	f1
NGC 2422 b	114.122352	−14.484943	476.4	35.86	−7.051	1.026	73	480.2	4.9	11.4	466	0.90
NGC 2451A b	115.754528	−38.264155	192.6	22.62	−20.954	15.333	58	182.2	4.9	8.3	311	0.82
NGC 2451B	116.079334	−38.016157	363.2	16.65	−9.575	4.846	50	276.1	8.0	9.1	400	0.54	f2
NGC 2516 b	119.492636	−60.786856	410.6	24.22	−4.620	11.203	123	1973.3	7.9	18.3	2690	0.88	t
NGC 2547 b	122.490342	−49.178877	387.2	13.42	−8.562	4.365	40	303.9	5.6	9.8	452	0.71	f2
NGC 3228	155.377537	−51.880903	482.2	16.49	−14.848	−0.573	63	84.4	6.3	6.1	81	0.50
NGC 3532	166.388592	−58.701886	478.3	5.44	−10.400	5.235	398	2210.5	7.4	18.2	2559	0.90	h
NGC 6405	265.095354	−32.227742	457.6	−8.12	−1.383	−5.859	79	596.9	5.3	11.8	654	0.87
NGC 6475	268.466733	−34.828885	279.5	−14.55	3.083	−5.403	186	1024.5	6.1	14.1	1157	0.86	h
NGC 6633 b	276.890553	6.680996	394.0	−28.11	1.242	−1.800	426	337.3	5.3	10.2	300	0.82
NGC 6774 b	289.091567	−16.317354	306.4	42.24	−0.889	−26.656	2650	152.6	4.7	7.8	154	0.76
NGC 7058	320.451323	50.821327	365.2	−19.21	7.483	2.772	80	126.8	3.4	7.0	190	0.83
NGC 7092	322.782191	48.247177	297.1	−5.29	−7.387	−19.718	350	193.1	6.5	8.1	303	0.61	t
Pleiades b	56.637062	24.132371	135.9	5.72	19.945	−45.365	125	744.0	4.7	12.7	1407	0.86	t
Praesepe	130.008023	19.601297	185.0	35.00	−35.904	−12.864	700	601.8	4.7	11.8	982	0.92	h
Roslund 5	302.659299	33.751082	540.5	−16.61	2.078	−1.183	98	191.5	6.8	8.1	207	0.57	f2
RSG 7	343.949646	59.772583	419.3	−8.55	4.760	−2.244	70	67.6	4.6	5.7	69	0.59
RSG 8	345.689844	59.113059	474.6	−8.68	5.665	−1.663	18	342.8	22.0	9.8	408	0.09	f2
Stephenson 1	283.406739	36.720713	358.2	−19.13	1.075	−3.020	46	263.4	7.0	9.0	407	0.58	f1
Stock 1	294.146449	25.145437	406.3	−19.75	6.105	0.205	470	136.3	4.6	7.2	125	0.77
Stock 12	353.914758	52.678590	437.7	−2.18	8.553	−1.920	112	122.0	4.3	6.9	124	0.73
Stock 23	49.134918	60.355183	606.7	−18.57	−4.313	−1.010	94	106.0	9.1	6.6	78	0.37	f1
UBC 7 b	106.935171	−37.677240	278.3	16.78	−9.730	7.004	40	192.3	7.6	8.1	336	0.54	f2
UBC 19	56.443686	29.936422	399.7	19.45	2.825	−5.217	7	42.2	6.1	4.9	57	0.38
UBC 31	60.707033	32.542298	365.0	20.89	3.698	−5.370	12	260.0	14.7	8.9	353	0.13	f1
UBC 31 gp1 a	63.563584	32.027203	339.2	20.90	3.766	−5.937	12	58.9	10.7	5.4	80	0.13	f1
UBC 31 gp2 a	57.959044	35.028000	381.4	23.12	3.385	−4.450	10	185.1	10.7	8.0	213	0.27	f1
UPK 82	298.265305	26.504003	542.2	−12.72	2.258	−2.117	81	57.6	4.8	5.4	61	0.56

Notes.

^a New hierarchical groups identified in this work. ^b Clusters members taken from Pang et al. (2021a, 2021b); Li et al. (2021).

Columns 2–7 list the median value of cluster members. R.A. and decl. are the R.A. and decl. Dist. is the distance after correction in Section 3. RV is the radial velocity. ${\mu }_{\alpha }\cos \delta$and μ_δ are the components of the PM. The age of the cluster is derived from PARSEC isochrone fitting. M_cl is the total mass of each star cluster. r_h and r_t are the half-mass and the tidal radii of each cluster, respectively. N is the total number of members, and f_m is the bound mass fraction, which is the ratio of mass within the tidal radius. The flag in Column 14 indicates the morphological type of the cluster substructures outside r_t: "f1" for filamentary, "f2" for fractal, "h" for halo, and "t" for tidal tail (see Section 4.1 for a further description). Group X is instead labeled as "d" for its total disruption state (Tang et al. 2019).

Download table as:  ASCIITypeset images: 1 2

This paper is organized as follows. In Section 2, we describe steps on the preprocessing of Gaia EDR 3 data that serve as input for cluster member selection. We then present the algorithm, StarGO, used to determine cluster members and identify hierarchical clustering groups. The reliability of the cluster membership is verified in Section 2.3. In Section 3, the distance correction for individual stars is carried out. The 3D morphologies and classification of the target open clusters are presented in Section 4. The analysis of their kinematic substructures is given in Section 5. The properties of the hierarchical groups are presented in Section 6. The morphology of the stars inside the tidal radius of each cluster is quantified in Section 7. In Section 8, we discuss the dependence of cluster morphology and its indication on star formation and evolution. Finally, we provide a summary in Section 9.

To increase the young open cluster sample size in Paper I, we apply the search for hierarchical stellar grouping in 50 open cluster regions in the solar neighborhood (up to a distance of about 650 pc). Unlike Paper I, where we only focus on studying the 3D morphology of the 13 clusters, and ignore stellar groups in the nearby sky region, we include all nearby closely related structures of the target region in this study, e.g., hierarchical clustering groups or neighboring clusters.

2.1. Gaia EDR 3 Data Preprocessing

The data preprocessing procedures prior to the member selection (Section 2.2) are similar to our previous studies (Tang et al. 2019; Pang et al. 2020, 2021a, 2021b). In short, we query the Gaia EDR 3 database for the regions around 50 target regions using a spherical spatial cut in the Cartesian coordinates. A target region can contain either just one open cluster or several clusters and/or stellar groups at the same time. The radius of the spherical cut is either 100 pc or 150 pc from the cluster center, depending on the size of the stellar groups present in the region. We then apply a PM cut for each cluster region based on a 2D PM density map (e.g., Figure 1 in Paper I). The PM cuts are done circularly to include as many potential members in the target region as possible. The average spatial coordinates and PMs of clusters used for the above selection are taken from Liu & Pang (2019) and Cantat-Gaudin et al. (2020). The data quality cut for the Gaia data follows the instruction in the Appendix C of Lindegren et al. (2018) to select stars with parallaxes and photometric measurements less than 10% uncertainty.

Zoom In Zoom Out Reset image size

2.2. Membership Determination

We use StarGO ⁸ (Yuan et al. 2018), an unsupervised machine-learning software based on the Self-Organizing-Map (SOM) algorithm, to select cluster members in data after the PM cut. StarGO has been successful in membership identification in both open clusters (Tang et al. 2019; Pang et al. 2020, 2021a, 2021b) and stellar streams (Yuan et al. 2020a, 2020b). Detailed descriptions about StarGO were given in the studies mentioned above; thus, here, we only present a brief summary of StarGO. The SOM algorithm starts with building a 2D neuron map with each neuron initially assigned with five random weight vectors. The dimension of the random weight vectors matches the number of the input parameters, which are X, Y, Z, ${\mu }_{\alpha }\cos \delta$, and μ_δ for this study. We adopt a network of either 100 × 100 or 150 × 150 neurons, depending on the number of input stars. During the learning process, target stars are fed to the 2D neuron map one by one so that the weight vectors of each neuron can be updated to become closer to the input vectors of a given star. One iteration of the training cycle is finished after the neurons are trained by every star in the sky region. For the weight vectors in all neurons to converge, the training cycle is set to iterate 400 times. After the training process of the 2D neuron map, each star will be assigned to a neuron where the five parameters of the star are closest to the neuron's weight vectors.

A cluster of stellar objects will appear as a grouping in the 5D input parameter space. In terms of the trained 2D neuron map, this grouping will show as a local minimum of the u values—the difference in weight vectors between adjacent neurons. The final cluster members are stars associated with neurons below a certain u value, which is determined by a ∼5% field star contamination rate (e.g., Figures 1(b) and (c) in Pang et al. 2020). The estimation of the field star contamination rate is computed with stars from the mock Gaia EDR 3 catalog (Rybizki et al. 2020). After applying the same selection criteria described in Section 2.1 to the mock stars, those remaining mock stars attached to the neurons same as the cluster members are considered as false positive, i.e., field star contaminants. Assuming the smooth Milky Way population has the same properties as the mock catalog, the contamination rate of each identified group can be estimated as the number of the field star contaminants over the selected cluster members (more details in Pang et al. 2020, Section 2.2).

For regions showing significant substructures in the 2D neural network, a sign of hierarchical clustering, we apply a top-down subgroup selection method developed in (Pang et al. 2021b, detailed in Section 2.2 and Figure 1(c)) to further disentangle them. In this top-down subgroup selection method, the top-level clustering is the 2D neural patch regions associated with u values corresponding to 5% contamination. To reveal the hierarchical substructures, we decrease the threshold value of u until the contamination rate reaches 1%–2%. The remaining patches are considered as the core regions of each subgroup. We then associate neurons in the top-level structures to these cores by calculating the minimum difference in weight vectors between a given neuron to each core patch.

An additional member selection process is carried out in the color–magnitude diagram. Similar to Section 2.3 in Pang et al. (2020), we first fit the PARSEC isochrone (version 1.2S, Bressan et al. 2012; Chen et al. 2015) to all clusters and stellar groups by eye with reddening and metallicity values adopted from the literature (Villanova et al. 2009; Zari et al. 2017; Gaia Collaboration et al. 2018a, 2018b; Bossini et al. 2019; Zhang et al. 2019; Carrera et al. 2019; Röser & Schilbach 2019; Casali et al. 2020; Pang et al. 2021a). Stars below the best-fit isochrone or an isochrone with older age for young clusters (<10 Myr) were then removed. The mass of each member star is then estimated from the best-fit isochrone using the K-D tree method (McMillan et al. 2007) to find the nearest point on the isochrone.

Finally, there are 65 clusters and stellar groups, with a total of 22,347 members, identified in the original 50 target regions. We provide basic information of these 65 clusters in Table 1 obtained in this study and their complete member lists in Table 2. Among the identified 65 clusters, 36 crossmatch with clusters in Liu & Pang (2019) and 52 crossmatch with clusters in Cantat-Gaudin et al. (2020). On average, the matched stars account for 70%–80% of members in Liu & Pang (2019) and Cantat-Gaudin et al. (2020). In general, we double the number of member stars in these matched clusters compare to previous studies (Liu & Pang 2019; Cantat-Gaudin et al. 2020). Among the 29 hierarchical groups identified in our samples, one-third (10 groups) are newly discovered; these are located in regions of Alessi 20, IC 348, LP 2373, and LP 2442.

Table 2. Columns of the Table of Individual Members of the 65 Clusters Identified in This Work

Column	Unit	Description
Cluster name	⋯	Name of the target cluster
Gaia ID	⋯	Object ID in Gaia EDR 3
ra	degree	R.A. at J2016.0 from Gaia EDR 3
er_RA	mas	Positional uncertainty in R.A. at J2016.0 from Gaia EDR 3
dec	degree	Decl. at J2016.0 from Gaia EDR 3
er_DEC	mas	Positional uncertainty in decl. at J2016.0 from Gaia EDR 3
parallax	mas	Parallax from Gaia EDR 3
er_parallax	mas	Uncertainty in the parallax
pmra	mas yr−1	PM with robust fit in $\alpha \cos \delta$ from Gaia EDR 3
er_pmra	mas yr−1	Error of the PM with robust fit in $\alpha \cos \delta$
pmdec	mas yr−1	PM with robust fit in δ from Gaia EDR 3
er_pmdec	mas yr−1	Error of the PM with robust fit in δ
Gmag	mag	Magnitude in the G band from Gaia EDR 3
BR	mag	Magnitude in the BR band from Gaia EDR 3
RP	mag	Magnitude in the RP band from Gaia EDR 3
Gaia_radial_velocity	km s−1	RV from Gaia DR 2
er_Gaia_radial_velocity	km s−1	Error of radial velocity from Gaia DR 2
Mass	M⊙	Stellar mass obtained in this study
X_obs	pc	Heliocentric Cartesian X coordinate computed via direct inverting Gaia EDR 3 parallax ϖ
Y_obs	pc	Heliocentric Cartesian Y coordinate computed via direct inverting Gaia EDR 3 parallax ϖ
Z_obs	pc	Heliocentric Cartesian Z coordinate computed via direct inverting Gaia EDR 3 parallax ϖ
X_cor	pc	Heliocentric Cartesian X coordinate after distance correction in this study
Y_cor	pc	Heliocentric Cartesian Y coordinate after distance correction in this study
Z_cor	pc	Heliocentric Cartesian Z coordinate after distance correction in this study
Dist_cor	pc	The corrected distance of individual member

Note. A machine-readable version of this table is available online.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

To further enlarge our sample size for a more robust statistical study on the 3D morphology and kinematic structures of nearby clusters and groups, we add 20 additional clusters with members from our previous studies. Twelve clusters are from Paper I: IC 2391, IC 2602, IC 4665, NGC 2422, NGC 2516, NGC 2547, NGC 6633, NGC 6774, NGC 2451A, NGC 2232, Blanco 1, and Coma Berenices, seven stellar groups (clusters) from the Vela OB2 region in (Pang et al. 2021b, Huluwa 1–5, Collinder 135, and UBC 7), and the final one, the Pleiades open cluster, from Li et al. (2021). Although a list of members of NGC 2415B was already provided in Paper I, we re-select its member in this study (increased by 10%) in a larger spherical space to better probe its extended substructures. Parameters of the abovementioned 20 extra clusters are obtained from previous studies and included in Table 1, except for the bound mass fraction (f_m ) and the morphological type (flag) that are determined in this work. Therefore, we use 85 open clusters for analysis in this study. Although some of the target cluster members come from other studies, they were all obtained with the same procedure of member identification in the current study, and therefore form a homogeneous data set that is ideal for statistical analysis.

2.3. Membership Verification

An apparent stretching along the line of sight in open clusters results from obtaining stellar distances from inverting the parallax; this phenomenon has been well studied (e.g., Carrera et al. 2019; Zhang et al. 2019; Pang et al. 2020, 2021a; Tarricq et al. 2022; Piecka & Paunzen 2021). Pang et al. (2020, 2021a) attempted to mitigate this problem via a Bayesian approach described in Carrera et al. (2019). This Bayesian approach is based on the assumption of a normal distribution for the individual distances to the cluster stars and an exponentially decreasing profile for the distances to field stars (Bailer-Jones 2015). To understand the limitation of this correction method, Paper I carried out a Monte Carlo (MC) simulation and showed that the uncertainty in the corrected distances, σ_dis, increases monotonically with distance. ¹⁰ In the MC simulations, the σ_dis reaches an accepted value of 3.0–6.3 pc at a distance of 500 pc under two initial conditions: (1) a spherical uniform distribution of cluster members, and (2) a cluster with members artificially elongated (perpendicular or along the line of sight). We adopt the same Bayesian procedure to correct the distances of all the members of the 65 clusters in this work. The corrected distances and the 3D positions are listed in Table 2, together with other parameters from Gaia EDR 3. Group X is excluded from distance correction because of its fragmented spatial distribution, and thus, the spatial positions of its members are directly computed from inverse Gaia parallax.

In Figure 2, we present the σ_dis (blue dots) for all 84 open clusters (12 taken from Paper I are highlighted with black open circles). The value of σ_dis resides mostly between the predicted value of uniform and elongated models (solid and dashed–dotted curves), depending on the quality of the parallax. Clusters at a further distance suffer from a reduced parallax quality. However, the value of the tidal radius of each cluster is typically three times larger than the uncertainty in the corrected distance σ_dis, ensuring an accurate determination of the intrinsic cluster morphology within the tidal radius.

Zoom In Zoom Out Reset image size

We apply the distance correction to the 84 simulated clusters in order to validate the reliability of the distance correction. Within the 84 simulated clusters, 18 develop into having tidal tails that stretch out beyond a few hundred parsec. Fewer than 10 stars remain in the bound region (i.e., within tidal tails) in those 18 clusters, a sign of ongoing dispersal. Therefore, we only correct for clusters with more than 10 members within r_t. The difference between the corrected position and the theoretical position from petar is defined as the σ_dis. Together with σ_dis from observation data, we present the σ_dis of simulated clusters (olive dots) in Figure 2. The σ_dis for simulated clusters is, in general, smaller than that of the observations, with a few outliers distributed in the region occupied by the observations. This result assures us that the distance correction is appropriate.

In Figure 3, we present the 3D spatial distributions for member stars in the 85 target clusters after distance correction. Our cluster samples have a larger variety in cluster total masses. The difference can be as high as two orders of magnitude, from the lowest-mass cluster, UBC 19, ∼42 M_⊙ (57 members), to the most massive cluster, NGC 3532, ∼2210 M_⊙ (2559 members). Each dot in Figure 3 represents a star, with the dot size proportional to the stellar mass (ranges from 0.1–13.7 M_⊙). An online interactive 3D demonstration is available for Figure 3 ¹¹ , through which stellar mass distribution for each cluster can be observed. The age of each cluster in Figure 3 is scaled by the color bar, with the blue color representing the youngest populations in our sample and red for the oldest. Most clusters (groups) in our sample are a few to a hundred million years old. These young clusters are closely associated with the center of the Local Arm (solid pink curve), where massive star formation takes place (Xu et al. 2013, 2016). On the contrary, older generations, are approaching the edge of the Local Arm (pink dashed curve).

Zoom In Zoom Out Reset image size

To further inspect the detailed spatial distribution of the target clusters, we project them onto X–Y, X–Z, and Y–Z planes in Figure 4, where we indicate the young cluster (group) regions and typical old clusters. The clustering effect of young clusters around the center of the Local Arm (solid pink curve in Figure 4(a)) is better observed. Younger clusters, especially those less than 30 Myr (blue), are distributed not only closer to the center of the Local Arm but also closer to the Galactic plane (Figure 4(b)). On the other hand, older clusters (yellow to red, >100 Myr), such as the two oldest NGC 6774 (2.65 Gyr) and LP 2429 (∼1.15 Gyr determined in this study), show a preference for being further away from the arm and disk. The spatial positions of young clusters on the Y–Z plane (Figure 4(c)) appear as a wave pattern. This wave pattern resembles to some extent the "Radcliffe" wave found by Alves et al. (2020) along the Y-axis, which is only associated with young clusters (Kounkel 2020). The distribution of young clusters (<100 Myr) in our samples on the Y–Z plane significantly differs from a uniform distribution (a p-value of 0.02 from a K-S test), implying the observed wave pattern is not a random event. Although our cluster sample in the solar neighborhood is the tip of the iceberg, the spatial preference of young and old clusters is consistent with that of full-sky cluster samples identified in earlier studies (Soubiran et al. 2018; Cantat-Gaudin et al. 2020; Castro-Ginard et al. 2020; Hao et al. 2022).

Zoom In Zoom Out Reset image size

4.1. Morphological Substructure Classification

In this work, we consider the extended features outside the tidal radius of each cluster as their morphological substructure. We qualitatively classify them into four types:

• 1.  
  Filamentary (f1 type): young clusters (groups) <100 Myr with unidirectional filaments that is mostly elongated along one direction.
• 2.  
  Fractal (f2 type): young clusters (groups) <100 Myr with multidirectional fractal substructures.
• 3.  
  Halo (h type): clusters with age >100 Myr having a compact core but show some low-density stars spread out in the outskirt.
• 4.  
  Tidal tail (t type): clusters older than 100 Myr show unidirectional tidal tails.

The classification results are listed in Column 14 of Table 1 and are discussed below. We label Group X as "d" in Table 1 to represent its disrupted state. Note that the classification does not apply to clusters (groups) without morphological features outside their tidal radii.

4.1.1. Filamentary and Fractal Types

Among the 60 open clusters in our sample that are younger than 100 Myr, 41 host extended substructures outside their tidal radii. These unidirectional elongated substructures are often referred to as filamentary (Jerabkova et al. 2019) or string-like structures (Kounkel & Covey 2019) in the literature. These types of structures are thought to be the relics of the star formation process that took place along the filaments in their parent molecular clouds (see Section 1). Coeval relic filaments have been found in the Orion cluster (Jerabkova et al. 2019), the Vela OB2 association (Cantat-Gaudin et al. 2019b; Beccari et al. 2020; Pang et al. 2021b; Wang et al. 2022), and other young (<100 Myr) clusters (Kounkel & Covey 2019).

We classify young clusters with a unidirectional filament-like substructure as the filamentary type (f1 type). Other young clusters show multidirectional fractal substructures are categorized as fractal type (f2 type). Examples of these two morphological types are shown in Figure 5 (in blue and red, respectively).

Zoom In Zoom Out Reset image size

According to the hierarchical formation scenario, if the filamentary substructures are formed in a low star formation efficiency (SFE) region, they will disperse after gas expulsion (Kruijssen 2012). However, if the filamentary substructures are formed in a high-SFE region, the filaments may have a chance to converge into the central hub and then merge with other filaments to form a dense cluster (Vázquez-Semadeni et al. 2017; Treviño-Morales et al. 2019; Ward et al. 2020). During this merger process, filaments interact and develop a fractal spatial configuration that might reflect an underlying dynamical interaction between subgroups of stars (Clarke 2010; Arnold et al. 2017; Fujii et al. 2021). As shown in Figure 5, there is an obvious overdensity outside the central core of BH 99 at (X_c , Z_c ) ∼ (100, 5) pc, resembling a subgroup in the merging process (Arnold et al. 2017; Vázquez-Semadeni et al. 2017).

4.1.2. Halo Type

Some clusters older than 100 Myr that still host a dense central core show a halo of low-density stars spread out in the outskirt. We considered this substructure as halo (h type). Only three clusters are classified as h type, and all are shown in Figure 5 (purple-dotted panels). Massive clusters are less susceptible to external disruption forces like the Galactic tide (Li et al. 2021); thus, as these h-type clusters all have a total mass of >600 M_⊙, their current internal kinematics is probably still dominated by the two-body relaxation. This relaxation process tends to result in low-mass stars to migrate toward the outskirts of the cluster. A previous study had found tidal tails ∼200 pc long in both directions for the Praesepe cluster (Röser & Schilbach 2019); however, because of the 100 pc limited size in the spatial cut, we can only recover its closest unbound halo substructures in this study. Praesepe is therefore classified as a halo type.

4.1.3. Tidal-tail Type

Unidirectional elongated substructures are detected in 10 clusters older than 180 Myr. We consider these substructures as tidal tails and designate them as t type. Besides four clusters that previously already been reported with tidal tails or early tidal-tail structure: Blanco 1, Coma Berenices, NGC 2516, and Pleiade (Tang et al. 2019; Zhang et al. 2019; Pang et al. 2021a; Lodieu et al. 2019; Li et al. 2021; Bouma et al. 2021), six more tidal tail clusters are detected in this study: Alessi 3, Collinder 350, IC 4756, LP 2429, Mamajek 4, and NGC 7092. The 3D spatial distribution of these six clusters is shown in Figure 6. The length of the tidal tails in these six clusters ranges from ∼20 to ∼50 pc (bottom panels in Figures 5 and 6). The elongation of the tidal tails is mostly parallel to the Milky Way disk, consistent with Paper I, reflecting that the tidal effects in these clusters have successfully overcome the effects of the internal dynamics.

Zoom In Zoom Out Reset image size

Stars in the tidal tails escape from the cluster through two Lagrangian points (e.g., Küpper et al. 2011). Because of the differential rotation of stars on the Galactic plane, theoretically, stars in the leading tail will orbit faster than the bulk motion of the cluster, while those in the trailing tail will stay behind the cluster as they orbit slower. Such stellar motions in the tails can lead to the observed inclination between tidal tails and the orbital direction of the star cluster (red arrows in the X–Y plane in Figure 6).

4.2. Radial Density Profile

To further quantify the four types of clusters hosting morphological substructures, we obtain the radial density profile (RDP) for all 85 target clusters and fit the Elson, Fall, and Freeman (EFF) model (Elson et al. 1987) to the RDP. The EFF model is characterized by an RDP given by

$$\begin{eqnarray}&&\rho (r)={\rho }_{0}{\left[1+{\left(\displaystyle \frac{r}{a}\right)}^{2}\right]}^{-\gamma /2},\end{eqnarray} \tag{ 1 }$$

where ρ is the number density, and r is the distance of a cluster member to the cluster center. The remaining parameters are fitted; these are the central number density, ρ₀; the power-law slope at large radii, γ; and the scale radius, a. A large value of γ corresponds to a steep radial distribution, and vice versa. The scale radius a is a measurement of the cluster's physical size. We also compute the core radius, r_c , with a and γ obtained from fitting the EFF profile using Equation 22 in Elson et al. (1987):

$$\begin{eqnarray}&&{r}_{c}\approx a\sqrt{{2}^{2/\gamma }-1}.\end{eqnarray} \tag{ 2 }$$

The fitted results of the EFF model for the example clusters shown in Figure 5 are displayed in Figure 7 (gray curves). Most of the clusters have the highest ρ at the center, except two hierarchical groups, Alessi 20 gp1 (f2 type) and UBC 31 gp2 (f1 type) having a value for ρ in the innermost 2 pc that is several times lower than the regions further away from cluster center. These noncentrally concentrated groups are excluded from further analysis, as the EFF profile cannot accurately describe them.

Zoom In Zoom Out Reset image size

The fitting of EFF profile becomes worse at larger radii (>20 pc), where substructures emerge and have a slight increase in ρ (e.g., filamentary, fractal, and halo types shown in Figure 7). The discrepancy becomes worst in the halo cluster Praesepe (purple), which has an extremely compact profile (γ = 7.62). The ρ value of Praesepe seems to remain constant at radius outside 20 pc; however, it is a result of low-number statistics (contributed by only 31 stars) so this does not imply any significant feature. This ρ plateau phenomenon becomes less pronounced in tidal-tail clusters, which are at a more advanced disruption state than halo type. In the tidal-tail clusters, a large number of tidally stripped stars probably have escaped to large distances from the cluster. Because of this significant mass loss, the density profile consequently becomes the shallowest among the tidal-tail clusters (green).

Considering the similarity in cluster ages of types f1 and f2, we combine these two samples to study the morphological substructures of young clusters (groups). In Figure 8, we show the distribution of ρ₀, r_c , and γ for three types of clusters (f1+f2, h, and t). Halo clusters with the highest central number density have the highest mean ρ₀ (panel (a) in Figure 8). The outlier with the highest ρ₀ value (>4 pc⁻³) is the distinct tidal-tail cluster, the Pleiades.

Zoom In Zoom Out Reset image size

The average values of r_c is very similar for the h-type and the t-type clusters, while the young f1+f2 type clusters tend to have larger cores. This result is consistent with the earlier study of Tarricq et al. (2022). As clusters grow older, mass segregation manifests. Lower-mass stars tend to obtain higher speeds over time, and migrate to the peripheral area, while massive objects sink to the center (e.g., Hillenbrand & Hartmann 1998; Pang et al. 2013; Tang et al. 2018). Over a longer timescale, this process will lead to the shrinking of core radius (Heggie & Hut 2003).

The value of γ decreases from the f1+f2 type to the h type and to the t type (Figure 8 (c)). A higher value of γ indicates more compact profile at large radius. A weak tendency of decreasing γ along the increasing cluster age can be seen in panel (d). Young (f1+f2) clusters have a wider range of profiles, both compact and shallow. Older t-type clusters only have small values of γ meaning a flat distribution, which should be attributed to their advanced disruption state. Therefore, the value of γ can be considered as an indicator of a clusters' dynamical state, which is consistent with the substructure classification based on 3D morphology. Our findings agree with previous results from globular clusters, that tidally affected clusters have a flatter profile (i.e., a small value of γ) than the unaffected ones (Carballo-Bello et al. 2012).

The unbound stars have been expelled either by stellar encounters or relaxation processes, and form the nearby halo populations. As clusters grow older, more stars will escape, which follow the Galactic shear and form the tidal tails. As a result, clusters with tidal tails have a more shallow and extended density profile than their younger counterparts. Halo clusters and tidal-tail clusters are undergoing similar dynamical disruption processes. Owing to the higher mass and higher density (panel (a) in Figure 8), halo clusters have a longer relaxation time. Therefore, it takes a longer time for the tidal tails to appear.

Expansion is often considered as a diagnostic of ongoing cluster dissolution. This signal is seen in both young and old clusters in a number of studies (Wright & Mamajek 2018; Cantat-Gaudin et al. 2019a; Pang et al. 2018; Tang et al. 2019; Pang et al. 2021b) using the 3D velocities or the 2D PMs. For old clusters hosting tidal tails, an expansion signal can be interpreted as the presence of the tidal disruption process; however, this is not always true for young clusters (groups). Only young stellar groups showing significant expansion with the speed of members increasing linearly toward the outer region are suggested to be undergoing infant disruption (Wright & Mamajek 2018; Kuhn et al. 2019; Pang et al. 2021b). Local expansion in other young groups is a dynamic relics of molecular gas cloud and will not affect the large-scale structure (Ward et al. 2020). When dynamical conditions are appropriate, young stellar groups interact and merge (Goodwin & Whitworth 2004; Kuhn et al. 2015). To trace the complexity of stellar motion, we use 3D velocities to search for global dispersion or interacting signatures in our target clusters (groups).

5.1. 3D Kinematic Substructures

In Figure 9, we present the 3D relative velocity (after bulk motion subtracted) computed from PMs and RVs for the example clusters shown in Figure 5. The median velocity value of cluster members is taken as the bulk motion of the cluster. In the filamentary (f1-type) group (blue arrows), stars of LP 2373 gp2 are moving perpendicularly away from the ridgeline of the filamentary structure. A similar trend is also detected in ASCC 127, but less significant. This kind of orthogonal expansion seen in young clusters or groups may reveal their destiny of total dispersal. Considering that the uncertainty of Gaia RV might bias the result, we require more accurate kinematic data to confirm the expansion hypothesis. For the fractal (f2-type) clusters (red arrows), simple expansion-like motion is evident in LP 2383, but stars move in a more complex way in NGC 2451B and BH 99. This complexity might be a reflection of turbulence in their parental molecular cloud (Elmegreen et al. 2000) or an aftermath of a merger event (Vázquez-Semadeni et al. 2017). Whether these expansions will impact the young clusters' large-scale structure or not needs to be confirmed with additional high-resolution spectroscopy. With the accuracy of Gaia DR2 RV, we cannot reach a quantitative conclusion.

Zoom In Zoom Out Reset image size

The star motions become complex in the bound region (within the tidal radius, gray dashed circle) of halo (h-type) clusters (purple arrows), with both radial expansion and tangential motions, which indicates the simultaneous interplay between two-body relaxation and Galactic tides. Halo clusters are transiting from early tidal disruption toward a more advanced disruption stage like the tidal-tail clusters. In the tidal-tail (t-type) clusters (green arrows), the outward motions are inclined to parallel to the Galactic Plane because of the stretch from the Galactic shear.

5.2. 2D Kinematic Substructures

Previous observations have shown that morphological substructures often relate to kinematic substructures in the PM space (Wright et al. 2016; Wright & Mamajek 2018; Pang et al. 2021b). We present the PM distributions for the example clusters (those in Figure 5) in Figure 10. To avoid the projection effect, we convert PMs into tangential velocities (V_t ;) and show in Figure 11. Several additional examples of each morphological type are added to explore the difference in the outcome of the kinematic evolution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The stellar group, Group X, is displayed as a reference in both Figures 10 and 11 (orange dots) to showcase the PM and V_t distribution for a disrupted stellar group (Tang et al. 2019). Like Group X, most filamentary (f1-type) clusters have elongated PM and V_t distributions, supporting their unstable and transient dynamical state. This result is consistent with the interpretation of their 3D motions (orthogonal expansion in Section 5.1). As kinematic substructures will be quickly erased by the fast dynamical relaxation or interaction (Parker et al. 2014), the elongated kinematic feature we are currently observing in these f1-type young clusters confirm that they are dynamically unevolved. The unbound kinematic substructures are likely relics of the primordial gas kinematics of their unbound parent molecular cloud (Heyer et al. 2001). In fractal (f2-type) clusters, the PM and V_t distributions have a higher degree of substructure but less elongation. This difference may be caused by two possibilities: (1) the primordial gas kinematics in f2-type clusters is different from f1-type clusters, owing to different SFE (Kruijssen 2012); (2) dynamical processing has taken place in these clusters so that the inherent kinematics of the molecular gas cloud is affected.

The 2D kinematic distribution of halo (h-type) clusters is more spherical with fine substructures (purple dots in Figures 10 and 11). Tail-like substructures appear in the Praesepe, referred to as kinematic tails (Li et al. 2021). The kinematic tails become more significant in tidal-tail (t-type) clusters, e.g., the Pleiades, Blanco 1, and NGC 7092. Most stars in the kinematic tails of t-type clusters are located along the tidal tails (Figure 7). Although the Pleiades only shows an early tidal-tail substructure outside the tidal radius in spatial space (Li et al. 2021) that is much shorter than Blanco 1's tidal tail, the remarkably long kinematic tails in both the PM and the V_t distributions may point out the uniqueness of this cluster, e.g., extremely high central density (blue shaded bin in Figure 8 (a)). The kinematic substructures in Coma Berenices and IC 4756, on the other hand, are less intriguing owing to a small number of tidal-tail members identified in Paper I.

The most complex morphological and kinematic substructures are seen in the hierarchical clustering groups around four clusters, Alessi 20, IC 348, LP 2373, and LP 2442. Together with Vela OB2, another sky region found to have hierarchical groups (Pang et al. 2021b), we highlight them in the large-scale spatial location of solar neighborhood in Figure 4, and provide a close-up look for the four identified in this study in Figure 12. Within these four sky regions, we in total identified 29 groups. Ten of the 29 (about one-third) groups are newly discovered, which we name after the target region.

Zoom In Zoom Out Reset image size

The hierarchical group relating to Alessi 20, is named Alessi 20 gp 1. A new isolated group found around Alessi 20 is named Alessi 20 isl 1. The cluster LP 2442, first discovered by Liu & Pang (2019), also named as UPK 640 in Cantat-Gaudin et al. (2020) host five hierarchical groups: LP 2442 gp 1–5. These nine groups mentioned above are all newly detected and none are cataloged in Liu & Pang (2019) or Cantat-Gaudin et al. (2020).

Six groups are identified in the region of LP 2373 (Liu & Pang 2019): LP 2373 gp 1–4, ASCC 16, and ASCC 19. The latter two are suggested to be a cluster pair in Soubiran et al. (2018). LP 2373 gp 1 is a newly discovered group that is not shown in Liu & Pang (2019) and Cantat-Gaudin et al. (2020). LP 2372 gp 2 is the biggest group in the LP 2372 region with 955 members. With this larger population, LP 2372 gp 2 have a fraction of members overlap with that in LP 2369 (240, Liu & Pang 2019), LP 2367 (117, Liu & Pang 2019), and ASCC 21 (82, Cantat-Gaudin et al. 2020). As for LP 2372 gp 3 (185 members), there are 66 members matches with LP 2370 (Liu & Pang 2019). LP 2372 gp 4 being the second massive subgroup (511 member), one-fifth of its members are matched with UBC 17a in (183, Cantat-Gaudin et al. 2020).

The majority of the hierarchical groups are younger than 30 Myr. They are either filamentary (f1) or fractal (f2) type (Section 4.1). We plot their 3D spatial distribution, 3D velocity, and PM distribution in Figure 12. The morphology of spatial substructures in these four regions varies significantly (Figure 12 panels (a) and (c)). Consistent with morphological substructures, a higher degree of kinematic substructures (Figure 12 panels (d)) are found in these hierarchical groups. We discuss the unique features of each hierarchical clustering region in the following sections.

6.1. Dispersing Filament Network: Alessi 20 and LP 2373 Regions

6.2. Infalling Motions: IC 348 Region

6.3. Well-mixed Structure: LP 2442 Region

As demonstrated in Section 2.3, cluster members identified within the tidal radius are most stringent and least affected by the contamination rate. In this section, we analyze the shape of all our target clusters within the bound region following the method developed in Paper I by performing the ellipsoid fitting. ¹² In Paper I, we found that, in general, the spatial distribution of member stars within the tidal radius of open clusters can be well described by oblate spheroids, prolate spheroids, or triaxial ellipsoids. During the ellipsoid fitting, we first center the fitted ellipsoid at the median position of all bound members, where we considered the cluster center. For the fitted semi-axes of the ellipsoid, a, b, and c, we make a the semimajor axis, b the semi-intermediate axis, and c the semiminor axis. We use the lengths of the semi-axes a, b, c, and axis ratios b/a and c/a to describe the morphology of the clusters. The direction of the semimajor axis a is defined as the direction of elongation of the fitted cluster. Smaller values of both the axis ratios, b/a and c/a, indicate a more elongated structure. The fitted values of the morphological parameters (a, b, c, b/a and c/a) are listed in Table 3. We also compute for each cluster the angle θ between the direction of a and the Galactic plane, and the angle ϕ between the direction of a and the line of sight. The values of these two angles are also given in Table 3. The morphological parameters of the 12 clusters from Paper I are taken from Table 3 in Paper I.

As can be seen in Figure 13, there is no dependence of the ratio b/a on ϕ (the angle between the elongated direction and the line-of-sight direction) with the Spearman coefficient: s = 0.09, and c/a on ϕ with s = 0.21. The slightly larger value of s = 0.21 in c/a versus ϕ does not truly imply any correlation. This evidence reassures that the distance correction is appropriate for members within the tidal radius. Surprisingly, the c/a ratio is also independent of θ (s = 0.13), the angle between the direction of elongation and the Galactic plane. As the long-term effect of the Galactic tides on clusters is to make their members show elongation parallel to the Galactic plane, this lie down effect is shown as the smaller θ goes, the smaller the c/a ratio. This trend is not seen in our cluster samples, we find that filamentary structures (small c/a) in young clusters tend to have a random angle relative to the Galactic plane (f1 type in Figure 5, Figure 14 (b)), which successfully dilutes the anticipated effect produced by the tidal field.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Based on the axis ratio, we define three kinds of shape describing the spatial distributions of stars inside tidal radius, based on the following criteria:

• 1.  
  Type I. oblate spheroid:

  $$\begin{eqnarray}&&\displaystyle \frac{c}{a}\lt \displaystyle \frac{b}{a},\quad \displaystyle \frac{b}{a}\geqslant 0.8,\quad \displaystyle \frac{c}{a}\lt 0.6.\end{eqnarray} \tag{ 3 }$$
• 2.  
  Type II. prolate spheroid:

  $$\begin{eqnarray}&&\left|\displaystyle \frac{c}{a}-\displaystyle \frac{b}{a}\right|\lt 0.1,\quad \displaystyle \frac{b}{a}\lt 0.6,\quad \displaystyle \frac{c}{a}\lt 0.6.\end{eqnarray} \tag{ 4 }$$
• 3.  
  Type III. triaxial ellipsoid:

  $$\begin{eqnarray}&&\left|\displaystyle \frac{c}{a}-\displaystyle \frac{b}{a}\right|\geqslant 0.2.\end{eqnarray} \tag{ 5 }$$

The shapes of stellar populations inside the tidal radius give 14 clusters that resemble an oblate ellipsoid, three as a prolate ellipsoid, and seven as a triaxial ellipsoid. Pleiades, the cluster with a prominent kinematics tail (Figures 10 and 11), exhibits a oblate shape. The most massive cluster in our sample, NGC 3532, resembles a triaxial ellipsoid. The tidal-tail cluster, Collinder 350 is oblate, with its a axis perfectly align with the Galactic plane. In all the hierarchical clustering groups, the LP 2442 region hosts the largest variety of morphologies. Three groups are oblate (LP 2442, LP 2442 gp 1, LP 2442 gp 4), two triaxial (LP 2442 gp 2, LP 2442 gp 5), and one prolate (LP 2442 gp 3).

The diversity of cluster morphology inside the tidal radius reveals different dynamical state among clusters. The discussion of individual cluster dynamical state requires a higher precision kinematic data and therefore is beyond the scope of this paper. However, the parameterization of 3D positional distribution of clusters do allow us to investigate the dependence of cluster morphology (e.g., axis ratio c/a) on other structural parameters (e.g., r_h) and properties (e.g., age), paving the way to understand how the internal and external dynamical processes reshape the cluster's spatial distribution over time.

The morphology of an open cluster evolves as it ages. It changes from inherent filamentary or fractal substructures to acquired halo or tidal tails. The variation of spatial distribution reflects different dynamical processes dominating the current cluster evolution. Inside the tidal radius, the influence of the cluster's gravitational potential is competing with Galactic tides. It is intriguing to investigate the similarity and differences in the morphology of young clusters' bound regions to explore the early dynamical evolution of clusters.

The relationship between clusters' spatial morphology and other properties can be linked to their evolution in space and time. We present the correlation between cluster age and other parameters in Figure 14. The mean ratio c/a seems to increase as the cluster grows older until 30 Myr (panel (a)), though a large scatter is observed. Most clusters younger than ∼30 Myr are filamentary or fractal. These young clusters are elongated in a direction with a random inclination to the Galactic plane (panel (b)). The mean value of θ tends to decrease from ∼30° at ∼10 Myr to ∼10°–15° in ∼30 Myr (panel (b)).

Many filamentary young clusters (groups) with small c/a are dispersing (Section 6.1) so that they were not able to grow older than 30 Myr. At the same time, other groups in the hierarchical clustering region follow the infalling flow and merger to form a denser cluster (Sections 6.2 and 6.3). During this merging event, the primordial filamentary morphology is significantly modified, and produces a more spherical morphology by the dynamical mixing. The inclination of the elongation direction of the cluster to the Milky Way disk, θ, as a result, is significantly reduced.

The apparent increasing trend of the mean value of c/a, M_cl, and r_h for the young clusters (groups) can be explained by filamentary group dissolution and group mergers in the hierarchical formation framework. Though the age correlation has a large scatter, this overall picture derived from it reaches an agreement with our morphological and kinematic analysis of clusters in previous sections (Sections 4–6).

The mean value of c/a declines after ∼30 Myr, and stays at ∼0.5 for age above ∼100 Myr (panel (a)). The inclination angle to the Galactic plane, θ, is generally smaller in old clusters (panel (b)). As Galactic tide becomes increasingly dominant in older clusters, clusters are elongated again (smaller c/a) with the development of tidal tails. This disruption process should go along with mass loss; however, we do not observe a significant decline of cluster mass in old clusters (panel (c)). This inconsistency shows that only the most massive and dense clusters can survive long enough under the influence of violent Galactic force. The decrease of r_h in older clusters is attributed to the two-body relaxation that generates more compact clusters (Heggie & Hut 2003), consistent with the finding of r_c in Section 4.2.

Observations have found that clusters younger than ∼30 Myr are mostly moving around the Galactic center in circular orbits (Tarricq et al. 2022). They follow much lower vertical velocities than the old clusters (Soubiran et al. 2018; Kounkel 2020). It seems that they are moving slowly in the dense star formation region in the spiral arms (Figure 4) and facing disruption and merger at the same time. On the other hand, old clusters run faster along the Z direction and manage to locate further away from the disk, which alleviates their disruption rate.

We obtain membership for 65 open clusters and stellar groups in the solar neighborhood via StarGO based on Gaia EDR3 data. Ten groups of stars in four hierarchical clustering regions (hosting 29 groups) are newly discovered. Those newly discovered groups in different sky regions are (1) Alessi 20 gp 1 and Alessi 20 isl 1 in the Alessi 20 region; (2) LP 2373 gp 1 in the LP 2373 region; (3) UBC 31 gp 1 and UBC 31 gp 2 in the IC 348 region; and (4) LP 2442 gp 1–5 in the LP 2442 region.

We verify the robustness of cluster membership by calculating the recovery rate of members from the simulated observation data set (combination of N-body simulated clusters and Gaia EDR3 mock field stars). The overall morphology of clusters can be strongly affected by a contamination rate >20%. An artificial halo-like substructure will appear when the contamination rate rises above 20%. The bound members (stars within tidal radius); however, unlike those unbound, are more stable in their membership. The recovery rate of bound members remains above 99% despite the overall contamination changing from 10% to 25%.

Besides the 65 clusters and groups, we adopt 12 clusters from Paper I, seven open clusters and groups in Vela OB2 from Pang et al. (2021b) and the Pleiades from (Li et al. 2021). The morphological parameters of ellipsoidal fitting a total of 72 open clusters (groups) are being computed for the first time by this work. Based on the morphological parameters derived from ellipsoidal fitting, the 3D morphology of the bound region of 14 clusters resemble oblate spheroid, three prolate spheroid, and seven triaxial ellipsoids.

We adopt the morphological parameters of 12 clusters from Paper I. In total, 85 clusters are analyzed for the 3D morphology investigation. We classify the morphological substructures detected in these 85 clusters into four types: (1) the filamentary (f1) type, for clusters younger than 100 Myr and having unidirectional elongated substructures; (2) the fractal (f2) type, for clusters that are as young as the filamentary type but with multidirectional and fluffy spatial substructures; (3) the halo (h) type, for clusters older than 100 Myr and harboring a diffuse halo around the massive core; and (4) the tidal-tail (t) type, for old clusters (>100 Myr) with extended tidal tails.

The RDP of 82 clusters is fitted with the EFF profile (Elson et al. 1987). Group X, Alessi 20 gp1 and UBC 31 gp2 are excluded in the fitting owing to irregular shapes. As a result of tidal disruption, t-type clusters have the smallest size of cluster core and have the most shallow RDP. Although clusters in the h type are as old as those in the t type, their high-mass and high-density nature slow down their disruption process.

Kinematic substructures are identified among all four types of morphological substructures. Clusters of f1 type have elongated kinematic substructures in the PM space, similar to the disrupted cluster Group X. This elongated feature seen in the velocity space for the f1-type clusters shows their unstable and short-lived dynamical state. Kinematic tail extended from the dense core in the PM diagram is most prominent in the t-type clusters, especially the Pleiades.

Morphology and kinematics of hierarchical clustering groups are highly complex. They form extended filament networks like in the Alessi 20, IC 348, and LP 2373 regions. Global orthogonal expansion seems to disperse the filaments in Alessi 20 and LP 2373. Tangential infalling motion along filament direction is detected in UBC 31, UBC 31 gp 1, and UBC 31 gp 2 in the IC 348 region. Hierarchical subgroups in the LP 2442 region are well mixed in positions and are distinctive in PM space. The complexity of the 3D motions in LP 2442 gp 1–5 might be attributed to an ongoing merger event.

Correlations between cluster age and the three structural parameters, the axis ratio c/a, the inclination of cluster elongated direction a to the Galactic plane θ, cluster mass M_cl, and half-mass–radius r_h, are observed. The mean value of c/a, M_cl, and r_h all follow an increasing trend with larger cluster ages when below ∼30 Myr. The mean θ, on the contrary, shows decreasing trend before ∼30 Myr. These signs imply two possible dynamical processes that young clusters (groups) may experience before 30 Myr. (1) Young stellar groups with low SFE (mostly along filaments in the molecular clouds) quickly dissolve and become unbound. (2) Filaments at high SFE, at the location of the hub or major filaments in the molecular cloud, have a chance to merge and form a dense cluster.

We thank the anonymous referee for advice that improved the paper. We are also grateful for the IT support of Lowell Observatory, and thank Xianhao Ye and Chaojie Hao for helpful discussions. This work is supported by grant No. 12173029 from the National Natural Science Foundation of China. X.Y.P. is grateful for the financial support of the Research Development Fund of Xi'an Jiaotong-Liverpool University (RDF-18–02–32). X.Y.P. gives thanks for the funding of Natural Science Foundation of Jiangsu Province, No. BK20200252. L.W. thanks the support from the one-hundred-talent project of Sun Yat-sen University and the National Natural Science Foundation of China through grant 12073090. M.B.N.K. is grateful for the support of the Continuous Support Fund (grant RDF-SP-93) of Xi'an Jiaotong-Liverpool University.

This work 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). This study also made use of the SIMBAD database and the VizieR catalog access tool, both operated at CDS, Strasbourg, France.

Software: Astropy (Astropy Collaboration et al. 2013, 2018), SciPy (Millman & Aivazis 2011), NumPy (Harris et al. 2020), galpy (Bovy 2015), TOPCAT (Taylor 2005), and StarGO (Yuan et al. 2018), Plotly (Plotly Technologies Inc 2015), petar (Wang et al. 2020a).