© ESO 2019

1. Introduction

The local velocity distribution of stars near the Sun has long been known to exhibit clear substructures most likely caused by a combination of resonances with multiple non-axisymmetric patterns (e.g. Dehnen 1998; Famaey et al. 2005) and incomplete phase-mixing (Minchev et al. 2009; Gómez et al. 2012) linked to external perturbations. The second data release from the Gaia mission (Gaia Collaboration 2018) has now revealed this rich network of substructures in unprecedented detail, displaying multiple clearly defined ridges in velocity space (Katz et al. 2018; Ramos et al. 2018) and even vertical velocity disturbances (Antoja et al. 2018; Quillen et al. 2018a; Monari et al. 2018), which have been associated with the perturbation of the disc by the Sagittarius dwarf galaxy (e.g. Laporte et al. 2019) or the buckling of the Galactic bar (Khoperskov et al. 2019). Regarding in-plane motion, it is also interesting to consider the distribution of stars in the space of actions, which are adiabatic invariant integrals of the motion constituting the natural coordinate system for Galactic dynamics and perturbation theory. Trick et al. (2019) produced similar plots in various volumes around the Sun. For local stars (d <  200 pc), they revealed several prominent ridges in the radial action distribution, among which a double-peak at the lowest border of the local azimuthal action (angular momentum) distribution corresponding to the well-known Hercules moving group at low azimuthal velocities, and one at high angular momentum corresponding to an arch at high velocities “covering” the velocity ellipsoid from above at V ∼ 40 km s⁻¹, where V is the heliocentric tangential velocity.

The Hercules moving group, in particular, has long been suspected to be associated with the perturbation of the potential by the central bar of the Galaxy (Dehnen 1999a, 2000). If the Sun is located just outside the bar’s outer Lindblad resonance (OLR), where stars make two epicyclic oscillations while making one retrograde rotation in the frame of the bar, the Hercules moving group is naturally generated by the linear deformation of the unperturbed background phase-space distribution function (DF). However, this implies a rather fast pattern speed for the bar, of the order of Ω_b = 55 km s⁻¹ kpc⁻¹ (see e.g. Fux 2001; Chakrabarty 2007; Minchev et al. 2007, 2010; Quillen et al. 2011; Antoja et al. 2014; Fragkoudi et al. 2019), which was favoured by past independent measurements (Englmaier & Gerhard 1999; Fux 1999; Debattista et al. 2002; Bissantz et al. 2003). By solving the linearized Boltzmann equation in the presence of the simple quadrupole bar potential of Dehnen (2000), we were able to show in Monari et al. (2016, 2017a,b) that the Hercules moving group was indeed naturally formed outside of the bar’s OLR, and that its observed position in velocity space was varying as predicted by such a model. However, recent measurements of both the three-dimensional density of red clump giants (Wegg et al. 2015) and the gas kinematics in the inner Galaxy (Sormani et al. 2015) indicate that the pattern speed of the Galactic bar could in fact be significantly slower, as hinted by some older studies (Weiner & Sellwood 1999; Rodriguez-Fernandez & Combes 2008; Long et al. 2013). From dynamical modelling of the stellar kinematics in the inner Galaxy, Portail et al. (2017, hereafter P17) recently deduced a pattern speed of Ω_b = 39 km s⁻¹ kpc⁻¹. Pérez-Villegas et al. (2017) then showed, with orbit simulations in the potential of P17, that stars trapped at the co-rotation resonance of such a bar could also reproduce the Hercules position in local velocity space. Using an update of the Tremaine & Weinberg (1984) method, which is similar to that of Debattista et al. (2002), on proper motion data from a combination of multi-epoch data from the VVV survey and Gaia DR2, Sanders et al. (2019) found a pattern speed of Ω_b = 41 ± 3 km s⁻¹ kpc⁻¹. Simultaneously, Clarke et al. (2019) also derived line-of-sight integrated and distance-resolved maps of mean proper motions and dispersions from the VVV Infrared Astrometric Catalogue combined with data from Gaia DR2, and found excellent agreement with a bar pattern speed of Ω_b = 37.5 km s⁻¹ kpc⁻¹. Studying in detail all the dynamical effects of such a slowly rotating bar on the local velocity field is thus extremely timely.

In Monari et al. (2017c, hereafter M17), we developed an analytical method for capturing the behaviour of the stellar phase-space DF at a resonance (see also Binney 2018), where the linearisation of the collisionless Boltzmann equation yields a divergent solution (problem of small divisors). This is a fundamental difference between a model where the Sun is located just outside of the bar’s OLR and one where it is outside of the bar’s corotation (CR) as in P17. While the Hercules moving group is located outside of the trapping region in the former case, and can be treated through linearisation of the Boltzmann equation (Monari et al. 2016, 2017a,b), it is precisely located in the resonant trapping region for a slow bar like that of P17. However, with the M17 method, the deformation of velocity space induced by Dehnen’s bar potential with a slow pattern speed was found to be rather minor. Here, by using the actual bar potential of P17 instead, and by computing the perturbed DF in the resonant regions, we confirm that the Hercules moving group can indeed be reproduced (Pérez-Villegas et al. 2017), and that many of the most prominent features in local action space can in principle be associated with the resonances of a slow bar model with Ω_b = 39 km s⁻¹ kpc⁻¹.

