Radial and azimuthal gradients of the moving groups in Gaia DR3: The slow/fast bar degeneracy problem

The structure and dynamics of the central bar of the Milky Way are still under debate whilst being fundamental ingredients for the evolution of our Galaxy. The recent Gaia DR3 offers an unprecedented detailed view of the 6D phase-space of the MW. We aim to characterise the dynamical moving groups across the MW disc, and use their large-scale distribution to help constrain the properties of the Galactic bar. We used wavelet transforms of the azimuthal velocity ($V_\phi$) distribution in bins of radial velocity to robustly detect the kinematic substructure in the Gaia DR3 catalogue. We then connected these structures across the disc to measure the azimuthal ($\phi$) and radial ($R$) gradients of the moving groups. We simulated thousands of perturbed distribution functions using Backwards Integration of feasible Galaxy models that include a bar, to compare them with the data and to explore and quantify the degeneracies. The radial gradient of the Hercules moving group ($\partial V_\phi/\partial R$ = 28.1$\pm$2.8 km$\,$s$^{-1}\,$kpc$^{-1}$) cannot be reproduced by our simple models of the Galaxy which show much larger slopes both for a fast and a slow bar. This suggests the need for more complex dynamics (e.g. spiral arms, a slowing bar, external perturbations, etc.). We measure an azimuthal gradient for Hercules of $\partial V_\phi/\partial \phi$ = -0.63$\pm$0.13$\,$km$\,$s$^{-1}$deg$^{-1}$ and find that it is compatible with both the slow and fast bar models. Our analysis points out that using this type of analysis at least two moving groups are needed to start breaking the degeneracies. We conclude that it is not sufficient for a model to replicate the local velocity distribution; it must also capture its larger-scale variations. The accurate quantification of the gradients, especially in the azimuthal direction, will be key for the understanding of the dynamics governing the disc. (ABR)


Introduction
The phase-space distribution function (DF) of stars in the Milky Way (MW) is a key element in understanding the structure and history of our Galaxy.A precise characterisation of this 6D phase-space (position and velocity R, ϕ, Z, V R , V ϕ , V Z ) of the MW has been made possible by the second and third Gaia data releases (Gaia DR2 andDR3, Gaia Collaboration et al. 2018a, 2023).Studying this Gaia data in various projections reveals the complex substructure within it; ridges in R − V ϕ (Antoja et al. 2018;Kawata et al. 2018;Fragkoudi et al. 2019), a wave in Friske & Schönrich 2019;Antoja et al. 2022), a bimodality in L Z − V Z (Gaia Collaboration et al. 2021;McMillan et al. 2022), the phase spiral in Z − V Z (Antoja et al. 2018(Antoja et al. , 2023;;Hunt et al. 2021Hunt et al. , 2022;;Laporte et al. 2019), and thin arches in V R − V ϕ (Gaia Collaboration et al. 2018b;Ramos et al. 2018).
Out of all of these projections, the V R − V ϕ distribution close to the solar neighbourhood (SN) is historically the most studied one (Strömberg 1946;Eggen 1965).It presents a complex configuration of overdensities -usually called moving groupsshaped as thin arches (Dehnen 1998;Famaey et al. 2005).Some of the moving groups, and also the ridges and the L Z − V R wave, have been related to the orbital resonances of the bar and spiral arms of the Galaxy (Kalnajs 1991;Raboud et al. 1998;Dehnen 2000;Antoja et al. 2011;Hunt & Bovy 2018;Hunt et al. 2018aHunt et al. ,b, 2019;;Monari et al. 2019b;Bernet et al. 2022) 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)).Despite using numerous analytical and numeric approaches to understand the structure of the phase space, we still face a significant challenge in explaining its origin in detail.
In particular, the origin of the Hercules moving group has been a subject of extensive debate.For the past two decades, it has been suggested that this moving group is connected to the resonant interactions between local stars and the central Galactic bar (e.g., Dehnen 2000;Antoja et al. 2014).According to this hypothesis, if the Sun is located just beyond the Outer Lindblad Resonance (OLR) of the bar, a Hercules-like overdensity can be generated, thus implying a bar pattern speed of Ω b ∼ 55 kms −1 kpc −1 (the short/fast bar scenario, Minchev et al. 2007;Chakrabarty 2007;Monari et al. 2017b).However, recent studies of the gas (Sormani et al. 2015) and stellar kinematics (Portail et al. 2017;Sanders et al. 2019) in the inner Galaxy point to a pattern speed of Ω b ∼ 40 kms −1 kpc −1 (the long/slow bar scenario).Following this new evidence, Pérez-Villegas et al. (2017) and Monari et al. (2019b) proved that orbits trapped in the co-rotation (CR) resonance of a long/slow bar can generate a Hercules-like overdensity in the local velocity space, although less pronounced than the one produced by the OLR (Monari et al. 2017a;Binney 2018;Hunt et al. 2018a;Fragkoudi et al. 2019).However, the combination of spiral arms and a bar (Hunt et al. 2018b), or a decelerating bar (Chiba et al. 2021;Chiba & Schönrich 2021) can create stronger Hercules-like overdensities even for a slow bar.
Indirectly measuring the pattern speed of the bar requires reproducing the moving groups at the SN through the bar resonances, something that has been mostly done qualitatively (e.g., Dehnen 2000;Pérez-Villegas et al. 2017;Hunt & Bovy 2018;Trick et al. 2021;Trick 2022, but see Asano et al. 2022;Clarke & Gerhard 2022).As the available data improved, first with RAVE (Steinmetz et al. 2006) and later on with Gaia, measurements of the pattern speed have been improved by including other regions of the galactic disk (Antoja et al. 2014;Monari et al. 2014Monari et al. , 2019b)).However, as explained above, there are a few combinations of pattern speeds and resonances that produce a similar local velocity distribution, i.e. degenerated solutions.Studying the position of the moving groups across the disc might have the potential to break this degeneracy.
In Bernet et al. (2022, hereafter B22) we presented a robust, large-scale characterisation of the moving groups across the MW disc.Our analysis revealed that the moving groups exhibit complex spatial changes, deviating from the expected lines of constant angular momenta along the radial direction and showing clear non-axisymmetries in azimuth.In particular, our study confirmed the azimuthal gradient of the Hercules moving group measured in Monari et al. (2019a).In their work, they favour the slow bar scenario using this measurement, as well as the significant azimuthal slope of the Horn, which is another moving groups at negative V R (Monari et al. 2013).Based on the predictions they obtained from perturbation theory, the gradient of L z with azimuth for a slow-bar Horn (created by the 6:1 resonance) should be non-zero, while the fast-bar Horn (created by OLR) should have a very small azimuthal gradient.
The goal of this study is to quantitatively characterise the moving groups in the phase space using the new Gaia DR3 data and improving on the methodology from B22 by including, in addition to V R and V ϕ , a novel set of measurables: the gradients of the moving groups at the SN in the radial (∂V ϕ /∂R) and azimuthal (∂V ϕ /∂ϕ) directions.The characterisation of the spatial variations of these groups allows for meaningful comparisons with theoretical models.As a particular case, we investigate the implications of kinematic substructures on the pattern speed of the MW bar, comparing the gradients obtained from the data with those of models with a simple bar, using backward integration simulations.We find that some measurables are compatible with both the fast and the slow bar models while others are incompatible with either.We also examine whether the gradients in the models truly differ between a slow and a fast bar scenario, and determine the minimum number of observables required to break the degeneracies.We observe a large disagreement between the measured gradients in the data and the simulations.
This disagreement leads us to the conclusion that simple models are insufficient to reproduce the current observations.This paper is organized as follows.In Section 2 we describe the data selection, summarize the methodology presented in B22, and explain the computation of the gradients.In Section 3, we introduce the simulations used in this analysis and how we apply our methodology to them.Section 4 describes the kinematic substructure observed in the data and the models, and compare the different trends shown by the moving groups and overdensities.In Section 5, we sweep the pattern speed Ω b and the slope of the rotation curve β in the models and compare the obtained measurables with the data.In Section 6 we study the change of the overdensities in time and azimuth in the models to understand the properties of the different structures.The implications of these results are discussed in Section 7. Finally, in Section 8 we summarise our results and list the main conclusions of this work.

