Kinematics with Gaia DR2: The Force of a Dwarf

TL;DR

This study uses Gaia DR2 data to analyze the 3D velocity distribution in the Milky Way disc, revealing complex vertical motions, radial velocity gradients, and potential influences from the Galactic bar or Sagittarius dwarf interactions.

Contribution

It provides new insights into the velocity structures and vertical perturbations in the Milky Way disc, highlighting the impact of recent interactions and the importance of accurate distance measurements.

Findings

Vertical motions are consistent with breathing and bending modes.

Radial velocity gradient varies with azimuth, indicating complex dynamics.

Recent perturbations, possibly from Sagittarius dwarf, influence disc kinematics.

Abstract

We use Gaia DR2 astrometric and line-of-sight velocity information combined with two sets of distances obtained with a Bayesian inference method to study the 3D velocity distribution in the Milky Way disc. We search for variations in all Galactocentric cylindrical velocity components ($V_{\phi}$, $V_R$ and $V_z$) with Galactic radius, azimuth, and distance from the disc mid-plane. We confirm recent work showing that bulk vertical motions in the $R\text{-}z$ plane are consistent with a combination of breathing and bending modes. In the $x\text{-}y$ plane, we show that, although the amplitudes change, the structure produced by these modes is mostly invariant as a function of distance from the plane. Comparing to two different Galactic disc models, we demonstrate that the observed patterns can drastically change in short time intervals, showing the complexity of understanding the origin of vertical perturbations. A strong radial $V_R$ gradient was identified in the inner disc, transitioning smoothly from $16$ km s$^{-1}$ kpc$^{-1}$ at an azimuth of $30^\circ<\phi<45^\circ$ ahead of the Sun-Galactic centre line, to $-16$ km s$^{-1}$ kpc$^{-1}$ at an azimuth of $-45^\circ<\phi<-30^\circ$ lagging the solar azimuth. We use a simulation with no significant recent mergers to show that exactly the opposite trend is expected from a barred potential, but overestimated distances can flip this trend to match the data. Alternatively, using an $N$-body simulation of the Sagittarius dwarf-Milky Way interaction, we demonstrate that a major recent perturbation is necessary to reproduce the observations. Such an impact may have strongly perturbed the existing bar or even triggered its formation in the last $1\text{-}2$ Gyr.

Figures

Equations

$$(U,V,W{)}_{\odot }=(11.1,12.24,7.25){\text{km s}}^{-1}$$

$$V_{\odot }+V_{\mathrm{LSR}}=4.74R_0{\mu }_{l_{\mathrm{Sgr}{\mathrm{A}}^{\ast }}},$$

$$\mathrm{RUWE}=\mathrm{UWE}/u_0(G,C)$$

$$\mathrm{UWE}=\sqrt{{\chi }^2/(N-5)}$$

$${\chi }^2=\text{astrometric\_chi2\_al}$$

$$N=\text{astrometric\_n\_good\_obs\_al},$$

$$G=\text{phot\_g\_mean\_mag},$$

$$C=\text{phot\_bp\_mean\_mag-phot\_rp\_mean\_mag}.$$

$$P(s|\varpi ,G_{RVS})\propto P(\varpi |s,{\sigma }_{\varpi })\times s^{2}P(\mathbf{r})P(M_{G_{RVS}}),$$

$$P(\mathbf{r})\propto {\mathrm{N}}_{\mathrm{a}}\exp \left(-\frac{R}{R_d^{\mathrm{a}}}-\frac{|z|}{z_d^{\mathrm{a}}}\right)+{\mathrm{N}}_{\mathrm{b}}\exp \left(-\frac{R}{R_d^{\mathrm{b}}}-\frac{|z|}{z_d^{\mathrm{b}}}\right)+{\mathrm{N}}_{\mathrm{c}}{r}^{-3.39}$$

$$\Delta V_{z,Breathing}=\left(V_{z,North}-V_{z,South}\right)$$

$$\Delta V_{z,Bending}=\frac{1}{2}\left(V_{z,North}+V_{z,South}\right)$$

Full text

††pubyear: 2018††pagerange: Kinematics with Gaia DR2: The Force of a Dwarf–Kinematics with Gaia DR2: The Force of a Dwarf

1 Introduction

It is now well established that the Milky Way disc contains non-axisymmetric structures, such as the Galactic bar and spiral arms. The effect of these asymmetries has been linked to substructure in the local stellar phase-space, first seen clearly as clumps in the velocity distribution (Dehnen, 1998) using data from the Hipparcos astrometric satellite (Perryman & ESA, 1997). Using only this small disc patch around the Sun ($d<200$ pc), it was already possible to establish the need for spiral (Quillen & Minchev, 2005; Pompéia et al., 2011) and bar structure (Dehnen, 2000; Minchev et al., 2007; Antoja et al., 2009).

Pre-Gaia, non-zero mean motions were found in the radial direction (e.g. Siebert et al., 2011; Williams et al., 2013) with origins usually attributed to internal perturbations due to the bar and spiral arms (e.g. Siebert et al., 2012; Monari et al., 2014). Non-axisymmetries were also found in the direction perpendicular to the Galactic disc (Widrow et al., 2012; Williams et al., 2013; Carlin et al., 2013) whose origins are still under debate and may result from both internal (e.g. Faure et al., 2014; Monari et al., 2015; Monari et al., 2016a) and external perturbations, such as the passage of a satellite galaxy or a dark matter subhalo (Gómez et al., 2013; Widrow et al., 2014; D’Onghia et al., 2016; Laporte et al., 2018).

The emerging Gaia era, where distances and proper motions are being obtained with unprecedented number, accuracy and precision, is changing our understanding of the velocity distribution in the Milky Way. Thanks to these data, we are now in a much better position to trace the origins of disc asymmetries. These features could be used as a dynamical diagnostic to model the history of interactions between the Milky Way and its satellite galaxies, as well as understand the importance of internal mechanisms.