In Sect. 2, we present the Galactic potential of P17 and our extraction of its Fourier modes. We then summarize in Sect. 3 the M17 method for treating the behaviour of the DF at resonances, and in Sect. 4 we apply this method to the main resonances of the m = 2, 3, 4, 6 Fourier modes of the P17 potential. We then validate our analytical treatment with numerical orbit integrations in Sect. 5, where we also compare the model to Gaia DR2 data. We conclude in Sect. 6.

2. The bar potential

We make use of the 2D Galaxy and bar potential Φ corresponding to the fit obtained by P17 in a range of radii between R = 0 and R = 12 kpc. As is evident from Fig. 1, this potential is very different from the pure m = 2 quadrupole that had been assumed in M17. In this potential model, the Sun is located at R = 8.2 kpc and at an inclined azimuth of 28° with respect to the long axis of the bar. We measure the azimuth ϕ from the long axis of the bar so that the Sun is at ϕ = −28°, i.e. in the direction opposite to the rotation of the Galaxy (and of the bar). The bar rotates at Ω_b = 39 km s⁻¹ kpc⁻¹, and the CR is located at R ∼ 6 kpc.

Fig. 1. Bar potential, i.e. the difference between the total P17 potential used in this work Φ and its axisymmetric m = 0 part Φ0. The Sun in this model is placed at (x, y) = (8.2 kpc, 0). The bar is rotating clockwise. Its long axis (blue negative contours of Φ − Φ0) is inclined at an angle of 28° from the line connecting the Sun to the centre of the Galaxy.

The potential Φ(R, ϕ) is defined on a grid, and the potential at any point in the Galaxy is obtained by spline interpolation. We also extract the Fourier components of the potential and in particular we use the m = 0 mode as the axisymmetric background potential, and the m = 2, 3, 4, and 6 modes in Sect. 3. In Fig. 2, we display these four Fourier modes. We note that the m = 3 mode should be zero for a model truly in equilibrium in the corotating frame, which is not exactly the case of the P17 model. We refer hereafter to the amplitude of the mth mode of the potential as Φ_m. In this way the background potential becomes Φ₀(R).

Fig. 2. Amplitude of the Fourier modes of the bar potential A(R), as a function of the Galactocentric radius R. In green m = 2, in blue m = 4, in purple m = 6, and in pink m = 3.

Using the background potential Φ₀, we can define the Galactic circular frequency Ω(R) and epicyclic frequency κ(R) in the usual way (see e.g. Binney & Tremaine 2008), and the circular velocity curve v_c(R) = RΩ(R).

For R ≳ 12 kpc, we assume an axisymmetric continuation of the form

(1)

i.e. a continuation corresponding to a flat circular velocity curve, and setting all non-axisymmetric modes Φ_{m >  0} = 0 in these regions. However, we do not set the boundary R_Φ precisely at 12 kpc, but slightly more inside, i.e. R_Φ = 10 kpc; in this way we find that the transition is smoother, and at this radius the main m >  0 components already have negligible amplitude.

3. Distribution function in resonant regions: action-angle formalism

Let (J_R, J_ϕ) be the radial and azimuthal actions and (θ_R, θ_ϕ) the canonically conjugated radial and azimuthal angles, defined in the background axisymmetric potential Φ₀. A star’s actions and angles (hereafter AA) are combinations of the star’s positions and velocities, and they are particularly convenient phase-space coordinates for several reasons (see Binney & Tremaine 2008). In particular, the equilibrium axisymmetric background DF can be written purely as a function of the actions from the Jeans theorem, and they are the most convenient coordinates for perturbation theory. In simple words, the actions identify a star’s orbit in phase-space, whilst the angles denote the phase of the star on that particular orbit. The larger the radial action J_R is, the more energetic its radial excursions are, and the more eccentric the orbit is. The azimuthal action J_ϕ represents the vertical component of the angular momentum L_z. Here, we approximate the true values of the AA using the epicyclic approximation (Binney & Tremaine 2008; M17). In further work we will extend the present model to more realistic AA variables obtained through a combination of the Torus Machinery (e.g. Binney & McMillan 2016) to go from AA variables to positions and velocities, and the “Stäckel fudge” (e.g. Binney 2012; Sanders & Binney 2016) for the reverse transformation. For this reason, the results obtained in this analytic approach are still not fully quantitative, and we confirm them in Sect. 5 with backward orbit integrations not making use of the AA variables in computing the response of the DF to the bar perturbation. However, our analytic approach offers a way to understand the physical mechanisms at play in the backward simulations.

Using the AA and the Galactic frequencies Ω and κ, we can define an unperturbed DF for the Galactic disc, i.e. a DF that would not change in time, because of the Jeans theorem, if there was no bar perturbation, which we denote f₀(J_R, J_ϕ). As in M17, we chose the quasi-isothermal DF defined by Binney & McMillan (2011), with a scale length of 2 kpc, a velocity dispersion scale-length of 10 kpc, and a local velocity dispersion σ_R(R₀) = 45 km s⁻¹ (slightly hotter than in M17).

The response of the unperturbed DF to the bar potential is strongest at the resonances, which happen at the locations of phase space where

(2)

where and are the orbital frequencies, and simply become, in the epicyclic approximation, ω_R = κ(R_g) and ω_ϕ = Ω(R_g)+[κ(R_g)/J_ϕ]J_R, where R_g(J_ϕ) is the guiding radius (e.g. Dehnen 1999b).

In Eq. (2), we consider the main resonances of the m = 2, 3, 4, 6 Fourier modes of the potential. For the m = 2 mode, we shall take care of the (l, m) = (0, 2) resonance, namely the CR, and the (l, m) = (1, 2), namely the OLR. For the other modes, we do not consider the CR in the analytic approach, in order to avoid overlap of resonances, but we treat the (l, m) = (1, 3),(1, 4), (1, 6) resonances of the m = 3, 4, 6 Fourier modes of the potential.