Data
We used data from the Gaia Data Release 3 (Gaia DR3, Gaia Collaboration et al. 2023).From it, we selected the approximately 34 million stars with position, proper motion, parallax, and line-of-sight velocity, and used the StarHorse (Anders et al. 2022) distances.We applied an astrometric quality selection (RUWE < 1.4), a selection in parallax quality (parallax_over_error > 5), and a selection of non-spurious solutions (fidelity_v2 > 0.5, Rybizki et al. 2022).As recommended for this sample, we corrected the line-of-sight velocity for stars with grvs_mag ≥ 11 and rv_template_teff < 8500 K using Eq. 5 in Katz et al. (2022).Blomme et al. (2023) discussed the need for another correction for stars with 8500 ≤ rv_template_teff ≤ 14, 500 K and 6 ≤ grvs_mag ≤ 12.However, after the correction, a residual bias of a few km s −1 remained.Since these stars are rare, we only kept stars with rv_template_teff < 8500 K in our sample.Our final sample has 25,397,569 stars.

Methodology
In B22, we presented a novel method to extract large kinematic substructures from a dataset of stellar positions and velocities.The first step of the method is to partition the data into small spatial volumes (∆R, ∆ϕ, ∆Z) and run a Wavelet Transform (WT) on the velocity distribution (V R , V ϕ ) of each individual volume to detect the overdensities.These overdensities -the moving groups-are known to form thin arches elongated around large ranges of V R , with a slight variation in V ϕ (Ramos et al. 2018).Because of this geometry, in the second step, we slice each (V R , V ϕ ) diagram in ∆V R bins, and run a 1D WT in each V ϕ histogram.The peaks obtained from these 1D WT are connected through the configuration space to form global substructures using a modification of the Breadth-First Search (BFS, Moore 1959) algorithm from graph theory.This method has proven to be effective in detecting large kinematic substructures in the Gaia EDR3 data and test particle simulations.For more details, we refer to the original paper.
We modified a few parameters of the method with respect to the execution in B22.Firstly, we used a larger scale of the wavelet, from 2 km s −1 in B22 to 6.25 km s −1 here.The main goal of the wavelet size modification was to robustly detect the largescale substructure in regions far from the Sun, where the statistical significance of the moving groups decreases.Secondly, we increased the resolution of the grids used in the configuration and velocity spaces.This reduced the degeneracies in the linking stage of the methodology and allowed for a more precise computation of the gradients (∂V ϕ /∂R and ∂V ϕ /∂ϕ).
The execution was run in the Cloud using resources granted by the Open Clouds for Research Environments (OCRE) from the European Union.In B22 we swept the entire disc volume, but in this analysis we focused on the radial (ϕ = 0 deg, Z = 0 kpc) and azimuthal (R = 8.277 kpc, Z = 0 kpc) directions.Our new binning to compute the positions of the moving groups was: -Radial direction (R): [5,14]  The possible biases introduced by a wrong distance estimation or the size of the bins were discussed in B22.We tested the bias of the bin size using a smaller bin, and estimated the impact to be below 2km s −1 , well below the WT size.Since the coordinate transformations and binning were very similar here, we refer to those analyses.
After the execution, we obtained independent detections of each moving group at each small spatial volume, and their link across the space in the form of large-scale ridge-like structures (Fig. 1), which we analyse later on.By construction, each detection corresponds to a single V R , and each moving group is in turn assigned a set of these detections that span the entirety of its arch in velocity space.In the rest of the paper, the term structure refers to the grouping for a single V R in the entire configuration space, either in the data or the simulations.Moving group refers to a set of structures covering an entire arch in velocity space in the data, and overdensity refers to the equivalent of a moving group in the simulations.
We want to quantify how the structures change across the MW disc.However, visualising how tens of structures evolve along 2 dimensions all at once is challenging.To reduce the dimensionality of the problem, we computed the radial (∂V ϕ /∂R) and azimuthal (∂V ϕ /∂ϕ) gradients of each structure at the SN.The computation of gradients is numerically unstable, small errors in the input data propagate throughout the computation, leading potentially to big errors.To produce a robust estimation of the gradient, we fitted a parabola to the structures in the radial direction and a straight line in the azimuthal direction (dashed lines in the top plots in Fig. 1), and computed the analytical derivative of these fittings at the SN.We can then plot these gradients in the SN (left panels in Fig. 2).We propagated the uncertainties in V ϕ for each measurement (WT size, σ = 6.25 km s −1 ) using Monte Carlo sampling.The uncertainties in the measurements of the gradients in the main moving groups (Hat, Sirius, Hyades, Horn, and Hercules) are below 1 kms −1 kpc −1 in ∂V ϕ /∂R, and below 0.05 kms −1 deg −1 in ∂V ϕ /∂ϕ.

Simulations
The Gaia data is unprecedented in its quality and quantity.As we explained in the previous section, it even allows us to compute the spatial gradients of the velocity structures.Running realistic simulations that match the quality of the data is complex and expensive.For this reason, we used the simpler, but many times faster, backwards integration (BI) method that we explain in Sect.3.1 to explore exhaustively the parameter space of different models.In Sect.3.2 we describe the specific method of detection of structures in the simulations.