Already with just early mission data from the first data release Gaia DR1 (Gaia Collaboration et al., 2016b) combined with data from the Radial Velocity Experiment (RAVE; Steinmetz et al., 2006) fifth data release (DR5; Kunder et al., 2017) and inferred distances by Astraatmadja & Bailer-Jones (2016) and McMillan et al. (2018), Carrillo et al. (2018, hereafter C18) studied the 3D velocity distribution in the extended solar neighbourhood, focusing on north-south differences and the observed vertical velocity asymmetries. In contrast to previous work showing a rarefaction-compression behaviour, in which the vertical velocity distribution has odd parity with respect to the Galactic plane and even parity in the density distribution, known as a breathing mode, C18 showed a vertical pattern with a combination of a breathing mode inside the solar radius ($R_{0}$) and a bending mode (even parity in $V_{z}$ with odd parity in the density distribution) outside $R_{0}$. This mode combination is physically intuitive as inwards of the solar radius the Galactic bar and spiral arms naturally induce breathing modes (Faure et al., 2014; Monari et al., 2015; Monari et al., 2016a), while outside $R_{0}$ bending modes are consistent with external perturbations such as a passing satellite galaxy or, on a different scale, dark matter subhaloes. C18 also showed that $V_{z}$ depends strongly on the adopted proper motions and that distance uncertainties can create artificial wave-like patterns.

The accuracy and precision achieved in the second Gaia data release (Gaia DR2) has brought a significant increase of scientific discoveries. A considerable effort mapping the kinematics of Gaia DR2 was made by Gaia Collaboration et al. (2018b) finding a variety of velocity substructures and seeing previously known ones with much more clarity. In comparison to the combination of breathing and bending mode observed by C18 with Gaia DR1, Gaia Collaboration et al. (2018b) showed that as the volume around the Sun extended to larger distances the bending mode outside $R_{0}$ reversed direction (as previously seen in the outer disc using LAMOST (Wang et al., 2018)), while a breathing mode inside $R_{0}$ was confirmed.

Gaia Collaboration et al. (2018b) and Ramos et al. (2018) showed that the solar neighborhood in-plane velocity distribution (known as the $U-V$ plane, with $U$ and $V$ the Galactocentric radial and tangential stellar velocity components) exhibits long arches separated by $\sim 20$ km s*-1* in $V$, in addition to the well known low-velocity moving groups (such as Hercules, Hyades, etc.). This confirmed predictions by Minchev et al. (2009), who used a semi-analytical model and test-particle simulations to show that phase wrapping following a perturbation from a recent merger could create such arches, linking the separation between them to an impact $\sim 2$ Gyr ago. Some of the arches found in Gaia DR2 $U-V$ plane could also be interpreted as the effect of spiral arm crossings (Quillen et al., 2018), transient spiral arms (Hunt et al., 2018), and the effect of a long bar (Monari et al., 2018).

Another relevant discovery using Gaia DR2 data was the phase-space ($z,V_{z}$) spiral found in the distribution of solar neighbourhood radial and azimuthal velocities by Antoja et al. (2018), interpreted as further evidence for ongoing phase mixing in the Milky Way disc, due to the Sagittarius dwarf galaxy (see e.g. Johnston et al., 2005; Law & Majewski, 2010; Laporte et al., 2018). This result reinforced earlier interpretations of observational signs in the Milky Way taken as evidence that the Galaxy was perturbed by external agents (see e.g. Minchev et al., 2009; Kazantzidis et al., 2009; Gómez et al., 2012; D’Onghia et al., 2016) and in particular Sagittarius (Quillen et al., 2009; Purcell et al., 2011; Gómez et al., 2013; de la Vega et al., 2015).

Although phase-mixing neglects self-gravity and structure dissipates with time, the patterns reproduced by Laporte et al. (2019) obtained pre-Gaia DR2 using a full $N$-body simulation of the impact of Sagittarius with the Milky Way are consistent with the timescale and structure seen by Antoja et al. (2018). Binney & Schönrich (2018) explained the phase-space spiral using toy-models of tracers defined by distribution functions reacting to a point mass perturber. Chequers et al. (2018) used $N$-body simulations of an isolated disc-bulge-halo system to investigate the continual generation of bending waves by a system of satellites or dark matter subhaloes. Their results suggest that the phase-space spiral may be due to long-lived waves of a continually perturbed disc rather than, or in addition to, single/recent satellite interactions. Similarly, Darling & Widrow (2018) used test-particles simulations and $N$-body models to show that disc bending perturbations can naturally create phase-space spirals. Bland-Hawthorn et al. (2019) used Gaia DR2 data complemented with spectroscopy from the GALAH survey (De Silva et al., 2015; Martell et al., 2017) to study the phase-space spiral in abundance and action spaces. On the other hand, Khoperskov et al. (2018) used a high-resolution $N$-body simulation to demonstrate that not only external perturbations can create structures qualitatively similar to those observed by Antoja et al. (2018). They showed that vertical oscillations driven by the buckling instability of the bar could also naturally create phase-space spirals. The different approaches used to explain the origins of the phase-space spiral suggests that this structure might result from a superposition of waves due to internal and external mechanisms and disc phase-wrapping.

Gaia DR2 has also made it possible to identify a number of ridges of negative slopes in the $V_{\phi}-R$ stellar distribution (Antoja et al., 2018; Kawata et al., 2018; Laporte et al., 2019; Fragkoudi et al., 2019; Khanna et al., 2019) separated roughly by $20-30$ km s*-1*, similarly to the arches found in the $U-V$ plane. This structure has been interpreted as phase-wrapping due to a perturbation with Sagittarius (Antoja et al., 2018), internal spiral arms (Kawata et al., 2018; Hunt et al., 2018), or the Galactic bar (Fragkoudi et al., 2019). Note that a ridge interpreted as the bar’s Outer Lindblad resonance in the $V_{\phi}-R$ plane was first identified in Gaia DR1 data by Monari et al. (2017).

