The signatures of the resonances of a large Galactic bar in local velocity space

TL;DR

This paper demonstrates that the Galactic bar alone can produce multiple prominent ridges in local velocity and action space, explaining complex structures observed by Gaia data through resonance analysis.

Contribution

The study introduces a method to analyze stellar phase-space at resonances and links six local action space ridges to specific bar resonances using a realistic Galactic bar model.

Findings

Six ridges in local action space are related to bar resonances.

Two ridges correspond to the corotation resonance, linked to the Hercules group.

Other ridges are associated with the outer Lindblad and higher-order resonances.

Abstract

The second data release of the Gaia mission has revealed a very rich structure in local velocity space. In terms of in-plane motions, this rich structure is also seen as multiple ridges in the actions of the axisymmetric background potential of the Galaxy. These ridges are probably related to a combination of effects from ongoing phase-mixing and resonances from the spiral arms and the bar. We have recently developed a method to capture the behaviour of the stellar phase-space distribution function at a resonance, by re-expressing it in terms of a new set of canonical actions and angles variables valid in the resonant region. Here, by properly treating the distribution function at resonances, and by using a realistic model for a slowly rotating large Galactic bar with pattern speed 39 km/s/kpc, we show that no less than six ridges in local action space can be related to resonances with the bar. Two of these at low angular momentum correspond to the corotation resonance, and can be associated to the Hercules moving group in local velocity space. Another one at high angular momentum corresponds to the outer Lindblad resonance, and can tentatively be associated to the velocity structure seen as an arch at high azimuthal velocities in Gaia data. The other ridges are associated to the 3:1, 4:1 and 6:1 resonances. The latter can be associated to the so-called 'horn' of the local velocity distribution. While it is clear that effects from spiral arms and incomplete phase-mixing related to external perturbations also play a role in shaping the complex kinematics revealed by Gaia data, the present work demonstrates that, contrary to common misconceptions, the bar alone can create multiple prominent ridges in velocity and action space.

Equations

$$\Phi (R,\varphi )={\Phi }_0(R_{\Phi })+v_{\mathrm{c}}(R_{\Phi }{)}^2\ln (R/R_{\Phi }),\mathrm{for}R\ge R_{\Phi },$$

$$l{\omega }_R(J_R,J_{\varphi })+m\left[{\omega }_{\varphi }(J_R,J_{\varphi })-{\Omega }_{\mathrm{b}}\right]=0.$$

$${\theta }_{\mathrm{s}}$$

$${\theta }_{\mathrm{f}}$$

$$f=⟨f_0(J_{\mathrm{f}},J_{\mathrm{s}}({\theta }_{\mathrm{p}},J_{\mathrm{p}}))⟩,$$

$$f=f_0(J_{\mathrm{f}},⟨J_{\mathrm{s}}({\theta }_{\mathrm{p}},J_{\mathrm{p}})⟩),$$

$$f(\mathbf{x},{\mathbf{v}}_i)=f_0({\mathbf{x}}_{i,1},{\mathbf{v}}_{i,1}).$$

Full text

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 of 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 et al., 2018) has now revealed this rich network of substructures with unprecedented details, 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., 2018b; Monari et al., 2018), which have been associated to 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). As far as in-plane motions are concerned, 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 such 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\sim 40\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{km}\leavevmode\nobreak\ \mathrm{s}^{-1}$, where $V$ is the heliocentric tangential velocity.

The Hercules moving group, in particular, has long been suspected to be associated to the perturbation of the potential by the central bar of the Galaxy (Dehnen, 1999b, 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). This however implies a rather fast pattern speed for the bar, of the order of $\Omega_{\mathrm{b}}=55\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{km}\leavevmode\nobreak\ \mathrm{s}^{-1}\leavevmode\nobreak\ \mathrm{kpc}^{-1}$ (see also, e.g., Fux, 2001; Chakrabarty, 2007; Minchev et al., 2007, 2010; Quillen et al., 2011; Antoja et al., 2014; Fragkoudi et al., 2019), which past independent measurements also favoured (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 could show in Monari et al. (2016, 2017b, 2017c) 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 be, in fact, 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 $\Omega_{\mathrm{b}}=39\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{km}\leavevmode\nobreak\ \mathrm{s}^{-1}\leavevmode\nobreak\ \mathrm{kpc}^{-1}$. 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. Recently, using an update of the Tremaine & Weinberg (1984) method – 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) also found a pattern speed of $\Omega_{\mathrm{b}}=41\pm 3\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{km}\leavevmode\nobreak\ \mathrm{s}^{-1}\leavevmode\nobreak\ \mathrm{kpc}^{-1}$. 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 $\Omega_{\mathrm{b}}=37.5\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{km}\leavevmode\nobreak\ \mathrm{s}^{-1}\leavevmode\nobreak\ \mathrm{kpc}^{-1}$. Studying in details all the dynamical effects of such a slowly rotating bar on the local velocity field is thus extremely timely.

