© The Authors 2023

1. Introduction

Young star clusters and star-forming regions are benchmarks where current theories of star formation and stellar evolution can be tested and validated. However, this validation requires statistically representative samples where biases are minimal. The most common biases affecting the parameters of nearby young clusters are observational and thus related to the survey characteristics, such as its completeness, limiting magnitude, and spatial extent. Methodological biases can also appear, for example, due to cuts in the observational space (Luri et al. 2018) or due to the membership selection. In young clusters, for example, activity related to youth can result in variable luminosity and colour indices, which can impact the membership analysis and derived properties. Moreover, in the analysis of young star clusters and star-forming regions, the remnants of dust and gas from the parent molecular cloud can also extinct the light of the newborn stars and introduce further biases in the inferred parameters of these young populations.

In this article we analyse the properties of the stellar groups in the nearby Perseus star-forming region using the Gaia Data Release 3 (DR3; Gaia Collaboration 2023) complemented with publicly available catalogues of radial velocity from the Apache Point Observatory Galactic Evolution Experiment (APOGEE; Abdurro’uf et al. 2022) and the Set of Identifications, Measurements and Bibliography for Astronomical Data (SIMBAD; Wenger et al. 2000), and photometry from the Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006) and the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS; Chambers et al. 2016). We mine these datasets with Bayesian methodologies developed for the Dynamical Analysis of Nearby Clusters (DANCe) and COSMIC-DANCe projects (P.I. H. Bouy). In particular, we focus on the phase-space (the joint space of 3D positions and 3D velocities), mass, and energy distributions of IC 348 and NGC 1333, which are two of the youngest, nearest, and richest star-forming clusters in the solar vicinity (Luhman et al. 2016). The phase-space structure, mass, and energy distributions of the young Perseus groups provide observational constraints to compare competing star formation scenarios, theories about the origin and variability of the initial mass function, and models for the star formation history of this region.

In Sect. 2 we review the recent publications about the membership, the spatial, kinematic, and mass distributions, as well as the star formation history, of young groups in the Perseus region. In Sects. 3 and 4 we introduce the dataset and methodologies we use. Then, in Sect. 5, we describe our results regarding the membership of the groups, their phase-space structure, mass, extinction, and energy distributions, and their dynamical state. Afterwards, in Sect. 6, we compare our results with those from the literature and discuss the differences and implications. In Sect. 7, we present our conclusions and future perspectives. Finally, Appendix A lists the assumptions that we take throughout this work.

2. Young groups in the Perseus region

2.1. Membership

IC 348 and NGC 1333 have been the subject of several literature studies. Concerning the membership status, we refer the reader to the excellent review by Luhman et al. (2016). These authors not only compiled and assessed the membership status of previously known candidate members but also identified new ones based on several indicators, in particular near-IR spectroscopy. With the updated lists of candidate members, these authors analysed the clusters’ ages, mass functions, disk fractions and spatial distributions. As a result, they identified 478 and 203 candidate members of IC 348 and NGC 1333, respectively (see Tables 1 and 2 of the aforementioned authors). Their survey is complete down to Ks < 16.8 mag and Ks < 16.2 mag, and in sky regions of 14′ and 9′, for IC 348 and NGC 1333 respectively. Shortly after, Esplin & Luhman (2017) obtained spectra of 11 members from Luhman et al. (2016) and confirmed that two and six are candidate members of IC 348 and NGC 1333, respectively. However, only 364 and 93 of Luhman et al. (2016) members of IC 348 and NGC 133, respectively, have Gaia DR3 parallax, proper motions, and photometry.

Cantat-Gaudin et al. (2018) used Gaia Data Release (DR2; Gaia Collaboration 2018) astrometry and a modified version of the Unsupervised Photometric Membership Assignment in Stellar Clusters algorithm (UPMASK; Krone-Martins & Moitinho 2014) to identify candidate members of hundreds of open clusters in the Milky Way. These authors found 144 and 50 candidate members of IC 348 and NGC 1333, respectively.

Castro-Ginard et al. (2018) utilised Gaia DR2 astrometry and the Density Based Spatial Clustering of Applications with Noise algorithm (DBSCAN; Ester et al. 1996) in combination with an artificial neural network to discover 31 new candidate open clusters. In the Perseus region, they found three cluster candidates: UBC4, UBC19, and UBC31, with 44, 34, and 84 candidate members, respectively. They proposed that these clusters are substructures of the Per OB2 complex, although they noticed that UBC4 is located farther away at 570 pc (see their Sect. 5.4).

Ortiz-León et al. (2018) combined observations of the Very Long Baseline Array (VLBA) with Gaia DR2 and measured the distance and kinematics of IC 348 and NGC 1333. Based on 3σ clipping in the independent spaces of parallax and proper motions these authors identified 133 and 31 members of IC 348 and NGC 1333, respectively (see their Table 7), of which 162 have Gaia DR3 parallax, proper motions and photometry.

Luhman & Hapich (2020) identified 12 new candidates for planetary-mass brown dwarfs in IC 348 based on infrared images obtained with the Wide Field Camera 3 of the Hubble Space Telescope. Their candidates have spectral types later than M8, while their faintest candidate reaches down to 4−5 M_Jup. Unfortunately, none of these sources has proper motions nor parallax in the Gaia DR3 catalogue as they are too faint.

Allers & Liu (2020) designed and implemented a medium band near-IR filter to detect low-mass stars and brown dwarfs. Using this filter, these authors survey 1.3 square degrees of IC 348 and NGC 1333 clusters, for which they identify 19 and 9 candidate members; however, only 13 and 3 of these sources, respectively, have Gaia DR3 parallax, proper motions and photometry.

Pavlidou et al. (2021) used the Gaia DR2 data to study the entire Perseus star-forming region. Through successive cuts and clusterings in the astrometric and photometric features, these authors report the discovery of several hundred members of five stellar groups with ages of 1−5 Myr. Their list of members recovers 50% and 78% of Luhman et al. (2016) members of NGC 1333 and IC 348, respectively. In addition, they identify 170, 27, 329, 85, and 302 candidate members of the groups Alcaeus, Autochthe, Electryon, Heleus, and Mestor, respectively. Although these authors claimed the discovery of the previous five groups, 29 members of Electryon and another 29 members of Heleus belong to the UBC 31 and UBC 19 clusters found by Castro-Ginard et al. (2018).

Kerr et al. (2021) used the Gaia DR2 data to identify ∼3 × 10⁴ young stars within a distance of 333 pc. Applying the Hierarchical Density Based algorithm (HDBSCAN; McInnes et al. 2017) to this sample, the authors recover young associations like Orion, Perseus, Taurus, and Sco-Cen. They analyse the star formation history of each group and find evidence of sequential star formation propagating at a speed of ∼4 km s⁻¹. In Perseus, they identified 264 candidate members that were broadly classified into groups 1A, 1B, 2A and 2B, based on cuts in the plane of the sky to separate the eastern populations (Per 1A and Per 1B) from the western ones (Per 2A and Per 2B), and in age to separate the youngest (Per 1A and Per 2A) from the oldest (Per 1B and Per 2B).

Pang et al. (2022) identified members of 65 open clusters using Gaia Early Data Release 3 (EDR3; Gaia Collaboration 2021) and an unsupervised algorithm based on the technique of self-organising maps after proper motions cuts. In the sky region of Perseus, these authors found 211 candidate members in IC 348, 353 in UBC 31, and 80 and 230 in two substructures related to UBC 31, which they claimed as new and called UBC 31 group 1 and UBC 31 group 2, respectively. However, comparing their candidate members with those of Pavlidou et al. (2021), we find that 177 (50%) of UBC 31 belong to Electryon, and 176 (82%) of UBC 31 group 2 belong to Mestor. Except for IC 348, all these groups fall outside the sky region analysed here.

Kounkel et al. (2022) identified 810 members in the Perseus region using Gaia EDR3 and HDBSCAN. In a previous application of HDBSCAN to other star-forming regions (Kounkel & Covey 2019), the authors normalised the data and used the same value for the HDBSCAN parameters of minimum sample size and minimum cluster size. However, in the Perseus regions, they did not normalise the data and used a minimum sample size of 10 stars and a minimum cluster size of 25 stars. As a result, they found 43 Perseus groups, out of which they selected nine based on their absolute colour-magnitude diagram (CMD); the remaining groups were deemed unrelated to the Perseus region. Out of their nine groups, seven correspond to the groups identified by Pavlidou et al. (2021), while the other two correspond to the California cloud and a new one called Cynurus and located outside the region analysed by Pavlidou et al. (2021).

Wang et al. (2022) identified 211 members in the Perseus region by applying astrometric, photometric, radial velocity, and quality cuts to the Gaia EDR3 data and confirming membership with spectroscopy from the Large Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST). In addition to IC 348 and NGC 1333, these authors also identified two subgroups corresponding to the Autochthe one from Pavlidou et al. (2021) and a new one associated with the Barnard 1 cloud.

Lalchand et al. (2022) searched for substellar objects in the Perseus cloud regions of IC 348 and Barnard 5 using narrow-band imaging centred on the water absorption feature at 1.45 μm. They confirm three brown dwarfs in IC 348 and discover the first one in Barnard 5. In addition, these authors used Gaia EDR3 to analyse the distance and proper motions of the Perseus regions of Barnard 1, L1448, and NGC 1333. They confirm that the western part of the region is closer than the eastern part, with Barnard 5 standing alone 100 pc away from the rest of the groups.

2.2. Spatial distribution

The cluster radius of IC 348 was measured by Scholz et al. (1999), Carpenter (2000), and Muench et al. (2003), with all finding consistent values of about 10 to 15 arcmin. However, Muench et al. (2003), studying the radial profile of the cluster, noticed that it showed two subregions (see their Fig. 4): a core within 5′, and a halo between 5′ and 10.33′, with this latter value corresponding to the limit of their survey. Nevertheless, the authors mentioned that they could not exclude the possibility that the cluster extends beyond their survey area coverage to larger radii. Transforming the previous measurements with the cluster distance (320 pc, Ortiz-León et al. 2018), we obtain cluster radii between 0.9 and 1.4 pc, while the core radius of Muench et al. (2003) corresponds to 0.5 pc. Muench et al. (2003) also warned about the possibility of halo sources within the core radius due to projection effects.

The core and halo populations of IC 348 were identified long ago by the pioneering work of Herbig (1998). He identified these populations based on their Hα emissions, age, and spatial location and concluded that the core of IC 348 lies in front of an older population of stars whose Hα has faded away. He also observed an age gradient in which older ages are located at larger distances. Moreover, he proposed that this extended and somewhat older population formed in the same molecular cloud that gave birth to IC 348 and NGC 1333.

Belikov et al. (2002a) created a Compiled Catalogue in the Perseus region with astrometry and photometry of about 30 000 stars distributed in an area of 10° radius. The identified ∼1000 members of the Perseus OB2 complex are distributed in two populations occupying the same sky region but differing in proper motions and distance (see Table 8 of those authors).

Kumar & Schmeja (2007) used deep infrared K-band data to analyse the spatial distribution of IC 348 with the minimum spanning tree (MST) method. They conclude that the stellar population displays a clustered distribution while the substellar one is homogeneously distributed in space within two times the cluster core radius. Although the substellar population is unbounded, it remains within the cluster limits. We notice that these authors establish the cluster limiting radius as the point where the density of the radial profile merges with that of the background, which is not an intrinsic property of the cluster but rather only of its contrast with the background.

Schmeja et al. (2008) analyse the spatial distribution of IC 348 and NGC 1333 by applying the nearest-neighbour and MST methods to infrared data. They found that both clusters are centrally concentrated and assembled from a hierarchical filamentary configuration that eventually built up to a centrally concentrated distribution. They also found that the stellar population of both clusters is mass segregated.

Kirk & Myers (2011) identified substructures in nearby star-forming regions (including IC 348) by applying the MST method to catalogues of young stellar objects from the literature. In all the groups, particularly in the two at IC 348, they found that the maximum mass member is typically more than five times more massive than the median mass member. Furthermore, these massive members are clustered in the central region and are associated with a relatively large density of low-mass stars. Given that all these groups are young, Kirk & Myers (2011) conclude that the observed configurations should be similar to the primordial ones.

Parker & Alves de Oliveira (2017) analysed IC 348 and NGC 133 using the membership lists of Luhman et al. (2016). They estimated the spatial structure, mass segregation and relative local surface density of these two systems using methods like the MST and the 𝒬-parameter. They found that both clusters are centrally concentrated with no evidence of mass segregation. They argue that the results of Kirk & Myers (2011) can be biased by the smaller list of members and their binning method. Afterwards, they compare their results with numerical simulations to estimate possible initial values for the density and velocities.

Ortiz-León et al. (2018) studied the distance and structure of the Perseus region by combining data from the VLBA and Gaia DR2. They estimated distances of 320 ± 26 pc and 293 ± 22 pc to IC 348 and NGC 1333, respectively, and concluded that given the large median value of the parallax uncertainties, which is larger than the parallax dispersion, then the depth of the groups cannot be extracted from these measurements.

Zucker et al. (2018) determined distances to major star-forming regions in the Perseus molecular cloud using a new method that combines Pan-STARRS and 2MASS stellar photometry, astrometric Gaia DR2 data, and ¹²CO spectral-line maps. For IC 348 and NGC 1333, they estimated distances of 295 ± 18 pc and 299 ± 17 pc, respectively. These uncertainties result from simple addition of the statistical and systematic ones reported in their Table 3.

González et al. (2021) analysed the spatial and kinematic structure of four star-forming regions, including IC 348. The authors applied their newly developed Small-scale Significant substructure DBSCAN Detection (S2D2) algorithm to Luhman et al. (2016) list of members (after removing binaries). They found that: (i) the densities of the six identified substructures indicate that these are Poisson fluctuations rather than imprints of star formation sites, and (ii) the members are centrally concentrated in a radial distribution with an exponent between one and two. The results of these authors are thus consistent with a single and centrally concentrated star formation event in IC 348.

2.3. Velocity distributions

