From ridges to manifolds: 3D characterization of the moving groups in the Milky Way disc

A&A 667, A116 (2022) https://doi.org/10.1051/0004-6361/202244070 c M. Bernet et al. 2022 Astronomy & Astrophysics From ridges to manifolds: 3D characterization of the moving groups in the Milky Way disc M. Bernet1,2,3 , P. Ramos1,2,3,4 , T. Antoja1,2,3 , B. Famaey4, G. Monari4, H. Al Kazwini4, and M. Romero-Gómez1,2,3 1Departament de Física Quàntica i Astrofísica (FQA), Universitat de Barcelona (UB), C Martí i Franquès, 1, 08028 Barcelona, Spain e-mail: [email protected] 2Institut de Ciències del Cosmos (ICCUB), Universitat de Barcelona (UB), C Martí i Franquès, 1, 08028 Barcelona, Spain 3Institut d’Estudis Espacials de Catalunya (IEEC), C Gran Capità, 2-4, 08034 Barcelona, Spain 4Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, 11 rue de l’Université, 67000 Strasbourg, France Received 20 May 2022 /Accepted 25 July 2022 ABSTRACT Context. The details of the effect of the bar and spiral arms on the disc dynamics of the Milky Way are still unknown. The stellar velocity distribution in the solar neighbourhood displays kinematic substructures, which are possibly signatures of these processes and of previous accretion events. With the Gaia mission, more details of these signatures, such as ridges in the Vφ−Rplane and thin arches in the Vφ−VRplane, have been revealed. The positions of these kinematic substructures, or moving groups, can be thought of as continuous manifolds in the 6D phase space, and the ridges and arches as specific projections of these manifolds. Aims. Our aim is to detect and characterize the moving groups along the Milky Way disc, sampling the galactocentric radial and azimuthal velocities of the manifolds through the three dimensions of the disc: radial, azimuthal, and vertical. Method. We developed and applied a novel methodology to perform a blind search for substructure in the Gaia EDR3 6D data, which consists in the execution of the wavelet transform in independent small volumes of the Milky Way disc, and the grouping of these local solutions into global structures with a method based on the breadth-first search algorithm from graph theory. We applied the same methodology to simulations of barred galaxies to validate the method and for comparison with the data. Results. We reveal the skeleton of the velocity distribution, uncovering projections that were not possible before. We sample nine main moving groups along a large region of the disc in configuration space, covering up to 6kpc, 60deg, and 2kpc in the radial, azimuthal, and vertical directions, respectively, extending significantly the range of previous analyses. In the radial direction we find that the groups deviate from the lines of constant angular momentum that one would naively expect from an epicyclic approximation analysis of the first-order effects of resonances. We reveal that the spatial evolution of the moving groups is complex and that the configuration of moving groups in the solar neighbourhood is not maintained along the disc. We also find that the azimuthal velocity of the moving groups that are mostly detected in the inner parts of the disc (Arcturus, Bobylev, and Hercules) is non-axisymmetric. For Hercules we measure an azimuthal gradient of −0.50 km s−1deg−1at R=8kpc. We detect a vertical asymmetry in the azimuthal velocity for the Coma Berenices moving group, which is not expected for structures originating from a resonance of the bar, supporting the previous hypothesis of the incomplete vertical phase mixing of the group. In our simulations we extract substructures corresponding to the outer Lindblad resonance and the 1:1 resonances and observe the same deviation from constant angular momentum lines and the non-axisymmetry of the azimuthal velocities of the moving groups in the inner part of the disc. Conclusions. This data-driven characterization is a starting point for a holistic understanding of the moving groups. It also allows for a quantitative comparison with models, providing a key tool to comprehend the dynamics of the Milky Way. Key words. Galaxy: disk – Galaxy: kinematics and dynamics – Galaxy: structure – Galaxy: evolution – methods: data analysis 1. Introduction The stellar velocity distribution in the solar neighbourhood (SN) has long been a key element in our understanding of the structure of the Milky Way (MW) (Dehnen & Binney 1998; Skuljan et al. 1999;Famaey et al. 2005;Antoja et al. 2008). Historically, several overdensities in this velocity distribution have been identified and discussed (Pleiades, Hyades, Sirius). These moving groups, as they are usually called, can be related to the orbital resonances of the bar and spiral arms of the Galaxy (Kalnajs 1991;Dehnen 2000;Antoja et al. 2011; Fragkoudi et al. 2019;Monari et al. 2019a) and/or attributed to ongoing phase mixing related to external perturbations (Minchev et al. 2009;Gómez et al. 2012;Antoja et al. 2018; Ramos et al. 2018;Hunt et al. 2018a;Khanna et al. 2019; Laporte et al. 2019,2020). The latest releases of the Gaia mission (Gaia Collaboration 2018a,2021a) have provided a full 6D phase-space catalogue of 7.2 million stars, increasing the size and precision of any previous survey by several orders of magnitude. This has been a game changer in many fields of astrophysics. In the SN the new high resolution velocity distribution has revealed a complex substructure with several thin arches never observed before (Gaia Collaboration 2018b). When extending the study to the entire disc, large ridges appeared in the R−Vφ(respectively, galactocentric radius and azimuthal velocity) diagram covering several kiloparsecs (Antoja et al. 2018;Kawata et al. 2018; Fragkoudi et al. 2019). 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. A116, page 1 of 16 A&A 667, A116 (2022) Orbits in a barred potential can be trapped into resonances (Weinberg 1994). Dehnen (2000) showed that in a short–fast bar scenario (i.e. Ωb=50 km s−1kpc−1) the transition between two types of non-axisymmetric orbital families across the bar’s outer Lindblad resonance (OLR) can explain the bi-modality formed by Hercules and the rest of the velocity distribution in the SN if the Sun is placed just outside the OLR of the bar (ROLR ≈ 7.2 kpc). This scenario was consistent with the gas dynamics measurements of the inner MW at the time. Later on, studies of star counts and the kinematics of the inner MW suggested that the bar might be longer and slower than previously thought (Portail et al. 2017). In this case, the OLR would be placed further out (ROLR ≈10.5 kpc, perhaps matching other groups such as the Arch/Hat instead of Hercules) and co-rotation (CR) would be closer to the SN (RCR ≈6 kpc). Pérez-Villegas et al. (2017) and Monari et al. (2019a) then explained Hercules as the overdensity formed by the orbits trapped at the CR, librating around the Lagrangian points of a long–slow bar. This moving group created by CR seems to be less pronounced than the one produced by the OLR (Binney 2018;Hunt et al. 2018b). However, Hunt et al. (2018a) showed that the addition of spiral structure in combination with the CR might create a strong distinct Hercules, and Chiba et al. (2021) showed that a decelerating Galactic bar could enhance the occupation on resonances, being able to reproduce Hercules through the CR resonance. This shows that the value of the pattern speed (Ωb) of the MW bar and the exact link between substructures and resonances are still a matter of debate, and more observables are needed to obtain a final answer. In this direction, Ramos et al. (2018, hereafter R18), used the wavelet transform (WT, Starck & Murtagh 2002;Chereul et al. 1999) to detect and characterize the kinematics of the moving groups along the disc. They matched the spatial evolution of the groups with the ridges in the R−Vφplane. They claimed that some of the arches follow lines of constant energy at a given volume, which could be related to phase mixing processes (Minchev et al. 2009;Gómez et al. 2012), and that others follow lines of common angular momentum in the radial direction, as expected approximately in the case of resonant kinematic substructures (e.g. Quillen et al. 2018a). They also claimed that the observed changes in the azimuthal direction for the Hercules moving group are consistent with being produced by the OLR of a short–fast bar (Dehnen 2000;Fux 2001;Antoja et al. 2014). The long–slow paradigm is relatively recent and there have been few analyses on the azimuthal variations of a substructure caused by CR. Monari et al. (2019b) found that the Hercules angular momentum changes significantly with azimuth as they predicted analytically for the co-rotation resonance of an old long–slow bar. They showed that the only way to obtain a similar change in azimuth for an OLR origin of Hercules would be if the orbits are still far from phase-mixed in the bar potential (bar perturbation younger than 2 Gyr; see also Trick et al. 2021). The link between the moving groups across the neighbourhoods in R18 was made visually, using a scatter plot of two variables and colour as a third variable. Therefore, the analysis of the moving groups link was restricted to three variables. Since VRand Vφare compulsory to select the moving groups, this limitation restricted the analysis to one dimension in space (either radial or azimuthal). The vertical direction was not explored. The correlation between the position of the moving groups (overdensities in VR−Vφ) and the ridges (overdensities in R−Vφ) indicate that both are projections of the same substructure in the 6D phase-space onto different planes. The positions of these kinematic substructures, also known as moving groups, can be described as continuous manifolds in 6D phase space, and the ridges and arches as specific projections of these manifolds. Our goal in this article is to extend the idea introduced in R18 by automatizing the n-dimensional link of the moving groups to avoid the limitation of projecting the data. Using this process, we intend to move from a ridge-moving group paradigm to a manifold paradigm, where we sample the position of these manifolds in the (R, φ, Z,VR,Vφ) space for each moving group. We present a novel methodology of detecting these manifolds in a dataset. It is based on the execution of the WT in independent small volumes, and the relation of these local solutions in global substructures with an algorithm based on the breadth-first search (BFS) algorithm from graph theory. With this methodology we processed the Gaia EDR3 6D data and detect the positions of the groups across the MW disc. We also sampled the manifolds of two test particle simulations with a fast (Ωb=50 km s−1kpc−1) and a slow (Ωb=30 km s−1kpc−1) bar, both to test our methodology and to compare it to the data. The Gaia DR3 (Gaia Collaboration 2022a) catalogue includes a larger and updated sample of radial velocities (33 M of stars, Gaia Collaboration 2022b), which covers a larger region of the MW disc and increases the resolution (number of stars and precision) of Gaia EDR3. This provides finer observables to untangle the different contributions in the complex dynamics of the Galaxy. To exploit these data in their totality new strategies must be developed (e.g. Contardo et al. 2022), such as the one we present here, to avoid the current limitations of the analysis. This paper is organized as follows. In Sect. 2we describe the observational data we used. In Sect. 3we introduce the methodology we developed. In Sect. 4we show the results of the application of the method to Gaia EDR3 data. In Sect. 5we present and analyse the simulations. In Sect. 6we compare the results from the data and the simulations, and with previous results in the literature. Finally, in Sect. 7we list the main conclusions of this work. 2. Data and sample preprocessing The Gaia Early Data Release 3 (EDR3) consists of an updated and enlarged source list, with improved astrometry and photometry. About 7.2 million stars have proper motion and a radial velocity measurements in the Gaia DR2, most of which are transferred to EDR3 (Seabroke et al. 2021;Torra et al. 2021). For this section we use a subset of these stars with photogeometric distances from Bailer-Jones et al. (2021), derived from a probabilistic approach including colour and apparent magnitude information. Distance is a critical parameter in the computation of the motion and position of the star in the 6D study of the MW, and a major source of uncertainty. This is why, in order to improve the quality of the sample, we also apply a cut in relative parallax error: $ σ$ >5.(1) The resulting sample contains 6 059 648 sources. In Appendix A, we study the errors that an overestimation or an underestimation of the distance could produce in our results. We determine that, in general, this error would be below 2 km s−1for a distance bias of ±10%, and that it would not affect the overall trend of the groups. We use a cylindrical galactocentric coordinate system, fixing the reference at the Galactic Centre with the radial direction A116, page 2 of 16 M. Bernet et al.: From ridges to manifolds: 3D characterization of the moving groups in the Milky Way disc (R) pointing outwards from it, the azimuthal (φ) negative in the direction of rotation, and the vertical component (Z) positive towards the north galactic pole. To transform Gaia observables to positions and velocities in this reference frame, we take the Sun to be at R=8.178 kpc (GRAVITY Collaboration 2019), φ=0◦and Z=0.0208 kpc (Bennett & Bovy 2019). For the solar motion, we use U=11.1,vcirc +V=248.5,W= 7.25 km s−1(Schönrich et al. 2010;Reid & Brunthaler 2020). 3. Method The data described in the previous section contains the 6D variables of position (R, φ, Z) and velocity (VR,Vφ,VZ) of the stars. Inside small volumes, cuts in (R, φ, Z), the moving groups appear as well-defined overdensities in the velocity distribution VR−Vφ (R18), which are easy to detect. However, we know that at large spatial scales the position of the overdensities in the velocity space changes (ridges in R−Vφ). Therefore, if we use larger volumes to construct the velocity distribution, the overdensities will blur and become undetectable. In this section we present the novel method that we developed to extract these large kinematic substructures from a dataset. It is divided into two steps: the execution of the WT in independent small volumes of the MW disc, and the relation of these local solutions in global substructures with an algorithm based on the BFS algorithm from graph theory (Moore 1959, described in Sect. 3.2). 3.1. Local wavelet transform We partition the data in a dense grid of small volumes (hereafter pixels) in the spatial coordinates. We construct it as follows: – Radial direction (Ri): [5,14] kpc in steps of 0.04 kpc, Rbin =±0.24 kpc around each centre; – Azimuthal direction (φj): [−34,34] deg in steps of 0.8 deg, φbin =±2.4 deg around each centre; – Vertical direction (Zk): [−1,1] kpc in steps of 0.08 kpc, Zbin =±0.24 kpc around each centre. This produces a dense grid of 2 700 000 pixels, with a maximum volume overlap between consecutive pixels of 83.3%. This bin size is bigger than that used in R18. When reaching regions far from the Sun, the statistical significance of the moving groups decreases and a bigger bin is the only way to obtain a robust determination of the group velocity. The drawback for this enlarged bin is that any inhomogeneity of the sample within a volume (e.g. due to the local extinction, the selection function, a change in the sampled populations) can bias the mean position of the stars, and therefore its kinematics. We tested this effect using a smaller bin, and we estimate the impact to be below 2 km s−1for the bin size used. For each pixel we construct the velocity distribution (VR,Vφ) diagram of the stars in it as a 2D histogram with bins of 1km s−1 (see background histogram in Fig. 1). R18 showed that the overdensities form thin arches elongated around large ranges of VR, with a small variation in Vφ. The use of 2D peak detection algorithms (like the one in R18) is sub-optimal for arch-like structure detection. When analysing regions with few observations, the search for peaks is translated into a very noisy determination of VR, and uncontrollable correlations between VRand Vφ (movement along the arch). To avoid this, we slice each VR−Vφ diagram in vertical columns (bins in VR), and run a 1D WT in the Vφhistogram of each column: radial velocity (VR): −100 – 100 km s−1in steps of 10 km s−1,VRbin =±15 km s−1around each centre. Since we are detecting each part of the arch separately, we avoid the movement along the arch of the overdensities, breaking the degeneracy between VRand Vφin the detection. To detect the peaks, we use the algorithm developed in Du et al. (2006) implemented in scipy (Virtanen et al. 2020) as find_peaks_cwt. This method performs the 1D WT in a range of length scales. A peak is then selected if it is present in enough scales consecutively. In our execution, we use a range of scales of [5,10] km s−1, with steps of 1 km s−1. We keep the peak if it is present in more than two scales consecutively. With this configuration of scales, we lose the thin resolution that we could extract in regions with a large number of sources, but we gain robustness in the detection of the large structures in poorly sampled regions. Since the scope of this work is the large-scale behaviour of the groups, we consider this approach to be better. At the end of the execution, the peak pinherits the spatial position from the pixel, VRfrom the position of the radial velocity bin, and Vφfrom the result of the WT detection. Therefore, a peak has the coordinates p=(R, φ, Z,VR,Vφ)p.(2) 3.2. Breadth-first search resolution with online interpolator We defined the pixels to have large overlaps; for example, two adjacent pixels will share 83.3% of their volume. Therefore, a given substructure in consecutive pixels will have an almost identical shape. We consider a pair of peaks from consecutive pixels to be adjacent if they are in the same VRbin and if their distance in Vφ is smaller than 4 km s−1. In R18, the maximum slope found in a moving group is 33 km s−1kpc−1. In our grid the step is 0.04 kpc, which translates into a maximum change of ≈1 km s−1between adjacent pixels. This 4kms−1limit in the adjacency is a compromise between including very steep groups (about 4 times that detected in R18) and reducing the number of adjacencies, which will determine the computational cost of the next step. It sometimes occurs, especially in poorly sampled regions, that a peak from one pixel can be exactly in the middle of two peaks in the other pixel. In these cases we do not consider any of the peaks adjacent. With this consideration, two adjacent peaks are always strong candidates to belong to the same substructure. This adjacency information constructs an enormous net of linked peaks. Ideally, substructures will be isolated subsets of peaks in this net. These are groups of peaks with no adjacencies to any peak outside their group. In graph theory these nets of linked points are called graphs and the isolated groups are the connected components of a graph. A very common algorithm to extract these connected components is the BFS algorithm, which proceeds as follows: 1. Add (enqueue) the initial peak pto the queue1Q; 2. Select (dequeue) the top peak ptop of the queue Q; 3. Visit all the peaks padj adjacent to ptop. For each adjacent peak, if we have already visited it, ignore the peak. If it is the first time we see the peak, enqueue it; 4. If there are still peaks in the queue, return to step 2; 5. If the queue is empty, our connected component is the list of visited peaks. 1A queue is a data structure similar to an array with limited access to the positions. One end is always used to insert data (enqueue) and the other is used to remove data (dequeue). Queue follows a first-in-first-out methodology, such that the data item stored first will be accessed first. A116, page 3 of 16 A&A 667, A116 (2022) Given an initial peak p0, Algorithm 1(see below) returns the entire substructure where it belongs. By repeating the process for all the non-matched peaks we can extract all the substructures. Algorithm 1 – Breadth-first search 1: queue Q 2: list V 3: add p0to Vand enqueue in Q 4: while Qnot empty do 5: pit ←dequeue Q (remove and assign) 6: for all pad j adjacent to pit do 7: if pad j not in Vthen 8: add pad j to V, enqueue in Q 9: end if 10: end for 11: end while 12: return V This solution would be enough in an ideal case, but in practice undersampling and Poisson noise especially in regions far from the Sun produce confusion and jumps between structures that a straightforward BFS implementation cannot filter out. In order to avoid these jumps between structures, we include an extra step in the algorithm. While the BFS is running, the peaks already matched give us information about the structure. Therefore, in order to accept a new peak, we require it to be consistent with the current structure. Let us suppose we have a group Vof already-visited peaks (the current substructure we are extracting). To see if a peak pis consistent with this substructure we select all the peaks in Vin a small subset S⊂Varound p. With this local sample we can compute a linear fit Vφ≈a0+a1R+a2φ+a3Z(3) of the subset Saround pand predict the expected Vφof the substructure in a given position. This works under the assumption that the manifold that follows the substructure is derivable and we can compute its first-order approximation locally. This prediction is already absorbing the offset in the structure position produced by its slope in a certain direction. Therefore, the criterion in the acceptance of a new peak should be stricter than that in the first adjacency step. Our limit in the resolution is the 1 km s−1bin in the Vφhistogram, and we include an extra tolerance of 0.5 km s−1. If the distance between the peak azimuthal velocity Vφ,p(Eq. (2)) and the prediction is smaller than 1.5 km s−1, we consider the peak to be consistent with the structure. We encapsulate this in the is_consistent_with function (Algorithm 2). We provide a summary of the final algorithm in pseudocode (Algorithm 3). Algorithm 2 – is_consistent_with(pad j,V) 1: S=V|Rpad j −RV|<Rf it & |φpad j −φV|< φf it & |Zpad j −ZV|<Zf it 2: f(R, φ, Z)=Vφ←linear_f it(S|pad j ) 3: if |f((R, φ, Z)pad j )−Vφ,pad j |<dthen 4: return True 5: else 6: return False 7: end if 8: Rf it =1 kpc, φf it =4 deg, Zf it =0.2 kpc, and d= 1.5 km s−1. Algorithm 3 – Breadth-first search with online interpolator 1: queue Q 2: list V 3: add p0to Vand enqueue in Q 4: while Qnot empty do 5: pit ←dequeue Q (remove and assign) 6: for all pad j adjacent to pit do 7: if pad j not in Vand is_consistent_with(pad j,V)then 8: add pad j to V, enqueue in Q 9: end if 10: end for 11: end while 12: return V 4. Results Within the 3D grid the method extracted hundreds of structures, covering 2.5 to 6kpc in Rand 30 to 60deg in φ. In R18 the arches in the SN and the radial direction were carefully characterized and matched to the groups previously studied in the literature. We rely on this matching to associate the results of our methodology with the different groups and arches. We end up with a sample of 99 structures, each one associated with one of the nine main moving groups: Arcturus, Bobylev, Hercules, Horn,Hyades,Sirius, Coma Berenices, Arch/Hat (R18, and references therein), and AC (see Anti-Centre newridge 1 in Gaia Collaboration 2021b). Each structure traces the position of the moving group along the space at a given VR. Therefore, each group of structures traces the manifold of the position of the moving groups in the (R, φ, Z,VR,Vφ) space. For each moving group we selected the largest structure as its representative. These representatives were used to study the behaviour of the groups in R−φand R−Zprojections, and are described in the following sections. In Fig. 1we show the selected groups in several neighbourhoods in the radial direction. In each arch we highlight its representative with a large black border. As explained in Sect. 3.1, our goal in this work is to be able to perform this analysis in a large extent of the disc, so we tuned the detection parameters to obtain a robust detection of the main structures in noisy regions (see R>10 kpc in Fig. 1), and this required the use of larger spatial bins. Because of this, the thin arches observed in the Gaia velocity distribution in the SN are slightly washed out. The methodology links the structures in a given VR, but does not provide the link of the arches in VR−Vφ. The complex nature of the arches formed by the moving groups at different positions in the disc and the high level of noise result in a sub-optimal global link of the arches in the velocity distribution. Therefore, we only provide a tentative manual link of the structures, based on the study of R18. In the rest of the paper we use this arch link as a qualitative tool in the analysis, but the main conclusions are based on the properties of the individual parts of the arches, which are determined by the described methodology. This linking procedure provides interesting results, different than in previous studies. For instance, the evident link of the Arch/Hat at R=9.4kpc creates an asymmetric arch shape for this structure at the SN. This could be an artefact of the detection of Arch/Hat at VR>60 km s−1or physical evidence of an unknown behaviour related to its origin. We note that tagging the groups in the SN and assuming that they remain united across the disc is a clear oversimplification. In the data we see that Sirius is formed by two arches at SN, but these arches merge into a single one at the outer parts of the disc (pannels R=8.16,9.4 kpc A116, page 4 of 16 M. Bernet et al.: From ridges to manifolds: 3D characterization of the moving groups in the Milky Way disc Fig. 1. Moving group detection in different neighbourhoods along the radial direction. For each moving group a parabolic fitting of the substructures associated with each group is included (highlighted in grey). Each moving group contains several structures, corresponding to different VRbins. The largest structure in each group is used as its representative (dots with larger black contours and the moving group name on top). There are two examples of bimodalities, which serve as evidence of the complex evolution of the arch morphology (Sirius at R=8.16 kpc and Arch/Hat at R=10.4 kpc). in Fig. 1). The same happens for Arch/Hat at R=10.4 kpc. In the simulation (Sect. 5) we observe the same behaviour for the overdensities related to the OLR. So far, this simplification is useful for the discussion and comparison to the state of the art. Moreover, the lack of data far away from the SN does not allow a robust arch characterization. In future releases an automatized arch detection will be needed to disentangle the complex orbit distribution. 4.1. Radial direction The first evidence of large-scale substructure in the dynamics of the disc was the presence of ridges in the R−Vφplane, directly related to the moving groups observed in the SN. Therefore, the first exercise we were able to do with the manifolds was to extract their subsets in the radial direction (i.e. φ=0 deg, Z=0 kpc) and plot them in R−Vφplane, coloured by VR(Fig. 2). In the top panel, we show the representative groups, tagged with their literature names. In the bottom panel, we show the rest of the structures as beams of lines that define the morphology of the corresponding moving groups. As expected, when tracing the different moving groups along Rwe observe diagonal lines in R−Vφ, matching the already known ridges. After comparison with R18, we detected the same moving groups, but we managed to extend their detection by several kiloparsecs. The results in the inner and outer part of the disc (R<6.5 kpc and R>10 kpc) are noisy due to Poisson noise. The Gaia DR3 release will improve the detection of groups in these regions, but even if we exclude this part the groups extend far beyond the range seen in other studies (see Fig. 6 in R18). This extension of the structures is due to a major improvement in the methodology and to the use of the updated astrometric Gaia EDR3 data. The lower error in proper motion and parallax increases the concentration of the moving groups in the undersampled regions. In Fig. 2we can see how the slope of the lines in radius is not constant across the different groups. In the top plot groups like Arcturus, Hercules, and Arch/Hat present slopes that are significantly steeper than Sirius and AC. However, this slope is not a common characteristic in all the parts of the arch of a group. For instance, in the bottom plot we can see that Arcturus is very steep at VR=30−40 km s−1(orange lines) and flattens for the negative part of the arch (VR=−20 km s−1, blue lines). We observe that secondary peaks have very different azimuthal and radial velocities when they start to show up at the inner radius, but end up having very similar azimuthal velocity at larger radii. This corresponds to a flattening of the arches in the velocity distribution as Rincreases. This is not observed in the simulations (Sect. 5), and could be an effect of the centroid of the distribution dominating in the case of undersampling in these regions. As for the global shape of the lines, the resonant effects of the bar and spiral arms are expected to create kinematic substructures that, from an epicyclic approximation analysis of the first-order effects, have an almost constant vertical angular A116, page 5 of 16 A&A 667, A116 (2022) Fig. 2. Azimuthal velocity of the kinematic substructures in the radial direction, φ=0◦,Z=0kpc as a function of the radius and coloured by their radial velocity. Each dashed grey line corresponds to the best fit of a structure’s constant angular momentum line. Top: structures corresponding to the main peak of a moving group, tagged with the name from the literature. Bottom: secondary peaks of the moving groups. The usual way to observe this projection is using the number of stars or the mean VRin each bin (see Fig. 1 in Fragkoudi et al. 2019). The different lines delineate the skeleton of the distribution and its complexity. momentum LZ=RVφ(Sellwood 2010;Quillen et al. 2018b). Thus, if the moving groups have a bar resonance origin, we might naively expect their azimuthal velocities to follow lines ∝R−1(dashed grey lines in Fig. 2top panel). In Ramos et al. (2018) (their Fig. 6), this trend is observed for Hercules and Hyades. When extending the analysis to a larger region, we find that the groups deviate from the lines of constant angular momentum. In Appendix Bwe compute the reduced chi-squared (χ2 ν) parameter for all the groups. The only group that is statistically well approximated globally by this Vφ∝R−1trend is Arch/Hat. We come back to this in Sect. 6. The ridges are usually studied in R−Vφdiagrams coloured by density or mean VR(see Fig. 1 in Fragkoudi et al. 2019). By doing these projections, the complexity of the moving groups (e.g. arch curvature, bi-modalities, arch disruption) is lost, offering only a partial understanding of the sample. With our A116, page 6 of 16 M. Bernet et al.: From ridges to manifolds: 3D characterization of the moving groups in the Milky Way disc methodology we can now visualize this complexity in a single plot. For instance, we can observe how the VR−Vφarch corresponding to the Arcturus moving group is a horizontal arch at R=8 kpc (structures with different VRand common Vφ), but that it curves towards the inner radius (the different structures fan out). This spreading clearly depends on the VR, which is a sign of a curved arch. Mixed with Arcturus, we are able to observe the morphology of Bobylev at VR>50 km s−1. In the previously mentioned projections, the visualization of both structures is not possible since the mean VRblends the two contributions. Beyond the detection and characterization of ridges in the radial direction, the main contribution of our method is the blind search of these kinematic structures in the three spatial dimensions at the same time. Next we focus on the representative of each moving group and study their kinematics in 3D space: the azimuth submanifold (Z=0 kpc) in Sect. 4.2 and the vertical submanifold (φ=0 deg) in Sect. 4.3. 4.2. Azimuth submanifold We first do a cut in the structures around Z=0 kpc (|Z|< 0.2 kpc) to observe the behaviour of the different moving groups as a function of Rand φ. We obtain surfaces covering up to ±25 deg (≈3.5 kpc at solar radius). This is the first time that the moving groups are traced along the Z=0 kpc plane with this completeness. In Fig. 3, apart from the already studied decrease in Vφwith R, we can see how the Arcturus, Bobylev, and Hercules moving groups present a slope in the azimuthal velocity along the azimuth, whereas the Horn,Sirius, and Arch/Hat moving groups present an axisymmetrical behaviour of the azimuthal velocity along the azimuth, as expected in an axisymmetric potential. It is interesting to quantify the variation of Vφwith φfor the structures. We evaluate this slope2∂Vφ/∂φ at a given R by restricting the structure to this Rvalue and doing a linear fitting of the surfaces in φ−Vφ. We compute the slope in the radius that minimizes the error in the fitting. The resulting values are −0.40 km s−1deg−1at R=7 kpc for Arcturus, −0.63 km s−1deg−1atR=7 kpcforBobylev,−0.50 km s−1deg−1 at R=8 kpc for Hercules, −0.04 km s−1deg−1at R=10 kpc for Sirius, and −0.01 km s−1deg−1at R=10 kpc for Arch/Hat. In Monari et al. (2019b) they study the mean angular momentum evolution in φfor Hercules in an analytical model. In the case of a co-rotation origin, they predict that the angular momentum Jφof Hercules at the solar radius must significantly decrease with increasing azimuth. Their model predicts the slope to be around −8 km s−1kpc deg−1, and they observe a similar trend in the Gaia DR2 data. Our equivalent value in angular momentum would be −4 km s−1kpc deg−1, which is smaller than the predicted value. In some parts of the disc the mean azimuthal velocity in the plane for the total 6D sample decreases with increasing azimuth at a constant radius (see Fig. 10 in Gaia Collaboration 2018b, R=8−10 kpc). This is the behaviour that we observe for Arcturus, Bobylev, and Hercules. It would be worth investigating the relative contribution of each moving group to the total sample to understand the relation between these individual groups and the total average motion, but we will postpone this to a future study. 2In our reference system a negative ∂Vφ/∂φ slope corresponds to an increase in |Vφ|with φ(a moving group with negative ∂Vφ/∂φ moves upwards in the velocity distribution with φ). Fig. 3. Mean azimuthal velocity of the representative groups in the R−φprojection, for |Z|<0.2 kpc. The contours of regions with the same velocity are shown for clarity (in white). The Arcturus, Bobylev, and Hercules moving groups present a constant slope in the variation of azimuthal velocity along azimuth, whereas the Horn,Sirius, and Arch/Hat moving groups present an axisymmetrical behaviour of the azimuthal velocity along azimuth. 4.3. Vertical submanifold We also studied the projection in the R−Zplane (|φ|<10 deg, see Fig. 4). In all the structures but Coma Berenices (see below), we observe decreasing |Vφ|values for increasing |Z|. In addition, Arcturus, Bobylev, and Hercules, the same structures that A116, page 7 of 16 A&A 667, A116 (2022) Fig. 4. Mean azimuthal velocity of the groups in the R−Zprojection, for |φ|<10 deg. The contours of regions with the same velocity are shown for clarity (in white). Coma Berenices clearly presents an increasing |Vφ|with Z, and thus strong vertical asymmetry. We measure a constant vertical slope of ∂Vφ/∂Z=−15 km s−1kpc−1. The rest of the structures show vertical symmetry. show steeper slope in azimuth, present a clear symmetry around Z=0 kpc. Arch/Hat is also very symmetric in Z. Instead, Horn,Hyades, and Sirius present a steeper decrease in |Vφ|for Z>0 kpc with respect to the other moving groups, and a more constant value for Z<0 kpc. Finally, AC does not have enough signal at this point to analyse it properly. Coma Berenices clearly presents an increasing |Vφ|with Z, and thus strong vertical asymmetry. We note that outside the range of R=[7,10] kpc, this moving group shows a change in behaviour in all the spatial projections (Figs. 2–4) possibly because our methodology is linking it to other close structures. Therefore, focusing only on the [7,10] kpc range, we measure a constant vertical slope of ∂Vφ/∂Z=−15 km s−1kpc−1, clearly different to the other groups and to the predictions from models with vertical symmetry. It would be interesting to obtain a measurement of the vertical curvature of the moving groups at each radius, as is done in the previous section with the slope in azimuth and with the vertical slope in Coma Berenices. With noisy data, each order of derivatives increases its uncertainty and the measurements we obtained were not significant enough. In the future, with better data and/or a robust analytical model to fit the curvature at all radii at the same time, this measurement could be produced. 5. Simulations In this section we apply the same methodology to the simulations. This has two main goals: evaluating the performance of our method in a case where there are no selection effects, and allowing a comparison of the data to a model where a particular and known perturbation is present. In our case we used a series of test particle simulations with 60M particles. The initial conditions and the Galactic potential are described in Romero-Gómez et al. (2015). In particular, the disc has a local radial velocity dispersion of σVR=30.3 km s−1at the radius of 8.5 kpc. We integrated the initial conditions, first in the axisymemtric potential of Allen & Santillan (1991) for 10 Gyr; then we introduced the Galactic bar potential adiabatically during 4 bar rotations (2.46 Gyr for the slow bar and 1.47 Gyr for the fast bar), to finally integrate another 4 bar rotations. The Galactic bar consists of the superposition of two aligned Ferrers ellipsoids (Ferrers 1877), one modelling the triaxial bulge with semi-major axis of 3.13 kpc and the second modelling the long thin bar with semi-major axis of 4.5 kpc. We used two simulations, where the bar rotates as a rigid body with a constant pattern speed of 30 and 50 km s kpc−1. For the slow bar, the CR was located at 7.3 kpc and the OLR at 12.2 kpc. For the fast bar, the CR was located at 4.3 kpc and the OLR at 7.6 kpc. We used the final snapshot of the simulations, and assumed that the bar is 30deg in azimuth with respect to the Sun in the direction of rotation, close to the estimations for the MW (Bland-Hawthorn & Gerhard 2016, references therein). In these final snapshots we execute the methodology described in Sect. 3and obtain an optimal detection of the moving groups in a large range of the sample. These robust results match the predictions from previous works, and validate the performance of our methodology in a known dataset. However, we found some substructures related to the centroid of the velocity distribution whose changes with azimuth, radius, and height are mostly related to the rotation curve of the model. In this section we only show the substructure related to bar resonances and ignore the rest of the groups extracted by the methodology. The selected moving groups are shown in the VR−Vφprojection in Fig. 5(as done in Fig. 1for the real data). In the simulation the moving groups also show arches in this projection, which we are able to detect at several radii. Again, we can show the groups projected in the radial direction (Fig. 6, compare Fig. 2). In Figs. 5and 6, the top panels show the structures of the fast bar model (depicting the effects of the OLR and 1:1 resonance); the bottom panels show the structures of the slow A116, page 8 of 16 M. Bernet et al.: From ridges to manifolds: 3D characterization of the moving groups in the Milky Way disc Fig. 5. Moving group detection in different neighbourhoods along the radial direction in the simulations (compare Fig. 1). Top row: fast bar model. Bottom row: slow bar model. Fig. 6. Azimuthal velocity of the kinematic substructures in the radial direction (φ=0◦,Z=0kpc) for the test particle simulations as a function of the radius, and coloured by their radial velocity. We include dashed grey lines corresponding to constant angular momentum lines as a guide. Top: structures for the fast bar model. Bottom: structures for the slow bar model. The detections of the OLR (for both the fast and slow model) and the 1:1 (only detected for the fast case) are marked on top of the lines. Here the complex morphology of the arches appears in a single image. bar model (only the effects of the OLR appear). In the following sections we analyse the fast bar model (Sect. 5.1) and the slow bar model (Sect. 5.2) in detail. As explained in the introduction, a fast bar model places Hercules near the OLR of the bar. In this model Arch/Hat can be explained as the 1:1 resonance. Instead, the slow bar model places Hercules in the CR of the bar and Arch/Hat in the OLR. Therefore, in this paper we use the term ‘Hercules-like’ for structures in the simulation that can be related to Hercules in the data (i.e. generated by OLR for the fast bar model and by CR for the slow bar model), and ‘Arch/Hat-like’ for the structures that can be related to the Arch/Hat (i.e. induced by the A116, page 9 of 16 A&A 667, A116 (2022) Appendix D: Extra projections of the simulations This is the first time that a simulation has been studied using the projection shown in Fig. 6. In general, these studies are done in projections of hVRi. In order to compare both results, in Fig. D.1 we show this projection for the simulation. In the top panel of Fig. D.1, in red, we see the Herculeslike overdensity, which steeply decreases in Vφaround R= 8 kpc. The upper part of the OLR bi-modality continues to decrease in a less steep trend, with negative hVRivalues above and positive values below the resonance. Finally, in the outer part of the disc we observe the effect of the 1:1 resonance, which shows a swap in hVRisign when crossing the rotation curve. In the bottom panel of Fig. D.1, CR should appear at R= 7.3 kpc, but we cannot see any significant structure in this region. In the outer parts of the disc we do observe the OLR placed at 12.2 kpc. Fig. D.1. Mean radial velocity of the kinematic substructures in the radial direction (φ=0◦,Z=0kpc) for the test particle simulations as a function of the radius and the azimuthal velocity. Top: Fast bar simulation. Bottom: Slow bar simulation. A116, page 16 of 16