In this work we follow up on the analysis of C18 using data from the second Gaia data release (Gaia DR2; Gaia Collaboration et al., 2018a) and two sets of distances obtained via Bayesian inference. The used priors consider the parallax and optical photometry from Gaia as well as multiband photometric information. These inferred distances, together with the outstanding improvement compared to the proper motions of Gaia DR1, allows for a deeper study of the radial and vertical motions by expanding the volume observed by Gaia Collaboration et al. (2018b). This could help us to further constrain the origins of such perturbations.

This paper is structured as follows: In Section 2, we introduce the data, our sample selection and present a comparison of our inferred distances and the direct use of the inverse parallax. A brief description of the Galactic disc models used in this work is given in Section 3. In Section 4, we present our kinematic analysis of the radial, azimuthal and vertical velocity distribution as a function of radius ($R$), azimuth ($\phi$), and height ($z$). Here we also briefly show the effects of systematic errors in the data and we further study the radial and vertical velocity distributions compared to simulations. Finally, Section 5 contains a summary of our results.

2 Data

2.1 Coordinate system

We used the sky positions, line-of-sight velocities, and proper motions from Gaia DR2, as well as two sets of distance estimates presented in Section 2.3, to compute the Galactocentric positions and velocities for stars in our sample. The method used to derive the velocities is described in detail by Johnson & Soderblom (1987). We computed the Galactocentric cylindrical velocity components ($V_{R}$, $V_{\phi}$ and $V_{z}$), following the coordinate transformation given in Appendix A of Williams et al. (2013).

We use the estimate of the peculiar velocity of the Sun obtained by Schönrich et al. (2010):

[TABLE]

and adopt $R_{0}=8.34$ kpc (Reid et al., 2014) for the solar Galactocentric distance. With these values and the proper motion of Sagittarius A*, $\mu_{l_{\mathrm{Sgr\,A^{*}}}}\,=6.379$ mas  yr*-1* (Reid & Brunthaler, 2004), we obtain:

[TABLE]

which yields a value of $V_{\mathrm{LSR}}\sim 240$ km s*-1* for the circular velocity of the local standard of rest (LSR). The selection of these parameters will affect the $V_{R}$ and $V_{\phi}$ values, but since $V_{z}$ is independent of $V_{\mathrm{LSR}}$ it will suffer no change.

2.2 Gaia DR2 sample selection

ESA’s mission Gaia (Gaia Collaboration et al., 2016a) is acquiring highly accurate parallaxes, proper motions, radial velocities, and astrophysical parameters for over a billion sources that will allows us to study the dynamics, the structure and the origins of our Galaxy. The current second data release Gaia DR2, contains median radial velocities for more than 7.2 million stars with a mean G magnitude between $\sim 4$ and 13 and an effective temperature in the range of $\sim 3550$ to 6900 K. The overall precision of the radial velocities (on the order of $\sim 200-300$ m s*-1* and $\sim 1.2$ km s*-1*, at the bright and faint ends, respectively) represents a considerable improvement to previous data. However, the biggest improvement comes from the parallaxes and proper motions with typical uncertainties of 0.04 mas and 0.06 mas yr*-1*, respectively. This will improve the velocity uncertainties up to a factor of $6$ compared to the results presented in C18.

In this work we select our Gaia DR2 sample based mainly on the recommended astrometric quality indicators such as the unit weight error (UWE) and the renormalised unit weight error (RUWE). The RUWE is computed following the recipe of Lindegren (2018):

[TABLE]

where ,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

$u_{0}(G,C)$ is a normalisation factor obtained from interpolating in $G$ and $C$ the provided table on the ESA Gaia DR2 known issues page.

We select stars with $\mathrm{RUWE}\leq 1.4$, which ensures astrometrically well-behaved sources. Further cuts include to select stars with VISIBILITY_PERIODS_USED $>8$, these cuts select stars with a better parallax determination and remove stars with parallaxes more vulnerable to errors. Finally, we select stars with an inferred distance uncertainty of $\sigma_{d}/d<0.2$.

2.3 Distance estimates

Over the last few years, it has become clear that using the inverse parallax as a distance estimator is not trivial (e.g. Brown et al., 1997; Arenou & Luri, 1999; Bailer-Jones, 2015; Astraatmadja & Bailer-Jones, 2016; Luri, Xavier et al., 2018). Parallaxes containing an error larger than $20\%$ result in biased distance estimates but it is also not recommended to select only stars with lower parallax uncertainties (although widely used), as this also creates biases. Bayesian inference therefore, has become a common technique to infer parallaxes from the observed data. In this Section, we briefly describe the methods used to estimate both sets of inferred distances and present a comparison of both samples with the inverse parallax of Gaia DR2.

McMillan (2018) estimated stellar distances, $s$, using the observed Gaia parallax, $\varpi$, and the radial velocity spectrometer magnitude, $G_{RVS}$, (calculated using the approximation given by Gaia Collaboration et al., 2018a). The Bayesian estimate is then given by:

[TABLE]

where $M_{G_{RVS}}$ is the absolute magnitude in the $G_{RVS}$ band. The prior $P(M_{G_{RVS}})$, which is an approximation to the distribution of $M_{G_{RVS}}$, was modelled with PARSEC isochrones (Marigo et al., 2017). The density model used to give $P(\mathbf{r})$ is the same as that used by McMillan et al. (2018), i.e., a three component model (thin disc, thick disc and halo, labelled a, b and c respectively),

[TABLE]

with $R_{d}^{\rm a}=2\,600\,\mathrm{pc}$, $z_{d}^{\rm a}=300\,\mathrm{pc}$, $R_{d}^{\rm b}=3\,600\,\mathrm{pc}$, $z_{d}^{\rm b}=900\,\mathrm{pc}$ and normalisations $\mathrm{N_{i}}$ (see McMillan et al., 2018 and McMillan, 2018 for further details).