Cottaar et al. (2015) analysed the velocity distribution of IC 348 utilising a subpopulation of 152 members with APOGEE radial velocity measurements. They fitted these measurements using Gaussian mixture models (GMMs; see Kuhn & Feigelson 2017, for examples of applications of mixture models in astronomy) with one and two components and found that the second component improved the fit. Their method allowed them to marginalise the contribution of binaries thanks to a binary model that incorporates a wide range of binary configurations into their expected radial velocities. The authors hypothesised that this second component could arise from: (i) contaminants from nearby groups, (ii) a halo of dispersed or evaporated stars, or (iii) the cluster has not yet relaxed to a single Gaussian distribution. The authors argue that the measured velocity dispersion (0.72 ± 0.07 km s⁻¹, or 0.64 ± 0.08 km s⁻¹ with two components) implies a super-virial state unless the gravitational potential has been underestimated by, for example, unaccounted gas. They accounted for the gas and dust mass by adding 40 M_⊙ and 210 M_⊙ as lower and upper limits. They found no evidence of a gradient in the velocity dispersion regarding distance to the cluster centre or stellar mass. However, they found evidence of convergence along the line of sight, which the small cluster rotation (0.024 ± 0.013 km s⁻¹ arcmin⁻¹) cannot explain.

Using APOGEE radial velocity measurements and a similar methodological analysis as that of Cottaar et al. (2015), Foster et al. (2015) determined the velocity distribution of NGC 1333 based on a sample of 70 members. They found a radial velocity dispersion of 0.92 ± 0.12 km s⁻¹, which is consistent with the virial velocity of the group.

Ortiz-León et al. (2018) analysed the velocity distribution of IC 348 and NGC 1333 using VLBA and Gaia DR2 data. They conclude that there is no evidence of expansion or rotations given that the velocities they measure (V_{Exp; IC 348} = −0.06 km s⁻¹, V_{Exp; NGC 1333} = 0.19 km s⁻¹, V_{Rot; IC 348} = [ − 0.16, 0.0, −0.10] km s⁻¹ and V_{Rot; NGC 1333} = [ − 0.10, 0.10, 0.19] km s⁻¹) are smaller than the observed velocity dispersions: σ_IC 348 = 2.28 km s⁻¹ and σ_NGC 1333 = 2.0 km s⁻¹.

Kuhn et al. (2019) studied the internal kinematics of 28 young (≲5 Myr) clusters and associations with Gaia DR2 data. Using proper motions and distances these authors computed transverse velocities in the plane of the sky. After correcting for perspective effects, they derived outward and azimuthal velocities, which are 2D proxies for the internal motions of rotation and expansion. In IC 348 and NGC 1333, these authors found no evidence of contraction or expansion. With respect to rotation, although they found non-zero azimuthal velocities of −0.45 ± 0.21 km s⁻¹ and −0.27 ± 0.23 km s⁻¹ in NGC 133 and IC 348, respectively, they deemed these values not significant under the large values of the observational uncertainties, and identified only one system, Tr 15, as having significant rotation.

2.4. Age distribution

Luhman et al. (2003) determined the age of IC 348 by comparing Baraffe et al. (1998) and Chabrier et al. (2000) theoretical isochrones of 1, 3, 10, and 30 Myr to their observational Hertzsprung-Russell diagram, which they obtained with spectral measurements of effective temperatures and bolometric luminosities. These authors found that the cluster members have ages compatible with 1 to 10 Myr, with a median value of 1−2 Myr (see their Fig. 9).

Bell et al. (2015) obtained an age of 6 Myr for IC 348 based on photometric data from a list of confirmed members from the literature. These authors derived ages comparing the extinction-corrected photometry to main sequence evolutionary models as well as their own semi-empirical pre-main sequence isochrones. They conclude that for star-forming regions younger than 10 Myr, the age estimates from the literature are younger by a factor of two.

Lada et al. (1996) estimated the age of NGC 1333 at 1−2 Myr based on its large fraction of members bearing a disk (61%). Later on, Gutermuth et al. (2008) measured a fraction of 83 ± 11% of members with a disk, confirming the young age of this group.

Pavlidou et al. (2021) derived the ages of the newly identified groups by comparing the photometry of their candidate members to those predicted by theoretical isochrones. According to their results, Authochte is coeval with NGC 1333, while Heleus and Mestor are slightly older, with some of their members close to or below the 1 Myr isochrone. Finally, Electryon and Alcaeus appear older with lower fractions of discs and their members being compatible with the 5 Myr isochrone. We notice that their age estimates were obtained without applying any extinction correction. These authors assumed that extinction was negligible based on the observation that the sequences of their candidate young stars were close to the theoretical isochrones of young stars.

Kerr et al. (2021) used PAdova TRieste Stellar Evolution Code (PARSEC; Bressan et al. 2012; Chen et al. 2015, and references therein), models and Gaia photometry to determine isochrones ages of 6.0 ± 3.2 Myr and 4.7 ± 0.5 for NGC 1333 (Per 1A) and IC 348 (Per 2A), respectively. These authors also derived ages ∼17 Myr for the older and eastern populations of Perseus (Per 1B and Per 2B).

Recently, several groups have estimated isochrone ages for the Perseus groups. Kounkel et al. (2022) estimated group and individual star ages and found that the age of the new Cynurus group is 7 Myr, whereas the rest of the Perseus groups have ages similar to those of Pavlidou et al. (2021). Also, Wang et al. (2022) used ancillary data from several photometric surveys to estimate ages of 5.4 Myr, 2.9 Myr, and 5.7 Myr for IC 348, NGC 1333, and the remaining cloud regions (i.e. Autochthe and Barnard 1 group), respectively. Moreover, Lalchand et al. (2022) estimated ages of 5 Myr for IC 348 and the Barnard 5 group.

2.5. Mass distribution

Muench et al. (2003) derived the mass distribution of IC 348 down to 10 M_Jup based on infrared (JHK) photometry. They found a mass distribution similar to that of the Trapezium, a brown-dwarf-to-stars ratio of 25%, and radial variations of the mass distribution on the parsec scale. They identified two distinct peaks in the mass distribution attributed to the core and halo populations. They mentioned that to reconcile these two peaks, the age of the halo needed to be 5−10 Myr and that although the age gradient reported by Herbig (1998) was in the correct direction, it was not big enough to account for this difference.

Thies & Kroupa (2007) analysed the mass distribution of the Trapezium, Taurus-Auriga, IC 348, and the Pleiades. They found evidence for correlated but disjoint populations of stars on the one hand and very low-mass stars and brown dwarfs on the other hand, which suggests different dynamical histories for both populations. They obtain a ratio of one brown dwarf for every five stars, although, in IC 348, this ratio reaches up to 30%.

Kirk & Myers (2012) analysed the mass distributions of Taurus, Lupus3, ChaI, and IC 348 based on their previous results on the spatial distributions of these star-forming regions (Kirk & Myers 2011). They found that massive stars are more easily located in regions with higher stellar surface density. Their results suggest strong evidence of this effect in Taurus and IC 348, where stars typically have 10−20% higher mean mass in the clustered environments.

Alves de Oliveira et al. (2013) performed a large survey of IC 348 to uncover its brown dwarf population, for this, they used deep optical and near-infrared images of MegaCam and Wide-Field Infrared Camera (WIRCam) to photometrically select candidate members. They also conducted a spectroscopic follow-up of their candidate members, of which 16 new members were confirmed, including 13 brown dwarfs. Five of these new members have L0 spectral types corresponding to masses of about 13 M_Jup. Combining their new members with those from the literature, they constructed the cluster mass distribution and found no significant differences with the mass distributions of other young clusters. Based on a Kolmogorov–Smirnov (KS) test, they conclude that the IC 348 mass distribution is well-fitted by a Chabrier (2003) log-normal distribution. However, we notice that their mass bin at log Mass [M_⊙] ∼ −1.2 shows a deficit with respect to the density predicted by Chabrier (2003) mass distribution. Interestingly, this feature is predicted by the mass distribution of Thies & Kroupa (2007), as can be observed when comparing Fig. 7 of Thies & Kroupa (2007) with Fig. 11 of Alves de Oliveira et al. (2013).

Scholz et al. (2013) determined the mass distributions of IC 348 and NGC 1333 for a wide interval of isochrone ages, models, extinction and distances. They warn about the strong dependence of the results on these parameters and point out the importance of comparing under similar assumptions. They found brown-dwarfs-to-stars ratios of 40% to 50% in NGC 1333 and 25% to 35% in IC 348. This latter value is in agreement with that of Thies & Kroupa (2007). Comparing these two clusters, they found differences in their cumulative distributions that resulted from a relative excess of low-mass stars in NGC 1333. They conclude that the environment plays an important role, with higher-density regions producing larger fractions of low-mass objects, as predicted by gravitational fragmentation models in which filaments fall into the cluster potential.

2.6. Extinction

Foster et al. (2013) investigated the shape of the extinction law in two regions of the Perseus molecular cloud. They combined red-optical and near-infrared images of Megacam and the UKIRT Infrared Deep Sky Survey (UKIDSS) to measure the colours of background stars. They developed a Bayesian hierarchical model to simultaneously infer the parameters of individual stars as well as those of the population and found a strong correlation between the extinction (A_v) and the slope of the extinction law (R_v), which they interpreted as evidence of grain growth. Later on, based on the correlation found by the previous authors, Zucker et al. (2018) adopted an R_v = 3.3 for moderate extinction values of A_v up to 4 mag. These latter authors determined distances to the Perseus molecular clouds using CO spectral-line maps, photometry and Gaia DR2 parallaxes.

Chen et al. (2016) obtained dust emissivity spectral index, dust temperature and optical depth maps of the Perseus molecular cloud from fitting spectral energy distributions to combined Herschel and James Clerk Maxwell Telescope (JCMT) data. They found that the distribution of the dust emissivity spectral index varies from cloud to cloud, which indicates grain growth. This effect was already reported by Strom et al. (1974) based on multiband photometry of 20 IC 348 bright members. Chen et al. (2016) also found evidence of heating from B stars and embedded protostars, as well as outflows.

Zari et al. (2016) derived optical depth and temperature maps of the clouds in the Perseus region based on Planck, Herschel, and 2MASS data. Their maps have resolutions from 36 arcsec to 5 arcmin and a dynamic range indicating that the extinction in this region reaches up to 20 mag in A_K.

Green et al. (2019) determined 3D maps of dust reddening based on Gaia parallaxes and stellar photometry from Pan-STARRS and 2MASS. Thanks to a spatial prior, they obtain smooth maps with isotropic clouds and small distance uncertainties. They made their map available online through the dustmaps package. Later on, Leike et al. (2020) used variational inference and a Gaussian process to infer highly resolved 3D dust maps up to 400 pc. A detailed comparison between the 3D dust maps of Leike et al. (2020) and Green et al. (2019) is shown in Figs. 9–11 of the former. We notice that although the maps of Leike et al. (2020) show better spatial resolution than those of Green et al. (2019), the latter have better 2D sky resolution (see Fig. 11 of Leike et al. 2020).

Doi et al. (2021) used optical and near-infrared stellar polarimetry in combination with Gaia DR2 parallaxes to study the magnetic field polarisation in the Perseus molecular cloud. They found a bimodal distribution in the polarisation angles that identified with foreground and background molecular clouds. The foreground cloud is located at ∼150 pc and has a contribution to the extinction of A_G ∼ 0.3 mag. On the other hand, the background cloud is at ∼300 pc and has a larger contribution to the extinction: A_G ∼ 1.6. Thus, these authors interpret these two clouds as the edges of an ellipsoidal HI shell of about 100−160 pc in size created by the Per OB2 association, with its foreground edge coinciding with the Taurus molecular cloud.

2.7. Star formation history

Lada et al. (1996) suggested that the ongoing star formation in IC 348 and NGC 1333, located at the opposite ends of the Perseus cloud complex, is produced by a similar physical mechanism, with the main difference between both clusters being that IC 348 has produced stars for a longer time than NGC 1333.

Herbig (1998) proposed the following four-point scenario for the star formation history of the Perseus region. First, the formation of stars in this region has been taking place for at least 10−20 Myr, with the latest episodes corresponding to IC 348 and NGC 1333. This first star formation episode created the OB stars o Per and ζ Per together with low-mass members that are expected to be spread over a large region. Most of this region has been emptied of molecular gas and dust except for the Perseus ridge. Second, there is a population of young stars with ages 1−12 Myr that is in and around IC 348. It is composed of the bright A, B, and F stars as well as ‘field’ ones entangled with IC 348 and at its boundaries. The Hα emission lines of these stars have decayed below the limit of detection. Third, within IC 348, there is a population of stars with Hα emission that are probably younger than the population of the previous point but entangled with it. Fourth, in the densest parts of the Perseus ridge, there is ongoing star formation, as suggested by the highly embedded sources.

de Zeeuw et al. (1999) compiled a comprehensive census of OB associations within 1 kpc of the Sun based on HIPPARCOS data. The results of their census are in qualitative agreement with a large-scale star formation scenario in which the Scorpio-Centaurs-Lupus-Crux, Orion, Perseus, and Lacerta associations formed ∼20 Myr ago out of the bubble blown by the high-mass stars of the Cas-Tau association.

Belikov et al. (2002b) considered that the two populations they identified in Belikov et al. (2002a) constitute an example of propagated star formation. It started in the Per OB2b region approximately 30 Myr, continued in Per OB2a 10 Myr ago, and is now in progress on the southern border of Per OB2, where IC 348 is located.

Pavlidou et al. (2021) suggested that the older groups (Alcaeus, Electryon, Heleus, and Mestor) are closer to the Galactic plane with low latitudes (> −19), while the younger ones (NGC 1333 and Autochthe) are at higher latitudes (< −19). They also point out that NGC 1333 and Autochthe are part of the same star formation event due to the similarity in their properties and their close location particularly. These two groups, together with IC 348, are the only ones with ongoing star formation, while the other older groups in the region have stopped forming stars.

Kerr et al. (2021), using Gaia and HDBSCAN, identified two populations in Per OB2: Per A and Per B, corresponding to the western and eastern regions in the sky. They subdivided each of these populations into two distance subgroups, Per-1 at 283 pc and Per-2 at 314 pc, with Per-1A and Per-2A corresponding to NGC 1333 and IC 348, respectively, and Per-1B and Per-2B their corresponding eastern extensions. While Per A is young (NGC 1333: 6.0 ± 3.2 Myr and IC 348: 4.7 ± 0.5) and concentrated, Per B is older (Per 1B: 17.5 ± 0.9 Myr and Per 2B: 17.1 ± 1.1 Myr) and sparser (35 pc). These authors suggest that due to their kinematic similarities, Per 1 and Per 2 most likely formed in the same star-forming process. However, they notice that a continuous star-forming process between the two populations seems unlikely due to the considerable time lag and lack of age gradient. Instead, they hypothesise that in a parallel fashion, the feedback from the first generation (Per B) dispersed the gas of the parent cloud but did not prevent the continuous in-falling flow of external material, which resulted in a new star formation burst (Per A). In the end, this process produced two distinct epochs of star formation in both Per 1 and Per 2.