Setup and bar potential
The BI technique (Vauterin & Dejonghe 1997;Dehnen 2000, also Hunt & Bovy 2018 for a more recent example) allows us to approximate the response of the DF to a bar perturbation by integrating the orbit of stars, all starting at a certain point (R, ϕ) in configuration space, but with different velocities on a grid in the (V R , V ϕ ) plane.According to the collisionless Boltzmann equation, the density of the DF inside an infinitesimal phase-space volume remains constant.However, the DF can only be measured in finite volumes.The BI technique assumes that the mean value within the local volume is the same as its central value, regardless of how the volume is deformed along the orbital evolution.This means that we trace back a single orbit at the central point of each volume to compute the DF before the perturbation.It is important to note, however, that these local volumes undergo significant stretching and kneading, potentially challenging the assumption that the mean is equal to the central value.In other words, the measurable DF, also sometimes called the coarse-grained DF, does not obey the collisionless Boltzmann equation.Moreover, when smoothly switching on the perturbation amplitude up to its plateau, the separatrixes will move in time whilst virtually always involving non-adiabatic behaviour on their surface, implying a dependence on the switch-on choice.Nevertheless, the assumptions made here are still extremely useful, especially for studying the dynamical effect of the perturbation for relatively few dynamical times.
In practice, we follow the model presented in Dehnen (2000).To validate our approach, we refer to existing literature: Monari et al. (2017a) modified the setup proposed by Dehnen (2000) to ensure complete phase-mixing and calculated the present-day DF using the pendulum approximation.Additionally, Trick et al. (2021) ran forward test particle simulations in an identical setup.Both studies produced results in agreement with the BI technique, affirming the reliability of our computations and the low impact of the assumptions presented in the previous paragraph.
We use the Dehnen (1999) distribution function to model the stellar disc before the bar formation where R E , Ω(R E ), and L c (E) are radius, circular frequency, and angular momentum of a circular orbit with energy E. The ana-lytical form of the DF allows for an efficient computation.Σ(R) and σ R (R) are defined as (3) where the values of the constants Σ 0 , σ 0 , R s , and R σ (Tab. 1) are the same as in Dehnen (2000).This DF is used extensively in the literature (e.g.Hunt et al. 2018a;Monari et al. 2019b), providing an initial setup consistent with previous analysis.Possible bias induced by the selection of the DF parameters are discussed later on.
The axisymmetric potential is assumed to be such that the circular velocity is given by where v 0 is the circular velocity at the Solar radius R 0 , and β characterises the shape of the rotation curve.
Also following Dehnen (2000), we model the bar potential as a simple quadrupole where Ω b is the pattern speed of the bar, ϕ b is angle of its long axis and R b its size, which we fix to be 80% of the corotation radius2 , given by where Ω 0 ≡ v 0 /R 0 denotes the circular frequency at the Sun.The bar amplitude is initialized at A b (t) = 0 for t = 0, increases smoothly for t < t 1 as and remains at a fixed A f for t > t 1 until the end of the simulation at t = t 2 .We express t in units of the bar rotation period T b ≡ 2πΩ b .Finally, the amplitude constant is defined as which depends on the dimensionless bar strength α, defined as the ratio of the radial forces from the bar's quadrupole and the axisymmetric power-law background at R 0 along the major axis of the bar (Dehnen 2000).
We present two fiducial models, which aim to reproduce the velocity distribution of the SN in a MW with a fast bar and a slow bar.In the Slow Bar Model (SBM) we included a bar with a pattern speed of Ω b = 39 kms −1 kpc −1 .The Fast Bar Model (FBM) has a pattern speed of Ω b = 56 kms −1 kpc −1 .The Sun is at R ⊙ = 8 kpc, and the bar is at ϕ b = −30 deg with respect to the Sun (Bland-Hawthorn & Gerhard 2016).The parameters of the fiducial simulations can be found in Table 1.The orbit integration was performed using AGAMA (Vasiliev 2019).
For consistency, we tested the effect of a significant modification of the parameters of the DF.We increased individually the fiducial values of Σ 0 , σ 0 , R s , and R σ by a 25% and assessed

Symbol
that the change on the position of the overdensities is below 2km s −1 .In addition, we tested the relative error in the gradients with the same increase of 25%.The relative error in radial gradients is below 2%, while the relative error in azimuthal gradients is below 5%.These shifts are even less significant on higher order resonances.With this analysis, we show that the measurements are robust to significant variations of the input parameters.We remark that the local velocity dispersion of the used model (σ R = 17.6 km s −1 ) would correspond to a relatively younger population of stars of the MW (1 − 2 Gyr, Robin et al. 2022).We tested the measurables obtained with a larger dispersion (σ R = 35 km s −1 ) which would correspond to an older population (6 − 7 Gyr, Robin et al. 2022) and the variations are well within the errors estimated in Table 2. Finally, we computed the gradients with a McMillan (2017) potential, a Ferrers (1877) bar and a Quasi-isothermal DF (Binney & McMillan 2011), and observed an error of 5% in radial gradient, and 20% in azimuthal gradient.Therefore, we conclude that the choice of DF has an impact on the final models, but the obtained measurables are still dominated by the global physics of the resonances.
We also note that the used models are 2D for simplicity, ignoring the contribution of the vertical motion.It is known that the position of the moving groups depends on the vertical component (B22).However, the studied sample is very dominated by stars close to the plane, and we tested that different vertical sizes in the selection induce biases below 1 km s −1 in the positions of the moving groups.Therefore, we conclude that the vertical substructure will not bias our measurements, thus allowing for a fair comparison with 2D models.

