# Is the Milky Way still breathing? RAVE–Gaia streaming motions

I. Carrillo,¹★ I. Minchev,¹ G. Kordopatis,¹,² M. Steinmetz,¹ J. Binney,³ F. Anders,¹ O. Bienaymé,⁴ J. Bland-Hawthorn,⁵ B. Famaey,⁴ K. C. Freeman,⁶ G. Gilmore,⁷ B. K. Gibson,⁸ E. K. Grebel,⁹ A. Helmi,¹⁰ A. Just,⁹ A. Kunder,¹,¹¹ P. McMillan,¹² G. Monari,¹,¹³ U. Munari,¹⁴ J. Navarro,¹⁵† Q. A. Parker,¹⁶,¹⁷ W. Reid,¹⁸,¹⁹ G. Seabroke,²⁰ S. Sharma,⁵ A. Siebert,⁴ F. Watson,²¹ J. Wojno,¹ R. F. G. Wyse²² and T. Zwitter²³

*Affiliations are listed at the end of the paper*

Accepted 2017 December 24. Received 2017 December 21; in original form 2017 October 11

## ABSTRACT

We use data from the Radial Velocity Experiment (RAVE) and the Tycho-Gaia astrometric solution (TGAS) catalogue to compute the velocity fields yielded by the radial ($V_R$), azimuthal ($V_\phi$), and vertical ($V_z$) components of associated Galactocentric velocity. We search in particular for variation in all three velocity components with distance above and below the disc mid-plane, as well as how each component of $V_z$ (line-of-sight and tangential velocity projections) modifies the obtained vertical structure. To study the dependence of velocity on proper motion and distance, we use two main samples: a RAVE sample including proper motions from the Tycho-2, PPMXL, and UCAC4 catalogues, and a RAVE–TGAS sample with inferred distances and proper motions from the TGAS and UCAC5 catalogues. In both samples, we identify asymmetries in $V_R$ and $V_z$. Below the plane, we find the largest radial gradient to be $\partial V_R/\partial R = -7.01 \pm 0.61 \text{ km s}^{-1} \text{ kpc}^{-1}$, in agreement with recent studies. Above the plane, we find a similar gradient with $\partial V_R/\partial R = -9.42 \pm 1.77 \text{ km s}^{-1} \text{ kpc}^{-1}$. By comparing our results with previous studies, we find that the structure in $V_z$ is strongly dependent on the adopted proper motions. Using the Galaxia Milky Way model, we demonstrate that distance uncertainties can create artificial wave-like patterns. In contrast to previous suggestions of a breathing mode seen in RAVE data, our results support a combination of bending and breathing modes, likely generated by a combination of external or internal and external mechanisms.

**Key words:** Galaxy: disc – Galaxy: kinematics and dynamics – Galaxy: structure.

## 1 INTRODUCTION