In Monari et al. (2017a, hereafter M17), we developed an analytical method to capture 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 indeed a fundamental difference between a model in which 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. Indeed, 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, 2017b, 2017c), 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 to the resonances of a slow bar model with $\Omega_{\mathrm{b}}=39\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{km}\leavevmode\nobreak\ \mathrm{s}^{-1}\leavevmode\nobreak\ \mathrm{kpc}^{-1}$.

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 to treat the behaviour of the DF at resonances, and we apply in Sect. 4 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 $\Phi$ corresponding to the fit obtained by P17 in a range of radii between $R=0$ and $R=12\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{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\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{kpc}$ and at an azimuth inclined of $28^{\circ}$ with respect to the long axis of the bar. We measure the azimuth $\phi$ from the long axis of the bar, so that the Sun is at $\phi=-28^{\circ}$, i.e. in the direction opposite to the rotation of the Galaxy (and of the bar). The bar rotates at $\Omega_{\mathrm{b}}=39\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{km}\leavevmode\nobreak\ \mathrm{s}^{-1}\leavevmode\nobreak\ \mathrm{kpc}^{-1}$, and the CR is located at $R\sim 6$ kpc.

The potential $\Phi(R,\phi)$ 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 $m$-th mode of the potential as $\Phi_{m}$. In this way the background potential becomes $\Phi_{0}(R)$.

Using the background potential $\Phi_{0}$, we can define the Galactic circular frequency $\Omega(R)$ and epicyclic frequency $\kappa(R)$ in the usual way (see e.g. Binney & Tremaine, 2008), and the circular velocity curve $v_{\mathrm{c}}(R)=R\Omega(R)$.

For $R\gtrsim 12\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{kpc}$, we assume an axisymmetric continuation of the form

[TABLE]

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

3 The distribution function in resonant regions: action-angle formalism

Let $(J_{R},J_{\phi})$ be the radial and azimuthal actions and $(\theta_{R},\theta_{\phi})$ the canonically conjugated radial and azimuthal angles, defined in the background axisymmetric potential $\Phi_{0}$. A star’s actions and angles (from now on AA) are combinations of the star’s positions and velocities, and they are particularly convenient phase-space coordinates for several reasons, listed in Binney & Tremaine (2008). In particular, the equilibrium axisymmetric background DF can be written purely as a function of the actions from 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_{\phi}$ 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 modelling 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 of 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 will 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. Our analytic approach however offers a way to understand the physical mechanisms at play in the backward simulations.

Using the AA and the Galactic frequencies $\Omega$ and $\kappa$, 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_{0}(J_{R},J_{\phi})$. As in M17, we choose 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 $\sigma_{R}(R_{0})=45\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{km}\leavevmode\nobreak\ \mathrm{s}^{-1}$ (slightly hotter than in M17).

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

[TABLE]

where $\omega_{R}=\dot{\theta}_{R}$ and $\omega_{\phi}=\dot{\theta}_{\phi}$ are the orbital frequencies, and simply become, in the epicyclic approximation, $\omega_{R}=\kappa(R_{\mathrm{g}})$ and $\omega_{\phi}=\Omega(R_{\mathrm{g}})+[\kappa(R_{\mathrm{g}})/J_{\phi}]J_{R}$, where $R_{g}(J_{\phi})$ is the guiding radius (e.g., Dehnen, 1999a).

In Eq. (2), we will 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 will not consider the CR in the analytic approach, to avoid overlap of resonances, but we shall 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_{0}$ 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_{\phi},\theta_{R},\theta_{\phi})$ to new ‘fast’ and ‘slow’ AA $(J_{\mathrm{f}},J_{\mathrm{s}},\theta_{\mathrm{f}},\theta_{\mathrm{s}})$. The canonical transformation is (Weinberg, 1994, M17):

