The wobbly Galaxy: kinematics north and south with RAVE red clump giants

As we are working at large distances from the solar position, we use cylindrical co-ordinates for the most part, with $V_{R},\,V_{\phi}$ and $V_{Z}$ defined as positive with increasing $R$, $\phi$ and $Z$, with the latter towards the North Galactic Pole (NGP). We also use a right-handed Cartesian co-ordinate system with $X$ increasing towards the Galactic centre, $Y$ in the direction of rotation and $Z$ again positive towards the NGP. The Galactic centre is located at $(X,Y,Z)=(R_{\odot},\,0,\,0)$. Space velocities, $UVW$, are defined in this system, with $U$ positive towards the Galactic centre.

To aid comparison with the predictions of Galaxia (see Section 2.3), we use mostly the parameter values that were used to make the model for which Galaxia makes predictions. Thus, the motion of the Sun with respect to the LSR is taken from Schönrich, Binney & Dehnen 2010, namely $U_{\odot}=11.1,\ V_{\odot}=12.24,W_{\odot}=7.25\,\ \mathrm{kms}^{-1}$. The LSR is assumed to be on a circular orbit with circular velocity $V_{\mathrm{circ}}=226.84\,\ \mathrm{kms}^{-1}$. Finally, we take $R_{\odot}=8\,\mathrm{kpc}$. The only deviation from Galaxia’s values is that they assumed that the Sun is located at $Z=+15\,\textrm{pc}$, while we assume $Z=0\,\textrm{pc}$.

S11 explored the variation of the observed $V_{R}$ gradient on the values of $V_{\mathrm{circ}}$ and $R_{\odot}$, finding that changing these parameters between variously accepted values could reduce but not eliminate the observed gradient in $V_{R}$. We do not explore this in detail in this paper but note that the trends that we observe are similarly affected by the Galactic parameters; amplitudes are changed but the qualitative trends remain the same.

The wide-field RAVE (RAdial Velocity Experiment) survey measures primarily line-of-sight velocities and additionally stellar parameters, metallicities and abundance ratios of stars in the solar neighbourhood (Steinmetz et al., 2006; Zwitter et al., 2008; Siebert et al., 2011; Boeche et al., 2011). RAVE’s input catalog is magnitude limited ($8<I<13$) and thus creates a sample with no kinematic biases. To the end of 2012 RAVE had collected more than 550,000 spectra with a median error of $1.2\,\ \mathrm{kms}^{-1}$ (Siebert et al., 2011b).

We use the internal release of RAVE from October 2011 which contains 434,807 RVs and utilizes the revised stellar parameter determination (see the DR3 paper, Siebert et al. 2011b for details). We applied a series of cuts to the data. Firstly, those stars flagged by the Matijevič et al. 2012 automated spectra classification code as having peculiar spectra were excluded. This removes most spectroscopic binaries, chromospherically active and carbon stars, spectra with continuum abnormalities and other unusual spectra. Further cuts restricted our sample to stars with signal-to-noise ratio $STN>20$ (STN is calculated from the observed spectrum alone with residuals from smoothing, see Section 2.2 of the DR3), Tonry & Davis 1979 cross-correlation coefficient $R>5$, $|\mu_{\alpha},\ \mu_{\delta}|<400\,\mathrm{masyr}^{-1}$, $e\mu_{\alpha},\ e\mu_{\delta}<20\,\mathrm{masyr}^{-1}$, $\mathrm{RV}<600\,\ \mathrm{kms}^{-1}$ (see Section 4.4) and stars whose SpectraQualityFLAG is null. Where there are repeat observations, we randomly select one observation for each star. With this cleaning, the data set has 293,273 unique stars with stellar parameters from which to select the red clump.

In some fields at $|b|<25^{\circ}$ the RAVE selection function includes a colour cut $J-K>0.5$ with the object of favouring giants. We imposed this colour cut throughout this region to facilitate the comparison to predictions by Galaxia. The selection of red clump stars is unaffected by this cut. Additional data cuts and selections were performed later in the analysis and are broadly a) the selection of red-clump giants (Section 3.3), b) removal of stars with large extinction/reddening (Section 4.1.5), and removal of data bins with large errors in measurements of mean velocity (Section 4.4).

As a sample with alternative distance determinations, we also used an internal release with stellar parameters produced by the pipeline that was used for the third Data Release (VDR3) and the method of Zwitter et al. 2008 from September 2011 with 334,409 objects. We cleaned the sample as above, leaving us with 301,298 stars.

A newer analysis of RAVE data is presented in the 4th RAVE data release, Kordopatis et al. (in preparation). This applies an updated version of the Kordopatis et al. 2011 pipeline to RAVE spectra. Selecting red clump stars using the stellar parameters from this pipeline does not produce significantly different results to those selected using the VDR3 pipeline above; the conclusions of this paper are unaffected by the pipeline used for the red-clump selection.

In Williams et al. 2011 we introduced the use of the Galaxy modelling code Galaxia (Sharma et al., 2011) to investigate the statistical significance of the new Aquarius stream found with RAVE. In this study we use Galaxia both to provide predictions with which to compare our results, and to investigate the effects of contamination of the red-clump sample by stars that are making their first ascent of the giant branch (see Section 4.1.1). Galaxia enables us to disentangle real effects from artifacts of our methodology, and further to understand the population we are examining.