To a first approximation, the Milky Way disc is assumed to be axisymmetric and in equilibrium (e.g. Reid et al. 2009; Gnedin et al. 2010; McMillan & Binney 2010). However, in the last two decades with the acquisition of high-quality spectroscopic and astrometric data, and expansion of the surveyed volume around the Sun, non-negligible deviations from axisymmetry have become apparent. One indication of such asymmetries came from the clear identification in data from the Hipparcos astrometric satellite (Perryman & ESA 1997) of overdensities in the velocity space of local stars (Chereul,

Creze & Bienaymé 1998; Dehnen 1998; Asiain et al. 1999). Subsequent analysis of the local velocity field has revealed that the most prominent moving groups (or streams) in the solar neighbourhood probably have a dynamical origin since within a given stream there is a mixture of chemical abundances and ages (e.g. Famaey et al. 2005; Antoja et al. 2008; Famaey, Siebert & Jorissen 2008). Dehnen (2000) and Fux (2001) showed that the Galactic bar can successfully reproduce the Hercules stream if the Sun is situated close to the bar’s outer Lindblad resonance, while a solar position near the inner ultra-harmonic 4:1 resonance of a two-armed spiral density wave can similarly create resonant structures consistent with the Pleiades/Hyades and Coma Berenices moving groups (Quillen & Minchev 2005; Pompéia et al. 2011). While internal perturbations can explain the low-velocity moving groups, high-velocity streams (such as Arcturus, which lags the local standard of rest by

<footnote>∗ E-mail: icarrillo@aip.de</footnote>

<footnote>† Senior CIfAR Fellow.</footnote>
~100 km s⁻¹; Williams, Freeman & Helmi 2009) have been related to external perturbations from an infalling satellite galaxy (Minchev et al. 2009; Gómez et al. 2012a,b; D'Onghia et al. 2016). A crucial diagnostic in understanding what causes these streaming motions is how structure of the velocity space varies with location in the disc (Dehnen 2001; Minchev et al. 2010; McMillan 2013; Antoja et al. 2014; Monari et al. 2017).

Another manifestation of non-axisymmetries in the Milky Way disc is in-plane stellar streaming motions detected in the extended solar neighbourhood. Using line-of-sight velocities measured in the Radial Velocity Experiment (RAVE; Steinmetz et al. 2006), Siebert et al. (2011b) measured a gradient in the mean Galactocentric radial velocity (|∂⟨V<sub>R</sub>⟩/∂R| ≥ 3 km s⁻¹ kpc⁻¹). Siebert et al. (2012) attributed this gradient to a two-armed spiral perturbation in which the Sun is again close to the inner ultraharmonic 4:1 resonance. More recently, Monari et al. (2014) used test-particle simulations to show that the radial gradient found in RAVE could alternatively be caused by the Galactic bar.

In addition to stellar bulk motions in the disc plane, structure has been found also in the direction perpendicular to the Galactic disc. Widrow et al. (2012) and Williams et al. (2013, hereafter W13) used SEGUE (Yanny et al. 2009) and RAVE, respectively, to detect a wave-like pattern in the mean vertical velocity of stars near the Sun. Similar vertical asymmetries were also found by Carlin et al. (2013) using data from the LAMOST survey (Cui et al. 2012) and Xu et al. (2015) with data from the Sloan Digital Sky Survey (Aihara et al. 2011).

While the radial-velocity gradient has been associated with internal perturbations, the origin of vertical velocity structure is debatable. Gómez et al. (2013) attributed the vertical patterns to the passing of the Sagittarius dwarf galaxy (Ibata, Gilmore & Irwin 1994), while Widrow et al. (2014) noted that they could also be caused by a dark matter subhalo. In addition to external perturbations, vertical streaming motions have now been shown to result also from internal mechanisms such as the Galactic bar (Monari, Famaey & Siebert 2015) and spiral arms (Faure, Siebert & Famaey 2014; Monari, Famaey & Siebert 2016). However, according to Monari et al. (2015), the Galactic bar is unlikely to induce mean vertical motions greater than ~0.5 km s⁻¹ in the outer disc, and therefore does not explain the observed vertical motions in the solar neighbourhood.

Another possible explanation of the observed vertical motions is the Galactic warp. The Galactic disc becomes warped beyond the solar circle: the H<sub>1</sub> disc at R ≥ 10 kpc is warped such that the Sun lies near the line of nodes and the z coordinate of the centre of the gas layer increases in the direction of Galactic rotation (Binney & Merrifield 1998; Levine, Blitz & Heiles 2006; Kalberla et al. 2007). Consequently, if (as is likely) the pattern speed of the warp is smaller than the circular frequency of the Sun, gas at R > R₀ and the Sun's azimuth should be moving downwards (V<sub>z</sub> < 0). A warp in the stellar disc has been detected from R₀ outwards (e.g. Dehnen 1998; Drimmel & Spergel 2001; Robin, Reylé & Marshall 2008). None the less, Poggio et al. (2017) used OB stars and proper motions from Hipparcos and Gaia DR1 (Gaia Collaboration et al. 2016b) to develop a simple model of a stable long-lived warp and found that vertical motions in the disc cannot be explained by the warp, suggesting that the warp's kinematic signal could be overwhelmed by other systematic motions.

More generally, the analysis of 16 high-resolution, fully cosmological simulations of the evolution of individual Milky Way-sized galaxies led Gómez et al. (2016) to conclude that vertical disc

asymmetries are ubiquitous in A cold dark matter (ACDM) cosmology. In fact, ~70 per cent of the analysed simulations exhibited strong vertical patterns, with amplitudes exceeding 2 kpc from the disc mid-plane. Hence, vertical motions are predicted by ACDM cosmology.

Using RAVE red-clump (RC) stars and a compilation of proper motions, W13 studied the 3D velocity distribution in the extended solar neighbourhood, focusing on north-south differences, and detected a rarefaction-compression behaviour in the vertical velocity pattern. Such a pattern, which has odd parity in the V<sub>z</sub> distribution and even parity in the density distribution, is known as a breathing mode. In contrast, even parity in V<sub>z</sub> with odd parity in the density distribution is known as a bending mode. By identifying the observed mode, we hope to constrain the nature of the exciting perturber. Breathing modes have been found to be induced naturally by the Galactic bar and spiral arms (Faure et al. 2014; Monari et al. 2015, 2016), whereas bending modes are attributed mainly to external perturbations. While more recently, Chequers & Widrow (2017) found that bending waves could also arise without excitation by a satellite/merging event, Widrow et al. (2014) used a toy-model simulation of disc-satellite interactions to show that a passing satellite galaxy could produce bending or breathing modes depending on the vertical velocity of the satellite as it passes through the Galactic disc. This shows that the picture can actually be rather complex, as external perturbations also generate internal spiral perturbations that can in turn excite breathing modes. While both bending and breathing modes can arise from external interactions, the observations so far have been found to be consistent with a breathing mode. Another interpretation was found by de la Vega et al. (2015), who showed that phase wrapping in the disc following the passage of a satellite can look very similar to a bending or a breathing mode.

ESA's mission Gaia (Gaia Collaboration et al. 2016a) is acquiring highly accurate parallaxes and proper motions for over a billion sources brighter than magnitude 20.7 in the G band. In this paper, we use early mission data from the first data release (Gaia DR1) combined with the line-of-sight velocities of the fifth data release from RAVE (DR5; Kunder et al. 2017) to study the three-dimensional velocity distribution of stars in the extended solar neighbourhood. We extend the analysis of W13 by increasing the number of stars through inclusion of stars that are not in the RC and using better proper motions and distances. Since uncertainties in velocity fields are dominated by errors in distance and proper motion (RAVE uncertainty in line-of-sight velocity <2 km s⁻¹), our goal is to use proper motions from Gaia to scrutinize with higher accuracy the vertical velocity pattern of the extended solar neighbourhood. In this way, we aim to determine whether the Milky Way exhibits a breathing or a bending mode.

This paper is structured as follows: in Section 2, we introduce our coordinate conventions and explain how we selected stars. In Section 3, we present the radial, azimuthal, and vertical velocity distribution obtained by W13 and the velocity distributions obtained from RAVE DR5 and different proper-motion catalogues. In Section 4, we study the composition of the vertical velocity pattern and how errors affect its structure. In Section 5, we examine the differences between RAVE DR5 and samples that use data from DR5 combined with the Tycho-Gaia astrometric solution (TGAS) catalogue (RAVE-TGAS samples) and present the north-south asymmetries obtained using the most accurate estimates of distance and proper motion. Finally, Section 6 contains a summary of our results and conclusions.
## 2 COORDINATE SYSTEMS AND DATA SELECTION

### 2.1 Coordinate systems

We compute the heliocentric rectangular components of the Galactic space velocity *U*, *V*, and *W* using the right-handed coordinate system, with *U* positive towards the Galactic Centre, *V* positive in the direction of Galactic rotation, and *W* positive towards the North Galactic Pole. The method used to derive the Galactic space velocities is described in detail in Johnson & Soderblom (1987). The Galactocentric cylindrical velocity components (*V*<sub>R</sub>, *V*<sub>φ</sub>, and *V*<sub>z</sub>) are computed following the coordinate transformation given in appendix A of W13.

For the solar Galactocentric distance, we adopt *R*<sub>0</sub> = 8 kpc (Reid 1993) and use the estimate of the peculiar velocity of the Sun obtained by Schönrich, Binney & Dehnen (2010):

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

With these values and the proper motion of Sagittarius A*, $\mu_{l_{Sgr A*}} = 6.379$ mas yr<sup>-1</sup> (Reid & Brunthaler 2004), we obtain

$$ V_{\odot} + V_{LSR} = 4.74 R_{0} \mu_{l_{Sgr A*}}, \quad (2) $$

which yields a value of $V_{LSR} \sim 230$ km s<sup>-1</sup> for the circular velocity of the local standard of rest (LSR). We also tested our results with $R_0 = 8.3$ kpc and $V_{LSR} = 240$ km s<sup>-1</sup>. The changes in $V_R$ and $V_\phi$ were negligible, and since $V_z$ is independent of $V_{LSR}$, it suffers no changes aside from the position in *R*.

### 2.2 RAVE data selection

RAVE is a spectroscopic magnitude-limited ($9 \lesssim I \lesssim 12$) Southern hemisphere survey, which collected data from 2003 to 2013. It acquired line-of-sight velocities and stellar parameters (effective temperature, surface gravity, overall metallicity) in order to probe the chemical and dynamical evolution of the Milky Way. In this paper, we use its latest data release, Data Release 5 (Kunder et al. 2017). DR5 contains 520 781 spectra of 457 588 unique stars. It includes observations that were previously discarded, resulting in ~30 000 more spectra than the fourth data release (DR4; Kordopatis et al. 2013a). DR5 also improves on the distance pipeline of DR4 (Binney et al. 2014), especially as regards metal-poor stars, and applies a new calibration of its stellar parameters that improves their accuracy by up to 15 per cent compared to DR4 (see Section 6 in Kunder et al. 2017 for details). The proper motions in DR5 are taken from the Tycho-2 (Høg et al. 2000), PPMXL (Roeser, Demleitner & Schilbach 2010), and UCAC4 (Zacharias et al. 2013) catalogues.

For our data selection (RAVE sample), we exclude stars without measured proper motions and stars with a proper-motion error $\sigma_{\mu_a}, \sigma_{\mu_b} > 15$ mas yr<sup>-1</sup> (catalogue independent) to remove outliers. Further cuts applied were to include only stars with a signal-to-noise ratio $SNR_K > 20$ and to remove stars whose spectral morphological flag indicates peculiar spectra (c1,c2,c3 ≠ 'n'); these cuts improve the quality of our data and remove most spectroscopic binaries, spectra with continuum abnormalities, and other unusual spectra. We also removed stars with an error of heliocentric radial velocity $eHRV > 10$ km s<sup>-1</sup>, cool dwarfs ($\log g > 3.5$ and $T_{eff} < 4000$ K), metal-poor stars ([M/H] < -1.2 dex), hot stars ($T_{eff} > 7400$ K), and spectra for which atmospheric parameters could not be determined reliably (Algocv = 1), since in that case, the derived distances

(obtained through isochrone fitting; see Binney et al. 2014) were not reliable either.

Finally, stars at Galactic latitude $|b| < 10^\circ$ were removed due to uncertain interstellar extinctions possibly affecting the derived distances. Our data selection applies to all proper-motion catalogues, with the exception of the Galactic latitude cut for the PPMXL catalogue, from which we remove all stars with $|b| < 20^\circ$. This cut was suggested by Vickers, Röser & Grebel (2016), who presented a correction to the PPMXL proper motions. Although we do not make use of the code provided to correct the proper motions, the Galactic latitude cut solves the high discrepancy between the velocity fields of PPMXL and all other catalogues. The applied cuts produce, respectively, samples of 173 162, 116 632, and 105 331 stars when UCAC4, PPMXL, and Tycho-2 proper motions are used. The adopted distance of each star is the inverse of the expectation of the star's parallax, which Binney et al. (2014) recommended as the most reliable distance estimate.

### 2.3 The TGAS catalogue

An important aspect of RAVE DR5 is that with almost 256 000 spectra that overlap with a TGAS star, this data release currently has a larger overlap with the TGAS catalogue than any other spectroscopic survey. TGAS is the primary astrometric data set included in Gaia DR1, which includes proper motions, sky positions, and parallaxes of 2057 050 stars from the Hipparcos (Perryman & ESA 1997; van Leeuwen 2007) and Tycho-2 catalogues. In addition to the TGAS proper motions and positions, we use two distance estimates for our RAVE–TGAS sample: the purely astrometric distances of Astraatmadja & Bailer-Jones (2016, ABJ), which employ an anisotropic prior based on a three-dimensional density model of the Milky Way, and the improved RAVE spectrophotometric distances from McMillan et al. (2017), which take advantage of both TGAS parallaxes and effective temperatures from the infrared flux method (IRFM; Blackwell, Shallis & Selby 1979; Casagrande et al. 2010).

The RAVE–TGAS sample with McMillan distances involves a cut additional to those listed in Section 2.2: we removed stars with $\log g < 2$, because distances to low-gravity stars are systematically overestimated (see McMillan et al. 2017). The data selection results in a sample of 68 477 stars.

For the RAVE–TGAS sample with ABJ distances, we use the selection criteria described in Section 2.2 less the cuts on heliocentric velocity error, effective temperature, metallicity, surface gravity, or the convergence of the atmospheric parameters, because the ABJ distance estimates do not rely on the accuracy of the atmospheric parameters. The RAVE–TGAS sample with ABJ distances comprises 184 954 stars.

To further study the proper-motion differences in the RAVE–TGAS samples, we include the newly obtained UCAC5 (Zacharias, Finch & Frouard 2017) proper motions, which were derived by combining UCAC data with Gaia DR1.

## 3 RAVE FIFTH DATA RELEASE

W13 studied the differences in the Galactocentric velocity distributions for stars above and below the Galactic plane. For this, they used a sample of 72 635 RC stars obtained from the RAVE internal third data release (Siebert et al. 2011a). For comparison purposes, we use a similar data set obtained from the internal fourth data release used in Kordopatis et al. (2013b). We apply the same cuts described in Section 2.2. The data selection results in a sample of
<figure><img src="image_1.png" /><figcaption>**Figure 1.** Median azimuthal velocity fields obtained using different proper-motion catalogues, as indicated on top of each panel. Each velocity field is shown in $(0.1 \text{ kpc})^2$ pixels covering a box of up to $(5 \text{ kpc})^2$ around the Sun, which we locate at $(R, z) = (8,0) \text{ kpc}$. All pixels are smoothed by computing the median of the velocities of all stars located in a square of size $(0.2 \text{ kpc})^2$ with a minimum of 50 stars. For comparison purposes, the bottom-right panel shows the velocity field obtained by Kordopatis et al. (2013b), which is very similar to the one obtained using only RC stars by W13. The velocity fields of all proper motions are consistent with a decrease of $V_\phi$ with increasing distance from the plane.</figcaption></figure>

127 722 stars (Williams/Kordopatis sample). The resulting velocity fields differ mainly from W13 in being based on the full RAVE sample instead of only RC stars. This allows for a better comparison to the velocity fields obtained with our RAVE samples.

In addition to the advantages of using RAVE DR5 data over RAVE DR4 (see Section 2.2), a further improvement to previous works is the choice of proper motions used to compute the velocity fields. Indeed, W13 and Kordopatis et al. (2013b) used an inhomogeneous compilation of proper-motion values, a remnant of RAVE DR3 approach; for a star present simultaneously in several proper-motion catalogues available at the time, the proper motion with the smallest reported uncertainty was adopted. However, it is to be feared that combining proper motions measured with different techniques imprinted unknown systematics on the data set. Consequently, in this work, we do not mix the catalogues and use individual proper motions. The data presented here, with more homogeneous proper motions and better distance estimates, therefore yield velocities that

are less systematically biased than the ones presented in W13 and Kordopatis et al. (2013b).

In Fig. 1, we present maps of the median azimuthal velocity for each proper-motion catalogue, as well as the one obtained by Williams/Kordopatis. The velocity maps are displayed in $(0.1 \text{ kpc})^2$ pixels covering up to 3 kpc from the Sun. Each pixel is smoothed by computing the median velocities over a box of $(0.2 \text{ kpc})^2$ with a minimum of 50 stars. The main characteristic observed in $V_\phi$ is the asymmetric drift. Fig. 1 displays a symmetrical distribution of velocity in the northern and southern Galactic hemispheres, with velocities lagging that of the LSR more at larger distance from the Galactic plane. In all proper-motion samples, $V_\phi$ behaves similarly to the results obtained by Williams/Kordopatis, which used a compilation of proper-motion catalogues: within $0 < |z| < 0.5 \text{ kpc}$, the median $V_\phi$ is of the order of 220 km s$^{-1}$ and decreases with distance from the plane to 190 km s$^{-1}$ at $|z| \approx 1 \text{ kpc}$ and up to 130 km s$^{-1}$ at $|z| \approx 2 \text{ kpc}$.
<figure><img src="image_2.png" /><figcaption>**Figure 2.** Same as Fig. 1 but for $V_R$. Independent of the proper-motion catalogue used, the velocity field above the plane is structured similarly, with mostly positive velocities everywhere. In contrast, below the plane, the structure varies significantly between proper-motion catalogues.</figcaption></figure>

The dependence of $V_R$ on proper motion is presented in Fig. 2. Inwards of the solar radius, independent of the proper-motion catalogue used, stars both above and below the plane move away from the Galactic Centre with median velocities up to $V_R = 30$ km s$^{-1}$. Below the plane in the range $-1 < z < 0$ kpc, we estimate the largest radial-velocity gradient in the Williams/Kordopatis velocity field to be $\partial V_R/\partial R = -7.98 \pm 0.87$ km s$^{-1}$ kpc$^{-1}$, a value similar to that found by Siebert et al. (2011b). Above the plane in the range $0 < z < 1$ kpc, $R \lesssim 8$ kpc, we estimate the largest value to be $\partial V_R/\partial R = -9.08 \pm 1.48$ km s$^{-1}$ kpc$^{-1}$. In Table 1, we report the largest $V_R$ gradients obtained from each proper-motion catalogue. The $(R, z)$ planes used to obtain these gradients are presented in Fig. 3.

The radial gradients in the northern and southern Galactic hemispheres are consistent with one another within $1\sigma$. Thus, with these subsamples, we find no north-south variation in the strength of the gradient within 8 kpc. Consistent with the gradient reported by Siebert et al. (2011b) and Anguiano et al. (2017), the gradients lie in the range $\sim 5-10$ km s$^{-1}$ kpc$^{-1}$.

<table><caption>**Table 1.** Largest gradient in $V_R$ obtained from the Williams/Kordopatis sample (W/K) and from the RAVE DR5 sample combined with three proper-motion catalogues. The values are computed using a least-squares fit (see Fig. 3). The uncertainties are taken from the covariance matrix.</caption><thead><tr><th>Catalogue</th><th>North<br>∂V<sub>R</sub>/∂R<br>(km s<sup>-1</sup> kpc<sup>-1</sup>)</th><th>South<br>∂V<sub>R</sub>/∂R<br>(km s<sup>-1</sup> kpc<sup>-1</sup>)</th></tr></thead><tbody><tr><td>PPMXL</td><td>-7.59 ± 1.15</td><td>-5.99 ± 0.81</td></tr><tr><td>Tycho-2</td><td>-8.53 ± 1.82</td><td>-5.75 ± 0.47</td></tr><tr><td>UCAC4</td><td>-5.53 ± 1.24</td><td>-4.10 ± 0.89</td></tr><tr><td>W/K</td><td>-9.08 ± 1.48</td><td>-7.98 ± 0.87</td></tr></tbody></table>

Unlike $V_R$, $V_z$ depends largely on proper motions, which are vulnerable to catalogue systematics and distance errors. Fig. 4 illustrates the strong dependence of the vertical velocity structure on the proper-motion catalogue used. Inside the solar radius, all velocity
<figure><img src="image_3.png" /><figcaption>Figure 3. Largest negative gradient in $V_R$ estimated in each proper-motion catalogue for stars above (left) and below the plane (right) computed in bins of 0.2 kpc using the shown (R, z) planes.</figcaption></figure>

fields feature a positive stellar motion above the plane. The velocity fields obtained with Tycho-2 and PPMXL show similar structure to that obtained by Williams/Kordopatis: inside the solar radius, we see upward stellar motion above the plane and downward motion below it. By contrast, at $R < R_0$, the UCAC4 velocity field shows upward motion both above and below the plane. Hence, inside the solar radius, the PPMXL, Tycho-2, and Williams/Kordopatis velocity fields have the signature of a breathing mode as discussed by Widrow et al. (2012, 2014), while the UCAC4 shows that of a bending mode.

Outside the solar radius, gradients in $V_z$ are weaker. PPMXL and Tycho-2 show no clear signature, while UCAC4 and Williams/Kordopatis show weak evidence for opposite patterns: in UCAC4, there appears to be upward motion above the plane and downward motion below it, while the Williams/Kordopatis field shows oppositely directed velocities. Either pattern would be characteristic of a breathing mode.

Although the combination of breathing and bending modes outside and inside $R_0$ shown by UCAC4 has not previously been reported in the literature, it could arise from a combination of external and internal excitations. For example, de la Vega et al. (2015) used non-interacting test-particle integrations to study the effects of a dwarf galaxy passage on the stellar epicyclic motions and the resulting streaming caused by the subsequent phase wrapping. Since their simulations lacked self-gravity, structures could not arise from bending or breathing modes. None the less, the velocity distribution observed at an azimuthal angle $\Theta = 270^\circ$ in their fig. 8 is consistent

with the streaming motions caused by a combination of breathing and bending mode.

However, the observed dependence of vertical velocity on proper motions indicates the need for more accurate data. The ESA Gaia mission is thus crucial in understanding the origins of the vertical streaming motions observed in the Milky Way.

## 4 DECOMPOSITION OF THE VERTICAL VELOCITY PATTERN

In this section, before studying how the more precise Gaia data affect the vertical velocity pattern as a whole, we study the composition of the vertical velocity and the contribution of each component in the observed structure. In Johnson & Soderblom (1987), the Galactocentric vertical velocity $V_z$ is obtained by computing the heliocentric velocity and adjusting for the solar motion afterwards. However, by using the transverse velocity in Galactic coordinates ($V_l, V_b$) and applying the solar motion correction before computing $V_z$, we are able to separate the vertical velocity in two components – one involving the line-of-sight velocity corrected by the solar motion, $V_{los}$, and the other dependent on the distance, $d$, and the proper motion, $\mu$. The Galactocentric vertical velocity is thus given by

$$V_z = V_{los} \sin b + V_b \cos b, \quad (3)$$

where

$$V_b = 4.74 d \mu_b + (-U_{\odot} \cos l \sin b - V_{\odot} \sin l \sin b + W_{\odot} \cos b). \quad (4)$$

Since the small uncertainty in $V_{los}$ from RAVE does not affect the $V_z$ pattern, we are in a position to study how distance and proper motion affect the vertical velocity. First, we want to avoid the effects of disc perturbations. For this, we used Galaxia (Sharma et al. 2011) as a front end for the (axisymmetric) Besançon Galaxy Model (Robin et al. 2003). Galaxia uses isochrones from the Padova data base to compute photometric magnitudes for the model stars (Bertelli et al. 1994; Marigo et al. 2008), and it quickly generates from models mock observational catalogues that cover a specified area on the sky with any specified selection function. We use the Galaxia sample that was generated by Wojno et al. (2017) to obtain a mock stellar sample based on the RAVE selection function.

Fig. 5 shows the median vertical velocity obtained from the mock RAVE sample and its components. The top panels display $V_z$ considering all components (left), only $V_{los} \sin b$ (middle), and $V_b \cos b$ [right, see equation (3)]. Here, we are able to identify, in dilute form, the main patterns in $V_z$, in which at $R < 8$ kpc $V_{los} \sin b$ is negative above the plane and positive below it, and at $R > 8$ kpc it is positive above the plane and negative below it. As expected, such a pattern alone does not match any of the $V_z$ patterns found in Fig. 4, indicating a stronger proper-motion/distance dependence. The top-right panel of Fig. 5 does display key aspects of the pattern visible in the W/K panel of Fig. 4, suggesting that in W13 errors may have caused the proper-motion component of $V_z$ to dominate the more precise contribution from $V_{los}$. The $V_z$ subcomponents shown in the top-middle and top-right panels have an inverse symmetry, which results in the structure-free velocity field when summed up to estimate $V_z$ (top-left panel).

In the bottom panels of Fig. 5, we show how the distance, proper motion, and solar motion affect the transverse velocity component $V_b \cos b$ (equation 4). In the bottom-left panel, we set $\mu_b = 1$ mas yr$^{-1}$ in the $V_b \cos b$ term and do not consider the solar motion. This displays the effects of distance on $V_b$, where higher velocities are found at larger distances. Similarly, in the bottom-middle
<figure><img src="image_4.png" /><figcaption>**Figure 4.** Same as Fig. 1 but for $V_z$. The trend in vertical velocity varies significantly with proper-motion catalogue, especially below the plane. The arrows in the bottom-right panel show the direction of vertical motion consistent with a breathing mode. While some catalogues agree with the breathing mode observed by Williams/Kordopatis, others display a bending mode or a combination of bending and breathing modes not previously seen (top-right panel).</figcaption></figure>

panel, we set $d = 1$ kpc and exclude the solar motion term, in order to see the effects of $\mu_b$ alone. Here, it is shown that the proper motions contribute mainly with negative velocities aside from the volume below the plane at $R > 8$ kpc. Finally, the bottom-right panel displays only the contribution to $V_b$ from the solar motion. In all panels, the obtained structure is better understood by looking at equations (3) and (4). The term with $\cos b$ in the transverse velocity component gives always a positive value, so the negative structure is given by $V_b$. Here, the velocity pattern dependent on distance is also positive since distances are inherently positive. In contrast, the proper motions and the solar motion have negative values, with the solar motion remaining relatively constant between samples. The different $V_b$ patterns are thus dependent mainly on the proper motions, with the distance increasing its amplitude. This can be seen in Fig. 6, which displays the components of $V_z$ in the Williams/Kordopatis sample. Here, although the pattern obtained from the distance has lower velocities than in Galaxia, the pattern of proper motions is strong enough to dominate the contribution

from $V_{\text{los}}$, breaking the balance previously shown in Galaxia and generating the pattern of a breathing mode.

In Fig. 7, we study the effects of introducing errors to the mock data from the axisymmetric and relaxed Galaxia model. The first panel shows the Galaxia sample with no errors. In the second panel, we introduce artificial errors and draw new velocity values, from a normal distribution centred on the actual velocity and $\sigma_\mu = 5$ mas yr$^{-1}$. Although these errors are larger than the typical error in our DR5 samples (up to $\langle\sigma_\mu\rangle \sim 4$ mas yr$^{-1}$), the assigned values produce some clumps but do not affect the general structure of the sample.

The third panel shows the effect of adding a random distance error with a standard deviation $\sigma_d/d = 0.8$ (we chose such large errors to show a clear pattern, although this is already seen with $\sigma_d/d = 0.5$). This random error increases the volume of our sample and, since all distance values are positive, increases the value of the $V_b$ component (equation 4). This results in a pattern similar to the one obtained by Williams/Kordopatis. The right-hand panel
<figure><img src="image_5.png" /><figcaption>**Figure 5.** Decomposition of the median value of $V_z$ obtained from the RAVE-like sample from Galaxia assuming a magnitude limit $0 < I < 13$, with no colour restriction, where RAVE-like errors were applied to stellar parameters. The velocity fields are shown in (0.1 kpc)² pixels covering a box of up to (4 kpc)² around the Sun. The top panels display $V_z$ and its two main components based on equation (3). The bottom panels show the contributions to $V_b \cos b$ (distance, proper motion, and solar motion) from equation (4). These panels indicate how each component could affect the map of $V_z$. The transverse velocity depends mainly on the proper motion with an amplitude that increases with distance.</figcaption></figure>

displays a stronger more noisy pattern as a combination of both random errors. The obtained patterns in this section allow us to understand better how each component affects $V_z$. In Section 5, we will analyse how this compares to our RAVE–TGAS sample.

# 5 RAVE–TGAS SAMPLE

The superior parallaxes and proper motions obtained by the ESA Gaia mission should yield marked improvements on the velocity fields presented in Section 3. In Section 5.1, we first compare the distance estimates obtained purely from Gaia DR1 with the improved RAVE DR5 spectrophotometric distances. Then in Section 5.2 we compare the impact of using proper motions from RAVE DR5, TGAS, and UCAC5. In Section 5.3, we introduce the velocity fields obtained using the best estimates in distance and proper motion available for our RAVE–TGAS sample.

## 5.1 Distance estimate

As discussed in Section 2.3, our RAVE–TGAS sample makes use of two distance estimates: those of Astraatmadja & Bailer-Jones (2016) and McMillan et al. (2017). We do this on account of the problems discussed by Bailer-Jones (2015) associated with use of inverse parallax as a distance estimator.

If we assume measured parallaxes $\omega$ to be normally distributed with a standard deviation $\sigma_{\omega}$, and adopt the inverse parallax as the distance estimator, we encounter two important issues: (i) the estimator fails for negative $\omega$, even though these are valid measurements (see Bailer-Jones 2015 for further details). (ii) For fractional parallax errors $f_{Obs} = \sigma_{\omega}/\omega > 0.2$, using the inverse parallax creates a skewed distribution that gives a biased distance estimator. Fig. 8 shows a cumulative histogram of $f_{Obs}$ for our RAVE–TGAS sample as well as the full TGAS sample. The parallaxes of Gaia DR1 suffer from a possible offset by $\pm 0.1$ mas of the parallax zero-point and $\pm 0.2$ mas due to position and colour-dependent
<figure><img src="image_6.png" /><figcaption>**Figure 6.** Same as Fig. 5 but for the map of $V_z$ obtained by Williams/Kordopatis.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 7.** Median $V_z$ obtained from the Galaxia RAVE-like sample together with the effects of adding random errors in *d* and *μ* as described in the text. The second panel shows that by adding a random error in *μ* the $V_z$ pattern displays some clumps, while the third panel shows that adding distance errors substantially changes the map. The combination of both errors in the last panel results in a pattern similar to the one obtained by Williams/Kordopatis.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 8.** Cumulative distribution of the fractional parallax errors $f_{\text{Obs}}$ of stars in the full TGAS (blue line) and RAVE–TGAS (red line) samples. The dotted curves represent the same samples taking into account systematic errors in the parallax uncertainties. Although there is a significant number of stars with $f_{\text{Obs}} < 0.2$, using the inverse parallax as a distance estimate produces biased distances for at least ~50 per cent of the stars in both samples.</figcaption></figure>

systematics. A detailed description of these systematics is given in Lindegren et al. (2016, appendix E). It is therefore recommended to add a systematic error of ±0.3 mas to the parallax uncertainties (Gaia Collaboration et al. 2016b). The dotted curves show the same TGAS data but consider the systematic errors. As can be seen, the full TGAS and our RAVE–TGAS sample have a relatively large number of stars (~50 per cent) within $f_{\text{Obs}} < 0.2$. Nevertheless, using the inverse parallax in the RAVE–TGAS sample would lead to potentially biased distance estimates for a minimum of ~50 per cent of the stars. If we consider the systematic errors, we see that the distance estimate would be biased for at least ~60 per cent of the stars.

The inverse parallax is therefore a poor distance estimate for our RAVE–TGAS sample. If we select only the stars with $f_{\text{Obs}} < 0.2$ including systematic errors, we would limit the observed volume of stars within just 0.5 kpc around the Sun. The volume obtained is then not big enough to recognize a distinctive mode in the Galaxy. Hence, another approach to estimate distances is required. Recently, Astraatmadja & Bailer-Jones (2016) used Bayesian inference to estimate distances for the full TGAS sample. For this, they used two different priors: an exponentially decreasing space density prior and an anisotropic prior based on a three-dimensional density model of the Milky Way. The consistency of both priors was tested by comparing their distance estimates with the determined distances of 105 Cepheid variable stars from Groenewegen (2013) cross-matched with Gaia DR1 (see fig. 4 of Astraatmadja & Bailer-Jones 2016). It was shown that for distances $d < 2$ kpc, the Milky Way prior comes closer to the more precise Cepheid distance estimate than the exponentially decreasing space density prior. Beyond this limit, the latter performs better. Due to the volume of our samples ($d \lesssim 2$ kpc), we choose the distance estimate obtained from the Milky Way prior for further analysis in this work.

Recently, McMillan et al. (2017) determined new distances to RAVE stars also using Bayesian inference, but considering both TGAS parallaxes and stellar parameters derived from the RAVE spectra, with $T_{\text{eff}}$ values derived using the IRFM. Therefore, it is recommended to use this new estimate when combining RAVE and TGAS data.

<figure><img src="image_9.png" /><figcaption>**Figure 9.** Comparison of the distance estimates of McMillan, which combine TGAS parallaxes with spectrophotometric data, and those of ABJ using the Milky Way prior. The red line indicates a perfect match between distances. Under 1 kpc, about ~66 per cent of the distance estimates are in good agreement.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 10.** Comparison between relative uncertainties $\sigma_d/d$ in the distance estimates of McMillan and ABJ. The red line denotes equal $\sigma_d/d$ values. For most stars, the McMillan distances have the smaller relative uncertainties.</figcaption></figure>

Fig. 9 presents a comparison between the McMillan distances and the ABJ distance estimates with the Milky Way prior. The residuals between the two distance estimates have a mean value of $\mu = 0.10$ kpc and a dispersion of $\sigma = 0.36$ kpc. Within ~1 kpc, the agreement is very good. A comparison of the relative uncertainties $\sigma_d/d$ between both estimates is presented in Fig. 10. As can be seen, most stars have smaller relative uncertainties when using the McMillan estimate. Due to the consistency between the distance estimates, we want to further study the ABJ distance uncertainties. Since the McMillan and the ABJ distance estimates are not totally independent – both estimates use TGAS parallaxes to infer distances, we are not able to compare their differences and learn about biases in one measurement versus another or whether the error estimates are accurate. Therefore, we compare ABJ to RAVE DR5 distances to further study the ABJ distance uncertainties.

Here, as in Binney et al. (2014) and Kunder et al. (2017), we divide the sample into three groups of stars: giants (log $g < 3.5$), hot dwarfs (log $g > 3.5$ and $T_{\text{eff}} > 5500$ K), and cool dwarfs (log $g > 3.5$)
<figure><img src="image_11.png" /><figcaption>**Figure 11.** Histograms of the difference between the distance of ABJ and RAVE DR5. Stars are divided into giants (log g < 3.5), hot dwarfs (log g > 3.5 and T<sub>eff</sub> > 5500 K), and cool dwarfs (log g > 3.5 and T<sub>eff</sub> < 5500 K). The differences are computed using equation (5); the solid red curves are Gaussians of zero mean and unit dispersion, representing perfect results. To remove outliers, we use only values with |Δ<sub>ABJ-DR5</sub>| < 4. The most characteristic feature in all three panels is the non-zero mean, indicating a clear difference between both distance estimates and a dispersion less than unity, corresponding to overestimated uncertainties of one or both measurements.</figcaption></figure>

and T<sub>eff</sub> < 5500 K). The histograms of Fig. 11 show the differences between the ABJ distance and the RAVE DR5 distances as follows:

$$ \Delta_{\text{ABJ-DR5}} = \frac{d_{\text{ABJ}} - d_{\text{DR5}}}{\sqrt{\sigma_{d,\text{ABJ}}^2 + \sigma_{d,\text{DR5}}^2}}. \qquad (5) $$

Ideally, Δ<sub>ABJ-DR5</sub> would have zero mean (no biases in one measurement versus another) and unit dispersion (consistent with the uncertainties being correctly estimated). The solid red curves are Gaussians with the desired mean and dispersion. For hot and cool dwarfs, the mean value of Δ<sub>ABJ-DR5</sub> is positive and deviates from zero, which shows that the distance bias is a significant fraction of the uncertainty, with the ABJ distances being larger than the RAVE DR5 distances. The negative Δ<sub>ABJ-DR5</sub> = -0.23 for giants indicates larger distances in DR5. The characteristic feature in all three panels is the dispersion σ < 1; this implies that the uncertainties of one or both measurements are overestimated. However, the agreement of distance estimates in Fig. 9 combined with the differences between relative uncertainties in Fig. 10 suggests that the ABJ distance uncertainties are overestimated. This could be due to the possible overestimation of TGAS uncertainties as discussed in McMillan et al. (2017).

## 5.2 TGAS proper motions

In RAVE, the UCAC4 catalogue is the last proper-motion catalogue that is independent of Gaia, and it contains the largest number of stars in RAVE DR5. So we use the RAVE–TGAS sample to compare UCAC4 proper motions with TGAS and UCAC5 proper motions. Fig. 12 displays the cumulative distribution of the uncertainties in proper motion. The upper panel is for the right-ascension component, σ<sub>μα</sub>, while the lower panel is for the declination component, σ<sub>μδ</sub>. In UCAC4 and UCAC5, very few stars have uncertainties smaller than 1 mas yr<sup>-1</sup>, whereas in TGAS more than 60 per cent of stars have smaller uncertainties in declination and only slightly fewer have smaller uncertainties in right ascension. On the other hand, less than 30 per cent of stars have UCAC5 uncertainties larger than 1 mas yr<sup>-1</sup> in either component while both TGAS and UCAC4 have longer tails of stars with larger uncertainties in right ascension. Since our sample has a mean/median σ<sub>μ</sub> ~ 1 mas yr<sup>-1</sup>, the TGAS proper motions are more reliable in the selected volume. The im-

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Cumulative distribution of uncertainties in proper motion in right ascension σ<sub>μα</sub> (top) and in declination σ<sub>μδ</sub> (bottom). The blue line is for UCAC4, the black line is for UCAC5, and the red line is for TGAS. Up to ~70 per cent of TGAS proper motions have uncertainties smaller than 1 mas yr<sup>-1</sup>. In UCAC5 and UCAC4, the same is true only for very few stars.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Comparison between relative uncertainties in proper motion in right ascension $\sigma_{\mu_{\alpha}}/\mu_{\alpha}$ (top) and declination $\sigma_{\mu_{\delta}}/\mu_{\delta}$ (bottom) in the UCAC5 and TGAS catalogues. In both panels, independent of direction, most of the stars have larger relative uncertainties in UCAC5 than in TGAS.</figcaption></figure>

