<figure><img src="image_1.png" /></figure>

HAL
open science

# Disruption of Hierarchical Clustering in the Vela OB2 Complex and the Cluster Pair Collinder 135 and UBC 7 with Gaia EDR3: Evidence of Supernova Quenching

Xiaoying Pang, Zeqiu Yu, Shih-Yun Tang, Jongsuk Hong, Zhen Yuan, Mario

Pasquato, M. Kouwenhoven

► To cite this version:

Xiaoying Pang, Zeqiu Yu, Shih-Yun Tang, Jongsuk Hong, Zhen Yuan, et al. Disruption of Hierarchical Clustering in the Vela OB2 Complex and the Cluster Pair Collinder 135 and UBC 7 with Gaia EDR3: Evidence of Supernova Quenching. The Astrophysical Journal, 2021, 923 (1), pp.20. 10.3847/1538-4357/ac2838. hal-03479259

HAL Id: hal-03479259

https://hal.science/hal-03479259v1

Submitted on 2 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Disruption of Hierarchical Clustering in the Vela OB2 Complex and the Cluster Pair Collinder 135 and UBC 7 with Gaia EDR3: Evidence of Supernova Quenching

Xiaoying Pang<fnref n="1" />,<fnref n="2" /> ID, Zeqiu Yu<fnref n="1" /> ID, Shih-Yun Tang³,<fnref n="4" /> ID, Jongsuk Hong⁵ ID, Zhen Yuan<fnref n="6" /> ID, Mario Pasquato⁷,<fnref n="8" /> ID, and M. B. N. Kouwenhoven<fnref n="1" /> ID

## Abstract

We identify hierarchical structures in the Vela OB2 complex and the cluster pair Collinder 135 and UBC 7 with Gaia EDR3 using the neural network machine-learning algorithm StarGO. Five second-level substructures are disentangled in Vela OB2, which are referred to as Huluwa 1 (Gamma Velorum), Huluwa 2, Huluwa 3, Huluwa 4, and Huluwa 5. For the first time, Collinder 135 and UBC 7 are simultaneously identified as constituent clusters of the pair with minimal manual intervention. We propose an alternative scenario in which Huluwa 1–5 have originated from sequential star formation. The older clusters Huluwa 1–3, with an age of 10–22 Myr, generated stellar feedback to cause turbulence that fostered the formation of the younger-generation Huluwa 4–5 (7–20 Myr). A supernova explosion located inside the Vela IRAS shell quenched star formation in Huluwa 4–5 and rapidly expelled the remaining gas from the clusters. This resulted in global mass stratification across the shell, which is confirmed by the regression discontinuity method. The stellar mass in the lower rim of the shell is $0.32 \pm 0.14 M_\odot$, higher than in the upper rim. Local, cluster-scale mass segregation is observed in the lowest-mass cluster Huluwa 5. Huluwa 1–5 (in Vela OB2) are experiencing significant expansion, while the cluster pair suffers from moderate expansion. The velocity dispersions suggest that all five groups (including Huluwa 1A and Huluwa 1B) in Vela OB2 and the cluster pair are supervirial and are undergoing disruption, and also that Huluwa 1A and Huluwa 1B may be a coeval young cluster pair. N-body simulations predict that Huluwa 1–5 in Vela OB2 and the cluster pair will continue to expand in the future 100 Myr and eventually dissolve.

*Unified Astronomy Thesaurus concepts: Star clusters (1567); Open star clusters (1160); Stellar kinematics (1608)*

*Supporting material: machine-readable table*

## 1. Introduction

Star clusters tend to form in complexes (Efremov 1978). A typical complex contains several OB associations and a few young clusters, with a star formation process lasting from a few to a dozen million years. Stellar complexes tend to be hierarchical, with smaller substructures (of roughly several to several tens of parsecs in size) embedded in larger structures (a few hundred parsecs in size). This feature is inherited from the parental molecular clouds. This is often accompanied by a hierarchical configuration shaped by supersonic turbulence (Elmegreen et al. 2000; Orkisz et al. 2017; Torniamenti et al. 2021). Many studies have identified hierarchies in stellar structures in the solar neighborhood (e.g., Megeath et al. 2016; Kounkel et al. 2018; Ballone et al. 2020; Kerr et al. 2021).

Two strategies are commonly used to identify different levels of hierarchical structures. The first strategy is the “bottom-up” approach, which initially identifies the lowest-level structures and merges these to form higher-level larger structures (Cantat-Gaudin et al. 2019a). The second strategy is the “top-down” method, which starts by initially identifying the highest-level structures and then moves down by fragmentation to locate lower-level substructures (Beccari et al. 2018; Cantat-Gaudin et al. 2019b; Kerr et al. 2021).

The Vela OB2 complex is a young stellar grouping located in the direction of the Vela and Puppis constellations, at a distance of 350–400 pc from the Sun (Pozzo et al. 2000). De Zeeuw et al. (1999) first made use of proper motions (PMs) and parallaxes from Hipparcos to identify members of Vela OB2. Subsequently, Pozzo et al. (2000) used X-ray observations and identified a group of pre-main-sequence (PMS) stars in this region, which is often referred to as Gamma Velorum (hereafter Gamma Vel), or Pozzo 1 (Dias et al. 2002). The formation of the Vela OB2 structure is thought to have been triggered by a supernova explosion of a $15 M_\odot$ star (Cantat-Gaudin et al. 2019b) that was located at the center of the IRAS Vela shell (Sahu 1992). A similar supernova scenario has also been suggested for the origin of the RCW 34 complex (Bik et al. 2010) and the Orion complex (Kounkel 2020).

Recent studies have identified the hierarchical clustering in Vela OB2, mostly using “top-down” approaches. Beccari et al. (2018) adopted the DBSCAN method to identify six clusters in Vela OB2 covering a distance range of 82 pc, including the Gamma Vel cluster and NGC 2547. Cantat-Gaudin et al. (2019b) used the UPMASK algorithm (Cantat-Gaudin et al. 2018) and increased the number of identified clusters in Vela OB2 to 11, which extends ~133 pc in space. Groups in Vela

<footnote n="1">¹ Department of Physics, Xi’an Jiaotong-Liverpool University, 111 Ren’ai Road, Dushu Lake Science and Education Innovation District, Suzhou 215123, Jiangsu Province, People’s Republic of China; Xiaoying.Pang@xjtu.edu.cn</footnote>

<footnote n="2">² Shanghai Key Laboratory for Astrophysics, Shanghai Normal University, 100 Guilin Road, Shanghai 200234, People’s Republic of China
³ Lowell Observatory, 1400 W. Mars Hill Road, Flagstaff, AZ 86001, USA</footnote>

<footnote n="4">⁴ Department of Astronomy and Planetary Sciences, Northern Arizona University, Flagstaff, AZ 86011, USA
⁵ Korea Astronomy and Space Science Institute, Daejeon 34055, Republic of Korea</footnote>

<footnote n="6">⁶ Université de Strasbourg, CNRS, Observatoire Astronomique de Strasbourg, UMR 7550, F-67000 Strasbourg, France
⁷ Center for Astro, Particle and Planetary Physics (CAP³), New York University Abu Dhabi, United Arab Emirates</footnote>

<footnote n="8">⁸ INFN- Sezione di Padova, Via Marzolo 8, I-35131 Padova, Italy

*Received 2021 June 11; revised 2021 August 20; accepted 2021 September 4; published 2021 December 8*</footnote>
OB2 are likely formed along the densest filaments of the same giant molecular clouds.

Survival of hierarchical clusters in such complexes mainly depends on the star formation efficiency (SFE) and the rate of gas expulsion. The SFE of clusters ranges from a few percent to 30% or above, depending on the cluster surface density of the region in which the clusters are born (Evans et al. 2009; Megeath et al. 2016). The remaining gas in the cluster is removed through stellar feedback. Violent feedback, such as the feedback caused by supernova explosions, removes the gas much more quickly than stellar winds or radiation. The faster the gas is expelled, the smaller the chance of cluster survival (Dinnbier & Kroupa 2020a, 2020b). Among clusters with an identical gas removal timescale, those with an SFE smaller than 33% will disrupt faster (Baumgardt & Kroupa 2007).

Besides disruption, another essential physical process in hierarchical clusters in the complex is the mutual gravitational interaction between its substructures. De la Fuente Marcos & de la Fuente Marcos (2009a) suggested that such interactions can be important when the physical separation between neighboring clusters is less than 30 pc, which is roughly three times the typical value of the tidal radius for open clusters in the Galactic disk (e.g., Binney & Tremaine 2008). When the relative velocity between two neighboring clusters is smaller than their internal dispersions (Gavagnin et al. 2016), clusters are strongly affected by the dynamical interaction with their neighbors. This dynamical interaction may foster the formation of gravitationally bound cluster pairs (e.g., Priyatikanto et al. 2016; Arnold et al. 2017).

Many cluster pairs have been observed in both the Milky Way (de la Fuente Marcos & de la Fuente Marcos 2009a, 2009b; Currie et al. 2010) and the Large Magellanic Cloud (Dieball et al. 2002; Dieball & Grebel 1998; Palma et al. 2016). About half of the Galactic cluster pairs in the work of de la Fuente Marcos & de la Fuente Marcos (2009a, 2009b) are younger than 25 Myr, and almost all these pairs are coeval. Since 3D kinematic data were unavailable in earlier studies, determining the physical status of these cluster pairs requires further investigation.

Located a projected distance of 15° from Vela OB2 is the 40 Myr old coeval cluster pair of Collinder 135 and UBC 7, which were first discovered by Castro-Ginard et al. (2018) and were later studied in further detail by Kovaleva et al. (2020). The distance to this cluster pair is ~300 pc (Kovaleva et al. 2020). In this work we will investigate the 3D spatial configuration and the kinematic relation between these two constituent clusters.

The Gaia Early Data Release 3 (EDR 3; Gaia Collaboration et al. 2021) provides a 30% higher precision in parallaxes and doubled the accuracy in PMs, when compared to Data Release 2 (DR 2; Gaia Collaboration et al. 2018). The 3D spatial positions and PMs allow for an accurate determination of cluster membership and enable measuring the physical separation between neighboring clusters. The hierarchical structure of Vela OB2 and the neighboring cluster pair provide an excellent laboratory to study the formation, dynamical evolution, and fate of hierarchical clusters that have formed from the same molecular cloud. Here, we present a detailed study of the Vela OB2 complex and the cluster pair Collinder 135 and UBC 7, using data from Gaia EDR3 and the Gaia-ESO survey (Gilmore et al. 2012). We aim to quantify the dynamical states of individual clusters and identify interactions between the clusters.

This paper is organized as follows. In Section 2.1 we discuss the quality and limitations of the Gaia EDR 3 data and describe

our input data set for member star identification. We then present the algorithm, StarGO, which is used to disentangle the hierarchical clustering in Section 2.2. The stellar membership of the hierarchical structures of Vela OB2 and the cluster pair is presented in Section 3, in which the age spread among different clusters in Vela OB2 is discussed. We discuss the kinematic characteristics of each structure in Section 4. The 3D spatial configuration and mass distribution of Vela OB2 and the cluster pair are presented in Section 5. In this section, we use the regression discontinuity method (Section 5.2) and the minimum spanning tree method (Section 5.3) to quantify the mass distribution. The dynamical state of the clusters is investigated in Section 6. We present the star formation history in the target regions in Section 7. In Section 8 we conduct N-body simulation to predict the future evolution of Vela OB2 and the cluster pair. Finally, we provide a brief summary of our findings in Section 9.

## 2. Identification of Hierarchical Clustering

### 2.1. Gaia EDR3 Data Reduction