Based on the Besançon Galaxy model, the Galaxia code creates a synthetic catalog of stars for a given model of the Milky Way. It offers several improvements over the Besançon model which increase its utility in modelling a large-scale survey like RAVE, the most significant of which is the ability to create a continuous distribution across the sky instead of discrete sample points. The elements of the Galactic model are a star formation rate, age-velocity relation, initial mass function and density profiles of the Galactic components (thin and thick disk, smooth spheroid, bulge and dark halo). The parameters and functional forms of these components are summarized in Table 1 of Sharma et al. 2011.

We follow a similar methodology to that described in Williams et al. 2011 to generate a Galaxia model of the RAVE sample, from which we select a red clump sample following the selection criteria applied to the real data. A full catalog was generated over the area specified by $0<l<360$, $\delta<2^{\circ}$ and $9<I<13$, with no under-sampling. As described in Section 3.2, the RAVE data were then divided into three signal-to-noise regimes and a sample drawn from the Galaxia model so that, for each regime, the $I$-band distribution was matched to that of the RAVE sample in $5^{\circ}\times 5^{\circ}$ squares. The Galaxia $I$-band was generated after correcting for extinction.

The $K$-band magnitude of the red clump, while being relatively unaffected by extinction, has also been shown observationally to be only weakly dependent on metallicity and age (Pietrzyński, Gieren & Udalski, 2003; Alves, 2000), so a single magnitude is usually assigned for all stars in the red clump. Studies such as Grocholski & Sarajedini 2002 and van Helshoecht & Groenewegen 2007 have shown that there is some dependence on metallicity and age which can be accounted for by the theoretical model of Salaris & Girardi 2002.

If models correctly predict the systematic dependence of $M_{K}$ of the red clump on metallicity and age, this would have implications for our study of kinematics in the solar suburb, for with increasing distance above the plane the metallicity will decrease and the population become older. From Burnett et al. 2011 we estimate that the population means change from $\mathrm{[M/H]}\sim 0$, $\mathrm{Age}=4\,\mathrm{Gyr}$ at $Z=0\,\mathrm{kpc}$ to $\mathrm{[M/H]}\sim-0.6$, $\mathrm{Age}=10\,\mathrm{Gyr}$ at $Z=|4|\,\mathrm{kpc}$. According to Salaris & Girardi 2002, these age and metallicity changes will change the $K$-band absolute magnitude from $M_{K_{s}}\textrm{(RC)}=-1.64$¹¹1Here $K$ is used to denote $K$-band magnitudes in the Bessell & Brett 1989 system, while $K_{s}$ is used to denote 2MASS values. The relations from Carpenter 2003 were used to convert between the two photometric systems. Note though that $J-K_{S}$ is shortened to $J-K$ in Section 3.2 and beyond. at $Z=0\,\mathrm{kpc}$ to $M_{K_{s}}\textrm{(RC)}=-1.39$ at $|Z|=4\,\mathrm{kpc}$. Hence, the distances to stars at higher $Z$ will be systematically underestimated by $\sim 10\>\mathrm{per\>cent}$.

Furthermore, there is uncertainty regarding the average value of $M_{K}$. Alves 2000 gives $M_{K}\textrm{(RC)}=-1.61\pm 0.03$ for local red clump giants with Hipparcos distances and metallicities between $-0.5\leq\mathrm{[Fe/H]}\leq 0.0$. In the 2MASS system the corresponding absolute magnitude is $M_{K_{s}}\textrm{(RC)}=-1.65\pm 0.03$. Grocholski & Sarajedini 2002 derives a similar value of $M_{K}\textrm{(RC)}=-1.61\pm 0.04$ for $-0.5\leq\mathrm{[Fe/H]}\leq 0.0$, $1.6\,\mathrm{Gyr}\leq\mathrm{Age}\leq 8\,\mathrm{Gyr}$ from 2MASS data of open clusters. Extending the sample of clusters, van Helshoecht & Groenewegen 2007 derived $M_{K}\textrm{(RC)}=-1.57\pm 0.05$ for $-0.5\leq\mathrm{[Fe/H]}\leq 0.4$, $0.3\,\mathrm{Gyr}\leq\mathrm{Age}\leq 8\,\mathrm{Gyr}$. Then using the new van Leeuwen 2007 Hipparcos parallaxes, Groenewegen 2008 found that the Hipparcos red clump giants now give $M_{K_{s}}\textrm{(RC)}=-1.54\pm 0.05$ ($M_{K}\textrm{(RC)}=-1.50\pm 0.05$) over $-0.9\leq\mathrm{[Fe/H]}\leq 0.3$, with a selection bias, whereby accurate K-magnitudes are only available for relatively few bright stars, meaning the actual value is likely brighter. The newer values hold better agreement with the theoretical results of Salaris & Girardi 2002, who derived an average value for the solar neighbourhood of $M_{K_{s}}\textrm{(RC)}=-1.58$, derived via modelling with a SFR and age-metallicity relation.

Given the disagreement over the $K$-band magnitude for the clump and possible metallicity/age variations, we investigated the use of the three normalizations for $M_{K_{s}}$ for our derivation of the RAVE red clump distances. Table 1 summarizes these values, where the Version A is the standard value from Alves 2000 and Grocholski & Sarajedini 2002, while Version B is the new value from Groenewegen 2008. Version C is derived from the theoretical models of Salaris & Girardi 2002 and attempts to take into account some of the possible systematics for lower metallicity stars, where we use the variation of age and metallicity of RAVE stars described above (i.e., $M_{K_{s}}(RC)=-1.64$ at $|Z|=0\,\mathrm{kpc}$ to $M_{K_{s}}(RC)=-1.39$ at $|Z|=4\,\mathrm{kpc}$) to develop a simple linear relation of $M_{K_{s}}$ with $Z$, where the value for each star is derived iteratively. In Section 4.1.2 we examine the effect of the different normalizations on our results.