While no colour information nor dust extinction was considered, the obtained distance nonetheless represents an important improvement over the naïve inverse parallax estimate.

Anders et al. (2019) derived distances using the StarHorse code (Santiago et al., 2016; Queiroz et al., 2018), which finds the posterior probability over a grid of stellar models, distances and extinctions, given a set of astrometric, spectroscopic, and photometric observations and a number of Galactic priors. For this, they combined parallaxes, adopting a fixed parallax zero-point shift of 0.05 mas (Zinn et al., 2018) and Gaia DR2 optical photometry ($G,G_{\rm BP},G_{\rm RP}$) with the photometric catalogues of 2MASS $JHK_{s}$, PANSTARRS and AllWISE $W1,W2$. Other priors used include an initial mass function as well as stellar density, metallicity, and age distributions for the main Galactic stellar components (see Queiroz et al., 2018 for further details). Due to the multi-band photometric data, the realistic Galactic density priors used and the treatment of dust extinction, an improvement in the distance precision is expected.

The left panel of Fig. 1 shows the distance uncertainty as a function of our two sets of inferred distances, as well as the distance obtained directly by computing the inverse parallax from Gaia DR2. Here it is clear, that StarHorse distances have smaller uncertainties at large ranges. The middle panel displays a comparison between our distance estimates. Both StarHorse and McMillan distances agree quite well up to $\approx 3$ kpc, beyond which we see a shift towards larger distances by McMillan compared to StarHorse. The shift may be due to dust extinction, which decreases apparent stellar brightness and if unaccounted for, leads to overestimate distances. The right panel of Fig. 1 displays a comparison of the relative distance uncertainties. Here, as previously expected, we see that the distance estimates from StarHorse have smaller uncertainties than the ones obtained by McMillan (2018).

As a result of the more precise estimates and the extinction treatment, we use the StarHorse distances derived by Anders et al. (2019) as the main distance set in this work. The data selection described in Section 2.2 yields a sample of 5,167,034 stars using the inferred distances of McMillan (2018) and 5,420,754 stars using StarHorse, which include additionally the cuts described in Anders et al. (2019): $\texttt{SH\_GAIAFLAG}="000"$ and $\texttt{SH\_OUTFLAG}="00000"$.

Another approach to obtain distance estimates was recently used by Lopez-Corredoira & Sylos Labini (2018), who used Lucy’s inversion method of the Fredholm integral equations of the first kind to estimate distances up to 20 kpc. This method, however, makes use of no Galaxy priors and does not consider extinction, which as discussed, overestimate distances.

3 Models

In Section 4.4 we compare our sample from Gaia DR2 to the results obtained with two different Galactic disc models with properties close to those of the Milky Way. The first one is selected from the suite of zoom-in simulations in the cosmological context presented by Martig et al. (2009, 2012) (Halo106). The radius was rescaled to 65% of its original value so as to match the Milky Way’s disc scale-length and rotation curve as described by Minchev et al. (2013), who used the same model for their chemo-dynamical model. We refer to this simulation hereafter as Model 1.

The second model is an $N$-body simulation presented by Laporte et al. (2018, 2019), which considered the interaction of a Sagittarius-like dSph with the Milky Way. This simulation has quantitatively demonstrated that Sagittarius can simultaneously account for the distribution of outer disc structures seen in the anticenter (Newberg et al., 2002; Price-Whelan et al., 2015; Bergemann et al., 2018; Sheffield et al., 2018) and most recently, the phase space spiral in ($z,V_{z}$) and ridges in the $V_{\phi}-R$ plane (Laporte et al., 2019). We use snapshots 690 and 648 of their L2 model and refer to it as Model 2. For Model 2 we rescale the distance to $80\%$.

4 6D phase space

4.1 Velocity maps in the $R-z$ plane

Figure 2 shows the median Galactocentric azimuthal ($V_{\phi}$), radial ($V_{R}$), and vertical ($V_{z}$) velocity fields obtained using the distances from StarHorse (top) and the ones of McMillan (2018) (bottom) in the $R$-$z$ plane. The StarHorse volume expands on the volume covered by the McMillan estimate because it has lower distance uncertainties at larger radii and it takes into account extinction, allowing for better estimates at lower latitude. The maps are shown in bins of (0.1 kpc)2, with a minimum of 50 stars per bin to avoid issues with low number statistics. Each map is computed by taking the mean of 100 realizations via Monte Carlo iterations. We do this by taking into account individual errors in distance, proper motion, line-of-sight velocity and drawing velocity values from a normal distribution centred on the originally estimated value. The black and white contours indicate significance larger than $2\sigma$ and $3\sigma$, respectively.

The $V_{\phi}$ maps (left column of Fig. 2) show a clearly defined wedge shape for contours representing $V_{\phi}\lesssim 220$  km s*-1* - as vertical distance increases at a given radius, we encounter slower rotating stars. This is related to the asymmetric drift effect, where stars with lower angular momentum reach outer radii at higher distance from the mid-plane.

The $V_{R}$ velocity fields (middle column of Fig. 2) present a more complex variation with $R$ and $z$. For example, as evident from the 2 and 3$\sigma$ contours (overlaid black and white lines) in the top panels, we can see a negative radial gradient in $V_{R}$ at $R\lesssim 10\,$kpc and $|z|\lesssim 2$ kpc, and a positive one at larger radii (similar results are seen in the bottom panels, where the negative $V_{z}$ inside $R\sim 5$ kpc is probably due to the larger uncertainties above the 2 and 3$\sigma$ levels). This is consistent with the results of C18 using Gaia DR1 and Gaia Collaboration et al. (2018b) using Gaia DR2. Radial variation of $V_{R}$ appears similar but shifted to larger radius at larger $|z|$, possibly due to the asymmetric drift effect shifting higher eccentric orbits to larger $R$, similar to the shift in resonances described by Mühlbauer & Dehnen (2003).

