© The Authors 2024

1. 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_ϕ, and V_Z) of the MW has been made possible by the second and third Gaia data releases (Gaia DR2 and DR3, Gaia Collaboration 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 L_Z − V_R (L_Z = RV_ϕ, Friske & Schönrich 2019; Antoja et al. 2022), a bimodality in L_Z − V_Z (Gaia Collaboration 2021; McMillan et al. 2022), the phase spiral in Z − V_Z (Antoja et al. 2018, 2023; Hunt et al. 2021, 2022; Laporte et al. 2019), and thin arches in V_R − V_ϕ (Gaia Collaboration 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 groups – shaped 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. 2018a,b, 2019; Monari et al. 2019a; 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. 2019, 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 km s⁻¹ kpc⁻¹ (the short and fast bar scenario, Minchev et al. 2007; Chakrabarty 2007; Monari et al. 2017a). 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 km s⁻¹ kpc⁻¹ (the long and slow bar scenario). Following this new evidence, Pérez-Villegas et al. (2017) and Monari et al. (2019a) proved that orbits trapped in the co-rotation (CR) resonance of a long and slow bar can generate a Hercules-like overdensity in the local velocity space, although it is less pronounced than the one produced by the OLR (Monari et al. 2017b; 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 one to reproduce 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 disc (Antoja et al. 2014; Monari et al. 2014, 2019a). However, as explained above, there are a few combinations of pattern speeds and resonances that produce a similar local velocity distribution, that is to say 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. (2019b). 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 investigated 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 (BI) 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 examined whether the gradients in the models truly differ between a slow and a fast bar scenario, and we determined the minimum number of observables required to break the degeneracies. We observed 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 organised as follows. In Sect. 2 we describe the data selection, summarise the methodology presented in B22, and explain the computation of the gradients. In Sect. 3, we introduce the simulations used in this analysis and explain how we applied our methodology to them. Section 4 describes the kinematic substructure observed in the data and the models, and we compare the different trends shown by the moving groups and overdensities. In Sect. 5, we explain how we swept the pattern speed Ω_b and the slope of the rotation curve β in the models and compare the obtained measurables with the data. In Sect. 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 Sect. 7. Finally, in Sect. 8 we summarise our results and list the main conclusions of this work.

2. Data and methodology

2.1. Data

We used data from the Gaia Data Release 3 (Gaia DR3, Gaia Collaboration 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. (2023). 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⁻¹ 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.

We transformed the observables into cylindrical phase space coordinates using R₀ = 8.277 kpc (GRAVITY Collaboration 2022)¹, Z_⊙ = 0.0208 kpc (Bennett & Bovy 2019), and U_⊙ = 9.3, V_⊙ + V_c(R₀) = 251.5, and W_⊙ = 8.59 km s⁻¹ from the proper motion of SagA* (Reid & Brunthaler 2020) and its radial velocity (GRAVITY Collaboration 2022). The reference system is right-handed with ϕ oriented contrary to the disc rotation and thus V_ϕ < 0 for most stars. We set the origin at the Sun’s azimuth (ϕ_⊙ = 0).

2.2. 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⁻¹ in B22 to 6.25 km s⁻¹ here. The main goal of the wavelet size modification was to robustly detect the large-scale 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] kpc in steps of 0.01 kpc,ΔR = ±0.24 kpc around each centre;

• Azimuthal direction (ϕ): [ − 50, 50] deg in steps of 0.1 deg,Δϕ = ±2.4 deg around each centre;

• Vertical direction (Z): Not swept, ΔZ = ±0.24 kpc around Z = 0 kpc;

• Radial velocity (V_R): [ − 120, 120] km s⁻¹ in steps of 2 km s⁻¹,ΔV_R = ±15 km s⁻¹ around each centre.

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 2 km s⁻¹, 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.

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.

We want to quantify how the structures change across the MW disc. However, visualising how tens of structures evolve along two 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⁻¹) 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 km s⁻¹ kpc⁻¹ in ∂V_ϕ/∂R, and below 0.05 km s⁻¹ deg⁻¹ in ∂V_ϕ/∂ϕ.

Fig. 2. Moving group and overdensities’ detection, coloured by their ∂Vϕ/∂R (top panels) and ∂Vϕ/∂ϕ (bottom panels) at each VR. 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.

3. Simulations

The Gaia data are 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, 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.

3.1. 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. (2017b) 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

$$\begin{array}{c}\hfill f(E,L)\propto \frac{\mathrm{\Sigma }(R_E)}{{\sigma }_R^2(R_E)}exp[\frac{\mathrm{\Omega }(R_E)(L-L_c(E))}{{\sigma }_R^2(R_E)}],\end{array}$$(1)

where R_E, Ω(R_E), and L_c(E) are radius, circular frequency, and angular momentum of a circular orbit with energy E. The analytical form of the DF allows for an efficient computation. Σ(R) and σ_R(R) are defined as

$$\begin{array}{cc}& \mathrm{\Sigma }(R)={\mathrm{\Sigma }}_{0}exp(-R/R_s),\hfill \end{array}$$(2)

$$\begin{array}{cc}& {\sigma }_R(R)={\sigma }_{0}exp(-R/R_{\sigma }),\hfill \end{array}$$(3)

where the values of the constants Σ₀, σ₀, R_s, and R_σ (Table 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. 2019a), 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

$$\begin{array}{c}\hfill v_c(R)=v_0(R/R_0{)}^{\beta },\end{array}$$(4)

where v₀ is the circular velocity at the Solar radius R₀, and β characterises the shape of the rotation curve.

Also following Dehnen (2000), we model the bar potential as a simple quadrupole

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_b(R,\varphi )=A_b(t)cos(2(\varphi -{\varphi }_b-{\mathrm{\Omega }}_{b}t))\times \{\begin{array}{cc}-{(R_b/R)}^3,\hfill & \mathrm{for}R\ge R_b,\hfill \\ {(R/R_b)}^3-2,\hfill & \mathrm{for}R\le R_b,\hfill \end{array}\end{array}$$(5)

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 radius², given by

$$\begin{array}{c}\hfill R_{\mathrm{CR}}\equiv R_0({\mathrm{\Omega }}_b/{\mathrm{\Omega }}_0{)}^{1/(\beta -1)},\end{array}$$(6)

where Ω₀ ≡ v₀/R₀ denotes the circular frequency at the Sun.

The bar amplitude is initialised at A_b(t) = 0 for t = 0, increases smoothly for t < t₁ as

$$\begin{array}{c}\hfill A_b(t)=A_f[\frac{3}{16}{\xi }^5-\frac{5}{8}{\xi }^3+\frac{15}{16}\xi +\frac{1}{2}],\xi =2\frac{t}{t_1}-1,\end{array}$$(7)

and it remains at a fixed A_f for t > t₁ until the end of the simulation at t = t₂. We express t in units of the bar rotation period T_b ≡ 2πΩ_b. Finally, the amplitude constant is defined as

$$\begin{array}{c}\hfill A_f=\alpha v_0^{2/3}(\frac{R_0}{R_b}{)}^3,\end{array}$$(8)

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₀ 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 km s⁻¹ kpc⁻¹. The Fast Bar Model (FBM) has a pattern speed of Ω_b = 56 km s⁻¹ kpc⁻¹. 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 Σ₀, σ₀, R_s, and R_σ by a 25% and assessed that the change on the position of the overdensities is below 2 km s⁻¹. 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⁻¹) 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⁻¹) 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⁻¹ 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.

3.2. 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. 2023). 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” resolution³, and can be considered twice-differentiable. 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