To establish the best selection method for the red-clump stars we first tried to match the observed distribution of the selected RC stars with the prediction of Galaxia. The errors in RAVE stellar parameters decrease with STN, so the red clump is more localized at higher STN values. Further, the errors in $T_{\mathrm{eff}}$ and $\log g$ are correlated. To account for this in the modelling of the clump, we split the data into three regimes: STN values between $20<\mathrm{STN}<40$, $40<\mathrm{STN}<60$ and $60<\mathrm{STN}$. We then generated Galaxia models for each STN regime for $5^{\circ}\times 5^{\circ}$ squares, matching the I-band distribution of Galaxia to that of RAVE in bins of 0.2 mag. The three models for each STN slice of the data were added together to give an overall Galaxia model, with an STN value of 20, 40 or 60 assigned to the Galaxia ‘stars’ to indicate which STN regime they were generated from. Note that, as stated in the DR2 and DR3 releases, the DENIS I-band photometry has large errors for a significant fraction of the stars. We therefore used Equation 24 from the DR2 paper to calculate I magnitudes for those RAVE stars which do not satisfy the condition $-0.2<(I_{\mathrm{DENIS}}-J_{\mathrm{2MASS}})-(J_{\mathrm{2MASS}}-K_{\mathrm{2MASS}})<0.6$, as well as those that have 2MASS photometry but no DENIS values. For the models the reddening at infinity is matched to that of the value in Schlegel map and to convert $E(B-V)$ to extinction in different photometric bands we used the conversion factors in Table 6 of Schlegel et al. 1998.

Errors were then added to the values of $T_{\mathrm{eff}}$, $\log g$ and $J-K$ from Galaxia to match the distributions in the RAVE data. The same seed was used for the random number generator for the temperature and gravity values to account for the error covariance. Table 2 gives the standard deviations of the errors that were added to the values from Galaxia. The average $J-K$ error, $0.01$, is smaller than the expected average observational error in $J-K$, $\sqrt{eJ^{2}+eK_{S}^{2}}=0.03$. An error of $0.01$ in $J-K$ corresponds to an error in $T_{\mathrm{eff}}$ of $\sim 40\,K$ (Alonso, Arribas & Martinez-Roger, 1996). Given this small value and that, unlike $T_{\mathrm{eff}}$, the spread in $J-K$ does not change with STN, we opted to select the RAVE red clump in $(J-K,\log g)$.

Figure 1 plots the $(J-K,\log g)$ plane for the RAVE giants within $8<I<13$ in the three STN regimes along with the corresponding Galaxia predictions including errors. Here we see the clump at $J-K\sim 0.65$ and $\log g\sim 2.2$. There are some notable differences between the models and the observed distributions. First, the position of the clump in $\log g$ decreases with decreasing STN for the RAVE stars. This change is not matched by Galaxia: the predicted distribution has mean $\log g\sim 2.4$. With increasing dispersion at lower STN it elongates but does not shift. Thus, the effect seen in the RAVE data appears not to be astrophysical in origin but rather a result of the parameter determination from RAVE spectra (for a discussion on systematic trends of stellar parameter estimates with SNR see DR3). Also, there is an absence of horizontal-branch stars below $J-K<0.55$ in the RAVE data relative to the prediction by Galaxia. This latter difference does not affect our results as they are not included in the selection region.

We select the red clump as those stars within $0.55\leq J-K\leq 0.8$, $1.8\leq\log g\leq 3.0$. We do not include a STN dependence as extending the selection box further up in gravity will mean an increased fraction of sub-giants in the sample and thus erroneous distances. The red box in Figure 1 displays the selection criteria, where we see that for $20<\mathrm{STN}<40$ some of the clump is cutoff at higher values of $\log g$. Applied to our entire cleaned RAVE set, this selection yields 78,019 stars. A sample of red clump stars were also selected from the Galaxia model in the same way, with a slightly smaller number of 73,594 ‘stars’.

Figure 2 gives the distributions of the entire RAVE data set and red clump stars in $I,\ V_{\mathrm{los}}$ and overall proper motion, $\mu=\sqrt{\mu_{\alpha}^{2}+\mu_{\delta}^{2}}$. The distributions from Galaxia are also given, with an additional spread of $2\,\ \mathrm{kms}^{-1}$ and $2.7\,\mathrm{masyr}^{-1}$ added to the Galaxia line-of-sight velocity and proper motion, respectively, to simulate the observational uncertainties. These values give the average value of the error in these quantities for the RAVE catalog stars. Here we see that Galaxia reproduces quite well the basic distributions of the two observational kinematic parameters which affirms our use of it in comparisons. The small difference in the proper motion distribution could suggest that the actual errors in proper motion are slightly larger than estimated.

The non-systematic sources of errors in the kinematics can be broadly broken into (i) the contribution of random uncertainties in the measurements and (ii) the finite sample size for each volume bin. We will examine each of these in turn, before discussing any kinematic cuts that were applied to the data.

