© The Authors 2023

1. Introduction

Coma Berenices (hereafter Coma Ber, 85.5 ± 0.4 pc, Gaia Collaboration 2017) is the second closest open cluster to the Sun after the Hyades (47.03 ± 0.20 pc, Lodieu et al. 2019; 45.5 ± 0.42 pc, Gaia Collaboration 2017). Given its proximity, mid-age (between 550 and 800 Myr, e.g. Martín et al. 2020; Tang et al. 2018), and relatively large population (∼200 members, according to Tang et al. 2018), it is an excellent benchmark cluster for studies of stellar evolution and dynamical interactions. Since the pioneering work of Trumpler (1938), Coma Ber has been the subject of several studies. Hereafter, we mention only those works done in the Gaia (Gaia Collaboration 2016) era and refer the reader to the thorough membership review of Casal et al. (2020).

Tang et al. (2018) identify 192 candidate members within a 5° radius of the cluster centre using a combination of the Two Microns All Sky Survey (2MASS, Skrutskie et al. 2006), the UKIRT Infrared Deep Sky Survey (UKIDSS, Lawrence et al. 2007), the first USNO Robotic Astrometric Telescope catalogue (URAT1, Zacharias et al. 2015) and Gaia Data Release 2 (DR2, Gaia Collaboration 2018). These authors find that the cluster is mass-segregated and elongated along the Galactic plane. In addition, they identify nine sub-stellar members.

Fürnkranz et al. (2019), using Gaia DR2, applied the density based spatial clustering of applications with noise algorithm (DBSCAN, Ester et al. 1996) to the space of positions and tangential velocities resulting in the discovery of a kinematically cold group of 177 stars located in the vicinity of Coma Ber. They propose that this new group is younger (∼400 Myr) than Coma Ber, which formed in another part of the Galaxy, and that it will reach a minimum separation of 25 pc with Coma Ber in the following ∼13 Myr.

Tang et al. (2019) identify the tidal tail structures of Coma Ber using Gaia DR2 data and find that they have 120 members, extending up to ∼50 pc from the cluster centre. In addition, they also identify 218 candidate members of the new group identified by Fürnkranz et al. (2019), which they call Group X. They find evidence that this group was composed of two substructures. As mention by these authors, groups 10, 81, and 1805 from Oh et al. (2017) partially overlap with this new Group X.

Martín et al. (2020) performed a photometric search for new faint objects in the core of Coma Ber. They used UKIDSS, Panoramic Survey Telescope and Rapid Response System (PanSTARRS, Chambers et al. 2016), Sloan Digital Sky Survey (SDSS, Alam et al. 2015 DR12; and AllWISE, Wright et al. 2010) data to identify a couple dozen candidate brown dwarfs within the tidal radius of the cluster. Using spectral types and radial velocities derived from optical spectra collected with the Gran Telescopio de Canarias (GTC), these authors confirm the membership of four L2−L2.5 dwarfs, two of which were already known in the literature (Tang et al. 2019). Given the large Lithium (Li) depletion in the spectra of some of these Coma Ber members, Martín et al. (2020) propose that the Li depletion boundary must be located at spectral types later than L2.5, which restricts the lower limit of the cluster age to 550 Myr. Combining this constraining information with other dating techniques, these authors provide an age estimate of 780 ± 230 Myr for Coma Ber.

Casal et al. (2020) combined their own BVI_c photometry with Gaia DR2 parallaxes and rotational periods of K and M members of Coma Ber. They conclude that background giants from the field with short surface rotational periods could mimic main sequence cluster member stars and, thus, contaminate membership lists. The latter may result in older age estimates when using the isochrone fitting technique. These authors then argue that gyrochronology age estimates, as those of Collier Cameron et al. (2009, 590 ± 41Myr), should produce less biased results.

Sapozhnikov & Kovaleva (2021) performed an unsupervised blind search for candidate members in Coma Ber and Group X using Gaia Early Data Release 3 (EDR3, Gaia Collaboration 2021). As a result, they identify ∼250 candidate members for both Coma Ber and Group X. Also, they fitted isochrones and determine ages of 700 ± 70 Myr and 450 ± 100 Myr, respectively.

Souto et al. (2021) inferred the metallicity of the literature members of Coma Ber based on high-resolution near-infrared spectra from the SDSS-IV APOGEE survey (Majewski et al. 2017). Their results show that the Coma Ber members have near solar and homogeneous metallicities ([Fe/H] = +0.04 ± 0.02 dex).

Messina et al. (2022) used data from the Transiting Exoplanet Survey Satellite (TESS) mission (Ricker et al. 2015) to measure the rotation periods of 168 Group X candidate members from Tang et al. (2019). Then, using these periods and the gyrochronology age-mass-rotation relations of Angus et al. (2015) and Mamajek & Hillenbrand (2008), they derive the mean group’s age at 300 ± 60 Myr.

Newton et al. (2022) performed a search of Group X’s members using Gaia EDR3. They queried the 25 pc volume around the planet-host star TOI-2048, a Group X candidate member according to Tang et al. (2019), and identify – as candidate members of Group X – the 208 stars whose tangential velocity differed from that of TOI 2048 by less than 5 km s⁻¹. From this list, only 69 are in common with the candidate members from Tang et al. (2019). The authors also identify a potential new group, which they named MELANGE-2, containing 81 stars that share the 20 pc radius and 4 km s⁻¹ volume around TIC 224606446. The latter is a non-member of Group X that differs 5 km s⁻¹ from TOI 2048. The rotation periods of this putative new group suggest an age similar to that of Praesepe (670 ± 67 Myr; Douglas et al. 2019).

He et al. (2022) recently searched for new clusters by applying DBSCAN on Gaia EDR3. They identify 99 candidate members in Coma Ber and 61 in each of the two Group X substructures. These authors also determine isochrone ages of 560 ± 65 Myr, 178 ± 20 Myr, and 141 ± 16 Myr for Coma Ber, Group Xa, and Group Xb, respectively. We notice that their age estimates and the numbers of candidate members are at the lower edge of the literature values.

1.1. Group X or Latyshev 2

In 1977, Latyshev (1977) reported a new open cluster candidate comprised of seven stars located in the Ursa Majoris constellation. Later on, this candidate open cluster was named Latyshev 2 by Archinal & Hynes (2003). A few years ago, Mamajek (2016) deem this candidate open cluster a non-physical one given the large velocity dispersion of its putative members. However, Faherty et al. (2018) note that three stars belonging to Group 10 of Oh et al. (2017 which contains 29 members) were already reported as candidate members of an unnamed open cluster by Latyshev (1977). It turns out that this unnamed open cluster is Latyshev 2.

As mentioned above, Fürnkranz et al. (2019) report the discovery of a moving group close to the Coma Ber open cluster. Tang et al. (2019) also find this group and (re)named it Group X. Given that this group corresponds to Group 10 of Oh et al. (2017) and thus to the open cluster candidate reported by Latyshev (1977), hereafter, we refer to it as Latyshev 2 instead of Group X.

Coma Ber and its neighbouring group will experience a flyby in 13 Myr, according to Fürnkranz et al. (2019), or between 10−16 Myr, according to Tang et al. (2019). Given their proximity, relatively young ages, and unrelated origin, the possibility of an encounter in this pair makes it an interesting benchmark to study the effects of dynamical interaction(s) and possible merger(s) of stellar populations.

Several interacting stellar cluster pairs have recently been discovered in the Galaxy (e.g. Ye et al. 2022; Piatti & Malhan 2022; Angelo et al. 2022). However, none of these pairs is as close and populated as Coma Ber and Latyshev 2. Moreover, thanks to their proximity to the Sun and large populations, it is possible to obtain high-precision estimates of the pair’s parameters and flyby time, and thus provide constraining information for future studies of the population of OC pairs.

In this work, we aim to analyse the membership, luminosity, mass, phase-space (defined as the joint space of Cartesian positions and velocities), and energy distributions of Coma Ber and Latyshev 2 based on Gaia Data Release 3 (DR3, Gaia Collaboration 2023b) and complementary data. These precise and updated distributions will allow us to re-examine the future encounter of Coma Ber and Latyshev 2 and its consequences on the stellar populations of both clusters.

Section 2 introduces our Gaia DR3 data set, complementary photometry, and radial velocity measurements. Then, in Sect. 3, we present the methods by which we identified the cluster members and inferred the luminosity, mass, phase-space, and energy distributions. In this section, we also describe the trace-forward method by which we analyse the time of the future encounter between Coma Ber and Latyshev 2. In Sect. 4, we present our results, and in Sect. 5, we compare them with previous ones from the literature and discuss their implications. Finally, in Sect. 6, we present our conclusions.

2. Data

We downloaded¹ the astrometry and photometry of ∼45 million Gaia DR3 sources with the following constraints: b > 30°, and proper motions within −400 < pmra/[mas yr⁻¹]<  + 200 and −400 < pmdec/[mas yr⁻¹]<  + 400. These proper motion filters are permissive enough to encompass most open clusters and star-forming regions in the solar neighbourhood. We do not apply parallax cuts in the Gaia DR3 data to avoid biasing our sample of low-mass stars and brown dwarfs. In addition, we complemented the Gaia DR3 data of our candidate members with photometry from 2MASS (JHK) and PanSTARRS (grizy), as well as radial velocity measurements from the Apache Point Observatory Galactic Evolution Experiment DR17 (APOGEE, Abdurro’uf et al. 2022), the Set of Identifications, Measurements and Bibliography for Astronomical Data (SIMBAD, Wenger et al. 2000), and our own radial velocity measurements. The catalogue cross-match and data query with external catalogues were done based on the Gaia identifiers using Astropy (Astropy Collaboration 2013, 2018). The next two sections describe our spectroscopic observations and the data processing schema that we use to combine the radial velocity measurements from the different catalogues.

2.1. Mid-resolution optical spectroscopy

We obtained optical spectra of 36 of our bright (V < 15.5 mag) candidate members (see Sect. 4.1) as well as spectrophotometric and radial velocity standards with the Intermediate Dispersion Spectrograph² (IDS, using the RED+2 CCD). IDS is mounted at the Cassegrain focus of the 2.5m Isaac Newton Telescope (INT) on the Roque de Los Muchachos Observatory in La Palma, Spain (P.I. J. Olivares, programme C77). We refer the reader to Appendix for a detailed description of our INT-IDS observations.

We reduced the spectra using the spextractor tool (Bouy et al., in prep.). Briefly, it corrects the raw images with the median-combined bias and flat fields. Next, it applies cosmetic corrections (e.g. cosmic rays) and corrects for geometric distortions to align the spectra with the columns and rows of the detector. Then, it extracts the spectra using a rectangular aperture (of 2 FWHM width and sky subtraction of a region between 5 and 7 FWHM from the maxima) together with Gaussian and Moffat aperture profiles. Afterwards, it does the wavelength calibration with an interactive approach that allows the user to match the lines of the arc lamps (Copper-Argon and Copper-Neon in the case of INT-IDS) with their corresponding templates. Finally, it corrects for the instrumental response based on the spectra of the spectrophotometric standard stars.

Radial velocities were obtained by cross-correlating the calibrated 1D spectra of our targets with templates created from our radial velocity standards (see Appendix ). We selected these standards from the catalogue of Soubiran et al. (2013) based on the criteria of similarity of spectral type and proximity on the sky. We measured the radial velocities of all our 1D spectra (i.e. targets and standards) using the iSpec tool (Blanco-Cuaresma 2019; Blanco-Cuaresma et al. 2014) by first removing telluric lines and then applying the cross-correlation function. The templates of the radial velocity standards were created by shifting the 1D spectra with the barycentric correction and the radial velocity reported by Soubiran et al. (2013). Zero-point offsets were applied after correcting for the measured radial velocity of the standard star. We notice that our radial velocity measurements were obtained after the membership selection procedure and thus did not affect it. The typical uncertainty of our resulting radial velocity measurements is 10 km s⁻¹. Although these measurements are less precise than the median ones from APOGEE (0.12 km s⁻¹) or Gaia (1.8 km s⁻¹) catalogues, we will use them for those sources with missing radial velocities. The impact that largely heterogeneous and heteroscedastic data have on the cluster parameters will be handled by deconvolution methods (see Sect. 3.2). These types of methods offer an optimal solution for the extraction of the underlying distribution of the population parameters under the presence of heterogenous, heteroscedastic, and incomplete data sets (see, for example, Bovy et al. 2011).