provement in precision provided by TGAS proper motions rather than UCAC5 is shown in Fig. 13. Both panels indicate smaller relative proper-motion uncertainties in TGAS than in UCAC5 in that the ridge lines of the distributions are steeper than the red lines for a perfect match between catalogue values. The mean relative uncertainties in right ascension for the most precisely observed 90 per cent of the stars are of the order of 6.0 per cent in UCAC5 and 4.3 per cent in TGAS for positive proper motions and 6.6 per cent in UCAC5 and 6.3 per cent in TGAS for negative proper motions.

Since even small proper-motion differences can significantly affect the space velocities obtained, the use of TGAS/UCAC5 data over UCAC4 takes us one step closer to solving the Milky Way mode discrepancy.

## 5.3 Gaia's wobbly Galaxy

We now present the velocity fields constructed from the most accurate estimates for our sample. When using the McMillan distances, we exclude stars with an error of $\sigma_d/d > 25$ per cent to exclude uncertain velocities from our velocity fields. We join these distances to line-of-sight velocities from RAVE and TGAS proper motions to form the RAVE–TGAS–McMillan data set with 58 972 stars. The top panels of Fig. 14 show the median values of $V_{\phi}$, $V_R$, and $V_z$

as a function of Galactic position for the RAVE–TGAS–McMillan set. The bottom panels show the corresponding velocity fields for the RAVE–TGAS–ABJ set. Since the ABJ distance uncertainties appear to be overestimated (Section 5.1), no distance cut has been applied.