The Monte Carlo method of error propagation enables an easy tracking of error covariances, so we employ it here by generating a distribution of 100 test particles around each input value of distance, proper motion, line-of-sight velocity²²2The positions of the stars are assumed to have negligible errors. and then calculating the resulting points in $(V_{R},\,V_{\phi},\,V_{Z})$. We assume that the proper motions and line-of-sight velocities have Gaussian errors with standard deviations given by the formal errors for each star. For the distribution in distance for each star however, we generate a double-Gaussian distribution in the parallax as in Figure 3, which we then invert to derive the distance distribution. This takes into account the fact that errors in distance are due both to the intrinsic width of the RC and the mis-classification of other giants. The ratio between the RC and first-ascent giants for Galaxia model stars changes with distance, where we found a smaller amount of contaminants for $d>1\,\mathrm{kpc}$. The double-Gaussian in parallax then has dispersions $(0.05,0.26)\,\mathrm{mas}$ at a ratio of $2:1$ for $d\leq 1\,\mathrm{kpc}$ and $4:1$ for $d>1\,\mathrm{kpc}$.

In each of the bins that are used to calculate the mean velocities, there are a finite number of stars $N$, so Poisson noise contributes to the errors in the derived kinematic properties. These errors scale as $1/\sqrt{N}$. To establish the error contribution from this source we used Bootstrap case resampling with replacement: for each kinematic quantity we derived a distribution in the values by randomly resampling from the distribution of values in the bin. The variance in each value could then be calculated from the resulting distribution. Poisson errors can dominate over measurement errors in bins at large distances that contain very few stars.

To ensure that Poisson errors do not dominate our plots, we only use bins which have $N_{\textit{stars}}>50$ and only include those points that have errors in the mean of less than $5\ \mathrm{kms}^{-1}$ (see Section 4.4).

There are several ways in which to prune the data to those that are deemed more reliable. However, pruning is liable to introduce kinematic biases. We investigated the effects on the measured velocities of trimming data via a) proper motion, b) error in proper motion, c) distance errors, d) magnitude of total velocity and e) total velocity error, $eV_{\mathrm{total}}$. In general, we found that as we increase the cut-off point there is a steady increase in velocity dispersion and fluctuations in the average velocity, asymptotically approaching a value as the number of stars approaches the full sample. It is therefore difficult to justify cuts that remove a significant proportion of stars. We therefore introduced cuts that remove only the outlier values, as given in Table 3, and we do not perform a cut on $eV_{\mathrm{total}}$.

Since we seek only to follow trends with $R,\,Z$, we do not distinguish between halo, thick and thin disk stars. The cuts in total velocity therefore aim to be inclusive of halo stars; the $600\,\ \mathrm{kms}^{-1}$ limit is $\sim 3\sigma$ the halo dispersion of $213\,\ \mathrm{kms}^{-1}$ (Vallenari et al., 2006). An additional cut in $V_{\phi}$ is also introduced to limit ourselves to stars that have plausible rotation velocities.

In our analysis we bin the data in physical space, calculating the means and dispersion in each bin. A further data cut was performed post-binning. For each bin we calculate the mean $\langle V_{R}\rangle,\ \langle V_{\phi}\rangle,\ \langle V_{Z}\rangle$. For each of these we also calculate an error in the mean, given by standard error propagation as

where $q=R,\ \phi,\ Z$ and $i=1,..,N$, with $N$ the number of stars in the average. The values $e_{V_{q},i}$ are given by the MC propagation described in Section 4.2. We remove bins with large errors, i.e., $e_{\langle V_{q}\rangle}>5\ \mathrm{kms}^{-1}$. This affects only peripheral points at large distances from the Sun.

The top two panels of Figure 10 show that, as we might expect, $V_{\phi}$ is largest in the plane and decreases fairly symmetrically both above and below the plane. At given values of $Z$, $V_{\phi}$ also decreases inwards. These trends are independent of the distances used (RC or Zwitter) and of the adopted proper-motion catalog, and they are exactly what the Jeans equations lead us to expect: $V_{\phi}$ decreases as the asymmetric drift increases and the asymmetric drift increases with the velocity dispersion (see e.g. Binney & Tremaine, 1998, eq. (4.34)). In Figure 11 the velocity dispersion trends in $V_{\phi}$ are shown for the RAVE results and the Galaxia model using the RC-assumption distances. The increase of velocity dispersion with increasing $Z$ and decreasing $R$ is seen in both the data and the model. Consequently, whether we move inwards or away from the plane to regions of increased velocity dispersion, $V_{\phi}$ will decrease.

Figure 8 shows that adopting a different proper-motion catalog does slightly change the predicted values of $V_{\phi}$: away from the plane slightly higher values of order $\Delta\simeq 5-10\,\ \mathrm{kms}^{-1}$ are obtained with the SPM4 proper motions than those in the RAVE catalog, while the UCAC3 proper motions give intermediate values ($\Delta\simeq 2-5\,\ \mathrm{kms}^{-1}$). However, these differences between results from various proper-motion catalogs are smaller than the difference between the observations and the predictions of Galaxia: the latter predicts a markedly flatter profile of $V_{\phi}$ with $Z$. In other words, the data imply that $V_{\phi}$ falls away as we move away from the plane significantly more rapidly than Galaxia predicts. The data give values of $V_{\phi}$ $5\,\ \mathrm{kms}^{-1}$ greater than Galaxia at $Z=0\,\mathrm{kpc}$, falling below the model at $|Z|=0.5\,\mathrm{kpc}$ to be $-10\,\ \mathrm{kms}^{-1}$ lower at $|Z|=1\mathrm{kpc}$.