To study the response of f₀ near a resonance, we have to define, in each of the five resonant zones considered here, a new set of AA variables. We have to go through a first canonical transformation of coordinates, from the old AA (J_R, J_ϕ, θ_R, θ_ϕ) to new “fast” and “slow” AA (J_f, J_s, θ_f, θ_s). The canonical transformation is Weinberg (1994; M17)

(3)

where (l, m) = [(0, 2),(1, 2),(1, 3),(1, 4),(1, 6)] for the CR and 2:1 OLR of the m = 2 mode, 3:1 resonance of the m = 3 mode, 4:1 of the m = 4 mode, and 6:1 of the m = 6 mode. The angle θ_s is called the slow angle because it evolves slowly near a resonance, as is evident from the definition of the frequencies and that of the resonance in Eq. (2). On the other hand, θ_f changes more rapidly and is called the fast angle. Physically, the slow angle θ_s represents the azimuth of the apocentre of the orbit in the reference frame where the unperturbed orbit would be close to the resonance. Therefore, θ_s represents the angle of precession of the orbit. At this point, expanding the potential in AA coordinates, it is possible to show that near the resonances J_f is almost constant along a star’s orbit, and that we can average the Hamiltonian over the fast angles. For each value of J_f, the motion in θ_s then becomes approximately like a pendulum, with J_s(t) the momentum of the pendulum. The steps to show this are explained in detail in M17. We can then define a pendulum energy parameter k depending on the phase-space coordinates that determines whether the pendulum is librating (i.e. the θ_s angle oscillates back and forth between a maximum and a minimum without ever covering the whole [0, 2π] range), or circulating (i.e. θ_s has a motion that covers the whole range of angles). Stars at a position of phase-space where the θ_s pendulum is librating are called “trapped to the resonances”. In Fig. 3, we display in local velocity space, at the position of the Sun, the five zones of trapping associated with the five resonances which we consider in the model of P17.

Fig. 3. Resonant zones in local velocity space (at the Sun’s position) for different Fourier modes and resonances of the P17 bar, within the epicyclic approximation used in our analytical model (see Sects. 3 and 4). The resonant zones are defined as the regions where the energy parameter k <  1 (see M17 and Sect. 3), i.e. the pendulum associated with each resonance is librating. Shown are the corotation and 2:1 (OLR) resonance for the m = 2 mode of the potential (green and red, respectively); the 3:1 resonance of the m = 3 mode (pink); the 4:1 resonance of the m = 4 mode (blue); and the 6:1 resonance of the m = 6 mode (purple).

For a pendulum, it is of course possible to make a new canonical transformation defining the actual AA of the pendulum itself: (θ_p, J_p). The way these depend on the phase-space coordinates changes according to whether the orbit is trapped or not. The pendulum AA also define the motion in J_s(θ_p, J_p) near the resonances. We can therefore rewrite the unperturbed DF as f₀(J_f, J_s(θ_p, J_p)).

In M17 (see also Binney 2018), the perturbed DF close to the resonances was then defined as the original DF phase-mixed over the angles θ_p

(4)

where the average is taken over the angle θ_p. Outside the zone of resonance, the DF is instead described as

(5)

where in this case the average represents the average J_s of the circulating motion. We note that these recipes used to obtain the perturbed DF are, in principle, different for every resonance also outside of the regions of trapping, while here we want to simultaneously describe five resonances of the Galactic bar. However, the value obtained outside the zone of resonance for f is very similar in the case of the five resonances because the J_s oscillations amplitudes decrease quickly going away from the resonant zone. It is therefore sufficient to take, outside the resonant regions of the phase space, the mean of the five DFs obtained for the five resonances.

4. Analytical results

As in M17, hereafter we concentrate on computing the perturbed DF at a few configuration space points, scanning through velocity space in order to define, at each phase-space point and for each resonance, the pendulum equation of motion for θ_s and its associated AA coordinates (J_p, θ_p), allowing us to compute the perturbed DF as described in the previous section and as detailed in M17.

We start by plotting the perturbed DF as a function of the coordinates of velocity space at the position of the Sun in the model of P17. In the left panel of Fig. 4, we plot the DF in (u, v) coordinates, namely the Galactocentric radial velocity in the direction of the Galactic centre, and the peculiar tangential velocity with respect to the circular velocity at the Sun. As can be seen in this figure, all five main resonances are present as deformations of the velocity distribution, which would be otherwise smooth. In particular, the effect of the m = 2 CR is to form a Hercules-like moving group (Pérez-Villegas et al. 2017) with a peak at (u, v) ≃ (−25, −50) km s⁻¹, and two asymmetric extensions at positive and negative u. At the same time, the 6:1 resonance forms a horn-like structure at larger v and positive u, while the OLR forms a structure similar to the arch present at high azimuthal velocities in Gaia data (although not as strong as in the data).

Fig. 4. Perturbed DF for stars located at the Sun’s position in the model, obtained with the analytical method described in Sect. 3. Left panel: velocity space (u, v). Right panel: perturbed DF in axisymmetric epicyclic action space (Jϕ − Jcirc, JR), where Jcirc is the angular momentum of a circular orbit at R = 8.2 kpc. The contours include (from darkest to lightest) 50, 68, 80, 90, and 95% of the stars. The radial velocity dispersion of the unperturbed DF f0 at the Sun is σR = 45 km s−1. The coloured lines correspond to each other in the two panels, and are plotted to identify the same features in velocity and action spaces. The green and yellow lines correspond to the corotation, the red line to the 2:1 (OLR) resonance for the m = 2 mode, the pink line to the 3:1 resonance of the m = 3 mode, the blue line to the 4:1 resonance of the m = 4 mode, and the purple line to the 6:1 resonance of the m = 6 mode.