Arch detection
In these simulations, we applied the methodology described in Sect.2.2.The first step of this methodology is to detect overdensities in the velocity distribution.For the data, the WT has proven to be optimal for this task (e.g.Chereul et al. 1998;Skuljan et al. 1999;Antoja et al. 2008;Kushniruk et al. 2017;Ramos et al. 2018;Lucchini et al. 2023b).However, for the simulated velocity distributions, we find that in the case of the SBM the detections with the WT are not robust enough.This is because the overdensity caused by CR is not as prominent as when observed in the data or the FBM (Binney 2018;Hunt et al. 2018a).Luckily, the velocity distributions produced with BI have potentially "infinite" resolution3 , and can be considered twicedifferentiable.Taking advantage of that, we used the statistic defined in Contardo et al. (2022) that estimates the "gappiness" of each point for differentiable distributions.We constructed the statistic as follows.
For each point in the V R − V ϕ grid, we computed the gradient vector ∇p and the Hessian matrix H of the density on the velocity distribution.Then we calculated the projected second derivative tensor ΠHΠ , with where I is the identity.The maximum eigenvalue of ΠHΠ is the "gappiness" estimator of the point.Notice that "gappiness" is the opposite of "overdensityness", therefore the minimum eigenvalue of ΠHΠ gives information on the position of the overdensities in the velocity distribution.For more details, we refer to the original paper (Contardo et al. 2022).We note that it would be ideal to apply the same methodology to the Gaia data.However, we aim to reach distant regions of the disc, where the Poisson noise is very significant and we do not have the differentiability conditions to apply this method.We tested the difference between both methods in the fiducial models, and found it to be below 1 km s −1 in the main overdensities, well below the WT size used in the data.When the overdensities are not sharp enough horizontally, the different methods may diverge.However, these overdensities are not studied in this paper, since we focus on the main structures.For less significant structures, a careful discussion on the bias induced by the method would be needed.
In the simulations, we applied the same methodology as in the data, described in 2.2, but using the "overdensityness" in V R − V ϕ .This reveals a complex velocity field, rich in overdensities, the equivalent in the simulations of the moving groups.To then measure the gradient of each overdensity, we first computed auxiliary velocity distributions in the radial direction (between 7 and 9 kpc, in steps of 0.1 kpc) and the azimuthal direction (between 20 and 40 deg in steps of 1 deg).For each BI-generated V R − V ϕ map, we detected the overdensities present and used the same methodology that we used for the data to find the set of structures, groupings across configuration space with a single V R , that compose them.
In Sect.5, we explore the parameter space of the models.For this particular case, the amount of data we generate would be too large to deal with.Instead, we fitted a parabola in V R − V ϕ for each overdentisty of interest and selected the maximum as their representative.We note, however, that sometimes one overdensity might get split into two independent ones after crossing a certain radius, for instance, through the mechanism explained in Dehnen (2000) for the OLR.To avoid confusion, we select the representative obtained from the larger of the resulting overdensities after the split.Overall, this representative selection is feasible thanks to the stability of the overdensities in the simulations, which remain bounded along the configuration space, but it is infeasible for the data due to the unstable shape of the moving groups once we move a few kpc away from the SN.

Data vs Fiducial Models
In this section, we explore the phase-space and the spatial gradients of the kinematic structures in the SN for the Gaia DR3 data and the models.

Global ridges in radius and azimuth
Figure 1 shows the structures obtained with our methodology, which trace the skeleton of the phase space.We show their V ϕ in the radial (ϕ = 0 deg, left panels) and azimuthal (R = 8.277 kpc, right panels) directions.This projection allows us to study all the structures at the same time.We see how the different lines organise in bundles related to each moving group.In the inner disc, the high kinematic temperature and the low number of stars observed blur these structures.Despite the usefulness of our methodology elsewhere, in this region it clearly struggles to obtain reliable detections.Thus, we ignore the detections with V ϕ > 80R − 740 km s −1 in this analysis (translucent white region in Fig. 1, left panels).
In the first row of Fig. 1, we show the V ϕ of the moving groups as a function of R and ϕ, selecting, from all the structures in each moving group at different V R , only the one that covers a largest area (accounting for both the radial and azimuthal directions).This subset selection is used only for visualization purposes.In the second row, we show all the structures that are obtained (about 1500).The new results show an improvement in the resolution and extent of the structures with respect to B22.Good examples of this are the extended range of detection of Hercules in the outer disc, the new structure below Sirius in the high |V ϕ | part of the inner disc, and the azimuthal extension of Hyades and Sirius up to ±40 deg.
In the second row, we observe how the different bundles of lines reveal the rich ridge-like pattern.In the outer parts of the disc, we observe two clear ridges (AC1, and AC2).These were first detected in Gaia Collaboration et al. (2021), where it was hypothesised that, due to the likely weaker effects of the bar resonances at these radii, they could be due to spiral arms and/or the interaction with external satellites.Below these ridges we observe the Hat, which in the radial direction, is traced robustly from the SN to R = 13.5 kpc.In the azimuth direction, we are only able to detect the Hat in a small angular range.
Sirius is the most dominant structure in the radial direction.It is detected in almost the entire studied range robustly, both in the radial and azimuthal directions.Below Sirius, we observe Hyades and the Horn, which are especially well detected in azimuth.Hyades shows a robust linear behaviour within ±40 deg, while the Horn shows a negative slope for ϕ < 0 and a flat profile for ϕ > 0.
Finally, we focus on Hercules, which is known to be formed by three thin arches (Gaia Collaboration et al. 2018b;Asano et al. 2020).In this execution, due to the use of a WT with a larger scale, two of these arches (A8 & A9 in Ramos et al. 2018) appear joined as Hercules1, while the other arch (A10 + A11) is detected and labelled as Hercules2.We note that the azimuthal velocity of Hercules1 has a clear linear increase with azimuth.As for the radial behaviour, with our methodology, both in B22 and here, we observe a ridge that flattens as we move to the inner parts of the disk.In B22, we discussed that this flattening could be related to the influence of the centroid of the distribution in the peak detection of the methodology.Interestingly, in this work, a new set of structures have appeared in the inner disk, which we label Hercules_In.This structure appears as a symmetric parabola (peak of |V ϕ | at V R = 0 km s −1 ), in good agreement with Fig. 4 of Dehnen (2000), where x 1 (2) orbits dominate (his Fig. 7).In the R − V ϕ and V R − V ϕ maps, we observe a potential continuity between Hercules1 and Hercules_In, thus, we manually linked both sets of structures at each V R (thick dashed line in Fig. 1, bottom left panel).In the rest of the paper, we will consider these grouped structures as Hercules1.We find that, in  doing so, the resulting structure is more coherent with the results obtained in the literature simulations (e.g.Fig. 13 in Hunt et al. 2019).