$$\begin{array}{c}\hfill \mathrm{\Pi }\equiv \mathbb{I}-\frac{\mathrm{\nabla }p{\mathrm{\nabla }p}^{\top }}{{\mathrm{\nabla }p}^{\top }\mathrm{\nabla }p},\end{array}$$(9)

where 𝕀 is the identity. The maximum eigenvalue of ΠHΠ is the “gappiness” estimator of the point. It is important to notice that “gappiness” is the opposite of “overdensityness”, and 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⁻¹ 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 Sect. 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.

4. Data versus the 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.

4.1. Global ridges in the radial and azimuthal directions

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⁻¹ 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 visualisation 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 (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 2018b; Asano et al. 2020). In this execution, due to the use of a WT with a larger scale, two of these arches (A8 and 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 disc. 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 disc, which we label Hercules_In. This structure appears as a symmetric parabola (peak of |V_ϕ| at V_R = 0 km s⁻¹), in good agreement with Fig. 4 of Dehnen (2000), where x₁(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).

4.2. 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), coloured by their gradients ∂V_ϕ/∂R (top) and ∂V_ϕ/∂ϕ (bottom). Thus, we are including information of the large-scale shape, that is, 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, that is, the structures shown in the top panel of Fig. 1, which we included for visualisation 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. 2019a). 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 (Table 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 km s⁻¹ deg⁻¹ (Table 2), compatible with both models within the error bars. In ∂V_ϕ/∂R, we detected two sub-segments for Hercules1: in the V_R < 50 km s⁻¹ part of the arch, we measured ∂V_ϕ/∂R ∼ 29 km s⁻¹ kpc⁻¹, while in the V_R > 50 km s⁻¹ part we obtained ∂V_ϕ/∂R ∼ 25 km s⁻¹ kpc⁻¹. The overall average slope of Hercules1, without separating into segments, ends up being ∂V_ϕ/∂R = 28.1 ± 2.8 km s⁻¹ kpc⁻¹. The corresponding values of ∂V_ϕ/∂R for the SMB (36.5 km s⁻¹ kpc⁻¹) and FBM (40.5 km s⁻¹ kpc⁻¹), however, are both about 10 km s⁻¹ kpc⁻¹ larger than the gradient measured in the data. In Sect. 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₂ = 14 bar laps, in contrast to the t₂ = 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 Sect. 6, we study these overdensities in the models to explore an alternative hypothesis on the origin of these moving groups.

5. 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 (slow-fast bar degeneracy). Other parameters can also shift the position of the overdensities in the models, thus leading to extra degeneracies, for example 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 with the bar’s effects. Possible links between other moving groups and the overdensities seen in the models are discussed in Sect. 6.

5.1. 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:

• 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 2022), so we maintain the fiducial model value R_⊙ = 8 kpc.