Bialy et al. (2021) used different indicators (e.g. X-rays, HI, and ²⁶Al) to identify a dust-free cavity between the Perseus and Taurus star-forming regions, which they call the Per-Tau shell. This most likely formed through one or multiple supernova episodes that created a super-bubble that swept up the interstellar medium and created today’s extended shell, with the age of this shell being ≃6−22 Myr. These authors hypothesise that the supernova that created the shell may have had its origin in: (a) a young (< 20 Myr) star cluster with a total mass between 800 and 3300 M_⊙, (b) a single supernova from a dynamically ejected O or B star, or (c) an ultra-luminous X-ray source. They mention that the most likely scenario is the first one, given that evidence supporting the existence of a young (< 20 Myr) population has been found in Taurus and Perseus. In the latter case, this young population corresponds to the Perseus groups identified by Pavlidou et al. (2021).

Zucker et al. (2022) found that almost all the star-forming complexes in the solar neighbour lie on the surface of the Local Bubble, with the young stars showing expansion with their motions being perpendicular to the Bubble surface. Their trace-back analysis supports a scenario in which the Local Bubble was formed by supernovae at approximately 14 Myr age. The only nearby star-forming complex that does not lie at the Local Bubble surface is the Perseus complex, which is related to the Local Bubble through the Taurus star-forming complex and the Per-Tau shell.

Kounkel et al. (2022) analysed three star formation scenarios to explain the observed kinematics of the Perseus region: a supernova explosion, a cloud-cloud interaction, and the first generation of stars from the Per-Tau shell, with the most likely scenario being the collision of two clouds. According to these authors, the evidence supporting the other two scenarios is not conclusive.

Wang et al. (2022) found that the Perseus region also follows the star formation scenario of the expanding Local Bubble (Zucker et al. 2022; Cox & Reynolds 1987), except that in this case, the Per-Tau shell interacts with it and with the Perseus molecular cloud and results in the elongated shape of the latter. The ages and positions of the Perseus groups suggest that Electryon, Heleus, and Mestor are far away from the Per-Tau shell and unrelated to it. On the contrary, Alcaeus, Autochthe, IC 348 and NGC 1333 are near the edge of the shell and may have formed during the same star formation event.

Lalchand et al. (2022) propose a star formation scenario in which a supernova explosion in the Perseus region triggered the star formation of the region through the driving of a HI super-shell. The latter has at its centre the eastern part of the Perseus region.

3. Data

3.1. Gaia Data Release 3

We downloaded¹ the astrometry and photometry of 164 502 Gaia DR3 sources within the sky region: 51° < RA < 59°, and 30° < Dec < 33°, and proper motions within −100 mas yr⁻¹ < pmra < 200 mas yr⁻¹ and −200 mas yr⁻¹ < pmdec < 100 mas yr⁻¹. From these, 163 178 sources (99.2%) have observed proper motions and photometry, which are necessary to apply our membership methodology. Our initial list of members comprises the 194 candidate members, with a probability > 0.5, found by Cantat-Gaudin et al. (2018) on the IC 348 and NGC 1333 open clusters. We use this sample of members due to its purely Gaia origin and the simplicity of its membership algorithm. Although the pre-Gaia sample of Luhman et al. (2016) is the most extensive one from the literature, we do not use it because it contains several contaminants in the form of astrometric outliers (see Sect. 6.1).

We processed the astrometric data by only applying a parallax zero point correction of −17 μas (see Sect. 2.2 of Gaia Collaboration 2021). As stated in Sect. 7 of Gaia Collaboration (2021) the current parallax bias correction as a function of magnitude, colour, and ecliptic latitude is only a tentative recipe.

3.2. Complementary data

We complement the Gaia DR3 data of our candidate members with APOGEE and SIMBAD radial velocities, as well as 2MASS and Pan-STARRS photometry. Out of our 1052 candidate members (see Sect. 5.1), 428, 407, and 149 have radial velocity entries in APOGEE, SIMBAD, and Gaia DR3, respectively. Whenever a source has multiple radial velocity entries, we select based on the following ordered preference: APOGEE, Gaia DR3, and SIMBAD. In the case of APOGEE, we used as radial velocity uncertainty the dispersion of several measurements (i.e. the VSCATTER column), except when it was zero, in which case we used the individual uncertainty (i.e. the VERR column). In some cases, the radial velocity catalogues report missing or zero value uncertainties for sources with non-missing radial velocity. Thus, we processed the radial velocities as follows. If the star has either a missing uncertainty but a non-missing value or an uncertainty larger than 50 km s⁻¹, then we set the uncertainty to 50 km s⁻¹. This large value diminishes the contribution of the source and avoids the need to discard it. Similarly, if the uncertainties are smaller than 0.01 km s⁻¹, we replace them with the latter value. This uncertainty soil avoids convergence issues in our kinematic inference methodology (see Sect. 4.2). After the previous processing, a total of 626 (60%) of our candidate members have radial velocity measurements, with a median uncertainty of 0.75 km s⁻¹.

In addition, we query the HIPPARCOS (Perryman et al. 1997) data in the same sky region as defined above, and we find 47 sources. Out of these, 46 have a cross-match in Gaia DR3, and the remaining one has parallax and proper motions (μ_α = −139.6 ± 0.89 mas s⁻¹, μ_δ = 20.19 ± 1.05 mas s⁻¹, ϖ = 18.32 ± 0.91 mas) clearly incompatible with the Perseus groups. Due to the previous reason, we assume that the Gaia DR3 of the Perseus region is complete on the bright side.

4. Methods

The following sections describe the methodology that we use to determine the properties of the stellar content of the Perseus region. We start by describing the membership methodology, and afterwards, we describe the steps to identify the distinct physical populations, as well as their properties: empirical isochrones, magnitude and mass distributions. We base the inference of the properties of a physical group on the obtained list of members (see Assumption 1).

4.1. Membership selection

We determined members of the Perseus star-forming region using the Miec code (Olivares et al. 2021), which is an improvement over the Bayesian hierarchical methodology developed by Olivares et al. (2018b) and focuses on the analysis of extincted nearby stellar clusters. Briefly, Miec is a statistical model that describes the observed astrometry and photometry of possible hundreds of thousands of sources in a sky region that encompasses an open cluster. It delivers candidate members of the open cluster as well as its astrometric (proper motions and parallax, if available) and photometric distributions (colour-index, and photometric bands). The likelihood of the data is a mixture of the field and cluster models, where the former consists of independent and multivariate GMMs in the astrometric and photometric spaces, and the latter is also made of a GMM in the astrometric space and a multivariate Gaussian in the photometric one. The median value of the latter corresponds to the cluster photometric sequence, in which each photometric band is modelled by a spline function of the colour index. The model is Bayesian because it infers the posterior distribution of the cluster parameters given the likelihood of the data and the prior distribution. This latter is constructed from the initial list of members. Once the posterior distributions of model parameters have been inferred (through Markov chain Monte Carlo methods) the cluster membership probability of each source in the dataset is computed using Bayes’ theorem and the cluster and field likelihoods as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}& \mathrm{Probability}(\mathrm{cluster}|\mathrm{data},{\mathcal{M}}_{\mathrm{cluster}},{\mathcal{M}}_{\mathrm{field}})=\hfill \\ \hfill & \frac{\mathcal{L}(\mathrm{data}|{\mathcal{M}}_{\mathrm{cluster}})·\mathcal{P}(\mathrm{cluster})}{\mathcal{L}(\mathrm{data}|{\mathcal{M}}_{\mathrm{cluster}})·\mathcal{P}(\mathrm{cluster})+\mathcal{L}(\mathrm{data}|{\mathcal{M}}_{\mathrm{field}})·\mathcal{P}(\mathrm{field})},\hfill \end{array}\end{array}$$(1)

where ℒ, ℳ, and 𝒫 stand for likelihood, model, and prior, respectively. We use as prior probabilities for the field and cluster their fractions of sources in the entire dataset.

The Miec code has been designed for open clusters, and thus, it models the astrometric features of the representation space using GMMs in which all Gaussian components share the same mean value, this is, they are concentric (for more details see Olivares et al. 2018b, 2021). However, the members of star-forming regions have proper motions and parallax distributions that are not necessarily concentric (see, for example, the parallax and proper motions of the Taurus candidate members depicted in Fig. 5 of Galli et al. 2019). For this reason, we modified the Miec code to deal with the dispersed populations present in star-forming regions by allowing non-concentric GMMs in the proper motions and parallax features.

In the Miec methodology, the representation space (i.e. set of observable features) is of paramount importance since it allows the disentanglement of the field and target populations. Given the known issue of the overestimated GaiaBP flux for faint red sources (see Sect. 8.1 of Riello et al. 2021), we used as colour index G − RP instead of BP − RP. Thus, our choice for the representation space comprises the following Gaia features: pmra, pmdec, parallax, G − RP, BP, and G. The Miec code requires that the spline functions describing the cluster photometric sequence be injective functions of the colour index (Olivares et al. 2018a, 2021). However, the previous representation space only allows for this condition to be fulfilled for G values fainter than 5−6 mag. For this reason, we searched for Perseus candidate members brighter than this magnitude limit using only their astrometric membership probabilities, which are also delivered by Miec. We classified sources brighter than G ∼ 5 mag as candidate members if their astrometric membership probability is larger than 3σ (0.997). We chose this highly conservative probability threshold given that, in these cases, the discriminant photometric information is not taken into account.

In addition to the membership probabilities, the code delivers the astrometric and photometric distributions of the target population. While the astrometric distributions are multivariate mixtures in the joint space of proper motions and parallax, the photometric ones are multivariate mixtures in the joint space of colour index and photometric bands. More details about methods for obtaining the astrometric, colour index, and magnitude distributions can be found in Olivares et al. (2018b).

Since the Perseus region contains several dust and gas clouds, we used the extinction module of Miec (see Sect. 2.2 of Olivares et al. 2021). Briefly, this module permits the extraction of the extinction-free population parameters (i.e. those that describe the cluster’s or group’s colour-index and magnitude distributions) by marginalising the possible extinction values, A_v ∈ [0, A_v, max] of each source. For the maximum extinction value, A_v, max, we used the 3D extinction map of Green et al. (2019)² at the group distance. We preferred the maps of Green et al. (2019) over those of Leike et al. (2020), given their better 2D sky resolution and because we are interested in individual stars rather than in the 3D structure of the dust clouds.

As explained in Olivares et al. (2021), the extinction module of the Miec code faces two main caveats: an increased contamination rate due to sources with missing values and a reduced recovery rate in sources with high extinction values. Given that our dataset comprises sources with complete astrometry and photometry (see Sect. 3), that the mean value of extinction to the Perseus region provided by the 3D extinction map is A_v, max ∼ 3 ± 2 mag, and that we classified members based on an optimum probability threshold that is optimised as a function of the G magnitude, then we expect that the performance of the Miec code in our specific conditions will be better than the extreme conditions reported in Table C.3 of Olivares et al. (2021). In other words, we expect a recovery rate better than 87% and a contamination rate less than 7%.

Although the physical members of distant open clusters can be identified using clustering methods working in the proper-motions-parallax space, nearby clusters and dispersed stellar populations extending several degrees on the sky are affected by projection effects that distort the proper-motions-parallax space and make difficult the identification of their members. Thanks to the improved non-concentric GMM, the Miec code is now flexible enough to accommodate possible distortions in the proper-motions-parallax space created by these projection effects. Nonetheless, these distortions increase the mixing of the physical groups in the proper-motions-parallax space and difficult their disentanglement. This effect can be seen in Fig. 10 of Pavlidou et al. (2021), where the proper motions of the Perseus groups are heavily mixed. To disentangle these populations, we iteratively ran the Miec and Kalkayotl (Olivares et al. 2020) codes. The first one identifies the candidate members in the astro-photometric space, while the second one separates the group in the phase-space (more below).

4.2. Phase-space structure

We analysed the phase-space distribution of the Perseus star-forming region using the Kalkayotl code (Olivares et al. 2020, a new version of the code is in prep.). This code implements a Bayesian hierarchical model that allows the joint inference of stellar positions, velocities, and population parameters without imposing a fixed prior. On the contrary, the code allows us to test different 1D (distance), 3D (positions) or 6D (positions+velocities) prior families and infer their parameters based on the Gaia data. Moreover, it corrects for the parallax and proper motions angular spatial correlations and zero points using the values provided by Lindegren et al. (2021). We notice that the output phase-space Cartesian coordinates returned by Kalkayotl are in the equatorial ICRS reference system rather than a Cartesian Galactic one. Throughout the rest of this work, unless stated otherwise, we use this equatorial reference system and the names X, Y, and Z for the 3D positions and U, V, and W for the 3D velocities.

As mentioned in the previous section, the identification of the physical groups is a mandatory step for the subsequent astrophysical analyses. Given that the region is known to host several populations (see Sect. 1), we use the Kalkayotl code to probabilistically disentangle the possible physical groups.

We modelled the stellar positions and velocities of the Perseus star-forming region using 6D GMMs (see Assumption 2). We classified the candidate members into the physical groups by maximising membership probability to each Gaussian component. It is important to notice that the modelling of the phase-space structure is computationally expensive because the number of inferred parameters, N_p, grows linearly with the number of sources: N_p = (6 × N_s)+(28 × N_c − 1), with N_s and N_c the number of sources and Gaussian components, respectively. In the latter equation, the first term corresponds to the parameters of the individual sources, with six phase-space coordinates for each one of them, and the second term corresponds to the global or population parameters. In these latter, each Gaussian component needs 28 parameters: 21 for the covariance matrix, six for the median, and one for the weight. Given that the components’ weights are restricted to add to one, there is one non-free weight.

Our iterative approach to identifying candidate members and physical groups proceeded as follows. Once the candidate members of the entire region were identified by the Miec code, we used the Kalkayotl code to fit their observables using 6D GMMs with one to six components. We considered that several Gaussian components pertained to the same physical group if their medians are mutually contained within one Mahalanobis distance (see Assumption 2). Otherwise, each Gaussian component corresponds to a physical group. We rejected as physical groups those Gaussian components for which: (i) the Hamiltonian Monte Carlo sampler did not converge (see Olivares et al. 2020), or (ii) the contribution to the mixture was ≲5%. We recursively fitted 6D GMM to each physical group until Assumption 2 (see Appendix A) was fulfilled. This recursive fit allowed us to reject non-physical members in the 3D space that, due to projection effects, had similar astrometric features as the bulk of the group. Once the physical groups were disentangled, we joined the list of members of each of them to the field population and ran again the Miec code independently on the resulting dataset. The independent run of Miec on each identified group ensures that the candidate members’ uncertainties (the photometric ones in particular) are propagated into the group’s empirical isochrone and mass distribution (see Sects. 4.3 and 4.4). We iterated the procedure until the number of candidate members of each physical group converged under Poisson uncertainties.