2.2. Data processing

We processed our Gaia DR3 data by only applying a parallax zero point correction of −0.017 mas. As stated in Sect. 7 of Gaia Collaboration (2021), the current parallax bias correction is only a tentative recipe to correct this systematic.

The radial velocities were processed as follows. If a candidate member has both Gaia DR3 and APOGEE measurements, we use those of APOGEE due to their better precision. We computed zero-point offsets for the different surveys and conclude that given their small values and large uncertainties they can be neglected (see Appendix ). Given that our INT-IDS observations were taken before the release of Gaia DR3, 26 of our INT-IDS targets also have DR3 radial velocity measurements. Thus, we used their Gaia DR3 radial velocities since they resulted in smaller uncertainties. Some sources, particularly those from APOGEE, have radial velocity uncertainties of a few m s⁻¹ that considerably increase the convergence time for our phase-space modelling methods (see Sect. 3.2). For this reason, we set radial velocity uncertainties smaller than 0.1 km s⁻¹ to this latter value.

2.3. Initial list of members

As will be explained in Sect. 3.1 our membership methodology requires an initial list of members for the cluster under analysis. As this input list, we choose the 214 and 177 candidate members that Fürnkranz et al. (2019) identified in Coma Ber and Latyshev 2, respectively. Although those authors identified 214 candidate members in Coma Ber, the cross-match with Gaia DR3 excluded one Gaia source (source_id 3960475762679374208) due to its lack of parallax and proper motions. An initial inspection of the 6D space of the input list members revealed that 15 Coma Ber members were outliers in the velocity space. These sources were located more than six standard deviations beyond the candidate members’ mean U and W values. Therefore, we decided to exclude them from our input list. In the case of Latyshev 2, we do not find outliers in the initial list of members.

3. Methods

The following sections describe the methodologies we used to determine the different properties of stellar systems under analysis. We start by describing our novel membership methodology. Afterwards, we outline the steps to identify the phase-space structure, mass and energy distributions of Coma Ber and Latyshev 2. Finally, the section concludes by presenting the trace-forward method that we use to estimate the time of encounter between the stellar systems.

3.1. Membership selection

Our novel membership methodology, which we call Mecayotl³, is a generic one and thus, can be applied to open clusters, stellar associations, and moving groups in the solar neighbourhood. It proceeds by iteratively applying the following steps.

– In the first step, we use the cluster’s input list of members (i.e. their astrometry, radial velocities, uncertainties, and correlations) and the Kalkayotl⁴ code (Olivares et al. 2020, see Sect. 3.2) to determine the cluster’s phase-space model. The latter consists of the cluster’s central coordinates and size, both in the joint space of physical Cartesian positions and velocities.

– In the second step, we take the previous phase-space model and sample 10⁶ synthetic stars whose coordinates we transform to the observed space (i.e. ra, dec, parallax, pmra, pmdec, and radial_velocity) with the PyGaia⁵ code.

– In the third step, we use Gaussian Mixture Models⁶ (hereafter GMMs) with up to 80 components to independently fit the density distribution of the synthetic cluster’s stars and that of an equal-size sample of real field stars randomly selected from the input catalogue. The flexibility of GMMs allows us to model the complex shapes of both the cluster and field stars distributions. We use the GMM algorithms of Olivares et al. (2018b), which enable uncertainty deconvolution and the treatment of sources with missing values. The latter is of paramount importance to include sources with missing radial velocities. We select the best number of clusters and field GMM components using the Akaike information criterion (AIC, Akaike 1974). However, we restrict the search to only those models with at least ten stars in each mixture component. This robust restriction avoids spurious components fitted to outliers.

– In the fourth step, we use the previous two best GMMs, one for the field and another for the cluster, to compute membership probabilities of the input data set (i.e. all the real stars) as

(1)

where ℒ, ℳ, and 𝒫 stand for likelihood, model and prior, respectively. We use a uniform prior probabilities for the field and cluster.

– In the fifth step, the candidate members are selected based on probability thresholds optimised for the photometric magnitude of the sources (see Sect. 3.1.1 below).

– In the sixth step, we mask as missing the radial velocities of candidate members laying one standard deviation beyond the mean radial velocity of the cluster. This masking prevents the radial velocities of binary stars from inflating the velocity dispersion of the cluster while conserving good astrometric candidates.

Finally, the candidate members resulting from the sixth step are used as the input list for the next iteration of the algorithm, which we iterated until the number of candidate members converges over successive iterations. We used as convergence tolerance the Poisson uncertainty in the number of added members with respect to the previous iteration.

We handled the generation of the cluster’s synthetic Gaia observables as well as the analysis of the classifier quality with the free open-source code Amasijo⁷ (see, for example, Casamiquela et al. 2022 for another application of this code). The complete membership methodology described here will be publicly available as a free open-source code⁸ (Olivares et al., in prep.).

3.1.1. Optimum probability thresholds

It is well known that the Gaia astrometric uncertainties heavily depend on the apparent magnitude of the sources (see, for example, Fig. 7 of Lindegren et al. 2021). Therefore, any classifier incorporating Gaia astrometric uncertainties will also show a correlation of its quality with the photometric magnitude of the sources to which it is applied. To minimise this correlation’s impact and obtain the best classification at different magnitude intervals, we estimated the optimum probability threshold at each g magnitude bin in the full Gaia domain (i.e. from 2 to 21 mag). We did this at each iteration of our membership algorithm by generating ten synthetic clusters with 10³ sources each. In these clusters, we randomly assigned masses from a uniform distribution and obtained their photometry using the isochrones⁹ code (Morton 2015), with the cluster’s observed age and metallicity (see Sect. 1). Afterwards, using the photometry and astrometry, we computed observational uncertainties with PyGaia and randomly sampled the observed values from these uncertainties. Finally, we applied the fourth step of our membership algorithm to the data sets resulting from combining the synthetic clusters with the real field stars. The obtained membership probabilities together with the true labels of the sources (i.e. cluster or field) allow us to compute classification confusion matrices¹⁰ and quality indicators, such as the true positive rate (TPR) or the contamination rate (CR), at each magnitude bin and probability threshold. We used magnitude and probability grids with 8 and 40 steps, respectively. The probability grid spans from one to five sigma, with each interval containing ten steps. At each grid step, we selected the optimum probability threshold as the one that maximises the Mathews correlation coefficient (MCC, Matthews 1975) of the classifier. This criterion measures the correlation coefficient between the true and observed classifications, and it is known to be the most objective metric to evaluate the quality of a classifier in the presence of unbalanced classes, as is the case here. We also tested other metrics (e.g. Accuracy and F1M), obtaining similar probability thresholds. We notice that due to our conservative approach for estimating contamination rates, which assumes as field stars all sources not classified as cluster members by either the literature or a previous iteration of the algorithm, possibly true cluster members will be considered contaminants and thus inflate our estimates of the contamination rate.

3.1.2. Minimising contamination

During the validation of our membership methodology, we noticed that the contamination rate of successive iterations of the algorithm tended to increase. We inspected the contaminants and found that these resulted from two main sources: radial velocity outliers and faint sources. On the one hand, the radial velocity outliers inflate the velocity dispersion of the cluster and thus result in accreted contaminants in the following iterations of the algorithm. On the other hand, at the faintest end of the cluster photometric sequence, the numbers and uncertainty of the sources increase exponentially. Both these factors result in large chances of confusing field sources as cluster members.

We minimised the contamination rate of our membership methodology and thus prevent possible biases by applying the following two filtering criteria. First, in the sixth step of our membership methods, we masked the radial velocities of outliers as missing. We opted for masking with a rather conservative 1σ filtering than discarding good astrometric candidate members. Second, we inferred the phase-space cluster model based only on sources within magnitudes bins with contamination rates < -pagination10%. This effectively discarded from the inference of the phase-space cluster model those sources fainter than 20 mag in g band. Although these sources were not used in the inference of the cluster phase-space model, they remain in the list of cluster candidate members. Due to our conservative approach to estimating the phase-space model, which filters out sources with G > 20 mag, our membership algorithm may have missed some of the faintest cluster members.

3.2. Phase-space structure

We analysed the phase-space structure of the stellar groups using the Kalkayotl code (Olivares et al. 2020). In its 6D version (Olivares et al., in prep.), the code simultaneously infers the Cartesian (ICRS) positions and velocities of the stellar system and its members, given the observed Gaia astrometry (including its uncertainties and correlations) and (possibly missing) radial velocities. The cluster or group-level parameters are modelled with a hierarchical prior, whereas the stellar or source-level parameters are drawn from this prior. Thanks to the hierarchical approach, the influence of the prior is minimised by inferring its parameters at the same time as those of the individual stars. Nonetheless, the prior family (i.e. the type of distribution) must be set in advance. The prior distribution should be selected based on the a priori knowledge about the properties of the stellar system. Given that we expected a certain degree of contamination in the list of members and possible substructure in the positions and velocities of these systems, we decided to test the three types of prior families provided by Kalkayotl. These are Gaussian, Gaussian mixture models (GMMs), and concentric Gaussian mixture models (CGMMs). They all have the joint 6D space of Cartesian positions and velocities as a domain.

It is expected that the velocity distributions of these relatively old (∼300 Myr for Latyshev 2 and 800 Myr for Coma Ber) stellar systems have already been relaxed and thus can be described with a single Gaussian distribution. However, the presence of structures, such as tidal tails or halos, requires more flexible models. While the CGMM imposes a common location for all the Gaussian components of the mixture, the GMM is the most flexible of our models with completely free weights (restricted to add to one), means, and covariance matrices. To choose among these, we used the convergence of the inference algorithm and non-negligible components weights as selection criteria. We tested GMM with one to four components, starting with the simplest one. If the model with additional components failed to converge or had weights lower than 10%, we step down to the previous simple one.

Our phase-space modelling is conceptually simple but allows constructing complex models based on Gaussian building blocks. It is not the first time that this rather simple multi-component approach has been successfully applied in the literature. For example, Tarricq et al. (2022) fitted the 3D spatial distribution of open clusters using CGMMs in which, similarly to what we did here, one Gaussian component describes the cluster’s core while the other two may be used to describe the halo and tails of the cluster.

We noticed the following advantages of our GMM phase-space modelling. First, it allowed the inference of the Cartesian positions and velocities for all the cluster members, even those without radial velocity measurements. This latter is accomplished assuming that the cluster members share the global cluster’s velocity vector but differ only by the cluster’s internal velocity dispersion. We refer the reader to the works of Madsen et al. (2002), Lindegren et al. (2000), Dravins et al. (1999) for the conceptual definition of astrometric radial velocities and its applications to nearby open clusters and moving groups. Second, GMMs can be used as classifiers; thus the identified candidate members can be partitioned into the substructures corresponding to each of the Gaussian building blocks of the GMM. This classification is based on the source’s multidimensional observed data and the GMM’s inferred parameters. Finally, our GMMs live in the joint space of positions and velocities, which allows them to capture correlations among positions and velocities. The latter represented an improvement over previous models living in independent spaces of positions and velocities (e.g. Lodieu et al. 2019; Jerabkova et al. 2021).

3.3. Luminosity and mass distributions

We determined the luminosity and mass of each candidate member using the PAdova TRieste Stellar Evolution Code (PARSEC, Pastorelli et al. 2020; Marigo et al. 2013) and BT-Settl (Baraffe et al. 2015; Allard 2014) theoretical isochrone models. We utilised a grid of ages spanning the literature values (see Sect. 1): Coma Ber ∈ [700 Myr–900 Myr] and Latyshev 2 ∈ [200 Myr–400Myr], both grids with steps of 100 Myr. We inferred the luminosity and mass distributions with the free and open-source code Sakam¹¹ (Olivares et al. 2019). It infers the joint posterior distribution of luminosity or mass together with A_v and R_v for each candidate member given its distance (see Sect. 3.2) and observed photometry (see Sect. 2). We worked with the following prior distributions. For the luminosity and mass, we used a uniform distribution in the logarithmic space and the Chabrier (2005) distribution, respectively. We used a Gaussian distribution for the total-to-selective extinction: R_v ∼ 𝒩(μ = 3.1, σ = 0.5). In the case of A_v, we used a uniform prior, with A_v ∈ [0, 0.1] mag. Appendix provides a detailed explanation to support this decision. Briefly, this is the most probable extinction interval according to the values reported by Planck Collaboration Int. XLVIII (2016) for the sky positions of our candidate members.