[TABLE]

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 respectively. The angle $\theta_{\mathrm{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, $\theta_{\mathrm{f}}$ changes more rapidly and is called the ‘fast angle’. Physically, the slow angle $\theta_{\mathrm{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, $\theta_{\mathrm{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_{\mathrm{f}}$ is almost constant along a star’s orbit, and that one can average the Hamiltonian over the fast angles. For each value of $J_{\mathrm{f}}$, the motion in $\theta_{\mathrm{s}}$ then becomes approximately a pendulum one, with $J_{\mathrm{s}}(t)$ the momentum of the pendulum. The steps to show this are explained in detail in M17. One can then define a pendulum energy parameter $k$ depending on the phase-space coordinates, that determines whether the pendulum is ‘librating’ (i.e. the $\theta_{\mathrm{s}}$ angle oscillates back and forth between a maximum and a minimum without ever covering the whole $[0,2\pi]$ range), or ‘circulating’ (i.e. $\theta_{\mathrm{s}}$ has a motion that covers the whole range of angles). Stars at a position of phase-space where the $\theta_{\mathrm{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 to the five resonances which we consider in the model of P17.

For a pendulum, it is of course possible to make a new canonical transformation defining the actual AA of the pendulum itself: $(\theta_{\mathrm{p}},J_{\mathrm{p}})$. The way these depend on the phase-space coordinates changes according to whether the orbit is trapped or not. The pendulum AA define also the motion in $J_{\mathrm{s}}(\theta_{\mathrm{p}},J_{\mathrm{p}})$ nearby the resonances. We can therefore rewrite the unperturbed DF as $f_{0}(J_{\mathrm{f}},J_{\mathrm{s}}(\theta_{\mathrm{p}},J_{\mathrm{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 $\theta_{\mathrm{p}}$

[TABLE]

where the average is done over the angle $\theta_{\mathrm{p}}$. Outside the zone of resonance, the DF is instead described as

[TABLE]

where, in this case, the average represents the average $J_{\mathrm{s}}$ of the circulating motion. Note that these recipes to obtain the perturbed DF are, in principle, different for every resonance also outside of the regions of trapping, while we want here to describe simultaneously 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_{\mathrm{s}}$ oscillations amplitudes decrease fast 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

{comment}

As in M17, we will concentrate hereafter 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 $\theta_{\mathrm{s}}$ and its associated AA coordinates $(J_{\mathrm{p}},\theta_{\mathrm{p}})$, allowing to compute the perturbed DF as described in the previous Sect. 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. On 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 center, and the peculiar tangential velocity with respect to the circular one at the Sun. As can clearly be seen on this figure, all five main resonances are present as deformations of the velocity distribution, that would be otherwise smooth. In particular the effect of the $m=2$ CR is to form an Hercules-like moving group (Pérez-Villegas et al., 2017) with a peak at $(u,v)\simeq(-25,-50)\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{km}\leavevmode\nobreak\ \mathrm{s}^{-1}$, 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).

To show even more clearly all these structures, on the right panel of Fig. 4, we show the DF in the $(J_{R},J_{\phi})$ action space of the unperturbed axisymmetric background. In this case, the deformation due to the CR appears clearly as two stripes at $J_{\phi}-J_{\text{circ}}\lesssim-500\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{km}\leavevmode\nobreak\ \mathrm{s}^{-1}\leavevmode\nobreak\ \mathrm{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_{\phi}-J_{\text{circ}}\simeq 500\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{km}\leavevmode\nobreak\ \mathrm{s}^{-1}\leavevmode\nobreak\ \mathrm{kpc}$. In total, the bar alone thus creates no less than six prominent ridges in velocity and action space, at high as well as 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 an 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.

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 from Gaia data as a mapping of 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 produce 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, one can note a splitting between the ridges associated to the 6:1 and CR resonances at $R\sim 9.5\,{\rm kpc}$, reminiscent of Fig. 13 of 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 to the 3:1 and 4:1 resonances are very close to each other, and create a thick ridge at $R\gtrsim 9.5\,{\rm 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 $\Omega_{\mathrm{b}}=50\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{km}\leavevmode\nobreak\ \mathrm{s}^{-1}\leavevmode\nobreak\ \mathrm{kpc}^{-1}$ (which is of course 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$ becomes inexistent, and the two ridges corresponding here to the CR and 6:1 resonances do not display their characteristic splitting around 9 kpc.

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, e.g., Trick et al. (2019). However, because of the epicyclic approximation used here, it will be useful to switch to orbit integrations in order to directly compare the effects of the bar to Gaia DR2 data.

5 Backward integrations and comparison to Gaia

{comment}

In our analytical model hereabove, we approximated the values of the AA variables using the epicyclic approximation. One should however keep in mind that the position of the resonant zones in this analytic approach will depend on the ability of the AA epicyclic approximation to represent the true AA variables, which is the case only close to the center of the $(u,v)$-plane. While this drawback of the analytical method will be cured in further work, we chose to confirm our analytical results here with another technique to obtain the response of the DF to the bar perturbation, the ‘backward integrations’ used for example by Vauterin & Dejonghe (1997) and Dehnen (2000). It was also recently 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 $\boldsymbol{x}$ of configuration space now, 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_{1}$ 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_{0}$ at $t=t_{1}$ for the $i$-th velocity grid point (i.e. an orbit at $(\boldsymbol{x},\boldsymbol{v}_{i})$ at $t=0)$ was $f_{0}(\boldsymbol{x}_{i,1},\boldsymbol{v}_{i,1})$, where $\boldsymbol{x}_{i,1}$ and $\boldsymbol{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 $\boldsymbol{x}$ and the velocity $\boldsymbol{v}_{i}$ is

[TABLE]

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 by P17, with its extension outside of $R=12\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{kpc}$ as in the other Sections. Due to some heating induced by the bar growth, we choose a colder initial DF, with a local radial velocity dispersion of $35\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{km}\leavevmode\nobreak\ \mathrm{s}^{-1}$, noting 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. 10. 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 on Fig. 4) associated to the Hercules stream. This Hercules group is centred at $(u,v)\simeq(-25,50)\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{km}\leavevmode\nobreak\ \mathrm{s}^{-1}$ for $R=8.2\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{kpc}$, and shifting towards lower (higher) $v$ as $R$ increases (decreases), as in the analytical prediction. In action space, one can also clearly identify the several ridges associated to the resonances described in Sect. 4, which are shifting in $J_{\phi}$ as a function of $R$.

In Fig. 11 we compare the actual velocity distribution of stars in the solar neighbourhood from the exquisite Gaia DR2 data (Gaia Collaboration et al., 2018, only stars with line-of-sight velocities in a spherical volume of 200 pc) with the one obtained from the backward integrations shown in the central panel, top row, of Fig. 10. 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 center 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_{\odot},V_{\odot})$ is known, which we discuss hereafter. The model’s kinematics is hotter than that of the Gaia stars, and therefore for the latter we show the contour containing 99 per cent of the stars as the most external contour, instead of 95 per cent 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 presence of the high $V$ arch (even more visible in action space in Fig. 12) which appears only in the 99 per cent contour of the Gaia stars at $V\sim 40\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{km}\leavevmode\nobreak\ \mathrm{s}^{-1}$ (but is admittedly more marked in the data), two ‘horn’-like features, especially clear in the 68 per cent contours, at $(U,V)\sim(50,-25)\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{km}\leavevmode\nobreak\ \mathrm{s}^{-1}$, and at $(U,V)\sim(40,-40)\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{km}\leavevmode\nobreak\ \mathrm{s}^{-1}$. These could be connected to the structures in the 50 and 68 per cent contours in the models, at $(u,v)\sim(60,-30)\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{km}\leavevmode\nobreak\ \mathrm{s}^{-1}$, and at $(u,v)\sim(50,-40)\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{km}\leavevmode\nobreak\ \mathrm{s}^{-1}$. Finally, we can notice the presence of another large arc-like feature in the data, from $(U,V)\sim(-80,0)\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{km}\leavevmode\nobreak\ \mathrm{s}^{-1}$ to $(U,V)\sim(60,0)\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{km}\leavevmode\nobreak\ \mathrm{s}^{-1}$ and which could be identified with the arc-like feature in the model between $(u,v)\sim(-70,0)\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{km}\leavevmode\nobreak\ \mathrm{s}^{-1}$ and $(u,v)\sim(70,0)\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{km}\leavevmode\nobreak\ \mathrm{s}^{-1}$.

We analyze in more detail these structures in Fig. 12, by associating them to features in epicyclic action space, using colored lines. 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\sim 0$ with a blue line, and the most prominent horn at smaller $V$ with a purple line.

To compare the colored lines from the $(U,V)$ plane of the data to the $(u,v)$ plane of the model, and convert both to action space, we have to assume a value for the Sun’s peculiar motion $(U_{\odot},V_{\odot})$. In this case we use $U_{\odot}=11\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{km}\leavevmode\nobreak\ \mathrm{s}^{-1}$ (Schönrich et al., 2010) and an unusually low value of $V_{\odot}=0$, which we will discuss below. As we have already seen 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 data111The action distribution for the data is slightly more spread in $J_{\phi}$, even at small $J_{R}$, than in the model. This is due to the fact that, while in the model we plot the action space corresponding to stars all at one point in configuration space, the stars in the data are distributed in a sphere of finite volume and have observational errors., using the epicyclic approximation, the background potential of P17, and $V_{\odot}=0$. The distribution in this space presents ridges, already observed by other authors (Trick et al., 2019), and we see how nicely these correspond to structures identified in velocity space and present also 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 is not standing 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_{\odot}=0$. The latter is indeed the value that allows the best superposition with all the structures in the data and in the model, but for some individual structures such as this high $V$ arch, the superposition implies another value of $V_{\odot}$: the high $V$ arch is better fit for the more usual value $V_{\odot}=10\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{km}\leavevmode\nobreak\ \mathrm{s}^{-1}$. 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_{\odot}=0$ is of course unusual, and 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_{\mathrm{c}}(R)$ also modifies $\Omega(R)$ and $\kappa(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_{\mathrm{c}}(R)$ and $V_{\odot}$ 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_{\mathrm{c}}(R)=v_{0}(R/R_{0})^{\alpha}$ where $R_{0}=8.2\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{kpc}$ is the Sun’s radius and $v_{0}=241\leavevmode\nobreak\ \mathrm{km}\leavevmode\nobreak\ \mathrm{s}^{-1}$ is the value of $v_{\mathrm{c}}(R_{0})$ in the background model of P17. This simple fit shows that, for $V_{\odot}=8\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{km}\leavevmode\nobreak\ \mathrm{s}^{-1}$ and a locally decreasing rotation curve ($\alpha=-0.1$), 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. 13, we present a further comparison of the model with the data, this time on a larger range of $R$, to show how the gap between the Hercules moving group and the main velocity mode changes with radius (something that was already studied by Antoja et al., 2014; Monari et al., 2017c; 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. 13, 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\leavevmode\nobreak\ \leavevmode\nobreak\ \mathrm{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 not reproduced by the model. Again such 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.

{comment}

6 Discussion and conclusion

In this work have used the analytical method recently developed in Monari et al. (2017a, 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 from this Galactic potential model the Fourier modes of the bar, 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 to features like the Hercules moving group (concave downwards 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 of the circular velocity curve in the model can result in agreement with the data and a peculiar velocity of the Sun $V_{\odot}=8\leavevmode\nobreak\ \mathrm{km}\leavevmode\nobreak\ \mathrm{s}^{-1}$ 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 to the possibility to describe higher eccentricity orbits using more precise AA and frequencies 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 the 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 (local and non-local) velocity space 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 to test whether such structures are present in the data at all radii and azimuths, allowing to test such large bar models with a low pattern speed. However, one should keep in mind that other perturbations, mostly external perturbations and spiral arms, are present too, 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 $\sqrt{u^{2}+v^{2}}\gtrsim 40\,\leavevmode\nobreak\ \mathrm{km}\leavevmode\nobreak\ \mathrm{s}^{-1}$ (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 have also 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. (2018a) 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 at explaining the strong Hyades overdensity at $(U,V)\sim(-35,-15)\leavevmode\nobreak\ \mathrm{km}\leavevmode\nobreak\ \mathrm{s}^{-1}$, 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.

{comment}

Acknowledgements.

We thank the anonymous referee for a constructive report which greatly improved the article. BF and OG acknowledge hospitality at the KITP, supported by the National Science Foundation under Grant No. NSF PHY-1748958, in the final stages of this work. 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.