Gradients of the moving groups and overdensities
Figure 2 shows the moving groups detected at the SN in the Gaia DR3 data (left panels), and the overdensities in the fiducial models (middle and right panels), colored by their gradients ∂V ϕ /∂R (top) and ∂V ϕ /∂ϕ (bottom).Thus, we are including information of the large-scale shape, i.e. how the velocities of the structures change with R and ϕ.In the left panel, the big dots correspond to the largest structure of each moving group, i.e. the structures shown in the top panel of Fig. 1, which we included for visualization purposes only.In the middle and right panels, the crosses correspond to the representatives of the overdensities in the simulations, described in Sect.3.2.These are the representatives used in Sect. 5.
As explained in the introduction, Hercules is a moving group commonly associated with bar resonances.In the slow bar scenario, the Hat has been related with the OLR (e.g.Monari et al. 2019b).In the fast bar scenario, it is associated with the 1:1 resonance (Dehnen 2000).In Fig. 2 we represent the regions of these moving groups with dashed black boxes.We computed the mean and standard deviation of the gradients of the structures within the boxes in the data and the models (Tab.2).We do this to capture the gradients of the entire moving group, not just at a specific V R , as it allows a more robust comparison with the simulations.
In Hercules1, we measured a similar azimuthal gradient throughout the entire moving group ∂V ϕ /∂ϕ ∼ − 0.63 kms −1 deg −1 (Table 2), compatible with both models within the error bars.In ∂V ϕ /∂R, we detected two subsegments for Hercules1: in the V R < 50 km s −1 part of the arch, we measured ∂V ϕ /∂R∼29 kms −1 kpc −1 , while in the V R > 50 km s −1 part we obtained ∂V ϕ /∂R∼25 kms −1 kpc −1 .The overall average slope of Hercules1, without separating into segments, ends up being ∂V ϕ /∂R= 28.1 ± 2.8 kms −1 kpc −1 .The corresponding values of ∂V ϕ /∂R for the SMB (36.5 kms −1 kpc −1 ) and FBM (40.5 kms −1 kpc −1 ), however, are both about 10 kms −1 kpc −1 larger than the gradient measured in the data.In Section 5 we explore this discrepancy in more detail.The gradients of Hercules2 are significantly different than Hercules1 but no counterpart is found in the models.
In the Hat, we observe different ∂V ϕ /∂R values within various parts of the same arch.It is important to note that the data detections in the Hat are not as reliable compared to the Hercules data, due to the low number of sources in the border of the distribution.Despite this, the gradients seen in both models are compatible with the data distribution.
The gradients in Sirius, Hyades, and Horn behave smoothly within large parts of the moving groups, with sudden breaks in some places.These breaks are a projection of the pattern already observed in Fig. 1, where on the large scale the ridges break and merge.
We also explored the dependence of these gradients with the integration time.In Fig.
A.1 and Table A.1 we show the equivalent of Fig. 2 and Table 2 for t 2 = 14 bar laps, in contrast to the t 2 = 4 bar laps in the fiducial models.We observe that the gradients in Hercules and Hat are very similar at short and long integration times, thus indicating that time is not playing a major role in the value of these measurables, similar to Dehnen (2000), where he pointed out that the velocity of Hercules does not strongly depend on time.In the FBM, the Horn structure is maintained for long integration times, while it is not seen in the SBM.
Finally, in the negative V R part of the distribution we note a resemblance, both in position and gradients, between the overdensities in the models and Hat, Sirius, and Horn.In Section 6, we study these overdensities in the models to explore an alternative hypothesis on the origin of these moving groups.

Model Parameter Exploration
In the literature, the resonant origin of the moving groups is usually tested by placing a given resonance (or set of resonances) in the SN, and comparing the resulting V R −V ϕ distribution with the data.We have already discussed in the introduction that different pattern speeds create resonances compatible with the data (slowfast bar degeneracy).Other parameters can also shift the position of the overdensities in the models, thus leading to extra degeneracies, e.g. the slope of the rotation curve.In this section we use the azimuthal and radial gradients as large-scale measures and sweep different parameters of the fiducial models, to test if we can indeed break some of these degeneracies.We start by studying how the properties of each simulated structure change with the model parameters and identify properties that may strongly depend on the model (Sect.5.1).Secondly, we compare the properties of the models with those of the data and try to identify the parameter region that offers a better match.
Apart from the model parameters that we want to explore, there is a bunch of kinematic structures from these models to compare to the data.For this analysis we focus on Hercules and Hat as possible moving groups associated to the bar's effects.Possible links between other moving groups and the overdensities seen in the models are discussed in Sect.6.

Parameter Exploration
The parameter space to explore is huge and multidimensional (Table 1).Below we detail each of the parameters and specify which ones we investigate:  (Dehnen 2000).Therefore, we leave it fixed to the value of the fiducial models.
-Ω p : We explored a parameter range of ±2 kms −1 kpc −1 around the slow and fast bar pattern speeds, within which the V ϕ of the overdensities in the models remain compatible with Hercules and Hat.
β (rotation curve slope): We explored this parameter around the flat rotation curve (0 ± 0.1).This range covers the realistic values for the MW rotation curve (but see discussion in Sect.7).
In Fig. 3, we show how the measurables V ϕ , V R , ∂V ϕ /∂ϕ, ∂V ϕ /∂R (different rows) of the Hercules-like (left) and Hat-like (right) overdensities behave, as a function of Ω p (horizontal axis in all panels) and β (vertical axis in all panels).Each pixel of each panel represents a different Galaxy Model.We cover the parameters around the two fiducial models: SBM for the first and third columns, and FBM for the second and fourth columns.
For the interpretation of some of these maps it is important to have in mind that, when a resonance moves inwards in the disc, the associated orbits reach the SN with a smaller V ϕ , and vice-versa.The first row of Fig. 3 shows that the azimuthal velocity |V ϕ | of both the Hercules-and Hat-like groups decreases as Ω b increases.This is expected, since the position of all resonances move inwards in the disc as the bar rotates faster.On the other hand, we see opposed behaviours with β: while the |V ϕ | of Hercules-like overdensities decreases with β, that of the Hat-like overdensities increases.This can be explained by the same argument as above once we consider the fact that changes in β affect the position of the resonances differently whether they are inside R 0 or outside.An increase in β decreases the slope of the azimuthal frequency curve, thus sending the resonances inside R 0 inwards, and the resonances outside R 0 outwards.Conversely, if β decreases the frequencies change faster with radius, which pulls the resonances closer to R 0 .
Regarding ∂V ϕ /∂ϕ (third row in Fig. 3), it shows an almost linear decrease with Ω b and β for both Hercules-like4 and the Hat-like overdensities across the parameter space for the slow bar (columns 1 and 3).In the fast bar case, the gradients of the Hat-like overdensity remain at a fixed value and shows little to no variation (rightmost column).
In contrast, the ∂V ϕ /∂R of Hercules-like structures increases almost linearly with both Ω b and β, while in Hat-like structures, ∂V ϕ /∂R increases with Ω b and decreases with β.At first glance, this difference in behaviour might be related to the effect β has on V ϕ , which we described above.In this sense, we can interpret all this information as follows: if the resonance moves inwards, the slope observed at the Sun increases, whether it moves inwards because of β (increasing for inner resonances and decreasing for the outer ones) or by increasing the pattern speed of the bar.Moreover, the slope seems to be proportional to Ω b regardless of the type of resonance or its resonant radius.It is difficult to obtain an intuition of these slopes beyond the empirical measurements, but these trends are a good milestone for an analytical modelling of these measurables.
A visual comparison between the slow and the fast bar models reveal that, to first order, the overdensities present a similar behaviour.Specially in the Hercules-like overdensities, the contours of all the measurables are parallel, indicating a consistent directional change in the parameter space.The only significant qualitative difference is the steepness of this change (distance between white contours) in the gradients.This similarity complicates the task of breaking degeneracies between the slow and fast bar scenarios.
In Appendix C we study the differences between the fiducial SBM and FBM measurable for Hercules and Hat (Tab.2) and the galaxy models around the FBM.The goal is to evaluate the minimal parameters we should consider to distinguish between the models, in case the data matched the models.We conclude that the best combination to constrain the bar parameters is to consider the measurements for both Hercules and Hat, and that including the azimuthal gradient of both structures brings the statistical significance over 8σ in almost the entire parameter space.