We notice the following points. First, in the case of unresolved binaries or multiple systems, the inferred mass for each candidate member will correspond to the system mass distribution. We had no method to infer the individual masses of possible unresolved binary or multiple systems. Second, none of the theoretical isochrone models that we used here spans the full magnitude interval of our candidate members. Therefore, we inferred masses using not only the BT-Settl and PARSEC models independently but also the join of these two, which we call the PB model. For the latter, we selected as the joining point a mass of 1.4 M_⊙, which is the one where the two isochrone models reach their maximum agreement (see, for example, Lodieu et al. 2019).

3.4. Energy distribution

The energy distributions will allow us to identify physical members gravitationally bound to the stellar system as well as possible escapers. In the following, we describe the method to quantify the total energy of each of our candidate members.

We obtained the energy distribution of each candidate member with respect to its parent stellar system by adding its kinetic and potential energies. These energies were computed from the candidate’s posterior distributions of mass, positions, and velocities inferred with the methods of the previous sections. Then, we took 1000 samples from these posterior distributions, translate the positions and velocities to the reference frame of the stellar system, and computed the sample’s total energy as the sum of its kinetic and potential energies. By taking samples of the posterior distribution, we avoided assuming a particular family distribution and thus the need for Jacobian computations.

The total energy of each sample was computed as

(2)

where r and v are the distance and speed in the reference frame of the stellar system, M is the total mass enclosed within the distance r, m is the sample’s mass, and G is the gravitational constant.

We notice the following assumptions. First, given the low photometric extinction of these two groups (A_v ≲ 0.75 mag, see Sect. 4.3), we assumed no contribution of the dust and gas to the total mass of these systems. Second, our mass estimates for individual candidate members are most likely underestimated in the case of unresolved binary or multiple systems. Thus, we assumed that the contribution of binary stars and unresolved systems to the total mass of the stellar group is 20% of its inferred mass. The estimated value for the fraction of binary systems in open clusters varies from 11% to 70% (Sollima et al. 2010), with the most typical value for unresolved binaries being ∼20% (e.g. Jadhav et al. 2021). Third, in the energy computation, see Eq. (2), we assumed that the stellar systems are self-gravitating, and thus we neglected the contribution of the Milky Way’s gravitational potential. This assumption was valid for establishing each candidate member’s energy relative to its parent group. Due to the negligible size of these systems compared to the Galactic scale, we assumed that the Milky Way’s gravitational potential can be treated as constant and thus removed from the energy computation. This assumption was dropped in the next section, where our interest lies in tracing our candidate members’ present-day positions and velocities forward in time.

3.5. Trace-forward

To determine the time of the flyby between Coma Ber and Latyshev 2, we traced the present-day phase-space coordinates inferred with Kalkayotl (see Sect. 3.2) forward into the future. We did this not only the group-level parameters (i.e. mean positions and velocities of the groups) but also the source-level parameters (i.e. positions and velocities of the candidate members). The uncertainties of these phase-space coordinates were taken into account by tracing 100 samples from the posterior distributions of the group- and source-level parameters forward into the future. We assumed that these samples are particles subject only to the Milky Way’s gravitational potential, and in the following, we refer to them as sample-particles.

We integrated the orbit of each of these sample-particle 20 Myr into the future using the Galpy code (Bovy 2015), a time grid with steps of 0.5 Myr and assuming as gravitational potential the Galpy’s MWPotential2014. The latter assumes a rotational velocity of 220 km s⁻¹ at a solar radius of 8 kpc (for specific details of this gravitational potential, see Bovy 2015).

Finally, for each group-level sample-particle, we determined the time of the flyby as the time at which the distance between the two groups reaches its minimum. We preferred to work with the group-level parameters to estimate the flyby time rather than the source-level parameters because the former have better precision than the latter. We report the time of encounter with statistics (mean and standard deviation) of the distribution of flyby times. Once the flyby time was determined, we use the source-level sample particles to analyse the degree of attachment of the individual stars to their parent group and the rest of the groups and thus test the hypothesis of mixing of stellar populations.

4. Results

This section presents the results we obtained after applying the methods of Sect. 3 to the data introduced in Sect. 2. The comparison of these results with those from the literature will be done in Sect. 5.

4.1. Membership

We applied our novel membership algorithm (see Sect. 3.1) to the ∼45 million Gaia DR3 sources (see Sect. 2), using as input lists of members those described in Sect. 2.3. In Coma Ber, our membership algorithm converged in six iterations to 302 candidate members, while in Latyshev 2, it converged to 332 candidate members in five iterations. Table D.1 contains the identifiers and group classification for all our candidate members together with their properties as inferred throughout the rest of this work.

Tables 1 and 2 show the quality indicators of the Coma Ber and Latyshev 2 classifiers at the final iteration of the membership algorithm. The first line of the tables (column Strategy All) shows, only for comparison purposes, the properties of the classifier if no magnitude bins were used for the optimisation. On the other hand, the remaining lines of the tables (Bins 1–9) show the classifier properties measured as a function of the magnitude bins. The edges of these bins were fixed at 2, 6, 8, 10, 12, 14, 16, 18, 20, and 22 magnitudes in the Gaiag band. The remaining columns show, for each magnitude bin and from left to right, the mean value of the sources’ g magnitude, the optimal probability threshold (truncated to two digits), the total number of sources, the entries of the confusion matrix (True Positives, TP, False Positives, FP, True Negatives, TN, and False Negatives, FN), and the quality indicators of True Positive Rate (TPR), Contamination Rate (CR) and MCC (see footnote¹⁰ and Sect. 3.1.1). We notice that the values shown correspond to the mean values computed from the ten synthetic clusters used to establish the optimum probability thresholds (see Sect. 3.1.1). As can be observed from these tables, the quality indicators of the membership methodology are excellent, with true positive rates ≳99% and contamination rates ≲5%. The only exception is at the faintest magnitude bin ( g > 20 mag), where the TPR of the Latyshev 2 classifier drops to 95%, and the CR of the Coma Ber classifier grows up to 44%.

The astrometric quality of our candidate members is briefly summarised as follows. The parallax fractional error is typically (mode) better than 1% with only one extreme value of 20% for a Coma Ber tails member. The mean proper motions uncertainties are 0.16 mas yr⁻¹, 0.08 mas yr⁻¹, 0.09 mas yr⁻¹, and 0.10 mas yr⁻¹ for Coma Ber’s tails, core, Latyshev 2, and Mecayotl 1, respectively. With respect to the RUWE parameter (see, for example, Lindegren et al. 2021), the typical (mode) values for ComaBer’s core, tails, Latyshev 2 and Mecayotl 1 are 1.05, 1.06, 1.05, and 1.03, respectively, while their cumulative distributions reach the 1.4 threshold at 85%, 80%, 86% and 83%, respectively. A RUWE value larger than 1.4 may indicate that the source is non-single or has a problematic astrometric solution¹². Concerning the astrometric excess noise parameter (astrometric_excess_noise, ϵ_i) and its significance (astrometric_excess_noise_sig, σϵ_i), which are measures of the disagreement between the observations of a source and its best fitting astrometric model (see Sects. 3.6 and 5.1.2 of Lindegren et al. 2012), the median values in Coma Ber’s core are ϵ_i = 0.13 mas and σϵ_i = 9.36, in Coma Ber’s tail are ϵ_i = 0.13 mas and σϵ_i = 2.58, in Mecayotl 1 are ϵ_i = 0.19 mas and σϵ_i = 3.70, and in Latyshev 2 are ϵ_i = 0.12 mas and σϵ_i = 3.46. However, the Gaia DR3 non-single star catalogue (NSS, see, for example, Gaia Collaboration 2023b) classifies only nine, two, and four of our candidate members of Coma Ber, Latyshev 2, and Mecayotl 1, respectively, as potential binaries (i.e. their non_single_star entry is non zero).

4.2. Phase-space structure

As explained in Sect. 3.2, we fitted three different models, Gaussian, CGMM, and GMM, to the observed astrometry and radial velocity (when available, see Sect. 2) of our 302 Coma Ber and 332 Latyshev 2 candidate members. The best phase-space models that describe these observables are a CGMM with two components for Coma Ber and a GMM with two components for Latyshev 2. As explained in Sect. 3.2, we selected the best model based on the criteria of the convergence of the inference algorithm and the non-negligible weights of the mixture components. Our attempts to fit models with three and four components returned negligible weights for the third component or failed to converge in the case of four components.

We observed that one candidate member of Coma Ber (source id: 4016309615972923648) and five of Latyshev 2 (source_id: 1608752435341695872, 1560784144636654976, 857806987370330624, 1560938316782309248, and 1512261012875442304) were located more than three sigma away from the mean UVW velocities of the clusters. To prevent these sources from biasing the internal velocity dispersion of the clusters, we masked their radial velocities as missing. We check the binary status of the previous sources and found that only source 1608752435341695872¹³ is a confirmed binary star according to SIMBAD, while sources 1560938316782309248¹⁴ and 1560784144636654976¹⁵ are most likely binary stars. However, the remaining three sources have RUWE values between 1.06 and 1.14 and astrometric_excess_noise between 0.03 and 0.16, which are not conclusive enough to claim binarity. None of these six stars appears in the Gaia DR3 NSS sample (i.e. their non_single_star entry is zero), which indicates that according to Gaia, they are neither astrometric, spectroscopic nor eclipsing binaries. However, due to the low completeness of the NSS sample (see Sect. 3 of Gaia Collaboration 2023a) and the discrepant radial velocities of these sources, we cannot rule out the possibility that they are spectroscopic binaries.

In the fitted CGMM of Coma Ber, we observe that one of the Gaussian components serves to describe the cluster core while the other one describes the cluster tails. In the fitted GMM of Latyshev 2, the two Gaussian components are mutually exclusive in the joint 6D space of positions and velocities (more below). Thus, we identify one of these components with the original candidate open cluster Latyshev 2 and the other with a new stellar association that we call Mecayotl 1.

Figures 1 and 2 show, with orange lines, 100 posterior samples from the inferred group-level parameters of Coma Ber’s CGMM and Latyhsev 2 + Mecayotl 1 GMM phase-space models, respectively. These figures show multidimensional projections of the phase-space, with the top panels depicting positions and the bottom ones velocities. The dots and error bars in the figures represent the mean and standard deviation of each candidate member’s posterior samples of its phase-space coordinates. The colour code of the individual sources indicates their probabilistically assigned parent structure, which, as mentioned in Sect. 3.2 corresponds to the classification into each of the Gaussian of the GMM. We notice that this classification is done in the 6D space and is based on the model parameters and the source’s observed data. This probabilistic classification results in 186 candidate members in Latyshev 2, 146 candidate members in Mecayotl 1, and 136 and 166 candidate members in the core and tails of Coma Ber, respectively.

Fig. 1. Positions (top panel) and velocities (bottom panel) of Coma Ber candidate members. The 100 orange lines show the posterior samples from the inferred group-level parameters while the dots and error bar depict the mean and standard deviation, respectively, of the inferred positions and velocities.

Fig. 2. Positions (top panel) and velocities (bottom panel) of the candidate members of Latyshev 2 and Mecayotl 1. Captions as in Fig. 1.

Figure 3 shows the distance distribution of our identified physical groups. These densities were computed by applying a kernel density estimate to the aggregation of 100 posterior samples of each candidate member. As can be observed, the core of Coma Ber is the closest group while the Mecayotl 1 one extends up to 140 pc away.

Fig. 3. Distance distributions of Coma Ber, Latyshev 2 and Mecayotl 1.