In Section 4.1.1 we saw that the assumption of a single RC magnitude can lead to an under-estimation of $V_{\phi}$ by $\sim 5\,\ \mathrm{kms}^{-1}$. This effect is also evident in Figure 8 in that the dotted green lines obtained for the Galaxia model using RC distances lie below the dashed blue lines obtained with the true distances, again by about $5\,\ \mathrm{kms}^{-1}$. Further, in Section 4.1.2 and Figure 5 we saw that for large values of $Z$, $M_{K}$ normalizations B and C give lower values of $V_{\phi}$ as large as $15\,\ \mathrm{kms}^{-1}$. The effect was only at the extremities of the data however. Since even the green lines lie above the data for $|Z|>0.5\,\mathrm{kpc}$, we conclude that the use of a single RC magnitude or the chosen $M_{K}$ normalization does not explain the offset between the data and the predictions of Galaxia and the Galaxy’s velocity field must differ materially from that input into Galaxia. Comparison of the top and bottom panels of Figures 10 reveals that, compared to the predictions of Galaxia, the observations show more clearly the expected tendency for $V_{\phi}$ to decrease as we move in at fixed $|Z|$.

In their study of the HTDC, Parker et al. 2004 measured the line-of-sight velocities of thick-disk and halo stars, finding that in Quadrant 1 the rotation of this body of stars lagged the LSR by 80-90$\,\ \mathrm{kms}^{-1}$, while in Quadrant 4 the corresponding lag was only $20\ \mathrm{kms}^{-1}$. Far from confirming this effect, the top-left panel of Figure 14 suggests if anything the opposite is true: at $X\simeq 1.5\mathrm{kpc}$ and $-1<Z<-0.5\,\mathrm{kpc}$, $V_{\phi}$ is lower in the Quadrant 4 than Quadrant 1.

The red curves of Figure 14 show the positions of spiral arms. Close to the plane, at $-0.5<Z<0.5$ there is some hint that there is an increasing lag in $V_{\phi}$ associated with the spiral features. Further from the plane, such an association is less clear, as one might expect.

The $(R,\,Z)$ dependence of $V_{\phi}$ is best described by the power-law:

with Table 4 giving the coefficients for the three different proper motion sources.

The form is not dissimilar to that given for high-metallicity stars in Ivezić et al. 2008 between $2<R<16\,\mathrm{kpc}$ and $0<|Z|<9\,\mathrm{kpc}$ ($\langle V_{Y}\rangle=20.1+19.2|Z/\mathrm{kpc}|^{1.25}\,\ \mathrm{kms}^{-1}$) and Bond et al. 2010 between $0<R<20\,\mathrm{kpc}$ and $0<|Z|<10\,\mathrm{kpc}$ ($\langle V_{\phi}\rangle=-205+19.2|Z/kpc|^{1.25}\,\ \mathrm{kms}^{-1}$). These results are not directly comparable to our fitting formula as we have an extra component that takes radial $R$ dependence into account, which may lead to interaction between the fit coefficients. Nevertheless, our range of values for the exponent (1.06-1.38) for the vertical $Z$ dependence is similar to their results (1.25).

In the top panels of Figure 12, contours of constant $V_{R}$ are by no means vertical, whether one adopts the proper motions in the RAVE or SPM4 catalogs. In fact, with the SPM4 proper motions the contours are not far from horizontal. Thus the radial gradient $\delta V_{R}/\delta R=-3\ \mathrm{kms}^{-1}/\mathrm{kpc}$ reported by S11 is at the very least just one aspect of a complex phenomenon. In fact, if the SPM4 proper motions are correct, the trend in $V_{R}$ is essentially that the further stars are from the plane, the more they are moving away from the Galactic centre.

The middle panel of Figure 8 shows this situation in a different way by showing that the red curves for the SPM4 data are horizontal within the errors at all distances from the plane, but they move down from $V_{R}\simeq 8\ \mathrm{kms}^{-1}$ at $Z\simeq-1.25\,\mathrm{kpc}$ to zero in the plane and back up to $\sim 8\ \mathrm{kms}^{-1}$ at $Z\sim 0.75\,\mathrm{kpc}$. In this figure the black curves that join the RAVE data-points tell a different story in that in the panels for $Z\la 0$ they slope firmly downwards to the right. In the two panels for $Z>0$ the data points from the RAVE proper motions are consistent with no trend in $V_{R}$ with $R$ with the exception of the innermost point in the panel for $Z\sim 0.25\,\mathrm{kpc}$, which carries a large errorbar. The present analysis is consistent with the value of $\delta V_{R}/\delta R$ reported by S11 in that $-3\ \mathrm{kms}^{-1}/\mathrm{kpc}$ is roughly the average of the gradient $\delta V_{R}/\delta R\simeq-7\>\mathrm{to}\>-8\ \mathrm{kms}^{-1}/\mathrm{kpc}$ given by the RAVE proper motions at $Z<0$ and the vanishing gradient at $Z>0$.

If the dominant gradient in $V_{R}$ is essentially in the vertical direction and an even function of $Z$ as the SPM4 proper motions imply, the suspicion arises that it is an artifact generated by the clear, and expected, gradient of the same type that we see in $V_{\phi}$. The gradient could be then seen to be caused by systematics in the proper motions creating a correlation between the measured value of $V_{R}$ and $V_{\phi}$ (see e.g. Schönrich, Binney & Asplund 2012). In Section 7 we re-explore the line-of-sight detection of the $V_{R}$, which is a proper-motion-free approach to observing the radial gradient, which however corroborates the existence of a gradient and North-South differences.