Both the $V_{\phi}$ velocity fields shown in the left-hand panels of Fig. 14 agree with those computed in Fig. 1 using other data sets and proper motions. Within 500 pc, stars both above and below the disc mid-plane exhibit velocities $V_{\phi} \sim 220$ km s$^{-1}$, which decrease with z. At z ∼ 1 kpc and despite the smaller observed volume, the azimuthal velocity in the bottom-left panel reaches $V_{\phi} \sim 170$ km s$^{-1}$, similar to the values observed with Tycho-2 in Fig. 1.

In the middle panels of Fig. 14, both maps of median $V_R$ differ from the maps obtained with the other proper-motion catalogues (see Fig. 2). At $0.5 < z < 1$ kpc, a negative radial-velocity gradient is evident above as well as below the mid-plane. When the McMillan distance is used, the peak gradient below the plane is $\partial V_R/\partial R = -7.01 \pm 0.61$ km s$^{-1}$ kpc$^{-1}$, consistent with previous values, while above the plane the peak gradient is $\partial V_R/\partial R = -9.42 \pm 1.77$ km s$^{-1}$ kpc$^{-1}$. The ABJ distances yield similar values: $\partial V_R/\partial R = -6.01 \pm 0.62$ km s$^{-1}$ kpc$^{-1}$ below the plane and $\partial V_R/\partial R = -9.24 \pm 2.28$ km s$^{-1}$ kpc$^{-1}$ above it. Both maps show similar radial-velocity structure, with a negative gradient in both hemispheres inside $R_0$ and a positive gradient at $R > R_0$. This radial structure could be related to the simulations of Faure et al. (2014) and the analytical results of Monari et al. (2016), who studied the response of stars to a stable spiral perturbation. They found that the mean value of $V_R$ is negative within the arms and positive in the interarm region.

As a comparison to the radial gradient reported by Siebert et al. (2011b), Fig. 15 shows the projection of line-of-sight velocity in the direction of the Galactic Centre ($|z| < 5^\circ$) and anticentre ($175^\circ < l < 185^\circ$) for the RAVE–TGAS–McMillan sample, the RAVE DR5 full sample, and the results obtained in their fig. 3. We compute the mean line-of-sight velocity in bins of 0.2 kpc with a minimum of 50 stars per bin. Each solid curve represents the expected velocities of a thin-disc model with the radial-velocity gradients as indicated. The dashed curve represents a thick-disc model with zero gradient (see Siebert et al. 2011b for further details). As can be seen, both our samples are consistent within $1\sigma$ with their results.

