The Milky Way Bar Pattern Speed using Hercules and Gaia DR3

The distribution of moving groups in the solar neighborhood has been used to constrain dynamical properties of the Milky Way for decades. The kinematic bimodality between the main mode (Hyades, Pleiades, Coma Berenices, and Sirius) and Hercules can be explained by two different bar models -- via the outer Lindblad resonance of a bar with a high pattern speed ($\sim$55 km s$^{-1}$ kpc$^{-1}$), or via the corotation resonance of a bar with a low pattern speed ($\sim$40 km s$^{-1}$ kpc$^{-1}$). Recent works directly studying the kinematics of bar stars and gas flows near the center of the Galaxy have converged on the low pattern speed model. In this paper, we independently confirm this result by using Gaia DR3 to directly study the variation of Hercules across Galactic azimuth. We find that Hercules increases in $V_\phi$ and becomes stronger as we move towards the minor axis of the bar, and decreases in $V_\phi$ and becomes weaker as we move towards the major axis of the bar. This is in direct agreement with theoretical predictions of a low pattern speed bar model in which Hercules is formed by the corotation resonance with stars orbiting the bar's L4/L5 Lagrange points.


INTRODUCTION
Being embedded within the Milky Way (MW) makes it difficult to observe its properties.We must use indirect methods to determine the distribution of mass, and corresponding nonaxisymmetric features in our Galaxy.Specifically, determining the properties of our own Galactic bar -an oblong stellar overdensity at the center of our Galaxy -have required several independent techniques to converge on a consensus.Three main methods that have been used are: direct observations of the kinematics of bar stars, gas flows in the inner Galaxy, and the kinematics of stars in the solar neighborhood (i.e."moving groups").
One of the first attempts to directly measure the bar pattern speed used a version of the Tremaine-Weinberg method (Tremaine & Weinberg 1984) modified for use with radial velocities (Debattista et al. 2002).Using OH/IR stars, they found a high pattern speed of Ω  ∼ 60 km s −1 kpc −1 .However recently, more advanced MW surveys have allowed us to directly observe the transverse kinematics of stars towards the center of our Galaxy (Portail et al. 2017;Sanders et al. 2019;Clarke & Gerhard 2022), and these studies all support at low pattern speed of Ω  ∼ 30 − 40 km s −1 kpc −1 .
Using the full 3D velocities of stars near the Sun, we can also investigate the inhomogeneities in the velocity distributions to learn about the nonaxisymmetric features of our Galaxy.How these stars cluster in velocity space (moving groups) contains a wealth of information about the evolution of our Galaxy's disk and the forces these stars are feeling.Specifically, the bimodality in the Galactocentric azimuthal velocity versus radial velocity plot has been one of the main features models have attempted to reproduce.The main mode (at   ∼ −230 km s −1 ) contains the moving groups of Hyades, Pleiades, Coma Berenices, and Sirius, while Hercules (at   ∼ −200 km s −1 ) is separated by a gap (a strong underdensity).This bimodality has been explained through resonances of the MW's bar.Since the bar contributes a nonaxisymmetric gravitational potential, it has an associated pattern speed, as discussed above.Over time, the bar perturbs the stars to align the frequencies of the stellar orbits with the bar's frequency.These stars are "in resonance" with the bar.A bar's strongest resonances are the corotation resonance (CR; in which the stellar and bar frequencies match exactly), and the inner and outer Lindblad resonances (ILR and OLR; in which the star completes two radial oscillations for every orbit around the galaxy).Unfortunately, depending on the pattern speed of the bar, these resonances fall at different locations, and multiple models are able to explain this bimodality.For a high pattern speed bar (Ω  ∼ 55 km s −1 kpc −1 ), the OLR falls just inside the solar neighborhood with Hercules being formed by stars moving outward in the disk coming from inside  OLR (Dehnen 2000;Minchev et al. 2010;Antoja et al. 2014;Fragkoudi et al. 2019).For a low pattern speed bar (Ω  ∼ 40 km s −1 kpc −1 ), the CR is able to create Hercules through stars orbiting the bar's Lagrange points (Pérez-Villegas et al. 2017;Hunt et al. 2018;Monari et al. 2019a;Asano et al. 2020;D'Onghia & L. Aguerri 2020).While most modern models are more consistent with a low pattern speed, a direct observation using stellar kinematics indicating Hercules having been formed via CR or via OLR has remained elusive.
Recent works explored various bar's resonances to discriminate between the models.Monari et al. (2019a) analytically explored the impact of the  = 3, 4, and 6 modes in addition to the  = 2 CR and OLR modes.Finding that many of these other resonances aligned with ridges and structures in the local velocity plane, they determined that a low pattern speed bar (Ω  = 39 km s −1 kpc −1 ) provides the best fit.Furthermore, analysis of local gas flows and metallicity gradients also indicate that Hercules has been formed by CR (Binney 2020;Chiba & Schönrich 2021).An in-depth study of the effect of resonance trapping on Hercules also concludes that a slowing bar with a current pattern speed of 35 km s −1 kpc −1 provides the best fit to the data (Chiba et al. 2021).Additionally, Antoja et al. (2014) was able to use RAVE data to explore the behavior of Hercules beyond the solar neighborhood varying radius.While they find that their data is more consistent with a high pattern speed bar model (Ω  = 54 km s −1 kpc −1 ), there is little difference in the R-dependence for the OLR and the CR.
While both the CR and the OLR have similar signatures in the   −   plane in the solar neighborhood, analytical calculations and numerical simulations have shown that the predictions of these two models differ as we move throughout the Galactic disk (Monari et al. 2019b;D'Onghia & L. Aguerri 2020).With Gaia's latest data release (Data Release 3; Gaia Collaboration et al. 2016, 2018, 2022a), we can begin to differentiate between these two models by looking beyond the solar neighborhood in azimuth.In this paper, we use Gaia's full 6D phase space information to explore how Hercules changes as we move around the disk in , independently confirming the pattern speed of the MW's bar.

METHODS
We follow the methods outlined in Lucchini et al. (2023) (hereafter Paper I) for Gaia selection and analysis using the wavelet transformation.Figure 1 shows the solar neighborhood velocity plane and its corresponding wavelet transformed image.The plus signs denote the five most significant classical moving groups as detected by local maxima in the wavelet image.Starting with the 33,653,049 stars with radial velocities and geometric distances computed by Bailer-Jones et al. (2021) (Bailer-Jones et al. 2020), we transformed the six-dimensional data into Galactocentric cylindrical coordinates 1 .See Figure 2 for a schematic defining our coordinate system with  increasing away from the Galactic center,  increasing counter to the direction of rotation, and  increasing towards the Galactic north pole.With the Sun located at  0 = 0 • , this means the major axis of the Milky Way's bar is at  ∼ −20 • (Gaia Collaboration et al. 2022b).Figure 2 also shows the extent of the Gaia DR3 data that we used: 6.5 <  < 10 kpc, and −15 • <  < 15 • .As in Paper I, we used the pyia2 code to propagate the errors (including correlations) from the source properties (right ascension, declination, proper motions, and radial velocities) to the final properties (Galactocentric cylindrical coordinates; Price-Whelan 2018).
We use the wavelet transform code, MGwave3 (described in Paper I), to analyze the velocity distributions of these different neighborhoods.The wavelet transformation is well suited to this problem since it allows us to isolate structures of a given size in an image (in this case, a histogram).Each velocity plane histogram (  −  ) is transformed using the Starlet transform (Starck & Murtagh 1994;Starck & Pierre 1998;Starck & Murtagh 2006) and the relative significance of the results are evaluated with respect to Poisson noise (Slezak et al. 1993) and errors in the source Gaia data (using Monte Carlo simulations).In this paper, we used a wavelet scale of 8 − 16 km s −1 to best identify the Hercules group.See Paper I for more details on the wavelet code and methodology.

RESULTS
We track Hercules across the Galactic disk using Gaia DR3.In particular, we are looking for variations in the size and strength of the Hercules moving group as we vary the azimuth.Figure 3 shows the velocity plane for five different neighborhoods:  = ±15 • , ±9 • , and 0 • (corresponding to the solar neighborhood shown in Figure 1).These angles were chosen to show the full extent of usable data from Gaia DR3 (±15 • , where the number of stars per bin reaches ∼ 10 4 ), while showing an intermediate region in which Hercules is just beginning to merge with the main mode (−9 • ).The upper panels are towards the major axis of the bar (−), and the lower panels are towards the minor axis of the bar (+).The left column shows the   −  histogram of all the stars in the specified neighborhood.The center column shows the wavelet transform of these data in which purple represents overdensities and green represents underdensities (as in Figure 1).The right-most column shows a plot of wavelet coefficient as a function of   obtained by summing the wavelet transformed image along   (the "1D summed wavelet histogram").From these plots, we are able to clearly see the location of Hercules in each neighborhood as the peak in this histogram around   ∼ −200 km s −1 .
Figure 4 summarizes the properties of Hercules as a function of azimuth.Figure 4a shows the azimuthal velocity of Hercules as a function of .This is identified by the location of the peak in the 1D summed wavelet histogram (right panels in Figure 3).A linear fit to   of Hercules versus azimuth gives  ,Herc = −0.74 − 199 (with  ,Herc in km s −1 and  in degrees) increasing from 190 km s −1 at  = −15 • to 211 km s −1 at  = +15 • .In Figure 4b, we show the percentage of stars that constitute Hercules in each neighborhood.This is determined by selecting all the stars within ±15 km s −1 of the   of Hercules (Figure 4a) and dividing by the total number of stars in the given bin.These 30 km s −1 regions are shown bounded by the dashed horizontal lines in Figure 3.In the solar neighborhood, Hercules constitutes 20.0% of stars.As we move towards the bar's major axis (upwards in Figure 3), this percentage decreases to 15.0% at  = −15 • .As we move towards the bar's minor axis (downwards in Figure 3), this percentage increases to 35.7% at  = +15 • .
Due to the increase in  ,Herc (Figure 4a), one would naturally expect the fraction of stars in Hercules (within ±15 km s −1 of  ,Herc ) to increase since it will get closer to the main mode (  ∼ −230 km s −1 at  =  0 ).However, by comparing against an unperturbed Dehnen distribution function (Dehnen 1999;Hunt et al. 2018) we see that  ,Herc only increases the Hercules fraction to 25% compared with 36% seen in the data.Therefore we conclude that the increase in  ,Herc cannot account for the increase in the Hercules fraction seen in Figure 4b and must be due to the resonance interactions.
If Hercules was formed through interaction with the Galactic bar's outer Lindblad resonance, we would expect its strength to be relatively constant across azimuth (Fragkoudi et al. 2019, see their Figure 12, left panel showing "Region 1" corresponding to Hercules).However, if the corotation resonance is responsible for Hercules, its member stars would be orbiting around the Galactic bar's L4/L5 Lagrange points (Pérez-Villegas et al. 2017; D'Onghia & L. Aguerri 2020) located along the bar's minor axis.Therefore, we should expect to see a larger fraction of stars in Hercules as we approach L4/L5 (in the + direction; see Figure 4 from D'Onghia & L. Aguerri 2020).The behavior identified here in Gaia DR3, mimics exactly the predictions of D'Onghia & L. Aguerri (2020) (see their Figure 8).Towards the bar's major axis, Hercules diminishes (for both the fast and slow bar scenarios).However, towards the minor axis, we are able to discriminate between these two models.For the high Ω  bar, Hercules should remain subdominant and separated from the main mode.While in the low Ω  bar model, Hercules should become extremely prominent, even merging with the main mode, which is what we see in the data.
Figure 4b shows that Hercules increases in strength as we move towards the bar's minor axis (+).This panel also shows the fraction of stars within 175 − 205 km s −1 at  = 12 kpc (corresponding to / OLR = 0.9 for a ∼40 km s −1 kpc −1 bar, to compare against Dehnen 2000) to show the angular dependence due to the OLR.As expected, the fraction of stars trapped by the OLR (at  = 12 kpc) is relatively constant in azimuth, while at the solar radius, there is a sharp increase in the + direction, indicating that Hercules is comprised of stars trapped at the Lagrange points.
Moreover, as shown in Monari et al. (2019b), we would expect the angular momentum of Hercules to vary with azimuth.Figure 3 shows that the Hercules overdensity in the data changes its value of Hercules is denoted by a horizontal line across each plot in a row, identified as a peak in the summed wavelet histogram (right column).Additionally, the percentage of stars that constitute Hercules (defined as all stars within 15 km s −1 of Hercules'   ) is printed in the top left of the left plots.The total number of stars in each bin, N, is also listed in the bottom left of the left plots.This corresponds to mean number densities of 0.2, 1.02, 11.69, 0.96, and 0.24 pc −2 in each bin from top to bottom.In the direction of rotation (i.e.− , towards the major axis of the MW bar, top panels) Hercules diminishes and the percentage drops, while moving towards the minor axis of the bar (+, towards the Lagrange points, bottom panels) Hercules grows and merges in with the main mode.We additionally see that   of Hercules increases as we approach the Lagrange points, as predicted by models (D'Onghia & L. Aguerri 2020; Monari et al. 2019b; see also Figure 4a).et al. (2019b).However, it is clear that the slope is < 0, consistent with Hercules being formed by corotation, not the OLR.
This was also explored in Bernet et al. (2022) where they found a slope of −0.5 km s −1 deg −1 for the Hercules group at  = 8 kpc.Their analysis focused on the detection of arches in kinematic space using a 1D wavelet transformation on data binned along   .Additionally, their technique to isolate the groups differs from ours which could lead to the slight discrepancy between our result of −0.74 km s −1 deg −1 and theirs.However, as with the result from Monari et al. (2019b), it is encouraging that we are consistently finding a significant negative slope.
Figure 3 shows that at  = ±15 • , the number of stars available from Gaia is ∼2% that of  = 0 • .Therefore the sample of stars will be very different between these regions due to the Gaia selection function.Since we are mainly observing the brightest stars at the large distances, the different stellar population observed may affect the properties of Hercules.However, as the fainter stars are usually older (and therefore have had more time to be perturbed by the angular momentum versus azimuth.In their model the Sun is located at 8.2 kpc giving us a predicted slope of −8/8.2 = −0.96km s −1 deg −1 .the Galaxy's gravitational potential), we would expect that Hercules would become stronger as we include more faint stars in our sample.This can be tested with Gaia DR4 to see how the inclusion of fainter stars at these limits affects Hercules.While the variation across azimuth is what allows us to discriminate between the low and high Ω  bar models, we also briefly explore the change in Hercules'   as a function of radius to compare with Antoja et al. (2014).Figure 5 shows our calculated   for Hercules as a function of Galactocentric radius compared against the four data points in Antoja et al. (2014) from RAVE.We followed the same process as in that paper using  = 6 • (closer to the bar minor axis; see their Figure 8), and determining  ,Herc by finding the minimum of the 1D summed wavelet histogram between the main mode and the Hercules mode.Our measured azimuthal velocities are consistent with those found in Antoja et al. (2014) (with the largest deviation being 1.2 at 8.3 kpc), however since their exploration exclusively looked at signatures of the OLR from a high Ω  bar, this doesn't rule out the low Ω  bar model.

CONCLUSIONS
Here we have used Gaia DR3 to track the properties of Hercules through Galactic azimuth.Varying Galactocentric  15 • either side of the Sun, we see that there is a strong variation in the azimuthal velocity and strength of Hercules.Hercules constitutes a larger fraction of stars as we move towards the minor axis of the bar (+).This is in direct agreement with predictions of a low pattern speed bar model in which Hercules is formed through stars trapped at corotation, orbiting the L4/L5 Lagrange points.This independently corroborates several recent works converging on the slow bar model (Monari et al. 2019b;Binney 2020;Chiba & Schönrich 2021;Chiba et al. 2021;Gaia Collaboration et al. 2022b).Future analytic work exploring the effect of the bar on Hercules across a full parameter space of models will also give us more direct constraints on its properties.
With the next release from Gaia, DR4, we can expect further significant improvements using this technique.While DR3 extended our view out to  = ±15 • , we can hope to get similar signal to noise results out to  = ±30 • and beyond.This will allow us to get a direct measurement of the bar angle with respect to the Sun's position, by looking for a minimum in Hercules' strength and   .

Figure 1 .Figure 2 .
Figure 1.Panel a shows the 2D histogram of all the solar neighborhood stars in the   −   kinematic plane using Gaia DR3 with a bin size of 0.5 km s −1 .Panel b shows the resultant wavelet transformed image at a scale of 8 − 16 km s −1 .The purple and green regions depict the relative strength of the positive and negative wavelet coefficients, respectively.Contours are shown at the −0.1% (dashed), 5%, 15%, 50%, and 90% levels.The plus signs denote the locations of the five major moving groups as detected by our wavelet transformation -Hyades, Pleiades, Coma Berenices, Sirius, and Hercules.

Figure 3 .
Figure 3.The velocity plane and its wavelet transformation for five neighborhoods located at  = −15 • , −9 • , 0 • , 9 • , and 15 • . = 0 • corresponds to the location of the Sun, and negative  is in the direction of rotation, while positive  is counter to the direction of rotation.The left column shows the   −   histogram for each region, the center column shows the wavelet transformed kinematic plane as in Figure 1, and the right column shows the wavelet histogram summed along   .In the left column, the number of bins in the histogram has been scaled to the total number of stars such that we use 200 bins in each dimension at  = ±15 • , and we use 600 bins at  = 0 • .Hercules is denoted by a horizontal line across each plot in a row, identified as a peak in the summed wavelet histogram (right column).Additionally, the percentage of stars that constitute Hercules (defined as all stars within 15 km s −1 of Hercules'   ) is printed in the top left of the left plots.The total number of stars in each bin, N, is also listed in the bottom left of the left plots.This corresponds to mean number densities of 0.2, 1.02, 11.69, 0.96, and 0.24 pc −2 in each bin from top to bottom.In the direction of rotation (i.e.− , towards the major axis of the MW bar, top panels) Hercules diminishes and the percentage drops, while moving towards the minor axis of the bar (+, towards the Lagrange points, bottom panels) Hercules grows and merges in with the main mode.We additionally see that   of Hercules increases as we approach the Lagrange points, as predicted by models (D'Onghia & L. Aguerri 2020;Monari et al. 2019b; see also Figure4a).

Figure 4 .
Figure 4. Properties of Hercules as a function of azimuth.The top panel shows the mean   for the Hercules group as detected using the location of the peak in the wavelet histogram.A linear fit with the 95% confidence interval is shown as a dashed line and grey shading.The best fit equation is shown in the top left.The bottom panel shows the percentage of stars in the given neighborhood that have   values within 15 km s −1 of Hercules.The dotted line shows the percentage of stars within 175 − 205 km s −1 at  = 12 kpc (such that / OLR = 0.9 for a ∼40 km s −1 kpc −1 bar, to compare against Dehnen 2000).At the OLR, we see no variation in intensity as a function of .These two plots clearly show that in the + direction (towards the bar's minor axis), Hercules sees an increase in angular momentum and becomes more dominant as predicted by the long, slow bar model (D'Onghia & L. Aguerri 2020).

Figure 5 .
Figure 5. Hercules azimuthal velocity as a function of Galactocentric radius.Following the technique outlined in Antoja et al. (2014), we determine   for the Hercules mode as a function of radius using our Gaia DR3 data set at  = 6 • .Our results are shown as a solid black line with the grey region marking the 1 error based on 1,000 monte carlo simulations varying the velocities of the stars within their Gaia errors.The data points with errors are the measurements fromAntoja et al. (2014).