The analysis described below is carried out for the Vela OB2 complex and for the cluster pair UBC 7 and Collinder 135. The spatial and kinematic structure of the target regions is investigated using Gaia EDR 3 data, within 100 pc of the position of Vela OB2 ($X = -53.9$ pc, $Y = -334.6$ pc, $Z = -33.4$ pc, in heliocentric Cartesian coordinates with the Sun located at the origin) and of the mean position between Collinder 135 and UBC 7 ($X = -103.7$, $Y = -266.4$ pc, $Z = -61.3$ pc), taken from Liu & Pang (2019). We remove artifacts and poor-quality sources from the sample by filtering out stars with relative uncertainties in the parallax or photometry above 10%, and with additional astrometric cuts described in Lindegren et al. (2018, their Appendix C). After applying this procedure, 327,950 stars remain for subsequent reduction.

To allow the algorithm to identify fine structures in the target region more efficiently, we construct a 2D density map of PMs to select stars around the overdensities of Vela OB2 and the cluster pair. Figures 1(a).1 and (a).2 show the 2D density maps of the PMs for the target regions. These maps only show bins with overdensities >3σ of the number of stars in all bins. There is a clear overdensity near the mean PM (Liu & Pang 2019) of the Gamma Vel cluster. The overdensity in the region of Collinder 135 and UBC 7 (indicated with a cross) matches the mean PMs of these two clusters measured by Liu & Pang (2019). We apply a circular cut (the black circles in Figures 1(a).1 and (a).2) to include the target clusters for further analysis. The radius of the circle is chosen to include as many potential members as possible. Although other overdensities of nearby clusters can be observed in Figures 1 (a).1 and (a).2, we only focus on the target clusters and ignore their neighbors in this work. Application of this circular cut reduces the number of candidate member stars in the two target regions to a total of 31,040. Stars in this sample have magnitudes ranging between $G \sim 3.8$ mag and $G \sim 20.3$ mag. All samples are complete for $G \lesssim 19-19.5$ mag.

In this study we use the 5D parameters of stars after PM selection (R.A., decl., $\varpi$, $\mu_o \cos \delta$, and $\mu_\delta$) from Gaia EDR 3. Radial velocities (RVs) from Gaia DR 2 are available for only a small fraction of the stars and have an uncertainty of ~2 km s⁻¹. We adopt the higher-accuracy RVs from Jeffries et al. (2014) as supplementary data for the stars in Vela OB2, to
<figure><img src="image_4.png" /><figcaption>**Figure 1.** (a) 2D density map of the PM vectors for the regions in Vela OB2 (top panels) and the cluster pair Collinder 135 and UBC 7 (bottom panels). The blue crosses indicate the mean overdensity of the Gamma Vel cluster and the mean between Collinder 135 and UBC 7 taken from Liu & Pang (2019). The location of NGC 2547 is also indicated, which is excluded from the circular PM cut. Each bin is smoothed by averaging neighboring eight bins. Only bins with a number count >3σ are shown, where σ is the standard deviation of all bins. The shade of gray indicates the number count in each bin. (b) Histogram of the distribution of *u*-matrix elements. The vertical dashed and solid lines denote the threshold values of *u* that respectively give a 5% and a 1% contamination rate for the identified groups. These correspond to the transparent colored patches and filled colored patches, respectively, in the 2D neural network (panel (c)). (c) 2D neural network resulting from SOM. The neurons with a *u*-threshold corresponding to 5% contamination rates are shown with transparent colors, which as a whole is considered as the top-level structure. The solid color patches having 1% contamination rate are the cores of second-level hierarchical substructures. Vela OB2 (the transparent patch) is divided into five second-level groups, indicated with red, blue, purple, cyan, and green. The cluster pair region (transparent patch) is fragmented into Collinder 135 (yellow) and UBC 7 (pink).</figcaption></figure>

aid the study of the dynamical state of the clusters. The RVs obtained from Jeffries et al. (2014) are part of the Gaia-ESO Survey (GES; Gilmore et al. 2012) and have an uncertainty of 0.4 km s⁻¹.

## 2.2. Disentangling Hierarchical Structures

StarGO<fnref n="9" /> (Yuan et al. 2018) has proven to be successful in membership determination of open clusters (Tang et al. 2019; Pang et al. 2020, 2021) and more diffuse stellar streams (e.g., Yuan et al. 2020a, 2020b). The method is based on the Self-Organizing-Map (SOM), which belongs to the domain of unsupervised learning. SOM can map high-dimensional data onto a 2D neural network, while preserving the topological structures of the input data.

We apply StarGO to the selected sample of stars after applying the circular cut in PM space (see Figures 1(a).1 and (a).2). The input space to the neural network is of the form (X, Y, Z, $\mu_\alpha \cos\delta$, $\mu_\delta$). Stars in clusters form overdensities in the input space and are associated with neuron groups that can be identified by StarGO. We adopt networks of 120 × 120 and 150 × 150 neurons (depending on the remaining number of stars after PM selection is applied) for Vela OB2 and the cluster pair, respectively (see Figures 1(c).1 and (c).2). Each neuron is assigned an initially randomized 5D weight vector, which has

the same dimensions as the input space. During the learning process, the weight vector of each neuron is updated to become closer to the input vector of a given star. One iteration of the training cycle is finished after the neurons are trained by every star from the entire sample. This process is iterated 400 times when the weight vectors converge.

After the self-supervised learning is completed, the difference in weight vectors between adjacent neurons is defined as the *u*-matrix. The distribution of element values in the *u*-matrix is shown in Figures 1 (b).1 and (b).2. The tail on the left side of the distribution consists of neurons with relatively low *u* values, which corresponds to cluster members grouped in the input space. The method for identification is to first associate each star from the entire sample with the neuron that has the closest weight vector to its input vector. Subsequently, we select neurons with *u* below a threshold value indicated by the dashed line in Figures 1(b).1 and (b).2, which forms the transparent colored patches in panels (c).1 and (c).2. This threshold value of *u* is selected to ensure an estimated contamination rate of ~5% among the identified members. Specifically, the contamination rate is estimated from the smooth Galactic disk population using the Gaia EDR3 mock catalog (Rybizki et al. 2020). An identical PM selection to that described in Section 2.1 is applied to the stars in the mock catalog in the same volume of the sky. Each of these mock stars is attached to the trained 2D neural network. We then count the mock stars associated with selected patches. Assuming that the Milky Way population

<footnote n="9">⁹ https://github.com/zyuan-astro/StarGO-OC</footnote>
has the same properties as the mock catalog, the contamination rate of each identified group can be estimated (see details in Pang et al. 2020, 2021). We adopt this *u*-threshold for the ~5% contamination rate and obtain the neuron groups shown with transparent colors in Figures 1(c).1 and (c).2.

We obtain one large transparent patch in each of the two runs in the target regions, which represents the top-level structure in Vela OB2 (panel (c).1) and the cluster pair (panel (c).2). To reveal the inner hierarchical substructures, we further decrease the threshold value of *u* until the contamination rate reaches 1%. The corresponding *u*-values are denoted by the solid lines in Figures 1(b).1 and (b).2, and the resulting patches are shown with solid colors in Figures 1(c).1 and (c).2. These patches represent the core regions of each cluster. We then associate neurons in the transparent outskirts with these cores by calculating the minimum difference in weight vectors between a given neuron and each core patch. For instance, neurons that have the smallest differences to the red core patch compared to the other four core patches are colored red. In this way, all the neurons in the transparent outskirt are assigned to a core patch, denoted by different transparent colors in panels (c).1 and (c).2. In total, Vela OB2 contains five second-level clusters (five transparent color patches in Figure 5(c).1), and the cluster pair is fragmented into two second-level structures that correspond to Collinder 135 and UBC 7 (the two transparent color patches in Figure 5(c).2). This is the first time that Collinder 135 and UBC 7 are simultaneously identified with minimal manual intervention. Previous studies (Castro-Ginard et al. 2018) did not detect Collinder 135 using the DBSCAN method since the structure is spatially very extended.

Our identification of hierarchical clustering is purely based on 5D phase-space information without pre-assumption about age, cluster size, and density threshold. Therefore, the extended distribution of both young and old generations can be identified independently, which enables us to investigate the progression of star formation inside the stellar complex. Top-level clustering is first extracted at a large scale, which includes significant substructures. We further fragment top-level structures into subclusters. This top-down method ensures that the identified hierarchical structures are truly coherent in both space and kinematics.

# 3. Stellar Membership in Hierarchical Structures

## 3.1. Sample Cleaning

We construct color–magnitude diagrams (CMDs) for the member candidates of the five second-level clusters in Vela OB2 (Figure 2) and for Collinder 135 and UBC 7 (Figure 4). Stars in the top-level structure of Vela OB2, together (gray dots in Figure 2(a)), track a locus of a main sequence (MS). The previous estimated age of Vela OB2 ranges from 10 Myr via isochrone fitting (Prisinzano et al. 2016; Cantat-Gaudin et al. 2019b) to 18–21 Myr via lithium abundance (Jeffries et al. 2017). We plot PARSEC isochrones of 10 and 20 Myr (black solid and dashed curves) with solar metallicity and $A_V = 0.35$ (Cantat-Gaudin et al. 2019b). A large number of stars are still in the PMS stage, while fewer are located in the MS. A significant number of upper PMS stars are older than 10 Myr, i.e., bluer than the 10 Myr old isochrone with $M_G$ in the 1.0–5.0 mag range (region I) in Figure 2(a), which were excluded from Cantat-Gaudin et al. (2019b), as they excluded the older population from NGC 2547. NGC 2547 has been

excluded from our analysis using the circular cut of PMs (see Figure 1(a).1). The older upper PMS stars in our sample therefore belong to the Vela OB2 complex. The blue edge of the upper PMS (region I) follows the 20 Myr isochrone. We overplot a binary sequence of the 10 Myr isochrone with mass ratio equal to unity (dashed–dotted curve). With binary stars and photometric uncertainties taken into account, the broadening of the PMS at regions I and II may still imply an age spread on the order of a million years in Vela OB2. Such an age spread is commonly seen in other young star clusters (e.g., NGC 3603, Beccari et al. 2010; Pang et al. 2013) as a consequence of the propagation of star formation.

A handful of candidate stars are located below the MS locus ($M_G > 4.3$ mag) and correspond to an age older than 50 Myr (black dotted curve in Figure 2 (a)). These stars are located in a more extended region than the majority of the stars in Vela OB2. These stars are likely field stars. As the cluster is located very close to the Galactic plane, field star contamination is significant, and additional photometric cleaning is required (see Section 2.3 in Pang et al. 2020). Therefore, we consider stars bluer than the 50 Myr old isochrone and fainter than $M_G > 4.3$ mag as field stars (black crossed symbols in panel (a) in Figure 2), and we exclude them from the analysis below.

We label the five second-level clusters in Vela OB2 as Huluwa<fnref n="10" /> 1–5, following a Chinese animation series, Huluwa (Calabash brothers). From Huluwa 1 to Huluwa 5, the number of stellar members decreases. Huluwa 1 is the most populous group, which is the Gamma Vel (Pozzo 1) cluster.

## 3.2. Age Spread in the Hierarchical Structures

We construct the CMD for each group in Vela OB2 after cleaning additional contamination in panels (b)–(f) in Figure 2. In order to identify a possible age difference between the five clusters, we enhance the upper PMS region (I) and lower PMS region (II) in the CMD (indicated with arrows in panel (a)). The former is the region where the PMS merges with the MS, and the latter is a region where a significant spread exists. The PMS isochrones are not well established in the lower PMS region (Kos et al. 2019). Previous studies have shown that coeval stars with different accretion rates will induce a spread in luminosity, especially in the lower-mass end of the PMS (Tout et al. 1999; Baraffe et al. 2009; Kunitomo et al. 2017). As can be seen in Figures 2(b)–(f), stars in the lower PMS region spread out to the 10 Myr binary sequence in all groups. Based on the above considerations, the age estimation from this lower PMS region will suffer from large uncertainties.

Stars in the upper PMS regions in Huluwa 1–3 populate the region between the 10 and 20 Myr isochrones (panels (b)–(d)), suggesting that these structures have ages between 10 and 20 Myr. On the contrary, stars in the upper PMS region in Huluwa 4 and Huluwa 5 are mostly younger than the 10 Myr isochrone. However, these two clusters suffer from low number statistics in this region.

