Galactic Bar Resonances Inferred from Kinematically Hot Stars in Gaia EDR3

Using a numerical simulation of an isolated barred disc galaxy, we first demonstrate that the resonances of the inner bar structure induce more prominent features in the action space distribution for the kinematically hotter stars, which are less sensitive to the local perturbation, such as the transient spiral arms. Then, we analyse the action distribution for the kinematically hotter stars selected from the Gaia EDR3 data as the stars with higher values of radial and vertical actions. We find several resonance features, including two new features, in the angular momentum distribution similar to what are seen in our numerical simulations. We show that the bar pattern speeds of about $\Omega_{\rm bar}\sim34$~km~s$^{-1}$~kpc$^{-1}$ and 42~km~s$^{-1}$~kpc$^{-1}$ explain all these features equally well. The resonance features we find correspond to the inner 4:1, co-rotation, outer 4:1, outer Lindblad and outer 4:3 (co-rotation, outer 4:1, outer Lindblad, outer 4:3 and outer 1:1) resonances, when $\Omega_{\rm bar}\sim34$ (42) km~s$^{-1}$~kpc$^{-1}$ is assumed.


INTRODUCTION
The central few kpc of the Milky Way show a prominent bar structure (e.g.Bland-Hawthorn & Gerhard 2016).The solidly rotating bar components affect the radial and rotational velocity distribution of the Galactic disc stars, and the presence of groups of stars moving with particular radial and rotational velocities in the Solar neighbourhood can be attributed to the bar (Dehnen 2000).Dehnen (1999) suggested that the Hercules stream, which is a group of stars rotating slower and moving outward in the disc, is caused by the outer Lindblad resonance (OLR hereafter) of the bar being just inside of the Sun's orbital radius.If true, this allows us to derive the pattern speed of the bar (Dehnen 1999;Monari et al. 2017a;Fragkoudi et al. 2019), and it is found to be fast (e.g.53 km s −1 kpc −1 ; Dehnen 1999).
Gaia data release 2 (DR2, Gaia Collaboration et al. 2018a) revolutionised our view of the kinematic structure of stars not only in the Solar neighbourhood but also for several kpc across the Galactic disc (Gaia Collaboration et al. 2018b;Antoja et al. 2018;Kawata et al. 2018;Friske & Schönrich 2019).It is well complimented by ★ E-mail: d.kawata@ucl.ca.uk near-infrared photometric surveys, such as VISTA Variables in the Via Lactea (VVV; Minniti et al. 2010), and ground-based spectroscopic surveys, such as Bulge Radial Velocity Assay (BRAVA; Kunder et al. 2012), the Abundances and Radial velocity Galactic Origins Survey (ARGOS; Freeman et al. 2013;Ness et al. 2013) and the Apache Point Observatory Galactic Evolution Experiment (APOGEE; Majewski et al. 2017), which revealed the detailed stellar structure and line-of-sight velocities within the Galactic bar itself.These observations were directly compared with theoretical models (e.g.Shen et al. 2010;Portail et al. 2015), which suggest a slower pattern speed of the bar than is found when assuming the OLR to be just inside of the Solar radius.Recently, both the gas dynamics (Sormani et al. 2015) and stellar dynamics from Gaia DR2, combined with VVV and APOGEE data (Sanders et al. 2019;Bovy et al. 2019) are converging on a value for the bar pattern speed of around 40 km s −1 kpc −1 .
Interestingly, the pattern speed of 40 km s −1 kpc −1 can explain the Hercules stream with the outer 4:1 resonance (Hunt & Bovy 2018).However, for such a pattern speed there should be a clear OLR feature in the kinematics just outside of the Solar neighbourhood (Hunt et al. 2019;Trick et al. 2021).Pérez-Villegas et al. (2017); Monari et al. (2019) suggested that the co-rotation (CR hereafter) is attributed to the Hercules stream, which explains the other features better (see a comprehensive discussion in Trick et al. 2021).Most models that attempt to create the Hercules stream from the CR alone show that the effect on the local velocity distribution is significantly weaker than the observed data (e.g.Binney 2018;Hunt & Bovy 2018).This is reconcilable with the addition of spiral structure (Hunt et al. 2018(Hunt et al. , 2019)).Alternatively, while previous studies assume a bar which rotates with a fixed pattern speed, Chiba et al. (2021); Chiba & Schönrich (2021) demonstrated that the observed detailed features of the stellar phase space distribution are better explained by the scenario that the Galactic bar is slowing down, and their CR reproduces the kinematics of the Hercules stream without requiring spiral structure.Hence, the pattern speed of the bar is still in-debate and requires more data and theoretical modelling studies.
Confronting the data with various theoretical models has made us realise that the kinematic features observed in the Solar neighbourhood can be explained by several different pattern speeds of the Galactic bar (e.g.Trick et al. 2019Trick et al. , 2021;;Hunt et al. 2019), because the bar can induce similar features with different resonances, such as the co-rotation resonance (e.g.Pérez-Villegas et al. 2017;D'Onghia & L. Aguerri 2020;Chiba et al. 2021), the outer 4:1 resonance (simply 4:1R hereafter; Hunt & Bovy 2018), the OLR, or other higher order resonances (e.g.Monari et al. 2019;Asano et al. 2020).While kinematic structure induced by different order resonances will vary significantly over Galactic scales, it is non-trivial to identify causation with Solar neighborhood data, and strong features like the Hercules stream can be explained in multiple ways.In addition, it is further complicated by the fact that similar phase space features can also be caused by transient spiral arms (De Simone et al. 2004;Quillen et al. 2011;Hunt et al. 2018;Fujii et al. 2019), and the influence of dwarf galaxies such as Sagittarius (e.g.Laporte et al. 2019;Khanna et al. 2019), and in many cases these will overlap in kinematic space.
In particular, transient spiral arms can have a systematic effect on the velocity field around the spiral arms (e.g.Stephens & Boesgaard 2002;Grand et al. 2012;Baba et al. 2013;Kawata et al. 2014) and radial migration causes nontrivial effects on the orbital phase (e.g.Grand et al. 2015).Hence, the velocity field and phase angles of orbits derived from an asymmetric potential are heavily affected by the transient spiral arms.In addition, the phase angles currently suffer from significant selection effects (Hunt et al. 2020;Trick 2020).Hence, in this paper we focus on kinematically 'hotter' stars in the Galactic disc, which can be defined as stars with larger actions.We consider that kinematically hotter stars are less affected by the transient spiral arms.On the other hand, the Galactic bar resonances are relatively strong and long lived effects, and can affect the kinematically hot stars as well (Binney 2018).
In Section 2 we first demonstrate that this working assumption is valid, based on the results of an -body/SPH simulation.Then, Section 3 shows the action distribution of stars in the recently released Gaia early data release 3 (Gaia EDR3; Gaia Collaboration et al. 2021a).Using the results from -body/SPH simulation as a prior, we discuss which observed action space features of stars correspond to which resonances of the Galactic bar.Summary and discussion of this study are presented in Section 4.

BAR RESONANCES FEATURES IN N-BODY/SPH SIMULATION
We utilise the -body/SPH simulation of a Milky Way-like galaxy presented in Baba & Kawata (2020) and Baba et al. (2021) to study the features expected to arise from the Galactic bar resonances.This simulation is an isolated disc galaxy, consisting of gas and stellar discs and a classical bulge, evolved self-consistently within a rigid dark matter halo.It includes gas radiative cooling, far-ultraviolet heating, star formation and stellar feedback (Saitoh et al. 2008;Baba et al. 2017).The gas and stellar particle masses are about 9.1 × 10 3 M and 3 × 10 3 M , respectively, and the softening length is set to be 10 pc.
We use a snapshot at  = 7 Gyr of this simulation, where there is a clear bar whose pattern speed is around Ω bar ∼ 40 km s −1 .The bar also has an X-shaped boxy inner bar/bulge, and there are several transient spiral arms (Baba & Kawata 2020).The bar pattern speed, Ω bar , is measured by the change of the phase angle of  = 2 Fourier mode, and the time evolution of Ω bar is shown in Fig. 3 of Baba et al. (2021).Interestingly, Ω bar fluctuates in the time scale of about 100 Myr, likely due to the interaction with the spiral arm features (Wu et al. 2016;Hilmi et al. 2020).We have selected the star particles with galactocentric radius 3 <  < 18 kpc to avoid analysing too many particles in the inner region and with height || < 0.5 kpc from the mid-plane of the disc.The relatively broad vertical region is selected to include kinematically hot disc particles.We select the wide radial range to cover many different resonances.Since we focus on the actions and orbital frequencies only, we use all the particles irrespective of their azimuthal angle position.
The number of particles are peaked around 8 kpc, because the density profile of the disc falls exponentially in the outer disc, and the area of the disc becomes smaller at smaller radii.To compensate for the decrease in the number of particles available for our analysis in the inner and outer disc, we weight the contribution of the particles from the inner and outer disc to make the weighted number of our particle sample at different radii constant.To compute the weight for each star, we first count the number of particles,  p ( i ), in the 64 radial bins within our sample radial range of 3 <  < 18 kpc, and compute the weight at the centre of each bin by  i = max( p )/ p ( i ), where max( p ) is the number of particles in the bin containing the maximum number of particles.Then, we compute the weight for each particle, depending on their radius by a linear interpolation of the weights at the centre of the radial bins.Although this has a negligible effect on our -body/SPH simulation results, we find that this is important for the observational data we analyse in the next section.
We compute the actions and orbital frequencies of the selected star particles using AGAMA (Vasiliev 2019) under the approximated gravitational potential of the -body/SPH simulation snapshot evaluated by AGAMA itself.In this paper, we focus on the radial action,  R , vertical action,  z , and azimuthal action, which is angular momentum,  z , and the radial frequency, Ω R , vertical frequency, Ω z , and azimuthal frequency, Ω  .
Fig. 1 shows the distribution of  z - R for the selected particles in our simulations.The actions are normalised by the circular velocity of  circ = 197 km s −1 at 8 kpc of the simulation.As shown with the Gaia DR2 data (e.g.Trick et al. 2019Trick et al. , 2021;;Hunt et al. 2019), our simulation also shows several strong ridge features.These features are considered to be caused by the resonances of the bar and transient spiral arms (e.g.Hunt et al. 2019;Trick et al. 2021), as discussed in Section 1.The resonances of the bar with pattern speed Ω bar are defined by the condition of where  and  are integer values (e.g.Binney & Tremaine 2008).
Fig. 2 demonstrates that the two ridges associated to the 4:1R and OLR are more prominent for the stars with the higher  R .As suggested from the cosmological simulations of the barred galaxies in Fragkoudi et al. (2020), the OLR ridge is most prominent.The upper panel of Fig. 2 shows the distribution of  z for the stars with 0.07 <  R ( z,0 ) < 0.15.The distributions of  z are computed with scikit-learn Kernel Density Estimation (KDE; Pedregosa et al. 2011) with a Epanechnikov kernel with the kernel size of 0.03.The two peaks are prominent and co-located with the 4:1R and OLR highlighted with the orange and red bands, respectively.On the other hand, the lower panel shows the  z distribution for kinematically colder, 0.01 <  R ( z,0 ) < 0.02, stars.The distribution is much flatter, because many features are overlapping with each other.Hence, we consider that the strong resonance features, such as the 4:1R and OLR, are easier to identify in the kinematically hot stars.
The mechanism causing the ridge features around the resonances is still in-debate (e.g.Trick et al. 2021).The resonance scattering (e.g.Lynden-Bell & Kalnajs 1972; Sellwood 2010) could be responsible for the strong features in the higher  R .The orbital trapping (e.g.Monari et al. 2017b;Binney 2018;Chiba et al. 2021;Chiba & Schönrich 2021) could also be responsible for the feature along the resonances, and the features in higher  R could be due to a higher number of stars being trapped in the resonance, and a higher  R tail being more clear.The velocity fields of stars are also expected to be affected by these mechanisms.However, as discussed in Section 1, we consider that the transient spiral arms can wash out these velocity features.In fact, although it is not shown here, we find that the velocity fields of our -body/SPH simulation are not similar to what are shown in test particle simulations without transient spiral arms.The velocity fields around the spiral arms show the systematic motions due to the transient spiral arms (e.g.Kawata et al. 2014;Grand et al. 2016) to be the dominant effects on the velocity field rather than the bar resonances.Hence, we do not look at the velocity fields or the orbital phase angles, but focus on the action distribution only.Selecting particles with higher action also helps reduce the effect of the transient arms (e.g.Solway et al. 2012).
To make sure that these high  R particles are influenced by the bar resonance, we analyse the orbit of the particles around CR, 4:1R and OLR.Fig. 3 shows the orbits of eight randomly selected particles around  R ( z,0 ) = 0.1 and around CR, 4:1R and OLR.The orbits are drawn by connecting the position of these particles in the bar rotation frame at the previous outputs.Typical orbits in these resonances are seen for these particles.There are some contaminants from 3:1 and 5:1-like orbits found in 4:1R.However, there are particles with clear 4:1 orbit.Hence, we think that these particles are affected by the bar resonances.
The rest of the vertical bands highlighted with blue, cyan, green and grey in Fig. 2 correspond to the i4:1R, CR, 4:3R and 1:1R, respectively.There is no obvious peak around these resonances, except subtle peaks at i4:1R and 1:1R in the upper panel of Fig. 2. To further focus on the kinematically hotter stars, we have selected the particles in the upper panel of Fig. 2, i.e. particles with 0.07 <  R ( z,0 ) < 0.15, and analysed the distribution of  z as a function of  z , which is shown in Fig. 4. It is interesting to see that there are ridge features toward higher  z around the 4:1R and OLR, which are highlighted with the orange and red bands in the upper panel.This is because more stars are in these resonances and the high  z tail becomes conspicuous at lower  z (see also Trick et al. 2021).In other words, kinematically hotter stars again show a clearer signal of the stellar number distribution around the resonances.
Hence, we have further selected high  z star particles, and the upper panel of Fig. 4 shows the  z KDE distribution of the star particles with 0.005 <  z ( z,0 ) < 0.05 (229,099 particles), which are within the region highlighted with the pink shaded region in the main panel.The coloured vertical bands again highlight different resonances.The  z distribution of higher  z star particles show clear peaks around the 4:1 (orange) and OLR (red).There are also some small peaks around i4:1R (blue), 4:3R (green) and 1:1R (grey), though the 4:3R and 1:1R are more tentative.Interestingly, at the CR, we find a small dip or no particular feature in the number of particles.It looks that the CR is unstable for the particles to stay, perhaps because it is where radial migration is efficient (e.g.Sellwood & Binney 2002) and many higher order resonances overlap, and/or its overlap with the transient spiral arms (Wu et al. 2016;Hilmi et al. 2020) may help to release particles from the CR (Baba et al. in prep.).Investigating mechanisms causing these features is not the aim of this paper.Rather we will use these features seen for kinematically hot star particles, to identify the resonance features in the real Galaxy in the next section.

BAR RESONANCES FEATURES IN GAIA EDR3
Similar to how we analyse our -body/SPH simulation data, we have selected stars in Gaia EDR3 and analyse their actions and or- bital frequencies.We first select the stars in Gaia EDR3 which have radial velocities from Gaia RVS (Cropper et al. 2018, Seabroke et al. in prep.).We then apply quality cuts, selecting stars with renormalised unit weight error, RUWE < 1.4, and /  > 4.0, where  and   are parallax and parallax uncertainty, respectively.We obtain the distance to the stars simply with the inverse of the reported parallax, after applying the zero-point correction suggested by Lindegren et al. (2021) using the python code provided by the Gaia collaboration 1 .We assume a distance to the Galactic centre from the Sun,  0 = 8.178 kpc (Gravity Collaboration et al. 2019) and the Sun's vertical height from the mid-plane,  0 = 20.8pc (Bennett & Bovy 2019).We assume the Sun's rotation speed to be  = 248.5 km s −1 and  , = 8.5 km s −1 , calculated from the combination of the assumed  0 and the proper motion measurement of Sgr A * (Reid & Brunthaler 2020).We also use the Sun's peculiar motion in the radial direction,  , = −12.9km s −1 (positive toward the outer disc) and  −  circ ( 0 ) = 12.32 km s −1 , where  circ ( 0 ) is the circular velocity at  0 .After transforming the data to Galactocentric cylindrical coordinates, we select stars with 4 <  < 12 kpc to minimise contamination from spuriously large distances, and || < 0.5 kpc.Using the same method as Hunt et al. (2020), we compute the actions and orbital frequencies of the selected stars using the actionAngleStaeckel (Binney 2012) function in galpy, assuming the MWPotential2014 potential, which is fit to various observational constraints (Bovy 2015).Note that we again normalise orbital frequencies and actions with Ω 0 =  circ ( 0 )/ 0 and  ,0 =  0  circ ( 0 ), respectively, as we did in the previous section.
The upper panel of Fig. 5 shows the distribution of  z and 1 https://gitlab.com/icc-ub/public/gaiadr3_zeropoint R for our selected Gaia EDR3 stars.As known from the previous studies with Gaia DR2, it is striking to see the various ridge features in this action space.We note that the action space distribution of Gaia EDR3 with the RVS data are similar to what is seen in Gaia DR2 with the Bayesian distances derived by Schönrich et al. (2019).Still, thanks to the superb astrometric accuracy of Gaia EDR3, we find tentative new ridge features, as discussed later in more detail.
The upper panel of Fig. 5 shows much finer structures than the result of our -body/SPH simulation in Fig. 1, because the Gaia data are tracing the phase space distribution of stars with much finer resolution, especially in the local volume.However, it also means that the sample is dominated by the stars near the Sun.The lower panel of Fig. 5 shows a similar image to the upper panel, but stars with | −  0 | < 0.2 kpc are excluded.We can see a clear parabola feature of an excluded zone centred at ( z ,  R ) = (1.0z,0 , 1.0 z,0 ), and some of the strong features seen in the upper panel, e.g. a feature extending from ( z ,  R ) ∼ (0.7 z,0 , 0.05 z,0 ) to (0.6 z,0 , 0.1 z,0 ), disappear.This demonstrates that these disappeared features are purely due to the dominance of stars close to the Sun.
To mitigate this effect, as we did in the previous section, we weight the contribution of stars to this distribution in the action space depending on the Galactocentric radius of the stars, so that the weighted number of stars at different radii becomes constant.The weight for each star is computed with the same method as described in Section 2, but using the 64 radial bins within 4 <  < 12 kpc, the different sample radial range in this section.Fig. 6 shows the  z and  R distribution of stars after weighted the stellar contribution to the distribution depending on their Galactocentric radius.Although we can not eliminate the entirety of the selection bias, we at least eliminate the spurious features caused by the overwhelming number of local stars.Hence, in this paper we show the radius weighted results.
As seen in our -body/SPH simulations, only a few strong features extend to high  R , i.e.  R > 0.03.As we did in the previous section, Fig. 7 displays the  z KDE distribution of stars with 0.03 <  R ( z,0 ) < 0.1 (upper panel) and 0.01 <  R ( z,0 ) < 0.2 (lower panel).Note that our selected Gaia EDR3 stars show relatively lower action values than our -body/SPH simulation, and therefore we select a lower  R range to pick up the kinematically hot stars.As shown in the previous section using the -body/SPH simulation, we find clearer peaks in the  z distribution for higher  R stars, while colder stars show many smaller peaks, like waves, as discussed in the previous papers with Gaia DR2 (e.g.Friske & Schönrich 2019;Hunt et al. 2019;Trick et al. 2019Trick et al. , 2021)).
Following the strategy of the previous section, we select the stars in the upper panel of Fig. 7, and analyse the distribution of  z and  z in Fig. 8.As seen in our -body/SPH simulation (Fig. 4), there are several vertically extended ridges for example at  z ∼  z,0 and  z ∼ 1.2 z,0 .The upper panel shows the  z KDE distribution, when we further restrict the stars with a relatively higher  z range of 0.005 <  z ( z,0 ) < 0.02, shaded in pink in the lower panel.We do not select stars with  z > 0.02 z,0 , because these high  z stars are dominated by stars with low angular momentum of  z 0.5 z,0 , which are likely thick disc stars, and which overwhelm the  z distribution, making it difficult to identify smaller peaks.As a result, the  z distribution of our selected stars show several weak, but clear peaks.Thanks to Gaia EDR3 where astrometry and radial velocity are precisely measured for a large number of stars, even after this strict selection, there are 163,838 stars which contribute to this distribution.
Using the prior from our -body/SPH simulation, we consider that these features in high  R and  z stars are due to the bar resonances.Selecting stars with higher actions is analogous to selecting older stars with a higher velocity dispersion, because the actions are correlated with the age of stars in the Milky Way (e.g.Beane et al. 2018;Ciucă et al. 2021).As demonstrated with the -body/SPH simulation, we consider that these relatively old (but not as old as the thick disc stars) stars are good tracers to identify the resonances caused by the Galactic bar.A remaining question is which features correspond to which resonances.
After trying different pattern speeds of the bar and comparing the resonant position with the features in the  z distributions.We find two pattern speeds which equally well explain these features.The first one is a pattern speed of the Galactic bar of Ω bar ∼ 1.16Ω 0 .Note that exact location of the resonance is sensitive to the shape of the Galactic potential.The pattern speed which we show here is only a rough fit by eye under our assumed potential shape, i.e.MWPotential2014 from galpy (Bovy 2015).The pattern speed value does not mean to be quantitatively accurate, but should only provide a rough estimate.The lines and bands highlighted with blue, cyan, orange, red and green in Figs. 6, 7 and 8 correspond to the location of i4:1R, CR, 4:1R, OLR, 4:3R, respectively, when we assume Ω bar ∼ 1.16Ω 0 .Here, we again use the orbital frequencies of Ω  and Ω R to identify these resonances, as done in the previous section.
In the upper panel of Fig. 8 the three weak peaks at  z ∼  z,0 , 1.2 z,0 and 1.4 z,0 nicely aligned with the 4:1R, OLR and 4:3R.A large and broad peak around  z ∼ 0.4 z,0 could potentially be explained by i4:1R.However, as mentioned above, this feature could just be the dominance of the thick disc stars in the inner disc.Interestingly, the CR corresponds to a dip around  z ∼ 0.75 z,0 , which is more prominent in the upper panel of Fig. 7.This is consistent with our -body/SPH simulation, which also shows a subtle dip rather than peak of the  z distribution at the CR.Hence, we regard this as a consistent result with our -body/SPH simulation expectation.In this case, there is no resonance to explain a peak around  z ∼ 0.8 z,0 in top panel of Figs. 7 and 8.However, it can be considered that this peak appears because this is next to the dip of the CR.
If this is the true pattern speed of the Galactic bar, the peak features corresponding to the i4:1R and the 4:3R are newly identified features in the action space, to our knowledge.We believe that the latter one corresponds to the ridge feature found in the radius and rotation velocity distribution in one of Gaia EDR3 performance verification papers of Gaia Collaboration et al. (2021b), which they call AC Newridge1.Hence, the peak at  z ∼ 1.4 in the upper panel of Fig. 8 is a confident detection, and this could be the furthest bar resonance feature we have newly identified.As a result, if the bar pattern speed of about 1.16Ω 0 is close to the true bar pattern speed, we find the i4:1R, CR, 4:3R, OLR and 4:3R.It is quite remarkable to find i4:1R, which is expected to exist from our simulation in the previous section, and Gaia EDR3 may be revealing the resonance inside the Galactic bar.
However, we also find that Ω bar = 1.45Ω 0 shows an equally good match to the features in the action distributions.Figs. 9 and 10 show the same results as Figs.7 and 8, but overlaid with the position of the CR, 4:1R, OLR, 4:3R and 1:1R, when a bar pattern speed of 1.45Ω 0 is assumed.In this case, we consider the peak at  z ∼ 0.8 z,0 in the top panels of Figs. 9 and 10 to be due to the 4:3R.Again, the CR corresponds to the dip around  z ∼ 0.6 z,0 .However, this dip is not as clear as the one which we associate with the CR, when assuming Ω bar = 1.16Ω 0 .With this pattern speed, the furthest resonance is associated to the 1:1R, which is also expected to be visible from the prediction of our -body/SPH simulation in the previous section.As a result, if the bar pattern speeds is 1.45Ω 0 , our identified features correspond to the CR, 4:1R, OLR, 4:3R and 1:1R.It is also remarkable to identify the 1:1R.