To prevent the radial velocity of unresolved binary stars from biasing the group-level parameters (e.g. internal velocity dispersion) of the Perseus groups, we made an additional run of Kalkayotl using the final members of each group as the input list. This time, however, we set as missing the radial velocities of sources lying more than 3σ away from any of the group’s mean U, V, and W space velocities. Setting the radial velocity as missing rather than removing the source means we do not have to discard the astrometric information of these sources.

Finally, we also analysed the internal kinematics of the identified physical groups by searching for evidence of their expansion or rotation. To do this, we computed the dot and cross product of the positions and velocity vectors of each candidate member on the reference system of its parent group. The average values of these vector products are proxies for the expansion and rotation rates of stellar groups (see Galli et al. 2019, 2021).

4.3. Empirical isochrones and age estimates

The empirical isochrones of the physical groups are inferred from the data and delivered by Miec as a by-product. These empirical isochrones are cubic spline functions that model the mean value of the BP and G magnitudes as functions of the colour index G − RP. We notice that thanks to the use of the extinction module (see Sect. 4.1), these empirical isochrones are free of extinction.

We estimated the age of each physical group by comparing its extinction-free empirical isochrone with the theoretical ones from the PARSEC (Pastorelli et al. 2020; Marigo et al. 2013), MESA Isochrones & Stellar Tracks (MIST; Dotter 2016; Choi et al. 2016), and BT-Settl (Allard 2014) models. We are aware that this dating method provides only a rough estimate of the group age due to well-known issues of theoretical isochrones to reproduce the observed colour-magnitude diagram of young clusters (e.g. Bell et al. 2015; Bouy et al. 2015; Binks et al. 2022).

4.4. Mass distributions

We inferred the mass distribution of each physical group using two different methods, but both based on the PARSEC, MIST and BT-Settl theoretical isochrones at each group’s estimated age (see Sect. 4.3). The first method uses the Sakam code (Olivares et al. 2019) to independently infer the mass of each candidate while the second method transforms the group’s magnitude distributions delivered by Miec into mass distributions.

None of the theoretical isochrone models that we used fully covers the magnitude interval of our candidate members. Thus, when used independently, these models introduce border artefacts in the resulting mass distributions. We overcame this problem by computing a unified theoretical model that we call PMB (which stands for PARSEC-MIST-BT-Settl). In this, we fitted cubic splines to the grid values of mass and magnitudes provided by the three theoretical models. We used cubic splines because they provide continuous derivatives of the magnitude-mass relations and thus avoid the typical problems of simple polynomials (i.e. Runge’s phenomenon).

Sakam is a Bayesian inference code that samples the joint posterior distribution of mass, A_v and R_v of individual stars based on theoretical isochrones and the star’s distance (see Sect. 4.2) and available photometry (see Sect. 3.2). As prior distributions of the mass, A_v and R_v we used the Chabrier (2005) distribution, a uniform distribution (A_v ∈ [0, 10] mag), and a Gaussian distribution (R_v ∼ 𝒩(3.1, 0.5)), respectively. Once the posterior distributions of all the candidate members were inferred, we computed the group’s mass distribution as a kernel density estimate on the aggregated mass samples of all the group members. In the second method, we used the mass-magnitude relations provided by the unified theoretical isochrones (at the group’s estimated age; see Sect. 4.3) together with the group distance to transform Miec’s magnitude distributions of each group into mass distributions. We notice that the theoretical mass-magnitude relations are not one-to-one and have abrupt changes of slope in the 1.5−2.5 M_⊙ mass interval, resulting in mass distributions with a large scatter in this region (see the discussion in Sect. 6.4).

Working with the previous two methods offered the following advantages. First, the Sakam method allowed us to study possible variations in the R_v value across the Perseus groups. These variations have already been suggested in the literature (see, for example, Foster et al. 2013; Zucker et al. 2018). On the other hand, the Miec method has two advantages over the Sakam one. First, it obtains the group’s magnitude distributions using the entire dataset, weighting each source by its membership probability to the group. This approach removes the sample bias introduced when working on the subsample of the most probable group members. Second, it does a full propagation of the model uncertainties and observational uncertainties to the magnitude and mass distributions, whereas the Sakam method only propagates the photometric and distance uncertainties. Although the Miec method offers a statistically more robust approach than that of the Sakam method, none of them is a perfect solution to the mass distribution inference problem. Nonetheless, until the arrival of a complete and spectroscopically confirmed and characterised list of the group’s members (see Assumption 1), the comparison of these two methods offered what we considered to be the best strategy to derive the mass distribution of the groups.

4.5. Dynamical analysis

We performed a dynamical analysis of the Perseus groups based on the source and group level parameters inferred with the methods presented in the previous subsections. We determined the dynamic state of each physical group with two methods.

The first method takes the mass, position, and velocity posterior distributions of each candidate member and computes its energy distribution with respect to its parent group, under the assumption that the latter are self-gravitating (see Assumption 8). We propagated uncertainties by taking samples from the posterior distributions of each candidate member and computing the energy of each sample as follows:

$$\begin{array}{c}\hfill E=\frac{1}{2}m·v^2-\frac{G·M·m}{r},\end{array}$$(2)

where r and v are the distance and speed in the reference system of the stellar group, M is the total group’s mass enclosed within the distance r from its centre, m is the sample’s mass, and G is the gravitational constant. To obtain r and v in the reference system of the stellar group, we used the population parameters delivered by Kalkayotl (see Sect. 4.2).

The second method compares the observed velocity dispersion of each group with the theoretical one expected if the stellar system were at virial equilibrium. To compute the latter, we followed the approach that Cottaar et al. (2015) took in the analysis of IC 348 (see their Sect. 4.3.3). Briefly, these authors assume that the velocity dispersion at virial equilibrium, σ_vir, can be estimated using the total mass of the cluster M, its half-mass radius, r_hm, and a structural parameter, η, (see Assumption 4) following Eq. (4) of Portegies Zwart et al. (2010). Moreover, we followed Cottaar et al. (2015) additional assumption that these parameters can be obtained by fitting an Elson et al. (1987) profile (hereafter EFF) to the 2D stellar number density of the system (see Assumption 5). We did this fitting for the Perseus physical groups with the free and open-source code PyAspidistra (Olivares et al. 2018a). In addition, given that our methods deliver mass and 3D positions for each member in the Perseus groups, we also estimated the half-mass radius by finding the radial distance at which the group’s mass reaches 50%.

Furthermore, we notice the following two aspects of the Perseus star-forming region. First, it is known that this region is still embedded in the dust and gas of its parent molecular cloud and that the contribution of this non-stellar mass to the total mass of the groups is non-negligible. Cottaar et al. (2015) account for this non-stellar mass (see Sect. 2.3), assuming that the dust and gas still follow the observed distribution of stars and that its total mass contribution has lower and upper limits equal to 65% (80 M_⊙) and 169% (210 M_⊙), respectively, of the total stellar mass of IC 348. Here, we also took the previous two assumptions and extended them to the rest of the Perseus groups (see Assumption 6). Second, it is known that the fraction of binary systems in open clusters varies between 11% and 70% (Sollima et al. 2010), with the fraction of unresolved binaries between 12% and 20% (e.g. Jadhav et al. 2021). However, our methodologies are unable to identify and infer the mass of these possibly unresolved binaries. Therefore, it follows that our gravitational potential will be underestimated due to the unaccounted mass of the unresolved binaries. Thus, we corrected for this bias by increasing by 20% the mass contribution of the individual stars when computing the gravitational potential of the groups (see Assumption 7).

5. Results

5.1. Membership

We iteratively apply the Miec and Kalkayotl codes (as described in Sects. 4.1 and 4.2) to the Perseus dataset (see Sect. 3). In the first iteration, the code recovered 920 candidate members, and after successive iterations utilising the extinction module, we recover 130 more candidate members. Our search for bright (G > 5 mag) members delivered only two astrometric candidate members: ζ Per and o Per, with astrometric membership probabilities of 0.99989 and 0.99988, respectively, which implies a ∼4σ discovery.

According to our iterative methodology, the final 1052 candidate members (see Table C.1) are distributed into eight statistical groups (see Sect. 5.2). Table 1 shows the names, number of members, mean distance, age, and mass estimates of these groups. In Sect. 5.2.1 we show that two of these statistical groups pertain to the same physical one, thus effectively reducing the number of physical groups to seven. We identify the well-known IC 348 (with its core and halo) and NGC 133 young clusters (see Fig. 1) and three of the recently discovered populations of Pavlidou et al. (2021): Heleus, Alcaeus and Autochthe. In addition, we discover a putatively new young physical group of ∼7 Myr and 191 candidate members that is composed of a core and halo populations. Following the nomenclature style of Pavlidou et al. (2021), we call this group Gorgophone. In Sect. 6.1, we present a detailed comparison between the candidate members that we find in this work and those from the literature.

Fig. 1. Sky coordinates of the Perseus candidate members. The colour code shows the probabilistic classification, and the background image shows the thermal dust emission (545 GHz) from Planck Collaboration Int. LVII (2020).

5.2. Phase-space structure

As explained in Sect. 4.2, we infer the phase-space structure of the Perseus groups by fitting 6D GMMs using the Kalkayotl code. We jointly inferred the parameters of all candidate members and chose the model with six components as the best one. We based this decision on the convergence properties of the sampler and the weights of the components. Models with more than six components resulted in inefficient sampling and negligible weights for the additional components. Figure 2 shows the inferred phase-space coordinates of the candidate members as well as the six-components GMM (the orange lines show samples from the posterior distribution of the one-sigma covariance matrices). The colour code shows the probabilistic classification of each of the components (i.e. IC 348 core, IC 348 halo, Alcaeus, Heleus, Gorgophone, and NGC 1333+).

Fig. 2. Cartesian equatorial (ICRS) positions (top panel) and velocities (bottom panel) of the 1052 Perseus candidate members. The colour code shows the probabilistic classification, and the orange ellipses depict samples from the posterior distribution of the group level parameters corresponding to the one-sigma covariance matrix, and are centred at the mean position and velocity of the groups.

The joint inference with the 1052 candidate members, our method was unable to completely disentangle the NGC 1333 members from those of Autochthe (see the group called ‘NGC 1333+’ in Fig. 2). To overcome this problem and to break any possible entanglement, we iteratively fitted two-component GMMs to the candidate members of each identified group. At the end of this hierarchical-tree exploration, we find out that all the groups except for Gorgophone required only one Gaussian component to describe their phase-space structure. In the rest of the groups, the additional components showed negligible weights ≤5% and convergence issues. In the two-component GMM of Gorgophone, the weights of the additional component were non-negligible (> 5%).

Our iterative methodology delivers phase-space parameters of the Perseus groups as well as posterior distributions of the positions and velocities of each candidate member. The group-level parameters (mean and standard deviation) of the identified groups are shown in Tables 2 and 3. As an example of the inferred group- and source-level parameters, Fig. 3 shows the Cartesian (ICRS) positions and velocities of the Alcaeus group. The dots and error bars show the mean and standard deviation of the posterior distributions of the source-level parameters (i.e. 6D Cartesian coordinates of each candidate member), while the orange ellipses show samples from the posterior distribution of the group-level parameters (i.e. the one-sigma covariance matrix centred at the mean position and velocity of the group). The total spatial and velocity dispersions of the identified groups are shown in the second and third columns of Table 4. As can be observed from this table, the most distant Heleus group is also the most spread in the XYZ space. On the contrary, the core of IC 348 is the most compact one, with only a 0.66 pc radius.

Fig. 3. Cartesian equatorial (ICRS) positions (top panel) and velocities (bottom panel) of the Alcaeus group. The colour code shows the speed (top panel) and the distance (bottom panel) in the radial direction and both relative to the group centre. The orange ellipses show samples from the posterior distribution of the group-level parameters (see Fig. 2).

5.2.1. Core and halo populations of IC 348 and Gorgophone

In IC 348 and Gorgophone, our methodology finds two Gaussian components that we call core and halo populations, with the core having the smallest dispersion in the 3D positional space and halo the largest (see the value of σ_XYZ in Table 3). In IC 348, the medians and covariance matrices of these two Gaussians result in Mahalanobis distances between them of 10.16 (halo with respect to the core) and 1.23 (core with respect to halo). Given that these distances are mutually farther away than one Mahalanobis distance, then we conclude that they correspond to independent physical groups (see Assumption 2). However, in the case of Gorgophone, the Mahalanobis distances are 1.6 (halo with respect to the core) and 0.96 (core with respect to halo), which prevent us from concluding that they pertain to independent physical groups (see Assumption 2). For historical reasons, we continue using Herbig (1998) nomenclature of core and halo populations for the two identified physical groups of IC 348 (see Sects. 2.2 and 2.7).

5.2.2. Internal kinematics

We analyse the internal kinematics of the groups using the inferred positions and velocities of both the sources and the groups. As an example, Fig. 3 shows the 3D positions (top panel) and 3D velocities (bottom panel) of the Alcaeus group. The colour code of this figure shows the distance (bottom panel) and speed (top panel), in the radial direction, both with respect to the group’s centre. As shown in this figure, there are no observable trends of expansion. Appendix B shows figures with the Galactic Cartesian positions and velocities of the candidate members in our eight statistical groups. As can also be seen in those figures, there are no observable trends of expansion in any of the Perseus groups.

As explained at the end of Sect. 4.2, to objectively quantify the internal kinematics of the groups, we computed the average magnitude of the dot and cross products of the radial distance and velocity vectors of all the group’s members, which are proxies of the group’s expansion and rotation, respectively. Columns fourth and fifth of Table 4 show the average values of the dot and cross product vectors, respectively. Although the previous values show some trend of contraction, particularly in NGC 1333 and the core of IC 348, the current uncertainties do not allow us to draw firm conclusions. Similarly, the uncertainties in the cross product show that the observed trends of rotation are significant only at the one-sigma level, but fail to exceed the two-sigma level.

5.3. Empirical isochrones and age estimates