The $V_{z}$ maps (right column of Fig. 2) also show rich variations with radius and distance from the disc mid-plane. The arrows indicate bulk stellar motions at two different radii, showing a switch from a breathing at $R\approx 6.5$ kpc to a positive bending mode at $R\approx 13$ kpc. A similar combination of modes was previously found in Gaia DR1 (C18), where a negative bending was observed at $R\approx 8.5$ kpc.

Our data sample obtained with StarHorse distances (top row) is not only more precise (Section 2.3) but also covers a larger volume than the sample using McMillan (2018) distances. We will, therefore, continue our analysis with the former dataset.

4.2 Velocity maps in the $x-y$ plane

In Fig. 3 we present colour maps in the ($x\text{-}y$) plane for all velocity components, as indicated. The Sun is located at $x,y=(-8.34,0)$ kpc, the Galactic centre is at $x,y=(0,0)$ kpc, the Galaxy rotation is clockwise and there is no cut in $z$. This orientation is the same as in Gaia Collaboration et al. (2018b), Fig. 10.

The strong radial gradient in $V_{\phi}$ apparent in Fig. 2 is also seen in the left panel of Fig. 3. For $R\gtrsim 8$ kpc some azimuthal variations may be present on the order of 15  km s*-1*, while at $R\lesssim-4$ kpc strong double peaks are seen. Away from the Sun, the drop in $V_{\phi}$ may be due to selection effects, where the density of stars close to the plane decreases with distance from the Sun.

The middle panel of Fig. 3 presents the Galactocentric radial velocity, $V_{R}$. This shows similar structure to the top panels of Fig. 10 by Gaia Collaboration et al. (2018b), but the radial range covered is about twice as large. The larger volume covered in the inner disc is particularly useful to study the bar/bulge region. At $R\lesssim 4\,$kpc with a significance of at least $3\sigma$ (white contours), we see a clear separation between stars lagging and stars ahead of the solar azimuth with negative and positive $V_{R}$, respectively. A more detailed analysis of this structure compared to our models will be done in Section 4.4.

The vertical velocity map in the right panel of Fig. 3 shows an interesting variation close to the solar radius. Stars outside the solar radius exhibit upward stellar motion and downward motion inside. This pattern is consistent with the observed kinematic warp seen by Poggio et al. (2018). The ridge of positive $V_{z}$ is, however, not exactly centred at $\phi=0$, as expected from this feature. Poggio et al. (2018) argued that this could be due to the Sun not being on the line of nodes. This is possible if the deviation is large, but probably we are just observing a combination of modes or it is just simply that the line of nodes has an inclination with respect to the plane as a function of radius. Such a non-straight line of nodes is consistent with recent results presented by Romero-Gómez et al. (2018), who used data from Gaia DR2 to study the global shape of the warp.

4.3 Variations with Galactic azimuth

Figure 3 showed that there exists significant azimuthal variations in all three velocity components, suggesting that some structure in the $R\text{-}z$ plane velocity maps (Fig. 2) may have been washed out since we integrated over the entire azimuthal range.

To study variation in velocity structure with Galactic azimuth ($\phi$), we divide our sample in 15∘ slices in the range $-45^{\circ}<\phi<45^{\circ}$. The results are shown in Fig. 4. The left column shows the residual $\Delta V_{\phi}$ in (0.1 kpc)2 bins, computed by subtracting the mean of the median $V_{\phi}$ at all heights within the same $R$, from the median $V_{\phi}$ of each $R$-$z$ bin. This removes the Galactic rotation and the effect of the asymmetric drift seen in Fig. 2, allowing to better see structure at a given radius. In the inner disc, inside $R\sim 3$ kpc, we see a peak in $\Delta V_{\phi}$ at $0^{\circ}<\phi<30^{\circ}$ and a second one just inside the solar radius at $R\sim 7$ kpc. These peaks decrease at large and small azimuths, likely indicating the crossing of spiral arms.

The middle column of Fig. 4 also shows strong azimuthal dependence in $V_{R}$ bulk motions. The strongest variations are found in the inner disc and close to the mid-plane, in agreement with the maps in Fig. 3 (which is dominated by stars close to the plane). Inside $R_{0}$, we see symmetrical radial velocities across $z=0$. However, at $R>R_{0}$, $V_{R}$ is asymmetric across the mid-plane for $\phi<0^{\circ}$ but mostly symmetric at $\phi>0^{\circ}$. This is interesting because if the radial streaming motions were solely due to spiral structure, they should be symmetric across the Galactic plane. This suggests external perturbations, e.g., from the Sagittarius dwarf galaxy, are at play here. In the next Section we will study the azimuthal variations in more detail and compare them to our two models.

The $V_{z}$ maps displayed in the right column of Fig. 4 agree with the combined $R\text{-}z$ map showing a combination of breathing and bending modes (Fig. 2, right column). Nonetheless, some interesting azimuthal variations are seen at $-30^{\circ}<\phi<-15^{\circ}$ within $10<R<12$ kpc where negative $V_{z}$ is seen at $z\sim-1$ and positive $V_{z}$ above, corresponding to a breathing mode. This contrasts with the strong positive $V_{z}$ seen at $0^{\circ}<\phi<15^{\circ}$ above and below the plane and the negative bending mode seen outside the solar radius in the range $-45^{\circ}<\phi<-30^{\circ}$. Both the breathing and the negative bending mode, however, suffer from large uncertainties (see white and black contours) probably due to a low number of stars at these azimuths.

4.4 Radial velocity gradient