SUMMARY AND DISCUSSION
Using an -body/SPH simulation of an isolated barred disc galaxy, we demonstrate that the resonances of the bar induce more prominent features in the action space distribution for the kinematically hot star particles, i.e. particles with relatively high actions, than the kinematically colder star particles.This is because kinematically colder stars are more affected by the weaker but local nonaxisymmetric strucures, such as transient spiral arms, than the kinematically hotter stars.Using this as a working assumption, we analyse the action distribution for the kinematically hotter stars identified in the recently provided Gaia EDR3 data with the radial velocities from Gaia RVS.After computing the actions of these stars, we find several features in the angular momentum,  z , distribution for  kinematically hot, relatively high  R and high  z stars.Due to the improved astrometry in Gaia EDR3, we find new ridge features extending from around  z = 0.6 z,0 and  z = 1.5 z,0 at  R = 0, although the features are tentative and close to the edge of the data.
Assuming these features correspond to the bar resonances as seen in our -body/SPH simulation, we find that bar pattern speeds of Ω bar = 1.16Ω 0 and Ω bar = 1.45Ω 0 both explain all these features well.With our assumed  0 and  circ ( 0 ), these pattern speeds correspond to 33.6 km s −1 kpc −1 and 42 km s −1 kpc −1 , respectively.When we adopt Ω bar = 1.16Ω 0 , the features correspond to the i4:1R, CR, 4:1R, OLR and 4:3R.When we adopt Ω bar = 1.45Ω 0 ,      the features can be explained by the CR, 4:1R, OLR, 4:3R and 1:1R.To our knowledge, these many resonance features have never been revealed in the Milky Way before.This demonstrates the power of the Gaia EDR3 data.
Interestingly, in both cases, the CR is identified as a dip rather than peak in the  z distribution.We find a similar deficit of the particles at the CR in the  z distribution of the -body/SPH simulation, and hence we regard this dip (or no peak) as a consistent feature of the CR.The mechanism responsible for the deficit of stars at the CR is not clear, and it requires more theoretical studies with both a bar and transient spiral arms, which our -body/SPH simulation has, and which are expected to impact the stellar motion.
The bar pattern speed of 1.16Ω 0 is consistent with what was suggested by Pérez-Villegas et al. ( 2017); Monari et al. (2019) to explain the Hercules stream with the CR and also consistent with the current pattern speed suggested by Chiba et al. (2021); Chiba & Schönrich (2021) with a slowing bar.This is also roughly consistent with the recently measured pattern speed of Ω bar = 35.4± 0.9 km s −1 kpc −1 from dynamical modelling of the proper motion of the stars in the bar region (Clarke & Gerhard 2021).This slow bar pattern speed puts the OLR on the so-called "Hat" phase space velocity feature identified in the Solar neighbourhood radial and rotational velocity distribution (Hunt et al. 2017;Hunt & Bovy 2018).Although we consider that the local velocity fields can be disturbed by the transient spiral arms (Hunt et al. 2018) and are not necessarily reliable indicators of the resonance features, Trick et al. (2021) discussed that at the OLR the radial motion of the stars flips from outward motion (inside the OLR) to inward motion (outside the OLR), and the Hat has this characteristic flip expected at the OLR.
The other pattern speed of 1.45Ω 0 is consistent with what is inferred from the kinematics of the stars in the Galactic bar, which are converging to around 1.45Ω 0 as suggested by Sanders et al. (2019); Bovy et al. (2019).Trick et al. (2021) discussed that this pattern speed places the OLR on the Sirius moving group, where they do not see the outward to inward motion flip.Ramos et al. (2018) discuss that the Sirius moving group is unlikely to be induced by the resonance, because their kinematic feature is more consistent with the constant kinetic energy rather than the constant angular momentum.However, we consider that the velocity fields could be distorted by perturbers, such as transient spiral arms or the impact of Sagittarius, and we do not require these kinematic features to occur at the resonances in this paper.Hence, we also consider this to be an acceptable pattern speed.In this case, the strong Hat ridge feature is explained with the 4:3R.Our conclusion of the two potential bar pattern speeds whose OLR corresponds to Hat or Sirius is consistent with the conclusion from the phase-angle analysis in Trick (2020).
We note that as mentioned in Section 3, the computed values of actions depend on the shape of the Galactic potential in the radial range of the orbits of our sampled stars.The results of this paper are under the assumption of the Galactic potential of MWPoten-tial2014 in galpy.To note the dependence on this assumption, we briefly discuss the results when we adopt the different shapes of the Galactic potential.To this end, we compute the actions of our sampled stars using McMillan17 (McMillan 2017) and Ir-rgang13III (Model III of Irrgang et al. 2013) Galactic potentials in galpy.We then applied the same selection of kinematically hot stars as the top panel of Fig. 8, i.e. 0.03 <  r ( z,0 ) < 1.0 and 0.005 <  z ( z,0 ) < 0.02.The angular momentum distribution of these stars are displayed in Figs.11 and 12. Here, we normalise the actions and angular momentum using ( 0 ,  circ ( 0 )) = (8.21kpc, 233.1 km s −1 ) and (8.33 kpc, 239.7 km s −1 ) for McMil-lan17 and Irrgang13III potentials, respectively.The vertical bands indicating the resonances are computed with the bar pattern speeds of Ω bar = 35.5 km s −1 kpc −1 ∼ 1.25Ω 0 (McMillan17) and 39 km s −1 kpc −1 ∼ 1.36Ω 0 (Irrgang13III), which are chosen to match the OLR with the Hat local kinematic feature as did in Fig. 8.
These results indicate that the required pattern speed of the bar to locate the OLR to the Hat feature depends on the shape of the Galactic potential.Also, the intervals of the resonance locations in  z are sensitive to the assumed Galactic potential, and the locations of i4:1R, CR and 4:3R with respect to the dips and peaks of the  z distribution are different from Fig. 8 with MWPotential2014.Both McMillan17 and Irrgang13III place the CR around the peak of the  z distribution, next to (lower side of) the dip which coincides with the CR in MWPotential2014 in Fig. 8. Also, the locations of i4:1R and 4:3R are not aligned well with the peaks, compared to the results in Fig. 8.Although this is beyond the scope of this paper, this may indicate that MWPotential2014 is a preferable shape of Galactic potential, because of better matching of the resonances with the peaks and dips of the  z distribution.In other words, alignment of the resonances with the  z distribution of kinematically hot stars could constrain the shape of the Galactic potential, if the features identified in this papers are truly induced by the resonances of the Galactic bar, as expected from our N-body simulation.
Although we tried to correct the bias due to the overwhelming number of local stars by weighting the contribution of stars in the analysis depending on their Galactocentric radius, we should still be careful with the observational selection bias which could influence our conclusions.Also, we must always remind ourselves that -body/SPH simulations can still be far from a true representation of the real Milky Way.Our simulation formed the bar at about 5.5 Gyr before the snapshot we used for this study.Also, the pattern speed of the bar does not change significantly since its formation, because of the rigid dark matter halo we used.Hence, the disc star particles were subject to the same resonance locations for a long time.The effect on the action distribution of the disc stars could be significantly different, if the bar of the Milky Way formed recently, or if the pattern speed of the bar was slowing down (Chiba et al. 2021;Chiba & Schönrich 2021).Also, if the Galactic disc recently experienced bar-buckling (e.g.Khoperskov et al. 2019) and/or perturbations from the satellite galaxies, such as the Sagittarius dwarf (e.g.Laporte et al. 2019), their effects could be significant enough to erase the resonance features.Further comparison between the observational data and the theoretical models taking into account all these potential effects would be necessary to recover the true Galactic bar pattern speed confidently, and ultimately understand the nature of the Galactic bar and its impact on the Galactic disc evolution.To this end, obtaining precise proper motions inside the Galactic bar will be crucial (e.g.Baba & Kawata 2020).The upcoming nearinfrared astrometry mission, Japan Astrometry Satellite Mission for INfrared Exploration (JASMINE; Gouda & Jasmine Team 2020) 2 , will be invaluable in providing proper motion of stars between the Sun and the Galactic centre.