As described in Sect. 4.3, the Miec code delivers, as a by-product, the extinction-free empirical isochrone of the group under analysis. Figure 4 shows the empirical isochrones of the identified groups. We notice that in the case of Gorgophone, both the core and halo have the same isochrone. In addition, Fig. 5 shows, for each of the groups, the absolute CMD of the candidate members (black dots) as well as their extinction-free empirical isochrones (solid black lines). The apparent magnitudes of both the candidate members and the empirical isochrones were transformed into absolute ones using the source- and group-level distances (see Sect. 5.2), respectively. The figures also show the theoretical isochrones from the PARSEC, MIST, and BT-Settl models for the ages of 1, 3, 5, 7, and 10 Myr. We notice that due to the scarcity of candidate members in the high-luminosity region (absolute G < 7 mag), the empirical isochrones do not follow the curvature of the theoretical isochrones but that of our prior, which is a simple linear regression (see Fig. B.3 of Olivares et al. 2021).

Fig. 4. Empirical isochrones of the Perseus groups as obtained by Miec.

Fig. 5. Absolute CMD of the Perseus groups candidate members (black dots). The black and coloured lines show the empirical isochrone and the theoretical ones of different evolutionary models, respectively.

We estimate the ages of the physical groups by comparing their extinction-free empirical isochrone to the theoretical ones (see Sect. 4.3) in the faint region (G > 8 mag) where the bulk of the candidate members are located. Our age estimates are shown in the fourth column of Table 1. Given that these estimates are based on a simple visual comparison, we do not provide uncertainties. We stress the fact that these ages highly depend on the isochrone models and extinction map and thus will benefit from further refinement.

5.4. Mass distributions

We infer the mass distributions of each of the physical groups using the two methods described in Sect. 4.4. Figure 6³ shows the result of these inferences. The orange lines depict a hundred realisations of the mass distribution that result from the propagation of the same number of samples taken from the posterior distributions of the Miec parameters that describe the G band magnitude distribution. The magnitude distributions are transformed into mass distributions using the theoretical mass-magnitude relations of the unified theoretical model PMB (see Sect. 4.4) at the group’s age and distance (see Table 1). The mass distributions inferred using the Sakam code with the theoretical (i.e. PARSEC, MIST, and BT-Settl) and unified (PMB) models are shown in the same figure as coloured lines. The grey line shows the Chabrier (2005) mass prior, and the grey area depicts the incompleteness region of the Gaia data, which corresponds to G = 19 mag (see Tables 4 and 5 of Lindegren et al. 2021) and has been extinction corrected.

Fig. 6. Mass distributions of the Perseus groups. The orange lines show 100 realisations from the mass distribution obtained after transforming Miec’s magnitude distributions (see text). The rest of the coloured lines (those of the legend) depict the mass distributions computed with Sakam. The grey area shows the Gaia incompleteness region.

As can be observed in Fig. 6, the mass distributions inferred with the Sakam code agree for all the theoretical models except at the borders of their mass domains. Indeed, the peaks observed at log Mass [M_⊙] ∼ − 1 and ∼0.2 correspond to border effects introduced by the lower limits of PARSEC and MIST models and the upper one of the BT-Settl model, respectively. As shown in the figure, the unified PMB model does not show these artefacts.

We observe that the uncertainty in the mass distributions obtained with the Miec code (dispersion of the orange lines) is proportional to the population size of the group. The uncertainty of IC 348 mass distribution is the smallest, while that of NGC 1333 and Autochthe are the largest ones.

Comparing the mass distributions inferred with the two methods, we observe that the largest discrepancies are observed in IC 348 while the smallest ones are in Gorgophone and Alcaeus. The extent of these discrepancies is proportional to the extinction value of the group. As will be shown in the next section, the mode of the inferred extinction distributions is the largest in IC 348 and the lowest in Gorgophone (see Fig. 7). The discrepancy in the mass distributions is explained by the difficulties that the Miec code has to infer the magnitude distributions of extincted regions under the presence of low-information-content datasets (see the discussion of Olivares et al. 2021), in this particular case, the visual bands of Gaia. Thus, we expect that in the heavily extinct groups, the mass distributions inferred with the Sakam method are more realistic than those of the Miec one because Sakam uses additional photometric bands (see Sect. 3), in particular the infrared ones that are less affected by extinction.

Fig. 7. Extinction distributions of the Perseus groups. The colour code indicates the results of the theoretical and unified models as well as the uniform prior.

5.5. A_v and R_v distributions

The mass inference done with the Sakam code also delivers samples from the posterior distributions of the A_v and R_v values of each source. Figures 7 and 8 show histograms and kernel density estimates of the posterior samples of A_v and R_v, respectively, of all the candidate members of each physical group. As in Fig. 6, the coloured lines indicate the theoretical and unified models, and the assumed prior distribution.

Fig. 8. Rv distributions of the Perseus groups. Captions are the same as in Fig. 7. The vertical dotted line at Rv = 3.1 depicts the mode of the Gaussian prior.

The figures show that the A_v extinction is highly variable, with the mode ranging from 0.8 mag in Gorgophone to 2.3 mag in IC 348. In addition, the within-group differential extinction has a large dispersion, with the exception of the Alcaeus group, in which the internal dispersion is only ∼0.5 mag.

We notice that the PARSEC models in all physical groups, except in Gorgophone, for which we use the 7 Myr isochrone, produces lower values of A_v compared to those obtained with the MIST and BT-Settl models. This effect is a consequence of the redder theoretical isochrones of the PARSEC model compared to the other two models, particularly at the faintest magnitudes of the infrared bands (absolute K > 4). However, at 7 Myr, the isochrones of the three models completely agree, which explains the negligible difference among the inferred A_v values of the three theoretical models in the Gorgophone group.

Concerning the inferred distributions of R_v, we observe that the three theoretical models result in similar distributions in all the groups. Thus, the results are all consistent. We notice, though, that in the case of NGC 1333, the mode of the distribution is shifted from that of the prior towards a larger value of 3.4. We further discuss the consequences of this latter value in Sect. 6.5.

5.6. Dynamical analysis

Using the group level parameters inferred in previous sections, we now analyse the dynamical state of each of the Perseus physical groups. As mentioned in Sect. 4.5, we make this analysis with two methods: the observed energy distributions and the comparison of the group’s internal velocity dispersions with that expected for virial equilibrium.

5.6.1. Energy distributions

We computed the energy distributions of the Perseus groups using Eq. (2) and samples from the posterior mass, position, and velocity distributions of each candidate member. As explained in Sect. 4.5, we corrected our mass estimates by accounting for possible unresolved binaries (see Assumption 7) and the non-stellar mass of dust and gas (see Assumption 6).

Figure 9 shows the cumulative energy distribution functions (CEDFs) of the Perseus physical groups. As can be observed, all the Perseus groups have positive energies except for Alcaeus and the core of IC 348. The fraction of gravitationally bound stars varies between 10% and 25% in Alcaeus and 15% to 40% in the core of IC 348. In this latter group, the fraction of bounded stars depends highly on the applied correction for the non-stellar mass fraction (dust and gas), with the larger fraction of bound stars corresponding to the most massive gravitational potential, as expected. We also observe that this non-stellar mass correction has a negligible impact on energy distributions of the other Perseus groups. Given that our non-stellar mass correction follows the distribution of observed stars (see Assumption 6) rather than being massive particles located at the centres of the groups, then its effect is most pronounced in the most compact group, which is the core of IC 348 (see the second column of Table 4).

Fig. 9. Cumulative energy distributions of the Perseus groups. The dotted and dashed lines depict the results of accounting for the dust and gas mass by increasing the stellar mass by 65% and 169%, respectively.

5.6.2. Virial state

We now estimate the virial state of the Perseus physical groups. First, we computed the velocity dispersion, σ_vir, that the groups would have if they were in virial equilibrium (see Sect. 4.5 and Assumptions 4 and 5). Then, we compare these σ_vir with the observed velocity dispersions σ_UVW (see Sect. 5.2) and estimate the dynamical state of each group. For completeness reasons, we also compute the total mass that would be needed for the groups to be in virial equilibrium given their observed velocity dispersion.

Table 5 shows, for each Perseus group, the EFF’s core radius (r_c) and γ parameters, the half-mass radius (r_hm), the virial (σ_vir) velocity dispersions, and the mass factors (π_vir) that should multiply the observed stellar mass for the group to be in virial equilibrium (assuming a binary mass fraction of 20%). In the virial velocity dispersions and mass factors, we show the values computed using the core and half-mass radii (sub-indices EFF and HM, respectively) as well as the original, lower, and upper limits of the cluster mass (sub-indices 0%, 65%, and 169%, respectively). As expected, the most massive clusters (i.e. those with the 169% non-stellar mass correction) have larger virial velocity dispersions. Comparing the observed velocity dispersions (third column of Table 4) with the virial ones, we observe that all Perseus groups are in a super-virial state. In other words, for the groups to be in virial equilibrium, their stellar mass should be increased in factors ranging from six times to a thousand times. We further discuss these results in Sect. 6.6.

6. Discussion

We now compare our results with those from the literature and discuss their differences. We proceed in the same order as in the previous sections.

6.1. Membership

When comparing lists of candidate members is important to highlight that the parallax is the most discriminant feature that Gaia DR3 provides to identify candidate members. Thus, we use it to exclude clear outliers from previous results of the literature. We consider that sources with parallaxes lower than 1.9 mas (farther than 526 pc) and higher than 6.5 mas (closer than 178 pc) belong to the field population. We chose this highly conservative parallax interval to ensure the inclusion of possible members within more than −100 pc and +200 pc around the traditional Perseus distance (320 pc, Ortiz-León et al. 2018). Figure 10 shows the representation space coordinates of our candidate members that are common with (dots) and rejected from the literature works. The colour code shows the membership probability of the sources. The following paragraphs present detailed comparisons of our candidate members with selected works from the literature.

Fig. 10. Representation space coordinates of the confirmed (dots) and rejected members from the literature after removal of clear outliers (ϖ ∉ [1.9, 6.5] mas). The membership probability is shown as a colour code.

Luhman et al. (2016) identified 478 and 203 members in IC 348 and NGC 133, respectively. Only 450 of these members have Gaia DR3 parallaxes, proper motions, and photometry. Of these, 16 are clear outliers according to our parallax limits. Our methodology recovers 380 and rejects 54 of their candidate members. As can be seen in Fig. 10, most of the rejected sources lay on the outskirts of the group’s loci.

From the list of members of Ortiz-León et al. (2018) 162 have a counterpart in our dataset. Out of these, we recover 158 as candidate members and reject one due to its discrepant parallax. We identify the remaining three sources as false negatives of our methodology, given that their coordinates are consistent with the locus of the Perseus groups.

Sixteen of Allers & Liu (2020) members are in our dataset, and four of them are clear parallax outliers. Of the remaining sources, nine are recovered by our methodology, and three were rejected. Of the latter, only one can be identified as a false negative of our methodology since the other two lay at the outskirts of the groups’ loci.

Pavlidou et al. (2021) identified 913 candidate members in their five groups, which do not include IC 348 or NGC 1333. Contrary to the previous authors, our main objective is the analysis of these latter two clusters, and we focus on a smaller region (see Sect. 3). This region contains 183 of their candidate members, and 2 of them are clear parallax outliers. Our methodology recovers 166 of their candidate members and rejects 15. Half of the rejected candidates lay slightly below our probability thresholds (blue and green diamonds) and can be considered false negatives by our methodology.

Kounkel et al. (2022) identified 810 candidate members in the Perseus region, but only 429 are in our sky region and 425 in our dataset⁴. From the latter, we recover 416 and reject nine. Only four of the rejected have membership probabilities above 0.5, which can be considered false negatives of our methodology.

Pang et al. (2022) found 211 candidate members to IC 348, 195 of these are within our sky region, and we recover 192 of them. The three rejected sources have probabilities larger than 0.5 and thus can be considered false negatives of our methodology.

Wang et al. (2022) found 207 candidate members, out of which 197 are in our sky region, and we recover 191 of them. Out of the six rejected, two have negligible membership probabilities, while the remaining four have membership probabilities > 0.5 but lower than our probability threshold. Unfortunately, these authors do not provide a classification of their candidate members into the Perseus groups. Therefore, we cannot make more detailed comparisons with our candidate members to the Perseus groups.

As can be observed in Fig. 10, there are sources with very low membership probabilities that lay at the outskirts of the cluster locus and that we consider field population. Nonetheless, there are 12 sources from Luhman et al. (2016), three of Ortiz-León et al. (2018), one of Allers & Liu (2020), five of Pavlidou et al. (2021), four of Kounkel et al. (2022) and three of Pang et al. (2022) with coordinates consistent with the cluster locus that lay below our probability thresholds, we consider these members as false negatives of our methodology. Under the conservative assumption that all the rejected members from the literature (excluding clear outliers with discrepant parallaxes, as described above) are indeed true Perseus members, then the true positive rate of our methodology is 93.2%, which is better than the 88% estimated by Olivares et al. (2021, see their Table C.3) for a cluster with A_v ∼ 6 mag. The better performance of the current application of the Miec code results from the use of the highly discriminant parallax feature and the exquisite precision of the Gaia DR3 astrometry. Finally, we highlight that our membership methodology recovers 267 candidate members not previously identified as such by the literature. This represents an increase of 31% with respect to the previous studies of the Perseus region.

6.1.1. Group’s classification

We now compare the list of members that we recover for IC 348 (core and halo) and NGC 1333 with those found in other works. Tables 6 and 7 show, for IC 348 and NGC 1333, respectively, the number of members from other works that are in our dataset and that we recover and reject. The third column of these tables shows the number of recovered members together with those sources that we still recover as members of the Perseus region but belonging to another group (shown with the first two letters of its name).

In addition to the well-known IC 348 and NGC 1333 groups, the Perseus region contains several substructures that have been recently identified thanks to the high-precision Gaia data. For example, Pavlidou et al. (2021) identified five Perseus groups, named Alcaeus, Autochthe, Electryon, Heleus, and Mestor. In the rest of this subsection, we discuss the recovery of the members of these groups, with the exception of the Mestor one, whose members are outside the sky region that we analyse.

Concerning the Electryon group, out of the 329 members of Pavlidou et al. (2021), only 21 are in our dataset. Our methodology recovers 15 of these, although they are classified as members of the Heleus group. The fact that our methodology classifies these Electryon members into the Heleus group is explained by the small fraction of the former group within our dataset (6%), its negligible contribution to the Perseus total number of members (< 1%), and the entanglement that these two group show in the parallax-proper-motions space. We notice that Pavlidou et al. (2021) were able to disentangle Heleus from Electryon due to its relative compactness in the sky coordinates; however, our methodology identifies members of this group spread over a larger sky region.