To further quantify the age spread among Huluwa 1–5, we adopt and modify the method used in Liu & Pang (2019) and Kos et al. (2019). We define the average distance of each

<footnote n="10">¹⁰ Huluwa means calabash brothers, which is one of the most popular animations in China in the 1980s. The legend begins when the calabashes sequentially fall off from the same stem and transform into seven boys (https://en.wikipedia.org/wiki/Calabash_Brothers). Similarly, the five clusters are sequentially born from the same molecular cloud.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** CMDs for identified structures in Vela OB2. (a) Gray dots represent all member candidates of the top-level structure in Vela OB2 (Huluwa 1–5). Black crosses are field stars. The solid, dashed, and dotted curves in panel (a) are PARSEC isochrones of 10, 20, and 50 Myr, respectively, with solar metallicity and $A_V$ = 0.35. The dashed–dotted curve represents the binary sequence of the 10 Myr isochrone assuming equal mass ratios. The arrows indicate the upper PMS (I) and lower PMS (II) regions. The formats of the curves in panels (b)–(f) are identical to those in panel (a). The two insets show zoomed-in versions of the upper PMS (I) and lower PMS (II) regions.</figcaption></figure>

member in the cluster to the isochrone in the CMD as

$$ \bar{d} = \frac{1}{n} \sum_{k=1}^{n} \sqrt{[\Delta(G_{BP} - G_{RP})]_{i,n}^2 + [\Delta M_{G_{i,n}}]^2}, \quad (1) $$

where $\Delta(G_{BP} - G_{RP})_{i,n}$ and $\Delta M_{G_{i,n}}$ are, respectively, the color difference and magnitude difference between the member star and its nearest point of the isochrone with the k-D tree method (Millman & Aivazis 2011). We study a set of isochrones with solar metallicity, $A_V$ = 0.35 (Cantat-Gaudin et al. 2019b), and age ranging from 0 to 30 Myr with a step size of 0.1 Myr. The isochrone with the smallest value of $\bar{d}$ (solid gray line in Figure 3) is considered as the best-fitted age of the cluster. The value of $\bar{d}$ shows a local minimum and increases toward older and younger ages (Figure 3). The uncertainty of the fitted age $\sigma_{age}$ of each group is taken as half of the age difference (between the two vertical gray lines) at which the average distance reaches $1\sigma$ (gray dashed line) above the minimum $\bar{d}$.

The age uncertainty is comparable to the best-fitted age, which further confirms the existence of an age spread found in Figure 2. Although an absolute age determination for young clusters using isochrone fitting may have an uncertainty of 10%

(Kos et al. 2019), a relative age spread of a few, up to ~10 Myr in Huluwa 1–5, still can be identified in Figure 3. Note that the effects of stellar binarity, the accretion of PMS, and the photometric uncertainty can also position stars above or below the isochrone and contribute to the average distance. Based on our current data, we are unable to disentangle these three contributions. Therefore, the estimated age spread from the uncertainty $\sigma_{age}$ should be treated as an upper limit for the age spread in each group.

From Huluwa 1 to Huluwa 5, the age becomes younger and the age spread decreases. A progression of star formation could have taken place in Huluwa 1–5. The best-fitted age of Huluwa 1 is 12.1 Myr with a spread of 10.3 Myr. Huluwa 5 is the youngest cluster with an age of 7.5 Myr and the smallest age spread of 7 Myr. This age estimation is consistent with the previous findings of Jeffries et al. (2017), who suggested an age of 18–21 Myr for Gamma Vel (Huluwa 1) based on lithium-abundance analysis, and is also consistent with the recent work by Franciosini et al. (2020), who used a lithium-depletion model of 18 Myr with 20% spot coverage and achieved a reasonable fit for Gamma Vel (Huluwa 1).

Considering the age spread in each group, we adopt an age of 10 Myr for Huluwa 1–5 in Vela OB2 and obtain the mass of
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Dependence of the average distance $\bar{d}$ on the cluster age. Colored curves show the average distance $\bar{d}$ between members of each group and an isochrone at a given age as measured in the CMD. The solid gray horizontal line indicates the minimum average distance at the best-fitted age. The horizontal gray dashed line corresponds to the $1\sigma$ deviation above the minimum average distance. The uncertainty of the fitted age is taken as half of the age difference between two vertical dotted gray lines.</figcaption></figure>

each member star from the nearest point on the isochrone via the k-D tree method (Millman & Aivazis 2011). The most massive group, Huluwa 1, contains 723.5 $M_{\odot}$. Huluwa 2 is the second most massive cluster with 467.4 $M_{\odot}$. The total mass for the five clusters in Vela OB2 is 1805 $M_{\odot}$, with an uncertainty of 104 $M_{\odot}$ when a different age of 20 Myr is taken into consideration.

The CMD locations of the members of cluster pair Collinder 135 and UBC 7 (Figure 4) are consistent with an age of 40 Myr (Kovaleva et al. 2020). A binary sequence is clearly visible in the CMDs of both Collinder 135 and UBC 7, which is within the equal-mass binary sequence (dashed-dotted curves). However, no obvious contamination sequence is observed. We do not apply additional cleaning for these two clusters. Adopting an age of 40 Myr for each of the clusters in the pair, the total mass of the cluster pair reaches 445.5 $M_{\odot}$ (Collinder 135: 253.2 $M_{\odot}$; UBC 7: 192.3 $M_{\odot}$).

### 3.3. Membership Cross-match

In total, 3074 members of Vela OB2 are selected after CMD cleaning, 1294 members for Huluwa 1, 743 members for Huluwa 2, 588 members for Huluwa 3, 347 members for Huluwa 4, and 102 members for Huluwa 5. Collinder 135 contains 377 members, and 336 members are identified for UBC 7. Parameters of each of these structures are listed in Table 1. We provide a detailed member list of all clusters in Table 2, with parameters obtained in this study.

Our member list reaches a good agreement with those of earlier studies but provides a more complete sample for Vela OB2. For the largest group, Huluwa 1, we extend the number of members from 437 in Liu & Pang (2019) (ID:2394, also known as LP 2394) and 341 in Cantat-Gaudin et al. (2020) to 1294, where 354 and 325 members are overlapped with the current study, respectively. There are 139 members in Huluwa 1 and 27 members in Huluwa 2 that cross-match with members identified in Jeffries et al. (2014).

Among the members of the cluster pair, we compare our member list to that of Cantat-Gaudin et al. (2020), in which 232 stars matched in Collinder 135 and 118 stars matched in UBC 7.

<figure><img src="image_7.png" /><figcaption>**Figure 4.** CMDs for Collinder 135 and UBC 7. A PARSEC isochrone of 40 Myr (black solid curve) is displayed in panels (g) and (h) with solar abundance and $A_V = 0.13$. The dashed-dotted curve represents the binary sequence of the 40 Myr isochrone, assuming equal mass ratios.</figcaption></figure>

## 4. Kinematical Characteristics of the Hierarchical Structures

### 4.1. Proper-motion Distribution

One characteristic of hierarchical structures in the complex is kinematic coherence between the stellar groups. Such coherence is expected since all substructures were formed from the same molecular cloud and have inherited the internal cloud kinematics (see, e.g., Elmegreen et al. 2000). The PM distribution and contour maps of Huluwa 1–5 and the cluster pair are presented in panels (a) and (b) of Figure 5. Each cluster occupies a distinct but adjacent region in the PM vector diagram, suggesting that their kinematical properties are related. Huluwa 1 exhibits an extended PM distribution. Two overdensities are observed within Huluwa 1 in the PM contour map (see Figure 5(b)). The PM distributions of Collinder 135 and UBC 7 are adjacent but still well separated, which supports a common formation origin for this coeval cluster pair (Castro-Ginard et al. 2018).
<table><caption>**Table 1**

General Parameters of Target Clusters</caption><thead><tr><th rowspan="2">Cluster</th><th>Age</th><th>X<sub>c</sub></th><th>Y<sub>c</sub></th><th>Z<sub>c</sub></th><th>r<sub>h</sub></th><th>r<sub>t</sub></th><th>M<sub>cl</sub></th><th>μ<sub>α</sub> cos δ</th><th>μ<sub>δ</sub></th><th>σ<sub>PMRA</sub></th><th>σ<sub>Decl.</sub></th><th>σ<sub>RV</sub></th><th>Memb.</th></tr><tr><th>(Myr)</th><th>(2)</th><th>(3)</th><th>(4)</th><th>(5)</th><th>(6)</th><th>(M<sub>☉</sub>)</th><th>(9)</th><th>(10)</th><th>(mas yr<sup>-1</sup>)</th><th>(12)</th><th>(km s<sup>-1</sup>)</th><th>(number)</th></tr></thead><tbody><tr><td>Huluwa 1</td><td>12.1–22.4</td><td>-46.3</td><td>-348.6</td><td>-46.1</td><td>14.6</td><td>13.1</td><td>724</td><td>-6.38</td><td>9.33</td><td>0.62<sup>+0.07</sup><sub>-0.05</sub></td><td>1.01 ± 0.10</td><td>1.34<sup>+1.01</sup><sub>-0.84</sub></td><td>1294</td></tr><tr><td>(Huluwa 1A)</td><td>12.1–22.4</td><td>-44.2</td><td>-342.6</td><td>-45.7</td><td>10.8</td><td>10.6</td><td>381.2</td><td>-6.55</td><td>9.72</td><td>0.46 ± 0.03</td><td>0.39 ± 0.03</td><td>0.35<sup>+0.16</sup><sub>-0.20</sub></td><td>675</td></tr><tr><td>(Huluwa 1B)</td><td>12.1–22.4</td><td>-48.2</td><td>-354.9</td><td>-47.0</td><td>13.6</td><td>10.2</td><td>342.8</td><td>-6.20</td><td>8.90</td><td>0.49 ± 0.07</td><td>0.38 ± 0.05</td><td>0.92<sup>+0.41</sup><sub>-0.49</sub></td><td>619</td></tr><tr><td>Huluwa 2</td><td>11.4–21.5</td><td>-45.4</td><td>-391.2</td><td>-63.2</td><td>15.3</td><td>11.3</td><td>467</td><td>-5.55</td><td>8.22</td><td>0.81<sup>+0.11</sup><sub>-0.09</sub></td><td>0.66 ± 0.08</td><td>2.35<sup>+1.08</sup><sub>-0.92</sub></td><td>743</td></tr><tr><td>Huluwa 3</td><td>10.6–19.6</td><td>-64.2</td><td>-383.1</td><td>-70.1</td><td>7.9</td><td>10.5</td><td>373</td><td>-4.74</td><td>8.94</td><td>0.77<sup>+0.17</sup><sub>-0.15</sub></td><td>0.75<sup>+0.13</sup><sub>-0.11</sub></td><td>1.80<sup>+1.20</sup><sub>-0.85</sub></td><td>588</td></tr><tr><td>Huluwa 4</td><td>10.3–20.0</td><td>-64.6</td><td>-335.4</td><td>-12.3</td><td>14.8</td><td>8.3</td><td>181</td><td>-7.15</td><td>10.02</td><td>0.68<sup>+0.24</sup><sub>-0.16</sub></td><td>0.50<sup>+0.18</sup><sub>-0.11</sub></td><td>...</td><td>347</td></tr><tr><td>Huluwa 5</td><td>7.5–14.5</td><td>-95.4</td><td>-341.7</td><td>12.8</td><td>4.7</td><td>5.7</td><td>61</td><td>-7.01</td><td>10.81</td><td>0.49<sup>+0.34</sup><sub>-0.15</sub></td><td>0.47<sup>+0.30</sup><sub>-0.15</sub></td><td>...</td><td>102</td></tr><tr><td>Collinder 135</td><td>40</td><td>-105.8</td><td>-277.7</td><td>-58.0</td><td>9.7</td><td>9.2</td><td>253</td><td>-10.05</td><td>6.17</td><td>0.73<sup>+0.10</sup><sub>-0.09</sub></td><td>0.69<sup>+0.09</sup><sub>-0.07</sub></td><td>1.13 ± 0.53</td><td>377</td></tr><tr><td>UBC 7</td><td>40</td><td>-98.0</td><td>-252.6</td><td>-63.7</td><td>7.6</td><td>8.4</td><td>192</td><td>-9.77</td><td>7.00</td><td>0.46<sup>+0.16</sup><sub>-0.12</sub></td><td>0.51<sup>+0.07</sup><sub>-0.05</sub></td><td>1.36<sup>+0.72</sup><sub>-0.76</sub></td><td>336</td></tr></tbody></table>

**Note.** Note that Huluwa 1A and 1B are two components of Huluwa 1. $X_c$, $Y_c$, and $Z_c$ are the median positions of members in each cluster, in heliocentric Cartesian coordinates after distance correction (Section 5.1). Parameters $r_h$ and $r_t$ are half-mass and tidal radii of each cluster, respectively. The tidal radius is computed using Equation (12) in Pinfield et al. (1998). The quantity $M_{cl}$ is the total mass of each star cluster. Parameters $\mu_\alpha$, $\cos \delta$ and $\mu_\delta$ are the mean PMs for each group. Parameters $\sigma_{pmra}$ and $\sigma_{pmdec}$ are the dispersions of the R.A. and decl. components of the PMs. Parameter $\sigma_{RV}$ is the RV dispersion after correction of binary stars and field star contamination. The latter is not available for Huluwa 4 and 5 owing to the small number of RV measurements. Errors in the velocity dispersion are obtained using the MCMC method.
<table><caption>**Table 2**

Columns for the Table of Individual Members of the Five Clusters in Vela OB2 and the Cluster Pair Collinder 135 and UBC 7</caption><thead><tr><th>Column</th><th>Unit</th><th>Description</th></tr></thead><tbody><tr><td>Cluster name</td><td></td><td>Name of the target cluster</td></tr><tr><td>Gaia ID</td><td></td><td>Object ID in Gaia EDR 3</td></tr><tr><td>ra</td><td>degree</td><td>R.A. at J2016.0 from Gaia EDR 3</td></tr><tr><td>er_R.A.</td><td>mas</td><td>Positional uncertainty in R.A. at J2016.0 from Gaia EDR 3</td></tr><tr><td>dec</td><td>degree</td><td>Decl. at J2016.0 from Gaia EDR 3</td></tr><tr><td>er_Decl.</td><td>mas</td><td>Positional uncertainty in decl. at J2016.0 from Gaia EDR 3</td></tr><tr><td>parallax</td><td>mas</td><td>Parallax from Gaia EDR 3</td></tr><tr><td>er_parallax</td><td>mas</td><td>Uncertainty in the parallax</td></tr><tr><td>pmra</td><td>mas yr<sup>-1</sup></td><td>Proper motion with robust fit in α cos δ from Gaia EDR 3</td></tr><tr><td>erpmra</td><td>mas yr<sup>-1</sup></td><td>Error of the PM with robust fit in α cos δ</td></tr><tr><td>pmdec</td><td>mas yr<sup>-1</sup></td><td>Proper motion with robust fit in δ from Gaia EDR 3</td></tr><tr><td>er_pmddec</td><td>mas yr<sup>-1</sup></td><td>Error of the PM with robust fit in δ</td></tr><tr><td>Gmag</td><td>mag</td><td>Magnitude in G band from Gaia EDR 3</td></tr><tr><td>BR</td><td>mag</td><td>Magnitude in BR band from Gaia EDR 3</td></tr><tr><td>RP</td><td>mag</td><td>Magnitude in RP band from Gaia EDR 3</td></tr><tr><td>Gaiaradial速度</td><td>km s<sup>-1</sup></td><td>Radial velocity from Gaia DR 2</td></tr><tr><td>er_Gaia_rational速度</td><td>km s<sup>-1</sup></td><td>Error of radial velocity from Gaia EDR 3</td></tr><tr><td>Jeffries_rational速度</td><td>km s<sup>-1</sup></td><td>Radial velocity from Gaia/ESO survey (Jeffries et al. 2014)</td></tr><tr><td>er_Jeffries_rational速度</td><td>km s<sup>-1</sup></td><td>Error of radial velocity from Gaia/ESO survey (Jeffries et al. 2014)</td></tr><tr><td>Mass</td><td>M<sub>☉</sub></td><td>Stellar mass obtained in this study</td></tr><tr><td>Xobs</td><td>pc</td><td>Heliocentric Cartesian X coordinate computed via direct inverting Gaia EDR 3 parallax ϖ</td></tr><tr><td>Yobs</td><td>pc</td><td>Heliocentric Cartesian Y coordinate computed via direct inverting Gaia EDR 3 parallax ϖ</td></tr><tr><td>Zobs</td><td>pc</td><td>Heliocentric Cartesian Z coordinate computed via direct inverting Gaia EDR 3 parallax ϖ</td></tr><tr><td>Xcor</td><td>pc</td><td>Heliocentric Cartesian X coordinate after distance correction in this study</td></tr><tr><td>Ycor</td><td>pc</td><td>Heliocentric Cartesian Y coordinate after distance correction in this study</td></tr><tr><td>Zcor</td><td>pc</td><td>Heliocentric Cartesian Z coordinate after distance correction in this study</td></tr><tr><td>Dist_cor</td><td>pc</td><td>The corrected distance of individual member</td></tr></tbody></table>

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

(This table is available in its entirety in machine-readable form.)

Previous studies have succeeded in disentangling Vela OB2 into different hierarchical groups. In panel (c) of Figure 5, we overplot the location of the 11 groups described in Cantat-Gaudin et al. (2019b, hereafter CG19) in the PM distribution. Huluwa 1–5 overlap with seven groups from CG19. Huluwa 1 includes groups A and B of CG19, which is consistent with the PM contour overdensity in panel (b). Groups I, J, and K of CG19 are more distant structures (more than 100 pc from Vela OB2) that are not included in our analysis. Additionally, Huluwa 1–3 relate to groups 4, 2, and 1 in Beccari et al. (2018). Note that 27 members in the Gamma Vel B component (Huluwa 1) from Jeffries et al. (2014; black crosses in the blue-dot region in panel (a)) do not actually belong to Gamma Vel (Huluwa 1). Instead, they are members of another cluster, Huluwa 2 (blue dots in panel (a)).

We follow the PM locations of groups A and B in CG19 and separate Huluwa 1 into Huluwa 1A and Huluwa 1B for the following kinematic analysis, with a straight line in panel (c) of Figure 5. Compared to previous studies, cross-matched members of Huluwa 1 (Gamma Vel) with Liu & Pang (2019; black diamonds) have a larger coverage in Huluwa 1B than Cantat-Gaudin et al. (2020; dark-green crosses). The adjacent distribution and the closeness of PMs in Huluwa 1A and Huluwa 1B resemble the configuration of the cluster pair Collinder 135 and UBC 7.

## 4.2. Radial Velocity Distributions

Earlier studies have found that the most massive cluster, Huluwa 1 (Gamma Vel), features two kinematically distinct coeval stellar populations (Jeffries et al. 2014; Sacco et al. 2015;

Franciosini et al. 2018), which were suggested to correspond to a bound core and an expanding halo (Jeffries et al. 2014), or they may have been generated by a past encounter between Gamma Vel and the nearby cluster NGC 2547 (Damiani et al. 2017). In previous studies, the member list of Gamma Vel B (Huluwa 1B) included members of Huluwa 2 (see Figure 5(a), black crosses); this may result in an overestimation of the computed velocity dispersion. We therefore rebuild the RV distribution of Huluwa 1 using RV measurements from GES (Figure 6(a)), based on members identified in the current paper. Huluwa 1A (~17.1 km s⁻¹; shaded cyan histogram) is distinct from Huluwa 1B (~18.2 km s⁻¹; shaded salmon histogram) with an RV distribution of mean RV offset by 1 km s⁻¹, which supports the earlier claims of two distinct kinematic components in Huluwa 1.

The RV distribution in Huluwa 2 (blue histogram, using RV from GES) generally is more dispersed, with a single peak at 21.0 km s⁻¹. For Huluwa 3–5 (black histogram, using RV from Gaia DR2), the RV profile is much more extended. Note that the number of RV measurements in Huluwa 3–5 is much smaller, and the RV error in Gaia DR2 is roughly a factor of 10 higher than that of GES. For the cluster pair, Collinder 135 has an RV distribution almost identical to that of UBC 7 (Figure 6(b)), with a mean RV of 16.7 and 17.3 km s⁻¹, respectively. The RV similarity between the cluster pair constituents is very similar to that of Huluwa 1A and Huluwa 1B. In both Vela OB2 and the cluster pair, we observe stars with RVs that deviate from the bulk of the population; these are likely stellar binary systems. Since RV information is not included in our member identification, we do not recover these two third-level structures in Huluwa 1 directly from StarGo. Note that groups A and B are identified in different parallax ranges in CG19.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** (a) PM diagram for the five clusters in Vela OB2 and the cluster pair Collinder 135 and UBC 7. Colored dots indicate members of different groups. Black crosses indicate members that have GES RV measurements from Jeffries et al. (2014). (b) Contour map of the PM distribution. (c) Modified PM diagram of the top left panel in Figure 5 of Cantat-Gaudin et al. (2019b). Members of Huluwa 1–5 are overplotted. Each ellipse represents the mean and standard deviation of 11 groups in Cantat-Gaudin et al. (2019b). The dotted straight line divides Huluwa 1 into the two subgroups, Huluwa 1A (red dots) and Huluwa 1B (yellow crosses). Members of Huluwa 1 that are cross-matched with the catalogs of Liu & Pang (2019) and Cantat-Gaudin et al. (2020) are indicated as black diamonds and dark-green crosses, respectively.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Histograms of RVs for (a) Huluwa 1–5 in Vela OB2 and (b) the cluster pair Collinder 135 and UBC 7. The bin size of the histogram is 2 km s⁻¹. Histograms for Huluwa 1 (red) and Huluwa 2 (blue) are obtained from members with GES RVs from Jeffries et al. (2014) only. We separate Huluwa 1 into Huluwa 1A (cyan-shaded histogram) and Huluwa 1B (salmon-shaded histogram) following the division criteria in Figure 5(c). RVs of Huluwa 3–5 and the cluster pair were obtained from Gaia DR 2.</figcaption></figure>

# 5. 3D Space and Mass Distribution

## 5.1. Spatial Configuration

We construct the 3D spatial distribution for each cluster by correcting the positions of individual stars for the pseudo-elongation generated by parallax errors, via the Bayesian method described in Carrera et al. (2019) and Pang et al. (2020, 2021). We

adopt a Gaussian distribution for the distances to the cluster stars (for extended structures, a Gaussian profile only provides a rough estimation) and an exponential profile (Bailer-Jones 2015) for the distances to field stars in the prior. We adopt the median value of cluster-centric distance of individual stars as the scale radius of the Gaussian distribution. The membership probability of each star is 95%, considering a 5% field star contamination rate resulting from member identification (see Section 2.2). All members are located at a distance of 278–399 pc. Therefore, the uncertainty in the corrected distance for stars is 2–3 pc, according to the simulations in Pang et al. (2021, see their Figure 4).

In Figure 7, panels (a)–(c), we show the projections of the members of Huluwa 1–5 in Vela OB2 and the pair Collinder 135 and UBC 7 onto the X-Y, X-Z, and Y-Z planes, after distance correction. Among the hierarchical substructures in Vela OB2, Huluwa 1 is the most massive group (red dots), extending up to 50 pc in the X-, Y-, and Z-directions. Huluwa 2 (blue dots) has an extent of ~60–70 pc especially along the Z-direction. Huluwa 3 (purple dots) is more compact than Huluwa 1 and 2. Huluwa 4 (cyan dots) has an elongated shape, connecting the largest cluster, Huluwa 1, and the smallest cluster, Huluwa 5 (green dots).

Collinder 135 and UBC 7 are separated along the Y-direction, with separate concentrations of stars. Collinder 135 has an extended tail (~50 pc long) along the Y-axis, while the tail of UBC 7 is oriented along the X-direction and is 40 pc in length. The centers of both clusters are separated by only 25 pc.

Similar to the cluster pair, the two components of Huluwa 1, Huluwa 1A and 1B, can be divided along the Y-direction at Y~ − 360 pc. Both Huluwa 1A (red dots in Figure 7) and Huluwa 1B (yellow crosses) have a diffuse distribution without any clear central concentration. The centers of Huluwa 1A and Huluwa 1B (see Table 1) are separated by merely 13 pc, which is significantly closer than the previous value of 38 pc reported by Franciosini et al. (2018). Their spatial separation of less than 30 pc (de la Fuente Marcos & de la Fuente Marcos 2009a), the similarity of their PM and RV distributions, and their almost identical mass suggest that Huluwa 1A and Huluwa 1B may be a coeval cluster pair candidate, resembling the cluster pair Collinder 135 and UBC 7.

Spatial distribution of the members of Huluwa 1–5 mimics a ring-like structure, which is located approximately along the
<figure><img src="image_10.png" /><figcaption>**Figure 7.** (a–c) Distribution of 3D spatial position of members in Huluwa 1–5 and the pair Collinder 135 and UBC 7 in heliocentric Cartesian coordinates, after distance correction using a Bayesian approach. The color-coding is identical to that in Figure 5. We separate Huluwa 1B as yellow crosses and Huluwa 1A as red dots in panel (a). In panels (b) and (c) the red dots represent Huluwa 1 instead. The dashed straight line in panel (c) is the proposed location for the rim of the shell, which is used to separate the upper shell and the lower shell. Yellow stars highlight stars with distance within 1 pc (along the Z-axis) from the rim of the shell (dashed line). Panel (d) is the CMD for Vela OB2 members (gray dots). The stars located at the upper shell region are highlighted as orange circles. Upper PMS stars with $M_G$ between 2.5 and 6 mag and older (bluer) than the 10 Myr old isochrone are indicated as black circles, which are mainly located at the lower shell region (panel (c)). The inset is a zoom-in version of the upper PMS region.</figcaption></figure>

rim of the IRAS Vela shell (Sahu 1992; Cha & Sembach 2000). In Figure 8 we display the 2D positions of the identified members, overlaid on the IRIS image (Miville-Deschênes & Lagache 2005), which shows the gas structures in this region. The ring-like structure is also apparent in 2D projection, along the dense gas structure, which might be at the background. CG19 suggested that star formation in Vela OB2 was triggered by a supernova explosion inside the cavity of the IRAS Vela shell. The supernova might be a member of a 30 Myr old generation of stars (group III in Cantat-Gaudin et al. 2019a).

However, the remnant of the hypothetical supernova that created the IRAS Vela shell has not been detected yet.

In panel (c), we highlight the positions of the higher-mass upper PMS stars with $M_G$ magnitudes between 2.5 and 6 mag and older (bluer) than the 10 Myr old isochrone (black circles in Figure 7(d)) as black circles. They are mostly located in the lower part of the shell structure. Following the spatial distribution of these upper PMS stars (black circles), we define the rim of the shell as a straight dashed line (panel (c) in Figure 7) and define the stars in Huluwa 1–5 located above the dashed line as the upper shell, and those below as the lower shell. All stars of the upper shell (the
<figure><img src="image_11.png" /><figcaption>**Figure 8.** IRAS-IRIS infrared image of the 60 µm band (Miville-Deschênes & Lagache 2005). The members of Huluwa 1–5 and the cluster pair Collinder 135 and UBC 7 are displayed as colored dots. Yellow stars highlight stars located along the rim of the shell, within 1 pc (along the Z-axis) from the rim (dashed line in Figure 7(c)). The color-coding is identical to that in Figure 7. Stars located in the upper shell region (Figure 7) are outlined with orange circles. The high-mass upper PMS stars older than 10 Myr shown in the inset of Figure 7(d) are highlighted with black circles.</figcaption></figure>

colored dots with orange circle outline in Figure 8) form a clear shell structure in the 2D projection. The rim of the shell (yellow stars in Figure 8) appears to agree with the location of the edge of the gas shell at the background.

By plotting stars in the upper shell in the CMD in Figure 7(d) as orange circles, we can see that they are primarily low-mass PMS. We find a deficiency of higher-mass upper PMS stars with $M_G$ magnitudes between 2.5 and 6 mag in the upper shell region (zoom-in inset in Figure 7(d)). It seems that most of these older and massive PMS stars are located at the lower shell region (Figure 7(c)). On the contrary, stars at the upper shell are mostly younger and of lower mass.

## 5.2. Mass Stratification along the Shell

To further quantify mass stratification across the shell, we adopt a regression discontinuity (hereafter RD) estimation. The RD method has been used in a number of fields, such as medicine and economics, and has only recently been applied in astrophysics (Pasquato & Matsiuk 2019).

The RD approach is motivated by the need to identify causal effects in observational studies (Imbens & Lemieux 2008). For example, the causal effect in this case is that the supernova explosion triggers mass stratification across the rim of the shell (red straight line in Figure 9(a)). If the spatial mass distribution is discontinuous at the rim of the shell, this causal effect is very probable. Alternatively, a continuous distribution at the rim means a low probability that the supernova explosion is the cause of observed mass distribution.

We adopt the rim of the shell in Figure 7 (dashed line in panel (c); solid red line in Figure 9(a)). Stars located in the upper rim (orange dots in Figures 7(c) and 9(a)) are located along the Vela IRAS gas shell (orange dots in Figure 8). We hypothesize that the masses of stars on either side of the shell should differ only owing to the influence of the supernova explosion. We carry out a local

linear regression for the mass dependence on the star's distance (D) to the rim at the upper (+D) side and lower (-D) side of the shell (solid curve in panels (b)–(d) in Figure 9). From comparing the intercepts of the local linear regression on both sides, when the confidence level of the intercepts do not overlap on both sides, we conclude that the causal effect is a likely explanation.

A critical issue in the RD approach is the determination of the optimal bandwidth for the regression, which determines the boundary of the kernel function included in the local linear regression. Small bandwidths may result in a large variance of the estimator because few data points are included, while large bandwidths may smooth out the discontinuity. We adopt a prescription to generate the optimal bandwidth by Imbens & Kalyanaraman (2012), which is implemented in the package for the R programming language (Dimmery 2016). We obtain an optimal bandwidth corresponding to 2.08 pc.

Figure 9(b) displays our results in a graphical form. The stellar mass at the lower side of the rim (-D; lower shell) is found to be $0.32 \pm 0.14 M_\odot$ higher than at the upper side (+D; upper shell), with an associated *p*-value of 1.8%. The *p*-value is a measure for the probability that this finding can be explained by a random occurrence. A *p*-value less of than 5% can be interpreted that the proposed causal effect is very likely. The estimated difference in the intercepts remains robust, even when we reduce the bandwidth to half (*p* = 0.45%) or double the value (*p* = 4.3%). When the position of the rim of the shell (solid red line in Figure 9(a)) is varied by moving it upward (green dashed line) or downward (blue dashed line) by 5 pc along the Z-axis, the resulting *p*-value becomes less significant and the intercepts on both sides of the rim overlap (in Figures 9(c) and (d)). We vary the location of the rim (with steps of 1 pc) from the proposed position (solid red line in Figure 9(a)) by 15 pc along the Z-direction and find that the mass discontinuity is most robust at the proposed rim position, where the *p*-value reaches a minimum (Figure 9(e)).
<figure><img src="image_12.png" /><figcaption>**Figure 9.** (a) Distribution of members in Huluwa 1–5 on the Y-Z plane. The size of each circle is proportional to the mass of a star. The red solid line indicates the rim of the shell in panel (c) of Figure 7, which is taken into consideration for the regression discontinuity estimation. Stars in the upper rim of the shell are colored orange, and stars in the lower rim are colored black. The dashed green and blue lines are the rim of the shell (red line) moved upward (+5 pc) and downward (−5 pc) along the Z-axis, respectively. (b) The mass of stars as a function of the distance of stars within 5 pc from the rim of the shell (red solid line in panel (a)). Panels (c) and (d) show the results for the case when the dashed blue and green lines in panel (a) are taken as the rim positions. The stellar mass is obtained in Section 3.2 via the k-D tree method by taking the mass from the nearest neighbors of the 10 Myr isochrone. Therefore, the stellar mass is discretely distributed. Positive distance values represent the upper shell (orange), and negative values represent the lower shell (black). The local linear regression curve is obtained from the RD estimation with an associated 95% confidence interval marked by dashed lines. When the confidence areas do not overlap at the rim (D = 0, −5, 5), it indicates a presence of a discontinuity in the mass distribution on both sides. (e) Dependence of the *p*-value on the distance within 15 pc from the proposed rim of the shell (D = 0; red solid line) along the Z-axis. The *p*-value is computed by moving the rim (the red solid line) along the Z-axis in increments of 1 pc. The green and blue dashed lines indicate positions of 5 pc from the rim. The horizontal line indicates the *p*-value of 5%. Color-coding is identical to that in panel (a).</figcaption></figure>

As Huluwa 1–3 mainly are located in the lower shell region, while Huluwa 4–5 are located in the upper shell, we search for evidence in the mass function in each group to cross-check the results obtained using the RD method. We plot the present-day

mass function for Huluwa 1–5 in Figure 10. The slope $\alpha$ in the mass function $dN/dm \approx m^{-\alpha}$ is constrained by using a linear least-squares fit. The slopes of the mass function in the five groups are in an agreement when the errors are taken into consideration.
<figure><img src="image_13.png" /><figcaption>**Figure 10.** The present-day mass function of Huluwa 1–5, with stellar masses derived from the nearest neighbors of the 10 Myr isochrone via the k-D tree method. The Y-axis is shifted manually with a constant so that the mass function of each cluster is visible. The black dashed line is the fitted mass function $dN/dm \approx m^{-\alpha}$ via linear least-squares fitting. The quantity $\alpha$ is the fitted slope and is indicated for each cluster in the upper right corner.</figcaption></figure>

There is an apparent deficiency of high-mass stars in Huluwa 4–5, which is expected from the Weidner et al. (2013) model, in which more massive clusters tend to have more high-mass stars. The observed maximum stellar mass is 11 $M_{\odot}$ in Huluwa 1. If an age of 20 Myr (upper limit) is adopted for Huluwa 1, the maximum stellar mass for a star remaining in the MS phrase is 12 $M_{\odot}$, while higher-mass stars would have evolved off the MS (Banerjee et al. 2020). This is consistent with the age spread in Huluwa 1–5. However, if the massive OB stars have been expelled from the cluster (Wang et al. 2019), or stars more massive than 20 $M_{\odot}$ have directly evolved into black holes (Banerjee et al. 2020), the observed 11 $M_{\odot}$ does not represent the maximum stellar mass.

Huluwa 1–3 contain more stars above 3 $M_{\odot}$ than Huluwa 4–5. This confirms that more high-mass stars are in the lower shell, where Huluwa 1–3 are located, than in the upper shell, which mainly contains Huluwa 4–5.

### 5.3. Mass Segregation

As high-mass stars are stratified along the rim of the shell (dashed lines in panel (c) in Figure 7), we adopt the Minimum Spanning Tree (MST) method (see Allison et al. 2009, for details) to further quantify the degree of mass distribution in Vela OB2 globally and in the cluster pair. 3D spatial positions of members (Pang et al. 2021) are adopted to avoid the dependence of the results on the viewing angle in 2D.

The MST method compares the minimum path length among the Nth most massive members of a stellar grouping to that of the minimum path length of N random members. The method is insensitive to massive stars with extended configuration like a “ring,” “cross,” or “clump,” unless a higher weight is assigned to the dominant configuration (Olczak et al. 2011). Although massive stars are stratified by the rim of the shell (Section 5.2) in Vela OB2, such extended configuration cannot be detected using the MST method (black histogram in Figure 11(a)).

<figure><img src="image_14.png" /><figcaption>**Figure 11.** The “mass segregation ratio” ($\Lambda_{MST}$) for the 60 most massive members, with a bin size of 12 stars in the whole Vela OB2 region (panel (a)), for each of the individual clusters Huluwa 1–5 (panels (b)–(f)), and for the cluster pair Collinder 135 (panel (g)) and UBC 7 (panel (h)). The dashed line ($\Lambda_{MST} = 1$) indicates an absence of mass segregation. An increasing value of $\Lambda_{MST}$ indicates a more significant degree of mass segregation. The error bars indicate the uncertainties obtained from 100 realizations of randomly selected stars.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 12.** The relative 3D velocity vectors for members, projected onto the X-Y, X-Z, and Y-Z planes. The red, blue, purple, cyan, and green vectors represent the velocities of member stars in Huluwa 1–5, respectively, relative to the mean motion of Huluwa 1. The yellow and pink vectors represent the velocities of member stars in Collinder 135 and UBC 7, relative to the mean motion of Collinder 135. The centers of Huluwa 1 and Collinder 135 are indicated with plus signs. The scale of the velocity vectors is indicated in the upper right corner of each panel.</figcaption></figure>

We again apply the MST method to individual group Huluwa 1–5 (panels (b)–(f) in Figure 11). We find mass segregation down to $1M_\odot$ in the most compact but lowest-mass group, Huluwa 5. The other more massive groups Huluwa 1–3 show no clear evidence of mass segregation. The degree of mass segregation in Huluwa 4, Collinder 135, and UBC 7 is marginal.

Mass segregation in a star cluster can be either primordial or dynamical. Dynamical mass segregation occurs as a consequence of the equipartition of energies of stars during two-body encounters. The timescale of dynamical mass segregation is proportional to the relaxation time, which in turn depends on the number of stars, the half-mass radius, and the total mass of the cluster. Therefore, the smallest group Huluwa 5 has the shortest relaxation timescale of ~ 3.8 Myr, computed using Equation (7.108) in Binney & Tremaine (2008, adopting the half-mass radius, the total mass, and the number of members from Table 1). The age of 10 Myr may already allow for a dynamical segregation taking place in Huluwa 5. Even though primordial segregation may have been present in these groups at early times, the violent binary interactions in the center that expel massive stars or a significant expansion can altogether erase the observable signatures of primordial mass segregation.

# 6. Dynamical State

## 6.1. Expansion in the Hierarchical Structures

Cantat-Gaudin et al. (2019a) report significant expansion within the Vela OB complex (group VII in their paper, corresponding to Huluwa 1–5 in the current study). The expansion rate of Vela OB2 is suggested to be the highest in the Vela-Puppis region (Cantat-Gaudin et al. 2019a). Expansion is also observed in the individual clusters that make up Vela OB2; for example, Franciosini et al. (2018) and Armstrong et al. (2020) detected expansion in population B of Gamma Velorum (Huluwa 1B).

We further investigate the expansion in the 3D velocity space for both Vela OB2 and the cluster pair. The 3D velocity vectors of member stars are projected onto the X-Y, X-Z, and Y-Z planes in panels (a)–(c) of Figure 12. We obtain RVs from Jeffries et al. (2014) for members in Huluwa 1 and Huluwa 2. RVs of the remaining stars are obtained from Gaia DR 2. The

velocity vectors of stars in Huluwa 1–5 are relative to the mean motion of Huluwa 1, $(U, V, W) = (-21.1, -14.7, -2.9) \text{ km s}^{-1}$ in heliocentric Cartesian coordinates. In the cluster pair, the velocities are relative to the mean motion of Collinder 135, $(U, V, W) = (-18.6, -8.6, -12.1) \text{ km s}^{-1}$. A clear sign of expansion is found in Huluwa 1–5, as the hierarchical clusters are seen to recede from the most massive group (Huluwa 1). Expansion is also observed in the cluster pair for the first time, with members moving away from Collinder 135.

To quantify the expansion, we show the relative 3D velocity $(U, V, W)$ of these stars versus the relative distance to the center of Huluwa 1 (left panels) or to the center of Collinder 135 (right panels) along the (X, Y, Z) axes in Figure 13. A positive/negative correlation indicates expansion/contraction, according to the linear expansion model by Wright & Mamajek (2018). Unlike previous studies that use all identified members to probe the signature of expansion, in this study we select stars with 3D velocity value between the 14th and 86th percentiles in each distance bin (filled colored circles) to fit a linear correlation, so as to avoid the influence of binary candidates. We observe a positive relationship between relative velocities and positions, with a 1D expansion rate $\kappa$ (i.e., the slope of the fitted line) ranging from $(6.9 \pm 0.3) \times 10^{-2} \text{ km s}^{-1} \text{ pc}^{-1}$ to $(7.9 \pm 0.5) \times 10^{-2} \text{ km s}^{-1} \text{ pc}^{-1}$ in Vela OB2, which is largest along the X-axis and smallest along the Z-axis. This value is more significant than the expansion rate in the 10 Myr old Scorpius-Centaurus OB2 region (on average $(2.0 \pm 3.5) \times 10^{-2} \text{ km s}^{-1} \text{ pc}^{-1}$; Wright & Mamajek 2018) and in the 6–9 Myr old Orion complex (on average $(1.0 \pm 1.3) \times 10^{-2} \text{ km s}^{-1} \text{ pc}^{-1}$; Swiggum et al. 2021). We find higher expansion rates than those obtained by CG19, but the findings are consistent within their uncertainties. The 1D expansion rate in the cluster pair is much smaller. The 3D expansion rate even indicates contraction, with a negative value for $\kappa$. The expansion is most pronounced along the X-direction, with $\kappa_X = (3.94 \pm 1.40) \times 10^{-2} \text{ km s}^{-1} \text{ pc}^{-1}$. Such anisotropic expansion in Vela OB2 and the cluster pair might be a result of spatially varied stellar feedback (Fujii et al. 2021a, 2021b).

Compared to the expansion in a young cluster with similar age (NGC 2232; see Pang et al. 2020), the expansion in Vela OB2 is more significant and is comparable to that of disrupting clusters with tidal tails (e.g., Coma Berenices; Tang et al. 2019; Pang et al. 2021). The degree of expansion experienced by a
<figure><img src="image_16.png" /><figcaption>**Figure 13.** The relative 3D and 1D velocities from Figure 12 vs. the relative 3D and 1D distances (along the X-, Y-, and Z-axes) to the center of Huluwa 1 (left panels) and to the center of Collinder 135 (right panels). The color-coding is identical to that of Figures 5 and 12. We fit a linear relation (solid black line) to the members with velocities within the 14th and 86th percentiles in each distance bin (filled colored circles). Stars outside this percentile cut (open circles) are potential binary candidates, or stars with large uncertainties in their velocities. The expansion rate $\kappa$ is the slope of the line, which is indicated in the upper left corner of each panel, in units of $10^{-2}$ km s$^{-1}$ pc$^{-1}$.</figcaption></figure>
cluster depends not only on how fast the natal gas is dispersed (impulsive vs. adiabatic) but also on the SFE (Baumgardt & Kroupa 2007; Dinnbier & Kroupa 2020a, 2020b). An SFE of at least 33% is required for the remaining cluster to reach equilibrium again after an impulsive gas expulsion event has occurred (Baumgardt & Kroupa 2007). The high expansion rate in Vela OB2 may indicate not only fast gas removal but also an SFE below ~33%.

## 6.2. Velocity Dispersions

To determine whether these hierarchical structures will remain gravitationally bound after such pronounced expansion, we compute the dispersion of RVs and PMs of each cluster, to quantify their dynamical states. We follow the procedure described in Pang et al. (2021), assuming that the likelihood function for the RV distribution can be modeled with two Gaussian components (cluster members and field stars; Equations (1) and (8) in Cottaar et al. 2012), with the effects of binary broadening and measurement uncertainties taken into account. The PM distribution follows a very similar likelihood function (Equations (1)–(3) in Pang et al. 2018) to the RV distribution, except for binary correction. For Huluwa 1–2 (including Huluwa 1A and 1B), we only use RVs from GES in the computation, while for the other groups we use Gaia DR 2. Due to the relatively small number of RV measurements in Huluwa 4 and 5, we do not calculate RV dispersion for these two clusters. We apply the Markov Chain Monte Carlo (MCMC) method to obtain the best-fit values and the corresponding uncertainties. The final dispersions of the PMs and RVs are presented in Table 1.

The dispersion of the RV distribution tends to be larger than that of the PMs. Though the global expansion after gas expulsion may induce a velocity anisotropy (Baumgardt & Kroupa 2007) in Huluwa 1–3 and the cluster pair, these different dispersions more likely have originated from smaller uncertainties in PM measurements with respect to those of the RVs. Higher-resolution spectroscopy is required to justify this trend.

Based on the velocity dispersion, we estimate the virial ratio (the ratio between the total kinetic energy and the total gravitational energy) for Huluwa 1–3 and for the cluster pair. The resulting virial ratio is larger than 10 (where a value of 0.5 corresponds to a virialized system). This requires at least an additional 20 times more mass to retain the expanding star. Even observational incompleteness cannot conceal such a large amount of mass. Incompleteness often accounts for a few percent of the mass of the cluster (Tang et al. 2019). This is not surprising, since many OB associations are known to be supervirial, with virial ratios ranging from 10 to 1000 (Mel'nik & Dambis 2017; Wright 2020). Over half of the groups have half-mass radii larger than their tidal radii (Table 1, Equation (12) in Pinfield et al. 1998), which further confirms their disrupting states.

We also measure the velocity dispersion for the two components that make up Huluwa 1 (i.e., Huluwa 1A and 1B). Huluwa 1B has a velocity dispersion that is almost three times higher than that of Huluwa 1A. The virial ratio is 3 for Huluwa 1A and 9 for Huluwa 1B. The dispersion error causes 10% uncertainty in the virial ratio. Having 50% of its mass already unbound (the half-mass radius is identical to the tidal radius; see Table 1), Huluwa 1A appears to be in an early stage of disruption, while Huluwa 1B (the half-mass radius is 30% larger than the tidal radius) is already in an advanced stage of disruption. Unlike in previous studies (Jeffries et al. 2014; Franciosini et al. 2018; Armstrong et al. 2020), we find

that Huluwa 1A may not be a bound core. The entire cluster Huluwa 1 is supervirial, although Huluwa 1A is probably dispersing at a slower rate than Huluwa 1B. The velocity dispersion and virial ratio confirm the dispersing state of the Vela OB2 structures and the cluster pair, which may be caused by impulsive gas expulsion with an SFE that is well below 33% (Baumgardt & Kroupa 2007; Dinnbier & Kroupa 2020a, 2020b).

## 7. Star Formation History in Vela OB2: the Supernova Quenching Scenario

Hierarchical structures can be formed in complexes through two channels. The first channel is mild and progressive and occurs via the propagation of star formation (Elmegreen et al. 2000; Kerr et al. 2021). The stellar feedback of the first generation of stars excited turbulence in the ambient molecular cloud and produces density fluctuations that drives the cloud to collapse and form stars (Mac Low & Klessen 2004; Fujii et al. 2021b). This turbulence-driven star formation often occurs along filamentary structures. The second formation channel, on the other hand, is violent, with supernova shocks compressing gas and generating overdensities. Star formation following this mode often takes place along the rim of a shell around the supernova (Kounkel 2020).

We propose an alternative scenario in which the age and spatial configuration of the stars in Vela OB2 (Huluwa 1–5) are a consequence of sequential star formation. The supernova inside the Vela IRAS shell (Cha & Sembach 2000) plays a destructive role in quenching star formation (Krause et al. 2016; Wang et al. 2020). The formation of the 20 Myr old generation in Huluwa 1–3 may have been triggered by the interplay of stellar feedback from surrounding older clusters that powered H II regions by pushing gas to the periphery (Fujii et al. 2021a), such as Collinder 135/UBC 7.

After about 3 Myr, the massive stars in Huluwa 1–3 progressively produced stellar feedback and pushed the gas away from the central clusters and excited turbulence in the molecular cloud. This generated local density fluctuations along the filamentary structures (Mac Low & Klessen 2004; Fujii et al. 2021a) and produced Huluwa 4–5. Simultaneously, the stellar population of Huluwa 1–3 started to expand after the gas expulsion phase; the majority of their members are now located in the lower region of the shell. When Huluwa 4–5 were forming, explosion of a supernova (maybe a member of a 30 Myr old generation of Population III stars in Cantat-Gaudin et al. 2019a) quenched the star formation by efficiently exhausting the gas in the upper shell region, which is primarily the location of Huluwa 4–5. Therefore, most low-mass protostars in this region stopped accreting gas and failed to grow into massive stars, which require a longer accretion time than their lower-mass counterparts (Walker et al. 2021). This violent removal of gas has substantially accelerated the disruption process of Huluwa 1–5.

To further investigate the sequential star formation and supernova quenching scenario, we date back the position of the member stars based on their current positions and 3D velocities (see Figure 14, assuming a linear motion for members; Wright & Mamajek 2018). Only the filled colored circles in Figure 13 are used in this calculation. If Huluwa 1–5 were formed from gas compressed by the expanding shell owing to the shock front of the hypothetical supernova, all stars should return to the central region of the cavity inside the Vela IRAS shell (Kounkel 2020). On the contrary, if stars of Huluwa 4–5 were
<figure><img src="image_17.png" /><figcaption>**Figure 14.** 3D configuration of Vela OB2 and the cluster pair, from the present to past 10 Myr. To trace back −5 and −10 Myr, we only use members with velocities within the 14th and 86th percentiles in each distance bin (filled colored circles in Figure 13) and assume that stars follow linear motions with their present-day observed 3D velocities. The panels in the first row show the current positions of these members. The color-coding is identical to that of Figure 13.</figcaption></figure>

formed from turbulence-driven overdensities, they would
follow the filaments and return to the source that empowered
the turbulence (Huluwa 1–3). As can be seen from Figure 14,
the filament-shape-like Huluwa 4–5 shrunk back to a small
space close to Huluwa 1 in the recent 10 Myr, like the other two
clusters Huluwa 2–3. This evidence supports the hypothesis
that sequential star formation occurred in Vela OB2. The stellar
feedback of the more massive group Huluwa 1–3 triggered
turbulence and formed the less massive clusters Huluwa 4–5.

However, considering the uncertainty in both the position
and the distance of the undetected supernova remnant, the
proposed scenario provides an alternative formation channel
for Vela OB2. We cannot exclude the scenario of supernova

triggering star formation (CG19) or nonuniform gas density
based on the current analysis. Further observational detection
of the supernova and observational constraints on its distance
and position will certainly aid in identifying the formation
history of the stellar structures in this region.

8. Future Predictions Using *N*-body Simulations

In order to investigate the dynamical fates of Vela OB2 and
the cluster pair Collinder 135 and UBC 7, we perform *N*-body
simulations using the NBODY6++GPU code (Wang et al.
2015, 2016). We aim to investigate the mutual interaction
among clusters Huluwa 1–5, and between Collinder 135 and
UBC 7, and determine whether they will eventually merge to become one large structure owing to gravitational interactions between subclusters as suggested in the literature (Goodwin & Whitworth 2004; Allison et al. 2009; Smith et al. 2011; Farias et al. 2018). We use the observed members of each group with their 3D spatial configuration after distance correction (Section 5.1) as the initial positions. To assign each star an individual 3D velocity, we make use of PMs and RVs from Gaia EDR3 and RVs from GES. However, only a small fraction of members have RV measurements.

To correct for the lack of RVs, we assume that the RV distributions of members in each group are Gaussian, with a mean value obtained from the observed average value and the dispersion based on PMs from Gaia EDR3 as a reference value. We evaluate models with three different RV dispersions. First, we study the case in which the RV dispersion is identical to that of the PMs, under the assumption of isotropic motion. We subsequently study models with a velocity dispersion that is half of, or twice, the isotropic value derived from the PM distribution. The RV of each star is randomly assigned for each of the three adopted RV distributions. In our simulations we ignore stellar evolution, stellar binaries, and the tidal effect from the Milky Way that will accelerate the disruption process of the clusters. We also assume that all stellar groupings have been cleared of their gas contents; we therefore do not include any remaining gas in the simulations.

The simulations are run from the present time into the future 100 Myr. Although the Galactic tide is absent, both Huluwa 1–5 and the cluster pair continue to expand and end up in a much more diffuse and eventually unbound structure, irrespective of the choice of adopted RV dispersion (Figures 15 and 16). We find that the relative velocities among Huluwa 1–5, or between Collinder 135 and UBC 7, are higher than their internal velocity dispersions, so that there is no chance for them to form a merged system after gravitational interactions (Gavagnin et al. 2016). Such significant expansion also prevents any groups in Huluwa 1–5 (including Huluwa 1A and Huluwa 1B) from interacting to become a physical pair.

According to the Kovaleva et al. (2020) simulations, the physical separation between Collinder 135 and UBC 7 was closer in the past than in the present. Our N-body simulations show that Collinder 135 and UBC 7 will continue to expand and that their separation will increase in the future 100 Myr. This cluster pair will eventually become dynamically unbound. The same fate seems to happen to the young cluster pair candidate Huluwa 1A and Huluwa 1B. Considering the coevality and the spatial and kinematic proximity between Collinder 135 and UBC 7 and between Huluwa 1A and Huluwa 1B, they were most likely formed from the same molecular cloud, i.e., via the fission of a giant molecular cloud induced by oblique cloud–cloud collisions (Fujimoto & Kumai 1997), but not through tidal capture (van den Bergh 1996) or resonance of Galactic structures (Dehnen 1998; de la Fuente Marcos & de la Fuente Marcos 2009a). The latter processes generally result in cluster pairs of different ages.

Therefore, Huluwa 1–5 in Vela OB2 and the cluster pair Collinder 135 and UBC 7 are all undergoing disruption. Such fast disruption within 10–20 Myr may imply a low SFE (less than 33%), especially in the region of Vela OB2. After the gas evacuation, which was facilitated by the supernova explosion, the shallow potential of the cluster provides no chance to maintain a bound cluster core.

## 9. Summary

1. We have identified hierarchical structures in the Vela OB2 complex and the star cluster pair consisting of Collinder 135 and UBC 7, using data from Gaia EDR3, via the machine-learning algorithm StarGO. Five second-level clusters of stars are kinematically disentangled in Vela OB2 (Huluwa 1, Huluwa 2, Huluwa 3, Huluwa 4, and Huluwa 5). Huluwa 1 is the Gamma Vel (or Pozzo 1) cluster from the literature. In total, 3074 member stars of the Vela OB2 complex are identified. The clusters Huluwa 1–3 (10–22 Myr) are generally older than Huluwa 4–5 (7–20 Myr). Collinder 135 and UBC 7 are of the same age, 40 Myr.

2. The five clusters in Vela OB2 have distinct but coherent PM and RV distributions. Consistent with previous studies, Huluwa 1 is found to consist of two components: Huluwa 1A and Huluwa 1B. However, unlike suggestions in earlier works, Huluwa 1A is not a bound core. Both Huluwa 1A and Huluwa 1B are dispersing structures. We suggest that Huluwa 1A and Huluwa 1B may be the components of a coeval dispersing binary cluster.

3. Collinder 135 and UBC 7 have adjacent PM distributions and nearly identical RV distributions. This confirms previous suggestions that both clusters in this cluster pair have fragmented from the same molecular cloud.

4. The 3D morphology of Huluwa 1–5 resembles a shell-like structure, which is elongated along the Z-axis, up to 100 pc. Stars at the upper part of the IRAS Vela shell are younger than those located at the lower shell. A mass stratification is observed and quantified via the RD method: more massive stars tend to be located in the lower shell, while low-mass stars tend to be more abundant in the upper shell. The stellar masses is $0.32 \pm 0.14 M_{\odot}$ higher in the lower shell region than in the upper shell region.

5. Mass segregation is detected only in the lowest-mass group Huluwa 5. This mass segregation is likely dynamical evolution, i.e., a consequence of two-body relaxation. Marginal mass segregation is found in Huluwa 4, Collinder 135, and UBC 7.

6. Significant expansion is observed in the five clusters Huluwa 1–5 in Vela OB2, with 1D expansion rates ranging from $(6.9 \pm 0.3) \times 10^{-2}$ km s$^{-1}$ pc$^{-1}$ to $(7.9 \pm 0.5) \times 10^{-2}$ km s$^{-1}$ pc$^{-1}$. Expansion in the cluster pair is moderate. 1D contraction in the cluster pair is observed along the Y-direction.

7. The velocity dispersions of Huluwa 1–5 and the cluster pair Collinder 135 and UBC 7 are consistent with a state of disruption. Vela OB2 was probably formed with an SFE that was less than 33%. After a phase of violent gas expulsion, the shallow potential is unable to retain the member stars in a gravitationally bound cluster.

8. We suggest that the formation of the Vela OB2 complex occurred through a process of sequential star formation, with a younger generation of clusters (Huluwa 4–5) formed from turbulence triggered by the older generation of clusters (Huluwa 1–3). The supernova in the cavity of the IRAS Vela shell quenched star formation in the young generation (Huluwa 4–5) along the upper part of the shell by sweeping away their gas content. This rapid phase of gas expulsion further contributed to the significant expansion observed in Huluwa 1–5. However, due to the uncertainty in the position of the supernova, the
<figure><img src="image_18.png" /><figcaption>**Figure 15.** 3D spatial positions for member stars in the Vela OB2 complex (Huluwa 1–5) from N-body simulations, which start from the present-day positions (panels 1) and evolve into the future for 20 Myr (panels 2), 50 Myr (panels 3), and 100 Myr (panels 4). Here we present the model in which the RV dispersions are inferred from those of the PMs, under the assumption of an isotropic velocity distribution. Vela OB2 continues to expand, and its components show no signs of mutual gravitational interactions. The color-coding is identical to that of Figure 5.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure 16.** 3D spatial position for member stars in the cluster pair Collinder 135 and UBC 7 from N-body simulations, which start from the present positions (panels 1) and evolve into the future 20 Myr (panels 2), 50 Myr (panels 3), and 100 Myr (panels 4). Here we present the model in which the dispersion of the RVs is obtained from that of the PMs under the assumption of isotropy. Collinder 135 and UBC 7 continue to expand and show no signs of merging or interaction. The color-coding is identical to that of Figure 5.</figcaption></figure>
supernova triggering star formation scenario cannot be excluded.

9. *N*-body simulations are carried out to predict the dynamical future of the Vela OB2 structure and the cluster pair. All clusters Huluwa 1–5 in Vela OB2 and the cluster pair (Collinder 135 and UBC 7) will expand for the future 100 Myr and will eventually dissolve into the Galactic field. None of the clusters in Vela OB2 will experience mutual interaction with the other components. Although the clusters Collinder 135 and UBC 7 and (possibly) Huluwa 1A and Huluwa 1B are currently a pair, they will not experience a merger event in the future.

We wish to express our gratitude to an anonymous referee for providing comments and suggestions that helped to improve the quality of this paper. X.Y.P. is grateful for the financial support of the research development fund of Xi'an Jiaotong-Liverpool University (RDF-18-02-32). This study is supported by the XJTLU Summer Undergraduate Research Fund (SURF-2021119). X.Y.P. acknowledges grants of the National Natural Science Foundation of China (No. 12173029) and the Natural Science Foundation of Jiangsu Province (No. BK20200252). M.B.N.K. expresses gratitude to the National Natural Science Foundation of China (grant No. 11573004) and the Research Development Fund (grant RDF-SP-93) of Xi'an Jiaotong-Liverpool University (XJTLU). Z.Y. acknowledges funding from the Agence Nationale de la Recherche (ANR project ANR-18-CE31-0017) and the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme (grant agreement No. 834148). M.P.'s contribution to this material is based on work supported by Tamkeen under the NYU Abu Dhabi Research Institute grant CAP³. We express gratitude to Prof. Dr. Pavel Kroupa (University of Bonn) for helpful discussion on the gas expulsion, and Prof. Dr. Simon Goodwin (University of Sheffield) for in-depth communication on the binary cluster evolution. We give thanks to Dr. Long Wang for discussion on the stellar evolution model and Dr. Chien-Cheng Lin for IRIS data download. We thank Prof. Dr. Youzhou Zhou (Xi'an Jiaotong-Liverpool University) for suggestions on 3D morphology interpretation. 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), *TOPCAT* (Taylor 2005), and *StarGO* (Yuan et al. 2018).