Figure 2 .
Figure 2. The KDE distribution of  z for the star particles with different ranges of  R , 0.07 <  R < 0.15 (upper) and 0.01 <  R < 0.02 (lower).The vertical bands highlighted with blue, cyan, orange, red, green and grey respectively indicate the range of i4:1R, CR, 4:1R, OLR, 4:3R and 1:1R in the  R range of each panel measured from Fig. 1.

Figure 3 .
Figure 3. Examples of orbits of the particles around CR (left), 4:1R (middle) and OLR (right) with high  R in the bar rotation frame, as the bar is highlighted with the grey horizontal bar.Particle's  z and  R are shown in each panel.

Figure 4 .
Figure 4.The distribution of  z and  z for the star particles of our  -body simulation which have 0.07 <  R < 0.15, corresponding to the top panel of Fig. 2. The upper panel shows the KDE distribution of  z for the star particles with 0.07 <  R < 0.15 and 0.005 <  z < 0.05, which are in the pink shaded region in the main panel.The vertical bands highlighted with blue, cyan, orange, red, green and grey respectively indicate the range of i4:1R, CR, 4:1R, OLR, 4:3R and 1:1R for 0.07 <  R < 0.15.

Figure 5 .
Figure 5. Upper panel: The distribution of the angular momentum,  z , and radial action,  R , for the selected stars in Gaia EDR3 without any radius weight.Lower Panel: Same as the upper panel, but also overplotting the stars within distance of 0.1 kpc with white dots.