The right-hand panels of Fig. 14 display maps of median vertical velocities. Inside $R_0$, we see a breathing mode with positive velocities above the plane and negative below; however, outside $R_0$, we see negative velocities above and below the plane, corresponding to a bending mode. This combination of breathing and bending modes is inverted relative to the combination discussed in Section 3 on the basis of the UCAC4 proper motions, and it also differs from the pattern found by W13. The mode change we observe at $R \approx R_0$ could be attributed to the Galactic warp as mentioned before, but a deeper analysis is required to understand this effect. By comparing both right-hand panels of Fig. 14 with the pattern obtained by Williams/Kordopatis, it is clear that in both our samples $V_{los}$ and $V_b$ contribute equally, whereas the $V_z$ pattern obtained by Williams/Kordopatis is consistent with a dominant $V_b$ contribution (see Fig. 6).

As in Section 4, we now study the effects of introducing errors to our sample and compare these to the results obtained with Galaxia using the same data selection. Fig. 16 shows this comparison. Although the RAVE–TGAS–McMillan sample exhibits vertical structure and Galaxia is axisymmetric, by adding a random error in distance of $\sigma_d/d = 0.8$, we increase the
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Maps of median values of each component of Galactocentric velocity. All panels employ RAVE line-of-sight velocities and TGAS proper motions. The top panels employ McMillan et al. (2017) distances with $\sigma_d/d < 25$ per cent, while the bottom panels employ the distances derived by Astraatmadja & Bailer-Jones (2016) using the Milky Way prior. In contrast to Fig. 2, the middle panels for $V_R$ display a velocity gradient above the plane. The right-hand panels for $V_z$ show the signatures of a breathing mode perturbation inside and a bending mode outside $R_0$.</figcaption></figure>

dependence of $V_z$ on $V_b$ and thus reproduce a velocity pattern similar to that found by Williams/Kordopatis. This experiment suggests that the different relative importance between the $V_{los}$ and $V_b$ patterns in the RAVE–TGAS–McMillan sample and the ones obtained by Williams/Kordopatis could arise from better distance estimates for our sample.

To study further the combination of breathing and bending modes, we now divide our subsamples into different radial bins from $R=6$ to 9 kpc in 0.5 kpc intervals and compute the median $V_z$ as a function of $z$ in bins of 0.3 kpc. Each bin contains a minimum of 50 stars.

From left to right, Fig. 17 shows the results obtained with (i) axisymmetric Galaxia model including random errors of $\sigma_d/d = 0.8$ and $\sigma_{\mu b} = 1$ mas yr$^{-1}$, (ii) UCAC4 proper motions and RAVE DR5 distances, (iii) UCAC5 proper motions and McMillan distances cut to $\sigma_d/d < 25$ per cent, (iv) the RAVE–TGAS–McMillan sample cut to $\sigma_d/d < 25$ per cent, and (v) the RAVE–TGAS–ABJ sample.

The leftmost column shows the pattern obtained from the mock RAVE sample including random errors in distance of $\sigma_d/d = 0.8$ and proper motion of $\sigma_{\mu b} = 1$ mas yr$^{-1}$ to match the errors in the RAVE–TGAS–McMillan sample. The obtained pattern in the Galaxia model shows amplitudes that are similar to the data; how-

ever, the signature of the pattern does not match any of the structures in our samples. Therefore, the patterns observed with the RAVE–TGAS subsamples may in fact be due to non-axisymmetries in the Milky Way disc.

In all RAVE–TGAS subsamples (third, fourth, and fifth columns), the top panels, with distances $6 < R < 7.5$ kpc, show positive vertical velocities in the northern Galactic hemisphere and negative in the Southern hemisphere. Thus, the stellar motions agree closely with those expected for a breathing mode. This is consistent with the velocity trend observed in fig. 4 of Widrow et al. (2012). However, the bottom panels covering the distance range $7.5 < R < 9$ kpc are mainly dominated by negative values of $V_z$ above and below the plane, corresponding to a bending mode. This signal is clear in the last column, based on ABJ distances, but weak in the other two columns. From top to bottom, we see a gradual change from a breathing to a bending mode.

In the second column of Fig. 17, we see with weak significance the opposite behaviour: a gradual change from a bending mode at distances $6 < R < 7$ kpc to a breathing mode at distances $7 < R < 9$ kpc. Given the weakness of this signal and the inferior quality of both the distances and the proper motions on which the leftmost column
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Mean of the component of line-of-sight velocity within the plane for stars that lie close to the direction to either the Galactic Centre ($|l| < 5°$) or anticentre ($175° < l < 185°$) versus projected position along the centre-anticentre line. The mean is computed in bins of 0.2 kpc with a minimum of 50 stars per bin from the RAVE–TGAS–McMillan sample and the RAVE DR5 full sample compared to the results obtained by Siebert et al. (2011b). The solid curves represent a thin disc with the indicated radial-velocity gradient, while a thick disc with zero motion is represented by the dashed curve. All three samples agree within their uncertainties.</figcaption></figure>

is based, we discount this combination of breathing and bending modes.

In contrast to previous results suggesting a breathing mode perturbation in the extended solar neighbourhood, our analysis supports a combination of breathing mode inside $R_0$ and bending mode outside $R_0$. Inwards of the solar radius, a bar and/or spiral perturbations could have induced the breathing mode; alternatively, this could be due to a satellite perturbation. The expected vertical varia-

tions due to the spiral structure are of the order of ≈4 km s⁻¹ (Faure et al. 2014) and just a small contribution due to the bar 0.5 km s⁻¹ (Monari et al. 2015); thus, the vertical amplitude we observe of ≈10 km s⁻¹ suggests that the vertical structure at $R < R_0$ may not have a purely internal origin but arise from the combination of internal and external mechanisms. On the other hand, the bending mode outside the solar radius is consistent with an external perturbation, for example caused by the Sagittarius dwarf galaxy or a dark matter subhalo passing through the disc. Gómez et al. (2013) studied a high-resolution simulation of the interaction between the Milky Way and Sgr and showed that the amplitude of the induced mean vertical velocity is ~8 km s⁻¹. Using a similar simulated setup, D’Onghia et al. (2016) reported vertical streaming motions of 10–20 km s⁻¹. Both of these are consistent with our findings. The structure we find is likely related to the Galactic warp extending outside the solar radius.

The possible combination of breathing and bending mode could thus be seen as a superposition of waves existing simultaneously in the Milky Way disc. However, further modelling is needed to understand the individual effects that create the breathing–bending mode perturbation. This will be the subject of a forthcoming paper.

## 6 SUMMARY AND CONCLUSIONS

We have used RAVE data combined with the TGAS catalogue to study the evolution of median velocities with Galactic radius and position relative to the plane. We studied two main samples: (i) a sample obtained from RAVE DR5 cross-matched with the Tycho-2, PPMXL, and UCAC4 catalogues (RAVE DR5 sample); (ii) samples comprising stars in RAVE DR5 that have TGAS astrometry using distances either inferred from TGAS trigonometric parallaxes by Astraatmadja & Bailer-Jones (2016) or distances obtained by McMillan et al. (2017) by combining TGAS parallaxes with spectrophotometric data. In agreement with recent studies, we identified asymmetries in $V_R$ and $V_z$ and less pronounced asymmetries in $V_\phi$.

<figure><img src="image_17.png" /><figcaption>**Figure 16.** Comparison of $V_z$ from the RAVE–TGAS–McMillan sample (left) and the RAVE-like sample from Galaxia (right). Similar to Fig. 7, the panels display the $V_z$ pattern obtained by adding a random error in distance with dispersion $\sigma_d/d = 0.8$. The result is a pattern similar to the one obtained with the Williams/Kordopatis sample.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure 17.** Median values of $V_z$ as a function of $z$ for different $R$ ranges for the axisymmetric Galaxia model with random errors of $\sigma_d/d = 0.8$ and $\sigma_{\mu_b} = 1$ mas yr$^{-1}$ (first column) and using proper motions from UCAC4 with RAVE DR5 distances (second column); UCAC5 with McMillan distances cut on $\sigma_d/d < 25$ per cent (third column), TGAS with McMillan distances cut on $\sigma_d/d < 25$ per cent (fourth column); and TGAS with ABJ distances (fifth column). Overlapped bins are computed every 0.3 kpc with a minimum of 50 stars per bin. The error bars correspond to the standard error of the median. The upper panels in the three right columns (which all benefit from TGAS proper motions) exhibit far from the plane positive velocities in the Northern hemisphere and negative in the Southern hemisphere, the signature of a breathing mode. The lower panels in these columns show the signature of a bending mode, even parity in the vertical velocity distribution. By contrast, the lower panels for UCAC4 show the signature of a breathing mode.</figcaption></figure>
Independently of the proper-motion catalogue used, we observed a somewhat symmetrical decline in $V_\phi$ with increasing $|z|$.