Table 3 shows the mean and standard deviation of the inferred group-level parameters of the CGMM of Coma Ber’s core and tail and the GMM of Latyshev 2 and Mecayotl 1. For completeness, the bottom five lines of this table show the total position and velocity dispersions of the inferred phase-space models (σ_XYZ and σ_UVW, respectively), two kinematic proxy indicators, one for expansion and the other for rotation, and the tidal radius (r_t). The latter corresponds to the distance from the group’s centre at which the gravitational pull exerted by the Galactic potential equals that of the group’s potential under the assumption that the latter follows a circular orbit around the centre of the former. We computed this tidal radius using Eq. (1) of Tang et al. (2019) with the Oort’s constants values of those authors. The total position and velocity dispersions were computed by adding the dispersion of each coordinate in quadrature. Finally, the proxy indicators for expansion, , and rotation, , are computed as the mean values of the dot and cross products of the unit radial vector, , and the velocity vector, v; these latter two computed in the cluster’s reference frame (see, for example, Galli et al. 2019).

Concerning Coma Ber, the top panel of Fig. 1 shows that its tail members reach larger distances, with the farthest candidate member reaching up to 95 pc from the centre, which is roughly two times distance of the latest result from the literature (Tang et al. 2019). In the velocity space (bottom panel of Fig. 1), we observe that the tail members have similar but not identical velocities to those in the core.

Concerning Latyshev 2 and Mecayotl 1, both panels of Fig. 2 show that they have distinct positions and velocities. Although their members appear to be entangled, their 6D models are not. Quantitatively, we measure the degree of separation between Latyshev 2 and Mecayotl 1 with the Mahalanobis metric. This metric is the multidimensional extension of the distance to the mean value in a one-dimensional Gaussian distribution when measured in units of the standard deviation. We notice that these distances are not commutative since the internal dispersions of the groups are not the same. Briefly, the Mahalanobis metric can be understood as the multi-dimensional generalisation of measuring the number of standard deviations to the mean value. The Mahalanobis distance that Latyshev 2 has with respect to Mecayotl 1 is 2.03, whereas that of Mecayotl 1 with respect to Latyhsev 2 is 5.5, thus indicating that these groups are mutually exclusive at a level > 2σ. This significance level indicates that our discovery of the two groups is statistically significant, with credibility > 95%. Furthermore, our claim that Latyshev 2 is an open cluster while Mecayotl 1 is a stellar association is based on their velocity dispersion. On the one hand, Latyshev 2 is a compact group in the velocity space with a total velocity dispersion of σ_UVW = 0.80 ± 0.08 km s⁻¹, which is similar to that of Coma Ber’s core (σ_UVW = 0.89 ± 0.10 km s⁻¹) and that of the archetypical Pleiades open cluster (0.8 ± 0.1 km s⁻¹ according to Galli et al. 2017; and 0.83 ± 0.07 km s⁻¹ according to Torres et al. 2021, as computed from their 1D radial velocity dispersion of 0.48 ± 0.04 km s⁻¹). On the other hand, Mecayotl 1 has a total velocity dispersion of σ_UVW = 1.85 ± 0.18 km s⁻¹, which is more than two times the velocity dispersion of Latyshev 2 or the core of Coma Ber, and similar to the 1.73 ± 0.3 km s⁻¹ of the βPic stellar association (Miret-Roig et al. 2020). Further arguments for these classifications will be provided in Sect. 4.4.

4.3. Luminosity and mass distributions

To determine the luminosity and mass distributions of the identified physical groups, we first verify that the ages proposed by the literature are indeed compatible with the photometry of our candidate members. Then we compute the luminosity and mass of the individual candidate members and combine them to obtain the corresponding distributions for the physical groups.

Figure 4 shows the G vs G − RP absolute colour-magnitude diagrams of our candidate members, together with the theoretical isochrones from the PARSEC and BT-Settl models at our chosen grids of ages (see Sect. 3.3). In addition, for Mecayotl 1 we also plot the PARSEC isochrones for up to 700 Myr. As can be observed from this figure, the age intervals proposed in the literature are compatible with the empirical isochrones traced by our candidate members.

Fig. 4. Absolute colour-magnitude diagrams of the candidate members of Coma Ber, Latyshev 2 and Mecayotl 1. The lines depict the theoretical isochrone models with colour-coded ages.

In the case of Coma Ber, our brightest candidate members support the ages of 700 to 800 Myr, as Martín et al. (2020) proposed based on lithium measurements and Tang et al. (2019, 2018) on isochrone fitting, respectively. In the case of Latyshev 2, the three brightest members originally proposed by Latyshev (1977) and confirmed here serve as anchors for the isochrone age of the system. These three members seem to indicate that this cluster is young, most likely aged 300 to 400 Myr, which supports the gyrochronology age of 300 ± 60 Myr proposed by Messina et al. (2022) rather than the isochrone fitting ages of (450 ± 100 Myr Sapozhnikov & Kovaleva 2021) and He et al. (2022 141 ± 16 Myr). On the other hand, the brightest candidate members of Mecayotl 1 show that its isochrone age is between 400 Myr and 600 Myr. Under the lack of constraining evidence, in the following, we assume an age of 400 Myr for this stellar system. As can be observed in Fig. 4, this age corresponds to the closest lower-envelope isochrone to the observed photometry of our candidate members.

As mentioned in Sect. 3.3, we compute the individual luminosity and mass distributions of the candidate members given the theoretical isochrone models and the observed photometry. Then, we combine the resulting posterior samples of all candidate members into the luminosity and mass distributions of the physical groups. Given the excellent recovery rate of our membership algorithm (> 95%, see TPR column in Tables 1 and 2), we do not correct the resulting luminosity and mass distributions for possible missing members. Similarly, given the low contamination rate of our membership algorithm (< 7%, see CR column in Tables 1 and 2) we do not correct the luminosity and mass distributions for contamination. However, the faintest (g > 20 mag) magnitude bin has large contamination rates (∼40%), which indicates that the luminosity and mass distributions at this magnitude bin can only be taken as upper limits to the true distributions.

Figure 5 shows the luminosity distributions of the identified groups as inferred from different models (colour coded). The differences that result from inferring the luminosity and mass distributions using the grid of ages mentioned in Sect. 3.3 are negligible. Thus, to avoid overcrowding the plots, we show only the most likely ages: 300 Myr for Latyshev 2, 400 Myr for Mecayotl 1, and 800 Myr for Coma Ber. The lines shown in Fig. 5 result from applying a kernel density estimate to 1000 aggregated samples from the posterior luminosity distribution of each of our candidate members. As can be observed, the luminosity distributions resulting from different theoretical models agree well, except at the limit of the BT-Settl model, where there is a concentration of samples from sources more luminous than the BT-Settl limit.

Fig. 5. Luminosity distributions inferred with different theoretical models. The grey line and grey region depict the uniform prior and the completeness limit of the Gaia data. This latter corresponds to apparent G mag ∼19 transformed to luminosity using the mean group distance and BT-Settl model.

Similarly to Figs. 5 and 6 shows the mass distributions inferred for Coma Ber, Latyshev 2 and Mecayotl 1 using the same models and ages as for the luminosity. The lines also show the kernel density estimate of the aggregated posterior samples. The mass prior (grey lines) corresponds to Chabrier (2005) mass distribution. As can be observed, the mass distributions resulting from different models have varying degrees of agreement. In the core of Coma Ber, the mass distributions agree well for all models whereas, in the rest of the cases, the peak of the mass distributions resulting from the PARSEC model are shifted 0.06 M_⊙ with respect to those of the BT-Settl and PB models. Given that the BT-Settl models have been specifically created for low-mass stars and brown-dwarfs and the bulk of our candidate members is made of low-mass stars, in the following, we use the mass estimates obtained from the unified PB model, this is BT-Settl for masses < 1.4M_⊙ and PARSEC for masses > 1.4 M_⊙ (see Sect. 3.3). We select this joint model because it is the only one covering the entire magnitude interval of our candidate members. Concerning the ages, we also observe negligible differences in the mass distributions resulting from using ±100 Myr from the selected ones.

Fig. 6. Mass distributions inferred with different theoretical models. Caption as in Fig. 5. In this case, the grey line depicts the Chabrier (2005) mass prior.

The luminosity distributions of Coma Ber, Latyshev 2 and Mecayotl 1 all show a dip at log L/L_⊙​ ≃ ​ − 0.89, which corresponds to the Wielen dip (Wielen 1974) and to a mass of ≃0.7 M_⊙ (Kroupa et al. 1990). This dip has been consistently observed in the solar neighbourhood (< 20 pc) and in young (< 50 Myr) open clusters (see Guo et al. 2021 and references therein), and even in the old (2.5 Gyr) open cluster Ruprecht 147 (Olivares et al. 2018b). We notice that this dip is less prominent in the luminosity distribution of Latyshev 2 as well as in its mass distribution. A similar situation occurs with the mass distribution of the core of Coma Ber, although in this case the dip is clearly observed in its luminosity distribution.

As mentioned in Sect. 3.3, when inferring the luminosity or mass distributions, the Sakam code also delivers posterior distributions of A_v and R_v. In the case of A_v, the inferred values are all consistent with the chosen prior (i.e. A_v < 0.1 mag, see Sect. 3.3). Concerning R_v, the posterior distributions are also consistent with the prior (i.e. R_v = 3.1) in spite of the different theoretical isochrones used. Furthermore, it shows no apparent correlation with A_v, the luminosity, or the mass of the candidate members.

4.4. Energy distribution

We now compute the energy distributions of our candidate members with the methods described in Sect. 3.4. Figure 7 shows each group’s cumulative energy distribution function (CEDF). In this figure, each orange line depicts the CEDF resulting from taking one sample of the posterior distributions of position, velocity and mass. We choose to plot the cumulative distribution function rather than the probability distribution function because it allows direct reading of the fraction of sources within a given energy value. For example, it allows us to directly read that only 20–30% of the sources in the core of Coma Ber have negative energies.

Fig. 7. Cumulative energy distribution functions. Each orange line shows the result obtained from a sample of the posterior distributions of mass, position and velocity. For visual aid, the vertical grey dashed lines show the zero energy.

We observe that the CEDF of Latyshev 2 resembles that of Coma Ber’s core while the one of Mecayotl 1 resembles that of Coma Ber’s tail. On the one hand, the CEDFs of Latyhsev 2 and the core of Coma Ber have steep slopes, with 80% of its sources having energies smaller than 0.3 M_⊙ km² s⁻² and 95% having energies smaller than 1.0 M_⊙ km² s⁻². On the other hand, Mecayotl 1 has a large energy dispersion, with only 40% of its sources having energies smaller than 0.3 M_⊙ km² s⁻² and 75% having energies smaller than 1.0 M_⊙ km² s⁻². This energy dispersion is even larger than that of the tails of Coma Ber, where 85% of the sources have energies smaller than 1.0 M_⊙ km² s⁻². The previous results provide additional evidence supporting our classification of Latyshev 2 as an open cluster and Mecayotl 1 as a stellar association. Nonetheless, we notice that the low fraction of bound sources in these open clusters indicates that they are disrupted, as in the case of Latyshev 2, or in the process of disruption, as in the case of the core of Coma Ber. We come back to these points in Sect. 5, where we attempt to link the properties of these systems with their possible origin.

Figure 8 shows the 2D kernel density estimate of energy versus distance to the group centre. This figure will help us to explore possible correlations between energy and position.

Fig. 8. Density of candidate members as a function of energy and distance to the group centre.

Concerning Coma Ber, we observe that the core candidate members show no clear correlation between the energy and distance to the centre, with only a slow decay in the density of sources as a function of distance to the centre. On the contrary, the tidal tail candidate members show a clear correlation between energy and distance, with more distant sources resulting in larger energies. Interestingly, the candidate members do not uniformly distribute along the distance-energy diagram but agglomerate at a distance of 10–15 pc.