Figure 6 .
Figure 6.The distribution of the angular momentum,  z , and radial action,  R , for the selected stars in Gaia EDR3 when weighting by radius (see the text for more detail).The blue, cyan, orange, red and green solid lines indicate the i4:1R, CR, 4:1R and OLR and 4:3R, respectively, when we assume Ω bar ∼ 1.16Ω 0 .

Figure 8 .
Figure 8.The distribution of  z and  z for the stars in Gaia EDR3 which have 0.03 <  R ( z,0 ) < 0.1, corresponding to the top panel of Fig. 7.The upper panel shows the KDE distribution of  z for the star particles with 0.03 <  r ( z,0 ) < 1.0 and 0.005 <  z ( z,0 ) < 0.02, which is highlighted in the pink shaded region in the main panel.The vertical bands highlighted with blue, cyan, orange, red and green respectively indicate the range of the i4:1R, CR, 4:1R, OLR and 4:3R for 0.03 <  R ( z,0 ) < 0.1, when we assume Ω bar ∼ 1.16Ω 0 .

Figure 10 .
Figure 10.The distribution of  z and  z for the stars in Gaia EDR3 which have 0.03 <  R ( z,0 ) < 0.1, corresponding to the top panel of Fig. 9.The upper panel shows the distribution of  z for the star particles with 0.03 <  r ( z,0 ) < 1.0 and 0.005 <  z ( z,0 ) < 0.02, which is highlighted in the pink shaded region in the main panel.The vertical bands highlighted with cyan, orange, red, green and grey respectively indicate the range of CR, 4:1R, OLR, 4:3R and 1:1R for 0.03 <  R ( z,0 ) < 0.1, when we assume Ω bar ∼ 1.45Ω 0 .