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

M. Bernet, P. Ramos, T. Antoja, G. Monari, B. Famaey

► To cite this version:

M. Bernet, P. Ramos, T. Antoja, G. Monari, B. Famaey. Radial and azimuthal gradients of the moving groups in Gaia DR3: The slow and fast bar degeneracy problem. Astronomy and Astrophysics - A&A, 2024, 686, 10.1051/0004-6361/202348410. insu-04608580

HAL Id: insu-04608580

https://insu.hal.science/insu-04608580

Submitted on 11 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Radial and azimuthal gradients of the moving groups in Gaia DR3: The slow and fast bar degeneracy problem

M. Bernet, P. Ramos, T. Antoja, G. Monari, B. Famaey

► To cite this version:

M. Bernet, P. Ramos, T. Antoja, G. Monari, B. Famaey. Radial and azimuthal gradients of the moving groups in Gaia DR3: The slow and fast bar degeneracy problem. Astronomy & Astrophysics - A&A, 2024, 686, 10.1051/0004-6361/202348410. insu-04608580

HAL Id: insu-04608580

https://insu.hal.science/insu-04608580v1

Submitted on 11 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Radial and azimuthal gradients of the moving groups in Gaia DR3: The slow and fast bar degeneracy problem

M. Bernet<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> Ⓓ, P. Ramos<fnref n="4" /> Ⓓ, T. Antoja<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> Ⓓ, G. Monari<fnref n="5" />, and B. Famaey<fnref n="5" /> Ⓓ

Received 27 October 2023 / Accepted 1 March 2024

## ABSTRACT

**Context.** The structure and dynamics of the central bar of the Milky Way (MW) 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, allowing for a better understanding of the complex imprints of the bar on the phase space.

**Aims.** We aim to identify and characterise the dynamical moving groups across the MW disc, and use their large-scale distribution to help constrain the properties of the Galactic bar.

**Methods.** We used 1D 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 $V_\phi$ of the moving groups. We simulated thousands of perturbed distribution functions using backward integration, sweeping a large portion of parameter space of feasible Galaxy models that include a bar, in order to compare them with the data and to explore and quantify the degeneracies.

**Results.** The radial gradient of the Hercules moving group ($\partial V_\phi / \partial R = 28.1 \pm 2.8 \text{ km s}^{-1} \text{ kpc}^{-1}$) cannot be reproduced by our simple models of the Galaxy that show much larger slopes both for a fast and a slow bar. This suggests the need for more complex dynamics (e.g. a different bar potential, spiral arms, a slowing bar, a complex circular velocity curve, external perturbations, etc.). We measured an azimuthal gradient for Hercules of $\partial V_\phi / \partial \phi = -0.63 \pm 0.13 \text{ km s}^{-1} \text{ deg}^{-1}$ and find that it is compatible with both the slow and fast bar models. Our analysis points out that in using this type of analysis, at least two moving groups are needed to start breaking the degeneracies.

**Conclusions.** 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.

**Key words.** methods: data analysis – Galaxy: disk – Galaxy: evolution – Galaxy: kinematics and dynamics – Galaxy: structure

## 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, \phi, Z, V_R, V_\phi$, 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_\phi$ (Antoja et al. 2018; Kawata et al. 2018; Fragkoudi et al. 2019), a wave in $L_Z - V_R$ ($L_Z = RV_\phi$, 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_\phi$ (Gaia Collaboration 2018b; Ramos et al. 2018).