Concerning Latyshev 2, similarly to the core of Coma Ber, it shows no clear correlation between distance to the centre and energy, with only a smooth decay in the density of sources with distance to the centre, which we interpret as further evidence of its similarity to the core of Coma Ber. On the other hand, Mecayotl 1 shows large dispersion in energy with over-densities in the distance to the centre. We spot at least two clear over-densities and hints of a possible third one. The two most prominent ones are uniformly spread at distance intervals of 5–15 pc and 30–35 pc, whereas the third and less prominent one is located at 45 pc. Furthermore, we observe that these three over-densities have energy dispersion larger than that of Coma Ber’s tails. We discuss the implications of these over-densities and their dispersion in Sect. 5.

4.5. Trace-forward

As mentioned in Sect. 3.5, we use the source- and group-level parameters inferred with Kalkayotl to trace the current positions and velocities of our candidate members forward in time. In the following, we first obtain the time of flyby between Coma Ber and the other two groups using the time integration of their group-level parameters. Then, we estimate the candidate members’ binding energies with respect to the other groups using the source-level parameters. These binding energies with respect to all the groups will allow us to investigate possible population mixing.

4.5.1. Flyby time

We integrate forward in time (see Sect. 3.5) 100 sample-particles drawn from the posterior distribution of the group central phase-space coordinates. Although we give these central coordinates to Galpy in the ICRS frame, it internally transforms them into a Galactic frame (R, Z, and ϕ; see Bovy 2015 for the definition of these Galactic coordinates). Figure 9 shows the time evolution (colour coded) of the central Galactic coordinates of Coma Ber, Mecayotl 1 and Latyshev 2. The figure shows that these latter two will experience an abrupt divergence in the following 10 Myr. This divergence will mainly occur in their galactocentric distance, R.

Fig. 9. Galactic coordinates (R, Z, ϕ) of Coma Ber, Mecayotl 1 and Latyshev 2 as functions of time. To improve visibility, only the orbits of 20 samples of each stellar system are shown.

To obtain the time of flyby, we compute for each sample-particle and, as a function of time, the euclidean distances among the three groups. We notice that the core and tails of Coma Ber have, by construction, the same central position; thus, there is no need to treat them as different groups. The closest distance between Coma Ber and Latyshev 2, 28.4 ± 1.5 pc, will occur in 11.3 ± 0.5 Myr, between Coma Ber and Mecayotl 1, 23.6 ± 2.0 pc, will occur in 14.0 ± 0.6 Myr, and that between Mecayotl 1 and Latyshev 2, 12.5 ± 2.0 pc, will occur in 8.1 ± 1.3 Myr.

4.5.2. Energies as a function of time

We now compute the total energy of each sample-particle with respect to its parent physical group and with respect to the rest of the groups. While the former helps us to analyse the future state of these systems, the latter helps us to spot sources that may be captured by a group different from its parent one. The presence or absence of these sources will allow us to test the hypothesis of population mixing.

In Fig. 10, we show the CEDF of each group as a function of time. In this figure, the energy of the sources was computed with respect to their parent group. As can be observed, the most significant change in energy in the following 20 Myr will occur in the core of Coma Ber followed by that of Mecayotl 1. In the case of Coma Ber’s core, this change points towards a more steep CEDF function, which indicates that, as time goes by, the group’s sources will homogenise their energies, and as a consequence, the fraction of gravitationally bound stars (energies < 0) will diminish. In the case of Mecayotl 1, we observe that its CEDF has a similar but less pronounced behaviour as that of Coma Ber’s core given that it has virtually no energetically bound sources. The slope of its CEDF will increase for near-zero energy sources and will remain similar for higher energy sources, which will result in an overall increase in the fraction of sources for any given energy value.

Fig. 10. Time evolution of each physical group’s cumulative energy distribution function. For visual aid, the vertical grey dashed lines show the zero energy.

We use the present-day values of the energy of each sample-particle with respect to the other groups as an additional validation to our membership classification. We find that none of our candidate members has present-day energies lower than those with respect to its parent group. It shows that our statistical classification methodology is robust.

We observe that at all time stamps of our orbit integration, the energies of each sample-particle with respect to the other groups are higher than those with respect to its parent group. Thus indicating that at the flyby distances the gravitational potentials of the groups are not strong enough to capture members from another group. Therefore, we conclude that, under our current assumptions and observational uncertainties, we find no evidence of population mixing.

5. Discussion

In this section we compare the results presented in the previous sections with the most recent ones from the literature and the relevant theories that describe them. Table 4 summarises the general properties of the stellar systems analysed here. The distance and total mass values correspond to the mean and standard deviation computed from the sample of candidate members.

5.1. Physical groups

We found that Coma Ber is composed of two subgroups: the core and the tails, with the tails being more populous than the core (166 vs 136) but with a lower total mass (74 ± 4 M_⊙ vs 87 ± 4 M_⊙, see also Fig. 6). This behaviour is expected when the less massive stars gain energy, for example thanks to dynamical relaxation (e.g. two-body encounters) or tidal effects (e.g. a constant tidal field or tidal shocks), and thus, move towards larger cluster radii until eventually they get expelled through the Lagrangian points L₁ or L₂, and end up in the leading or trailing tidal tails, respectively (see, for example, Küpper et al. 2012 and references therein).

We found that the moving group simultaneously discovered by Fürnkranz et al. (2019) and Tang et al. (2019), and called Group X by the latter, is composed of two independent physical groups. One corresponds to the candidate open cluster Latyshev 2, while the other is a sparse and kinematically loose stellar association that we call Mecayotl 1.

Latyshev 2 has a compact velocity dispersion (0.80 ± 0.08 km s⁻¹) similar to that of the Pleiades (0.8 ± 0.1 km s⁻¹, Galli et al. 2017) and the core of Coma Ber (0.89 ± 0.1 km s⁻¹). Furthermore, its CEDF and energy-distance distribution (see Figs. 7 and 8) resemble those of Coma Ber’s core. The previous results allow us to confirm that Latyshev 2 is an open cluster as originally proposed by Latyshev (1977). We nonetheless notice that its large spatial dispersion (σ_XYZ = 21.04 ± 1.29 pc) and the negligible fraction of energy-bound sources (≲5%, see Fig. 7) indicate that it is in a disrupted state, which may have resulted from past encounters with molecular clouds or other stellar systems.

Mecayotl 1 has a large velocity dispersion (1.85 ± 0.18 km s⁻¹), similar to that of the βPic stellar association (1.73 ± 0.3 km s⁻¹Miret-Roig et al. 2020), and even larger than that of the tidal tails of Coma Ber (1.43 ± 0.1 km s⁻¹). Furthermore, its CEDF and energy-distance distribution (see Figs. 7 and 8) have larger dispersions than those of the tails of Coma Ber. Concerning its spatial dispersion, 23.5 ± 1.43 pc, we observe that it is three times larger than the core of Coma Ber (6.66 ± 0.64 pc), slightly larger than that of Latyshev 2 (21.04 ± 1.29 pc), but smaller than that of Coma Ber’s tails (34.34 ± 1.89 pc). Therefore, we conclude that Mecayotl 1 resembles a stellar association rather than an open cluster. We notice that its large (space, velocity, and energy) dispersions may be primordial, the result of dynamical encounters that disrupted the system, or artefacts of averaging unidentified and unrelated substructures.

Concerning this latter point, the energy vs distance diagram of Mecayotl 1 shows over-densities (see Fig. 8) that could be related to possible unidentified substructures. However, the present observational uncertainties do not provide evidence supporting these substructures. Therefore we conclude that Mecayotl 1 is one single physical group.

5.2. Membership

We now compare our membership lists with previous ones from the literature, in particular with the most recent ones obtained based on Gaia data. Figures 11–13 compare the astrometry and photometry of our candidate members with those previously known in the literature (i.e. Tang et al. 2019; Fürnkranz et al. 2019; He et al. 2022; and Newton et al. 2022; these latter authors split their sample into those having velocity offsets with respect to TOI-2048 larger and smaller than 2 km s⁻¹). As can be observed, our candidate members extend farther in the sky and further in photometric magnitudes than in previous searches from the literature. In particular, Fig. 11 shows that our search for candidate members extends to half the northern Galactic hemisphere (b > 30° without parallax cuts). The latter represents a considerable improvement with respect to previous works that restrict their searches in space volume, for example, to radii of 85 pc (Tang et al. 2019), 75 pc (Fürnkranz et al. 2019) or 25 pc (Newton et al. 2022).

Fig. 11. Galactic coordinates of our candidate members (olive stars) and those from the literature. The candidate members of Latyshev 2 and Mecayotl 1 are located to the left bottom of the image, while those of Coma Ber are at the centre-right.

Fig. 12. Gaia astrometry and photometry of the Coma Ber candidate members from the literature and this work.

Fig. 13. Gaia astrometry and photometry of Latyshev 2 (+Mecayotl 1) candidate members from the literature and this work.

We compare our Latyshev 2 list of candidate members with those proposed by Latyshev (1977). We recovered only three¹⁶ out of the seven candidate members originally proposed in that work. These three sources are shown as green diamonds in Figs. 11 and 13. The remaining four candidate members were discarded due to their small parallaxes (< 7 mas) and highly discrepant proper motions.

Table 5 shows the number of Coma Ber’s candidate members from the literature that our membership algorithm recovers and rejects. Similarly, Table 6 shows the number of recovered and rejected candidate members from the literature, but this time for what is known as Group X. In this latter case, we aggregated our lists of candidate members for Latyshev 2 and Mecayotl 1 into a single one corresponding to what is known as Group X in the literature works.

As can be observed in Tables 5 and 6, our novel membership methodology recovers between 42% to 305%, and between 52% and 275% more candidate members in Coma Ber and Group X (Latyshev 2 + Mecayotl 1), respectively, than the previous works from the literature. In particular, we increase the number of Coma Ber and Group X (Latyshev 2 + Mecayotl 1) candidate members by 42% and 87% with respect to those of Fürnkranz et al. (2019), by 53% and 52% with respect to Tang et al. (2019), and by 305% and 275% with respect to He et al. (2022), respectively. The work of Newton et al. (2022) focused only on Group X and split their list of candidate members into those having velocity offsets larger and smaller than 2 km s⁻¹ with respect to TOI-2048. Comparing our list of Group X (Latyshev 2 + Mecayotl 1) candidate members with those of the previous authors, we increase their numbers by 60% and reject 11% of their Δv_tan < 2 km s⁻¹ sample and 90% of their Δv_tan > 2 km s⁻¹ sample.

In Coma Ber, our new candidate members are primarily located in the farthest regions of the tidal tails and in the faintest magnitudes of the cluster photometric sequence (see Fig. 12). In the case of Group X (Latyshev 2 + Mecayotl 1), our new candidate members are primarily located in the regions of large parallaxes (ϖ > 11 mas) and faintest magnitudes (see Fig. 13).

Due to the large recovery rates (> 95%) and low contamination rates (≤5%) of our classifiers (see Sect. 4.1), we expect that our lists of candidate members have larger purity than those from the literature. In particular, our recovery and contamination rates are better than those reported by Tang et al. (2019 90% of recovery rate and 5-6% of contamination). Unfortunately, Fürnkranz et al. (2019) only reported an estimated contamination rate of a few per cent, and the most recent works of Newton et al. (2022) and He et al. (2022) do not report recovery or contamination rates.

To provide a quantitative comparison of our classifiers’ quality with respect to those from the literature, we use the numbers of recovered and rejected candidate members from those works together with our new candidate members. We compute the literature classifier’s recovery, contamination, and missing rates as the fractions of recovered, rejected and missed sources, with respect to our lists of candidate members, respectively, and assuming that these latter are complete. These quality indicators are shown in Table 7. As previously mentioned, the classifier’s quality indicators of Tang et al. (2019) agree well with the ones they report. The 4% to 10% contamination rates that we measure for Fürnkranz et al. (2019) classifiers are larger than their estimated few per cent contamination rate. In the case of Newton et al. (2022), we observe that the quality of their classifier for the sample with Δv_tan < 2 km s⁻¹ is similar to that of Fürnkranz et al. (2019), whereas the sample with Δv_tan > 2 km s⁻¹ has worse quality indicators than those of a random classifier.