## ORCID iDs

Xiaoying Pang https://orcid.org/0000-0003-3389-2263

Zeqiu Yu https://orcid.org/0000-0001-6980-2309

Shih-Yun Tang https://orcid.org/0000-0003-4247-1401

Jongsuk Hong https://orcid.org/0000-0002-5097-8707

Mario Pasquato https://orcid.org/0000-0003-3784-5245

M. B. N. Kouwenhoven https://orcid.org/0000-0002-1805-0570

## References

Allison, R. J., Goodwin, S. P., Parker, R. J., et al. 2009, ApJL, 700, L99

Armstrong, J. J., Wright, N. J., Jeffries, R. D., & Jackson, R. J. 2020, MNRAS, 494, 4794

Arnold, B., Goodwin, S. P., Griffiths, D. W., & Parker, R. J. 2017, MNRAS, 471, 2498

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Astropy Collaboration, Robitaille, T. P., Tollerud, J., et al. 2013, A&A, 558, A33

Bailer-Jones, C. A. L. 2015, PASP, 127, 994

Ballone, A., Mapelli, M., Di Carlo, U. N., et al. 2020, MNRAS, 496, 49

Banerjee, S., Belczynski, K., Fryer, C. L., et al. 2020, A&A, 639, A41

Baraffe, I., Chabrier, G., & Gallardo, J. 2009, ApJL, 702, L27