Concerning the Heleus group, our methodology recovers the 29 members of Pavlidou et al. (2021) present in our dataset, and we identify 95 more candidate members. With respect to Kounkel et al. (2022), out of their 66 candidate members, our dataset contains 22, of which we recover 21, all of them belonging to this group. In addition, we identify 103 more candidate members than these authors.

With respect to Alcaeus, from the 170 members of Pavlidou et al. (2021), 108 are in our dataset, and our methodology recovers 98 of these, with 88 classified into Alcaeus and 10 in Gorgophone. Kounkel et al. (2022) identified 124 candidate members, 83 of which are in our dataset. We recover 81 of these, although 20 of them are classified into Gorgophone. In addition, our method discovers 64 new candidate members of this group.

Finally, from the 27 members of Autochthe identified by Pavlidou et al. (2021), 25 are in our dataset, and our method recovers 24 of them, with two of them classified into Gorgophone. Kounkel et al. (2022) found 25 candidate members, with 23 of these in our dataset. Our methodology recovers 22 of these, although with one belonging to NGC 1333 and 2 to Gorgophone.

As shown in the previous paragraphs, further disentanglement of the Perseus substructures remains a difficult task in which more precise data, particularly radial velocities, and further methodological developments will still be needed.

6.2. Phase-space structure

We now compare the distribution of positions and velocities of the identified physical groups with those from the literature. The most striking difference with respect to previous pre-Gaia works is the number of identified substructures. Classically, the region comprised IC 348 and NGC 1333; however, the arrival of Gaia unravelled multiple populations. Compared to the latest analyses from the literature (Wang et al. 2022; Kounkel et al. 2022; Pavlidou et al. 2021), we identified one additional physical group that we call Gorgophone, which is composed of a core and halo populations.

In IC 348, our kinematic criterion to identify physical groups (see Assumption 2) allow us to conclude that the core and halo of IC 348 are distinct physical groups. The existence of these two populations is further supported by the age and mass distribution features of these populations, with the halo population being older and having fewer low-mass stars and more high-mass stars than the core. In addition, the most massive stars of the halo population are located off-centre while those in the core are centrally concentrated. These findings support a distinct halo population that formed earlier than the core population and possibly quenched the formation of high-mass stars in this latter.

6.2.1. 3D distribution

Our results show that the core of IC 348 has a (1σ) radius of 0.65 ± 0.20 pc, with small correlations (< 0.3) among the X, Y, and Z coordinates. These negligible correlations are in agreement with previous works that found a centrally concentrated distribution (see Sect. 2.2). On the contrary, the halo extends over a (1σ) radius of 6.89 ± 0.75 pc, with non-negligible correlations of 0.3 between X and Y, 0.5 between X and Z, and 0.7 between Y and Z.

Concerning IC 348, our values of the halo radius are larger than the 0.9−1.4 pc radius reported by Scholz et al. (1999), Carpenter (2000), and Muench et al. (2003). On the other hand, our 3D core radius estimate is in agreement with the 2D radius of 0.5 pc reported by Muench et al. (2003). These latter authors also warned about the possibilities of a larger cluster radius and possible entanglement between the halo sources and the core ones due to projection effects. Evidence of a larger cluster size was also found by Cottaar et al. (2015), who identified 63 of their candidate members beyond a radius of 1.8 pc radius (20′). Our results confirm the following suggestions proposed by Muench et al. (2003): a larger cluster size extending into larger areas and elongated in the north-south direction, halo sources entangled within the core radius (see Fig. 1), and a density profile that is far from trivial. On the other hand, we do not find evidence of substructures within the halo populations, as suggested by the aforementioned authors.

Contrary to IC 348, NGC 1333 is well described by a single Gaussian. It shows an elongation twice as large in the X coordinate that in the Y and Z, directions, although with negligible correlations and unrelated to the line-of-sight elongation. Given the cluster’s young age and the low number of stars (see Table 1), the observed elongation can be primordial.

Concerning distances to IC 348 and NGC 1333, our estimates are in agreement with those of Ortiz-León et al. (2018) and Zucker et al. (2018). The surprising agreement between the early distance estimate of 316 ± 22 pc made by Strom et al. (1974) and our 315 ± 1 pc Gaia value for the core of IC 348 further supports the hypothesis proposed by those authors about dust grain growth in the dense Perseus clouds (see Sects. 2.6 and 6.5).

We notice that Ortiz-León et al. (2018) warned about the impossibility of disentangling the cluster structure along the line of sight due to the fact that the parallax dispersion was dominated by individual uncertainties. However, this is no longer the case, thanks to the high precision of the Gaia DR3 and the uncertainty deconvolution applied by Kalkayotl.

6.2.2. Velocity distribution

Our results on the velocity distributions of the Perseus groups are shown in Tables 3 and 4. As can be observed, the total velocity dispersion (third column of Table 4) of these groups varies from 0.76 km s⁻¹ in Alcaeus to 3.79 km s⁻¹ in the halo of Gorgophone.

Comparing the halo velocity dispersions of IC 348 and Gorgophone with respect to their cores, we observe that in IC 348, these values are similar (and compatible within 2σ uncertainties), whereas, in Gorgophone, the halo velocity dispersion is significantly more than twice that of its core. The relatively large velocity dispersion and size of Gorgophone’s halo to that of IC 348, in combination with their similar ages, suggests that different mechanisms have formed these halo populations. We hypothesise that the core and halo of IC 348 formed out of the same molecular cloud but at different ages. Thus, although physically independent, they inherit similar velocity dispersions as those of the parent molecular cloud. On the contrary, the core and halo of Gorgophone pertain to the same physical group, but the halo-core dichotomy formed as a result of the dynamical interactions. Thus, the halo is dynamically hotter, and the core is relatively depleted of low-mass stars (see Fig. 6).

Cottaar et al. (2015) fitted the radial velocity profile of IC 348 with a GMM with two components (see Sect. 2.3). They argue that the presence of the second component may have three possible origins: contaminants, a halo population, or lack of relaxation into a Gaussian distribution due to the young age of the cluster. Unfortunately, given the 1.8 pc radius cut that these authors applied to their sample, their results cannot be directly compared with those of the IC 348 core that we found here. Indeed, at that radius, we find a mixture of both core and halo populations. Our results confirm the hypothesis of the halo population, which we found in the 6D space of positions and velocities rather than only in the 1D radial velocity, as those authors hypothesised. On the other hand, the contamination scenario can be ruled out thanks to our robust membership methodology, in which the expected contamination rate of this group is ≲7%.

In the case of NGC 1333, Foster et al. (2015) used the same methodology as Cottaar et al. (2015) and found that the radial velocity distribution can be well described by a single Gaussian component. These results agree with our single Gaussian model for the phase-space structure of NGC 1333.

6.2.3. Internal kinematics

Concerning IC 348, in its core, we find a 3D velocity dispersion of 1.44 ± 0.08 km s⁻¹, which is similar than the 3D 1.2 ± 0.1 km s⁻¹ value (0.72 ± 0.07 km s⁻¹ for the 1D radial velocity dispersion) found by Cottaar et al. (2015), and smaller to the 2.28 km s⁻¹ reported by Ortiz-León et al. (2018). On the one hand, the smaller value found by Cottaar et al. (2015) can be explained by their better treatment of binaries, which in our analysis are not distinguished from single stars. On the other hand, the larger value found by Ortiz-León et al. (2018) can be explained by their use of a single Gaussian to model stars in both the core and halo populations. In NGC 1333, the 3D velocity dispersion of 1.6 ± 0.2 km s⁻¹ (0.92 ± 0.12 km s⁻¹ for the 1D radial velocity dispersion) and 2.0 km s⁻¹ reported by Foster et al. (2015) and Ortiz-León et al. (2018), respectively, are smaller than and compatible with our 2.37 ± 0.23 km s⁻¹ value, respectively (see Table 4). We notice that the 3D velocity ellipsoid of NGC 1333 is 33% wider in the Y direction than in the Z direction, which makes it far from spherical, as assumed by Foster et al. (2015). On the other hand, the velocity dispersions measured by us and by Ortiz-León et al. (2018) are also probably inflated by unresolved binaries. Future work is warranted upon the arrival of more precise radial velocity measurements.

With respect to expansion and contraction, our results show that the observed precision is not enough to claim any clear trend. Nonetheless, we notice that the largest values of contraction are those of the core of IC 348, Alcaeus, and NGC 1333. Concerning the signature of internal rotation, our results indicate that it is significant only at the one-sigma level. In IC 348, Cottaar et al. (2015) have already reported a small but significant radial velocity rotation gradient in the plane of the sky: 0.024 ± 0.013 km s⁻¹ arcmin⁻¹ or 0.26 ± 0.14 km s⁻¹ pc⁻¹. Scaling their 1D gradient to 3D values at the cluster core and halo radii, we obtain 0.3 ± 0.1 km s⁻¹ and 3.1 ± 0.1 km s⁻¹, respectively. Although the rotation values we find for the core and halo of IC 348 are larger and smaller, respectively, than those reported by Cottaar et al. (2015), they still are not significant enough to rule out the null hypothesis of no internal rotation. Moreover, Ortiz-León et al. (2018) also found non-zero values for the rotation vectors of IC 348 and NGC 1333. However, they concluded that their values were negligible when compared to the internal velocity dispersion. Similarly, Kuhn et al. (2019) found non-zero azimuthal velocities in NGC 1333 and IC 348, but they also judged these as non-significant when compared to the measured uncertainties. In the light of the previous works and current evidence, future efforts are still needed, both on the observational and modelling sides, to test the hypothesis of internal rotation in the Perseus groups.

6.2.4. The core and halo of IC 348

As mentioned in Sect. 2.2 the core and halo populations were first found by Herbig (1998) on the base of the spatial distribution of stars with Hα emission. He also found an age gradient of 2.4 Myr pc⁻¹ in which older ages are located at larger radii (1.45 Myr at 4′ and 2.8 Myr at 10′). Here, we confirm the existence of these two IC 348 populations and observe a smaller but still non-negligible age gradient of 0.32 Myr pc⁻¹ computed on the base of the ages and radii of the core and halo.

We notice that Cottaar et al. (2015) found a correlation between the radial velocity and the reddening in IC 348, for which they proposed three different scenarios: a systemic offset in the radial velocities, the contraction of the cluster as a whole, or the convergence of two sub-clusters aligned in the line of sight. Their third scenario naturally coincides with the core and halo populations proposed by Herbig (1998) and confirmed here. These two populations naturally solve the previous issue because the halo population is red-shifted 1.6 km s⁻¹ with respect to the core population and has a slightly larger spread in extinction (see Fig. 7).

Furthermore, we observe that the halo contains slightly more massive stars than the core (see Figs. 5 and 6). This effect was already predicted by Herbig (1998) when he mentioned that the brightest stars from the halo formed first, and then the molecular cloud retreated behind them, where now the youngest population is formed.

6.3. Empirical isochrones and age estimates

Our age estimates for NGC 1333 and IC 348 are compatible with the literature values. In the case of IC 348, our 3 Myr estimate for the core is slightly larger than the 2 Myr median value of Luhman et al. (2003). Similarly, the 5 Myr age of IC 348 halo is similar to the 6 Myr reported by Bell et al. (2015), and in agreement with the recent 5 Myr age determination of Pang et al. (2022) and Lalchand et al. (2022). Both our core and halo age estimates fall within the age interval reported by Luhman et al. (2003). In the case of NGC 133, our age estimate of 3 Myr is older but still compatible with the 1−2 Myr derived by Lada et al. (1996). Furthermore, our age estimates for the halo of IC 348 and NGC 1333 are in clear agreement with those of Wang et al. (2022, see our Sect. 2.4), who determined isochrone ages of 5.4 Myr and 2.9 Myr for IC 348 and NGC 1333, respectively. We notice that this agreement is based not only on the theoretical isochrone models but on the similar extinction distribution, as observed when comparing their Fig. 7 to our Fig. 7.

Comparing the ages estimated by Pavlidou et al. (2021) for the newly identified Perseus groups with those obtained here, we observe that the relative internal ages are in agreement. Furthermore, we confirm the coevality of NGC 1333 and Authochte. However, the absolute ages differ by a factor of roughly two, with our estimates being twice older than theirs. In the particular case of the absolute age estimates obtained by Pavlidou et al. (2021), we conclude that their underestimation is a direct consequence of their lack of extinction correction. According to the extinction maps of Green et al. (2019), the members of Perseus groups have extinctions as large as 7 mag in A_v, with a median value of A_v = 2.6 ± 1.9 mag. Neglecting these extinction values when estimating isochronal ages results in younger age determinations. The same reasoning applies to the more recent isochrone age determinations by Kounkel et al. (2022).

6.4. Mass distributions

As mentioned in Sect. 5.4, the mass distributions obtained with our two methods are compatible within the uncertainties, except in those groups with the largest peaks in extinction, which are the core and halo of IC 348, NGC 1333, and Autochthe (see Fig. 7). The differences in the inferred mass distributions arise from the difficulties that the Miec method has to recover luminosity distribution of extincted regions based on datasets lacking infrared bands. The Sakam method complements Gaia data with infrared photometry and thus is not affected by this problem.

Comparing Sakam mass distributions of all groups within Gaia completeness limits, we observe that there is a general agreement with Chabrier (2005) initial mass distribution. This comes as no surprise given the young ages of the groups. Nonetheless, we notice a depletion of stars more massive than ∼5 M_⊙ in Authochte, NGC 1333, Heleus and the cores of IC 348 and Gorgophone compared to the values expected from Chabrier (2005) distribution. Unexpectedly, the halo populations of IC 348 and Gorgophone seem to have more massive stars than the core populations. In almost all groups, there is a visible peak at masses ranging from 1.5 M_⊙ to 2.5 M_⊙, which results from abrupt changes in the slope of the mass-luminosity relation of the PARSEC and MIST models at this particular mass interval (see Sect. 4.4).

Figure 11 shows the cumulative mass distribution and incompleteness regions (shaded areas) of the core and halo populations of IC 348 as obtained from the unified isochrone model (PMB; see Sect. 4.4). As can be observed in the figure, the core and halo populations are not identical. This observation is confirmed by a two-sample KS test. We applied this test to one hundred samples with bootstrap realisations of the core and halo stars. The mean value of these hundred KS tests results in the rejection (p-value < 1%) of the null hypothesis that both populations come from the same parent distribution.

Fig. 11. Cumulative mass distribution of IC 348 and NGC 1333.