From our previous estimates, we notice the following points. First, the StarGo classifier of Tang et al. (2019) misses only one-third of our candidate members. This classifier performs better in the more sparse case of Group X (Latyshev 2 + Mecayotl 1) than in the case of Coma Ber. Second, the classifier of Fürnkranz et al. (2019) performs better in Coma Ber than in Group X (Latyshev 2 + Mecayotl 1), which is not surprising given that this classifier was designed to identify extended structures like tidal tails and stellar streams (see Meingast et al. 2019). Third, the classifier of He et al. (2022) works better in Coma Ber than in Latyshev 2 and Mecayotl 1, which comes as no surprise since it was tuned for more compact open clusters. Fourth, the DBSCAN classifier of He et al. (2022) and the FindFriends query method used by Newton et al. (2022) miss up to two-thirds of our candidate members.

We highlight that the different performances that the literature classifiers have when applied to stellar systems in distinct dynamical states point towards intrinsic differences in the machine-learning algorithms they use and their fine-tuned and manually chosen parameters, particularly in the number of neurons used by the self-organising map of Tang et al. (2019) and the ϵ and minPts DBSCAN parameters used by Fürnkranz et al. (2019) or He et al. (2022). Contrary to the previous methods, our novel membership methodology offers clear physical and statistical interpretability of both the membership algorithm and its resulting phase-space model. Furthermore, our new method returns similar contamination and recovery rates (see Tables 1 and 2) in spite of the dynamical state of the stellar system to which it is applied. The latter are examples of the advantages that physical forward modelling offers over classical machine-learning algorithms imported from other domains. In spite of the previous advantages, our methodology still faces difficulties, particularly in recovering the tail’s most distant and low-mass members (see discussion in Sect. 5.4). Future improvements of the tail’s phase-space model are warranted.

5.3. Phase-space structure

Our comprehensive phase-space modelling allows us, on the one hand, to reveal that Latyshev 2 and Mecayotl 1 are two physically distinct stellar systems and, on the other hand, to extend the tidal tails of Coma Ber up to distances of 95 pc from the cluster centre. We now discuss the position and velocity distributions of our candidate members.

In Fig. 14, we show the distributions of offsets, both in distance and in velocity, that our candidate members have with respect to the centre of their parent group (histogram). These offsets are shown in Galactic heliocentric coordinates, for which we follow the PyGaia reference frame definitions and transformations. The figure also shows the groups’ Gaussian models obtained in Sect. 4.2 after being transformed from the ICRS to the Galactic reference frame. The uncertainties in both the individual stars and the groups’ parameters were incorporated by taking 100 samples from their posterior distributions. This figure allows us to confirm that although our model building blocks are Gaussian, the inferred position and velocity distributions of the individual sources are not necessarily Gaussian, as can be observed by the long tails in both positions and velocities of Coma Ber’s tails. For comparison purposes, the figures depict the group’s total dispersions (σ_XYZ and σ_UVW from Table 3) with black lines. In the case of Coma Ber’s tail, we also add the dispersion of Coma Ber’s core in grey lines.

Fig. 14. Distribution of position (top) and velocity (bottom) offsets with respect to the group centre in Galactic heliocentric coordinates. The histograms and smooth lines show samples from the posterior distribution of the source- and group-level parameters, respectively.

The 3D spatial distribution of Latyshev 2, Mecayotl 1, and the core and tails of Coma Ber show smaller dispersion in the Galactic Z direction than in the other two. This phenomenon has been observed in several open clusters and galactic structures (see, for example, Moranta et al. 2022; Meingast et al. 2019; Kraus et al. 2014), which indicates that it results from the interaction with the Galactic potential. Furthermore, the tidal tails in the Galactic Z direction have the same typical size as the cluster core, thus indicating that the escaped stars remain collimated in this direction. On the other hand, the Galactic X and Y distributions show the following properties. In the case of Latyshev 2, the X and Y distributions appear Gaussian, although with wider tails. Similarly, the X and Y distributions of Mecayotl 1 have wider tails than a Gaussian distribution, but in this case, they present the over-densities described in Sect. 4.4, particularly seen in the Y distribution. In Coma Ber, the core has two over-densities in the Y distribution roughly corresponding with the spatial dispersion and tidal radius of the core (6.66 ± 0.64 pc and 6.32 ± 0.22 pc, respectively, see Table 3). In the X distribution, these over-densities are shifted out at −9 pc and 11 pc. We interpret these over-densities as the result of the accumulation of escaping stars at the Lagrangian points, where they are most likely to escape. The tails of Coma Ber also have two over-densities in both X and Y located at −14 pc and 9 pc, which we interpret as the accumulation of escaped stars still remaining in the cluster’s vicinity due to their low escaping velocities. We notice that, as pointed out by Küpper et al. (2012), the tidal radius (r_t, see Table 3 and Sect. 4.2) is time-dependent for clusters with eccentric orbits and can significantly change between peri- and apogalacticon. Therefore, there can be previously escaped stars within the current tidal radius and future-bound stars outside of it. This breathing of the tidal radius will result in the accumulation of stars in the Lagrangian regions.

The 3D velocity distributions of Latyshev 2, Mecayotl 1, and the tails of Coma Bar have fewer features than their 3D spatial counterparts. On the other hand, the velocity distribution of Coma Ber’s core has symmetric over-densities located at the wings of the V distribution. Given that these over-densities are only present in the Galactic plane and in the V direction, we interpret them as resulting from either potential escapers or misclassified sources belonging to the tails. The Coma Ber’s tail velocity dispersion in the U direction lacks the central peak present in the other directions (V and W) and thus resembles a flat-top Gaussian. We interpret the lack of this central peak as resulting from escapers gradually accelerating in this direction and thus effectively moving them away from the central zero-velocity peak. Although significant only at the one-sigma level, the observed expansion of the tails (0.7 ± 0.7 km s⁻¹; see Table 3) also supports this scenario. Furthermore, while the tidal tail velocity dispersion in the U and V directions increases with respect to those of the core, the velocity dispersion in the W direction remains similar. This effect indicates that the main Galactic forces acting upon the tails remain in the Galactic plane.

Our results show that the stellar systems analysed here have non-zero rotation velocities (see Table 3), with similar values to those observed in open clusters and star-forming regions. For example, the rotation velocities of Coma Ber’s core and Latyshev 2 are consistent with the Praesepe value of 0.4 km s⁻¹ reported by Loktin & Popov (2020), and the rotation velocity of Mecayotl 1, 1.32 ± 0.74 km s⁻¹, is also consistent with the 1.5 ± 0.1 km s⁻¹ reported by Galli et al. (2019) for the Taurus complex. In addition, the total rotation of Coma Ber (1.13 ± 0.75 km s⁻¹, computed by adding in quadrature the core and tail values) is similar to that of Mecayotl 1 and compatible with that of the Taurus complex. Although our rotational velocities are non-zero and compatible with those from other stellar systems in the solar neighbourhood, their significance still not reaches the two-sigma level (i.e. 95%) needed to discard the null hypothesis of no rotational velocity. Future work is still needed both in the modelisation and the observational sides.

5.4. Tail overdensity discussion

We investigate the radial distribution of Coma Ber’s candidate members in the two tails. To identify the candidate members in the leading and trailing tails we use the inner product of their radial vector (with respect to the cluster centre) with the vector of the cluster’s orbit (see Sect. 4.5). If the resulting inner product is positive we associate the candidate to the leading tail, otherwise to the trailing tail. We identify 101 and 65 candidate members in the leading and trailing tails, respectively. This represents a leading-to-trailing tail number ratio of 1.55 ± 0.34¹⁷. A similarly large ratio of 1.76 was also reported for the Hyades open cluster (Jerabkova et al. 2021).

Figure 15 shows the distribution of the distance to the cluster centre for the Coma Ber’s leading and trailing tails. In this figure, the orange solid lines depict the kernel density estimate of samples from each candidate member’s radial distance to the cluster centre, as computed from samples of their posterior distributions of the Cartesian spatial coordinates (X, Y and Z). The black dots and grey bars depict the mean and standard deviations obtained after binning the mean radial distance of the candidate members in bins of 2 pc width. This width represents the maximum uncertainty in the radial distance determination (the typical one being 0.5 pc). As can be observed, the extensions of the tails differ, with the leading one reaching up to 95 pc while the trailing one only 80 pc. The radial distributions of both tails show over-densities. We investigate the statistical significance of these over-densities by comparing them to a model of the radial distribution. For such, we use a Gamma distribution whose parameters we infer from 100 bootstrap realisations of each tail candidate members. This figure also depicts, at each radial distance, the median value of the model (black solid lines) as well as one-sigma (blue dashed lines), two-sigma (orange dashed lines) and three-sigma (green dashed lines) confidence intervals. As can be observed, only two radial bins in each tail have densities (median minus standard deviation) exceeding the two-sigma (orange arrows) and three-sigma (green arrows) levels of the distribution. In the leading tail these correspond to the over-densities located at 35 ± 1 pc and 45 ± 1 pc, while in the trailing tail, these correspond to the over-densities at 29 ± 1 pc and 33 ± 1 pc. We notice that these over-densities are still significantly detected for bin-widths of up to ∼4 pc and up to ∼3 pc in the leading and trailing tails, respectively.

Fig. 15. Distance-to-centre distribution of Coma Ber’s leading and trailing tail members. The orange lines depict the kernel density estimate of samples from the radial distance while the black dots and grey bars depict radial bins mean and standard deviation of 100 bootstrap realisations from the tail’s members. A Gamma distribution was fitted to these bootstrap realisation and its median value (black solid lines) and one-sigma (blue dashed lines), two-sigma (orange dashed lines) and three-sigma (green dashed lines) confidence intervals are also plotted. The vertical arrows indicate the radial bin with densities exceeding the two-sigma (orange arrows) and three-sigma (green arrows) levels.

Clumps in the tidal tails are expected to appear due to gravitational shocks (e.g. due to encounters with molecular clouds or passages through the Galactic disk) or by the epicyclic motion of stars escaping from a cluster in a circular galactic orbit and in a static tidal field (see Küpper et al. 2008 and references therein). The latter are known as epicyclic over-densities and correspond to the places where escaped stars spend more time in their epicyclic motion away from the cluster (Küpper et al. 2008). These over-densities are expected to depend on several factors, like the mass of the cluster, the eccentricity of its orbit, or the velocity of the escaped stars. Nonetheless, the location of these over-densities can be roughly approximated to be y_C ≃ 3πr_t (see Eq. (4) of Küpper et al. 2012) with r_t as the tidal radius. As mentioned before, this tidal radius is expected to vary as a function of time for clusters in excentric orbits (see for example, Eq. (2) of the aforementioned authors).

Assuming that the observed over-densities correspond to the epicyclic ones, their distances would result in tidal radii of 4 ± 1 pc and 5 ± 1 pc for the leading tail, and of 3 ± 1 pc and 4 ± 1 pc for the trailing tail. These values are compatible within their uncertainties and smaller than the present-day value of the tidal radius, 6.32 ± 0.21 pc (see Table 3). A possible explanation for these small tidal radii is that their associated over-densities were produced at the perigalacticon and thus the cluster has an eccentric orbit. Future work will be needed to clarify if the over-densities found here are indeed of epicyclic origin.

5.5. Luminosity and mass distributions

Concerning the luminosity and mass distributions (see Figs. 5 and 6), we observe that in the case of Coma Ber, they show clear differences between the core and the tails, in particular at the low-mass domain. We interpret these differences as resulting from the ejection and evaporation of low-mass stars. We notice that the fact that our luminosity and mass distributions at the faint magnitude bin represent only upper limits to the true distributions has no impact on the previous conclusion. On the other hand, Latyshev 2 and Mecayotl 1 have similar luminosity and mass distributions with peak masses of ∼0.33 M_⊙ (0.36 M_⊙ for PARSEC and 0.30 M_⊙ for BT-Settl and PB) and ∼0.25 M_⊙ (0.28 M_⊙ for PARSEC and 0.22 M_⊙ for BT-Settl and PB), respectively. This relative shift towards higher masses in the peak of Latyshev 2 further supports our older age estimate for this system with respect to that of Mecayotl 1.