Baumgardt, H., & Kroupa, P. 2007, MNRAS, 380, 1589

Beccari, G., Boffin, H. M. J., Jerabkova, T., et al. 2018, MNRAS, 481, L11

Beccari, G., Spezzi, L., De Marchi, G., et al. 2010, ApJ, 720, 1108

Bik, A., Puga, E., Waters, L. B. F. M., et al. 2010, ApJ, 713, 883

Binney, J., & Tremaine, S. 2008, Galactic Dynamics (2nd ed.; Princeton, NJ: Princeton Univ. Press)

Cantat-Gaudin, T., Anders, F., Castro-Ginard, A., et al. 2020, A&A, 640, A1

Cantat-Gaudin, T., Jordi, C., Vallenari, A., et al. 2018, A&A, 618, A93

Cantat-Gaudin, T., Jordi, C., Wright, N. J., et al. 2019a, A&A, 626, A17

Cantat-Gaudin, T., Mapelli, M., Balaguer-Núñez, L., et al. 2019b, A&A, 621, A115

Carrera, R., Pasquato, M., Vallenari, A., et al. 2019, A&A, 627, A119

Castro-Ginard, A., Jordi, C., Luri, X., et al. 2018, A&A, 618, A59

Cha, A. N., & Sembach, K. R. 2000, ApJS, 126, 399