We observe that the halo region of IC 348 has more massive stars than the core region despite its lower number density. Furthermore, according to the unified PMB model, the core mass distribution peaks at 0.12 ± 0.03 M_⊙, whereas the halo peaks at 0.25 ± 0.07 M_⊙. This result is in disagreement with the observations of Muench et al. (2003), who found that the core and halo mass distribution peak at 0.56 ± 0.18 M_⊙ and 0.1 ± 0.02 M_⊙, respectively. As can be observed in the top right panel of Fig. 6, our halo mass distribution starts to decrease at 0.2 M_⊙, within the Gaia completeness limits (> 0.11 M_⊙). Muench et al. (2003) observed that the alignment of the peaks of their core and halo mass distributions would require a halo age between 5 and 10 Myr. Our age determination for the halo population, 5 Myr, agrees with the lower limit proposed by those authors. However, they mentioned that the hypothesis of two distinct populations with an older halo was difficult to accept. They base their conclusion on the following argument. If the two-population model were correct, then the similar sizes of their core and halo regions will imply a truncation of the halo population below 0.3 M_⊙, which was clearly not observed in their mass distributions. Here, we observe that the number of halo stars is half that of the core and that the halo mass distribution peaks at 0.25 ± 0.07 M_⊙ with a clear cut at 0.2 M_⊙. Based on our previous results, we have no compelling evidence to reject Herbig’s (Herbig 1998) hypothesis of distinct core and halo populations. Furthermore, the kinematic evidence further supports the hypothesis of core and halo populations.

Comparing the cumulative distribution of NGC 1333 (also shown in Fig. 11) with the core and halo populations of IC 348, we observe the following. On the one hand, comparing the mass distributions of the halo of IC 348 and that of NGC 1333, we observe that NGC 1333 has a general overabundance of stars less massive than 2 M_⊙. This overabundance of low-mass stars and brown dwarfs was already reported by Scholz et al. (2013), who observed it in several of their mass scenarios (combinations of age, distance and extinction). On the other hand, comparing the mass distributions of NGC 1333 and the core of IC 348, we observe that although in the mass interval of 0.3 − 1.0 M_⊙ the core of IC 348 appears to have slightly more stars than NGC 1333, the two distributions are remarkably similar with a KS p-value of 0.53, which prevents us from rejecting the null hypothesis that the two distributions are random realisation of the same underlying distribution. However, due to the Gaia completeness limits, we can only claim the previous findings in the domain of masses > 0.16 M_⊙. Therefore, we conclude that, within the completeness limits of the Gaia data, we observe no difference in the mass distributions of the core of IC 348 and NGC 1333, and thus, reject Scholz et al. (2013) hypothesis of variation in the initial mass distribution with respect to the crowdedness of the environment.

6.5. Extinction and R_v distributions

The A_v and R_v values that we derive here confirm previous findings from the literature. We observe cloud-to-cloud and intra-cloud extinction variations, evidence of which was already reported by Choi et al. (2016) for the dust emissivity spectral index. Furthermore, in the extinction interval of the Perseus group analysed here, we confirm the correlations reported by Foster et al. (2013) in the B5 and West-End clouds (see Fig. 9 of those authors). As can be observed from Fig. 12, there are no significant correlations between R_v and A_v for A_v < 2 mag in any of the groups. However, for A_v > 3 mag these correlations are clearly observed in NGC 1333 and Authochte, mildly in Gorgophone, almost perceptible in IC 348, and clearly absent in Alcaeus and Heleus. Interpreting these correlations as evidence of dust grain growth, and given that Alcaeus and Heleus are off the clouds (see Fig. 1), we conclude that there is an east-west gradient in the growth of the dust size across the stars of Perseus ridge.

Fig. 12. Av and Rv 2D histograms for stars in the different Perseus groups.

Comparing the A_v values inferred by Sakam with those reported by Bayestar19 (Green et al. 2019; see our Fig. 13), we observe that they are, in all cases, compatible (i.e. the zero value is covered by the 95% credible interval of the distribution). However, NGC 1333 and the cores of IC 348 and Gorgophone show the largest discrepancies (with modes at −2.2, −1.16, and −0.91 mag, respectively), while Heleus, Alcaeus, and the halos of IC 348 and Gorgophone display the smallest differences (with modes at 0.35, 0.26, 0.02, and −0.06 mag, respectively). The largest discrepancies observed in NGC 1333 and the cores of IC 348 and Gorgophone can be explained by the lack of stars in the 2MASS and Pan-STARRS catalogues of these embedded regions that lay in the Perseus ridge (see Fig. 1). The lack of stars for anchoring the extinction values causes the Bayestar19 algorithm to overestimate them.

Fig. 13. Distributions of Av difference between the values inferred with Sakam and those reported by Bayestar19 (Green et al. 2019). Captions as in Fig. 7.

6.6. Dynamical analysis

In Sect. 5.6, we presented the CEDFs of each Perseus physical group, and a comparison of its observed velocity dispersion with that expected for virial equilibrium under varying values of the mass and radius parameters. All these results showed that the majority of the stars in the Perseus groups are gravitationally unbound and in a super-virial state. Even in the core of IC 348, which is the most massive (146 ± 5 M_⊙) and compact (0.65 ± 0.20 pc) group, the upper limit in the fraction of gravitationally bound stars reaches only 40%. The velocity dispersions at virial equilibrium computed from varying radii and total masses are all similar, although with the more massive clusters assumption resulting in larger virial velocity dispersion (see Sect. 5.6). In the case of the core of IC 348, the observed velocity dispersion is 1.6 times larger than the most conservative estimate for the virial velocity dispersion (0.89 ± 1.71 km s⁻¹ computed from the core radius and a non-stellar mass correction of 169%), whereas in the rest of the groups, the difference is even larger. For the core of IC 348 to be in virial equilibrium, a total mass six times larger than the stellar mass will be needed. This value for the dust-and-gas mass exceeds by more than three times our most conservative estimate for this non-stellar mass. Moreover, recent studies have shown that the molecular cloud seen in projection towards IC 348 is at 280−300 pc (Zucker et al. 2022) and thus in front of the cluster. This observation then implies that the core of IC 348 has an even larger super-virial state than our conservative estimates. Thus, we conclude that it is highly unlikely that the young Perseus groups formed in virial state.

Our independent estimates of the dynamical state of the Perseus groups confirm Cottaar et al. (2015) findings about the super-virial state of IC 348 and extend it to the rest of the groups. We interpret the super-virial states and large fractions of gravitationally unbound stars in the young (≲10 Myr) Perseus groups as evidence that these groups formed through a hierarchical star formation mechanism (e.g. Kruijssen 2012) rather than a monolithic one (e.g. Lada & Lada 1991). More details about the differences between these star formation scenarios are available in Ward & Kruijssen (2018) and the references therein.

Alcaeus is an interesting physical group with 10−30% of its stars in an energetically bounded state with a steep CEDF (see Fig. 9). It also has tight velocity dispersion (0.76 ± 0.21 km s⁻¹; see Table 4) similar to that of the well-known Pleiades (0.8 km s⁻¹Galli et al. 2017) and Coma Berenices (0.89 ± 0.10 km s⁻¹Olivares et al. 2022) open clusters. Therefore, we conclude that this group is a young open cluster in the process of disruption. We notice that its kinematic indicators show hints of contraction. Although this contraction is not significant, if confirmed, it may be explained by either unidentified substructures or contaminants. Future work is warranted.

6.7. Star formation history

In this section we use our results to discuss some of the star formation histories proposed in the literature. Afterwards, we propose our own.

The results that we have collected in this work can be categorised into three major topics: identified physical groups, isochrone ages, and dynamical states. The similar, although not exact same, values of the ages, distances, and kinematics of the Perseus groups suggest that it is unlikely that all these groups randomly formed out of the same star formation event with this display of properties. The observed spatial (see Fig. 14), temporal (see Fig. 4), energetic (see Fig. 9), and mass distribution (see Fig. 6) gradients all indicate that a latent (hidden) process links the star formation history of the Perseus group. Whether this latent process was a supernova explosion, the feedback from nearby OB associations, triggering from within itself, or another unknown phenomenon remains to be explained.

Fig. 14. Galactic Cartesian positions and velocities of the Perseus groups (cones) together with the positions of individual members (dots). The interactive version of this figure will be available at http://www.project-dance.com/.

We notice that our results are insufficient to link the star formation history of the Perseus groups to that of other nearby star-forming regions, as suggested in the literature. Thus, we cannot confirm nor reject the star formation histories proposed by de Zeeuw et al. (1999), Bialy et al. (2021), Zucker et al. (2022), Wang et al. (2022), or Lalchand et al. (2022). On the contrary, the star formation histories proposed by Herbig (1998), Belikov et al. (2002a), Pavlidou et al. (2021), Kerr et al. (2021), and Kounkel et al. (2022) are related, up to a certain point, to the intrinsic properties of the region and thus can be scrutinised under the light of the evidence gathered here. We now discuss these works.

Kounkel et al. (2022) discuss two scenarios for the formation of the Perseus regions that are related to the internal properties of the groups (a third scenario associated with the Per-Tau shell is also discussed): a past supernova and cloud-cloud interaction. Those authors concluded that although the region shows evidence of expansion, it is unlikely that it resulted from the triggering of a supernova, given that the observed velocities are not consistent with it. Our results show (see Table 2) that the velocities of the Per OB2a (i.e. the core and halo of IC 348 and Heleus) and Per OB2b (i.e. Alcaeus, Authochte, NGC 1333, and Gorgophone) super-groups have distinct directions but similar magnitudes, with the maximum relative difference amongst the velocities of all the groups being only 4 km s⁻¹. Concerning the cloud-cloud interaction scenario, the aforementioned authors argue that the super-groups Per OB2a and Per OB2b may have originated from two distinct clouds that independently started forming stars when they were in close proximity, a fact that is compatible with the spatial and kinematics derived here. Those authors continue mentioning that after the initial burst, the two super-groups continue forming stars resulting in today’s similar ages for all the groups. Although our age estimates are systematically older than those of the previous authors (see Sects. 5.3 and 6.3), we observe clear age gradients in the groups of Per OBa and Per OB2b. Finally, after highlighting the difficulty of measuring the degree of possible mutual triggering (e.g. Dale et al. 2015) between these two super-groups, they stress that there appears to be at least some degree of mutual influence.

Kerr et al. (2021) suggested that given the kinematic similarities between Per 1 and Per 2 (which roughly correspond to Per OB2a and Per OB2b of Kounkel et al. 2022), these super-groups most likely formed in the same star-forming process (see Sect. 2.7). The spatial and kinematic results that we gather here support this conclusion. We consider it unlikely that two unrelated star-forming events produced the observed similar spatial, kinematic, and age properties of the Perseus region. Then, the aforementioned authors argue that given the lack of age gradient and the considerable time lag between the older (Per B) and younger populations (Per A) of the two super-groups (Per 1 or Per 2): ∼17 Myr versus ∼5 Myr, then they consider unlikely that a continuous star formation process was at work. Instead, they proposed a scenario in which the formation of the older generation (Per B) dispersed part of the gas but not entirely prevent the formation of the younger generation (Per A). Our results contradict this scenario because we observe age gradients in both super-groups. On the one hand, in Per 2, we have one intermediate age generation with Heleus, and the halo of IC 348, both formed at 5 Myr, and a younger generation constituted by the core of IC 348, which formed 3 Myr ago. On the other hand, in Per 1, there are three generations with a clear age gradient: Alcaeus at 10 Myr, Gorgophone at 7 Myr, and NGC 1333 and Autochthe both at 3 Myr.

Pavlidou et al. (2021) proposed that the older groups (Alcaeus and Heleus) formed closer to the Galactic plane while the younger ones did it at higher latitudes. They also point out that NGC 1333 and Autochthe are part of the same star formation event. Both these conclusions are supported by our spatial and age results. These authors also mentioned that only NGC 1333, Authochte, and IC 348 show evidence of continuing star formation. For this latter conclusion, we have no further evidence.

Belikov et al. (2002a) considered that the Perseus region constitutes an example of propagated star formation, which started in the edge of Auriga 30 Myr ago, continued to Per OB2a 10 Myr ago, and is now in progress in the southern region where IC 348 is located. Our results support these two latter steps, with Alcaeus at 10 Myr and the core of IC 348, NGC 1333 and Autochthe at 3 Myr.

Finally, Herbig (1998) proposed a four steps scenario (see Sect. 2.7), which is fully supported by the results we gather here. We confirm the following points: (a) continuous star formation in the region for at least the last 10 Myr, with Autochthe, NGC 1333, and the core of IC 348 being the youngest examples; (b) ζ Per and o Per are (4σ astrometric) members of Alcaeus, a group formed 10 Myr ago; (c) there are low-luminosity members of Alcaeus spread over a large region; (d) IC 348 is composed of two physical groups: an old 5 Myr halo and a young 3 Myr core.

We now propose a star formation history that closely follows that outlined by Herbig (1998). We highlight the fact that our history is descriptive rather than phenomenological, and it is based on the ages and kinematics of the identified physical groups. To arrive at a phenomenological description of the Perseus star formation history, we still require precise age determinations and the application of the methods similar to the ones used here to a wider sky region encompassing the other actors in the scene: the Taurus star-forming region and its connecting bridge with Perseus. Future steps will be taken to address these points.

6.7.1. The first generation

Alcaeus and Gorgophone constitute the first generation of stars in the Perseus region. Alcaeus formed approximately 10 Myr ago, followed by the core and halo of Gorgophone 7 Myr ago. The brightest members of Alcaeus, ζ Per and o Per are B-type stars with ages of 12.6 ± 1.5 Myr and 11.1 ± 0.5 Myr, respectively, according to Tetzlaff et al. (2011). These ages are compatible with our 10 Myr age estimate for the group. Most likely, the winds of this first generation sweep up part of the gas and dust from the molecular cloud resulting in the observed lowest extinction values of about A_v ∼ 1 mag (see Fig. 7). The birth of this generation is an example of high-mass formation in a clustered environment (i.e. Alcaeus) together with low-mass formation in a sparse one (i.e. the core and halo of Gorgophone), which has also been observed in other star-forming regions (e.g. Kuhn et al. 2022; Lim et al. 2022) and is predicted by the hierarchical star formation scenario (Ward & Kruijssen 2018).

We notice that ζ Per has been related to the Per OB2 association as a runaway member (e.g. de Zeeuw et al. 1999; Mirzoyan et al. 1999), a member with discrepant radial velocity (Steenbrugge et al. 2003) or a member of the trapezium type system ADS 2843 (Allen et al. 2004, 2018), which is located at the central region of Alcaeus.