• v₀: We keep this parameter fixed but we checked that changing it even by a considerable amount (±15 km s⁻¹) does not change significantly our results (variations in the radial slopes of about 1 km s⁻¹ kpc⁻¹).

• ϕ: 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 (Dehnen 2000). Therefore, we leave it fixed to the value of the fiducial models.

• Ω_p: We explored a parameter range of ±2 km s⁻¹ kpc⁻¹ 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.

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: VR of the overdensity. Third row: ∂Vϕ/∂ϕ of the overdensity. In black, we provide the statistical deviation to the measurements in Gaia DR3 (Table 2). Fourth row: Same as the third row but with ∂Vϕ/∂R.

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₀ or outside. An increase in β decreases the slope of the azimuthal frequency curve, thus sending the resonances inside R₀ inwards, and the resonances outside R₀ outwards. Conversely, if β decreases the frequencies change faster with radius, which pulls the resonances closer to R₀.

Regarding ∂V_ϕ/∂ϕ (third row in Fig. 3), it shows an almost linear decrease with Ω_b and β for both Hercules-like⁴ and the Hat-like overdensities across the parameter space for the slow bar (Cols. 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 (Table 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.

5.2. Comparison with Gaia data

The black contours in Fig. 3 show the statistical deviation between each galaxy model and the data, that is, the difference with regard to 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 Sect. 4.2. In some parts of the parameter 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 km s⁻¹ kpc⁻¹ and β ∼ 0.025. Opposite to that, the azimuthal gradient of the fast bar Hercules-like overdensity (second column) matches the data for Ω_b ∼ 54 km s⁻¹ kpc⁻¹ 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 km s⁻¹ kpc⁻¹, 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.

6. Azimuth and time exploration

In Sect. 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.

6.1. 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⁻¹. In Fig. 4, we present the azimuthal velocity (V_ϕ) of these structures across the azimuthal range of ±90° 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 links⁵. 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.

Fig. 4. Vϕ of the overdensities at VR = −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 t2 = 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 (dashed red lines). Middle: Structures for t2 = 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 t2 = 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.

Since this analysis is done at very short timescale (four 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₂ = 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⁻¹ and they are located in the opposite site of the major axis. Therefore, at V_R = 0 km s⁻¹ 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 stabilise 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.

6.2. 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⁻¹ 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 km s⁻¹ deg⁻¹. 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.

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

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.

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 (Fig. 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).

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 et al. 2020; Sanders et al. 2022). In the following section, we discuss the implications of this short timescale prediction.

7. 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 km s⁻¹ kpc⁻¹ is compatible with recent results such as the one in Ramos et al. (2018), of ∂V_ϕ/∂R = 26.5 km s⁻¹ kpc⁻¹, 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 Hercules⁶ 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 km s⁻¹ kpc⁻¹) 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 non-negligible 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.

Azimuthal gradient of Hercules and Horn. In this analysis, we measure an azimuthal slope for Hercules of ∂V_ϕ/∂ϕ = −0.63 ± 0.13 km s⁻¹ deg⁻¹. Using also Gaia DR3⁷, in Lucchini et al. (2024) they obtain a measurement of −0.74 ± 0.04 km s⁻¹ deg⁻¹. In this work, the V_R is marginalised, obtaining a smaller error than ours, but loosing the V_R information. Similarly, in B22 we computed ∂V_ϕ/∂ϕ = −0.5 km s⁻¹ deg⁻¹ using Gaia EDR3 data.

Using a combination of perturbation theory and the pendulum approximation (Monari et al. 2019a), in Monari et al. (2019b) they found a slope of ∂L_z/∂ϕ = −8 km s⁻¹kpc deg⁻¹ at a radius of 8.2 kpc, which corresponds to ∂V_ϕ/∂ϕ = −0.96 km s⁻¹ deg⁻¹, 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. (2019b), 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 cannot 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 (Table 2). Finally, we note that the measure in Lucchini et al. (2024), 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 observed 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⁻¹. For this same velocity separation, from our Fig. 6 we would obtain times of around 1 Gyr, that is to say 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 idealised 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 (Table 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.

8. 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 the following:

• We detected 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 observed a robust slope in azimuth for Hyades and Horn, which were previously unexplored.

• The radial gradient of Hercules is ∂V_ϕ/∂R = 28 ± 2.3 km s⁻¹ kpc⁻¹, and the azimuthal gradient is ∂V_ϕ/∂ϕ = −0.63 ± 0.13 km s⁻¹ deg⁻¹.

• 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, self-gravity, 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 an estimate of the bar age that is too low (∼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, especially 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 cannot capture the full complexity of our Galaxy, especially on a large scale. It emphasises the need of exploring more advanced models, and reminds us that an optimal model will have to go beyond matching local patterns.

Acknowledgments

We thank the anonymous referee for her/his helpful comments. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. We acknowledge the grants PID2021-125451NA-I00 and CNS2022-135232 funded by MICIU/AEI/10.13039/501100011033 and by “ERDF A way of making Europe”, by the “European Union” and by the “European Union Next Generation EU/PRTR”, and the Institute of Cosmos Sciences University of Barcelona (ICCUB, Unidad de Excelencia ‘María de Maeztu’) through grant CEX2019-000918-M. This work was partially supported by the OCRE awarded project Galactic Research in Cloud Services (Galactic RainCloudS). OCRE receives funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement no. 824079. M.B. acknowledges funding from the University of Barcelona’s official doctoral program for the development of a R+D+i project under the PREDOCS-UB grant. T.A. acknowledges the grant RYC2018-025968-I funded by MCIN/AEI/10.13039/501100011033 and by “ESF Investing in your future”. B.F. and G.M. acknowledge funding from the Agence Nationale de la Recherche (ANR project ANR-18-CE31-0006 and ANR-19-CE31- 0017) and from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 834148).

References

Appendix A: Time dependence of the arches

In order to understand how does the velocity distribution evolve in time, we integrate it 14 bar laps and compute the radial and azimuthal gradient of each overdensity. We show these results in Fig. A.1 (equivalent of Fig.2) and Table A.1 (equivalent of Tab. 2). We see that the both Hercules-like and Hat-like structures are located in the same position at all times. The gradients of the main structures are also very similar as the fiducial ones.

As for the presence of overdensities, the only features that are maintained are the Hercules-like and Hat-like for both pattern speeds, and the Horn-like for the fast bar. In the slow bar regime, a Horn-like structure can be produced by the 6:1 resonance, but our bar is a simple quadrupole that has no high order moments. Therefore, it is expected that these high order resonances are not appearing.

To gain an intuition in how do the overdensities evolve from one integration time to the other, we selected the structures in V_R = −60 km s⁻¹ and traced them along 50 bar laps (Fig. A.2). We observe that the structures are created in given positions of the phase space (green lines) and tend asymptotically to the boundaries of a resonance (red lines), following exponential decays in time. Understanding this behaviour could bring new intuitions to understanding the phase-mixing of this system.

Fig. A.1. Overdensities’ detection in the simulations, for t = 14 bar laps. They are coloured by their ∂Vϕ/∂R (top panels) and ∂Vϕ/∂ϕ (bottom panels) at each VR. 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. A.2. Time evolution of the structures at VR = −60 km s−1. The structures are created in given positions of the phase space (green lines) and tend asymptotically to the boundaries of a resonance (red lines), following exponential decays in time. Top: Slow bar model. Bottom: Fast bar model.

Appendix B: Correlation in other azimuths

In this section, we complement the parameter exploration done in Sect. 5 by exploring the same diagrams at different ϕ_b. We conclude that the trends observed at ϕ_b = −30 deg (Fig. 3) are maintained when looking at other azimuths. The only noticeable change is the steepness of the ∂V_ϕ/∂R and ∂V_ϕ/∂ϕ panels according to the distance to the major axis of the bar. In Fig. B.1 (ϕ_b = −20 deg) the contour lines are close to each other, indicating a rapid change of the properies. Opposite to that, in Fig. B.2 (ϕ_b = −40 deg), the contour lines are more separated, indicating a solution which is more stable.

Fig. B.1. Same as Fig. 3 but computed with ϕb = −20 deg.

Fig. B.2. Same as Fig. 3 but computed with ϕb = −40 deg.

Appendix C: Model versus model degeneracy

The computation of the Galactic models links a set of parameters (Ω_b, β, ϕ_b...) to a set of measurables (V_R, V_ϕ, ∂V_ϕ/∂R, ∂V_ϕ/∂ϕ). For a set of target measurables (e.g. Gaia DR3 measurables), a set of parameters will be optimal if the measurables it produces match the target measurables. Given a set of parameters, we define their statistical difference (σ, dimensionless) with respect to a target as the difference between the measurables in units of the measure error in the target. Notice that in this work we have not been able to obtain an optimal set of parameters to reproduce the target Gaia DR3 measurables. In this section we will use the fiducial models as targets.

Lets suppose that our target measurables are the ones from the fiducial SBM (Tab. 2). In this case, by definition the fiducial SBM will be an optimal set of parameters. To test the degeneracies, we want to check if there are other sets of parameters that also reproduce the SBM measurables. This is what we explore in Figure C.1, where we show the statistical differences between the target measurables (assuming those of the model SBM) and the models around the FBM using various combinations of measurables. The top panel exclusively utilises the measurables corresponding to the Hercules-like overdensity, while the bottom panel incorporates measurables for both Hercules and Hat. The black contours represent the statistical differences between the model and the target when considering only V_ϕ. When restricted to Hercules and its local value of V_ϕ alone (black contour, top panel), we observe a substantial region of parameter space where the solution is degenerate, with multiple combinations of Ω_b and β from the FBM models yielding the same V_ϕ. However, as we introduce additional measurables, such as ∂V_ϕ/∂ϕ (orange contour) and ∂V_ϕ/∂R (green contour), the statistical difference becomes more pronounced, scaling to over 4σ. When we further include a second structure like Hat (bottom panel), the degeneracy in the solution can be partially resolved using only V_ϕ. Once gradients are factored in, the distinction between models is once again emphasised.

To complement this study, we repeat the analysis assuming that the target measurables are the ones of the fiducial FBM (Table 2), and study the statistical difference with the models around it (Fig. C.2). The goal in this exercise is to characterise the local degeneracy. The results are very similar to the ones in Fig. C.1: we need to include the gradients to obtain a statistical difference.

This analysis suggests that the best combination of measurables to break the degeneracy is to consider both Hercules and Hat (bottom panel). Using only V_ϕ of both structures we already see that we are able to constrain a small region of Ω_b and β. In addition, including the azimuthal gradient brings the statistical significance over 8σ in almost the entire parameter space.

Fig. C.1. Differences between the fiducial SBM and the galaxy models around the FBM for Hercules (top) and Hercules + Hat (bottom). Black contours show the disagreement considering only Vϕ measurements. Orange contours are for the disagreement considering Vϕ and ∂Vϕ/∂ϕ. Green contours show the disagreement considering Vϕ, ∂Vϕ/∂ϕ, and ∂Vϕ/∂R.

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 theHercules- and Hat-like one (bottom).

All Tables

All Figures

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.

In the text

	Fig. 2. Moving group and overdensities’ detection, coloured by their ∂Vϕ/∂R (top panels) and ∂Vϕ/∂ϕ (bottom panels) at each VR. 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.
In the text

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: VR of the overdensity. Third row: ∂Vϕ/∂ϕ of the overdensity. In black, we provide the statistical deviation to the measurements in Gaia DR3 (Table 2). Fourth row: Same as the third row but with ∂Vϕ/∂R.

In the text

Fig. 4. Vϕ of the overdensities at VR = −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 t2 = 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 (dashed red lines). Middle: Structures for t2 = 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 t2 = 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.

In the text

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

	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 (Fig. 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).
In the text

	Fig. A.1. Overdensities’ detection in the simulations, for t = 14 bar laps. They are coloured by their ∂Vϕ/∂R (top panels) and ∂Vϕ/∂ϕ (bottom panels) at each VR. 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.
In the text

	Fig. A.2. Time evolution of the structures at VR = −60 km s−1. The structures are created in given positions of the phase space (green lines) and tend asymptotically to the boundaries of a resonance (red lines), following exponential decays in time. Top: Slow bar model. Bottom: Fast bar model.
In the text

	Fig. B.1. Same as Fig. 3 but computed with ϕb = −20 deg.
In the text

	Fig. B.2. Same as Fig. 3 but computed with ϕb = −40 deg.
In the text

	Fig. C.1. Differences between the fiducial SBM and the galaxy models around the FBM for Hercules (top) and Hercules + Hat (bottom). Black contours show the disagreement considering only Vϕ measurements. Orange contours are for the disagreement considering Vϕ and ∂Vϕ/∂ϕ. Green contours show the disagreement considering Vϕ, ∂Vϕ/∂ϕ, and ∂Vϕ/∂R.
In the text

	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 theHercules- and Hat-like one (bottom).
In the text