Cottaar, M., Meyer, M. R., & Parker, R. J. 2012, A&A, 547, A35

Currie, T., Hernandez, J., Irwin, J., et al. 2010, ApJS, 186, 191

Damiani, F., Prisinzano, L., Jeffries, R. D., et al. 2017, A&A, 602, L1

de la Fuente Marcos, R., & de la Fuente Marcos, C. 2009a, A&A, 500, L13

de la Fuente Marcos, R., & de la Fuente Marcos, C. 2009b, ApJ, 700, 436

de Zeeuw, P. T., Hoogerwerf, R., de Bruijne, J. H. J., Brown, A. G. A., & Blaauw, A. 1999, AJ, 117, 354

Dehnen, W. 1998, AJ, 115, 2384

Dias, W. S., Alessi, B. S., Moitinho, A., & Lépine, J. R. D. 2002, A&A, 389, 871

Dieball, A., & Grebel, E. K. 1998, A&A, 339, 773

Dieball, A., Müller, H., & Grebel, E. K. 2002, A&A, 391, 547

Dimmery, D. 2016, rdd: Regression Discontinuity Estimation, https://CRAN.R-project.org/package=rdd

Dinnbier, F., & Kroupa, P. 2020a, A&A, 640, A84

Dinnbier, F., & Kroupa, P. 2020b, A&A, 640, A85

