Issue 
A&A
Volume 675, July 2023



Article Number  A28  
Number of page(s)  25  
Section  Galactic structure, stellar clusters and populations  
DOI  https://doi.org/10.1051/00046361/202244703  
Published online  29 June 2023 
The cosmic waltz of Coma Berenices and Latyshev 2 (Group X)
Membership, phasespace structure, mass, and energy distributions^{⋆}
^{1}
Departamento de Inteligencia Artificial, Universidad Nacional de Educación a Distancia (UNED), c/Juan del Rosal 16, 28040 Madrid, Spain
email: jolivares@dia.uned.es
^{2}
Instituto de Astrofísica de Canarias, 38205 La Laguna, Tenerife, Spain
^{3}
Departamento de Astrofísica, Universidad de La Laguna (ULL), 38205 La Laguna, Tenerife, Spain
^{4}
Núcleo de Astrofísica Teórica, Universidade Cidade de São Paulo, R. Galvão Bueno 868, Liberdade, 01506000 São Paulo, SP, Brazil
Received:
5
August
2022
Accepted:
17
April
2023
Context. Open clusters (OCs) are fundamental benchmarks where theories of star formation and stellar evolution can be tested and validated. Coma Berenices (Coma Ber) and Latyshev 2 (Group X) are the second and third OCs closest to the Sun, making them excellent targets to search for lowmass stars and ultracool dwarfs. In addition, this pair will experience a flyby in 10–16 Myr, making it a benchmark to test pair interactions of OCs.
Aims. We aim to analyse the membership, luminosity, mass, phasespace (i.e. positions and velocities), and energy distributions for Coma Ber and Latyshev 2 and test the hypothesis of the mixing of their populations at the encounter time.
Methods. We developed a new phasespace membership methodology and applied it to Gaia data. With the recovered members, we inferred the phasespace, luminosity, and mass distributions using publicly available Bayesian inference codes. Then, with a publicly available orbit integration code and members’ positions and velocities, we integrated their orbits 20 Myr into the future.
Results. In Coma Ber, we identified 302 candidate members distributed in the core and tidal tails. The tails are dynamically cold and asymmetrically populated. The stellar system called Group X is made of two structures: the disrupted OC Latyshev 2 (186 candidate members) and a loose stellar association called Mecayotl 1 (146 candidate members), and both of them will fly by Coma Ber in 11.3 ± 0.5 Myr and 14.0 ± 0.6 Myr, respectively, and each other in 8.1 ± 1.3 Myr.
Conclusions. We study the dynamical properties of the core and tails of Coma Ber and also confirm the existence of the OC Latyshev 2 and its neighbour stellar association Mecayotl 1. Although these three systems will experience encounters, we find no evidence supporting the mixing of their populations.
Key words: open clusters and associations: individual: Coma Berenices / stars: kinematics and dynamics / solar neighborhood / methods: statistical / astrometry / parallaxes
Full Table D.1 is only available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/vizbin/cat/J/A+A/675/A28
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
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, midage (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 masssegregated and elongated along the Galactic plane. In addition, they identify nine substellar 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 highresolution nearinfrared spectra from the SDSSIV 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 agemassrotation 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 planethost star TOI2048, 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^{−1}. 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 MELANGE2, containing 81 stars that share the 20 pc radius and 4 km s^{−1} volume around TIC 224606446. The latter is a nonmember of Group X that differs 5 km s^{−1} 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 nonphysical 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 highprecision 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, phasespace (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 reexamine 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, phasespace, and energy distributions. In this section, we also describe the traceforward 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^{1} the astrometry and photometry of ∼45 million Gaia DR3 sources with the following constraints: b > 30°, and proper motions within −400 < pmra/[mas yr^{−1}]< + 200 and −400 < pmdec/[mas yr^{−1}]< + 400. These proper motion filters are permissive enough to encompass most open clusters and starforming regions in the solar neighbourhood. We do not apply parallax cuts in the Gaia DR3 data to avoid biasing our sample of lowmass 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 crossmatch 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. Midresolution 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^{2} (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 INTIDS observations.
We reduced the spectra using the spextractor tool (Bouy et al., in prep.). Briefly, it corrects the raw images with the mediancombined 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 (CopperArgon and CopperNeon in the case of INTIDS) 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 crosscorrelating 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 (BlancoCuaresma 2019; BlancoCuaresma et al. 2014) by first removing telluric lines and then applying the crosscorrelation 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). Zeropoint 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^{−1}. Although these measurements are less precise than the median ones from APOGEE (0.12 km s^{−1}) or Gaia (1.8 km s^{−1}) 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 zeropoint offsets for the different surveys and conclude that given their small values and large uncertainties they can be neglected (see Appendix ). Given that our INTIDS observations were taken before the release of Gaia DR3, 26 of our INTIDS 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^{−1} that considerably increase the convergence time for our phasespace modelling methods (see Sect. 3.2). For this reason, we set radial velocity uncertainties smaller than 0.1 km s^{−1} 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 crossmatch 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 phasespace structure, mass and energy distributions of Coma Ber and Latyshev 2. Finally, the section concludes by presenting the traceforward 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^{3}, 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^{4} code (Olivares et al. 2020, see Sect. 3.2) to determine the cluster’s phasespace 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 phasespace model and sample 10^{6} synthetic stars whose coordinates we transform to the observed space (i.e. ra, dec, parallax, pmra, pmdec, and radial_velocity) with the PyGaia^{5} code.
– In the third step, we use Gaussian Mixture Models^{6} (hereafter GMMs) with up to 80 components to independently fit the density distribution of the synthetic cluster’s stars and that of an equalsize 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
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 opensource code Amasijo^{7} (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 opensource code^{8} (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^{3} sources each. In these clusters, we randomly assigned masses from a uniform distribution and obtained their photometry using the isochrones^{9} 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^{10} 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 phasespace cluster model based only on sources within magnitudes bins with contamination rates < pagination10%. This effectively discarded from the inference of the phasespace cluster model those sources fainter than 20 mag in g band. Although these sources were not used in the inference of the cluster phasespace model, they remain in the list of cluster candidate members. Due to our conservative approach to estimating the phasespace model, which filters out sources with G > 20 mag, our membership algorithm may have missed some of the faintest cluster members.
3.2. Phasespace structure
We analysed the phasespace 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 grouplevel parameters are modelled with a hierarchical prior, whereas the stellar or sourcelevel 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 nonnegligible 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 phasespace modelling is conceptually simple but allows constructing complex models based on Gaussian building blocks. It is not the first time that this rather simple multicomponent 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 phasespace 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 BTSettl (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 opensource code Sakam^{11} (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 totaltoselective 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 BTSettl 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
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 selfgravitating, 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’ presentday positions and velocities forward in time.
3.5. Traceforward
To determine the time of the flyby between Coma Ber and Latyshev 2, we traced the presentday phasespace coordinates inferred with Kalkayotl (see Sect. 3.2) forward into the future. We did this not only the grouplevel parameters (i.e. mean positions and velocities of the groups) but also the sourcelevel parameters (i.e. positions and velocities of the candidate members). The uncertainties of these phasespace coordinates were taken into account by tracing 100 samples from the posterior distributions of the group and sourcelevel 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 sampleparticles.
We integrated the orbit of each of these sampleparticle 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^{−1} at a solar radius of 8 kpc (for specific details of this gravitational potential, see Bovy 2015).
Finally, for each grouplevel sampleparticle, 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 grouplevel parameters to estimate the flyby time rather than the sourcelevel 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 sourcelevel 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^{10} 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%.
Confusion matrix and quality indicators of the Coma Ber classifier as a function of photometric magnitude.
Confusion matrix and quality indicators of the Latyshev 2 classifier as a function of photometric magnitude.
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^{−1}, 0.08 mas yr^{−1}, 0.09 mas yr^{−1}, and 0.10 mas yr^{−1} 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 nonsingle or has a problematic astrometric solution^{12}. 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 nonsingle 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. Phasespace 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 phasespace 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 nonnegligible 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^{13} is a confirmed binary star according to SIMBAD, while sources 1560938316782309248^{14} and 1560784144636654976^{15} 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 grouplevel parameters of Coma Ber’s CGMM and Latyhsev 2 + Mecayotl 1 GMM phasespace models, respectively. These figures show multidimensional projections of the phasespace, 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 phasespace 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 grouplevel 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 grouplevel 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 phasespace 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).
Phasespace grouplevel parameters of the identified substructures.
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 onedimensional 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 multidimensional 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^{−1}, which is similar to that of Coma Ber’s core (σ_{UVW} = 0.89 ± 0.10 km s^{−1}) and that of the archetypical Pleiades open cluster (0.8 ± 0.1 km s^{−1} according to Galli et al. 2017; and 0.83 ± 0.07 km s^{−1} according to Torres et al. 2021, as computed from their 1D radial velocity dispersion of 0.48 ± 0.04 km s^{−1}). On the other hand, Mecayotl 1 has a total velocity dispersion of σ_{UVW} = 1.85 ± 0.18 km s^{−1}, 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^{−1} of the βPic stellar association (MiretRoig 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 colourmagnitude diagrams of our candidate members, together with the theoretical isochrones from the PARSEC and BTSettl 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 colourmagnitude diagrams of the candidate members of Coma Ber, Latyshev 2 and Mecayotl 1. The lines depict the theoretical isochrone models with colourcoded 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 lowerenvelope 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 BTSettl model, where there is a concentration of samples from sources more luminous than the BTSettl 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 BTSettl 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 BTSettl and PB models. Given that the BTSettl models have been specifically created for lowmass stars and browndwarfs and the bulk of our candidate members is made of lowmass stars, in the following, we use the mass estimates obtained from the unified PB model, this is BTSettl 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^{2} s^{−2} and 95% having energies smaller than 1.0 M_{⊙} km^{2} s^{−2}. 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^{2} s^{−2} and 75% having energies smaller than 1.0 M_{⊙} km^{2} s^{−2}. 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^{2} s^{−2}. 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 distanceenergy 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 overdensities in the distance to the centre. We spot at least two clear overdensities 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 overdensities have energy dispersion larger than that of Coma Ber’s tails. We discuss the implications of these overdensities and their dispersion in Sect. 5.
4.5. Traceforward
As mentioned in Sect. 3.5, we use the source and grouplevel 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 grouplevel parameters. Then, we estimate the candidate members’ binding energies with respect to the other groups using the sourcelevel 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 sampleparticles drawn from the posterior distribution of the group central phasespace 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 sampleparticle 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 sampleparticle 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 nearzero 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 presentday values of the energy of each sampleparticle with respect to the other groups as an additional validation to our membership classification. We find that none of our candidate members has presentday 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 sampleparticle 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.
Number of members and mean properties of the stellar systems.
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. twobody 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_{1} or L_{2}, 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^{−1}) similar to that of the Pleiades (0.8 ± 0.1 km s^{−1}, Galli et al. 2017) and the core of Coma Ber (0.89 ± 0.1 km s^{−1}). Furthermore, its CEDF and energydistance 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 energybound 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^{−1}), similar to that of the βPic stellar association (1.73 ± 0.3 km s^{−1}MiretRoig et al. 2020), and even larger than that of the tidal tails of Coma Ber (1.43 ± 0.1 km s^{−1}). Furthermore, its CEDF and energydistance 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 overdensities (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 TOI2048 larger and smaller than 2 km s^{−1}). 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 centreright. 
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^{16} 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.
Number of Coma Ber members from the literature that we recover and reject.
Number of Group X (Latyshev 2 + Mecayotl 1) members from the literature that we recover and reject.
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^{−1} with respect to TOI2048. 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^{−1} sample and 90% of their Δv_{tan} > 2 km s^{−1} 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 56% 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^{−1} is similar to that of Fürnkranz et al. (2019), whereas the sample with Δv_{tan} > 2 km s^{−1} has worse quality indicators than those of a random classifier.
Recovery, contamination, and missing rates of the literature classifiers.
From our previous estimates, we notice the following points. First, the StarGo classifier of Tang et al. (2019) misses only onethird 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 twothirds 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 machinelearning algorithms they use and their finetuned and manually chosen parameters, particularly in the number of neurons used by the selforganising 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 phasespace 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 machinelearning 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 lowmass members (see discussion in Sect. 5.4). Future improvements of the tail’s phasespace model are warranted.
5.3. Phasespace structure
Our comprehensive phasespace 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 grouplevel 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 overdensities described in Sect. 4.4, particularly seen in the Y distribution. In Coma Ber, the core has two overdensities 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 overdensities are shifted out at −9 pc and 11 pc. We interpret these overdensities 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 overdensities 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 timedependent 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 futurebound 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 overdensities located at the wings of the V distribution. Given that these overdensities 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 flattop 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 zerovelocity peak. Although significant only at the onesigma level, the observed expansion of the tails (0.7 ± 0.7 km s^{−1}; 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 nonzero rotation velocities (see Table 3), with similar values to those observed in open clusters and starforming 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^{−1} reported by Loktin & Popov (2020), and the rotation velocity of Mecayotl 1, 1.32 ± 0.74 km s^{−1}, is also consistent with the 1.5 ± 0.1 km s^{−1} reported by Galli et al. (2019) for the Taurus complex. In addition, the total rotation of Coma Ber (1.13 ± 0.75 km s^{−1}, 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 nonzero and compatible with those from other stellar systems in the solar neighbourhood, their significance still not reaches the twosigma 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 leadingtotrailing tail number ratio of 1.55 ± 0.34^{17}. 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 overdensities. We investigate the statistical significance of these overdensities 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 onesigma (blue dashed lines), twosigma (orange dashed lines) and threesigma (green dashed lines) confidence intervals. As can be observed, only two radial bins in each tail have densities (median minus standard deviation) exceeding the twosigma (orange arrows) and threesigma (green arrows) levels of the distribution. In the leading tail these correspond to the overdensities located at 35 ± 1 pc and 45 ± 1 pc, while in the trailing tail, these correspond to the overdensities at 29 ± 1 pc and 33 ± 1 pc. We notice that these overdensities are still significantly detected for binwidths of up to ∼4 pc and up to ∼3 pc in the leading and trailing tails, respectively.
Fig. 15. Distancetocentre 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 onesigma (blue dashed lines), twosigma (orange dashed lines) and threesigma (green dashed lines) confidence intervals are also plotted. The vertical arrows indicate the radial bin with densities exceeding the twosigma (orange arrows) and threesigma (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 overdensities 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 overdensities 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 overdensities 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 overdensities 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 presentday 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 overdensities were produced at the perigalacticon and thus the cluster has an eccentric orbit. Future work will be needed to clarify if the overdensities 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 lowmass domain. We interpret these differences as resulting from the ejection and evaporation of lowmass 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 BTSettl and PB) and ∼0.25 M_{⊙} (0.28 M_{⊙} for PARSEC and 0.22 M_{⊙} for BTSettl 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 lowmass 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 lowmass stars in the most distant edges of the tidal tails. Identifying this distant population of lowmass 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 phasespace 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 phasespace coordinates and masses sampled randomly from the King (1966) and Kroupa (2001) distributions, respectively. These clusters have halfmass radii of 1 pc and concentration parameters, W_{0}, varying between 1 and 12. We choose the total mass and halfmass radii of these synthetic clusters to resemble those of Come Ber; the chosen W_{0} 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 colourcoded with the value of their W_{0} 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 phasespace 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 timescale is of the order of the orbital timescale of the cluster (Küpper et al. 2012).
5.7. Traceforward
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 Nbody 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 traceback 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 Nbody traceback 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 traceback simulations.
6. Conclusions
Using public data from 2MASS, APOGEE, Gaia, and PanSTARRS, in combination with our own INTIDS observations, we studied the wellknown 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 lowmass 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^{−1}), 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 ultracool dwarfs and exoplanets.
Mecayotl 1 is a relatively populous (146 candidate members), young (400600 Myr), and sparse stellar association. It has a large velocity dispersion (1.85 ± 0.18 km s^{−1}), a total mass of 75 ± 5 M_{⊙}, and an energy distribution with more dispersion than that of the tails of Coma Ber.
Our traceforward 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 traceback 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 phasespace models for the tidal tails, traceback and traceforward numerical Nbody 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 ultracool dwarfs in the Euclid (Laureijs et al. 2011) era.
The confusion matrix of a binary classifier is a four entries matrix with the number of true positives (TP) and true negatives (TN) in the diagonal, and the number of false positives (FP) and false negatives (FN) in the antidiagonal. The classifier quality can then be evaluated using this matrix and a suitable metric.
The Latyshev (1977) sources that we recovered as candidate members are 81 UMa, HD 119765, and 84 UMa, with the following Gaia source_id: 1562168842092340352, 1559292485314843008 and 1561439694084279552, respectively. We notice that 83 UMa and 86 UMa could be ejected members that, due to their discrepant distance and large negative radial velocities, were not identified as members by our membership algorithm.
Given the large TPR and low FPR of our classifier (see Sect. 4.1) we expect the uncertainties in the number of sources in the leading and trailing tails to be Poisson dominated.
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 (AEIMCINN) under grant PID2019109522GBC53. 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/125188 and 2021/117789. Cofunded 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 highperformance 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 DeimosDiva 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 PanSTARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the PanSTARRS Project Office, the MaxPlanck 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 HarvardSmithsonian 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. AST1238877, 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. SDSSIV acknowledges support and resources from the Center for High Performance Computing at the University of Utah. The SDSS website is www.sdss.org. SDSSIV 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), MaxPlanckInstitut für Astronomie (MPIA Heidelberg), MaxPlanckInstitut für Astrophysik (MPA Garching), MaxPlanckInstitut 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 GodoyRivera for the useful discussions that took place during the Cool Stars 21 and greatly improved the quality of this work.
References
 Abdurro’uf, Accetta, K., Aerts, C., et al. 2022, ApJS, 259, 35 [CrossRef] [Google Scholar]
 Akaike, H. 1974, IEEE Trans. Autom. Control, 19, 716 [Google Scholar]
 Alam, S., Albareti, F. D., Allende Prieto, C., et al. 2015, ApJS, 219, 12 [Google Scholar]
 Allard, F. 2014, in Exploring the Formation and Evolution of Planetary Systems, eds. M. Booth, B. C. Matthews, & J. R. Graham, 299, 271 [Google Scholar]
 Andrae, R., Fouesneau, M., Sordo, R., et al. 2023, A&A, 674, A27 [CrossRef] [EDP Sciences] [Google Scholar]
 Angelo, M. S., Santos, J. F. C., Maia, F. F. S., & Corradi, W. J. B. 2022, MNRAS, 510, 5695 [CrossRef] [Google Scholar]
 Angus, R., Aigrain, S., ForemanMackey, D., & McQuillan, A. 2015, MNRAS, 450, 1787 [Google Scholar]
 Archinal, B. A., & Hynes, S. J. 2003, Star Clusters (Richmond, VA: WillmannBell) [Google Scholar]
 Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Astropy Collaboration (PriceWhelan, A. M., et al.) 2018, AJ, 156, 123 [Google Scholar]
 Baraffe, I., Homeier, D., Allard, F., & Chabrier, G. 2015, A&A, 577, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 BlancoCuaresma, S. 2019, MNRAS, 486, 2075 [Google Scholar]
 BlancoCuaresma, S., Soubiran, C., Heiter, U., & Jofré, P. 2014, A&A, 569, A111 [CrossRef] [EDP Sciences] [Google Scholar]
 Bovy, J. 2015, ApJS, 216, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Bovy, J., Hogg, D. W., & Roweis, S. T. 2011, Ann. Appl. Stat., 5, 1657 [Google Scholar]
 Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245 [Google Scholar]
 Casal, E., Fernández, M., Alfaro, E. J., Casanova, V., & Tobaruela, Á. 2020, MNRAS, 497, 2562 [NASA ADS] [CrossRef] [Google Scholar]
 Casamiquela, L., Olivares, J., Tarricq, Y., et al. 2022, A&A, 664, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chabrier, G. 2005, in The Initial Mass Function: From Salpeter 1955 to 2005, eds. E. Corbelli, F. Palla, & H. Zinnecker, 327, 41 [Google Scholar]
 Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, ArXiv eprints [arXiv:1612.05560] [Google Scholar]
 Chen, Y., Girardi, L., Bressan, A., et al. 2014, MNRAS, 444, 2525 [Google Scholar]
 Claydon, I., Gieles, M., Varri, A. L., Heggie, D. C., & Zocchi, A. 2019, MNRAS, 487, 147 [NASA ADS] [CrossRef] [Google Scholar]
 Collier Cameron, A., Davidson, V. A., Hebb, L., et al. 2009, MNRAS, 400, 451 [NASA ADS] [CrossRef] [Google Scholar]
 Douglas, S. T., Curtis, J. L., Agüeros, M. A., et al. 2019, ApJ, 879, 100 [Google Scholar]
 Dravins, D., Lindegren, L., & Madsen, S. 1999, A&A, 348, 1040 [NASA ADS] [Google Scholar]
 Ester, M., Kriegel, H. P., Sander, J., & Xu, X. 1996, in Proceedings of the Second International Conference on Knowledge Discovery and Data Mining, KDD’96 (AAAI Press), 226 [Google Scholar]
 Faherty, J. K., Bochanski, J. J., Gagné, J., et al. 2018, ApJ, 863, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Feiden, G. A., & Chaboyer, B. 2012, ApJ, 757, 42 [Google Scholar]
 Fürnkranz, V., Meingast, S., & Alves, J. 2019, A&A, 624, L11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaia Collaboration (Prusti, T., et al.) 2016, A&A, 595, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaia Collaboration (van Leeuwen, F., et al.) 2017, A&A, 601, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaia Collaboration (Brown, A. G. A., et al.) 2021, A&A, 649, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaia Collaboration (Arenou, F., et al.) 2023a, A&A, 674, A34 [CrossRef] [EDP Sciences] [Google Scholar]
 Gaia Collaboration (Vallenari, A., et al.) 2023b, A&A, 674, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Galli, P. A. B., Moraux, E., Bouy, H., et al. 2017, A&A, 598, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Galli, P. A. B., Loinard, L., Bouy, H., et al. 2019, A&A, 630, A137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Green, G. 2018, J. Open Source Softw., 3, 695 [NASA ADS] [CrossRef] [Google Scholar]
 Guo, D., de Koter, A., Kaper, L., Brown, A. G. A., & de Bruijne, J. H. J. 2021, A&A, 655, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 He, Z., Wang, K., Luo, Y., et al. 2022, ApJS, 262, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Jadhav, V. V., Roy, K., Joshi, N., & Subramaniam, A. 2021, AJ, 162, 264 [NASA ADS] [CrossRef] [Google Scholar]
 Jerabkova, T., Boffin, H. M. J., Beccari, G., et al. 2021, A&A, 647, A137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Katz, D., Sartoretti, P., Guerrier, A., et al. 2023, A&A, 674, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 King, I. 1962, AJ, 67, 471 [Google Scholar]
 King, I. R. 1966, AJ, 71, 64 [Google Scholar]
 King, I. R. 2008, in Dynamical Evolution of Dense Stellar Systems, eds. E. Vesperini, M. Giersz, & A. Sills, 246, 131 [NASA ADS] [Google Scholar]
 Kraus, A. L., Shkolnik, E. L., Allers, K. N., & Liu, M. C. 2014, AJ, 147, 146 [Google Scholar]
 Kroupa, P. 2001, MNRAS, 322, 231 [NASA ADS] [CrossRef] [Google Scholar]
 Kroupa, P., Tout, C. A., & Gilmore, G. 1990, MNRAS, 244, 76 [NASA ADS] [Google Scholar]
 Küpper, A. H. W., MacLeod, A., & Heggie, D. C. 2008, MNRAS, 387, 1248 [Google Scholar]
 Küpper, A. H. W., Kroupa, P., Baumgardt, H., & Heggie, D. C. 2010, MNRAS, 407, 2241 [NASA ADS] [CrossRef] [Google Scholar]
 Küpper, A. H. W., Lane, R. R., & Heggie, D. C. 2012, MNRAS, 420, 2700 [Google Scholar]
 Küpper, A. H. W., Maschberger, T., Kroupa, P., & Baumgardt, H. 2011, MNRAS, 417, 2300 [Google Scholar]
 Latyshev, I. N. 1977, Astronomicheskij Tsirkulyar, 969, 7 [NASA ADS] [Google Scholar]
 Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, ArXiv eprints [arXiv:1110.3193] [Google Scholar]
 Lawrence, A., Warren, S. J., Almaini, O., et al. 2007, MNRAS, 379, 1599 [Google Scholar]
 Lindegren, L., Madsen, S., & Dravins, D. 2000, A&A, 356, 1119 [NASA ADS] [Google Scholar]
 Lindegren, L., Lammers, U., Hobbs, D., et al. 2012, A&A, 538, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lodieu, N., Smart, R. L., PérezGarrido, A., & Silvotti, R. 2019, A&A, 623, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lindegren, L., Klioner, S. A., Hernández, J., et al. 2021, A&A, 649, A2 [EDP Sciences] [Google Scholar]
 Loktin, A. V., & Popov, A. A. 2020, Astron. Nachr., 341, 638 [NASA ADS] [CrossRef] [Google Scholar]
 Madsen, S., Dravins, D., & Lindegren, L. 2002, A&A, 381, 446 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Majewski, S. R., Schiavon, R. P., Frinchaboy, P. M., et al. 2017, AJ, 154, 94 [Google Scholar]
 Mamajek, E. E. 2016, in Young Stars& Planets Near the Sun, eds. J. H. Kastner, B. Stelzer, & S. A. Metchev, 314, 21 [NASA ADS] [Google Scholar]
 Mamajek, E. E., & Hillenbrand, L. A. 2008, ApJ, 687, 1264 [Google Scholar]
 Marigo, P., Bressan, A., Nanni, A., Girardi, L., & Pumo, M. L. 2013, MNRAS, 434, 488 [Google Scholar]
 Martín, E. L., Lodieu, N., & Béjar, V. J. S. 2020, A&A, 640, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Matthews, B. 1975, Biochimica et Biophysica Acta (BBA)  Protein Structure, 405, 442 [CrossRef] [Google Scholar]
 Meingast, S., Alves, J., & Fürnkranz, V. 2019, A&A, 622, L13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Messina, S., Nardiello, D., Desidera, S., et al. 2022, A&A, 657, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 MiretRoig, N., Galli, P. A. B., Brandner, W., et al. 2020, A&A, 642, A179 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Moranta, L., Gagné, J., Couture, D., & Faherty, J. K. 2022, ApJ, 939, 94 [NASA ADS] [CrossRef] [Google Scholar]
 Morton, T. D. 2015, Isochrones: Stellar Model Grid Package [Google Scholar]
 Newton, E. R., Rampalli, R., Kraus, A. L., et al. 2022, AJ, 164, 115 [Google Scholar]
 Oh, S., PriceWhelan, A. M., Hogg, D. W., Morton, T. D., & Spergel, D. N. 2017, AJ, 153, 257 [NASA ADS] [CrossRef] [Google Scholar]
 Olivares, J., Moraux, E., Sarro, L. M., et al. 2018a, A&A, 612, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Olivares, J., Sarro, L. M., Moraux, E., et al. 2018b, A&A, 617, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Olivares, J., Bouy, H., Sarro, L. M., et al. 2019, A&A, 625, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Olivares, J., Sarro, L. M., Bouy, H., et al. 2020, A&A, 644, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pastorelli, G., Marigo, P., Girardi, L., et al. 2020, MNRAS, 498, 3283 [Google Scholar]
 Piatti, A. E., & Malhan, K. 2022, MNRAS, 511, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Planck Collaboration Int. XLVIII. 2016, A&A, 596, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, J. Astron. Telescopes Instrum. Syst., 1, 014003 [Google Scholar]
 Sapozhnikov, S., & Kovaleva, D. 2021, Open Astron., 30, 191 [NASA ADS] [CrossRef] [Google Scholar]
 Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163 [Google Scholar]
 Sollima, A., CarballoBello, J. A., Beccari, G., et al. 2010, MNRAS, 401, 577 [NASA ADS] [CrossRef] [Google Scholar]
 Soubiran, C., Jasniewicz, G., Chemin, L., et al. 2013, A&A, 552, A64 [CrossRef] [EDP Sciences] [Google Scholar]
 Souto, D., Cunha, K., & Smith, V. V. 2021, ApJ, 917, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Tang, S.Y., Chen, W. P., Chiang, P. S., et al. 2018, ApJ, 862, 106 [CrossRef] [Google Scholar]
 Tang, S.Y., Pang, X., Yuan, Z., et al. 2019, ApJ, 877, 12 [Google Scholar]
 Tarricq, Y., Soubiran, C., Casamiquela, L., et al. 2022, A&A, 659, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Torres, G., Latham, D. W., & Quinn, S. N. 2021, ApJ, 921, 117 [Google Scholar]
 Trumpler, R. J. 1938, Lick Obs. Bull., 494, 167 [NASA ADS] [Google Scholar]
 Webb, J. J., Hunt, J. A. S., & Bovy, J. 2023, MNRAS, 521, 3898 [NASA ADS] [CrossRef] [Google Scholar]
 Wenger, M., Ochsenbein, F., Egret, D., et al. 2000, A&AS, 143, 9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wielen, R. 1974, Highlights Astron., 3, 395 [NASA ADS] [CrossRef] [Google Scholar]
 Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868 [Google Scholar]
 Ye, X., Zhao, J., Oswalt, T. D., Yang, Y., & Zhao, G. 2022, AJ, 164, 132 [NASA ADS] [CrossRef] [Google Scholar]
 Zacharias, N., Finch, C., Subasavage, J., et al. 2015, AJ, 150, 101 [Google Scholar]
 Zhong, J., Chen, L., Jiang, Y., Qin, S., & Hou, J. 2022, AJ, 164, 54 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: INTIDS observations
We use INTIDS 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^{−1}) or nonexisting 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^{−1}, with a typical value of ∼10 km s^{−1}. 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 zeropoints by median combining the standards of each night assuming Soubiran et al. (2013) values. These zeropoints were used to correct our target’s radial velocity, and its uncertainty was added in quadrature.
Observing log of our 36 radial velocity targets.
Observing log of the four radial velocity standards.
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 onetoone 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 zeropoints
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 zeropoint 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 zeropoint offset of SIMBAD, APOGEE and IDS are 0.02 ± 1.46 km s^{−1}, 0.44 ± 0.81 km s^{−1}, and 0.32 ± 8.86 km s^{−1}, 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: dwarfturnoff 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^{−1}) 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, BTSettl, and PB models (colourcoded 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 A_{v} as inferred using a uniform prior (A_{v} ∈ [0, 10] mag, grey lines) and different theoretical models (colour coded legend). The black and blue dashed lines show the A_{v} 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 BTSettl and the joint PB models deliver bimodal distributions with two overdensities, 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 BTSettl models, which is visible in Fig. 4 at the colour interval G − RP ∈ [0.8, 1.2] mag where the BTSettl 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 nearinfrared photometric sequences of lowmass stars in metalpoor Globular clusters and the solarmetallicity 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 BTSettl 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 lowmass 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 lowmass 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
Properties of the 634 candidate members.
All Tables
Confusion matrix and quality indicators of the Coma Ber classifier as a function of photometric magnitude.
Confusion matrix and quality indicators of the Latyshev 2 classifier as a function of photometric magnitude.
Number of Group X (Latyshev 2 + Mecayotl 1) members from the literature that we recover and reject.
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 grouplevel 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 colourmagnitude diagrams of the candidate members of Coma Ber, Latyshev 2 and Mecayotl 1. The lines depict the theoretical isochrone models with colourcoded 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 BTSettl 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 centreright. 

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 grouplevel parameters, respectively. 

In the text 
Fig. 15. Distancetocentre 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 onesigma (blue dashed lines), twosigma (orange dashed lines) and threesigma (green dashed lines) confidence intervals are also plotted. The vertical arrows indicate the radial bin with densities exceeding the twosigma (orange arrows) and threesigma (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 A_{v} as inferred using a uniform prior (A_{v} ∈ [0, 10] mag, grey lines) and different theoretical models (colour coded legend). The black and blue dashed lines show the A_{v} distributions of our candidate members as reported by Gaia DR3 and Planck Collaboration Int. XLVIII (2016), respectively. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.