Out of all of these projections, the $V_R - V_\phi$ 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 $\Omega_b \sim 55 \text{ km s}^{-1} \text{ kpc}^{-1}$ (the short and fast bar scenario,

<footnote n="1">¹ Departament de Física Quàntica i Astrofísica (FQA), Universitat de Barcelona (UB), C Martí i Franquès, 1, 08028 Barcelona, Spain
e-mail: mbernet@fqa.ub.edu</footnote>

<footnote n="2">² Institut de Ciències del Cosmos (ICCUB), Universitat de Barcelona (UB), C Martí i Franquès, 1, 08028 Barcelona, Spain</footnote>

<footnote n="3">³ Institut d'Estudis Espacials de Catalunya (IEEC), Edifici RDIT, Campus UPC, 08860 Castelldefels (Barcelona), Spain</footnote>

<footnote n="4">⁴ National Astronomical Observatory of Japan, Mitaka-shi, Tokyo 181-8588, Japan</footnote>

<footnote n="5">⁵ Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, 11 rue de l'Université, 67000 Strasbourg, France</footnote>
# Radial and azimuthal gradients of the moving groups in Gaia DR3: The slow and fast bar degeneracy problem

M. Bernet<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> Ⓓ, P. Ramos<fnref n="4" /> Ⓓ, T. Antoja<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> Ⓓ, G. Monari<fnref n="5" />, and B. Famaey<fnref n="5" /> Ⓓ

Received 27 October 2023 / Accepted 1 March 2024

## ABSTRACT

**Context.** The structure and dynamics of the central bar of the Milky Way (MW) 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, allowing for a better understanding of the complex imprints of the bar on the phase space.

**Aims.** We aim to identify and characterise the dynamical moving groups across the MW disc, and use their large-scale distribution to help constrain the properties of the Galactic bar.

**Methods.** We used 1D 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 $V_\phi$ of the moving groups. We simulated thousands of perturbed distribution functions using backward integration, sweeping a large portion of parameter space of feasible Galaxy models that include a bar, in order to compare them with the data and to explore and quantify the degeneracies.

**Results.** The radial gradient of the Hercules moving group ($\partial V_\phi / \partial R = 28.1 \pm 2.8 \text{ km s}^{-1} \text{ kpc}^{-1}$) cannot be reproduced by our simple models of the Galaxy that show much larger slopes both for a fast and a slow bar. This suggests the need for more complex dynamics (e.g. a different bar potential, spiral arms, a slowing bar, a complex circular velocity curve, external perturbations, etc.). We measured an azimuthal gradient for Hercules of $\partial V_\phi / \partial \phi = -0.63 \pm 0.13 \text{ km s}^{-1} \text{ deg}^{-1}$ and find that it is compatible with both the slow and fast bar models. Our analysis points out that in using this type of analysis, at least two moving groups are needed to start breaking the degeneracies.

**Conclusions.** 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.

**Key words.** methods: data analysis – Galaxy: disk – Galaxy: evolution – Galaxy: kinematics and dynamics – Galaxy: structure

## 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, $\phi$, Z, $V_R$, $V_\phi$, 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_\phi$ (Antoja et al. 2018; Kawata et al. 2018; Fragkoudi et al. 2019), a wave in $L_Z - V_R$ ($L_Z = RV_\phi$, 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_\phi$ (Gaia Collaboration 2018b; Ramos et al. 2018).

Out of all of these projections, the $V_R - V_\phi$ 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 $\Omega_b \sim 55 \text{ km s}^{-1} \text{ kpc}^{-1}$ (the short and fast bar scenario,

<footnote n="1">¹ Departament de Física Quàntica i Astrofísica (FQA), Universitat de Barcelona (UB), C Martí i Franquès, 1, 08028 Barcelona, Spain
e-mail: mbernet@fqa.ub.edu</footnote>

<footnote n="2">² Institut de Ciències del Cosmos (ICCUB), Universitat de Barcelona (UB), C Martí i Franquès, 1, 08028 Barcelona, Spain</footnote>

<footnote n="3">³ Institut d'Estudis Espacials de Catalunya (IEEC), Edifici RDIT, Campus UPC, 08860 Castelldefels (Barcelona), Spain</footnote>

<footnote n="4">⁴ National Astronomical Observatory of Japan, Mitaka-shi, Tokyo 181-8588, Japan</footnote>

<footnote n="5">⁵ Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, 11 rue de l'Université, 67000 Strasbourg, France</footnote>
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 $\Omega_b \sim 40 \text{ km s}^{-1} \text{ kpc}^{-1}$ (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_\phi$, a novel set of measurables: the gradients of the moving groups at the SN in the radial ($\partial V_\phi/\partial R$) and azimuthal ($\partial V_\phi/\partial \phi$) 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 $\Omega_b$ and the slope of the rotation curve $\beta$ 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 ≤ 14500 K and 6 ≤ grvs_mag ≤ 12. However, after the correction, a residual bias of a few km s⁻<fnref n="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.

We transformed the observables into cylindrical phase space coordinates using $R_0 = 8.277 \text{ kpc}$ (GRAVITY Collaboration 2022)<fnref n="1" />, $Z_\odot = 0.0208 \text{ kpc}$ (Bennett & Bovy 2019), and $U_\odot = 9.3$, $V_\odot + V_c(R_0) = 251.5$, and $W_\odot = 8.59 \text{ km s}^{-1}$ from the proper motion of SagA* (Reid & Brunthaler 2020) and its radial velocity (GRAVITY Collaboration 2022). The reference system is right-handed with $\phi$ oriented contrary to the disc rotation and thus $V_\phi < 0$ for most stars. We set the origin at the Sun’s azimuth ($\phi_\odot = 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 ($\Delta R$, $\Delta\phi$, $\Delta Z$) and run a Wavelet Transform (WT) on the velocity distribution ($V_R$, $V_\phi$) 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_\phi$ (Ramos et al. 2018). Because of this geometry, in the second step, we slice each ($V_R$, $V_\phi$) diagram in $\Delta V_R$ bins, and run a 1D WT in each $V_\phi$ histogram. The peaks obtained from these 1D WT are connected through the configuration space to form global substructures

<footnote n="1">¹ Clarke & Gerhard (2022) used an alternative measure of $R_0$, the centre of the bulge. In their study they found a great agreement between this measurement and the distance to Sag A*, confirming the hypothesis that Sag A* is at rest at the centre of the bulge.</footnote>
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 $\Omega_b \sim 40 \text{ km s}^{-1} \text{ kpc}^{-1}$ (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_\phi$, a novel set of measurables: the gradients of the moving groups at the SN in the radial ($\partial V_\phi/\partial R$) and azimuthal ($\partial V_\phi/\partial \phi$) 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 $\Omega_b$ and the slope of the rotation curve $\beta$ 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 ≤ 14500 K and 6 ≤ grvs_mag ≤ 12. However, after the correction, a residual bias of a few km s⁻<fnref n="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.

We transformed the observables into cylindrical phase space coordinates using $R_0 = 8.277 \text{ kpc}$ (GRAVITY Collaboration 2022)<fnref n="1" />, $Z_\odot = 0.0208 \text{ kpc}$ (Bennett & Bovy 2019), and $U_\odot = 9.3$, $V_\odot + V_c(R_0) = 251.5$, and $W_\odot = 8.59 \text{ km s}^{-1}$ from the proper motion of SagA* (Reid & Brunthaler 2020) and its radial velocity (GRAVITY Collaboration 2022). The reference system is right-handed with $\phi$ oriented contrary to the disc rotation and thus $V_\phi < 0$ for most stars. We set the origin at the Sun’s azimuth ($\phi_\odot = 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 ($\Delta R$, $\Delta\phi$, $\Delta Z$) and run a Wavelet Transform (WT) on the velocity distribution ($V_R$, $V_\phi$) 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_\phi$ (Ramos et al. 2018). Because of this geometry, in the second step, we slice each ($V_R$, $V_\phi$) diagram in $\Delta V_R$ bins, and run a 1D WT in each $V_\phi$ histogram. The peaks obtained from these 1D WT are connected through the configuration space to form global substructures

<footnote n="1">¹ Clarke & Gerhard (2022) used an alternative measure of $R_0$, the centre of the bulge. In their study they found a great agreement between this measurement and the distance to Sag A*, confirming the hypothesis that Sag A* is at rest at the centre of the bulge.</footnote>
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.

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φ/∂φ.

## 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

$$f(E, L) \propto \frac{\Sigma(R_E)}{\sigma_R^2(R_E)} \exp\left[\frac{\Omega(R_E)(L - L_c(E))}{\sigma_R^2(R_E)}\right], \quad (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

$$\Sigma(R) = \Sigma_0 \exp(-R/R_s), \quad (2)$$

$$\sigma_R(R) = \sigma_0 \exp(-R/R_\sigma), \quad (3)$$
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<sub>φ</sub>/∂R and ∂V<sub>φ</sub>/∂φ).

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<sub>R</sub>): [-120, 120] km s<sup>-1</sup> in steps of 2 km s<sup>-1</sup>, ΔV<sub>R</sub> = ±15 km s<sup>-1</sup> 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<sub>R</sub>, 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<sub>R</sub> 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 two dimensions all at once is challenging. To reduce the dimensionality of the problem, we computed the radial (∂V<sub>φ</sub>/∂R) and azimuthal (∂V<sub>φ</sub>/∂φ) 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<sub>φ</sub> for each measurement (WT size, σ = 6.25 km s<sup>-1</sup>) 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<sub>φ</sub>/∂R, and below 0.05 km s⁻¹ deg⁻¹ in ∂V<sub>φ</sub>/∂φ.

## 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<sub>R</sub>, V<sub>φ</sub>) 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

$$f(E, L) \propto \frac{\Sigma(R_E)}{\sigma_R^2(R_E)} \exp\left[\frac{\Omega(R_E)(L - L_c(E))}{\sigma_R^2(R_E)}\right], \quad (1)$$

where R<sub>E</sub>, Ω(R<sub>E</sub>), and L<sub>c</sub>(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 σ<sub>R</sub>(R) are defined as

$$\Sigma(R) = \Sigma_0 \exp(-R/R_s), \quad (2)$$

$$\sigma_R(R) = \sigma_0 \exp(-R/R_\sigma), \quad (3)$$
where the values of the constants $\Sigma_0$, $\sigma_0$, $R_s$, and $R_\sigma$ (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

$$v_c(R) = v_0(R/R_0)^\beta, \quad (4)$$

where $v_0$ is the circular velocity at the Solar radius $R_0$, and $\beta$ characterises the shape of the rotation curve.

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

$$\Phi_b(R, \phi) = A_b(t) \cos(2(\phi - \phi_b - \Omega_b t)) \times \begin{cases} -(R_b/R)^3, & \text{for } R \ge R_b, \\ (R/R_b)^3 - 2, & \text{for } R \le R_b, \end{cases} \quad (5)$$

where $\Omega_b$ is the pattern speed of the bar, $\phi_b$ is angle of its long axis and $R_b$ its size, which we fix to be 80% of the corotation radius<fnref n="2" />, given by

$$R_{CR} = R_0(\Omega_b/\Omega_0)^{1/(\beta-1)}, \quad (6)$$

where $\Omega_0 \equiv v_0/R_0$ 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_1$ as

$$A_b(t) = A_f \left[ \frac{3}{16} \xi^5 - \frac{5}{8} \xi^3 + \frac{15}{16} \xi + \frac{1}{2} \right], \quad \xi = 2 \frac{t}{t_1} - 1, \quad (7)$$

and it 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 \equiv 2\pi\Omega_b$. Finally, the amplitude constant is defined as

$$A_f = \alpha v_0^{2/3} \left( \frac{R}{R_b} \right)^3, \quad (8)$$

which depends on the dimensionless bar strength $\alpha$, 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 $\Omega_b = 39 \text{ km s}^{-1} \text{ kpc}^{-1}$. The Fast Bar Model (FBM) has a pattern speed of $\Omega_b = 56 \text{ km s}^{-1} \text{ kpc}^{-1}$. The Sun is at $R_\odot = 8 \text{ kpc}$, and the bar is at $\phi_b = -30 \text{ 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 $\Sigma_0$, $\sigma_0$, $R_s$, and $R_\sigma$ by a 25% and assessed that the change on the position of the overdensities is below $2 \text{ km 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

<table><caption>**Table 1.** Parameters of the fiducial models.</caption><thead><tr><th>Symbol</th><th>Units</th><th>SBM</th><th>FBM</th></tr></thead><tbody><tr><td>Ω<sub>b</sub></td><td>km s<sup>-1</sup> kpc<sup>-1</sup></td><td>39</td><td>56</td></tr><tr><td>t<sub>1</sub></td><td>Gyr</td><td>2T<sub>b</sub> = 0.32</td><td>2T<sub>b</sub> = 0.22</td></tr><tr><td>t<sub>2</sub></td><td>Gyr</td><td>4T<sub>b</sub> = 0.64</td><td>4T<sub>b</sub> = 0.44</td></tr><tr><td>R<sub>0</sub></td><td>kpc</td><td>8</td><td></td></tr><tr><td>v<sub>0</sub></td><td>km s<sup>-1</sup></td><td>220</td><td></td></tr><tr><td>α</td><td></td><td>0.01</td><td></td></tr><tr><td>β</td><td></td><td>0</td><td></td></tr><tr><td>Σ<sub>0</sub></td><td>km s<sup>-1</sup></td><td>1</td><td></td></tr><tr><td>σ<sub>0</sub></td><td>km s<sup>-1</sup></td><td>48</td><td></td></tr><tr><td>R<sub>s</sub></td><td>kpc</td><td>2.64</td><td></td></tr><tr><td>R<sub>σ</sub></td><td>kpc</td><td>8</td><td></td></tr></tbody></table>

measurements are robust to significant variations of the input parameters. We remark that the local velocity dispersion of the used model ($\sigma_R = 17.6 \text{ 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 ($\sigma_R = 35 \text{ 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 \text{ 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.

## 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<fnref n="3" />, 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.

<footnote n="2">² This value is the same used by Dehnen (2000) and is close to the reference value of $R = R_{CR}/R_{bar} = 1.2$ (Athanassoula 1992, but see Font et al. 2017).</footnote>

<footnote n="3">³ In practice, to have reasonable computing costs, we limit our velocity grid to a resolution of $0.2 \text{ km s}^{-1}$, which is of the order of the velocity uncertainties at the SN.</footnote>
where the values of the constants $\Sigma_0$, $\sigma_0$, $R_s$, and $R_\sigma$ (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

$$v_c(R) = v_0(R/R_0)^\beta, \quad (4)$$

where $v_0$ is the circular velocity at the Solar radius $R_0$, and $\beta$ characterises the shape of the rotation curve.

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

$$\Phi_b(R, \phi) = A_b(t) \cos(2(\phi - \phi_b - \Omega_b t)) \times \begin{cases} -(R_b/R)^3, & \text{for } R \ge R_b, \\ (R/R_b)^3 - 2, & \text{for } R \le R_b, \end{cases} \quad (5)$$

where $\Omega_b$ is the pattern speed of the bar, $\phi_b$ is angle of its long axis and $R_b$ its size, which we fix to be 80% of the corotation radius<fnref n="2" />, given by

$$R_{CR} = R_0(\Omega_b/\Omega_0)^{1/(\beta-1)}, \quad (6)$$

where $\Omega_0 \equiv v_0/R_0$ 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_1$ as

$$A_b(t) = A_f \left[ \frac{3}{16} \xi^5 - \frac{5}{8} \xi^3 + \frac{15}{16} \xi + \frac{1}{2} \right], \quad \xi = 2 \frac{t}{t_1} - 1, \quad (7)$$

and it 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 \equiv 2\pi\Omega_b$. Finally, the amplitude constant is defined as

$$A_f = \alpha v_0^{2/3} \left( \frac{R}{R_b} \right)^3, \quad (8)$$

which depends on the dimensionless bar strength $\alpha$, 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 $\Omega_b = 39 \text{ km s}^{-1} \text{ kpc}^{-1}$. The Fast Bar Model (FBM) has a pattern speed of $\Omega_b = 56 \text{ km s}^{-1} \text{ kpc}^{-1}$. The Sun is at $R_\odot = 8 \text{ kpc}$, and the bar is at $\phi_b = -30 \text{ 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 $\Sigma_0$, $\sigma_0$, $R_s$, and $R_\sigma$ by a 25% and assessed that the change on the position of the overdensities is below $2 \text{ km 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

<table><caption>**Table 1.** Parameters of the fiducial models.</caption><thead><tr><th>Symbol</th><th>Units</th><th>SBM</th><th>FBM</th></tr></thead><tbody><tr><td>Ω<sub>b</sub></td><td>km s<sup>-1</sup> kpc<sup>-1</sup></td><td>39</td><td>56</td></tr><tr><td>t<sub>1</sub></td><td>Gyr</td><td>2T<sub>b</sub> = 0.32</td><td>2T<sub>b</sub> = 0.22</td></tr><tr><td>t<sub>2</sub></td><td>Gyr</td><td>4T<sub>b</sub> = 0.64</td><td>4T<sub>b</sub> = 0.44</td></tr><tr><td>R<sub>0</sub></td><td>kpc</td><td>8</td><td></td></tr><tr><td>v<sub>0</sub></td><td>km s<sup>-1</sup></td><td>220</td><td></td></tr><tr><td>α</td><td></td><td>0.01</td><td></td></tr><tr><td>β</td><td></td><td>0</td><td></td></tr><tr><td>Σ<sub>0</sub></td><td>km s<sup>-1</sup></td><td>1</td><td></td></tr><tr><td>σ<sub>0</sub></td><td>km s<sup>-1</sup></td><td>48</td><td></td></tr><tr><td>R<sub>s</sub></td><td>kpc</td><td>2.64</td><td></td></tr><tr><td>R<sub>σ</sub></td><td>kpc</td><td>8</td><td></td></tr></tbody></table>

measurements are robust to significant variations of the input parameters. We remark that the local velocity dispersion of the used model ($\sigma_R = 17.6 \text{ 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 ($\sigma_R = 35 \text{ 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 \text{ 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.

## 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<fnref n="3" />, 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.

<footnote n="2">² This value is the same used by Dehnen (2000) and is close to the reference value of $R = R_{CR}/R_{bar} = 1.2$ (Athanassoula 1992, but see Font et al. 2017).</footnote>

<footnote n="3">³ In practice, to have reasonable computing costs, we limit our velocity grid to a resolution of $0.2 \text{ km s}^{-1}$, which is of the order of the velocity uncertainties at the SN.</footnote>
For each point in the $V_R - V_\phi$ grid, we computed the gradi- ent vector $\nabla p$ and the Hessian matrix $H$ of the density on the velocity distribution. Then we calculated the projected second derivative tensor $\Pi H \Pi$, with

$$ \Pi \equiv \mathbb{I} - \frac{\nabla p \nabla p^\top}{\nabla p^\top \nabla p}, \qquad (9) $$

where $\mathbb{I}$ is the identity. The maximum eigenvalue of $\Pi H \Pi$ 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 $\Pi H \Pi$ gives information on the posi- tion 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 mod- els, and found it to be below $1 \text{ km s}^{-1}$ in the main overdensities, well below the WT size used in the data. When the overden- sities 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 signifi- cant 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_\phi$. This reveals a complex velocity field, rich in overden- sities, 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_\phi$ 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_\phi$ for each overdentisty of interest and selected the maximum as their representative. We note, however, that sometimes one over- density 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 over- densities after the split. Overall, this representative selection is feasible thanks to the stability of the overdensities in the sim- ulations, 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 gradi- ents 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_\phi$ in the radial ($\phi = 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 else- where, in this region it clearly struggles to obtain reliable detec- tions. Thus, we ignore the detections with $V_\phi > 80R - 740 \text{ 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_\phi$ of the moving groups as a function of $R$ and $\phi$, selecting, from all the struc- tures in each moving group at different $V_R$, only the one that cov- ers 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_\phi|$ part of the inner disc, and the azimuthal extension of Hyades and Sirius up to $\pm 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 hypoth- esised that, due to the likely weaker effects of the bar resonances at these radii, they could be due to spiral arms and/or the interac- tion 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 azimuthal 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 $\pm 40$ deg, while the Horn shows a negative slope for $\phi < 0$ and a flat profile for $\phi > 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 detec- tion of the methodology. Interestingly, in this work, a new set of structures have appeared in the inner disc, which we label Her- cules_In. This structure appears as a symmetric parabola (peak of $|V_\phi|$ at $V_R = 0 \text{ 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_\phi$ and $V_R - V_\phi$ 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 litera- ture 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
For each point in the $V_R - V_\phi$ grid, we computed the gradient vector $\nabla p$ and the Hessian matrix $H$ of the density on the velocity distribution. Then we calculated the projected second derivative tensor $\Pi H \Pi$, with

$$ \Pi \equiv \mathbb{I} - \frac{\nabla p \nabla p^\top}{\nabla p^\top \nabla p}, \qquad (9) $$

where $\mathbb{I}$ is the identity. The maximum eigenvalue of $\Pi H \Pi$ 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 $\Pi H \Pi$ 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 \text{ 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 Sect. 2.2, but using the “overdensityness” in $V_R - V_\phi$. 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_\phi$ 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_\phi$ 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_\phi$ in the radial ($\phi = 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_\phi > 80R - 740 \text{ 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_\phi$ of the moving groups as a function of $R$ and $\phi$, 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_\phi|$ part of the inner disc, and the azimuthal extension of Hyades and Sirius up to $\pm 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 azimuthal 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 $\pm 40$ deg, while the Horn shows a negative slope for $\phi < 0$ and a flat profile for $\phi > 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_\phi|$ at $V_R = 0 \text{ 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_\phi$ and $V_R - V_\phi$ 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
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Azimuthal velocity of the kinematic substructures as a function of radius ($\phi = 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, $\sigma = 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.</figcaption></figure>

models (middle and right panels), coloured by their gradients $\partial V_\phi / \partial R$ (top) and $\partial V_\phi / \partial \phi$ (bottom). Thus, we are including information of the large-scale shape, that is, how the velocities of the structures change with $R$ and $\phi$. 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 $\partial V_\phi / \partial \phi$ ~

$-0.63$ km s$^{-1}$ deg$^{-1}$ (Table 2), compatible with both models within the error bars. In $\partial V_\phi / \partial R$, we detected two sub-segments for Hercules1: in the $V_R < 50$ km s$^{-1}$ part of the arch, we measured $\partial V_\phi / \partial R \sim 29$ km s$^{-1}$ kpc$^{-1}$, while in the $V_R > 50$ km s$^{-1}$ part we obtained $\partial V_\phi / \partial R \sim 25$ km s$^{-1}$ kpc$^{-1}$. The overall average slope of Hercules1, without separating into segments, ends up being $\partial V_\phi / \partial R = 28.1 \pm 2.8$ km s$^{-1}$ kpc$^{-1}$. The corresponding values of $\partial V_\phi / \partial R$ for the SMB (36.5 km s$^{-1}$ kpc$^{-1}$) and FBM (40.5 km s$^{-1}$ kpc$^{-1}$), however, are both about 10 km s$^{-1}$ kpc$^{-1}$ 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 $\partial V_\phi / \partial 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
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Azimuthal velocity of the kinematic substructures as a function of radius ($\phi = 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, $\sigma = 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.</figcaption></figure>

models (middle and right panels), coloured by their gradients $\partial V_\phi / \partial R$ (top) and $\partial V_\phi / \partial \phi$ (bottom). Thus, we are including information of the large-scale shape, that is, how the velocities of the structures change with $R$ and $\phi$. 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 $\partial V_\phi / \partial \phi$ ~

$-0.63$ km s$^{-1}$ deg$^{-1}$ (Table 2), compatible with both models within the error bars. In $\partial V_\phi / \partial R$, we detected two sub-segments for Hercules1: in the $V_R < 50$ km s$^{-1}$ part of the arch, we measured $\partial V_\phi / \partial R \sim 29$ km s$^{-1}$ kpc$^{-1}$, while in the $V_R > 50$ km s$^{-1}$ part we obtained $\partial V_\phi / \partial R \sim 25$ km s$^{-1}$ kpc$^{-1}$. The overall average slope of Hercules1, without separating into segments, ends up being $\partial V_\phi / \partial R = 28.1 \pm 2.8$ km s$^{-1}$ kpc$^{-1}$. The corresponding values of $\partial V_\phi / \partial R$ for the SMB (36.5 km s$^{-1}$ kpc$^{-1}$) and FBM (40.5 km s$^{-1}$ kpc$^{-1}$), however, are both about 10 km s$^{-1}$ kpc$^{-1}$ 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 $\partial V_\phi / \partial 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
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Moving group and overdensities' detection, coloured by their $\partial V_\phi / \partial R$ (top panels) and $\partial V_\phi / \partial \phi$ (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 overdensities.</figcaption></figure>

<table><caption>**Table 2.** Mean slope and the dispersion within the box (see Fig. 2) in the data and both fiducial models.</caption><thead><tr><th></th><th></th><th>∂V<sub>φ</sub>/∂R (σ)<br>[km s<sup>-1</sup> kpc<sup>-1</sup>]</th><th>∂V<sub>φ</sub>/∂φ (σ)<br>[km s<sup>-1</sup> deg<sup>-1</sup>]</th></tr></thead><tbody><tr><td rowspan="3">Hercules</td><td>Data</td><td>28.1 (2.8)</td><td>-0.63 (0.13)</td></tr><tr><td>39</td><td>36.5 (1.1)</td><td>-0.49 (0.05)</td></tr><tr><td>56</td><td>40.5 (0.9)</td><td>-0.76 (0.01)</td></tr><tr><td rowspan="3">Arch/Hat</td><td>Data</td><td>28.7 (6.8)</td><td>0.07 (0.13)</td></tr><tr><td>39</td><td>25.7 (1.4)</td><td>-0.20 (0.02)</td></tr><tr><td>56</td><td>22.7 (0.1)</td><td>-0.20 (0.01)</td></tr></tbody></table>

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 over-

densities 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_\phi$ 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.
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Moving group and overdensities' detection, coloured by their $\partial V_\phi / \partial R$ (top panels) and $\partial V_\phi / \partial \phi$ (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 overdensities.</figcaption></figure>

<table><caption>**Table 2.** Mean slope and the dispersion within the box (see Fig. 2) in the data and both fiducial models.</caption><thead><tr><th></th><th></th><th>∂V<sub>φ</sub>/∂R (σ)<br>[km s<sup>-1</sup> kpc<sup>-1</sup>]</th><th>∂V<sub>φ</sub>/∂φ (σ)<br>[km s<sup>-1</sup> deg<sup>-1</sup>]</th></tr></thead><tbody><tr><td rowspan="3">Hercules</td><td>Data</td><td>28.1 (2.8)</td><td>-0.63 (0.13)</td></tr><tr><td>39</td><td>36.5 (1.1)</td><td>-0.49 (0.05)</td></tr><tr><td>56</td><td>40.5 (0.9)</td><td>-0.76 (0.01)</td></tr><tr><td rowspan="3">Arch/Hat</td><td>Data</td><td>28.7 (6.8)</td><td>0.07 (0.13)</td></tr><tr><td>39</td><td>25.7 (1.4)</td><td>-0.20 (0.02)</td></tr><tr><td>56</td><td>22.7 (0.1)</td><td>-0.20 (0.01)</td></tr></tbody></table>

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 over-

densities 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_\phi$ 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 $\partial V_{\phi}/\partial R$.

- *R*: It can be fitted by external measurements (GRAVITY Collaboration 2022), so we maintain the fiducial model value $R_{\odot} = 8$ kpc.

- $v_0$: 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⁻¹).

- $\phi$: In each model, we sweep $\phi$ between -10 and 10 deg, in steps of 1 deg to compute $\partial V_{\phi}/\partial \phi$.

- $\phi_b$: It can also be fitted by external measurements, though there is less consensus. In this section, we keep the value of $\phi_b = -30$ deg used in the fiducial model. In Appendix B we repeat the analysis for $\phi_b = -20, -40$ deg and see that, to first order, all the trends are the same.

- $\alpha$ (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.

- $\Omega_p$: We explored a parameter range of ±2 km s⁻¹ kpc⁻¹ around the slow and fast bar pattern speeds, within which the $V_{\phi}$ of the overdensities in the models remain compatible with Hercules and Hat.

- $\beta$ (rotation curve slope): We explored this parameter around the flat rotation curve ($0 \pm 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_{\phi}$, $V_R$, $\partial V_{\phi}/\partial \phi$, $\partial V_{\phi}/\partial R$ (different rows) of the Hercules-like (left) and Hat-like (right) overdensities behave, as a function of $\Omega_p$ (horizontal axis in all panels) and $\beta$ (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_{\phi}$, and vice versa. The first row of Fig. 3 shows that the azimuthal velocity $|V_{\phi}|$ of both the Hercules- and Hat-like groups decreases as $\Omega_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 $\beta$: while the $|V_{\phi}|$ of Hercules-like overdensities decreases with $\beta$, 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 $\beta$ affect the position of the resonances differently whether they are inside $R_0$ or outside. An increase in $\beta$ 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 $\beta$ decreases the frequencies change faster with radius, which pulls the resonances closer to $R_0$.

Regarding $\partial V_{\phi}/\partial \phi$ (third row in Fig. 3), it shows an almost linear decrease with $\Omega_b$ and $\beta$ for both Hercules-like<fnref n="4" /> and the Hat-like overdensities across the parameter space for the slow bar (Colrs. 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 $\partial V_{\phi}/\partial R$ of Hercules-like structures increases almost linearly with both $\Omega_b$ and $\beta$, while in Hat-like structures, $\partial V_{\phi}/\partial R$ increases with $\Omega_b$ and decreases with $\beta$. At first glance, this difference in behaviour might be related to the effect $\beta$ has on $V_{\phi}$, 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 $\beta$ (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 $\Omega_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\sigma$ 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\sigma$ 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$\sigma$ curves). In the slow bar Hercules-like overdensity (leftmost column) the azimuthal gradient matches the data for $\Omega_b \sim 40$ km s⁻¹ kpc⁻¹ and $\beta \sim 0.025$. Opposite to that, the azimuthal gradient of the fast bar Hercules-like overdensity (second column) matches the data for $\Omega_b \sim 54$ km s⁻¹ kpc⁻¹ and $\beta \sim -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 $\Omega_b$ and $\beta$ that reproduces all the measures in the data. To obtain a match in radial gradient we would need to decrease the parameters under $\Omega_b \ll 37$ km s⁻¹ kpc⁻¹, away from the match in $V_{\phi}$, and $\beta \ll -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\sigma$ when comparing models and data. This confirms that our models are sub-optimal

<footnote n="4">⁴ The Hercules-Slow Bar (leftmost column) presents an irregularity around $\beta = 0.02$ in the panels of $V_R$, $\partial V_{\phi}/\partial \phi$, and $\partial V_{\phi}/\partial R$ due to the split in the arch (see Sect. 3.2). The global trends are maintained, so this does not affect our final conclusions.</footnote>
## 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 $\partial V_{\phi}/\partial R$.

- *R*: It can be fitted by external measurements (GRAVITY Collaboration 2022), so we maintain the fiducial model value $R_{\odot} = 8$ kpc.

- $v_0$: 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⁻¹).

- $\phi$: In each model, we sweep $\phi$ between -10 and 10 deg, in steps of 1 deg to compute $\partial V_{\phi}/\partial \phi$.

- $\phi_b$: It can also be fitted by external measurements, though there is less consensus. In this section, we keep the value of $\phi_b = -30$ deg used in the fiducial model. In Appendix B we repeat the analysis for $\phi_b = -20, -40$ deg and see that, to first order, all the trends are the same.

- $\alpha$ (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.

- $\Omega_p$: We explored a parameter range of ±2 km s⁻¹ kpc⁻¹ around the slow and fast bar pattern speeds, within which the $V_{\phi}$ of the overdensities in the models remain compatible with Hercules and Hat.

- $\beta$ (rotation curve slope): We explored this parameter around the flat rotation curve ($0 \pm 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_{\phi}$, $V_R$, $\partial V_{\phi}/\partial \phi$, $\partial V_{\phi}/\partial R$ (different rows) of the Hercules-like (left) and Hat-like (right) overdensities behave, as a function of $\Omega_p$ (horizontal axis in all panels) and $\beta$ (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_{\phi}$, and vice versa. The first row of Fig. 3 shows that the azimuthal velocity $|V_{\phi}|$ of both the Hercules- and Hat-like groups decreases as $\Omega_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 $\beta$: while the $|V_{\phi}|$ of Hercules-like overdensities decreases with $\beta$, 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 $\beta$ affect the position of the resonances differently whether they are inside $R_0$ or outside. An increase in $\beta$ 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 $\beta$ decreases the frequencies change faster with radius, which pulls the resonances closer to $R_0$.

Regarding $\partial V_{\phi}/\partial \phi$ (third row in Fig. 3), it shows an almost linear decrease with $\Omega_b$ and $\beta$ for both Hercules-like<fnref n="4" /> and the Hat-like overdensities across the parameter space for the slow bar (Colrs. 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 $\partial V_{\phi}/\partial R$ of Hercules-like structures increases almost linearly with both $\Omega_b$ and $\beta$, while in Hat-like structures, $\partial V_{\phi}/\partial R$ increases with $\Omega_b$ and decreases with $\beta$. At first glance, this difference in behaviour might be related to the effect $\beta$ has on $V_{\phi}$, 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 $\beta$ (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 $\Omega_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\sigma$ 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\sigma$ 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$\sigma$ curves). In the slow bar Hercules-like overdensity (leftmost column) the azimuthal gradient matches the data for $\Omega_b \sim 40$ km s⁻¹ kpc⁻¹ and $\beta \sim 0.025$. Opposite to that, the azimuthal gradient of the fast bar Hercules-like overdensity (second column) matches the data for $\Omega_b \sim 54$ km s⁻¹ kpc⁻¹ and $\beta \sim -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 $\Omega_b$ and $\beta$ that reproduces all the measures in the data. To obtain a match in radial gradient we would need to decrease the parameters under $\Omega_b \ll 37$ km s⁻¹ kpc⁻¹, away from the match in $V_{\phi}$, and $\beta \ll -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\sigma$ when comparing models and data. This confirms that our models are sub-optimal

<footnote n="4">⁴ The Hercules-Slow Bar (leftmost column) presents an irregularity around $\beta = 0.02$ in the panels of $V_R$, $\partial V_{\phi}/\partial \phi$, and $\partial V_{\phi}/\partial R$ due to the split in the arch (see Sect. 3.2). The global trends are maintained, so this does not affect our final conclusions.</footnote>
<figure><img src="image_4.png" /><figcaption>**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 $\beta$ and $\Omega_b$. First row: $V_\phi$ of the overdensity. Second row: $V_R$ of the overdensity. Third row: $\partial V_\phi/\partial\phi$ 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 $\partial V_\phi/\partial R$.</figcaption></figure>

to reproduce the large-scale velocity distribution for any param-
eter 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 \text{ km s}^{-1}$. In Fig. 4, we present the
azimuthal velocity ($V_\phi$) of these structures across the azimuthal
range of $\pm 90^\circ$ 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)
<figure><img src="image_4.png" /><figcaption>**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 $\beta$ and $\Omega_b$. First row: $V_\phi$ of the overdensity. Second row: $V_R$ of the overdensity. Third row: $\partial V_\phi/\partial\phi$ 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 $\partial V_\phi/\partial R$.</figcaption></figure>

to reproduce the large-scale velocity distribution for any param-
eter 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 \text{ km s}^{-1}$. In Fig. 4, we present the
azimuthal velocity ($V_\phi$) of these structures across the azimuthal
range of $\pm 90^\circ$ 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)
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** $V_\phi$ of the overdensities at $V_R = -60 \text{ 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 (dashed red 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.</figcaption></figure>

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<fnref n="5" />. 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 (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_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 $\phi$, 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_\phi$ 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 \text{ km s}^{-1}$ and they are located in the opposite site of the major axis. Therefore, at $V_R = 0 \text{ 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 stabilise relatively quickly (~8 bar laps) and remain fixed, as previously anticipated in Dehnen (2000). The transient arches are joined through an entire turn around the disc forming a global structure with helical shape. In addition, this helical shape 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 \text{ 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

<footnote n="5">⁵ Notice that we are assuming a 180-deg symmetry. Therefore we are linking each structure with its counterpart.</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** $V_\phi$ of the overdensities at $V_R = -60 \text{ 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 (dashed red 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.</figcaption></figure>

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<fnref n="5" />. 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 (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_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 $\phi$, 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_\phi$ 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 \text{ km s}^{-1}$ and they are located in the opposite site of the major axis. Therefore, at $V_R = 0 \text{ 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 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 \text{ 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

<footnote n="5">⁵ Notice that we are assuming a 180-deg symmetry. Therefore we are linking each structure with its counterpart.</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Detection of the 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.</figcaption></figure>

azimuth (grey line on the right panel). We obtain a global slope for this structure of $\partial V_\phi/\partial \phi = -0.177 \text{ km s}^{-1} \text{ 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 $\phi > 0$.

Back to the models, we observed that the space between the transient arches decreases with time (Fig. 4). That is, the $\partial V_\phi/\partial \phi$ 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 $\Omega_\phi - \phi$ 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 \approx 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 $\partial V_\phi/\partial R = 28.1 \pm 2.8 \text{ km s}^{-1} \text{ kpc}^{-1}$ is compatible with recent results such as the one in Ramos et al. (2018), of $\partial V_\phi/\partial R = 26.5 \text{ km s}^{-1} \text{ 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 Hercules<fnref n="6" /> 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 ($\partial V_\phi/\partial R \sim 33 \pm 10 \text{ km s}^{-1} \text{ 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

<figure><img src="image_7.png" /><figcaption>**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).</figcaption></figure>

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\sigma$. 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 \approx 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 $\partial V_\phi/\partial \phi = -0.63 \pm 0.13 \text{ km s}^{-1} \text{ deg}^{-1}$. Using also Gaia DR3⁷, in Lucchini et al. (2024) they obtain a measurement of $-0.74 \pm 0.04 \text{ km s}^{-1} \text{ deg}^{-1}$. 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 $\partial V_\phi/\partial \phi = -0.5 \text{ km s}^{-1} \text{ deg}^{-1}$ 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 $\partial L_z/\partial \phi = -8 \text{ km s}^{-1} \text{ kpc deg}^{-1}$ at a radius of 8.2 kpc, which corresponds to $\partial V_\phi/\partial \phi = -0.96 \text{ km s}^{-1} \text{ deg}^{-1}$, when looking at the mean $V_R$ of the model of a slow bar, which roughly agrees with the observations. They

<footnote n="6">⁶ Technically, the gap between Hercules and Hyades-Pleiades.</footnote>

<footnote n="7">⁷ They use Bailer-Jones et al. (2021) distances and a $R_0 = 8.15$ kpc instead of $R_0 = 8.277$ kpc used here. This can modify the value of the measurement, but the global trends are not expected to change significantly.</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Detection of the 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.</figcaption></figure>

azimuth (grey line on the right panel). We obtain a global slope for this structure of $\partial V_\phi/\partial \phi = -0.177 \text{ km s}^{-1} \text{ 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 $\phi > 0$.

Back to the models, we observed that the space between the transient arches decreases with time (Fig. 4). That is, the $\partial V_\phi/\partial \phi$ 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 $\Omega_\phi - \phi$ 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 \approx 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 $\partial V_\phi/\partial R = 28.1 \pm 2.8 \text{ km s}^{-1} \text{ kpc}^{-1}$ is compatible with recent results such as the one in Ramos et al. (2018), of $\partial V_\phi/\partial R = 26.5 \text{ km s}^{-1} \text{ 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 Hercules<fnref n="6" /> 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 ($\partial V_\phi/\partial R \sim 33 \pm 10 \text{ km s}^{-1} \text{ 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

<figure><img src="image_7.png" /><figcaption>**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).</figcaption></figure>

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\sigma$. 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 \approx 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 $\partial V_\phi/\partial \phi = -0.63 \pm 0.13 \text{ km s}^{-1} \text{ deg}^{-1}$. Using also Gaia DR3⁷, in Lucchini et al. (2024) they obtain a measurement of $-0.74 \pm 0.04 \text{ km s}^{-1} \text{ deg}^{-1}$. 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 $\partial V_\phi/\partial \phi = -0.5 \text{ km s}^{-1} \text{ deg}^{-1}$ 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 $\partial L_z/\partial \phi = -8 \text{ km s}^{-1} \text{ kpc deg}^{-1}$ at a radius of 8.2 kpc, which corresponds to $\partial V_\phi/\partial \phi = -0.96 \text{ km s}^{-1} \text{ deg}^{-1}$, when looking at the mean $V_R$ of the model of a slow bar, which roughly agrees with the observations. They

<footnote n="6">⁶ Technically, the gap between Hercules and Hyades-Pleiades.</footnote>

<footnote n="7">⁷ They use Bailer-Jones et al. (2021) distances and a $R_0 = 8.15$ kpc instead of $R_0 = 8.277$ kpc used here. This can modify the value of the measurement, but the global trends are not expected to change significantly.</footnote>
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_\phi - \phi$ 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_\phi - \phi$ space. Thus, it should contain a maximum and a minimum in $-90 \le \phi \le 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_\phi$) 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 $\partial V_\phi / \partial R = 28 \pm 2.3$ km s⁻¹ kpc⁻¹, and the azimuthal gradient is $\partial V_\phi / \partial \phi = -0.63 \pm 0.13$ km s⁻¹ deg⁻¹.

- The radial gradient of Hercules in the BI models disagrees with the data with a significance above $3\sigma$ for all the explored parameters ($\Omega_b, \beta$, and $\phi_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.
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_\phi - \phi$ 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_\phi - \phi$ space. Thus, it should contain a maximum and a minimum in $-90 \le \phi \le 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_\phi$) 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 $\partial V_\phi / \partial R = 28 \pm 2.3$ km s⁻¹ kpc⁻¹, and the azimuthal gradient is $\partial V_\phi / \partial \phi = -0.63 \pm 0.13$ km s⁻¹ deg⁻¹.

- The radial gradient of Hercules in the BI models disagrees with the data with a significance above $3\sigma$ for all the explored parameters ($\Omega_b, \beta$, and $\phi_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 assess 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.

Acknowledgements. 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-125451NA100 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

- Anders, F., Khalatyan, A., Queiroz, A. B. A., et al. 2022, A&A, 658, A91

- Antoja, T., Figueras, F., Fernández, D., & Torra, J. 2008, A&A, 490, 135

- Antoja, T., Figueras, F., Romero-Gómez, M., et al. 2011, MNRAS, 418, 1423

- Antoja, T., Helmi, A., Dehnen, W., et al. 2014, A&A, 563, A60

- Antoja, T., Helmi, A., Romero-Gómez, M., et al. 2018, Nature, 561, 360

- Antoja, T., Ramos, P., López-Guitart, F., et al. 2022, A&A, 668, A61

- Antoja, T., Ramos, P., García-Conde, B., et al. 2023, A&A, 673, A115

- Asano, T., Fujii, M. S., Baba, J., et al. 2020, MNRAS, 499, 2416

- Asano, T., Fujii, M. S., Baba, J., et al. 2022, MNRAS, 514, 460

- Athanassoula, E. 1992, MNRAS, 259, 345

- Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Demleitner, M., & Andrae, R. 2021, AJ, 161, 147

- Bennett, M., & Bovy, J. 2019, MNRAS, 482, 1417

- Bernet, M., Ramos, P., Antoja, T., et al. 2022, A&A, 667, A116

- Binney, J. 2018, MNRAS, 474, 2706

- Binney, J., & McMillan, P. 2011, MNRAS, 413, 1889

- Bland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529

- Blomme, R., Frémat, Y., Sartoretti, P., et al. 2023, A&A, 674, A7

- Bolzano, B. 1817, Rein analytischer Beweis des Lehrsatzes daß zwischen je zwey Werthen, die ein entgegengesetzes Resultat gewähren, wenigstens eine reelle Wurzel der Gleichung liege (Gottlieb Haase)

- Chakrabarty, D. 2007, A&A, 467, 145

- Chereul, E., Creze, M., & Bienayme, O. 1998, A&A, 340, 384

- Chiba, R., & Schönrich, R. 2021, MNRAS, 505, 2412

- Chiba, R., Friske, J. K. S., & Schönrich, R. 2021, MNRAS, 500, 4710

- Clarke, J. P., & Gerhard, O. 2022, MNRAS, 512, 2171

- Contardo, G., Hogg, D. W., Hunt, J. A. S., Peek, J. E. G., & Chen, Y.-C. 2022, AJ, 164, 226

- Darling, K., & Widrow, L. M. 2019, MNRAS, 484, 1050

- Darragh-Ford, E., Hunt, J. A. S., Price-Whelan, A. M., & Johnston, K. V. 2023, ApJ, 955, 74

- Dehnen, W. 1998, AJ, 115, 2384

- Dehnen, W. 1999, AJ, 118, 1201

- Dehnen, W. 2000, AJ, 119, 800

- Eggen, O. J. 1965, in Galactic Structure, eds. A. Blaauw, & M. Schmidt (University of Chicago Press), 111

- Famaey, B., Jorissen, A., Luri, X., et al. 2005, A&A, 430, 165

- Ferrers, N. 1877, Q. J. Pure Appl. Math., 14, 1

- Font, J., Beckman, J. E., Martínez-Valpuesta, I., et al. 2017, ApJ, 835, 279

- Fragkoudi, F., Katz, D., Trick, W., et al. 2019, MNRAS, 488, 3324

- Frankel, N., Bovy, J., Tremaine, S., & Hogg, D. W. 2023, MNRAS, 521, 5917

- Friske, J. K. S., & Schönrich, R. 2019, MNRAS, 490, 5414

- Gaia Collaboration (Brown, A. G. A., et al.) 2018a, A&A, 616, A1

- Gaia Collaboration (Katz, D., et al.) 2018b, A&A, 616, A11

- Gaia Collaboration (Antoja, T., et al.) 2021, A&A, 649, A8

- Gaia Collaboration (Vallenari, A., et al.) 2023, A&A, 674, A1

- Gómez, F. A., Minchev, I., Villalobos, Á., O'Shea, B. W., & Williams, M. E. K. 2012, MNRAS, 419, 2163

- Grady, J., Belokurov, V., & Evans, N. W. 2020, MNRAS, 492, 3128

- GRAVITY Collaboration (Abuter, R., et al.) 2022, A&A, 657, L12

- Hunt, J. A. S., & Bovy, J. 2018, MNRAS, 477, 3945

- Hunt, J. A. S., Bovy, J., Pérez-Villegas, A., et al. 2018a, MNRAS, 474, 95

- Hunt, J. A. S., Hong, J., Bovy, J., Kawata, D., & Grand, R. J. J. 2018b, MNRAS, 481, 3794

- Hunt, J. A. S., Bub, M. W., Bovy, J., et al. 2019, MNRAS, 490, 1026

- Hunt, J. A. S., Stelea, I. A., Johnston, K. V., et al. 2021, MNRAS, 508, 1459

- Hunt, J. A. S., Price-Whelan, A. M., Johnston, K. V., & Darragh-Ford, E. 2022, MNRAS, 516, L7

- Kalnajs, A. J. 1991, in Dynamics of Disc Galaxies, ed. B. Sundelius (Gäteborgs: Gäteborgs University and Chalmers University of Technology), 323

- Katz, D., Sartoretti, P., Guerrier, A., et al. 2023, A&A, 674, A5

- Kawata, D., Baba, J., Ciucà, I., et al. 2018, MNRAS, 479, L108

- Khanna, S., Sharma, S., Tepper-Garcia, T., et al. 2019, MNRAS, 489, 4962

- Kushniruk, I., Schirmer, T., & Bensby, T. 2017, A&A, 608, A73

- Laporte, C. F. P., Minchev, I., Johnston, K. V., & Gómez, F. A. 2019, MNRAS, 485, 3134

- Laporte, C. F. P., Famaey, B., Monari, G., et al. 2020, A&A, 643, L3

- Li, H., & Widrow, L. M. 2021, MNRAS, 503, 1586

- Lucchini, S., Pellett, E., D'Onghia, E., & Aguerri, J. A. L. 2023, MNRAS, 519, 432

- Lucchini, S., D'Onghia, E., & Aguerri, J. A. L. 2024, MNRAS, 531, L14

- McMillan, P. J. 2017, MNRAS, 465, 76

- McMillan, P. J., Petersson, J., Tepper-Garcia, T., et al. 2022, MNRAS, 516, 4988

- Michtchenko, T. A., Lépine, J. R. D., Pérez-Villegas, A., Vieira, R. S. S., & Barros, D. A. 2018, ApJ, 863, L37

- Minchev, I., Nordhaus, J., & Quillen, A. C. 2007, ApJ, 664, L31

- Minchev, I., Quillen, A. C., Williams, M., et al. 2009, MNRAS, 396, L56

- Monari, G., Antoja, T., & Helmi, A. 2013, arXiv e-prints [arXiv:1306.2632]

- Monari, G., Helmi, A., Antoja, T., & Steinmetz, M. 2014, A&A, 569, A69

- Monari, G., Kawata, D., Hunt, J. A. S., & Famaey, B. 2017a, MNRAS, 466, L113

- Monari, G., Famaey, B., Fouvry, J.-B., & Binney, J. 2017b, MNRAS, 471, 4314

- Monari, G., Famaey, B., Siebert, A., Wegg, C., & Gerhard, O. 2019a, A&A, 626, A41

- Monari, G., Famaey, B., Siebert, A., et al. 2019b, A&A, 632, A107

- Moore, E. F. 1959, Proc. Int. Symp. Switching Theory, 1959, 285

- Pérez-Villegas, A., Portail, M., Wegg, C., & Gerhard, O. 2017, ApJ, 840, L2

- Portail, M., Gerhard, O., Wegg, C., & Ness, M. 2017, MNRAS, 465, 1621

- Raboud, D., Grenon, M., Martinet, L., Fux, R., & Udry, S. 1998, A&A, 335, L61

- Ramos, P., Antoja, T., & Figueras, F. 2018, A&A, 619, A72

- Reid, M. J., & Brunthaler, A. 2020, ApJ, 892, 39

- Robin, A. C., Bienaymé, O., Salomon, J. B., et al. 2022, A&A, 667, A98
- 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 assess 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.

Acknowledgements. 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-100 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

- Anders, F., Khalatyan, A., Queiroz, A. B. A., et al. 2022, A&A, 658, A91

- Antoja, T., Figueras, F., Fernández, D., & Torra, J. 2008, A&A, 490, 135

- Antoja, T., Figueras, F., Romero-Gómez, M., et al. 2011, MNRAS, 418, 1423

- Antoja, T., Helmi, A., Dehnen, W., et al. 2014, A&A, 563, A60

- Antoja, T., Helmi, A., Romero-Gómez, M., et al. 2018, Nature, 561, 360

- Antoja, T., Ramos, P., López-Guitart, F., et al. 2022, A&A, 668, A61

- Antoja, T., Ramos, P., García-Conde, B., et al. 2023, A&A, 673, A115

- Asano, T., Fujii, M. S., Baba, J., et al. 2020, MNRAS, 499, 2416

- Asano, T., Fujii, M. S., Baba, J., et al. 2022, MNRAS, 514, 460

- Athanassoula, E. 1992, MNRAS, 259, 345

- Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Demleitner, M., & Andrae, R. 2021, AJ, 161, 147

- Bennett, M., & Bovy, J. 2019, MNRAS, 482, 1417

- Bernet, M., Ramos, P., Antoja, T., et al. 2022, A&A, 667, A116

- Binney, J. 2018, MNRAS, 474, 2706

- Binney, J., & McMillan, P. 2011, MNRAS, 413, 1889

- Bland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529

- Blomme, R., Frémat, Y., Sartoretti, P., et al. 2023, A&A, 674, A7

- Bolzano, B. 1817, *Rein analytischer Beweis des Lehrsatzes daß zwischen je zwey Werthen, die ein entgegengesetzes Resultat gewähren, wenigstens eine reelle Wurzel der Gleichung liege* (Gottlieb Haase)

- Chakrabarty, D. 2007, A&A, 467, 145

- Chereul, E., Creze, M., & Bienayme, O. 1998, A&A, 340, 384

- Chiba, R., & Schönrich, R. 2021, MNRAS, 505, 2412

- Chiba, R., Friske, J. K. S., & Schönrich, R. 2021, MNRAS, 500, 4710

- Clarke, J. P., & Gerhard, O. 2022, MNRAS, 512, 2171

- Contardo, G., Hogg, D. W., Hunt, J. A. S., Peek, J. E. G., & Chen, Y.-C. 2022, AJ, 164, 226

- Darling, K., & Widrow, L. M. 2019, MNRAS, 484, 1050

- Darragh-Ford, E., Hunt, J. A. S., Price-Whelan, A. M., & Johnston, K. V. 2023, ApJ, 955, 74

- Dehnen, W. 1998, AJ, 115, 2384

- Dehnen, W. 1999, AJ, 118, 1201

- Dehnen, W. 2000, AJ, 119, 800

- Eggen, O. J. 1965, in *Galactic Structure*, eds. A. Blaauw, & M. Schmidt (University of Chicago Press), 111

- Famaey, B., Jorissen, A., Luri, X., et al. 2005, A&A, 430, 165

- Ferrers, N. 1877, Q. J. Pure Appl. Math., 14, 1

- Font, J., Beckman, J. E., Martínez-Valpuesta, I., et al. 2017, ApJ, 835, 279

- Fragkoudi, F., Katz, D., Trick, W., et al. 2019, MNRAS, 488, 3324

- Frankel, N., Bovy, J., Tremaine, S., & Hogg, D. W. 2023, MNRAS, 521, 5917

- Friske, J. K. S., & Schönrich, R. 2019, MNRAS, 490, 5414

- Gaia Collaboration (Brown, A. G. A., et al.) 2018a, A&A, 616, A1

- Gaia Collaboration (Katz, D., et al.) 2018b, A&A, 616, A11

- Gaia Collaboration (Antoja, T., et al.) 2021, A&A, 649, A8

- Gaia Collaboration (Vallenari, A., et al.) 2023, A&A, 674, A1

- Gómez, F. A., Minchev, I., Villalobos, Á., O'Shea, B. W., & Williams, M. E. K. 2012, MNRAS, 419, 2163

- Grady, J., Belokurov, V., & Evans, N. W. 2020, MNRAS, 492, 3128

- GRAVITY Collaboration (Abuter, R., et al.) 2022, A&A, 657, L12

- Hunt, J. A. S., & Bovy, J. 2018, MNRAS, 477, 3945

- Hunt, J. A. S., Bovy, J., Pérez-Villegas, A., et al. 2018a, MNRAS, 474, 95

- Hunt, J. A. S., Hong, J., Bovy, J., Kawata, D., & Grand, R. J. J. 2018b, MNRAS, 481, 3794

- Hunt, J. A. S., Bub, M. W., Bovy, J., et al. 2019, MNRAS, 490, 1026

- Hunt, J. A. S., Stelea, I. A., Johnston, K. V., et al. 2021, MNRAS, 508, 1459

- Hunt, J. A. S., Price-Whelan, A. M., Johnston, K. V., & Darragh-Ford, E. 2022, MNRAS, 516, L7

- Kalnajs, A. J. 1991, in *Dynamics of Disc Galaxies*, ed. B. Sundelius (Gäteborgs: Gäteborgs University and Chalmers University of Technology), 323

- Katz, D., Sartoretti, P., Guerrier, A., et al. 2023, A&A, 674, A5

- Kawata, D., Baba, J., Ciucà, I., et al. 2018, MNRAS, 479, L108

- Khanna, S., Sharma, S., Tepper-Garcia, T., et al. 2019, MNRAS, 489, 4962

- Kushniruk, I., Schirmer, T., & Bensby, T. 2017, A&A, 608, A73

- Laporte, C. F. P., Minchev, I., Johnston, K. V., & Gómez, F. A. 2019, MNRAS, 485, 3134

- Laporte, C. F. P., Famaey, B., Monari, G., et al. 2020, A&A, 643, L3

- Li, H., & Widrow, L. M. 2021, MNRAS, 503, 1586

- Lucchini, S., Pellett, E., D'Onghia, E., & Aguerri, J. A. L. 2023, MNRAS, 519, 432

- Lucchini, S., D'Onghia, E., & Aguerri, J. A. L. 2024, MNRAS, 531, L14

- McMillan, P. J. 2017, MNRAS, 465, 76

- McMillan, P. J., Petersson, J., Tepper-Garcia, T., et al. 2022, MNRAS, 516, 4988

- Michtchenko, T. A., Lépine, J. R. D., Pérez-Villegas, A., Vieira, R. S. S., & Barros, D. A. 2018, ApJ, 863, L37

- Minchev, I., Nordhaus, J., & Quillen, A. C. 2007, ApJ, 664, L31

- Minchev, I., Quillen, A. C., Williams, M., et al. 2009, MNRAS, 396, L56

- Monari, G., Antoja, T., & Helmi, A. 2013, arXiv e-prints [arXiv:1306.2632]

- Monari, G., Helmi, A., Antoja, T., & Steinmetz, M. 2014, A&A, 569, A69

- Monari, G., Kawata, D., Hunt, J. A. S., & Famaey, B. 2017a, MNRAS, 466, L113

- Monari, G., Famaey, B., Fouvry, J.-B., & Binney, J. 2017b, MNRAS, 471, 4314

- Monari, G., Famaey, B., Siebert, A., Wegg, C., & Gerhard, O. 2019a, A&A, 626, A41

- Monari, G., Famaey, B., Siebert, A., et al. 2019b, A&A, 632, A107

- Moore, E. F. 1959, Proc. Int. Symp. Switching Theory, 1959, 285

- Pérez-Villegas, A., Portail, M., Wegg, C., & Gerhard, O. 2017, ApJ, 840, L2

- Portail, M., Gerhard, O., Wegg, C., & Ness, M. 2017, MNRAS, 465, 1621

- Raboud, D., Grenon, M., Martinet, L., Fux, R., & Udry, S. 1998, A&A, 335, L61

- Ramos, P., Antoja, T., & Figueras, F. 2018, A&A, 619, A72

- Reid, M. J., & Brunthaler, A. 2020, ApJ, 892, 39

- Robin, A. C., Bienaymé, O., Salomon, J. B., et al. 2022, A&A, 667, A98
Rybizki, J., Green, G. M., Rix, H.-W., et al. 2022, MNRAS, 510, 2597

Sanders, J. L., Smith, L., & Evans, N. W. 2019, MNRAS, 488, 4552

Sanders, J. L., Matsunaga, N., Kawata, D., et al. 2022, MNRAS, 517, 257

Skuljan, J., Hearnshaw, J. B., & Cottrell, P. L. 1999, MNRAS, 308, 731

Sormani, M. C., Binney, J., & Magorrian, J. 2015, MNRAS, 454, 1818

Steinmetz, M., Zwitter, T., Siebert, A., et al. 2006, AJ, 132, 1645

Strömberg, G. 1946, ApJ, 104, 12

Trick, W. H. 2022, MNRAS, 509, 844

Trick, W. H., Fragkoudi, F., Hunt, J. A. S., Mackereth, J. T., & White, S. D. M. 2021, MNRAS, 500, 2645

Tsoutsis, P., Kalapotharakos, C., Efthymiopoulos, C., & Contopoulos, G. 2009, A&A, 495, 743

Vasiliev, E. 2019, MNRAS, 482, 1525

Vauterin, P., & Dejonghe, H. 1997, MNRAS, 286, 812
Rybizki, J., Green, G. M., Rix, H.-W., et al. 2022, MNRAS, 510, 2597

Sanders, J. L., Smith, L., & Evans, N. W. 2019, MNRAS, 488, 4552

Sanders, J. L., Matsunaga, N., Kawata, D., et al. 2022, MNRAS, 517, 257

Skuljan, J., Hearnshaw, J. B., & Cottrell, P. L. 1999, MNRAS, 308, 731

Sormani, M. C., Binney, J., & Magorrian, J. 2015, MNRAS, 454, 1818

Steinmetz, M., Zwitter, T., Siebert, A., et al. 2006, AJ, 132, 1645

Strömberg, G. 1946, ApJ, 104, 12

Trick, W. H. 2022, MNRAS, 509, 844

Trick, W. H., Fragkoudi, F., Hunt, J. A. S., Mackereth, J. T., & White, S. D. M. 2021, MNRAS, 500, 2645

Tsoutsis, P., Kalapotharakos, C., Efthymiopoulos, C., & Contopoulos, G. 2009, A&A, 495, 743

Vasiliev, E. 2019, MNRAS, 482, 1525

Vauterin, P., & Dejonghe, H. 1997, MNRAS, 286, 812
## 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

<table><caption>**Table A.1.** Mean slope and the dispersion within the box (see Fig. 2) in both models for *t* = 14 bar laps.</caption><thead><tr><th></th><th></th><th>∂V<sub>φ</sub>/∂R (σ)</th><th>∂V<sub>φ</sub>/∂φ (σ)</th></tr></thead><tbody><tr><td rowspan="2">Hercules</td><td>39</td><td>36.2 (1.8)</td><td>-0.43 (0.11)</td></tr><tr><td>56</td><td>40.5 (5.7)</td><td>-0.62 (0.13)</td></tr><tr><td rowspan="2">Arch/Hat</td><td>39</td><td>24.5 (0.9)</td><td>-0.13 (0.06)</td></tr><tr><td>56</td><td>24.6 (1.1)</td><td>-0.10 (0.03)</td></tr></tbody></table>

$V_R = -60 \text{ km s}^{-1}$ 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.

<figure><img src="image_8.png" /><figcaption>**Fig. A.1.** Overdensities' detection in the simulations, for *t* = 14 bar laps. They are coloured by their $\partial V_{\phi}/\partial R$ (top panels) and $\partial V_{\phi}/\partial \phi$ (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 overdensities.</figcaption></figure>
## 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

<table><caption>**Table A.1.** Mean slope and the dispersion within the box (see Fig. 2) in both models for *t* = 14 bar laps.</caption><thead><tr><th></th><th></th><th>∂V<sub>φ</sub>/∂R (σ)</th><th>∂V<sub>φ</sub>/∂φ (σ)</th></tr></thead><tbody><tr><td rowspan="2">Hercules</td><td>39</td><td>36.2 (1.8)</td><td>-0.43 (0.11)</td></tr><tr><td>56</td><td>40.5 (5.7)</td><td>-0.62 (0.13)</td></tr><tr><td rowspan="2">Arch/Hat</td><td>39</td><td>24.5 (0.9)</td><td>-0.13 (0.06)</td></tr><tr><td>56</td><td>24.6 (1.1)</td><td>-0.10 (0.03)</td></tr></tbody></table>

$V_R = -60 \text{ km s}^{-1}$ 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.

<figure><img src="image_8.png" /><figcaption>**Fig. A.1.** Overdensities' detection in the simulations, for *t* = 14 bar laps. They are coloured by their $\partial V_{\phi}/\partial R$ (top panels) and $\partial V_{\phi}/\partial \phi$ (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 overdensities.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Fig. A.2.** Time evolution of the structures at $V_R = -60 \text{ 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.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Fig. A.2.** Time evolution of the structures at $V_R = -60 \text{ 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.</figcaption></figure>
## 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 $\phi_b$. We conclude that the trends observed at $\phi_b = -30$ deg (Fig. 3) are maintained when looking at other azimuths. The only notice-

able change is the steepness of the $\partial V_\phi/\partial R$ and $\partial V_\phi/\partial \phi$ panels according to the distance to the major axis of the bar. In Fig. B.1 ($\phi_b = -20$ deg) the contour lines are close to each other, indicating a rapid change of the properties. Opposite to that, in Fig. B.2 ($\phi_b = -40$ deg), the contour lines are more separated, indicating a solution which is more stable.

<figure><img src="image_10.png" /><figcaption>Fig. B.1. Same as Fig. 3 but computed with $\phi_b = -20$ deg.</figcaption></figure>
**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 $\phi_b$. We
conclude that the trends observed at $\phi_b = -30$ deg (Fig. 3) are
maintained when looking at other azimuths. The only notice-

able change is the steepness of the $\partial V_{\phi}/\partial R$ and $\partial V_{\phi}/\partial \phi$ panels
according to the distance to the major axis of the bar. In Fig. B.1
($\phi_b = -20$ deg) the contour lines are close to each other, indicat-
ing a rapid change of the properties. Opposite to that, in Fig. B.2
($\phi_b = -40$ deg), the contour lines are more separated, indicating
a solution which is more stable.

<figure><img src="image_10.png" /><figcaption>Fig. B.1. Same as Fig. 3 but computed with $\phi_b = -20$ deg.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. B.2. Same as Fig. 3 but computed with $\phi_b = -40$ deg.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. B.2. Same as Fig. 3 but computed with $\phi_b = -40$ deg.</figcaption></figure>
## Appendix C: Model versus model degeneracy

The computation of the Galactic models links a set of parameters ($\Omega_b, \beta, \phi_b$...) to a set of measurables ($V_R, V_\phi, \partial V_\phi/\partial R, \partial V_\phi/\partial \phi$). 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 ($\sigma$, 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_\phi$. When restricted to Hercules and its local value of $V_\phi$ alone (black contour, top panel), we observe a substantial region of parameter space where the solution is degenerate, with multiple combinations of $\Omega_b$ and $\beta$ from the FBM models yielding the same $V_\phi$. However, as we introduce additional measurables, such as $\partial V_\phi/\partial \phi$ (orange contour) and $\partial V_\phi/\partial R$ (green contour), the statistical difference becomes more pronounced, scaling to over $4\sigma$. When we further include a second structure like Hat (bottom panel), the degeneracy in the solution can be partially resolved using only $V_\phi$. 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_\phi$ of both structures we already see that we are able to constrain a small region of $\Omega_b$ and $\beta$. In addition, including the azimuthal gradient brings the

<figure><img src="image_12.png" /><figcaption>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_\phi$ measurements. Orange contours are for the disagreement considering $V_\phi$ and $\partial V_\phi/\partial \phi$. Green contours show the disagreement considering $V_\phi$, $\partial V_\phi/\partial \phi$, and $\partial V_\phi/\partial R$.</figcaption></figure>

statistical significance over $8\sigma$ in almost the entire parameter space.
# Appendix C: Model versus model degeneracy

The computation of the Galactic models links a set of parameters ($\Omega_b, \beta, \phi_b$...) to a set of measurables ($V_R, V_\phi, \partial V_\phi/\partial R, \partial V_\phi/\partial \phi$). 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 ($\sigma$, 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_\phi$. When restricted to Hercules and its local value of $V_\phi$ alone (black contour, top panel), we observe a substantial region of parameter space where the solution is degenerate, with multiple combinations of $\Omega_b$ and $\beta$ from the FBM models yielding the same $V_\phi$. However, as we introduce additional measurables, such as $\partial V_\phi/\partial \phi$ (orange contour) and $\partial V_\phi/\partial R$ (green contour), the statistical difference becomes more pronounced, scaling to over $4\sigma$. When we further include a second structure like Hat (bottom panel), the degeneracy in the solution can be partially resolved using only $V_\phi$. 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_\phi$ of both structures we already see that we are able to constrain a small region of $\Omega_b$ and $\beta$. In addition, including the azimuthal gradient brings the

<figure><img src="image_12.png" /><figcaption>**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_\phi$ measurements. Orange contours are for the disagreement considering $V_\phi$ and $\partial V_\phi/\partial \phi$. Green contours show the disagreement considering $V_\phi$, $\partial V_\phi/\partial \phi$, and $\partial V_\phi/\partial R$.</figcaption></figure>

statistical significance over $8\sigma$ in almost the entire parameter space.
<figure><img src="image_13.png" /></figure>

**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 the Hercules- and Hat-like one (bottom).
<figure><img src="image_13.png" /></figure>

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