Here we investigate in more detail the variation of the $V_{R}$ radial gradient with Galactic azimuth and distance from the disc mid-plane. The $V_{R}$ map shown in Fig. 3 is not very informative as the colorbar range of $\pm 15$ km s*-1* was chosen in order to show structure in the outer disc, which resulted in saturation in the inner one. Hence, in Fig. 5 we show the median $V_{R}$ as a function of $R$ for $15^{\circ}$ azimuthal slices in the range $-60^{\circ}<\phi<60^{\circ}$, as indicated by the colorbar and divide our sample between $0<\,\left|\,z\,\right|<\,0.5$ kpc (top) and $0.5<\left|z\right|<1$ kpc (bottom). The left column shows the results obtained for stars above and the middle column for stars below the plane. The error bars correspond to the standard error of the median. The black solid lines represent the solar position and $V_{R}=0$ km s*-1*.

At distances closer to the disc mid-plane (top panels of Fig. 5), variations with azimuth are on the order of 10 km s*-1* near the solar radius, increasing outside $R_{0}$ to 30 km s*-1* and especially toward the bulge region, where values of up to $\pm 100$ km s*-1* are seen. In the inner disc, we observe a strong radial $V_{R}$ gradient, transitioning smoothly from $\approx+16$ km s*-1* kpc*-1* in a range between $30^{\circ}<\phi<45^{\circ}$ ahead of the Sun-Galactic centre line to $-16$ km s*-1* kpc*-1* in $-45^{\circ}<\phi<-30^{\circ}$ lagging the solar azimuth. In contrast, at larger distances from the disc plane (bottom panels), the radial velocity amplitude decreases to approximately half, while the trends remain similar.

The right column of Fig. 5 displays the $V_{R}$ residuals between the northern and southern Galactic hemispheres $\Delta V_{R}=\left|V_{R,\,North}\right|-\left|V_{R,\,South}\right|$. In the top right panel the residuals show at positive azimuths larger $V_{R}$ above the Galactic plane and at negative azimuths larger $V_{R}$ below. This type of shearing is not observed for stars farther away from the mid-plane, where larger $V_{R}$ are seen at $z>0$ for most azimuthal bins.

We are interested in understanding what may cause the azimuthal gradient in $V_{R}$. Therefore, we will focus our attention mostly to the inner disc. It is remarkable that the fanning of $V_{R}$ towards the Galactic centre shifts systematically from $\sim-50$ km s*-1* for $-30^{\circ}<\phi<-15^{\circ}$ (light blue curves in Fig. 5) to $\sim+50$ km s*-1* for $15^{\circ}<\phi<30^{\circ}$ (yellow curves in Fig. 5). It should be noted that wavy behaviour is present for each azimuth, with $V_{R}$ peaks shifting systematically with azimuth, as expected for a spiral pattern.

4.4.1 The effect of a steady-state bar

In order to understand the origins of such a pattern we now use the models described in Section 3. The left column of Fig. 6 shows face-on density maps for 5 different snapshots of Model 1, separated by $37.5$ Myr. During this time interval there are no significant mergers disturbing the disc. The dashed white circle represents the solar neighbourhood, the dotted white circle the bar region, and the colour lines indicate the $\phi$ angles as in Fig. 5. For this model, we fixed the Galactic bar at roughly the angle that the Sun lags the bar ($\phi=30\,^{\circ}$), while we let the galaxy rotate clockwise. The constant density contour levels are spaced by $10\%$ from 90 to $30\%$. The $V_{R}$ maps due to the interaction of the spiral arms with the bar are shown in the middle column of Fig. 6. The changes between different evolutionary stages are more evident outside the solar radius. This can be seen better in the right column, showing $V_{R}$ as a function of radius for different azimuths (as in Fig. 5). In the bar region (inside the red dotted line), however, it can be seen that, independently of time, $V_{R}$ exhibits a similar pattern with velocity decreasing with increasing $\phi$. This pattern remains stable for later evolutionary stages and is the expected velocity structure due to the bar rotation. These negative/positive streaming at the leading/trailing bar side can also be seen in the Auriga simulations (e.g. Grand et al., 2016, Figure 1), as well as test particle models (Monari et al., 2016b, Figure 4). It is therefore very interesting, that the data shown in Fig. 5 shows exactly opposite trends inside $R\sim 5$ kpc.

4.4.2 The effect of systematics

To see if this discrepancy can be explained by uncertainties in the data, we add random errors to one of the snapshots (t=12.90 Gyr). We use the data typical errors in proper motion $\sigma_{\mu}=0.07$ mas  yr*-1* and line-of-sight velocity $\sigma_{\mathrm{\,los}}=1.8$ km s*-1* together with the maximum accepted distance error of $\sigma_{d}=20\%$. The results can be seen in the second panel of Fig. 7. Compared to the snapshot without errors in the first panel, the inclusion of random errors is enough to affect the pattern inside the solar radius, but the systematic shift as a function of azimuth is not as clear as in the data. The third and fourth panels of Fig. 7 show the effect of possible systematics in our distances. For this, we multiply the distance by 0.8 and 1.2 and recompute $V_{R}$ (this is a simplified attempt to account for the problem as in reality the parallax shift might vary a lot with position and even colour of stars creating larger uncertainties). In the fourth panel, scaling the distance to $120\%$ increases the amplitude of the original pattern. In contrast, scaling the distance to $80\%$ (third panel) modifies the velocity amplitudes and inverts the trends for different azimuth systematically similar to the data.