At distances $|z| < 1$ kpc, we confirmed the gradient in Galactocentric radial velocity $V_R$ previously observed. In our RAVE DR5 sample, using PPMXL proper motions, we found the largest gradient in the southern Galactic hemisphere to be $\partial V_R/\partial R = -5.99 \pm 0.81 \text{ km s}^{-1} \text{ kpc}^{-1}$, consistent with the gradient reported by Siebert et al. (2011b). In the Northern hemisphere at $R < R_0$, we identified similar gradients, with the largest value being $\partial V_R/\partial R = -8.53 \pm 1.82 \text{ km s}^{-1} \text{ kpc}^{-1}$ obtained with the Tycho-2 proper motions. Using our more accurate RAVE–TGAS–McMillan sample, we measured $\partial V_R/\partial R = -7.01 \pm 0.61 \text{ km s}^{-1} \text{ kpc}^{-1}$ in the southern Galactic hemisphere and $\partial V_R/\partial R = -9.42 \pm 1.77 \text{ km s}^{-1} \text{ kpc}^{-1}$ in the Northern hemisphere. As seen in Fig. 14, middle panels, the gradient seems to reverse sign at $R > R_0$, consistent with the expected effect of spiral structure (Faure et al. 2014; Monari et al. 2016; this, however, does not mean that the bar plays no role in shaping the gradient). In both the RAVE DR5 and RAVE–TGAS samples, we used the projection of the line-of-sight velocity on the Galactic plane in the direction of the Galactic Centre and anticentre to compare our results with the radial gradient previously obtained. Our samples proved to be consistent within $1\sigma$. While our sample and those of many previous works are a mixture of stars with different ages and chemistry, the radial gradient can vary strongly for stellar populations with narrow ranges in ages and/or metallicity (Wojno et al., in preparation).

The $V_z$ velocity fields display a more complex structure, the origin of which is under debate. Spiral arms (Faure et al. 2014; Monari et al. 2016), and to a lesser extent the bar (Monari et al. 2015), naturally induce a pattern with odd parity in the $V_z$ distribution of stars with respect to $z$, associated with a breathing mode. In contrast, bending modes (even parity in the $V_z$ distribution) are attributed mainly to external perturbations such as a satellite galaxy or dark matter subhalo crossing the Galactic plane. Widrow et al. (2014), however, showed that a passing satellite galaxy could produce both bending and breathing modes depending on the vertical velocity of the satellite.

We have shown that maps of $V_z$ depend strongly on the adopted proper motions and distances. W13 and Kordopatis et al. (2013b) employed a compilation of proper motions and found that inside the solar radius, there was upward motion above the plane and downward motion below it. Outside the solar radius, this velocity pattern reversed. This pattern is consistent with that of a breathing mode.

Using our RAVE DR5 sample, with more accurate proper motions from the UCAC4 catalogue, we confirmed the observed breathing mode outside the solar radius. However, inside $R_0$, we found upward motion both above and below the plane, consistent with a bending mode. Thus, using the UCAC4 proper motions with RAVE DR5 changes the breathing mode observed by Williams/Kordopatis to a possible combination of bending and breathing modes (compare the right-hand panels of Fig. 4).

Our most accurate velocity field, obtained with TGAS proper motions and McMillan distances, which exploits both TGAS parallaxes and the IRFM, supports a combination of bending and breathing modes.

After studying the components of the vertical velocity, we found the Williams/Kordopatis breathing mode to be due to the contribution of the transverse velocity component dominating that of line-of-sight velocities (Fig. 6). The structure resulting from the transverse velocity alone depends mainly on the proper motion, while the distance increases its amplitude. By adding artificial dis-

tance errors to our RAVE–TGAS sample and to an axisymmetric mock RAVE-like sample from Galaxia, we were able to increase the transverse velocity in both samples and artificially reproduce the breathing mode observed by W13 (see Fig. 16).

The combination of modes found in the RAVE–TGAS sample is the inverse of that found with UCAC4 proper motions. Inside $R_0$, we identified patterns that appear consistent with a breathing mode and outside $R_0$ they seem to agree with a bending mode. Similar structure is found with UCAC5 proper motions (see Fig. 17). This combination of breathing and bending modes could be seen as a superposition of waves existing simultaneously in the Milky Way disc. Unlike the pattern found with UCAC4, the TGAS/UCAC5 pattern makes physical sense; inwards of the solar radius, a bar and/or spiral perturbations could induce the observed breathing mode, while outside, bending modes would likely be generated by external perturbations such as a passing satellite galaxy or a dark matter subhalo. Further modelling is needed to understand what combination of perturbations to the Milky Way disc can induce a bending/breathing mode outside/inside the solar circle.

To confirm our results that the Milky Way exhibits both bending and breathing modes, a larger disc volume must be probed. Compared to Gaia DR1, Gaia DR2 will cover a significantly larger volume of the Milky Way disc and improve significantly the data systemsatics. For stars within the RAVE magnitude, Gaia DR2 will have a preliminary parallax error of $\sigma_{\omega} = 0.03$ mas (Katz & Brown 2017) and proper-motion uncertainties of $\sigma_{\mu} \approx 0.04$ mas yr$^{-1}$ (Marchetti et al. 2017). This will improve the median $\sigma_{V_z}$ and $\sigma_{V_R}$ by a factor of 6. Thus, using data from Gaia DR2 may solve the question whether the Milky Way is still just breathing.

## ACKNOWLEDGEMENTS

IC is grateful to Ralf-Dieter Scholz for the many discussions about proper motions, to Mary E. K. Williams for comments that improved the quality and clarity of this paper, to Lawrence M. Widrow and Matthew H. Chequers for useful discussions about breathing and bending modes, to Jorrit Hagen for discussions on the effects of distance errors, and to Haifeng Wang for discussions on velocity asymmetries. Funding for RAVE has been provided by the Australian Astronomical Observatory; the Leibniz-Institut für Astrophysik Potsdam (AIP); the Australian National University; the Australian Research Council; the French National Research Agency; the German Research Foundation (SPP 1177 and SFB 881); the European Research Council (ERC-StG 240271 Galactica); the Istituto Nazionale di Astrofisica at Padova; the Johns Hopkins University; the National Science Foundation of the USA (AST-0908326); the W. M. Keck foundation; the Macquarie University; the Netherlands Research School for Astronomy; the Natural Sciences and Engineering Research Council of Canada; the Slovenian Research Agency (research core funding No. P1-0188); the Swiss National Science Foundation; the Science & Technology Facilities Council of the UK; Opticon; Strasbourg Observatory; and the Universities of Groningen, Heidelberg, and Sydney. The RAVE website is https://www.rave-survey.org. EKG acknowledges support by Sonderforschungsbereich 'The Milky Way System' (SFB 881) of the German Research Foundation (DFG), particularly through subproject A5.