Efremov, Y. N. 1978, PAZh, 4, 125

Elmegreen, B. G., Efremov, Y., Pudritz, R. E., & Zinnecker, H. 2000, in Protostars and Planets IV, ed. V. Mannings, A. Boss, & S. S. Russell (Tucson, AZ: Univ. Arizona Press), 179

Evans, N. J. I., Dunham, M. M., Jørgensen, J. K., et al. 2009, ApJS, 181, 321

Farias, J. P., Fellhauer, M., Smith, R., Domínguez, R., & Dabringhausen, J. 2018, MNRAS, 476, 5341

Franciosini, E., Sacco, G. G., Jeffries, R. D., et al. 2018, A&A, 616, L12

Franciosini, E., Tognelli, E., Degl'Innocenti, S., Prada Moroni, P. G., & Randich, S. 2020, Mem. Soc. Astron. Italiana, 91, 80

Fujii, M. S., Saitoh, T. R., Hirai, Y., & Wang, L. 2021a, PASJ, 73, 1057

Fujii, M. S., Saitoh, T. R., Hirai, Y., & Wang, L. 2021b, PASJ, 73, 1074

Fujimoto, M., & Kumai, Y. 1997, AJ, 113, 249

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1

Gavagnin, E., Mapelli, M., & Lake, G. 2016, MNRAS, 461, 1276