To show even more clearly all these structures, in the right panel of Fig. 4 we show the DF in the (J_R, J_ϕ) action space of the unperturbed axisymmetric background. In this case, the deformation due to the CR appears clearly as two stripes at J_ϕ − J_circ ≲ −500 km s⁻¹ kpc, corresponding to the deformation of the DF at positive and negative u in the Hercules region of velocity space. In action space, the OLR is the pointy structure at J_ϕ − J_circ ≃ 500 km s⁻¹ kpc. In total, the bar alone thus creates no less than six prominent ridges in velocity and action space, at high and at low angular momenta, contrary to common misconceptions that the bar can only create one or two ridges, which is the case only when the bar is a pure m = 2 mode (see e.g. Khanna et al. 2019).

In Fig. 5, we reproduce the same plots of the distribution in velocity and action spaces as a function of radius at the azimuth of the Sun. All structures shift with R towards lower angular momenta as expected. In particular the effect of the CR is to form a Hercules-like moving group (Pérez-Villegas et al. 2017) that shifts to lower v and becomes less populated as R increases, like the Hercules moving group does in the real data.

Fig. 5. As in Fig. 4, but for stars at different positions in the Galaxy (always aligned with the Sun and the Galactic centre). Top row: DF in local velocity space. Bottom row: DF in local action space. The panels from left to right represent stars at R = 7.4, 7.8, 8.2, 8.6, 9 kpc.

It is also interesting to study how the different resonant ridges tend to move in terms of their radial velocity u as a function of R. In Laporte et al. (2019) and Fragkoudi et al. (2019), this has been illustrated using the Gaia data to map the radial velocity as a function of position in the R − v space (where v refers here to the azimuthal velocity, even outside of the solar neighbourhood). In Fig. 6, we show a similar plot from our analytical model. Given the epicyclic approximation used in this approach, this should not serve as a direct quantitative comparison, but does capture the behaviour of the different resonances as a function of radius in the P17 model. In particular, there is a splitting between the ridges associated with the 6:1 and CR resonances at R ∼ 9.5 kpc, reminiscent of Fig. 13 in Laporte et al. (2019), as well as a ridge at high v corresponding to the 2:1 OLR also present in the data. Perhaps more subtly, the ridges associated with the 3:1 and 4:1 resonances are very close to each other; they create a thick ridge at R ≳ 9.5 kpc, which becomes thinner at smaller radii once the effect of the 3:1 resonance vanishes in this space, a feature which is also tentatively seen in the data. We note that, when the same P17 bar is given a pattern speed of Ω_b = 50 km s⁻¹ kpc⁻¹ (which is unphysical in the case of this particular bar model, because it would make the bar extend beyond its corotation radius), the ridge at high v disappears, and the two ridges corresponding here to the CR and 6:1 resonances do not display their characteristic splitting around 9 kpc.

Fig. 6. Mean −⟨u⟩ velocity (see colour bar) predicted by the analytical model in the (R, v) space, for stars aligned with the Sun and the Galactic centre (i.e. at the solar azimuth). The bins are ΔR = 0.4 kpc and Δv = 5 km s−1. The coloured lines correspond to the positions in this space of the OLR (red), 3:1 (pink), 4:1 (blue), 6:1 (purple), and CR (green). The velocity −⟨u⟩ was chosen to facilitate the qualitative visual comparison with the figures in Laporte et al. (2019) and Fragkoudi et al. (2019).

These plots, and in particular the one at R = 8.2 kpc, thus allow us to physically understand the possible bar-related origin of many of the action space ridges identified by Trick et al. (2019), among others. However, because of the epicyclic approximation used here, it is useful to switch to orbit integrations in order to directly compare the effects of the bar to the Gaia DR2 data.

5. Backward integrations and comparison to Gaia

In the analytical model we present above, we approximate the values of the AA variables using the epicyclic approximation. It should be kept in mind, however, that the position of the resonant zones in this analytic approach depends on the ability of the AA epicyclic approximation to represent the true AA variables, which is only the case close to the centre of the (u, v)-plane. While this drawback of the analytical method will be resolved in future work, we chose to confirm our analytical results here with another technique for obtaining the response of the DF to the bar perturbation, the backward integration method, used by Vauterin & Dejonghe (1997) and Dehnen (2000), among others. It has also recently been used by Hunt & Bovy (2018) to explore the possibility that the Hercules moving group could be caused by the m = 4 mode of a bar, with an amplitude similar to that found in N-body simulations. This method consists in integrating backwards the orbit of stars all at a certain point x of configuration space, but having different velocities on a grid in the (u, v) plane. The method is based on the conservation of the phase-space density. Let us imagine that, at some time t = t₁ in the past, the amplitude of the bar was null, and then it started to grow with time until the current configuration at t = 0. The value of the unperturbed DF f₀ at t = t₁ for the ith velocity grid point (i.e. an orbit at (x, v_i) at t = 0) was f₀(x_i, 1, v_i, 1), where x_i, 1 and v_i, 1 are the initial position and velocity of the orbit. Then, because of the conservation of phase-space density, its current value at the point x and the velocity v_i is