This work has made use of data from the European Space Agency (ESA) mission Gaia (http://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, http://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 Multi-
lateral Agreement. This work was supported by the European
Research Council under the European Union’s Seventh Framework
Programme (FP7/2007-2013)/ERC grant agreement no. 321067.

REFERENCES

Aihara H. et al., 2011, ApJS, 193, 29

Anguiano B., Rebassa-Mansergas A., García-Berro E., Torres S., Freeman
K. C., Zwitter T., 2017, MNRAS, 469, 2102

Antoja T., Figueras F., Fernández D., Torra J., 2008, A&A, 490, 135

Antoja T. et al., 2014, A&A, 563, A60

Asiain R., Figueras F., Torra J., Chen B., 1999, A&A, 341, 427

Astraatmadja T. L., Bailer-Jones C. A. L., 2016, ApJ, 833, 119

Bailer-Jones C. A. L., 2015, PASP, 127, 994

Bertelli G., Bressan A., Chiosi C., Fagotto F., Nasi E., 1994, A&AS, 106, 275

Binney J., Merrifield M., 1998, Princeton Series in Astrophysics, Galactic Astronomy. Princeton Univ. Press, Princeton, NJ

Binney J. et al., 2014, MNRAS, 437, 351

Blackwell D. E., Shallis M. J., Selby M. J., 1979, MNRAS, 188, 847

Carlin J. L. et al., 2013, ApJ, 777, L5

Casagrande L., Ramírez I., Meléndez J., Bessell M., Asplund M., 2010, A&A, 512, A54

Chequers M. H., Widrow L. M., 2017, MNRAS, 472, 2751

Chereul E., Creze M., Bienayme O., 1998, A&A, 340, 384

Cui X.-Q. et al., 2012, Res. Astron. Astrophys., 12, 1197

D'Onghia E., Madau P., Vera-Ciro C., Quillen A., Hernquist L., 2016, ApJ, 823, 4

de la Vega A., Quillen A. C., Carlin J. L., Chakrabarti S., D'Onghia E., 2015, MNRAS, 454, 933

Dehnen W., 1998, AJ, 115, 2384

Dehnen W., 2000, AJ, 119, 800

Dehnen W., 2001, in Deiters S., Fuchs B., Just A., Spurzem R., Wielen R., eds, ASP Conf. Ser. Vol. 228, Dynamics of Star Clusters and the Milky Way. Astron. Soc. Pac., San Francisco, p. 413

Drimmel R., Spergel D. N., 2001, ApJ, 556, 181

Famaey B., Jorissen A., Luri X., Mayor M., Udry S., Dejonghe H., Turon C., 2005, A&A, 430, 165

Famaey B., Siebert A., Jorissen A., 2008, A&A, 483, 453

Faure C., Siebert A., Famaey B., 2014, MNRAS, 440, 2564

Fux R., 2001, A&A, 373, 511

Gaia Collaboration et al., 2016a, A&A, 595, A1

Gaia Collaboration et al., 2016b, A&A, 595, A2

Gnedin O. Y., Brown W. R., Geller M. J., Kenyon S. J., 2010, ApJ, 720, L108

Gómez F. A., Minchev I., Villalobos Á., O’Shea B. W., Williams M. E. K., 2012a, MNRAS, 419, 2163

Gómez F. A. et al., 2012b, MNRAS, 423, 3727

Gómez F. A., Minchev I., O’Shea B. W., Beers T. C., Bullock J. S., Purcell C. W., 2013, MNRAS, 429, 159

Gómez F. A., White S. D. M., Marinacci F., Slater C. T., Grand R. J. J., Springel V., Pakmor R., 2016, MNRAS, 456, 2779

Groenewegen M. A. T., 2013, A&A, 550, A70

Hög E. et al., 2000, A&A, 355, L27

Ibata R. A., Gilmore G., Irwin M. J., 1994, Nature, 370, 194

Johnson D. R. H., Soderblom D. R., 1987, AJ, 93, 864

Kalberla P. M. W., Dedes L., Kerp J., Haud U., 2007, A&A, 469, 511

Katz D., Brown A. G. A., 2017, in Reylé C., Di Matteo P., Herpin F., Lagadec E., Lançon A., Meliani Z., Royer F., eds, SF2A: Proc. Annu. Meeting French Soci. Astron. Astrophys., p. 259

Kordopatis G. et al., 2013a, AJ, 146, 134

Kordopatis G. et al., 2013b, MNRAS, 436, 3231

Kunder A. et al., 2017, AJ, 153, 75

Levine E. S., Blitz L., Heiles C., 2006, ApJ, 643, 881

Lindegren L. et al., 2016, A&A, 595, A4

McMillan P. J., 2013, MNRAS, 430, 3276

McMillan P. J., Binney J. J., 2010, MNRAS, 402, 934

McMillan P. J. et al., 2017, MNRAS, preprint (arXiv:1707.04554)

Marchetti T., Contigiani O., Rossi E. M., Albert J. G., Brown A. G. A.,
Sesana A., 2017, MNRAS, preprint (arXiv:1711.11397)

Marigo P., Girardi L., Bressan A., Groenewegen M. A. T., Silva L., Granato
G. L., 2008, A&A, 482, 883

Minchev I., Quillen A. C., Williams M., Freeman K. C., Nordhaus J., Siebert
A., Bienaymé O., 2009, MNRAS, 396, L56

Minchev I., Boily C., Siebert A., Bienayme O., 2010, MNRAS, 407, 2122

Monari G., Helmi A., Antoja T., Steinmetz M., 2014, A&A, 569, A69

Monari G., Famaey B., Siebert A., 2015, MNRAS, 452, 747

Monari G., Famaey B., Siebert A., 2016, MNRAS, 457, 2569

Monari G., Kawata D., Hunt J. A. S., Famaey B., 2017, MNRAS, 466, L113

Perryman M. A. C., ESA, eds, 1997, in ESA SP-1200: The Hipparcos and Tycho Catalogues. ESA, Noordwijk, p. L49

Poggio E., Drimmel R., Smart R. L., Spagna A., Lattanzi M. G., 2017, A&A, 601, A115

Pompéia L. et al., 2011, MNRAS, 415, 1138

Quillen A. C., Minchev I., 2005, AJ, 130, 576

Reid M. J., 1993, ARA&A, 31, 345

Reid M. J., Brunthaler A., 2004, ApJ, 616, 872

Reid M. J. et al., 2009, ApJ, 700, 137

Robin A. C., Reylé C., Derrière S., Picaud S., 2003, A&A, 409, 523

Robin A. C., Reylé C., Marshall D. J., 2008, Astron. Nachr., 329, 1012

Roeser S., Demleitner M., Schilbach E., 2010, AJ, 139, 2440

Schönrich R., Binney J., Dehnen W., 2010, MNRAS, 403, 1829

Sharma S., Bland-Hawthorn J., Johnston K. V., Binney J., 2011, ApJ, 730, 3

Siebert A. et al., 2011a, AJ, 141, 187

Siebert A. et al., 2011b, MNRAS, 412, 2026

Siebert A. et al., 2012, MNRAS, 425, 2335

Steinmetz M. et al., 2006, AJ, 132, 1645

van Leeuwen F. ed. 2007, Astrophysics and Space Science Library, Vol. 350, Hipparcos, the New Reduction of the Raw Data. Springer-Verlag, Berlin

Vickers J. J., Röser S., Grebel E. K., 2016, ApJ, 151, 99

Widrow L. M., Gardner S., Yanny B., Dodelson S., Chen H.-Y., 2012, ApJ, 750, L41

Widrow L. M., Barber J., Chequers M. H., Cheng E., 2014, MNRAS, 440, 1971

Williams M. E. K., Freeman K. C., Helmi A., RAVE Collaboration, 2009, in Andersen J., Nordström B., Bland-Hawthorn J., eds, Proc. IAU Symp. 254, The Galaxy Disk in Cosmological Context. Cambridge Univ. Press, Cambridge, p. 139

Williams M. E. K. et al., 2013, MNRAS, 436, 101 (W13)

Wojno J. et al., 2017, MNRAS, 468, 3368

Xu Y., Newberg H. J., Carlin J. L., Liu C., Deng L., Li J., Schönrich R., Yanny B., 2015, ApJ, 801, 105

Yanny B. et al., 2009, AJ, 137, 4377

Zacharias N., Finch C. T., Girard T. M., Henden A., Bartlett J. L., Monet D. G., Zacharias M. I., 2013, AJ, 145, 44

Zacharias N., Finch C., Frouard J., 2017, AJ, 153, 166

<footnote n="1">¹Leibniz Institut für Astrophysik Potsdam (AIP), An der Sterwarte 16, D-14482 Potsdam, Germany</footnote>

<footnote n="2">²Laboratoire Lagrange, Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Bd de l'Observatoire, CS 34229, F-06304 Nice cedex 4, France</footnote>

<footnote n="3">³Rudolf Peierls Centre for Theoretical Physics, Keble Road, Oxford OX1 3NP, UK</footnote>

<footnote n="4">⁴Observatoire astronomique de Strasbourg, Université de Strasbourg, CNRS, UMR 7550, 11 rue de l'Université, F-67000 Strasbourg, France</footnote>

<footnote n="5">⁵Sydney Institute for Astronomy, School of Physics A28, University of Sydney, NSW 2006, Australia</footnote>

<footnote n="6">⁶Research School of Astronomy & Astrophysics, Australian National University, Cotter Rd., Weston, ACT 2611, Australia</footnote>
This paper has been typeset from a TEX/LATEX file prepared by the author.

<footnote n="7">⁷Institute of Astronomy, University of Cambridge, Madingley Road, Cam-
bridge CB3 0HA, UK</footnote>

<footnote n="8">⁸E.A. Milne Centre for Astrophysics, University of Hull, Hull HU6 7RX, UK</footnote>

<footnote n="9">⁹Astronomisches Rechen-Institut, Zentrum für Astronomie der Universität Heidelberg, Mönchhofstr. 12–14, D-69120 Heidelberg, Germany</footnote>

<footnote n="10">¹⁰Kapteyn Astronomical Institute, University of Groningen, PO Box 800, NL-9700 AV Groningen, the Netherlands</footnote>

<footnote n="11">¹¹Saint Martin's University, 5000 Abbey Way SE, Lacey, WA 98503, USA</footnote>

<footnote n="12">¹²Lund Observatory, Department of Astronomy and Theoretical Physics, Lund University, Box 43, SE-22100 Lund, Sweden</footnote>

<footnote n="13">¹³The Oskar Klein Centre for Cosmoparticle Physics, Department of Physics, Stockholm University, AlbaNova, SE-10691 Stockholm, Sweden</footnote>

<footnote n="14">¹⁴INAF National Institute of Astrophysics, Astronomical Observatory of Padova, I-36012 Asiago, Italy</footnote>

<footnote n="15">¹⁵Department of Physics and Astronomy. University of Victoria. Victoria, BC V8P5C2, Canada</footnote>

<footnote n="16">¹⁶Department of Physics, University of Hong Kong, Hong Kong SAR, China</footnote>

<footnote n="17">¹⁷The Laboratory for Space Research, University of Hong kong, Hong Kong SAR, China</footnote>

<footnote n="18">¹⁸Department of Physics and Astronomy, Macquarie University, Sydney, NSW 2109, Australia</footnote>

<footnote n="19">¹⁹Western Sydney University, Locked Bag 1797, Penrith South, NSW 2751, Australia</footnote>

<footnote n="20">²⁰Mullard Space Science Laboratory, University College London, Holmbury St Mary, Dorking RH5 6NT, UK</footnote>

<footnote n="21">²¹Australian Astronomical Observatory, PO Box 915, North Ryde, NSW 1670, Australia</footnote>

<footnote n="22">²²Department of Physics and Astronomy, Johns Hopkins University, 3400 N. Charles St, Baltimore, MD 21218, USA</footnote>

<footnote n="23">²³Faculty of Mathematics and Physics, University of Ljubljana, Jadranska 19, 1000 Ljubljana, Slovenia</footnote>