Gilmore, G., Randich, S., Asplund, M., et al. 2012, Msngr, 147, 25

Goodwin, S. P., & Whitworth, A. P. 2004, A&A, 413, 929

Imbens, G., & Kalyanaraman, K. 2012, The Review of economic studies, 79, 933

Imbens, G. W., & Lemieux, T. 2008, Journal of econometrics, 142, 615

Jeffries, R. D., Jackson, R. J., Cottaar, M., et al. 2014, A&A, 563, A94

Jeffries, R. D., Jackson, R. J., Franciosini, E., et al. 2017, MNRAS, 464, 1456

Kerr, R., Rizzuto, A. C., Kraus, A. L., & Offner, S. S. R. 2021, ApJ, 917, 23

Kos, J., Bland-Hawthorn, J., Asplund, M., et al. 2019, A&A, 631, A166

Kounkel, M. 2020, ApJ, 902, 122

Kounkel, M., Covey, K., Suárez, G., et al. 2018, AJ, 156, 84

Kovaleva, D. A., Ishchenko, M., Postnikova, E., et al. 2020, A&A, 642, L4

Krause, M. G. H., Charbonnel, C., Bastian, N., & Diehl, R. 2016, A&A, 587, A53

Kunitomo, M., Guillot, T., Takeuchi, T., & Ida, S. 2017, A&A, 599, A49

Lindegren, L., Hernández, J., Bombrun, A., et al. 2018, A&A, 616, A2

Liu, L., & Pang, X. 2019, ApJS, 245, 32

Mac Low, M.-M., & Klessen, R. S. 2004, RvMP, 76, 125

Megeath, S. T., Guttermuth, R., Muzerolle, J., et al. 2016, AJ, 151, 5

Mel'nik, A. M., & Dambis, A. K. 2017, MNRAS, 472, 3887
Millman, K. J., & Aivazis, M. 2011, CSE, 13, 9

Miville-Deschênes, M.-A., & Lagache, G. 2005, ApJS, 157, 302

Olczak, C., Spurzem, R., & Henning, T. 2011, A&A, 532, A119

Orkisz, J. H., Pety, J., Gerin, M., et al. 2017, A&A, 599, A99

Palma, T., Gramajo, L. V., Clariá, J. J., et al. 2016, A&A, 586, A41

Pang, X., Grebel, E. K., Allison, R. J., et al. 2013, ApJ, 764, 73

Pang, X., Li, Y., Tang, S.-Y., Pasquato, M., & Kouwenhoven, M. B. N. 2020, ApJL, 900, L4

Pang, X., Li, Y., Yu, Z., et al. 2021, ApJ, 912, 162

Pang, X., Shen, S., & Shao, Z. 2018, ApJL, 868, L9

Pasquato, M., & Matsiuk, N. 2019, RNAAS, 3, 179

Pinfield, D. J., Jameson, R. F., & Hodgkin, S. T. 1998, MNRAS, 299, 955

Pozzo, M., Jeffries, R. D., Naylor, T., et al. 2000, MNRAS, 313, L23

Prisinzano, L., Damiani, F., Micela, G., et al. 2016, A&A, 589, A70

Priyatikanto, R., Kouwenhoven, M. B. N., Arifyanto, M. I., Wulandari, H. R. T., & Siregar, S. 2016, MNRAS, 457, 1339

Rybizki, J., Demleitner, M., Bailer-Jones, C., et al. 2020, PASP, 132, 074501

Sacco, G. G., Jeffries, R. D., Randich, S., et al. 2015, A&A, 574, L7

Sahu, M. S. 1992, PhD thesis, Kapteyn Institute, Postbus 800 9700 AV Groningen, The Netherlands

Smith, R., Fellhauer, M., Goodwin, S., & Assmann, P. 2011, MNRAS, 414, 3036

Swiggum, C., D'Onghia, E., Alves, J., et al. 2021, ApJ, 917, 21

Tang, S.-Y., Pang, X., Yuan, Z., et al. 2019, ApJ, 877, 12

Taylor, M. B. 2005, in ASP Conf. Ser. 347, Astronomical Data Analysis Software and Systems XIV, ed. P. Shopbell, M. Britton, & R. Ebert (San Francisco, CA: ASP), 29

Torniamenti, S., Pasquato, M., Di Cintio, P., et al. 2021, arXiv:2106.00684

Tout, C. A., Livio, M., & Bonnell, I. A. 1999, MNRAS, 310, 360
van den Bergh, S. 1996, ApJL, 471, L31

Walker, D. L., Longmore, S. N., Bally, J., et al. 2021, MNRAS, 503, 77

Wang, L., Kroupa, P., & Jerabkova, T. 2019, MNRAS, 484, 1843

Wang, L., Kroupa, P., Takahashi, K., & Jerabkova, T. 2020, MNRAS, 491, 440

Wang, L., Spurzem, R., Aarseth, S., et al. 2015, MNRAS, 450, 4070

Wang, L., Spurzem, R., Aarseth, S., et al. 2016, MNRAS, 458, 1450

Weidner, C., Kroupa, P., & Pflamm-Altenburg, J. 2013, MNRAS, 434, 84

Wright, N. J. 2020, NewAR, 90, 101549

Wright, N. J., & Mamajek, E. E. 2018, MNRAS, 476, 381

Yuan, Z., Chang, J., Banerjee, P., et al. 2018, ApJ, 863, 26

Yuan, Z., Chang, J., Beers, T. C., & Huang, Y. 2020a, ApJL, 898, L37

Yuan, Z., Myeong, G. C., Beers, T. C., et al. 2020b, ApJ, 891, 39