(6)

This method has to assume an integration time and a law for the growth of the bar. In our case we assume the default values used by Dehnen (2000), i.e. 4 bar rotation for the integration time, a growth of the bar regulated by the polynomial law in Eq. (4) of Dehnen (2000), and a bar growth time of 2 bar rotations. We assume in this case the full numerical potential developed by P17, with its extension outside of R = 12 kpc as in the other sections. Due to some heating induced by the bar growth, we chose a colder initial DF, with a local radial velocity dispersion of 35 km s⁻¹; we note that the location of the structures in velocity space are independent of the velocity dispersion.

The results of the backward integrations for the same points in configuration space as in Fig. 5 are presented in Fig. 7. From this figure it is clear that the structures with this method are more complex for two main reasons: the ongoing phase-mixing, which adds noise to the backward integrations, especially for long integration times (Fux 2001), and the effects of resonances not considered above, in particular the CR of the modes m >  2, which tend to boost in action space the second ridge (green line in Fig. 4) associated with the Hercules stream. This Hercules group is centred at (u, v)≃(−25, 50) km s⁻¹ for R = 8.2 kpc, and shifts towards lower (higher) v as R increases (decreases), as in the analytical prediction. In action space, we can also clearly identify the several ridges associated with the resonances described in Sect. 4, which shift in J_ϕ as a function of R.

Fig. 7. As in Fig. 5, but for the DF obtained using the backward integration technique described in Sect. 5. In this case the initial value of the radial velocity dispersion of the unperturbed DF f0 at the Sun is σR = 35 km s−1.

In Fig. 8 we compare the actual velocity distribution of stars in the solar neighbourhood from the exquisite Gaia DR2 data (Gaia Collaboration 2018, only stars with line-of-sight velocities in a spherical volume of 200 pc) with the distribution obtained from the backward integrations shown in the central panel, top row, of Fig. 7. For the kinematics of the solar neighbourhood, we show the (U, V) velocities, i.e. the Cartesian velocities in the rest frame of the Sun, with U positive towards the Galactic centre and V positive towards the Galactic rotation. The (u, v) velocities are related to the (U, V) velocities, once the velocity of the Sun with respect to the local circular speed (or the local standard of rest) (U_⊙, V_⊙) is known, which we discuss below. The model’s kinematics is hotter than that of the Gaia stars, and therefore for the latter we show the contour containing 99% of the stars as the most external contour, instead of 95% for the model. The comparison of the model with the data shows many similarities in the velocity distributions: the already cited Hercules moving group; the high V arch (even more visible in action space in Fig. 9), which appears in the 99% contour of the Gaia stars at V ∼ 40 km s⁻¹ (but is admittedly more marked in the data); two horn-like features, especially clear in the 68% contours, at (U, V) ∼ (50, −25) km s⁻¹, and at (U, V) ∼ (40, −40) km s⁻¹. These could be connected to the structures in the 50% and 68% contours in the models, at (u, v) ∼ (60, −30) km s⁻¹, and at (u, v) ∼ (50, −40) km s⁻¹. Finally, we see another large arc-like feature in the data, from (U, V) ∼ (−80, 0) km s⁻¹ to (U, V) ∼ (60, 0) km s⁻¹, which could be identified with the arc-like feature in the model between (u, v) ∼ (−70, 0) km s⁻¹ and (u, v) ∼ (70, 0) km s⁻¹.

Fig. 8. Left panel: observed local velocity UV plane for Gaia DR2 stars inside a spherical volume of radius 0.2 kpc, with the distances estimated by Bailer-Jones et al. (2018). The contours include (from darkest to lightest) 50, 68, 80, 90, and 99% of the stars. Right panel: (u, v) velocity distribution at the Sun obtained with the backward integrations discussed in Sect. 5 for the P17 model, i.e. the central panel (top row) of Fig. 7. In this case the contours include (from darkest to lightest) 50, 68, 80, 90, and 95% of the stars.

Fig. 9. Top row, left panel: Gaia (U, V) kinematics in the solar neighbourhood from Fig. 8; the contours contain 50, 68, 80, 90, and 99% of the stars. Top row, right panel: as in the left panel, but here the stars are plotted in the space of actions (Jϕ − Jcirc, JR); the actions are computed assuming the P17 background potential, the epicyclic approximation, and (U⊙, V⊙) = (11.1, 0) km s−1; Jcirc is the angular momentum of a circular orbit at the Sun. Bottom row, left panel: (u, v) distribution at the Sun obtained via backward integration from Fig. 8; the contours contain 50, 68, 80, 90, and 95% of the stars. Bottom row, right panel: as in the left panel, but here the stars are plotted in the space of actions (Jϕ − Jcirc, JR). In all the plots the coloured lines represent the same and corresponding positions in the (u, v), (Jϕ − Jcirc, JR), and (U, V) spaces. Shown is the effect of the CR (Hercules, green line), 6:1 (horn, purple line), and 4:1 (blue line) on both the model and data. The 2:1 resonance (red line) is also present in both, and most prominently seen in action space, but with a slight position mismatch between model and data, linked to the choice of V⊙ and the shape of the circular velocity curve of the axisymmetric background (see text). Also shown is the prominent deformation in the velocity plane of the data at (U, V) ∼ (−35, −15) km s−1, namely the Hyades moving group, which is not reproduced by our bar-only model, and which has long been suspected to be related to a spiral perturbation.