The mass distribution of Coma Ber’s tails also shows a lack of low-mass stars, ≲0.2 M_⊙, (see Fig. 6) for what is expected from the initial mass distribution of Chabrier (2005). Given that this lack is well within the Gaia completeness limits, we conclude that our membership algorithm still misses low-mass stars in the most distant edges of the tidal tails. Identifying this distant population of low-mass stars and brown dwarfs will be a significant future challenge given the low contrast and large uncertainties of these faint sources. Future steps will be taken to improve our current forward model of the cluster’s tails.

5.6. Energy distribution

In Sect. 4.4, we presented the energy distributions of Coma Ber, Latyshev 2 and Mecayotl 1 and described their observed properties. To the best of our knowledge, this is the first time that such observationally derived energy distributions are presented for open clusters in the literature. We now analyse these distributions and compare them with each other and theoretical models.

Classically, the phase-space distribution of open clusters has been described with King’s profile (King 1962) and model (King 1966). These were derived based on either the radial density profiles of globular clusters (King 1962) or the assumption that their velocity distribution is Gaussian and reaches a zero value at the cluster surface (King 1966). Küpper et al. (2010) provide a clear disambiguation between the latter two. The literature has shown that King (1962) profile faces difficulties when fitting the radial density profiles of open clusters (e.g. Zhong et al. 2022; Tarricq et al. 2022; Olivares et al. 2018a), particularly when these have halos or tails. Küpper et al. (2010) analysed the profiles of dynamically evolved simulated clusters and found that although the King (1962) profile is a reasonable tool to measure the Jacobi radius it fails to capture the distribution of escaped stars in the form of the tidal debris. Concerning the King (1966) model, the literature has shown that it fails to describe the observed properties of globular clusters (see, for example, Claydon et al. 2019) and several works have created improved models of these stellar systems (see, for example, Webb et al. 2023 and references therein). However, little attention has been paid to the modelling of open clusters despite the challenges they pose (see, Sect. 9.1 of King 2008), and King’s profile and model continue to be widely used by the open cluster community.

To compare our observational energy distribution with that expected from the King (1966) model, we simulated six stellar clusters using the McLuster code (Küpper et al. 2011). Each synthetic cluster has a total mass of 160 M_⊙ and between 280 and 330 stars, with phase-space coordinates and masses sampled randomly from the King (1966) and Kroupa (2001) distributions, respectively. These clusters have half-mass radii of 1 pc and concentration parameters, W₀, varying between 1 and 12. We choose the total mass and half-mass radii of these synthetic clusters to resemble those of Come Ber; the chosen W₀ values cover the full spectrum of this parameter. We compute the CEDFs of these synthetic clusters in the same way as our observed ones. Figure 16 shows the empirical CEDFs of ComaBer’s core (orange lines) together with the CEDFs of these synthetic clusters colour-coded with the value of their W₀ parameter. The solid lines show the original CEDFs and the dashed lines show the CEDFs that result from subtracting a total mass of 80 M_⊙. The latter is similar to the total mass estimated for the tidal tails (see Table 4). As can be observed, the original CEDFs of these simulated clusters show that 90% of its stars have negative energies, which is in clear disagreement with our observational CEDF of Coma Ber’s core. On the contrary, the CEDFs resulting from subtracting 50% of the cluster mass produce CEDFs resembling the observed ones. Therefore, we conclude that the King (1966) model fails to capture the energy distribution of open clusters in late dynamical stages in which they have lost considerable fractions of their initial mass, as already shown by dynamical simulations (Küpper et al. 2010). Future work will be needed to analyse the energy distributions of other open clusters and stellar associations at different stages of their dynamical evolution and to develop more realistic models of their phase-space and energy distributions.

Fig. 16. CEDFs of synthetic clusters with positions and velocities following the King (1966) profile and mimicking the properties of Coma Ber (see text). The orange lines depict the empirical CEDFs of Coma Ber’s core (see Fig. 7).

We observe that the main characteristic that distinguishes our observational CEDFs is their slope. On the one hand, low positive slope values result in dispersed energy distribution, as in the case of Mecayotl 1 and the tails of Coma Ber. On the other hand, large positive slopes result in concentrated energy distributions, as in the core of Coma Ber and Latyshev 2, where most of the sources have small but positive energies. The CEDFs of these two open clusters agree with previous results from the literature, in which energetically unbound stars remain within the cluster for several dynamical time scales, and are preferentially located at the cluster’s edge (as mentioned in Sect. 5.3), where the dynamical time-scale is of the order of the orbital time-scale of the cluster (Küpper et al. 2012).

5.7. Trace-forward

In Sect. 4.5, we estimated the times at which Mecayotl 1 and Latyshev 2 will have their closest distance to Coma Ber and between themselves. These flyby times have a precision between 5% and 16%, which represents a major improvement with respect to the 37% precision of Tang et al. (2019). The latter reported a single flyby time at 10–16 Myr, while Fürnkranz et al. (2019) also reported a single flyby time of 13 Myr but with no associated uncertainty.

We stress the fact that our flyby times and test for population mixing are computed under the assumption that the clusters and their stars move as free particles in the gravitational potential of the Milky Way. However, this assumption neglects both the gravitational potential of the clusters as well as possible encounters between the cluster’s members. Numerical N-body simulations will be needed to overcome this assumption. However, these simulations are beyond the scope of the present work and will be carried out in a following paper of this series.

In an attempt to disentangle the possible origin of Mecayotl 1 and Latyshev 2, we also performed a trace-back orbit integration down to −50 Myr and found no other past encounters among these three systems. In addition, the distances of Mecayotl 1 and Latyshev 2 to Coma Ber rapidly diverge in the first −10 Myr, thus confirming that Latyshev 2 and Mecayotl 1 are unrelated to Coma Ber. Furthermore, the distance between Latyshev 2 and Mecayotl 1 also diverges in time but at a slower rate, reaching separations of 70 pc in the past 50 Myr. This result suggests that if the current approach between Latyshev 2 and Mecayotl 1 is the origin of the large energy dispersion observed in the latter (see Fig. 8), then numerical N-body trace-back simulations will be needed to clarify this scenario. However, the unravelling of the origin of these two stellar systems still awaits the arrival of more precise radial velocity measurements, precise age determinations, and numerical trace-back simulations.

6. Conclusions

Using public data from 2MASS, APOGEE, Gaia, and PanSTARRS, in combination with our own INT-IDS observations, we studied the well-known Coma Ber open cluster and its neighbour moving group called Group X by Tang et al. (2019). Our novel forward model membership method confirmed, with more than 95% confidence level, that Group X is composed of two independent physical structures as Tang et al. (2019) originally suggested. These structures correspond to the disrupted open cluster Latyshev 2 (an open cluster proposed by Latyshev 1977) and a newly identified stellar association that we call Mecayotl 1.

Our novel membership method unravels ∼50% more candidate members than recent studies from the literature. Moreover, thanks to carefully crafted synthetic data sets that mimic the real stellar systems and have the uncertainty properties of the Gaia data, we measured recovery and contamination rates ≳90% and ∼5%, respectively.

Based on the most complete samples of candidate members for Coma Ber, Latyshev 2, and Mecayotl 1 to date, we analyse their phase space, mass, and energy distributions and we integrated their Galactic orbits 20 Myr into the future and 50 Myr into the past. From the previous set of results, we draw the following conclusions.

The tidal tails of Coma Ber are dynamically cold and formed through mass segregation and evaporation. The leading tail is more populous (101 sources) and elongated (95 pc) than the trailing one (65 sources, 80 pc). Although we identified tidal tails’ members at up to 95 pc from the cluster centre, their mass distribution indicates that our membership method still misses low-mass stars within the Gaia completeness limit.

Latyshev 2 is an open cluster, as originally proposed by Latyshev (1977), that is in a disrupted stage. Its mass (92  ±  5 M_⊙), velocity dispersion (0.80 ± 0.08 km s⁻¹), and energy distribution (see Fig. 7) resemble those of Coma Ber’s core, although with a larger population (186 candidate members). Its relatively young age (∼400 Myr) and proximity to the Sun make it an excellent target for studies of cluster evolution and searches for ultra-cool dwarfs and exoplanets.

Mecayotl 1 is a relatively populous (146 candidate members), young (400-600 Myr), and sparse stellar association. It has a large velocity dispersion (1.85 ± 0.18 km s⁻¹), a total mass of 75 ± 5 M_⊙, and an energy distribution with more dispersion than that of the tails of Coma Ber.

Our trace-forward analysis shows that Latyshev 2 and Mecayotl 1 will experience flyby encounters with Coma Ber in 11.3 ± 0.5 Myr and 14.0 ± 0.6 Myr, respectively, and between them in 8.1 ± 1.3 Myr. The precision in the time of these encounters (5% and 16%) represents a considerable improvement with respect to the 37% precision of Tang et al. (2019). Although the encounters’ distances will be smaller to the groups’ size, we find no evidence for the mixing of the groups’ populations. On the other hand, our trace-back analysis shows that these two stellar systems are most likely unrelated in origin. Despite our previous results, the dynamical analysis of Coma Ber, Latyshev 2, and Mecayotl 1 still requires further methodological improvements, such as better phase-space models for the tidal tails, trace-back and trace-forward numerical N-body simulations, together with precise and accurate absolute age determinations.

Finally, we notice that thanks to their proximity to the Sun, their high Galactic latitude, and their different dynamical stages, the three stellar systems analysed here will be fundamental for the search and characterisation of ultra-cool dwarfs in the Euclid (Laureijs et al. 2011) era.

Acknowledgments

We thank the anonymous referee for the useful comments that improved the quality of this work. JO acknowledges financial support from “Ayudas para contratos postdoctorales de investigación UNED 2021”. JO, NL, VB, and ELM acknowledge support from the Agencia Estatal de Investigación del Ministerio de Ciencia e Innovación (AEI-MCINN) under grant PID2019-109522GB-C53. NL and MŽ acknowledge support from the Consejería de Economía, Conocimiento y Empleo del Gobierno de Canarias and the European Regional Development Fund (ERDF) under grant with reference PROID2020010052. P.A.B. Galli acknowledges financial support from São Paulo Research Foundation (FAPESP) under grants 2020/12518-8 and 2021/11778-9. Co-funded by the European Union (ERC, SUBSTELLAR, project number 101054354). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. We acknowledge Lucía Suarez, Olga Zamora, George Hume, and the rest of the INT staff for their support during the observations, as well as their kind help afterwards. This article is based on observations made in the Observatorios de Canarias del IAC with the Isaac Newton telescope, under Spanish Service Time, operated on the island of La Palma by the Isaac Newton Group of Telescopes in the Observatorio Roque de Los Muchachos. The Isaac Newton Telescope is operated on the island of La Palma by the Isaac Newton Group of Telescopes in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias. 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 made use of computing time available on the high-performance computing systems at the Instituto de Astrofísica de Canarias. The author thankfully acknowledges the technical expertise and assistance provided by the Spanish Supercomputing Network (Red Española de Supercomputacion), as well as the computer resources used: the Deimos-Diva Supercomputer, located at the Instituto de Astrofísica de Canarias. 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. 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. Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. 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. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. We thank Diego Godoy-Rivera for the useful discussions that took place during the Cool Stars 21 and greatly improved the quality of this work.

References

Appendix A: INT-IDS observations

We use INT-IDS spectrograph to acquire the spectra of 36 targets and four radial velocity standards in visitor mode on the nights of April 12th to April 18th 2022. We selected our IDS targets from our list of candidate members (see Sect. 4.1) based on two criteria: poor quality (uncertainties > 10km s⁻¹) or non-existing radial velocity measurements in Gaia EDR3, APOGEE or SIMBAD, and V band magnitude brighter than 15.5 mag. The spectra were collected with the R1200R grating centred at 6500Å, covering the 5300–7100Å wavelength range and yielding a spectral resolving power of 6731. We used the IDS default detector, RED+2 CCD, which has pixels of 15.0 micron, a plate scale of 66.7 pixels/mm and 0.44″/pixel, and 2200 unvignetted pixels.