Schönrich 2012 found a larger value of $U_{\odot}=14\,\ \mathrm{kms}^{-1}$ compared to that used in this study and S11, and suggest that this larger value would reduce the gradient in $V_{R}$. In Figure 9 we plot $V_{R}^{s}$, the Galactocentric radial velocity with $U_{\odot}=14\,\ \mathrm{kms}^{-1}$. The qualitative trends are unaffected by the higher value of $U_{\odot}$, with the values only shifted down, bringing the values in the $0.5<Z<1.0\,\mathrm{kpc}$ bin closer to those predicted by Galaxia.

In Figure 14 the location of the spiral arms is overlaid for reference on the full $XYZ$ plots of $V_{R}$ trends. There is some indication that structure in $V_{R}$ is coincident with spiral arms in the $-1.0<Z<-0.5\,\mathrm{kpc}$ region, with Quadrant 3 and Quadrant 4 showing the largest gradient. With a scale height of the thin disk of 300 pc (Gilmore &Reid, 1983), at these large $Z$ values we would not expect however that spiral resonances would play a role. However, the recent results of S12 which explained the gradient in terms of such resonances suggests that they do exert an influence at large distances from the Galactic plane. The diminishment of the gradient at positive $Z$ values is also evident in this figure. Further modelling in 3D is required to understand if the North-South differences in the gradient can be explained in terms of resonances due to spiral structure, which one would expect to produce effects symmetric with respect to $Z$.

The bottom panels of Figure 12 shows the corresponding Galaxia model for the red clump. Here we see that there is an absence of the observed North-South asymmetry, and further, we confirm that the asymmetry is not produced by using the red clump distance estimate. We have explored the influence of changing the distances using the Zwitter distances. We found that these alternative distance calculations also show the same overall $R,\ Z$ variations for $V_{R}$, while exhibiting some differences in various regions. This indicates that while the distances may influence the quantitative results, qualitatively the overall trends remain the same.

The top panels of Figure 13 show that both proper-motion catalogs yield a similar dependence of $V_{Z}$ on $(R,Z)$. We see a ridge of enhanced $V_{Z}$ that slopes across the panel, making an angle of approximately $40^{\circ}$ with the $R$ axis and cutting that axis at approximately the solar radius. This feature is most pronounced in the bottom-left panel, for proper motions in the RAVE catalog. If we assume the LSR $V_{Z}$ value is zero, and thus the top of the ridge has $V_{Z}=0$, we could interpret the stars interior to this point showing an overall rarefaction: stars below the plane move downwards, while those above the plane, move upwards. Exterior to $R=8\,\mathrm{kpc}$ the behaviour is reversed and shows a compression; below the plane stars move upwards, while above the plane stars move downwards. The amplitude of these variations is large; up to $|V_{Z}|=17\,\ \mathrm{kms}^{-1}$ as seen from Figure 8.

The results for the SPM4 and UCAC3 proper motions do not exhibit quite so large an amplitude in the $V_{Z}$ variations, but as shown in Figure 13 for SPM4, exhibit a similar overall behaviour, with the line of higher $V_{Z}$ running diagonal across the $RZ$ plane. However, in Figure 8 we see that below the plane, the alternative proper motion results do not show the large dip in $V_{Z}$ at $R=9\,\mathrm{kpc}$, though this is mostly due to the smaller sampling area in the $XY$ plane of the alternative proper motion sources. Above the plane, the behaviour of the results from the various proper motion sources is similar, with clear deviation from the Galaxia model’s predictions in the bottom panels of Figure 13. The Galaxia model itself shows a random pattern of high and low $V_{Z}$ of the order of $5\,\ \mathrm{kms}^{-1}$. There is some smoothing spatially of these fluctuations with the use of RC-magnitude derived distances, though the observed pattern is not generated with this distance method. In general, the proper motions dominate the calculated $V_{Z}$ values, and the differences in observed structure reflect this. Nonetheless, it is reassuring that the overall pattern is recovered in Figure 13 with each proper motion catalog source.

Recently, Widrow et al. 2012 found with main-sequence SDSS stars that the vertical number-density and $V_{Z}$ profiles suggest vertical waves in the Galactic disk excited by a recent perturbation. A direct comparison with their results and ours is not possible as their stars are outside the RAVE spatial sampling region. However, our results support this proposition: the rarefaction and compression behaviour seen in the top panels of Figure 13 is indicative of wave-like behaviour. This could then be seen as further “ringing” behaviour of the disk caused by a recent accretion event, as in Minchev et al. 2009; Gomez et al. 2012. Indeed, Gomez et al. 2013 further suggest that such vertical waves may have been excited by the Sagittarius dwarf as it passed through the disk. Their simulations find deviations of up to $8\,\ \mathrm{kms}^{-1}$ in $V_{Z}$ in a heavy-Sgr scenario.

Turning now to Figure 14 we find some suggestion that the $V_{Z}$ field shows some alignment with spiral-arms feature, even at significant distances from the plane. Below the plane, the Sagittarius-Carina arm at $X\sim 1.5\,\mathrm{kpc}$ is aligned with the low $V_{Z}$ velocities seen with the RAVE catalog proper motions for $Y<0$. In this representation, the SPM4 and UCAC3 proper motion results also show a dip in their $V_{Z}$ velocities around this area, though the feature is not as pronounced as for the RAVE catalog proper motions. The Perseus arm is also aligned with an area of lower $V_{Z}$ at $X=-1.25,\ Y<0,\ Z<-0.5\,\mathrm{kpc}$. Across all values of $Z$, the alignment of the $V_{Z}$ features with the spiral arms is less clear for $Y>0$, which below the plane may be explained by the proximity of the Hercules thick disk cloud.Siebert et al. 2012 showed how spiral density waves could be used to explain the $V_{R}$ streaming motion. Siebert et al. (in preparation) goes further to investigate how they could also create the $V_{Z}$ structures, coupled to $V_{R}$.