We analyse these structures in more detail in Fig. 9 by associating them with features in epicyclic action space. We identify the high V arch with a red line, which becomes particularly prominent in action space, the Hercules moving group with a green and a yellow line (for the negative and positive u ridges respectively), the arch close to V ∼ 0 with a blue line, and the most prominent horn at smaller V with a purple line.

To compare the coloured lines from the (U, V) plane of the data to the (u, v) plane of the model, and to convert both to action space, we have to assume a value for the Sun’s peculiar motion (U_⊙, V_⊙). In this case we use U_⊙ = 11 km s⁻¹ (Schönrich et al. 2010) and an unusually low value of V_⊙ = 0, which we discuss below. As we show in Sect. 4, the position of the structures in the (u, v) plane in the bottom left panel correspond to ridges in action space, already identified in Sect. 4. In the top right panel we plot the action distribution for the data¹, using the epicyclic approximation, the background potential of P17, and V_⊙ = 0. The distribution in this space presents ridges already observed by other authors (Trick et al. 2019), and we see how nicely they correspond to structures identified in velocity space and also found in the model. Interestingly, we noted that the m = 3 mode should be zero for a model truly in equilibrium in the corotating frame, and indeed this resonance does not stand out in the local data, even though it might have an effect at larger radii (see Fig. 6). On the other hand, the high angular momentum arch is prominent in action space both in the data and model, but its position is not perfectly reproduced. This is due to the choice of V_⊙ = 0, which is the value that allows the best superposition with all the structures in the data and in the model. However, for some individual structures, such as this high V arch, the superposition implies another value of V_⊙: the high V arch is better fit for the more usual value V_⊙ = 10 km s⁻¹. Another possibility is that the P17 barred model is displaying deficiencies in explaining the data in the regions of velocity space where other effects such as phase-wrapping due to external perturbers (e.g. Laporte et al. 2019) become important.

The value V_⊙ = 0 is unusual, and is at odds with most estimates of the same parameter in the literature (e.g. Dehnen & Binney 1998; Schönrich et al. 2010). However, this unusual value could simply reflect the fact that the circular velocity curve of the Milky Way is somewhat different from that of the P17 model. Modifying only sightly the circular velocity curve v_c(R) also modifies Ω(R) and κ(R), and consequently the relative position of the resonances in local velocity space (as well as the extension of the “gap” between the Hercules moving group and the main velocity ellipsoid, as was already noted by Dehnen 2000). We performed a simple fit to see how we could modify v_c(R) and V_⊙ to minimize the distance between the resonances and the position of the structures defined in the (U, V) plane of the data. In this toy model, we assumed for the background potential v_c(R) = v₀(R/R₀)^α, where R₀ = 8.2 kpc is the Sun’s radius and v₀ = 241 km s⁻¹ is the value of v_c(R₀) in the background model of P17. This simple fit shows, for V_⊙ = 8 km s⁻¹ and a locally decreasing rotation curve (α = −0.1), that we get the best agreement between the position of the structures in the (U, V) plane and the position of the resonances from the bar model of P17. However, all this could also depend on the precise value of the pattern speed, which might be slightly lower than assumed here (e.g. Clarke et al. 2019).

Finally, in Fig. 10, we present a further comparison of the model with the data, this time on a wider range of R, to show how the gap between the Hercules moving group and the main velocity mode changes with radius (already studied by Antoja et al. 2014; Monari et al. 2017b; Pérez-Villegas et al. 2017, before the Gaia DR2). We confirm the clear identification of the gap in the data and how it shifts with R in the top row of Fig. 10, following the vertical red line that we overplot on both the data and model. This line corresponds, at each R, to the velocity v of stars that have the same angular momentum (i.e. the same guiding radius) as those that are at the gap velocity at R = 8.2 kpc. In the bottom row, we show how the gap is also clear in the backward integration model, although less marked than in the data, in part due to a larger velocity dispersion, but also because of uncertainty on the unperturbed DF we assumed, and because the data must also be affected by other physical mechanisms. In particular, two bumps are especially clear in the data at 7.8 kpc, and are not reproduced by the model. Again these limitations appear precisely in the regions of velocity space where other effects such as the inner resonances of outer spirals, or phase-wrapping due to external perturbers, become important.

Fig. 10. Top row: observed V distribution for Gaia DR2 stars inside spherical volumes (of radius 0.2 kpc) centred at different Galactocentric radii and along the line connecting the Galactic centre to the Sun; the distances of the stars used here are those provided by Bailer-Jones et al. (2018). Bottom row: v distribution obtained from the backward integration already shown in Fig. 7 and marginalizing along the u-axis (assuming V⊙ = 0). The red lines correspond to the same velocities in both rows, and to the same angular momentum as that of the gap at R = 8.2 kpc.

6. Discussion and conclusion

In this work used the analytical method recently developed in Monari et al. (2017c, M17) to study the response of the Galactic disc DF to the large bar Galactic potential model developed by Portail et al. (2017, P17). This method is based on the use of action-angle (AA) coordinates and perturbation theory.

Extracting the Fourier modes of the bar from this Galactic potential model, we have shown that in particular the m = 2, m = 4, and m = 6 modes deform the disc DF in ways that resemble those shown by the second data release of the Gaia satellite around the Sun. The m = 2 modes CR and OLR can be tentatively associated with features like the Hercules moving group (downward concave in velocity space, as shown by the green and yellow lines on Figs. 4 and 9) and the high azimuthal velocity arch in the local velocity distribution of stars (albeit slightly less marked than in the data). Interestingly, the 6:1 resonance of the m = 6 mode corresponds to the so called horn feature of local velocity space.