Comparison with Gaia Data
The black contours in Fig. 3 show the statistical deviation between each galaxy model and the data, i.e. the difference w.r.t the mean value in Table 2 in units of standard deviation.For instance, a 1σ deviation means that the difference between the model and the mean gradient of the corresponding moving group equals the standard deviation of the data reported in Table 2.
Since we are exploring the parameter space around the fiducial models, to first order the results are similar to what we have described in Section 4.2.In some parts of the parame- Third row: ∂V ϕ /∂ϕ of the overdensity.In black, the statistical deviation to the measurements in Gaia DR3 (Table 2).Fourth row: same as third row, with ∂V ϕ /∂R.
ter space, the azimuthal gradient is matching the one seen in the data (0σ curves).In the slow bar Hercules-like overdensity (leftmost column) the azimuthal gradient matches the data for Ω b ∼40 kms −1 kpc −1 and β∼0.025.Opposite to that, the azimuthal gradient of the fast bar Hercules-like overdensity (second column) matches the data for Ω b ∼54 kms −1 kpc −1 and β∼ − 0.05.
The radial gradients, on the other hand, do not match the values of the data in any part of the parameter space explored.Thus, in general, we do not find any combination of Ω b and β that reproduces all the measures in the data.To obtain a match in radial gradient we would need to decrease the parameters under Ω b ≪ 37 kms −1 kpc −1 , away from the match in V ϕ , and β ≪ −0.1, well below realistic models of the MW.Similar analysis in the other azimuths (see Figs. B.1 and B.2) lead to similar differences in radial gradients, always above 3σ when comparing models and data.This confirms that our models are sub-optimal to reproduce the large-scale velocity distribution for any parameter combination.

Azimuth and Time Exploration
In Section 4.2, we pointed out the resemblance, both in position and gradients, between the overdensities in the models and the moving groups Hat, Sirius, and Horn in the data.The goal of this section is to understand these structures in the models and the potential implication on the origin of the mentioned moving groups.

Origin of the overdensities in the models
We aim to study the global shape of these overdensities in the negative V R region of our models, which is the region of the Horn that at the same time includes Sirius and Hat.To simplify our analysis, we focus on a subset of structures, arbitrarily those with a radial velocity V R = −60 km s −1 .In Fig. 4, we present the azimuthal velocity (V ϕ ) of these structures across the azimuthal range of ±90 degrees with respect to the long axis of the bar.The models have 180-deg symmetry, so we just need to cover half of the disc.In the top panels, we study the structures in the fiducial models.We see that all the overdensities in the SN (black vertical line) belong to a few global structures, that link the different overdensities in the SN in an entire turn around the galaxy (forming double helical shapes).The dashed red lines in the top panels of Fig. 4 represent these links5 .These connections could already be suspected in Fig. 2 of Dehnen (2000), but to our knowledge, it is the first time that it is characterised in detail.Since this analysis is done at very short timescale (4 bar laps, i.e. 0.64 Gyr for the slow bar, and 0.44 Gyr for the fast bar) in the second and third rows of Fig. 4 we do the same analysis for t 2 = 10 and 40 bar laps, respectively.We observe two phenomena: the predicted resonant regions appear (CR and OLR for the slow bar and OLR and 1:1 for the fast bar), and the global helix winds, flattening in ϕ, except at the resonance boundaries, where the lines get squeezed.To understand the evolution from one row to another, in Fig . A.2 we show the velocity of these structures with time.We observe that the V ϕ of most of the structures evolves in time, asymptotically approaching the boundaries of one of the resonant regions.
In the last row of Fig. 4 we can clearly see the nodes of the resonant regions, which are not aligned with the bar major axis as one would expect.We computed the position of the nodes for positive V R = 50 km s −1 and they are located in the opposite site of the major axis.Therefore, at V R = 0 km s −1 the nodes will be aligned with the bar as expected, indicating that the displacement of these nodes is probably due to the eccentricity of the orbits.A more detailed description of these nodes is beyond the scope of this work.In summary, in the models we have two types of overdensities: resonances, and transient arches.The position of the resonance boundaries stabilize relatively quickly (∼8 bar laps) and remain fixed, as previously anticipated in Dehnen (2000).The transient arches are joined trough an entire turn around the disc forming a global structure with helical shape.In addition, this helical shapes ends at the resonance boundaries, joining all the overdensities.These transient arches are the result of introducing the bar potential in an axisymmetric DF, and thus contain temporal information of the model.

Comparison with Gaia Data
As an exercise, we now assume that, in the MW, the Hat, Horn, and Sirius moving groups are indeed related with the same global helical structure.In this scenario, the Hat would be formed in a resonance boundary and Sirius and Horn would in turn be the imprints of the ongoing phase-mixing.In Fig. 5, we show the detection of structures in the data for V R = −60 km s −1 in the radial and azimuthal direction.Following the assumption that they are related, we compute the fitting of a linear helix in azimuth (grey line on the right panel).We obtain a global slope for this structure of ∂V ϕ /∂ϕ= −0.177 kms −1 deg −1 .Although the fitting results are poor, especially for Sirius (B), it does match the slope of Hat (A) and the parts of Horn (C) at ϕ > 0.
Back to the models, we observed that the space between the transient arches decreases with time (Fig. 4).That is, the ∂V ϕ /∂ϕ of the transient arches in azimuth decreases with time.We note that the diagrams of this figure show aspects similar to typical phase-mixing in a frequency-angle plot (in our case equivalent to Ω ϕ -ϕ plane), where the lines get closer and have smaller slopes with time (e.g.Li & Widrow 2021;Frankel et al. 2023;Darragh-Ford et al. 2023, for the case of the phase spiral).In Fig. 6 we compute this slope of the transient arches for the slow (green points) and fast (orange points) bar model at different integration times.We see that, if t is expressed in bar laps, the evolution of the slope with time is the same for both pattern speeds.It seems therefore that in this case the phase-mixing times are modulated not only by the orbital frequencies but also by the bar as well.
Again, assuming that the structures in the model are the same as the moving groups, we see that the measured slope for the data (see above) gives a timescale estimate of the MW bar of t ≈ 0.6 Gyr.This is evidently a very short timescale compared with the usual MW bar age estimations (≳ 8 Gyr, e.g.Grady  ).In the following section, we discuss the implications of this short timescale prediction.