Tables A.1 and A.2 show the identifiers and characteristics of our 36 targets and four radial velocity standards, respectively. As can be observed, our targets have radial velocity uncertainties with values ranging from 5 to 40 km s⁻¹, with a typical value of ∼10 km s⁻¹. We notice that the radial velocities that we measured for our standard stars depart from the values reported by Soubiran et al. (Soubiran et al.) shown in the last two columns of Table A.2). For this reason, we computed radial velocity zero-points by median combining the standards of each night assuming Soubiran et al. (2013) values. These zero-points were used to correct our target’s radial velocity, and its uncertainty was added in quadrature.

A.1. Radial velocity comparison

With the aim of investigating the quality of the radial velocities measured here, we compare them with those reported in APOGEE, SIMBAD, Gaia DR2 and Gaia DR3 for the candidate members of Coma Ber, Latyshev 2 and Mecayotl 1, and using as reference the Gaia DR3 values due to their precision and large completeness. Figure A.1 shows one-to-one scatter plots of the subsets of candidate members with radial velocity measurements in Gaia DR3 and each of the other surveys. In the case of our IDS measurements (red dots), the figure shows the radial velocities of our 26 targets with Gaia DR3 values as well. As can be observed, our IDS radial velocity measurements compare well with the Gaia DR3 values and with those of APOGEE, SIMBAD, and Gaia DR2, although with larger uncertainties. Nonetheless, there are four IDS outliers that have a similar dispersion as the five SIMBAD outliers. We notice that in all these cases, the Gaia DR3 values were used instead.

Fig. A.1. Comparison of the radial velocities measured here (red dots) and reported in APOGEE (green), SIMBAD (blue) and Gaia DR2 (orange) with respect to those of Gaia DR3 (horizontal axis). The values shown correspond to the subsets of Coma Ber, Latyshev 2 and Mecayotl 1 candidate members having radial velocities from both Gaia DR3 and the other source. The grey dashed line shows the identity relation.

Appendix B: Radial velocity zero-points

We evaluate the consistency of the radial velocity measurements originating from the different input catalogues using as reference the Gaia DR3 radial velocities. We compute zero-point offsets for APOGEE, SIMBAD, and IDS using the subsets of our candidate members simultaneously having radial velocities from Gaia DR3 and the other survey. The zero-point offset of SIMBAD, APOGEE and IDS are 0.02 ± 1.46 km s⁻¹, 0.44 ± 0.81 km s⁻¹, and 0.32 ± 8.86 km s⁻¹, respectively. We notice that given their uncertainties (computed as standard errors) the latter values are compatible with zero and thus we assume they can be neglected.

Katz et al.(2023, Sect. 6.1) compared the Gaia DR3 radial velocities to those measured by other radial velocity surveys including APOGEE. In their comparison, these authors used two groups of stars: dwarf-turn-off stars and red giants and clump stars and found that the median radial velocity residuals (subtracting the Gaia radial velocities) show a trend in magnitude grvs_mag. They recommend correcting the Gaia radial velocities of stars cooler than 8500 K and fainter than 11 mag using the relation given by their Eq. 5. Given that the two samples of stars that those authors use are not representative of our main sequence dwarf candidate members, we explored possible trends in the residuals of the APOGEE, SIMBAD, and IDS radial velocities of our candidate members as functions of grvs_mag. Figure B.1 shows the radial velocity residuals of APOGEE, SIMBAD, and IDS with respect to Gaia DR3. As can be observed there is no clear trend in magnitude except by an increase of the dispersion for sources fainter than 11 mag. In addition, the figure also shows the median radial velocity residual (dashed lines) and its standard error (shaded regions). As mentioned above, these median values are consistent with zero. Furthermore, they are more than four times smaller than the median value of Gaia DR3 radial velocity uncertainties of our candidate members (1.8 km s⁻¹) and thus, we assume that they can be neglected.

Fig. B.1. APOGEE, SIMBAD, and IDS radial velocities residuals with respect to Gaia DR3 as a function of grvs_mag. The dashed lines and shaded regions depict the median value and its standard error, respectively.

Appendix C: Extinction distributions inferred under a uniform prior: A_v ∈ [0, 10] mag

In a first attempt to infer the luminosity and mass distributions of our candidate members, we used a wide and uniform extinction prior, with A_v ∈ [0, 10] mag. Figure C.1 shows the A_v posterior distributions resulting from the use of this prior (grey lines) and the PARSEC, BT-Settl, and PB models (colour-coded lines). In addition, the black and blue dashed lines show the Gaia DR3 and Planck Collaboration Int. XLVIII (2016) values reported for our candidate members, respectively. The GaiaA_v values were obtained by transforming the Gaia DR3 azero_gspphot to A_v with the Cardelli et al. (1989) relation assuming Rv = 3.1. The Planck Collaboration Int. XLVIII (2016) extinction values were obtained by querying the reddening values with the dustmaps package (Green 2018) at the sky coordinates of our candidate members and transforming them to A_v values assuming Rv = 3.1.

Fig. C.1. Distributions of Av as inferred using a uniform prior (Av ∈ [0, 10] mag, grey lines) and different theoretical models (colour coded legend). The black and blue dashed lines show the Av distributions of our candidate members as reported by Gaia DR3 and Planck Collaboration Int. XLVIII (2016), respectively.

As can be observed in Fig. C.1, the PARSEC model results in distributions with modal values A_v < 0.1 mag for all groups. Instead, the BT-Settl and the joint PB models deliver bimodal distributions with two over-densities, one located at A_v​ ≲ ​0.1 mag that coincides with the mode of the PARSEC model, and another one at A_v​ ∼ ​0.6 − 0.9 mag. In the case of the core of Coma Ber, there is only one peak at A_v​ ≲ ​0.1 mag. The disagreement between the observed photometry and the BT-Settl models, which is visible in Fig. 4 at the colour interval G − RP ∈ [0.8, 1.2] mag where the BT-Settl isochrone is bluer than the PARSEC one, resulted in larger A_v values for this model and thus for the PB one. Nonetheless, we notice that the PARSEC models have been empirically corrected to reproduce the optical and near-infrared photometric sequences of low-mass stars in metal-poor Globular clusters and the solar-metallicity open clusters M67 and Praesepe (Chen et al. 2014).

We further investigated the origin of the previous discrepancy by comparing our inferred A_v distributions with those reported by Gaia DR3 and Planck Collaboration Int. XLVIII (2016). As explained above, these latter are also shown in Fig. C.1. As can be observed, the Gaia extinction distributions agree well with those inferred with the BT-Settl and PB models, having similar modal values and shape, except in the core of Coma Ber, where its decay is more pronounced. On the other hand, the extinction distributions delivered by Planck Collaboration Int. XLVIII (2016) are contained within A_v < 0.2 mag and have modal values that agree with those obtained with the PARSEC model: A_v < 0.1 mag.

It is well known that stellar models face problems to reproduce the observed photometry of low-mass stars (e.g. Chen et al. 2014; Feiden & Chaboyer 2012), with the empirical correction applied to the PARSEC models (Chen et al. 2014) offering a practical solution to this issue. Moreover, it is also known that Gaia DR3 extinction values tend to be overestimated for low-mass dwarfs (see Sect. 3.6.1 of Andrae et al. 2023). Due to the previous well known issues of the stellar models, in Sect. 3.3, we decided to infer stellar masses with an extinction uniform prior fixed to A_v ∈ [0.0 − 0.1] mag. We choose this interval because of the following reasons. First, it contains the peak (the mode) values inferred using the PARSEC models with the wider prior ( A_v ∈ [0, 10] mag). Second, it is the most probable extinction interval according to the Planck Collaboration Int. XLVIII (2016) values, which are based on the diffuse thermal dust emission, and therefore are independent of the assumptions and issues of to stellar models.

Appendix D: List of candidate members

All Tables

All Figures

	Fig. 1. Positions (top panel) and velocities (bottom panel) of Coma Ber candidate members. The 100 orange lines show the posterior samples from the inferred group-level parameters while the dots and error bar depict the mean and standard deviation, respectively, of the inferred positions and velocities.
In the text

	Fig. 2. Positions (top panel) and velocities (bottom panel) of the candidate members of Latyshev 2 and Mecayotl 1. Captions as in Fig. 1.
In the text

	Fig. 3. Distance distributions of Coma Ber, Latyshev 2 and Mecayotl 1.
In the text

	Fig. 4. Absolute colour-magnitude diagrams of the candidate members of Coma Ber, Latyshev 2 and Mecayotl 1. The lines depict the theoretical isochrone models with colour-coded ages.
In the text

	Fig. 5. Luminosity distributions inferred with different theoretical models. The grey line and grey region depict the uniform prior and the completeness limit of the Gaia data. This latter corresponds to apparent G mag ∼19 transformed to luminosity using the mean group distance and BT-Settl model.
In the text

	Fig. 6. Mass distributions inferred with different theoretical models. Caption as in Fig. 5. In this case, the grey line depicts the Chabrier (2005) mass prior.
In the text

	Fig. 7. Cumulative energy distribution functions. Each orange line shows the result obtained from a sample of the posterior distributions of mass, position and velocity. For visual aid, the vertical grey dashed lines show the zero energy.
In the text

	Fig. 8. Density of candidate members as a function of energy and distance to the group centre.
In the text

	Fig. 9. Galactic coordinates (R, Z, ϕ) of Coma Ber, Mecayotl 1 and Latyshev 2 as functions of time. To improve visibility, only the orbits of 20 samples of each stellar system are shown.
In the text

	Fig. 10. Time evolution of each physical group’s cumulative energy distribution function. For visual aid, the vertical grey dashed lines show the zero energy.
In the text

	Fig. 11. Galactic coordinates of our candidate members (olive stars) and those from the literature. The candidate members of Latyshev 2 and Mecayotl 1 are located to the left bottom of the image, while those of Coma Ber are at the centre-right.
In the text

	Fig. 12. Gaia astrometry and photometry of the Coma Ber candidate members from the literature and this work.
In the text

	Fig. 13. Gaia astrometry and photometry of Latyshev 2 (+Mecayotl 1) candidate members from the literature and this work.
In the text

	Fig. 14. Distribution of position (top) and velocity (bottom) offsets with respect to the group centre in Galactic heliocentric coordinates. The histograms and smooth lines show samples from the posterior distribution of the source- and group-level parameters, respectively.
In the text

Fig. 15. Distance-to-centre distribution of Coma Ber’s leading and trailing tail members. The orange lines depict the kernel density estimate of samples from the radial distance while the black dots and grey bars depict radial bins mean and standard deviation of 100 bootstrap realisations from the tail’s members. A Gamma distribution was fitted to these bootstrap realisation and its median value (black solid lines) and one-sigma (blue dashed lines), two-sigma (orange dashed lines) and three-sigma (green dashed lines) confidence intervals are also plotted. The vertical arrows indicate the radial bin with densities exceeding the two-sigma (orange arrows) and three-sigma (green arrows) levels.

In the text

	Fig. 16. CEDFs of synthetic clusters with positions and velocities following the King (1966) profile and mimicking the properties of Coma Ber (see text). The orange lines depict the empirical CEDFs of Coma Ber’s core (see Fig. 7).
In the text

Fig. A.1. Comparison of the radial velocities measured here (red dots) and reported in APOGEE (green), SIMBAD (blue) and Gaia DR2 (orange) with respect to those of Gaia DR3 (horizontal axis). The values shown correspond to the subsets of Coma Ber, Latyshev 2 and Mecayotl 1 candidate members having radial velocities from both Gaia DR3 and the other source. The grey dashed line shows the identity relation.

In the text

	Fig. B.1. APOGEE, SIMBAD, and IDS radial velocities residuals with respect to Gaia DR3 as a function of grvs_mag. The dashed lines and shaded regions depict the median value and its standard error, respectively.
In the text

	Fig. C.1. Distributions of Av as inferred using a uniform prior (Av ∈ [0, 10] mag, grey lines) and different theoretical models (colour coded legend). The black and blue dashed lines show the Av distributions of our candidate members as reported by Gaia DR3 and Planck Collaboration Int. XLVIII (2016), respectively.
In the text