We also performed backward integrations using the whole Galactic bar model (no extraction of modes), showing that the same features obtained with the analytical method are also present in velocity and action space (but slightly offset because of the epicyclic approximation used for computing actions), further complicated by the presence of other modes of the bar and by ongoing phase-mixing.

The features obtained with the models can be seen most clearly in action-space, both in the model and in the Gaia data. The comparison between the model and the heliocentric velocity features in the data favours very small peculiar tangential velocities of the Sun. However, we showed that a small change in the circular velocity curve in the model can result in agreement with the data and a peculiar velocity of the Sun V_⊙ = 8 km s⁻¹ that is much more in line with other estimates in the literature.

One caveat of using our method, as implemented in this work, is the use of the epicyclic approximation to estimate AA and orbital frequencies, which is valid only for orbits that deviate little from circularity. In the future we will update these results and extend them in order to describe higher eccentricity orbits using more precise AA and frequency estimation methods like the Stäckel fudge (Sanders & Binney 2016) and the torus machinery (Binney & McMillan 2016; Binney 2018). It would also be extremely important to extend perturbation theory methods like the one used here (or linear theory to study the perturbations of the DF far from the resonances; Monari et al. 2016) to the temporal evolution of the perturbations, and to be able to include non-equilibrium phase-mixing effects (Agobert et al., in prep.).

The purpose of this paper was to demonstrate which structures of velocity space (local and non-local) can be created by the resonances of a large Galactic bar alone. Future fits to larger volumes in configuration space, exploiting the full capabilities of the Gaia data, as well as future spectroscopic surveys such as WEAVE and 4MOST, will allow us to test whether such structures are present in the data at all radii and azimuths, and will allow us to test such large bar models with a low pattern speed. However, it should be kept in mind that other perturbations, mostly external perturbations and spiral arms, are also present, and that they cannot be ignored in a more thorough modelling. Importantly, their effects could explain some of the phase-space structures that the present bar-only model cannot. In particular, phase-wrapping due to external perturbations affect velocity space by creating ridges locally at km s⁻¹ (Minchev et al. 2009; Gómez et al. 2012; Antoja et al. 2018; Laporte et al. 2019), while the bar model studied here does not create any structures beyond the region of velocity space delineated by its CR at low v and OLR at high v. Some features identified in Ramos et al. (2018) and Laporte et al. (2019) are clearly outside of this zone. Moreover, Laporte et al. (2019) identified structures in vertical velocities in the R − v space, which cannot be explained by the in-plane resonances of the bar studied here. The same wave-like features have also been identified in the study of the mean vertical and radial velocities of the stars in Gaia DR1 and DR2 as a function of angular momentum (or guiding radius) by Schönrich & Dehnen (2018) and Friske & Schönrich (2019). These wave-like features are partly aligned in mean radial and vertical velocity, and also have an azimuthal dependence compatible with m = 2 and m = 4 symmetries, and are reminiscent of the bending modes formed by satellite interaction on a galactic disc. This is also the favoured explanation that Carrillo et al. (2019) propose for many of the observed radial and vertical motions observed in Gaia DR2. On the other hand, Quillen et al. (2018b) have shown that some ridges and arcs seen in the local velocity distribution could be caused by multiple spiral features with different pattern speeds. While our results remove, in principle, the need for an alternative explanation for some observed ridges that can be created by the bar resonances, it certainly does not exclude that other features are related to spiral arms. For instance, the Hercules moving group appears as a double-structure in the UV-plane of the data (see also Katz et al. 2018; Li & Shen 2019), and this splitting is not present in the uv-plane of the P17 bar-only model, thus asking for additional effects in order to split Hercules in two. One other obvious deficiency of the P17 bar model is its inability to explain the strong Hyades overdensity at (U, V) ∼ (−35, −15) km s⁻¹, which has long been suspected to be related to a spiral perturbation (e.g. Quillen & Minchev 2005; Pompéia et al. 2011; McMillan 2013). The shape of the median velocity field as a function of position uncovered by Katz et al. (2018) is also highly suggestive of a spiral perturbation (see also Siebert et al. 2012). Nevertheless, we conclude that the bar model studied here has some merits. Indeed, the most recent proper motion data within the inner Galaxy point towards a large bar with a low pattern speed (e.g. Clarke et al. 2019; Sanders et al. 2019). We showed here that the strength of the m = 2, 4, and 6 modes of the P17 large bar model produces resonances with a noticeable amplitude (see also Hunt & Bovy 2018), and that a low pattern speed puts those resonances at about the right positions in order to explain some of the most prominent ridges observed in local velocity and action spaces. Such a large bar model could thus serve as a basis model for future detailed study of other internal and external perturbations of the Galactic disc.

Acknowledgments

We thank the anonymous referee for a constructive report which greatly improved the article. BF and OG acknowledge the hospitality at the KITP, supported by the National Science Foundation under Grant No. NSF PHY-1748958, in the final stages of this work. BF acknowledges support from the ANR project ANR-18-CE31-0006. CW acknowledges funding from the European Union?s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 798384. 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.

References

All Figures

	Fig. 1. Bar potential, i.e. the difference between the total P17 potential used in this work Φ and its axisymmetric m = 0 part Φ0. The Sun in this model is placed at (x, y) = (8.2 kpc, 0). The bar is rotating clockwise. Its long axis (blue negative contours of Φ − Φ0) is inclined at an angle of 28° from the line connecting the Sun to the centre of the Galaxy.