6.7.2. The second generation

Heleus and the halo of IC 348 constitute the second generation of Perseus stars. Although these two groups are separated by about ∼54 pc, they have a similar age and direction of their velocity vector (see Fig. 14 and Table 2), which indicate that probably these two groups were born from the same star formation event. The origin of Heleus may also be related to the ∼5 Myr, elongated (∼100 pc) but more distant (∼400 pc) Barnard 5 group identified by Lalchand et al. (2022). However, we currently lack evidence to relate the origin of these two groups.

6.7.3. The third generation

Finally, Autochthe, NGC 1333, and the core of IC 348 make the third and most recent generation of stars. While NGC 1333 and IC 348’s core have some hints that are still contracting (see Table 4) the value of this contraction is not significant enough to claim this effect. NGC 1333 and Autochthe have similar space velocities (see Fig. 14 and Table 2) and ages. However, the space velocity of IC 348’s core is different from that of the latter two groups (see Fig. 14 and Table 2). Thus, we conclude that the core of IC 348 formed in a coeval but rather different, possibly unrelated, or parallel (as suggested by Kerr et al. 2021) star formation episode than that giving birth to NGC 1333 and Autochthe, with these latter two most likely representing the youngest star formation episode.

6.7.4. Internal triggering?

Internal triggering, or self-propagated star formation, has been proposed as a mechanism to explain the origin of the Perseus groups (e.g. Belikov et al. 2002a). Although a clear demonstration of triggered star formation is extremely complex and elusive (e.g. Dale et al. 2015), the Perseus region gathers the following common indicators of internal triggering. Evidence of heating from B stars, embedded protostars, and outflows was found by Choi et al. (2016). Outflows have also been found by Arce et al. (2010) and Plunkett et al. (2013). These outflows affect the regulation of the star formation process, given that they can intensify or diminish it, particularly in the environment surrounding the active star formation regions. Our results show (see Fig. 14) that the super-group Per OB2b (i.e. Alcaeus, Gorgophone, NGC 1333 and Autochthe) follows a clear spatio-temporal vector gradient that points away from Alcaeus (10 Myr) and towards NGC 1333 and Autochthe (3 Myr), passing by Gorgophone (7 Myr). Per OB2a (Heleus and the core and halo of IC 348) also follows a clear age gradient but lacks the neat spatial relation observed in Per OB2b. Although the previous evidence is not conclusive enough to claim internal trigger in Per OB2a, Per OB2b, or the entire Perseus region, it nonetheless represents a step forward in understanding the star formation history of this region.

7. Conclusions

We have applied Bayesian methodologies to the public Gaia, APOGEE, 2MASS, Pan-STARRS, and SIMBAD catalogues to identify and characterise the physical groups of the Perseus star-forming region. We find 1052 candidate members distributed into seven physically distinct populations: the core and halo of IC 348, Autochthe, Alcaeus, Gorgophone, Heleus, and NGC 1333. Gorgophone is a new kinematic group composed of a core and a halo. Our new list of candidate members provides 31% more with respect to those from the literature in the Perseus region analysed here.

7.1. Most important conclusions

• Kinematics. We propose that Alcaeus is an open cluster in a state of disruption. On the other hand, the internal velocities of NGC 1333 and the core of IC 348 show the largest, although not significant, signals of contraction. In all the Perseus groups, the evidence of internal rotation is not conclusive.

• Age. Our age estimates are compatible with those from the literature. We notice that neglecting or underestimating the contribution of the extinction in this partly embedded star-forming region may result in twice lower age estimates than those reported here (see, for example, Pavlidou et al. 2021; Kounkel et al. 2022).

• Mass. The mass distribution of NGC 1333 shows no over-abundance of low-mass stars compared to that of the core of IC 348. This contradicts previous claims about the environmental differences of the initial mass function (e.g. Scholz et al. 2013). On the contrary, the mass distributions of the young Perseus groups are broadly compatible with Chabrier (2005), although with distinct features.

• Extinction. We find evidence of dust grain growth in NGC 1333 and Autochthe, as previously reported by Chen et al. (2016) for the Perseus clouds.

• Dynamical state. All of the Perseus groups are in a super-virial state, with CEDFs showing large fractions of energetically unbound stars. These findings, together with the lack of clear expansion patterns, support the hierarchical star formation scenario (e.g. Ward & Kruijssen 2018) over the monolithic one.

• Star formation history. The Perseus region contains stars from at least three generations, which supports the star formation scenario proposed by Herbig (1998).

7.2. Caveats and future perspectives

Thanks to the unprecedented quality and abundance of Gaia data in combination with our comprehensive Bayesian methodologies, the traditional uncertainty sources in the mass determinations, such as distance, extinction, completeness, and membership status, are now minimised. However, the following issues remain.

Age continues to be the largest source of uncertainty for the analysis of the initial mass function and star formation history of the Perseus region. Precise age determinations, such as those provided by dynamical trace-back analyses (Miret-Roig et al. 2020, 2022b), are the next mandatory step before attempting to construct a phenomenological star formation history of the region.

A phenomenological model for the star formation history of the Perseus region also demands a joint analysis of its neighbour structures (i.e. Taurus, Auriga, and Pleiades). This analysis will be the next step to linking the origin of the Perseus region to large-scale structures such as the Per-Tau shell. However, the disentanglement of physically coherent structures in large sky regions is a methodological challenge for which our team has made only the first steps (e.g. Olivares et al. 2022).

The incompleteness limits of the Gaia data still prevent sound comparisons of the brown dwarf ratios in the Perseus groups, particularly between IC 348 and NGC 1333. Deep, wide, and complete astrophotometric surveys, such as those provided by the COSMIC-DANCe project (e.g. Bouy et al. 2013, 2015; Miret-Roig et al. 2022a), are still needed. Our team is currently working on collecting and curating these data.

Acknowledgments

We thank the anonymous referee for the useful comments. J.O. acknowledges financial support from “Ayudas para contratos postdoctorales de investigación UNED 2021”. P.A.B. Galli acknowledges financial support from São Paulo Research Foundation (FAPESP) under grants 2020/12518-8 and 2021/11778-9. This research has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 682903, P.I. H. Bouy), and from the French State in the framework of the “Investments for the future” Program, IdEx Bordeaux, reference ANR-10-IDEX-03-02. We gratefully acknowledge the support of NVIDIA Corporation with the donation of one of the Titan Xp GPUs used for this research. We acknowledge Anthony Brown, the Gaia Project Scientist Support Team and the Gaia Data Processing and Analysis Consortium (DPAC) for providing the PyGaia code. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. The Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, the Queen’s University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation Grant No. AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the US Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High Performance Computing at the University of Utah. The SDSS website is www.sdss.org. SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics | Harvard & Smithsonian, the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU)/University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional/MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

References

Appendix A: Assumptions

In computing the properties of a physical group using its list of candidate members, we assume that this list provides a perfect classification. This assumption is optimistic, given that we have access only to a probabilistic classification in which a certain degree of entanglement is expected.

The positions and velocities of stars belonging to the same physical group are Gaussian distributed or can be modelled as a mixture of Gaussian distributions. Several Gaussian distributions pertain to the same physical group if their medians are mutually contained within one Mahalanobis distance. This latter is the multi-dimensional extension of the idea of measuring how many standard deviations a point is away from the mean of a distribution.

Gaia observes sources with a perfect selection function within its photometric completeness limit (G< 19 mag; see Tables 4 and 5 of Lindegren et al. 2021). This is a simplistic assumption given that Gaia misses sources in its expected magnitude coverage according to complex functions of its observables (see, for example, Boubert & Everall 2020).

Following Cottaar et al. (2015), we assume that the velocity dispersion at virial equilibrium, σ_vir can be computed as

$$\begin{array}{c}\hfill {\sigma }_{\mathit{vir}}=\sqrt{3\frac{G·M}{\eta ·r_{\mathit{hm}}}},\end{array}$$(A.1)

where G is the gravitational constant, η is a structural parameter, r_hm is the half-mass radius, and M is the total mass contained within r_hm. We added the factor $\sqrt{3}$ to obtain the total velocity dispersion rather than the 1D value computed by Cottaar et al. (2015).

Following Cottaar et al. (2015) and Portegies Zwart et al. (2010) we assume that r_hm and η in Eq. A.1 can be estimated by fitting an Elson et al. (1987) profile, and that η ∼ 10 for values of the EFF parameter γ > 2.5.

We assume that the dust and gas mass in the Perseus groups follows the observed 3D distribution of the stars and that the lower and upper limits to this mass contribution correspond to 65% and 169% of the stellar mass, respectively. Although this assumption will be generally correct for the young groups (i.e. IC348 and NGC1333), where the gas and stars are still coupled, it may not be entirely correct for the oldest groups like Alcaeus, where the stellar winds may have already cleared out some of its gas and dust. Nonetheless, this is the best we can do until the arrival of a 3D dust-and-gas mass map of the Perseus region.

We assume that the contribution of unresolved binaries to the total gravitational potential amounts to 20% of the inferred stellar mass. Generally, an unresolved binary is more luminous than a single star with the same mass as the primary of the binary. Then, it follows that our inferred masses, which are estimated based on the observed photometry of the stars, already partially mitigate the unaccounted mass of unresolved binaries. However, given the lack of information about the mass ratio of possible unresolved binaries within our list of candidate members, we can but make the conservative approach of increasing the inferred stellar mass.

When computing the internal energy distribution of the Perseus groups, we assume that the stellar system is self-gravitating. Given the negligible size of the Perseus groups when compared to the Galactic scale, it can be safely assumed that the Gravitational potential of the Galaxy is constant and thus can be removed from Eq. 2.

Appendix B: 3D velocities of the candidate members

In this Appendix, we show figures with the positions and velocities of our eight statistical Perseus groups. The figures show 2D projections of the Galactic Cartesian positions and velocities of the candidate members in each group. The velocity vectors are shown with respect to the group’s central velocity. For visual aid, their magnitude is also shown with the colour code.

Fig. B.1. 2D projections of the Galactic Cartesian positions (dots) and velocities (arrows) of the candidate members to the core of IC348. The velocities are computed with respect to the group’s central velocity. For visual aid, the colour code shows the magnitude of the 3D velocity.

Fig. B.2. Captions as in Fig. B.1 but for the halo of IC348.

Fig. B.3. Captions as in Fig. B.1 but for Heleus.

Fig. B.4. Captions as in Fig. B.1 but for Alcaeus.

Fig. B.5. Captions as in Fig. B.1 but for the core of Gorgophone.

Fig. B.6. Captions as in Fig. B.1 but for the halo of Gorgophone.

Fig. B.7. Captions as in Fig. B.1 but for NGC1333.

Fig. B.8. Captions as in Fig. B.1 but for Autochthe.

Appendix C: List of candidate members

All Tables

All Figures

	Fig. 1. Sky coordinates of the Perseus candidate members. The colour code shows the probabilistic classification, and the background image shows the thermal dust emission (545 GHz) from Planck Collaboration Int. LVII (2020).
In the text

	Fig. 2. Cartesian equatorial (ICRS) positions (top panel) and velocities (bottom panel) of the 1052 Perseus candidate members. The colour code shows the probabilistic classification, and the orange ellipses depict samples from the posterior distribution of the group level parameters corresponding to the one-sigma covariance matrix, and are centred at the mean position and velocity of the groups.
In the text

	Fig. 3. Cartesian equatorial (ICRS) positions (top panel) and velocities (bottom panel) of the Alcaeus group. The colour code shows the speed (top panel) and the distance (bottom panel) in the radial direction and both relative to the group centre. The orange ellipses show samples from the posterior distribution of the group-level parameters (see Fig. 2).
In the text

	Fig. 4. Empirical isochrones of the Perseus groups as obtained by Miec.
In the text

	Fig. 5. Absolute CMD of the Perseus groups candidate members (black dots). The black and coloured lines show the empirical isochrone and the theoretical ones of different evolutionary models, respectively.
In the text

	Fig. 6. Mass distributions of the Perseus groups. The orange lines show 100 realisations from the mass distribution obtained after transforming Miec’s magnitude distributions (see text). The rest of the coloured lines (those of the legend) depict the mass distributions computed with Sakam. The grey area shows the Gaia incompleteness region.
In the text

	Fig. 7. Extinction distributions of the Perseus groups. The colour code indicates the results of the theoretical and unified models as well as the uniform prior.
In the text

	Fig. 8. Rv distributions of the Perseus groups. Captions are the same as in Fig. 7. The vertical dotted line at Rv = 3.1 depicts the mode of the Gaussian prior.
In the text

	Fig. 9. Cumulative energy distributions of the Perseus groups. The dotted and dashed lines depict the results of accounting for the dust and gas mass by increasing the stellar mass by 65% and 169%, respectively.
In the text

	Fig. 10. Representation space coordinates of the confirmed (dots) and rejected members from the literature after removal of clear outliers (ϖ ∉ [1.9, 6.5] mas). The membership probability is shown as a colour code.
In the text

	Fig. 11. Cumulative mass distribution of IC 348 and NGC 1333.
In the text

	Fig. 12. Av and Rv 2D histograms for stars in the different Perseus groups.
In the text

	Fig. 13. Distributions of Av difference between the values inferred with Sakam and those reported by Bayestar19 (Green et al. 2019). Captions as in Fig. 7.
In the text

	Fig. 14. Galactic Cartesian positions and velocities of the Perseus groups (cones) together with the positions of individual members (dots). The interactive version of this figure will be available at http://www.project-dance.com/.
In the text

	Fig. B.1. 2D projections of the Galactic Cartesian positions (dots) and velocities (arrows) of the candidate members to the core of IC348. The velocities are computed with respect to the group’s central velocity. For visual aid, the colour code shows the magnitude of the 3D velocity.
In the text

	Fig. B.2. Captions as in Fig. B.1 but for the halo of IC348.
In the text

	Fig. B.3. Captions as in Fig. B.1 but for Heleus.
In the text

	Fig. B.4. Captions as in Fig. B.1 but for Alcaeus.
In the text

	Fig. B.5. Captions as in Fig. B.1 but for the core of Gorgophone.
In the text

	Fig. B.6. Captions as in Fig. B.1 but for the halo of Gorgophone.
In the text

	Fig. B.7. Captions as in Fig. B.1 but for NGC1333.
In the text

	Fig. B.8. Captions as in Fig. B.1 but for Autochthe.
In the text