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

Context. The details of the e ﬀ ect 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 φ − R plane and thin arches in the V φ − V R plane, 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 speciﬁc 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-ﬁrst 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.


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. 2019b) 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. 2018b;Khanna et al. 2019;Laporte et al. 2019Laporte et al. , 2020. The latest releases of the Gaia mission (Gaia Collaboration et al. 2018a, 2021b 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 et al. 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 Kawata et al. 2018;Fragkoudi et al. 2019).
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 kms −1 kpc −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 (R OLR ≈ 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 the bar might be longer and slower than previously thought . In this case, the OLR would be placed further out (R OLR ≈ 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 (R CR ≈ 6 kpc). Pérez-Villegas et al. (2017) and Monari et al. (2019b) 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. 2018a). However, Hunt et al. (2018b) 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. (2019a) 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 V R and 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 V R −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, V R , 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 kms −1 kpc −1 ) and a slow (Ω b = 30 kms −1 kpc −1 ) bar, both to test our methodology and to compare it to the data.
The Gaia DR3 (Gaia Collaboration et al. 2022b) catalogue includes a larger and updated sample of radial velocities (33 M of stars, Gaia Collaboration et al. 2022a), 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 Section 2 we describe the observational data we used. In Section 3 we introduce the methodology we developed. In Section 4 we show the results of the application of the method to Gaia EDR3 data. In Section 5 we present and analyse the simulations. In Section 6 we compare the results from the data and the simulations, and with previous results in the literature. Finally, in Section 7 we list the main conclusions of this work.

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 photo-geometric 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: 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 kms −1 for 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 (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 et al. 2019), φ = 0 o and Z = 0.0208 kpc (Bennett & Bovy 2019). For the solar motion, we use U = 11.1, v circ + V = 248.5, W = 7.25 km s −1 (Schönrich et al. 2010;Reid & Brunthaler 2020).

Method
The data described in the previous section contains the 6D variables of position (R, φ, Z) and velocity (V R , V φ , V Z ) of the stars. Inside small volumes, cuts in (R, φ, Z), the moving groups appear as well-defined overdensities in the velocity distribution V R -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 Section 3.2).

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: 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 kms −1 for the bin size used. For each pixel we construct the velocity distribution (V R , V φ ) diagram of the stars in it as a 2D histogram with bins of 1 kms −1 (see background histogram in Fig. 1). R18 showed that the overdensities form thin arches elongated around large ranges of V R , 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 V R , and uncontrollable correlations between V R and V φ (movement along the arch). To avoid this, we slice each V R − V φ diagram in vertical columns (bins in V R ), and run a 1D WT in the V φ histogram of each column: radial velocity (V R ): −100 -100 kms −1 in steps of 10 kms −1 , V R bin = ±15 kms −1 around 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 V R and 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] kms −1 , with steps of 1 kms −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 p inherits the spatial position from the pixel, V R from the position of the radial velocity bin, and V φ from the result of the WT detection. Therefore, a peak has the coordinates (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 V R bin and if their distance in V φ is smaller than 4 kms −1 . In R18, the maximum slope found in a moving group is 33 kms −1 kpc −1 . In our grid the step is 0.04 kpc, which translates into a maximum change of ≈ 1 kms −1 between adjacent pixels. This 4 kms −1 limit 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 p to the queue 1 Q; 2. Select (dequeue) the top peak p top of the queue Q; 3. Visit all the peaks p ad j adjacent to p top . 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.
Given an initial peak p 0 , 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 p 0 to V and enqueue in Q 4: while Q not empty do 5: p it ← dequeue Q (remove and assign)

6:
for all p ad j adjacent to p it do 7: if p ad j not in V then 8: add p ad 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 V of already-visited peaks (the current substructure we are extracting). To see if a peak p is consistent with this substructure we select all the peaks in V in a small subset S ⊂ V around p. With this local sample we can compute a linear fit of the subset S around p and 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 kms −1 bin in the V φ histogram, and we include an extra tolerance of 0.5 kms −1 . If the distance between the peak azimuthal velocity V φ,p (Eq. 2) and the prediction is smaller than 1.5 kms −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).