These results suggest that there are large-scale vertical movements of stars, with various cohorts at the same $Z$ moving in opposing directions. In Casetti et al. 2011, it was proposed that while the Galactic warp starts at the solar radius, the small elevations of the warp at these distances ($70-200\,\textrm{pc}$ for $R=8-10\,\mathrm{kpc}$ (López-Corredoira et al., 2002)) and the large distances of the RAVE stars from the plane would mean that kinematics of stars would be little affected by the warp. In Russiel 2002, the nearby Sagittarius-Carina arm as traced by star-forming complexes is found to lie mostly below the plane by $100-200\textrm{pc}$, which suggests that nearby spiral arms can be influenced by some warping. Whether the warp is long-lived or a transient feature would have implications for the associated velocities. The complexity of the vertical velocity distribution suggested by our results would tend to point towards transient features: they suggest a non-equilibrium state. A multi-mode travelling wave caused by a recent accretion event, or structure associated with the disk’s spiral arms, would be the most likely explanations for the observed velocity structure.

The second moments of the velocity components increase with increasing $Z$ and decreasing $R$ in a smooth fashion as expected, with Figure 11 showing the results for $\sigma(V_{\phi})$. As for $V_{\phi}$, they can be fit by a simple parametric function.

Figure 15 displays the trends in the dispersion of $V_{\phi},\ V_{R},\ V_{Z}$ as a function of $Z$ for $0.5\,\mathrm{kpc}$ thick slices in $R$ using $0.5\,\mathrm{kpc}$ bins in $Z$. Also plotted are the results based on UCAC3 and SPM4 proper motions, and the Galaxia model results using the “true” model distances plus RC-calculated distances. Unlike the mean-value trends, the differences between the various proper motion sources is minimal. We see however that the Galaxia model results indicate that the use of RC-calculated distances can lead to an increase in the measured dispersion values, particularly for $V_{R}$. Interestingly however, there is a reasonable agreement between the results of Galaxia and that observed with RAVE.

We note that the dispersions were corrected for errors following the same methodology of Casetti et al. 2011. This method adjusts a guess of the true dispersion, $\sigma_{V}$, in velocity component $V$ until

where $\sigma_{T}^{2}=\sigma_{V}^{2}+eV_{i}^{2}$, with $eV_{i}$ the MC-derived error in $V$ for $i$th star and $\bar{V}=\mathrm{mean}(V_{i})$ is the mean over the sample.

The velocity dispersions exhibit linear behaviour in $R$ and parabolic behaviour in $Z$. Thus, a fit to the dispersion trends was found to be provided by the simple form

with Table 5 giving the coefficients for the three velocity components for the three different proper motion sources. The similar coefficients for each velocity component underlines the similarity of the trends between the three sources.

Figure 3 from S11 gives the projection onto the plane of the mean $V_{\mathrm{los}}$ in bins $200\,\textrm{pc}$ wide as a function of $d\cos{l}\cos{b}$, with the latter a proxy for $X$ in their Galactic co-ordinates, with $X_{\textrm{S11}}\sim R_{\odot}-d\cos{l}\cos{b}$. For $|l|<5\deg$ and $175\deg<l<185\deg$ the values of the radial component, $V_{R}$ is essentially the same as the Cartesian value $-U$. As most stars in these narrow cones have $\cos{b}/\sin{b}\sim 1.5$, $U$ (and so $-V_{R}$) is dominated by the term $V_{\mathrm{los}}\cos{b}$. In the Appendix we list the relevant equations which shows how this follows. S11b compared the mean velocities to those expected for a thin disk in circular rotation and adding a radial gradient of $\partial(V_{R})/\partial R=-3,-5$ and $-10\ \ \mathrm{kms}^{-1}/\mathrm{kpc}$ (plus a thick disk lagging the LSR by $36\,\ \mathrm{kms}^{-1}$), finding that the results were consistent with a radial gradient of $\partial(V_{R})/\partial R=-3\ \mathrm{to}-5\,\ \mathrm{kms}^{-1}/\mathrm{kpc}$.

From Equations 1 and 7 in the Appendix, we see that if we are to consider the $V_{\mathrm{los}}$ as representative of $U$ (and so $-V_{R}$), technically it would be better to first subtract the solar motion from the line-of-sight velocity, to obtain $V_{\mathrm{los}}^{\prime}$. Otherwise the components above and below the plane will be shifted from each other, which can be problematic if the sample shifts from above-the-plane to below-the-plane, potentially producing spurious trends. This was omitted in the analysis of Figure 3 from S11, however, it did not affect the conclusions as the models were similarly shifted.

In Figure 16 we re-examine the trends in $V_{\mathrm{los}}\cos{b}$ for the RC giants and stars with Zwitter distances, where we plot against $R$ rather than $d\cos{l}cos{b}$ as this is easier to interpret and there is little difference between the two values. We plot the results with and without the LSR correction, with the curve being flattened in the latter with the removal of the $U_{\odot}$ component. We also plot the trends for a thin-disk as in S11, using 100 Monte Carlo realizations. We see here that not correcting for the solar motion does not affect the conclusions. Note that with the updated data sets the $V_{R}$ streaming motion is slightly less apparent in both the RC and Zwitter results than it was in S11: the model with $\partial(V_{R})/\partial R=-3\,\ \mathrm{kms}^{-1}/\mathrm{kpc}$ is a better fit to the data. Also, the differences between the RC and Zwitter trends around the solar radius in Figure 16 is a result of the different geometric sampling of the two.