To further study the effect of systematics, we create an idealized experiment of a Galaxy composed only of circular orbits by setting $V_{R}=\,0$ and $V_{\phi}=240$ km s*-1* for all the stars. We continue by introducing systematic errors in distance and proper motion and recomputing the velocities. The velocity maps shown in Fig. 8 and Fig. 9 were computed with the same cuts as described in Section 2.2, but do not include the $\sigma_{d}/d<0.2$ cut, as this would vary for each map. The top panels of Fig. 8 show the effects of including systematic errors of $20\%$. As can be seen, distance and proper motion systematics, both create a systematic shift as a function of $\phi$. While proper motion systematics are unlikely the main feature creating the $V_{R}$ dipole due to the $200\%$ systematic error needed to invert the pattern and the affected area not being in the bar region, a $20\%$ systematic error in distance creates a trend that is quite similar to the dipole observed in the data at similar amplitudes. This would imply that our distances could be overestimated. Such systematic errors could be the effect of the assumed parallax zero-point shift. However, StarHorse assumes a parallax shift of 0.05 mas, in agreement with recent works (e.g. Zinn et al., 2018; Graczyk et al., 2019; Schönrich et al., 2019). If the distances were overestimated this would require a larger parallax shift. Furthermore, applying these errors to the data (bottom panels of Fig. 8) shows that the pattern still remains even with a $20\%$ systematic error in distance. If systematics are indeed the main factor creating this $V_{R}$ trend, there should be a systematic distance or proper motion factor that inverts the pattern to match that expected from a typical bar model (Model 1). In Fig. 9, we further study the effect of distance systematics applied to the data. The maps show how a dipole is gradually created in the outer disc with decreasing distance and the one in the inner disc decreases in amplitude. The trends in the inner disc are never reversed (as seen in the second top panel of Fig. 8), even if we consider a $40\%$ systematic error. This indicates that the observed systematic shift in $V_{R}$ as a function of $\phi$ could be a real kinematic feature.

4.4.3 The effect of Sagittarius

We next study Model 2 to find out how an external perturbation from the Sagittarius dwarf galaxy modifies the inner disc velocity field.

The top row of Fig. 10 shows the variation of radial velocity with Galactic radius for different azimuthal slices (similar to Fig. 5 but for the range $3<R<7$ kpc) for Model 1 (left), Gaia DR2 (middle), and Model 2 (right). The bottom row shows the same samples but as maps in the $x\text{-}y$ plane.

We now see more clearly that Model 1, which lacks a significant external perturbation at the time period we consider, cannot reproduce the Galactocentric radial velocity structure in the bar vicinity. While the data reveal a region of high negative bulk velocity at azimuths lagging the bar and high positive velocity ahead of the bar, Model 1 showed the opposite trends, which we argued earlier are typical for a steady-state bar. This indicates that the Galactic bar and spiral arms cannot reproduce the azimuthal $V_{R}$ gradient alone.

The right column in Fig. 10, however, shows that the disc perturbed by Sagittarius matches the data in the trends, including the concave behaviour of the largest azimuthal bins (orange and red curves). In other words, the radial velocity field streaming motions on each side of the closer bar half are completely reversed compared to a steady-state model. We note that these trends in the inner disc can also be caused by distance systematics as discussed in Section 4.4.2 but assuming that is not the case, we find a good match in Model 2.

Looking at the bottom panels of Fig. 10, we find a remarkable match of the entire $V_{R}$ field between Gaia DR2 and Model 2, all the way out to 14 kpc. This includes the bridge between 8 and 10 kpc with $V_{R}\approx 0$ km s*-1*, also recently showed for the same simulation by Laporte et al. (2019), but note that here we have rescaled radius down to $80\%$. At a radius of about 14 kpc we find an arch of outward streaming motions in data and both models, although a better match again results from Model 2, except for the amplitude being larger by a factor of $\sim 2$.

4.4.4 A recently formed/evolved bar

The morphology of the Model 2 bar region at this time output can be seen in the bottom left panel of Figure 4 by Laporte et al. (2019). Note that the bar formation in this model started about 1 Gyr ago, triggered by the impact. This time is not sufficient for the disc to respond fully to the bar perturbation. Minchev et al. (2010) showed that, up to 2 Gyr after bar formation the $x_{1}(1)$ and $x_{1}(2)$ orbits at the bar OLR precessed in the bar reference frame thus giving rise to two low-velocity streams consistent with Pleiades and Coma Berenices or Pleiades and Sirius, depending on the time since bar formation (in addition to Hercules). This is unlike the expected behaviour of a steady-state bar, where these two orbital families are fixed in the bar reference frame along the bar major axis ($x_{1}(1)$) or perpendicular to it ($x_{1}(2)$). It is possible that such effects are present in the bar region as well, which could explain the inversion in the radial velocity field. A closer inspection of the central regions in both models shows that the inner $1\text{-}2$ kpc are consistent in the expected trends from a steady-state bar. However, at radii nearing the bar end ($\sim 4$ kpc in Model 1 and $\sim 3.5$ kpc in Model 2) Model 1 never reproduces the large area of positive motions needed to match the data, for any time output.

It is also interesting to note that the bar farther half does not show the same $V_{R}$ trends as the near one for Model 2, further strengthening our conclusion that the bar is time evolving. This can be contrasted to the near-perfect symmetry found in the central 2 kpc of the Model 1 disc. Future work will be dedicated to more detailed analyses.

4.5 Breathing and bending modes

In this Section, we focus our attention to the vertical velocity distribution, which has a complex structure with origins still under debate. In particular, we study the $V_{z}$ distribution by computing the breathing and bending modes in our sample and comparing these to our Models. Disentangling the observed mode, will helps us to put some constraints on the origins of the perturbation, as bending modes are attributed mainly to external perturbations and breathing modes are mainly induced by internal bar or spiral perturbations.

We define a breathing and bending mode based on the convention of Widrow et al. (2014):

[TABLE]

and

[TABLE]

We note that the structure resulting from the application of the above equations to the data and models will not necessarily signify breathing and bending modes, which need self gravity to propagate, but could be the result of phase mixing due to an external perturbation (such as Sagittarius) as discussed by de la Vega et al. (2015).