Discussion
Radial gradients.The main disagreement between the models and the data is observed in the radial gradients of the moving groups.Our measurement of the radial gradient of Hercules ∂V ϕ /∂R= 28.1 ± 2.8 kms −1 kpc −1 is compatible with recent results such as the one in Ramos et al. (2018), of ∂V ϕ /∂R= 26.5 kms −1 kpc −1 , but does not match the ones in the models by a significant amount (Table 2).In previous studies, they found a good agreement between the data and simulations in the case of the FBM.For instance, in Antoja et al. (2014) the detection of Hercules6 in the RAVE data at different radius matched the empirical predictions of test particle models.Revisiting this work, we note that the radial slope measured with RAVE (∂V ϕ /∂R= ∼33±10 kms −1 kpc −1 ) was significantly larger than the one found in this work but with a large measurement error.In our case, even when modifying the slope of the rotation curve (row 4 in Fig. 3) or the angle of the bar (row 4 in Figs.B.1 and B.2) we are still far from obtaining a match between the models and the data, with a disagreement above 3σ.However, it could be possible that the MW circular velocity curve is not well approximated by a single power-law.In fact, as discussed in Antoja et al. (2022) where we examined various models and observations for the MW, not all cases could be well fitted by a single power law model.For instance, the model by McMillan (2017) shows raising and decreasing trends with the transition at R ≈ 8 kpc.Quantifying this effect with the current set-up requires changing some of its pieces, like the DF used.However, we can use the formalism of Ramos et al. (in prep.) to assess the impact of different rotation curves on the slopes.
It may be that a more sophisticated axisymmetric model, the effect of spiral arms (Hunt et al. 2018b), self-gravity, Giant Molecular Clouds (GMCs), external perturbations, and/or, the scenario that is currently gaining the most traction, a decelerating bar (Chiba et al. 2021;Chiba & Schönrich 2021) are nonnegligible effects that need to be taken into account.A clear next step using the framework that we have constructed is to compute the radial gradients of the moving groups in a decelerating bar and check if we are able to reproduce the trends measured in Gaia.
Using a combination of perturbation theory and the pendulum approximation (Monari et al. 2019b), in Monari et al. (2019a) they found a slope of ∂L z /∂ϕ = −8 kms −1 kpc deg −1 at a radius of 8.2 kpc, which corresponds to ∂V ϕ /∂ϕ= −0.96 kms −1 deg −1 , when looking at the mean V R of the model of a slow bar, which roughly agrees with the observations.They also observe the Horn variation with azimuth, which we confirmed in Figs. 1 and 2. In Monari et al. (2019a), they claim that in the fast bar regime, a Horn moving group with an OLR origin would show a flat azimuthal structure.In contrast, in the slow bar regime, a Horn created by the 6:1 resonance would have a significant slope.In the case of Hercules, they predict a distinct azimuthal slope in the slow bar regime (CR origin) and a less pronounced slope in the fast bar regime (OLR origin).
Our measurements in the fiducial models show evidence that both the SBM and the FBM present a significant gradient in azimuth for the Hercules-like overdensity (Table 2) for short and long integration times.We also confirm that the azimuthal structure of the fast bar Horn-like structure is flatter than the Hercules-like.However, due to the absence of high-order moments in our simple bar model, we can not generate a Horn-like overdensity (6:1 resonance) in the slow bar regime.
We observe that for Hercules-like overdensities, both azimuthal gradients are still within the error bar of the data (even when exploring the parameter space around the fiducial models, see row 3 in Fig. 3).We would need to reduce the observational errors at least by a half to obtain a significant similitude (or difference) between the models and the data (Tab.2).Finally, we note that the measure in Lucchini et al. (2023a), with a smaller error, is very compatible with our models with a fast bar.
Azimuthal shape of Hercules.In an ideal, long-lived barred galaxy, the resonances form steady state closed islands in the V ϕ − ϕ diagrams (Fig. 4) with a certain azimuthal periodicity, depending on the resonance order (e.g.Tsoutsis et al. 2009;Michtchenko et al. 2018).The features in the velocity distribution related to these islands will then follow periodic curves in azimuth.Therefore, either they have null slope or they must contain a minimum and a maximum within a period (Bolzano's theorem, Bolzano 1817).If Hercules is related to a resonance, it should follow one of these periodic curves in the V ϕ − ϕ space.Thus, it should contain a maximum and a minimum in −90 ≤ ϕ ≤ 90 deg.
In our measurements, we observe Hercules to follow a linear trend covering up to 60 deg and we do not observe a maximum nor a minimum.There are two possibilities: the maximum and the minimum are located in a part of the disc that we have not ob-served yet, or the resonance regions are not stationary (closed), indicating that the MW disc is far from phase-mixed.Future spectroscopic surveys (WEAVE, 4MOST, SDSS-V/MWM) covering larger fractions of the disc will allow us to detect Hercules in a larger azimuthal range, ideally going beyond the long axis of the bar, and confirm its connection to this non-axisymmetric component of the Galaxy.
Origin of Hat, Sirius and Horn.In Sect.6, we used the moving groups Hat, Sirius and Horn (in particular, their azimuthal gradient) to time the stage of phase-mixing induced by the growth of the bar.Our estimate of ∼0.6 Gyr contradicts by far the current measurement of the age of the MW bar, ≳ 8 Gyr (Grady et al. 2020;Sanders et al. 2022).Thus, given the initial assumption, it is unlikely that our hypothesis is true.However, these timescales below one Gyr would be compatible with more recent events, such as the Sagittarius (Sgr) impact, and a paradigm of ongoing phase-mixing of the MW disc.In fact, other studies have determined the phase-mixing timescales in the planar velocities (V R -V ϕ ) following a perturbation from a galactic satellite.These studies find similar but slightly larger times.For instance, in Antoja et al. (2022) we found times of 0.8−2.1 Gyr based on the separation between peaks of the wave in V R -L Z .Prior to that, Minchev et al. (2009) found a start of phase-mixing time of 2 Gyr based on a separation between moving groups of 20 km s −1 .For this same velocity separation, from our Fig.6 we would obtain times of around 1 Gyr, i.e. smaller than in Minchev et al. (2009), and which also depend on the pattern speed of the bar, thus indicating a different underlying phase-mixing mechanism.Although the studies mentioned are more specific for an external perturbation than the present work, they do not explore the azimuthal dimension and are still simplified models.More tests must be done to establish the relation between certain moving groups and Sagittarius, both with idealized models (Antoja et al. 2022) and taking into account the effect of self-gravity (Darling & Widrow 2019).

Possible caveats of this work
The main contribution of this work is the development of novel techniques to quantify the large-scale kinematic substructure in both the data and models.The final goal of this quantification is to perform a direct discrimination between the likelihood of the models.However, this task has a clear obstacle: a proper modelling of the data.The models used in this work, despite being intentionally simple, depend on a large set of parameters.Throughout the present work, we tested the impact on the measurables of each of the parameters of the model (Tab. 1) and the 'observational' methodology (vertical size of the selection, overdensity detection, gradient fitting. . .).We assessed that the bias produced by each of the tested factors independently is well below the precision we can measure in the actual data.However, there is a non-negligible possibility that a combination of biases in the model and/or the methodology is the cause of the clear disagreement with the data.Regardless, as mentioned above, the simplicity of the model itself is most likely the main contributing factor, as the interaction with a satellite, spiral arms, and/or a slowing bar would play a major role that has not been tested in this work.

Summary and conclusions
In this work we provide a detailed analysis of the kinematic substructure of the MW disc in the Gaia DR3 data across 7 kpc in radius and 80 deg in azimuth.We compared the results in the data to a suite of BI simulations, and we found that the overdensities detected in the fiducial models, as well as their gradients to a lesser degree, are compatible between them, and that this similarity is maintained even in the parameter space around the fiducial models.Our main findings and conclusions are: -We detect new ridges in the outer disc, propose a new radial profile for Hercules, and precisely characterise Hyades, Horn, and Hercules ∼70 deg along the disc.-We observe a robust slope in azimuth for Hyades and Horn, previously unexplored.-The radial gradient of Hercules is ∂V ϕ /∂R= 28 ± 2.3 kms −1 kpc −1 , and the azimuthal gradient is ∂V ϕ /∂ϕ= −0.63 ± 0.13 kms −1 deg −1 .
-The radial gradient of Hercules in the BI models disagrees with the data with a significance above 3σ for all the explored parameters (Ω b , β, and ϕ b ).We conclude that more complex models, which take into account spiral arms, selfgravity, external perturbations, and/or a slowing-down bar, are required to reproduce the observations.Despite this, we acknowledge that this problem has a high complexity and we could be underestimating some biases in the analysis.-The azimuthal gradient of Hercules in the data is compatible with both the slow and the fast bar models.For these simple models, we require at least half of the error in the data to disentangle between them.-Our analysis points out that a robust determination of the pattern speed using our type of analysis requires fitting the azimuthal velocity and the azimuthal gradient of at least two moving groups.-We study the azimuthal slope of the phase-mixing structures induced by the growth of the bar in the models, and asses that for identical Galaxy models, the evolution of the slope only depends on the pattern speed of the bar.-We explore the possibility of a phase-mixing origin for Sirius and Horn in the context of the growth of the bar, which leads to a too low estimate of the bar age (∼0.6 Gyr).We conclude that a different hypothesis is required to explain these structures.
With this work we confirm the potential of our methodology and of the Gaia data to provide a quantitative description of the MW kinematic substructure.In the near future, the use of spectroscopic surveys to extend the range of the Gaia 6D sample, specially towards the tip of the MW bar, will allow us to disentangle the resonant origin of Hercules.Finally, this project makes it clear that simple models can not capture the full complexity of our Galaxy, specially in the large-scale.It emphasizes the need of exploring more advanced models, and reminds that an optimal model will have to go beyond matching local patterns.

Fig. 1 .
Fig. 1.Azimuthal velocity of the kinematic substructures as a function of radius (ϕ = 0 • , left panels) and azimuth (R = 8.24 kpc, right panels), and coloured by their radial velocity.The vertical error in these measurements is the size of the WT, σ = 6.25 km s −1 .Top: main structure for each moving group, tagged with the name from the literature.The dashed black lines illustrate the fitting used to compute the gradients in the SN.Bottom: all the structures detected by the methodology.

Fig. 2 .
Fig. 2. Moving group and overdensities detection, colored by their ∂V ϕ /∂R(top panels) and ∂V ϕ /∂ϕ(bottom panels) at each V R .The dashed boxes represent the region of the velocity distribution associated with Hercules and Hat.Left: Gaia DR3 SN moving groups, tagged by their literature names.Middle: SBM detection of the overdensities.Right: FBM detection of the overdensites.

-
R: In each model, we sweep R between 7 and 9 kpc, in steps of 0.1 kpc to compute ∂V ϕ /∂R.-R ⊙ : It can be fitted by external measurements (Gravity Collaboration et al. 2022), so we maintain the fiducial model value R ⊙ = 8 kpc.v 0 : We keep this parameter fixed but we checked that changing it even by a considerable amount (±15 km s −1 ) does not change significantly our results (variations in the radial slopes of about 1 kms −1 kpc −1 ).-ϕ: In each model, we sweep ϕ between −10 and 10 deg, in steps of 1 deg to compute ∂V ϕ /∂ϕ.-ϕ b : It can also be fitted by external measurements, though there is less consensus.In this section, we keep the value of ϕ b = −30 deg used in the fiducial model.In Appendix B we repeat the analysis for ϕ b = −20, −40 deg and see that, to first order, all the trends are the same.-α (bar strength).The strength of the bar hardly changes the position of the resonances

Fig. 3 .
Fig.3.Parameters of the Hercules-like (first and second columns) and Hat-like (third and fourth columns) overdensities representatives in the models.The parameter space is swept around the fiducial models β and Ω b .First row: V ϕ of the overdensity.Second row: V R of the overdensity.Third row: ∂V ϕ /∂ϕ of the overdensity.In black, the statistical deviation to the measurements in Gaia DR3 (Table2).Fourth row: same as third row, with ∂V ϕ /∂R.

Fig. 4 .
Fig. 4. V ϕ of the overdensities at V R = −60 km s −1 for an entire turn around the model of the galaxy.Left: Results for the SBM.Right: Results for the FBM.Top: Structures for t 2 = 4 bar laps.We observe how the arches seen in the fiducial models are part of the same global structure, linked through a 180 deg symmetry (red dashed lines).Middle: Structures for t 2 = 10 bar laps.The resonant regions start to be clear, and the structures in the middle are flatter and closer between them.Bottom: Structures for t 2 = 40 bar laps.The resonant regions are clearly defined, and the middle region is crowded with flat structures.The vertical dashed lines mark the location of the Sun.

Fig. 5 .
Fig. 5. Detection of structures in the data at V R = −60 km s −1 (extracted from the bottom panel of Fig. 1).Left: Radial position of the structures.Right: Azimuthal position of the structures.

Fig. 6 .
Fig. 6.Azimuthal slope of the transient arches at each integration time, for a Slow and a Fast bar.The black line shows the curve −0.72t −1 , which approximates both curves.The dashed black line corresponds to the slope fitted in the data 5), and the intersection returns the time estimation.The dashed red line marks the end of the growth of the bar (t = 2 bar laps).
Fig. A.1.Overdensities detection in the simulations, for t = 14 bar laps.Colored by their ∂V ϕ /∂R(top panels) and ∂V ϕ /∂ϕ(bottom panels) at each V R .The dashed black boxes represent the region of the velocity distribution associated with Hercules and Hat.Left: SBM detection of the overdensities.Right: FBM detection of the overdensites.

Fig. C. 2 .
Fig. C.2. Similar to Fig. C.1 but showing the differences between the fiducial FBM and the galaxy models around it when using only the Hercules-like overdensity (top) or Hercules-and Hat-like (bottom).

Table 1 .
Parameters of the fiducial models.

Table 2 .
Mean slope and the dispersion within the box (see Fig.2) in the data and both fiducial models.