In the text

	Fig. 2. Amplitude of the Fourier modes of the bar potential A(R), as a function of the Galactocentric radius R. In green m = 2, in blue m = 4, in purple m = 6, and in pink m = 3.
In the text

Fig. 3. Resonant zones in local velocity space (at the Sun’s position) for different Fourier modes and resonances of the P17 bar, within the epicyclic approximation used in our analytical model (see Sects. 3 and 4). The resonant zones are defined as the regions where the energy parameter k <  1 (see M17 and Sect. 3), i.e. the pendulum associated with each resonance is librating. Shown are the corotation and 2:1 (OLR) resonance for the m = 2 mode of the potential (green and red, respectively); the 3:1 resonance of the m = 3 mode (pink); the 4:1 resonance of the m = 4 mode (blue); and the 6:1 resonance of the m = 6 mode (purple).

In the text

Fig. 4. Perturbed DF for stars located at the Sun’s position in the model, obtained with the analytical method described in Sect. 3. Left panel: velocity space (u, v). Right panel: perturbed DF in axisymmetric epicyclic action space (Jϕ − Jcirc, JR), where Jcirc is the angular momentum of a circular orbit at R = 8.2 kpc. The contours include (from darkest to lightest) 50, 68, 80, 90, and 95% of the stars. The radial velocity dispersion of the unperturbed DF f0 at the Sun is σR = 45 km s−1. The coloured lines correspond to each other in the two panels, and are plotted to identify the same features in velocity and action spaces. The green and yellow lines correspond to the corotation, the red line to the 2:1 (OLR) resonance for the m = 2 mode, the pink line to the 3:1 resonance of the m = 3 mode, the blue line to the 4:1 resonance of the m = 4 mode, and the purple line to the 6:1 resonance of the m = 6 mode.

In the text

	Fig. 5. As in Fig. 4, but for stars at different positions in the Galaxy (always aligned with the Sun and the Galactic centre). Top row: DF in local velocity space. Bottom row: DF in local action space. The panels from left to right represent stars at R = 7.4, 7.8, 8.2, 8.6, 9 kpc.
In the text

Fig. 6. Mean −⟨u⟩ velocity (see colour bar) predicted by the analytical model in the (R, v) space, for stars aligned with the Sun and the Galactic centre (i.e. at the solar azimuth). The bins are ΔR = 0.4 kpc and Δv = 5 km s−1. The coloured lines correspond to the positions in this space of the OLR (red), 3:1 (pink), 4:1 (blue), 6:1 (purple), and CR (green). The velocity −⟨u⟩ was chosen to facilitate the qualitative visual comparison with the figures in Laporte et al. (2019) and Fragkoudi et al. (2019).

In the text

	Fig. 7. As in Fig. 5, but for the DF obtained using the backward integration technique described in Sect. 5. In this case the initial value of the radial velocity dispersion of the unperturbed DF f0 at the Sun is σR = 35 km s−1.
In the text

Fig. 8. Left panel: observed local velocity UV plane for Gaia DR2 stars inside a spherical volume of radius 0.2 kpc, with the distances estimated by Bailer-Jones et al. (2018). The contours include (from darkest to lightest) 50, 68, 80, 90, and 99% of the stars. Right panel: (u, v) velocity distribution at the Sun obtained with the backward integrations discussed in Sect. 5 for the P17 model, i.e. the central panel (top row) of Fig. 7. In this case the contours include (from darkest to lightest) 50, 68, 80, 90, and 95% of the stars.

In the text

Fig. 9. Top row, left panel: Gaia (U, V) kinematics in the solar neighbourhood from Fig. 8; the contours contain 50, 68, 80, 90, and 99% of the stars. Top row, right panel: as in the left panel, but here the stars are plotted in the space of actions (Jϕ − Jcirc, JR); the actions are computed assuming the P17 background potential, the epicyclic approximation, and (U⊙, V⊙) = (11.1, 0) km s−1; Jcirc is the angular momentum of a circular orbit at the Sun. Bottom row, left panel: (u, v) distribution at the Sun obtained via backward integration from Fig. 8; the contours contain 50, 68, 80, 90, and 95% of the stars. Bottom row, right panel: as in the left panel, but here the stars are plotted in the space of actions (Jϕ − Jcirc, JR). In all the plots the coloured lines represent the same and corresponding positions in the (u, v), (Jϕ − Jcirc, JR), and (U, V) spaces. Shown is the effect of the CR (Hercules, green line), 6:1 (horn, purple line), and 4:1 (blue line) on both the model and data. The 2:1 resonance (red line) is also present in both, and most prominently seen in action space, but with a slight position mismatch between model and data, linked to the choice of V⊙ and the shape of the circular velocity curve of the axisymmetric background (see text). Also shown is the prominent deformation in the velocity plane of the data at (U, V) ∼ (−35, −15) km s−1, namely the Hyades moving group, which is not reproduced by our bar-only model, and which has long been suspected to be related to a spiral perturbation.

In the text

Fig. 10. Top row: observed V distribution for Gaia DR2 stars inside spherical volumes (of radius 0.2 kpc) centred at different Galactocentric radii and along the line connecting the Galactic centre to the Sun; the distances of the stars used here are those provided by Bailer-Jones et al. (2018). Bottom row: v distribution obtained from the backward integration already shown in Fig. 7 and marginalizing along the u-axis (assuming V⊙ = 0). The red lines correspond to the same velocities in both rows, and to the same angular momentum as that of the gap at R = 8.2 kpc.

In the text