Figure 11 displays the breathing (top) and bending (bottom) face-on maps for different $z$ ranges, obtained from 100 Monte Carlo iterations as in Fig. 3. The breathing pattern observed close to the Galactic plane ($0<\left|z\right|<0.3$ kpc, top 1st column) is consistent with a median $\Delta V_{z}\approx 0$. This indicates bulk vertical motions pointing in the same direction and with similar amplitude above and below the plane and thus no breathing mode. This is in agreement with the results found by Bennett & Bovy (2019) within $\left|z\right|<1$ kpc with an amplitude $<1$ km s*-1* (half the amplitude compared to the breathing mode definition used in this work). However, at higher $z$ (top 2nd and 3rd column), the breathing mode seems to take shape, with mostly positive velocities at larger distances. This becomes more evident on the top 4th column showing amplitudes beyond $5$ km s*-1* in the overall map.

Unlike the breathing mode, all the bending mode patterns in the bottom panels of Figure 11 display the same trend, showing clear signatures of negative velocities at $R\lesssim 9$ kpc and positive outside. These patterns are quite similar to the $V_{z}$ map shown in the right panel of Fig. 3, implying that the bending mode is the dominant mode in our sample.

In Fig. 12 we compare two snapshots from Model 1 with our data. The data volume is illustrated by the black squares on each panel. Aside from the breathing mode at $0.7<\left|z\right|<1.5$ kpc, the breathing and bending modes at t=12.86 Gyr seem to follow a similar pattern as the data. However, just $37.5$ Myr later at t=12.90 Gyr, the bending pattern strongly changes. Similar results are seen in Fig. 13 with Model 2. At t= 6.9 Gyr, after virial radius crossing and at all heights, we observe patterns which follow similar trends as the data. In contrast, at t=6.48 Gyr, roughly 0.42 Gyr before, the disc is strongly bent and does not match the observations. These strong changes show that understanding the origin of the $V_{z}$ patterns can be rather complex. Beyond the scope of this work, a detailed analysis of individual perturbations of the Galactic bar, the spirals and external perturbations is required to understand the observed vertical motions.

5 Summary and Conclusions

In this work we have used the high accuracy proper motions and line-of-sight velocity obtained from Gaia DR2 combined with two sets of distances to study the three-dimensional velocity distribution of stars across a large portion of the Milky Way disc. The distances were obtained via Bayesian inference with two different approaches. Using the observed Gaia parallax $\varpi$ and the radial velocity spectrometer magnitude $G_{RVS}$ as priors (McMillan, 2018) and using the StarHorse code, which combines the Gaia astrometric information with multiband photometric information and a number of Galactic priors (Anders et al., 2019). In Section 2.3, we demonstrated how both estimates improve greatly the volume and accuracy of stellar distance compared to the naïve use of the inverse parallax as distance estimator: the observed volume was roughly doubled compared to previous authors using the inverse parallax and the mean uncertainties at $R>3$ kpc are at the level of $\sim 20\%$ instead of $\sim 60\%$ when using the inverse parallax (see left panel of Fig. 1).

We summarize our main results as follows:

• •

In both subsamples we identified asymmetries in the stellar velocity distribution. We studied the $V_{R}$ velocity distribution by dividing our sample in different azimuthal slices and heights from the plane. In the outer bulge/inner disc, we found variations up to $\pm 100$ km s*-1*.

• •

Inside the solar circle, we identified a strong radial $V_{R}$ gradient, transitioning smoothly from $-16$ km s*-1* kpc*-1* at an azimuth of $-45^{\circ}<\phi<-30^{\circ}$ lagging the Sun-Galactic centre line to $+16$ km s*-1* kpc*-1* at an azimuth of $30^{\circ}<\phi<45^{\circ}$ ahead of the solar azimuth.

• •

To understand the origins of such a gradient, we analysed the effect of introducing systematic errors in distance and proper motion. We found that distance systematics can create a similar pattern in the idealized case of a Galaxy composed only of circular orbits. However, when applied to the data, even correction for a large systematic shift did not remove the observed gradient in the inner disc.

• •

Assuming the gradient to be a real kinematic feature, we compared the data with two Milky Way models: a zoom-in simulation in the cosmological context, lacking recent massive external perturbers (Model 1) and an $N$-body simulation, which considers a Sagittarius-like interaction with the Milky Way (model 2). We found, that we were not able to reproduce the observed $V_{R}$ azimuthal gradient by considering mainly the internal effects of the Galactic bar and spiral arms. Instead, we demonstrated that introducing a major perturbation, such as the impact of Sagittarius, could reproduce a similar velocity field as in the observations (see Fig. 10).

• •

The velocity field can be reversed in the bar region only if the bar were currently forming or a previously existing one has been recently strongly perturbed, both of which could have been caused by Sagittarius. This suggest that the reversed $V_{R}$ field in the inner disc is a strongly time-dependent phenomenon, as a steady-state bar must exhibit the $V_{R}$ streaming of Model 1. As no external perturbations on the scale of Sagittarius are expected in the next several Gyr (except for the Sagittarius itself, who is constantly losing mass), we predict that the $V_{R}$ field should reverse in the next $1\text{-}2$ Gyr to match Model 1.

• •

Consistent with the results first shown by C18, we confirmed that the vertical velocity structure is a combination of vertical modes with a breathing mode inside the solar circle ($R\approx 6.5$ kpc) and a bending mode outside ($R\approx 13$ kpc). However, as Gaia Collaboration et al. (2018b) previously showed, the bending mode differed from the one found at $R\approx 8.5$ kpc by C18. With a larger dataset expanding to larger distances, we observed a bending mode with positive $V_{z}$ instead of negative.

• •

In the $x\text{-}y$ plane, we showed that, although the amplitude varies, the structure of the breathing and bending modes does not show any significant variations with distance from the Galactic plane. We showed that in both our models the breathing and bending patterns can drastically change in short time intervals, showing the complexity of understanding the origin of vertical perturbations. Thus, further and more detailed modelling is necessary to disentangle how internal and external perturbers individually affect the vertical velocity distribution in the Milky Way.

Acknowledgements

I. Carrillo is grateful to E. Poggio for discussions concerning the Galactic warp.

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.