In Figure 17 we examine the $V_{\mathrm{los}}^{\prime}\cos{b}$ results split into North and South components. In this we open up the criteria a bit to within $7\deg$ of the centre and anti-centre direction to reduce the effects of Poisson noise and make the trends clearer. Here we see that that differences can be discerned between the North and South trends, however, they are opposite to that seen in the previous section; the gradient is in the North ($-5\ \mathrm{to}-10\,\ \mathrm{kms}^{-1}/\mathrm{kpc}$) is stronger though absent in the South ($0\,\ \mathrm{kms}^{-1}/\mathrm{kpc}$)) for $R<8\,\mathrm{kpc}$. Plotting the corresponding results for $V_{R}$ in Figure 18, we find that the including the proper motion results shifts the overall values so that the northern trends are more in line with $\delta V_{Z}/\delta R=-10\,\ \mathrm{kms}^{-1}/\mathrm{kpc}$ and the southern $\delta V_{Z}/\delta R=-3\ \mathrm{to}-5\,\ \mathrm{kms}^{-1}/\mathrm{kpc}$. Despite the disparity with the actual numbers, the cause of which is discussed below, both $V_{\mathrm{los}}^{\prime}\cos{b}$ and the $V_{R}$ nonetheless exhibit differences between the north and south trends, opposite to those found in Section 6.2.

To understand the reversal of the trends note first that the selected sample cuts across a range of $Z$ as we change $R$ so these plots combine the $R$ and $Z$ trends seen previously. Second, the necessary restriction on $l$ means that we are sampling a very narrow beam and thus do not see the global patterns, but only those along that beam. In Figure 14 Quadrants 3 and 4 show the largest gradient, which we miss with the beams. Hence, these plots emphasize the 3D nature of the $V_{R}$ values in the solar neighbourhood: what you measure very much depends on where you look, be it north or south, and at different $R$ and $Z$ values.

From Equation 6 in the Appendix we can see that $V^{\prime}_{b}\cos{b}$ can be used as a proxy for $W=V_{Z}$ if $|\cos{b}|>|\sin{b}|$. This condition is for the most part met in the low latitudes sampled along the $|l|<7\deg$ and $173\deg<l<187\deg$ cones used above. So in Figure 19 we examine the trends above and below the plane for $V^{\prime}_{b}\cos{b}$ in these directions, with Figure 20 giving the corresponding trends in the total $V_{Z}$. Note that we use RC and Zwitter distances here to provide the abscissa, supplementing the proper motion data. We also plot both positive and negative trends in $V_{Z}$ for a thin disk model as above, with values from $\partial(V_{Z})/\partial R=-10,-5,-3,0,3,5,10\ \ \mathrm{kms}^{-1}/\mathrm{kpc}$.

Both Figures 19 and 20 show that the trends above and below the plane are markedly different, with the $V_{Z}$ decreasing with R above the plane at a rate of $\partial(V_{Z})/\partial R\sim-10\ \ \mathrm{kms}^{-1}/\mathrm{kpc}$ according to the $V^{\prime}_{b}\cos{b}$ plot and $-5\ \ \mathrm{kms}^{-1}/\mathrm{kpc}$, to the $V_{Z}$ plot. Below the plane, there is a positive trend of $\partial(V_{Z})/\partial R\sim+5\,\ \mathrm{kms}^{-1}/\mathrm{kpc}$ in both plots. This positive trend in $V_{Z}$ appears to stop just beyond the solar circle at $R\sim 8.5\,\mathrm{kpc}$. In contrast to $V_{R}$ above, this behaviour is consistent with what was found in Section 6.3.

Note that the differences between the observed magnitude of the trends in $V_{\mathrm{los}}^{\prime}\cos{b}$ and $V_{R}$, plus $V^{\prime}_{b}\cos{b}$ and $V_{Z}$ can be explained by the fact that both $V_{\mathrm{los}}^{\prime}$ and $V^{\prime}_{b}$ contain components of $U$ and $W$. For small $l$, $V_{\mathrm{los}}^{\prime}\cos{b}=U^{\prime}\cos^{2}{b}+W^{\prime}\cos{b}\sin{b}$ and $V^{\prime}_{b}\cos{b}=-U^{\prime}\sin{b}\cos{b}+W^{\prime}\cos^{2}{b}$. The cross-terms mean that there is some ‘leakage’ of the $V_{Z}$ trends into $V_{\mathrm{los}}^{\prime}\cos{b}$ and $V_{R}$ trends into $V^{\prime}_{b}\cos{b}$, which work to either diminish or enhance the observed trend in the single component. Nonetheless, the fact that there are differences between the north and south for the single components at all is further evidence of the 3D variations of the velocity values.

The $V_{R}$ and $V_{Z}$ trends, as seen via the line-of-sight velocities and proper motions respectively, are unaffected by potential systematics in the distances: a change of the distance scale would not affect the fact that a trend is seen at all. This is particularly so for the differences observed between the northern and southern samples. Furthermore, both the SPM and UCAC3 proper motions give similar results as in Figures 19, with a stark negative trend above the plane and a positive trend below for the region $R<8.5\,\mathrm{kpc}$. Thus, neither the distances nor the proper motions introduce large systematics into this method of detection.