Results
Within the 3D grid the method extracted hundreds of structures, covering 2.5 to 6 kpc in R and 30 to 60 deg 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, Algorithm 3 -Breadth-first search with online interpolator 1: queue Q 2: list V 3: add p 0 to V and enqueue in Q 4: while Q not empty do 5: p it ← dequeue Q (remove and assign) 6: for all p ad j adjacent to p it do 7: if p ad j not in V and is_consistent_with(p ad j , V) then 8: add p ad j to V, enqueue in Q 9: end if 10: end for 11: end while 12: return V Bobylev, Hercules, Horn, Hyades, Sirius, Coma Berenices, Arch/Hat (R18, and references therein), and AC (see Anti-Centre newridge 1 in Gaia Collaboration et al. 2021a). Each structure traces the position of the moving group along the space at a given V R . Therefore, each group of structures traces the manifold of the position of the moving groups in the (R, φ, Z, V R , 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 − Z projections, and are described in the following sections.
In Fig. 1 we 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 Section 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 V R , but does not provide the link of the arches in V R -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.4 kpc creates an asymmetric arch shape for this structure at the SN. This could be an artefact of the detection of Arch/Hat at V R > 60 kms −1 or 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 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.

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 V R (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 R we 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. 2 we 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 V R = 30−40 kms −1 (orange lines) and flattens for the negative part of the arch (V R = −20 kms −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 R increases. 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 momentum L Z = 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. 2 top 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 B we compute the reduced chi-squared (χ 2 ν ) parameter for all the groups. The only group that is statistically well approximated globally by this V φ ∝ R −1 trend 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 V R (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 methodology we can now visualize this complexity in a single plot. For instance, we can observe how the V R -V φ arch corresponding to the Arcturus moving group is a horizontal arch at R = 8 kpc (structures with different V R and common V φ ), but that it curves towards the inner radius (the different structures fan out). This spreading clearly depends on the V R , which is a sign of a curved arch. Mixed with Arcturus, we are able to observe the morphology of Bobylev at V R > 50 kms −1 . In the previously mentioned projections, the visualization of both structures is not possible since the mean V R blends 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.

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 R and φ. 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 Figure 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 slope 2 ∂V φ /∂φ at a given R by restricting the structure to this R value 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 kms −1 deg −1 at R = 7 kpc for Arcturus, −0.63 kms −1 deg −1 at R = 7 kpc for Bobylev, −0.50 kms −1 deg −1 at R = 8 kpc for Hercules, −0.04 kms −1 deg −1 at R = 10 kpc for Sirius, and −0.01 kms −1 deg −1 at R = 10 kpc for Arch/Hat.
In Monari et al. (2019a) 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 kms −1 kpc deg −1 , and they observe a similar trend in the Gaia DR2 data. Our equivalent value in angular momentum would be −4 kms −1 kpc 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 et al. 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.

Vertical submanifold
We also studied the projection in the R − Z plane (|φ| < 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 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, 3, 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 kms −1 kpc −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.

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 σ V R = 30.3 kms −1 at 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 semimajor 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 kms 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 30 deg 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 Section 3 and 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 V R -V φ projection in Fig. 5 (as done in Fig. 1 for 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. 5 and 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 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 ex-plained 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 1:1 for the fast bar model and by the OLR for the slow bar model).
In Fig. 6 it is clear that the global position of the overdensities does not follow exactly the lines of constant angular momentum predicted for small epicyclic amplitudes (V φ ∝ R −1 ). In the radial direction the curvature of the OLR_bott structures is opposite in the two models. In addition, the fact that the structures merge at some radius is clear evidence that they follow a trend that is different from ∝ R −1 . The first-order prediction for resonances is suboptimal for V φ values far from the circular velocity. In Fig. 6 we observe how the radial slope of the structures differs more from the prediction at small and large V φ values.

Fast bar model
In the fast rotating bar simulation (Ω b = 50 kms −1 kpc −1 ) we are able to detect substructures related to the OLR and the 1:1 resonance with our methodology. As expected, both structures present an arch shape in the V R -V φ diagram (Fig. 5, top). The complete sampling and the simplicity of the simulation allows us to trace these arches unequivocally and observe their morphology in ranges up to 6 kpc in the radial and azimuthal directions, characterizing with strong robustness not only their position in the velocity space, but also the trend of the central and outer part of the arch separately. We note that the simulation has been integrated for 4 bar rotations, which places this model in the regime of a young fast bar, as defined in Monari et al. (2019a).
In the top part of Fig. 6, the bi-modality of the OLR and the 1:1 resonance are clearly observed. For the OLR, at R < 7 kpc we observe a Hercules-like arch (OLR_bott), extending to a maximum V R of 40 kms −1 (orange line at top). In this same region we also observe a symmetric arch at the top of the distribution labelled OLR_top, with a maximum in V R = 0 kms −1 (green line at top). At R between 7 and 8 kpc, the V R -negative part of the bottom arch (V R around −10 kms −1 ) continues to decrease in |V φ |, and the positive part of the arch (V R > 20 kms −1 ) merges with the top arch (OLT_top) to form a unique structure covering the whole radial velocity range. This unique arch, labelled OLR_comm, has its maximum at V R = −40 kms −1 (blue line at top). The arch configurations at R = 7, 9, and 11 kpc can also be seen in the top row of Fig. 5.
The signature of the 1:1 resonance (Arch/Hat-like structure) is less prominent than in the case of the OLR. In Fig. 5 (top row, R = 11 kpc) we observe a bi-modality in the histogram, with the upper part of the distribution being more prominent at negative V R values and the lower part more prominent in the positive V R range. In the R-V φ diagram (Fig. 6, top panel) we see that at the inner radius we mostly detect negative V R , and at the outer parts mostly positive V R , giving more evidence of this bi-modality. Due to the small prominence of the resonance, we are not able to detect it when it is located in the centre of the distribution.
As we did with the real data, we exploit the 3D extent of the structures and show their trends in the azimuth submanifold. In Fig. 7 we show the R − φ projection of the azimuthal velocity for different parts of the arches (compare Fig. 3). The black lines in these panels show the azimuthal slope of V φ in each radius (right vertical axes). The OLR_comm structures present a discontinuity around R = 8 kpc. In Fig. 6 this can be seen as the sudden drop in V φ in the green-turquoise lines at R ∼ 8 kpc (V R = −20, −10, 0 km s −1 ). A similar discontinuity is present in two panels of the OLR_bott structure in Fig. 7. The peak detection algorithm guarantees a minimum distance between peaks for robustness. Therefore, when two arches merge, the algorithm stops detecting one of them a short time before the merge. Even so, a few spurious detections can act as a bridge for the algorithm and and make it join the two structures. Since the arches are merging, the algorithm described in Sect. 3.2 (Algorithm 3) detects a continuity and the structures are detected together.
It is important to remember that Fig. 7 shows a young fast bar, which is known to introduce azimuth correlations, especially around the OLR (see Fig. A1 in Trick et al. 2021). For most of the structures, we see an azimuth slope of V φ that depends on radius. Interestingly, some structures have negative ∂V φ /∂φ, such as Hercules in the real data, and others have a positive slope. We observe three main structures with different patterns: -Rows 1-4, and 7. Upper part of the OLR bimodality. Before crossing the rotation curve (R ≈ 8 kpc), the ∂V φ /∂φ slope of the structure is constant along all the V R values. For R > 8 kpc, a correlation appears between ∂V φ /∂φ and V R . This is a sign that the OLR_comm arch shown in Fig. 5 (top row, R = 9 kpc) is moving to negative V R values in the V R -V φ diagram as φ increases (when V R is more positive, its |V φ | decreases faster with φ). -Rows 5, 6, and 7 (R < 8 kpc). Hercules-like part of the distribution. The slope in azimuth decreases as R increases. The contour lines merge when approaching the bar's long axis (φ b = −30 deg). -Rows 8 and 9. Arch/Hat-like part of the distribution. The slope in azimuth is constant along the different detected parts of the arch, and tends to flatten for increasing R values.

Slow bar model
In the slow rotating bar simulation (Ω b = 30 kms −1 kpc −1 ), we are able to detect the Arch/Hat-like overdensity caused by the OLR along a wide range of radius values, covering up to 6 kpc (Figs. 5 and 6, bottom panels). In the radial direction we observe three different structures related with the OLR. Two of them have negative V R and the other has positive V R . The negative V R structures (OLR_top and OLR_bott) present a clear bi-modality in V φ . In the V R -V φ diagram (Fig. 5, bottom), the OLR_top forms a flat arch and the OLR_bott forms a curved arch. In the positive V R regime (orange-red lines), we observe one single structure (OLR_comm) compatible with OLR_top at R = 10 to 11 kpc, which continues to decrease its azimuthal velocity with R in a constant slope and merges OLR_bott from R = 12 kpc to the outer parts of the disc. Again, we can exploit the 3D extension of the manifolds to study how these groups are evolving in the R − φ projection. In Fig. 8 (as in Fig. 7, but for the slow bar model) we observe two main structures, corresponding to the different parts of the OLR bi-modality (whose upper part corresponds now to the Arch/Hat-like group). Rows 1-3 are the top part of the OLR bi-modality. It presents a constant slope in azimuth (∂V φ /∂φ = −0.1 kms −1 deg −1 at R = 12 kpc), flattening as R increases, similar to the Arch/Hat-like part observed in Fig. 7 for the fast bar. Rows 4-7 show the bottom part of the bi-modality. The slope is less constant, and we observe the same decrease in R as in the Hercules-like part of the fast bar model.
We do not detect overdensities caused by CR. In general, the expected Hercules-like moving group caused by CR is less prominent than when formed by the OLR (Binney 2018;Hunt et al. 2018a), especially if higher order modes are missing from the bar model. In addition, if the velocity dispersion of the disc is too large in this region, we could have less trapping and less resolution. Another explanation could be that the length of the  Fig. 7). In this simulation the shape of the arches is maintained along φ (same slope for all V R values), with a constant displacement to bigger |V φ | as φ increases. bar in our models does not favour CR trapping. In the future, we aim to perform the same analysis with other models to analyse this effect (e.g. Sormani et al. 2022).
In Fig. D.1 we show the V φ -R projection coloured by mean V R . In this figure we observe a very faint red overdensity in the CR region, but much weaker compared to the rest of the detected resonances in the fast and the slow bar models.

Discussion
Simulations that contain only a bar perturber offer us the possibility to characterize its effects in a robust way. With this characterization, we can then compare the manifolds extracted from the simulation with those seen in the data. As explained in the introduction, one of our goals with the characterization of the manifolds is to search for kinematic observables that distinguish resonances from short/fast or long/slow bars.
Radial direction. First, we focus on the slope, shape, and evolution of the moving groups in the data and the simulation along the radial direction (Figs. 2 and 6).
The global gradient of the structures in the radial direction deviates from the first-order prediction for resonances (curves following V φ ∝ R −1 ) in the data and the simulations, as expected in more realistic cases, even if they have a resonance origin. In the simulations we see that the lines are less curved than the simple prediction and some parts of the OLR_bott have even the opposite curvature. This is also observed in the data, where all the moving groups are less steep in the inner parts of the disc than the prediction, and Hercules and Arcturus present the same opposite curvature pattern as the OLR_bott.
We see in the V φ -R projection coloured by mean V R of the slow bar scenario of our simulations (Fig. D.1) that the corotation overdensity is not significant enough. This is discussed in the Introduction and in the previous section, and could be related to the specificities of the model considered here. In the fast bar scenario a Hercules-like structure is formed at R = [6, 8] kpc, which shows an arch shape with a maximum in V R = 40 kms −1 (Fig. 5, top left panel), consistent with the Hercules moving group in the data. In this same radial region the top part of the OLR forms an arch in negative V R , which can be related to Horn. Finally, from R = 8 kpc towards the outer parts of the disc, the OLR forms a single arch with a maximum in V R = −20 kms −1 . A similar arch, present at all the radial velocities, was already observed in the models by Bovy (2010). In the data this region (R > 9 kpc, V φ < −170 kms −1 ) contains very few stars, and we did not detect any group there.
The other feature commonly associated with bar resonances is the Arch/Hat moving group. For a fast bar, it can be explained by the 1:1 resonance trapping region, and for a slow bar it matches the position of the OLR. In the data, the number of stars with R > 10 kpc is low. Therefore, the quality of the shape characterization of the Arch/Hat moving group is poor. Even so, in the panel R = 10.4 kpc in Fig. 1 our method finds a tentative arch split around V R = 10 kms −1 . This would match the split we detect between the OLR_comm and OLR_top structures in the slow bar model (Fig. 6, bottom panel at R = 11 kpc). In the fast bar simulation, the prominence of the 1:1 resonance is low and we can only detect it in regions far from the centroid of the distribution. Even with this limited detection, we do not observe a bi-modality in the negative V R region, which would be a way to discriminate between resonances.
Azimuth submanifold. In both simulation models, for the Arch/Hat-like moving groups we observe a constant positive slope of |V φ | in azimuth which tends to flatten (get closer to 0 kms −1 deg −1 ) as R increases. In these Arch/Hat-like groups, the slope of the azimuthal velocity in azimuth does not depend on V R (Figs. 7 and 8). In the data (Fig. 3), the Arch/Hat group is very noisy. Therefore, as discussed in the previous paragraph it is still complicated to use the Arch/Hat velocities for a final relation to a specific resonance.
In Fig. 7 we can observe the different parts of the OLR resonance in the fast bar. The Hercules-like overdensity (R = [6, 8] kpc) shows a constant positive slope of |V φ | in azimuth common along all V R . In opposition, the trapping region of the resonance OLR_comm has a slope that depends on the V R . This can be interpreted as the arch moving along V R in the V R -V φ diagram when moving in azimuth, which is observed at the OLR in other simulations of a young bar far from phase-mixed (Dehnen 2000;Bovy 2010;Trick et al. 2021). Since we have the complete information of the moving groups in the data, we can reproduce the same analysis of the moving group slope at different V R for some moving groups in the data. In Figures C.1 and C.2, we show this slope for Hercules and Hyades, respectively. We do not observe this dependence in V R for any of the groups. In the velocity distribution this can be observed as the moving group arches increasing their |V φ | along the azimuth, but no displacement along V R .
Vertical submanifold. For the resonces of the bar, the vertical displacement of V φ should be dominated to first order by the vertical potential of the Galaxy (Al Kazwini et al. 2022). In our results the vertical curvature of all the moving groups (except for Coma Berenices) at a common radius is very similar, thus matching the analytical prediction. This means that this data is a good candidate to constrain the 3D structure of the potential. To second order, we could try to distinguish different curvatures of different resonances. In Al Kazwini et al. (2022) the displacement of the resonances at Z = 1 kpc with respect to the galactic plane is measured to be 8 kms −1 kpc −1 for the corotation, 6 kms −1 kpc −1 for the OLR, and 4 kms −1 kpc −1 for the 1:1 resonance. Our maximum resolution, given by the V φ histogram at each pixel, is 1 kms −1 kpc −1 . Therefore, disentangling which resonances create each moving group with the vertical information is beyond our current capabilities.
In the vertical behaviour of the moving groups, a clear outlier is Coma Berenices. In Quillen et al. (2018b) (also Monari et al. 2018;Laporte et al. 2019), the Coma Berenices moving group is observed to be present only at negative galactic latitudes, showing evidence of incomplete vertical phase mixing. With the new EDR3 data and our methodology, we are able to detect the group in a larger extension and at positive and negative galactic latitudes, but it does show a clear vertical asymmetry in its azimuthal velocity. We measure a constant vertical slope of ∂V φ /∂Z = −15 kms −1 kpc −1 . In the plots of a phase spiral coloured by V φ , it is seen that there must be a correlation between Z and V φ (higher |V φ | values at positive Z). This slope in Coma Berenices could be a projection of this correlation.

Conclusions
We have sampled, with the Gaia EDR3 6D data, the manifolds tracing the main moving groups in the SN along the (R, φ, Z, V R , V φ ) space, in an automatic way. We have revealed the skeleton of the velocity distribution in a multidimensional space that we can then explore along the radial direction, and characterize in the azimuth and vertical submanifolds. This methodology was successfully tested with two simulations of the effects of a dynamically young bar. We were able to observe and quantify the spatial evolution of the observed moving groups in a large range of about 3 kpc around the sun. Our main results and conclusions are the following: -The azimuthal velocity of the moving groups in the radial direction does not follow lines of constant angular momentum, deviating from our naive first-order prediction for resonances. In the simulations, resonant structures also deviate from this simple prediction, demanding more complex analytic predictions. -The spatial evolution of the moving groups is complex. The moving group configuration observed in the SN is not maintained throughout the disc. The relative position between the arches and their curvature changes across space, and the different moving groups split and merge several times. This is expected in a context of bifurcating orbital families, for example in the case of resonances. -In our slow bar simulation we observe a bi-modality created by the OLR in the outer parts of the disc. This bi-modality is also observed in the Arch/Hat moving group in the data. This intriguing agreement could favour the slow bar scenario, and opens the possibility of a test with future data. -The Arcturus, Bobylev, and Hercules moving groups present a positive slope of their V φ location with the azimuth. We measure this slope to be −0.50 kms −1 deg −1 at R = 8 kpc for Hercules. -The azimuthal velocity of the Horn, Sirius and Arch/Hat moving groups presents an axisymmetrical behaviour. In both our simulations we observed a small azimuthal gradient in V φ in the Arch/Hat-like structures, although it approches 0 as R increases. This could be related to the young bar model we are using. -The vertical curvature of the moving groups is similar at the same R. These curvatures are dominated by the gravitational potential to first order, independently of the observed resonance. However, we find that the Coma Berenices group deviates from this behaviour, which points to a different dynamical origin that deserves further investigation. -In the fast bar simulation a correlation between ∂V φ /∂φ and V R is observed for the OLR trapping region. The region where this correlation is observed in the simulation (R > 9 kpc, V φ < −170 kms −1 ) is poorly sampled in Gaia EDR3, but this could potentially be used to give information on the pattern speed of the bar with better data.
Spiral arms, resonances with the bar, accretion events, and possibly other effects can contribute to the present phase space distribution from which we obtain our observables. Disentangling all the contributions of these dynamical processes is difficult. In this work we have shown the complexity of the phasespace structure that even a single mechanism (namely the bar) can produce. Our methodology allows a quantitative and robust measurement of the observed phase space substructure to be extracted, which can be then compared and/or fit to different models.
In this paper we have developed the methodology for the study of the disc with the Gaia data, and its formulation is easily generalizable. The same approach can be exported to other substructures in astrophysics (e.g. blind search for streams and shells in the halo) and even to datasets outside this field.
In the following months Gaia DR3 will revolutionize once again our field of research. The new 6D sample contains about 33M stars, covering a significantly larger region of the disc. During the review process, this new data was already being used in Lucchini et al. (2022) to reveal four new moving groups candidates in the SN. With this work, we are ready to process the Gaia DR3 data and extract its full potential. The 3σ error of the slope is shown in a translucent region around the line. The group has the same slope at all V R . This corresponds to a vertical displacement of the moving group in the V R -V φ diagram.
with an optimal fitting. Finally, the shape of the rest of the structures is absolutely different from the lines of constant angular momentum.

Appendix C: Azimuth submanifold along a group
In Figure 3 we show the azimuth submanifold of a representative for each group. In this appendix we extend this analysis and show how each part of the Hercules and Hyades arches is evolving. In addition, since we have already shown these plots, we include the extra information on the slope of V φ in azimuth at each radius. With this information, we can compare these groups to the resonances studied in the simulations (Figs. 7 and 8).
The Hercules moving group (Fig. C.1) presents a stable and constant slope of V φ in azimuth at the centre of the sample (R = [6.5, 8] kpc). In the limits of the sample, where the significance of the group is smaller compared to the sample, this slope tends to flatten for the regions with V R ≈ 0 and to become more steep in the large V R parts of the arch.
The Hyades moving group (Fig. C.2) also presents a stable behaviour at all V R . In the negative V R end of the arch (top two rows in the figure) the group has a very low significance, thus leading to a noisy detection.