The tilt of the velocity ellipsoid in the Milky Way with Gaia DR2

TL;DR

This study analyzes Gaia DR2 data to understand the orientation of stellar velocity distributions in the Milky Way, revealing how tilt angles vary with galactic radius and highlighting the impact of parallax measurement errors.

Contribution

It introduces a detailed analysis of the velocity ellipsoid tilt angles across the Galaxy and models their radial dependence, considering Gaia DR2 systematic errors.

Findings

Tilt angles change from spherical to cylindrical alignment with increasing radius.

Systematic parallax errors significantly influence tilt angle measurements.

Tilt angles are consistent across different stellar populations and galactic azimuths.

Abstract

The velocity distribution of stars is a sensitive probe of the gravitational potential of the Galaxy, and hence of its dark matter distribution. In particular, the shape of the dark halo (e.g. spherical, oblate, or prolate) determines velocity correlations, and different halo geometries are expected to result in measurable differences. Here we explore and interpret the correlations in the $(v_R, v_z)$-velocity distribution as a function of position in the Milky Way. We selected a high-quality sample of stars from the Gaia DR2 catalogue and characterised the orientation of the velocity distribution or tilt angle over a radial distance range of $[4-13]~$kpc and up to $3.5~$kpc away from the Galactic plane while taking into account the effects of the measurement errors. We find that the tilt angles change from spherical alignment in the inner Galaxy ($R\sim4~$kpc) towards more cylindrical alignments in the outer Galaxy ($R\sim11~$kpc) when using distances that take a global zero-point offset in the parallax of $-29~\mu$as. However, if the amplitude of this offset is underestimated, then the inferred tilt angles in the outer Galaxy only appear shallower and are intrinsically more consistent with spherical alignment for an offset as large as $-54~\mu$as. We further find that the tilt angles do not seem to strongly vary with Galactic azimuth and that different stellar populations depict similar tilt angles. Therefore we introduce a simple analytic function that describes the trends found over the full radial range. Since the systematic parallax errors in Gaia DR2 depend on celestial position, magnitude, and colour in complex ways, it is not possible to fully correct for them. Therefore it will be particularly important for dynamical modelling of the Milky Way to thoroughly characterise the systematics in astrometry in future Gaia data releases.

Equations

$$\tan (2\gamma )=\frac{2\mathrm{cov}(v_R,v_z)}{\mathrm{var}(v_R)-\mathrm{var}(v_z)},$$

$$\tan (2\alpha )\equiv \frac{2\mathrm{cov}(v_r,v_{\theta })}{\mathrm{var}(v_r)-\mathrm{var}(v_{\theta })}.$$

$$\tan (2\gamma )=-\tan (2\theta +2\alpha ).$$

$$\text{@vec}{\Sigma }_{\mathrm{error},\mathrm{i}}=\left[\begin{array}{cc}\mathrm{var}(v_{R,i})& \mathrm{cov}(v_{R,i},v_{z,i})\\ \mathrm{cov}(v_{R,i},v_{z,i})& \mathrm{var}(v_{z,i})\end{array}\right],$$

$$\begin{array}{rl}{L}_i& =L_i[⟨v_R⟩,\sigma (v_R{)}_{\text{intr}},⟨v_z⟩,\sigma (v_z{)}_{\text{intr}},\text{cov}(v_R,v_z{)}_{\text{intr}}]\\ & =\frac{1}{\sqrt{\text{det}(2\pi \text{@vec}{\Sigma }_i)}}\exp \left[-\frac{1}{2}(\text{@vec}{x}_i-\text{@vec}\mu {)}^{⊺}\text{@vec}{{\Sigma }_i}^{-1}(\text{@vec}{x}_i-\text{@vec}\mu )\right],\end{array}$$

$$\hat{\varpi }=\varpi +0.029\mathrm{mas},\mathrm{and}{\widehat{ϵ}}_{\varpi }=\sqrt{{ϵ}_{\varpi }^2+0.043^2}.$$

$$\begin{array}{rl}\alpha (R,z)=\alpha (R_0,z_0)& +a_1(R-R_0)+a_2(z-z_0)\\ & +a_3(R-R_0)(z-z_0)\\ & +a_4(R-R_0{)}^2+a_5(z-z_0{)}^2+...,\end{array}$$

$$\alpha (R,z)\approx a_3(R-R_0)z.$$

$${\chi }^2=\sum _{j=1}^{N_{\mathrm{bins}}}{\left(\frac{\alpha (R_j,z_j{)}_{\mathrm{model}}-\alpha (R_j,z_j{)}_{\mathrm{obs}}}{ϵ[\alpha (R_j,z_j{)}_{\mathrm{obs}}]}\right)}^2.$$

$$\mathrm{var}(s^2)=\frac{1}{N}\left({\mu }_4-\frac{N-3}{N-1}\mathrm{var}(v{)}^2\right)$$

$$\mathrm{var}(S_{xy})=\frac{1}{N}\left[{\mu }_{22}-\frac{N-2}{N-1}\mathrm{cov}(x,y{)}^2+\frac{1}{N-1}\mathrm{var}(x)\mathrm{var}(y)\right],$$

$$\left(\begin{array}{c}U\\ V\\ W\end{array}\right)=\left(\begin{array}{ccc}\cos (l)\cos (b)& -\sin (l)& -\cos (l)\sin (b)\\ \sin (l)\cos (b)& \cos (l)& -\sin (l)\sin (b)\\ \sin (b)& 0& \cos (b)\end{array}\right)\left(\begin{array}{c}{v}_{\mathrm{los}}\\ \frac{k}{\varpi }{\mu }_{l^{\star }}\\ \frac{k}{\varpi }{\mu }_b\end{array}\right).$$

$${\left(\begin{array}{c}U\\ V\\ W\end{array}\right)}_1={\left(\begin{array}{c}U\\ V\\ W\end{array}\right)}_0+{\frac{\partial }{\partial \varpi }\left(\begin{array}{c}U\\ V\\ W\end{array}\right)|}_0\Delta \varpi +O(\Delta {\varpi }^2).$$

$$\frac{\partial }{\partial \varpi }\left(\begin{array}{c}U\\ V\\ W\end{array}\right)=-\frac{1}{\varpi }\left(\begin{array}{c}U-\cos (l)\cos (b)v_{\mathrm{los}}\\ V-\sin (l)\cos (b)v_{\mathrm{los}}\\ W-\sin (b)v_{\mathrm{los}}\end{array}\right),$$

$$v_{\mathrm{los}}=\cos (b)\cos (l)U+\cos (b)\sin (l)V+\sin (b)W.$$

$$\tan (2\delta )=\frac{2\mathrm{cov}(U,W)}{\mathrm{var}(U)-\mathrm{var}(W)}.$$

$$\tan (2{\delta }_1)\simeq \frac{2⟨U_0{W}_0⟩+{ϵ}_{\mathrm{A}}}{⟨U_0^2⟩-⟨W_0^2⟩+{ϵ}_{\mathrm{B}}},$$

$${ϵ}_A$$

$${ϵ}_B$$

$${\delta }_1$$

$$\Delta \delta$$

$$\tan (2{\delta }_1)$$

$${ϵ}_{\mathrm{C}}$$

$${ϵ}_{\mathrm{D}}$$

$$x\simeq {ϵ}_{\mathrm{D}}-{ϵ}_{\mathrm{C}}.$$

$$b=0^{\circ }:$$

$$|b|=90^{\circ }:$$

$$b=+45^{\circ }:$$

$$b=-45^{\circ }:$$

$${\delta }_1\simeq \pm \frac{1}{2}\sin (2b)\left(\frac{⟨U_0^2⟩+⟨W_0^2⟩}{⟨U_0^2⟩-⟨W_0^2⟩}\right)\left(\frac{\Delta \varpi }{{\varpi }_0}\right),$$

Full text

1 Introduction

The second data release of the Gaia space mission (Gaia Collaboration et al. 2018a) contains more than $1.3$ billion stars with measured proper motions and positions and a subset of over $7$ million stars with full six-dimensional (6D) phase-space information. The availability of the motions and positions of stars in the Milky Way and its satellite galaxies has already led to new insights about the Galaxy (e.g. Antoja et al. 2018; Belokurov et al. 2018; Helmi et al. 2018; Poggio et al. 2018; Price-Whelan & Bonaca 2018), and many more discoveries will likely follow before Gaia’s next data release.

Studies of the Galaxy provide insight about the formation and evolution of galaxies in general, and hence about elements of the cosmological paradigm. For example, detailed dynamical modelling of the Milky Way and its satellites, and in particular their mass distribution, provide critical constraints on the nature of dark matter (e.g. Bonaca et al. 2018). Mass models of the Galaxy, such as those by McMillan (2011), Piffl et al. (2014), and McMillan (2017), have been developed to fit many different observational constraints simultaneously, although this is very challenging. Therefore many works often focus on a specific aspect such as the characterisation of the velocity distribution across the Galaxy.

The in-plane velocity distribution $f(v_{R},v_{\phi})$ in the Solar vicinity has long been known to be complex, and many moving groups are known to exist (e.g. Proctor 1869; Eggen 1965; Dehnen 1998; Antoja et al. 2008). With Gaia DR2 the level of detail visible in the velocity distribution of stars has increased immensely (see e.g. Gaia Collaboration et al. 2018c; Antoja et al. 2018), and a plethora of substructures have become apparent. On the other hand, the 2D velocity distribution describing the radial and vertical velocity components, $f(v_{R},v_{z})$, shows significantly less substructure and the traditional velocity moments can still describe the data well to first order.

Such velocity moments and thus the axial ratios of the velocity ellipsoid, however, depend on the stellar distribution function and are different for different populations of stars. In contrast, its orientation (or better known as alignment or tilt) is directly related to (the shape of) the underlying gravitational potential in which the stars move (e.g. van de Ven et al. 2003; Binney & Tremaine 2008; Binney & McMillan 2011; An & Evans 2016) and is the focus of this paper.

Nearly spherically aligned velocity ellipsoids were found for the halo (Smith et al. 2009; Bond et al. 2010; King et al. 2015; Evans et al. 2016) by mainly using data from the Sloan Digital Sky Survey (York et al. 2000). Similar findings were obtained by Posti et al. (2018) for dynamically selected nearby halo stars. These authors obtained full 6D phase-space information by combining radial velocity measurements from the RAdial Velocity Experiment (RAVE DR5, Kunder et al. 2017) to the 5D subset of the Gaia DR1 catalogue (Gaia Collaboration et al. 2016). Most recently, Wegg et al. (2018) used $15,651$ RR Lyrae halo stars with accurate proper motions from Gaia DR2 and also inferred a nearly spherically aligned velocity ellipsoid over a large range of distances between $1.5$ kpc and $20$ kpc from the Galactic centre. When fed into the Jeans equations, this result seems to imply a spherical dark matter distribution.

Studies focusing on the orientation of the velocity ellipsoid in local samples of the Milky Way disk have also been consistently reporting (close to) spherical alignment. Siebert et al. (2008) have used RAVE DR2 and found a tilt angle $\gamma$ equal to $7.3^{\circ}\pm 1.8^{\circ}$ for red clump stars at $R=R_{\odot}$ and $z=1$ kpc, where $\gamma_{\textrm{sph}}=7.1^{\circ}$ would be expected for spherical alignment at this location. Casetti-Dinescu et al. (2011) found $8.6^{\circ}\pm 1.8^{\circ}$ for a sample of stars with heights between $0.7$ kpc and $2.0$ kpc and representative of the metal-rich thick disk, which can be compared to $\gamma_{\textrm{sph}}=8.0^{\circ}$ given the mean location of the sample. Subsequently, Smith et al. (2012) reinforced these findings using data from the Sloan Digital Sky Survey DR7 (SDSS, Abazajian et al. 2009). Binney et al. (2014) using RAVE data, and Büdenbender et al. (2015), using Sloan Extension for Galactic Understanding and Exploration (SEGUE, Yanny et al. 2009), characterised the tilt angle around the Galactic radius of the Sun up to $z\sim 2.0$ kpc by $\gamma(z)\approx a_{0}\arctan(z/R_{\odot})$. They found $a_{0}\sim 0.8$ and $a_{0}=0.9\pm 0.04$ respectively, values close to, but significantly different from, spherical alignment for which $a_{0}$ = 1.0. Recently, Mackereth et al. (2019) have analysed the kinematics of mono-age, mono-[Fe/H] populations for both low and high [$\alpha$/Fe] samples. They have cross matched the Apache Point Observatory Galactic Evolution Experiment (APOGEE DR14, Majewski et al. 2017) with Gaia DR2 to obtain a sample of $65,719$ red giant stars located between $4$ kpc and $13$ kpc in Galactic radius and up to $2$ kpc from the Galactic plane. Mackereth et al. (2019) report that the tilt angles found are consistent with spherical alignment for all populations, although they note that the uncertainties are very large.

In this work we characterise the orientation of the velocity ellipsoid over a larger section of the Milky Way by using a dataset of more than 5 million stars from Gaia DR2. The paper is organised as follows. In Sect. 2 the dataset is introduced as well as the selection criteria applied. In Sect. 3 we characterise the velocity distribution and the measurement errors. The results are presented in Sect. 4. In that section we also explore differences with azimuth, investigate trends with stellar populations, and put forward a fit that reproduces the variation of the tilt angle with position in the Galaxy. In Sect. 5 we explore the effect of systematic errors on our measurements and show that the systematic parallax errors present in Gaia DR2 have a significant impact on the tilt angles found. In that Section we therefore also discuss our findings in the context of Galactic models. We summarise in Sect. 6.

2 Data

We used the subset of Gaia DR2 with full 6D information (Gaia Collaboration et al. 2018c). We use the Bayesian distance estimates $\hat{d}$ provided by McMillan (2018) who uses the Gaia DR2 parallaxes $\varpi$ and $G_{\textrm{RVS}}$ magnitudes as input. McMillan (2018) takes into account Gaia DR2’s overall parallax offset of $-29\mu$as with a RMS error of $43\mu$as (Lindegren et al. 2018).

To construct a high-quality sample we select stars with at most 20% relative distance errors, that is $\hat{d}/\hat{\epsilon}(\hat{d})>5$, and $\hat{d}<5$ kpc. The sample contains $5,796,226$ stars. Stars with $\hat{d}<1$ kpc, typically have distances better than 5% (median 2.8%) and for stars at $4<\hat{d}<5$ kpc the relative distance errors are in between 12% and 20% (median 17.1%).

In Fig. 1 we show the extent of our sample in a number density map. To compute the Galactocentric cylindrical coordinates $(R,z,\phi)$, we assume111Use of the value of $R_{\odot}=8178\pm 13_{stat.}\pm 22_{sys.}$ pc, as determined by Gravity Collaboration et al. (2019), does not affect the main conclusions of this paper. $R_{\odot}=8.3$ kpc (Schönrich 2012) and $z_{\odot}=0.014$ kpc (Binney et al. 1997, and $\phi_{\odot}=180^{\circ}$) for the position of the Sun. Because of the imposed maximum distances to the stars, the sample extends from $R\sim 4$ kpc up to $R\sim 13$ kpc and reaches up to $z=\pm 4$ kpc. The white contours in Fig.1 indicate the location of bins containing $2,000$ and $100$ stars respectively. This shows that Galactic heights up to $\sim 3.5$ kpc are still covered with a statistically significant number of stars.

We derive the velocities of the stars in our sample in a Galactocentric cylindrical coordinate system ($v_{R}$, $v_{z}$, $v_{\phi}$). For the motion of the Local Standard of Rest (LSR), that is the velocity of a circular orbit at $R=R_{\odot}$, we assume $v_{\mathrm{c}}(R_{\odot})=240$ km/s (Piffl et al. 2014; Reid et al. 2014). The peculiar motion of the Sun with respect to the LSR is taken to be $(U,V,W)_{\odot}=(11.1,12.24,7.25)$ km/s (Schönrich et al. 2010), where $U$ denotes motion radially inwards and $V$ in the direction of Galactic rotation (both in the Galactic plane), and $W$ perpendicular to the Galactic plane and in the direction of the Galactic north pole. We propagate the errors and correlations in the observables to determine the errors on the velocities (and their correlations). Here we assume that the Bayesian distances are not correlated with the remaining astrometric parameters. The velocity errors for the stars in our sample at $\hat{d}<1$ kpc are typically smaller than $2$ km/s with a median value of $\sim 1$ km/s for the $v_{R}$-, $v_{z}$-, and $v_{\phi}$-components. At $4<\hat{d}<5$ kpc the median errors are in the range from $\sim 3$ km/s to $\sim 8$ km/s and generally smaller than $15$ km/s.

The characterisation of the kinematics, in terms of the mean motions and velocity dispersions, of a large part of the Milky Way disk have been presented in Gaia Collaboration et al. (2018c) using the 6D dataset from Gaia DR2. This characterisation has put on firm ground the evidence of the presence of streaming motions in all velocity components (Siebert et al. 2008; Williams et al. 2013; Tian et al. 2017; Carrillo et al. 2018) and revealed a large amount of substructure in the velocity distributions. In this paper we proceed to focus on the correlation between the radial and vertical velocity components across a large fraction of the Milky Way galaxy.

3 Methods

The 3D velocity distribution of stars $f(v_{\phi},v_{R},v_{z})$ at a given point in the Galaxy may be characterised by its various moments. As described in the Introduction, the tilt of the velocity ellipsoid refers to the orientation of the 2D velocity distribution $f(v_{R},v_{z})$, which would be obtained by integrating over $v_{\phi}$. As shown in Smith et al. (2009) and Büdenbender et al. (2015), this is equivalent to taking the moments of the 3D velocity distribution and neglecting the cross terms with $v_{\phi}$. These cross-terms are interesting in their own right, as they reveal also other physical mechanisms at work, such as for example the presence of substructures associated to resonances induced by the rotating Galactic bar (Dehnen 1998), but are not the focus of this work.

3.1 The tilt angle: the orientation of the velocity ellipse

In the Galactocentric cylindrical coordinate system we define the tilt angle $\gamma$, following for instance Smith et al. (2009), as:

[TABLE]

which therefore takes values from $-45$ degrees to $+45$ degrees, and is measured counterclockwise (i.e. from the $v_{R}$-axis towards the positive $v_{z}$-axis). For exact cylindrical alignment $\gamma_{\textrm{cyl}}=0^{\circ}$ and the major and minor axis align with the Galactocentric cylindrical coordinates.

It is also possible to define a tilt angle $\alpha$ with respect to the spherical coordinate system $(r,\theta,\phi)$, where $\tan(\theta)\equiv R/z$, that is:

[TABLE]

The tilt angle $\alpha$ thus measures directly the deviation from spherical alignment, which corresponds to $\alpha=0^{\circ}$. In such a case one of the principal axes of the ellipse points to the Galactic centre. The relation between $\alpha$ and $\gamma$ at every $(R,z)$ is

[TABLE]

From now on, we always refer to the tilt angle $\gamma$, thus as defined in the cylindrical coordinate system, unless stated otherwise. To explore the spatial variation of the tilt angle we measure the intrinsic moments of Eq. 1 after projecting all stars onto the $(R,z)$-plane, thus ignoring in the first stage the Galactic azimuthal angle of the stars (although this is considered in Sect. 4.2). We bin the meridional plane as in Fig. 1 and always require at least $100$ stars per bin.

3.2 Accounting for measurement errors

Measurement errors affect the observed velocity moments and can therefore have a significant effect on the inferred tilt angles (Siebert et al. 2008). To establish their effect we here explore two ‘methods’ to account for the errors and for recovering the (intrinsic) velocity moments.

Method 1. We assume that the stars in a given spatial bin have similar measurement errors. This assumption is reasonable because the measurement errors in a particular bin are usually much smaller than the intrinsic velocity dispersion. If the measurement errors were exactly the same for all stars in a bin, the intrinsic velocity covariance matrix can be recovered by subtracting the error covariance matrix from the observed velocity covariance matrix. This follows from the fact that convolving a Gaussian distribution with Gaussian distributed measurement errors again results in a Gaussian with covariance matrix $\@vec{\Sigma_{\mathrm{obs}}}=\@vec{\Sigma_{\mathrm{intr}}}+\@vec{\Sigma_{\mathrm{error}}}$, where $\@vec{\Sigma_{\mathrm{obs}}}$ and $\@vec{\Sigma_{\mathrm{intr}}}$ are the observed and intrinsic covariance matrix of the velocity distribution respectively. In our approximation $\@vec{\Sigma_{\mathrm{error}}}\approx\mathrm{median}\left(\@vec{\Sigma_{\mathrm{error,i}}}\right)$ for

[TABLE]

in which the diagonal terms denote the variance error of the corresponding velocity component of star $i$. Similarly $\mathrm{cov}(v_{R,i},v_{z,i})$ denotes the error covariance for the $(v_{R},v_{z})$ measurements of star $i$. For the required typical errors we take the relevant median errors of the stars in the bin. The recovered intrinsic velocity moments are then used to characterise the velocity distribution. The errors on these moments are analytically estimated and then propagated into uncertainties on the recovered tilt angles. More details can be found in Appendix A.

Method 2. We perform Markov Chain Monte Carlo (MCMC) modelling (Foreman-Mackey et al. 2013) for bins with a smaller number of stars (with $100<N<2,000$). This aims to solve for the intrinsic velocity dispersions $\sigma(v_{R})_{\textrm{intr}}$ and $\sigma(v_{z})_{\textrm{intr}}$, the mean velocities $\langle v_{R}\rangle$ and $\langle v_{z}\rangle$, and the covariance term $\textrm{cov}(v_{R},v_{z})_{\textrm{intr}}$ in each bin. This is done by maximizing the bivariate Gaussian likelihood function $L=\prod^{N}_{i=1}L_{i}$, where

[TABLE]

in which $\@vec{x_{i}}=[v_{R,i},v_{z,i}]$, $\@vec{\mu}=[\langle v_{R}\rangle,\langle v_{z}\rangle]$ and $\@vec{\Sigma_{i}}=\@vec{\Sigma_{\mathrm{intr}}}+\@vec{\Sigma_{\mathrm{error},i}}$. Whereas in Method 1 $\@vec{\Sigma_{\mathrm{error},i}}$ was assumed to be the same for each star, we here use $\@vec{\Sigma_{\mathrm{error},i}}$ for each star separately. We add priors to the model that only allow for positive velocity dispersions in $v_{R}$ and $v_{z}$ and that restrict the correlation coefficient between $v_{R}$ and $v_{z}$ always to be within [$-1$,$1$]. For a given bin, the samples drawn by the MCMC run translate into a distribution of tilt angles. We take the median as the best estimate of the tilt angle. For its error we take half the difference between the tilt angles corresponding to the $16^{\mathrm{th}}$ and $84^{\mathrm{th}}$ percentile.

In general we find that the effect of the measurement errors on the recovered moments is small. Moreover, for most bins the velocity measurement errors are sufficiently similar and small that we may use the computationally much faster Method 1 instead of the MCMC-based deconvolution. We have also compared the results to the case in which we simply compute the variances of the observed stellar velocities in the bins of interest, and take these at face value, meaning that we do not take into account the measurement errors. The results are again rather similar, see for example, Fig. 12 of Appendix A which shows the distributions of the measurement errors for the bin located at $R=11.5$ kpc and $z=1.5$ kpc. In what follows, we use the results from Method 1 unless stated otherwise.

4 Results

We present our measurement of the tilt angles by showing velocity ellipses in the meridional plane. At each position $(R,z)$, we define a set of axes with $v_{R}$ into the $R$-direction and $v_{z}$ in the $z$-direction. The centre of each velocity ellipse is always placed at its position $(R,z)$. The size of the major and minor axis of each ellipse scale with the intrinsic velocity dispersions along these directions. The $R$- and $z$-axis are both scaled by the same constant $c_{x}$. Similarly, all $v_{R}$- and $v_{z}$-axes are scaled by a constant $c_{v}$, thus both sets of axes have an aspect ratio of $1$. As a consequence, the velocity ellipses drawn will actually point to the Galactic centre when there is spherical alignment. As a reference, the inset in the figures shows a velocity distribution aligned in cylindrical coordinates and with $\sigma(v_{R})=100$ km/s and $\sigma(v_{z})=50$ km/s (unless stated otherwise).

4.1 Tilt angles projected onto the $(R,z)$-plane

Fig. 2 shows the velocity ellipses colour-coded by their angular misalignment with respect to spherical alignment. For $z\geq 0$ kpc we define this misalignment as $\gamma-\gamma_{\textrm{sph}}$, whereas for $z<0$ kpc the misalignment is $\gamma_{\textrm{sph}}-\gamma$. Steeper tilt angles result in positive misalignment (from light to dark red), shallower tilt angles in negative misalignment (from light to dark blue). Ellipses that are consistent with spherical alignment are greyish. At the midplane it is however not possible to distinguish between spherical and cylindrical alignment, since both $\gamma_{\textrm{sph}}=\gamma_{\textrm{cyl}}=0^{\circ}$ at $z=0$ kpc, thus here consistency with spherical alignment also implies consistency with cylindrical alignment. Only away from the midplane it is possible to differentiate between these types of alignment.

We further add contours of constant formal statistical error values on the recovered tilt angles in Fig. 2. We have drawn contours for errors reaching $0.5$, $1.0$, $2.0$, and $4.0$ degrees. These contours show the great quality of our dataset over the distance range explored.

From this figure it is evident that there are just a few bins that have tilt angles much steeper than spherical alignment (i.e. there are just two dark red ellipses). These are however located in the inner regions of the Galaxy and at those positions where the error on the tilt angle is also large.

In general, however, the following trend is apparent: for Galactocentric spherical radius $r\sim 4$ kpc, the orientations of the velocity ellipses seem to be slightly steeper than spherical alignment. For $r\sim 7$ kpc they seem fully consistent with spherical alignment. For larger radii, that is $R>8$ kpc and $|z|\gtrsim 1$ kpc, the ellipses have a negative misalignment, meaning that the orientations of the ellipses become shallower compared to prediction for spherical alignment. Here the orientation thus changes into the direction of cylindrical alignment and is no longer consistent with spherical alignment.

To be able to assess whether the tilt angles found are more consistent with spherical or cylindrical alignment we show them with error bars in Fig. 3 as a function of height for four Galactic radii, namely $R=[5,7,9,11]$ kpc. The red squares (without error bars; labelled ‘Raw data’) follow from computing the moments directly from the data, and the green diamonds (‘Analytic’) and blue crosses (‘MCMC’) are derived using Method 1 and Method 2 respectively, thus accounting for the measurement errors (see Sect. 3.2). They give consistent results given the error bars, although the MCMC-method seems to result in slightly steeper tilt angles.

The black curve in Fig. 3 shows the expectation in the case of spherical alignment. At $R=5$ kpc (left panel) the recovered tilt angles are in agreement with spherical alignment for the heights explored. At $R=7$ kpc (left centre panel) the data is consistent with spherical alignment up to $|z|\sim 2$ kpc. For larger heights the tilt angles are only mildly shallower. For $R=9$ kpc and $R=11$ kpc, however, the tilt angles are becoming increasingly shallower with respect to spherical alignment. In fact, for $R=12$ kpc (see Fig. 2) the orientation of the ellipses become more consistent with cylindrical alignment for the heights probed.

4.2 Tilt angles for different azimuthal angles

Since the Galaxy is not axisymmetric we now investigate whether the tilt angles vary with azimuth by taking into account the 3D location of the individual stars in our dataset. We bin the data into Cartesian bins $(x,y,z)$ whose volume is fixed to 1$\times$1$\times$1 kpc3, which implies that the different azimuthal cones we explore contain independent data for $R>4$ kpc. These cones are centred on three different angles $\phi=[165^{\circ},180^{\circ},195^{\circ}]$.

The resulting maps are shown in Fig. 4. Since the data is effectively sliced in $\phi$, the number of stars at a given $(R,z)$ is lower and as a consequence the spatial bins cover a lower spatial extent in comparison to Sect. 4.1. A coarse comparison of the different panels in this figure suggests that the variations with azimuth are relatively small compared to the global trend that is still apparent in each panel: the misalignment changes from positive to negative when moving outwards in Galactic radius.

The most prominent differences are seen for the bins at $R\sim 4$ kpc and $z\sim 1$ kpc. The $\phi=180^{\circ}$-slice indicates much steeper tilt angles than the $\phi=195^{\circ}$-slice. The statistical errors on these tilt angles are however large. In fact, most of these bins have consistent tilt angles given their error bars.

For a more direct comparison we show in Fig. 5, for specific radii $R=[6,8,10]$ kpc, the tilt angles for the different Galactic azimuths as a function of Galactic height. Here the different symbols, namely red squares, green diamonds, and blue crosses correspond to the measurements for $\phi=[165^{\circ},180^{\circ},195^{\circ}]$, respectively. The black starred symbols show the measurements from all stars at the given $R$ and $z$ and irrespective of azimuth (as in Sect. 4.1). At $R=10$ kpc the tilt angles for the different azimuths are less consistent with spherical alignment than those at at $R=6$ kpc, especially at positive Galactic heights.

Even though some bins reveal slight differences in the tilt angles when varying Galactic azimuth, the overall qualitative trends are similar to the case in which we projected all stars onto the $(R,z)$-plane, thus justifying the approach used in Sect. 4.1. These results also suggest that the degree of non-axisymmetry, in terms of the tilt angles, is modest over the azimuthal range explored.

4.3 Variations with stellar populations

In this section we explore whether different populations of stars follow similar trends in tilt angle. To this end we have cross matched the full Gaia DR2 catalogue with three spectroscopic datasets: the Large Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST DR4, Cui et al. 2012), RAVE DR5, and APOGEE DR14. If a star has radial velocity measurements from more than one survey, we take the measurement with the smallest quoted error. As for the spectroscopic sample delivered as part of Gaia DR2 (Arenou et al. 2018), we only consider stars whose radial velocity errors have been estimated to be smaller than $20$ km/s. By adding radial velocities from these other surveys the number of stars with full phase-space information is increased by over $30\%$.

To explore dependences on populations, we only use metallicities from LAMOST DR4 since this survey probes a much larger region than either RAVE or APOGEE. We refrain from merging the metallicity information from the different surveys to avoid possible offsets between metallicity scales. Finally, only stars with metallicity uncertainties up to $0.2$ dex are considered in our analysis.

A downside of extending our sample is that Bayesian distances are missing for the newly added stars to our sample. Since the purpose of this section is to inspect variations between different populations, we here approximate the distances to the stars by $\hat{d}=1/\hat{\varpi}$, where

[TABLE]

For the following analysis, we select those stars with at most 20% relative distance errors, that is $\hat{\varpi}/\hat{\epsilon}_{\varpi}>5$, and $\hat{d}<5$ kpc. We proceed to classify the stars according to a halo population as those with $\mathrm{[M/H]}<-1.0$ dex, a thick disk population for $-1.0<\mathrm{[M/H]}<-0.5$ dex, and a thin disk population for $\mathrm{[M/H]}>-0.4$ dex. With these criteria, our sample contains $\sim 23,000$ halo stars, $\sim 260,000$ thick disk stars, and $\sim 2$ million thin disk stars.

Fig. 6 shows the velocity ellipsoids and tilt angles as a function of position in the meridional plane for the halo (left), thick disk (middle), and thin disk (right) subsamples. The different spatial coverage of the subsets reflect differences in the number of stars (recall that to reliably measure a tilt angle we require at least $100$ stars in a spatial bin). In addition the ellipses for the halo population are much larger compared to those of the thick and thin disks. In fact, we have had to use different scales for the panels: the insets in the bottom right of each panel show ellipses whose semi-major and semi-minor axes correspond to dispersions of $\sigma(v_{R})=200$ km/s and $\sigma(v_{z})=100$ km/s for the halo and thick disk populations, and to $\sigma(v_{R})=100$ km/s and $\sigma(v_{z})=50$ km/s for the thin disk.

As in previous sections, the colours in Fig. 6 represent the misalignment of the tilt angles with respect to spherical alignment. The same trends as found earlier are visible for the populations independently: at $R\lesssim 7$ kpc the alignment is closer to spherical, while outwards from $R\sim 9$ kpc the misalignment becomes negative, which means that the tilt angles become shallower. This can be seen more easily when comparing the tilt angles derived for each population at specific radii, as shown in Fig. 7.

There are also some differences seen. For example, at $R=8.5$ kpc, the halo sample seems to be more consistent with spherical alignment than both disk samples. For $R=9.5$ kpc, however, the differences between the populations are minor, except for the flatter thin disk tilt angles at $z\sim 2.5$ kpc. Therefore we may conclude that the results shown in Sect. 4.1 are not strongly dependent on the different populations present throughout the volume probed by our dataset.

4.4 Quantifying the degree of spherical alignment

Because the trends seen in the tilt angles are not strongly dependent on Galactic azimuth nor on stellar population, we here aim to provide a simple description of their variation with radius $R$ and height $z$ as found in Sect. 4.1. Since we infer near spherical alignment for $R\sim 6$ kpc, we consider expanding $\alpha$ around a point $(R_{0},z_{0})$:

[TABLE]

where $a_{i}$ are constants and both $R$ and $z$ in kpc222We prefer to quantify the deviation from spherical symmetry directly on the spherical tilt angle $\alpha$ than to use the purely geometric parametrisation by Binney et al. (2014) of the cylindrical tilt angle $\gamma\prime=a_{0}\arctan(z/R)=a_{0}(\pi/2-\theta)$ where $\theta$ indicates the spatial location of the bin (see also Eq. 3). Although $a_{0}=1$ implies spherical alignment and $a_{0}=0$ cylindrical alignment, it is not intuitively clear what the quantitive meaning of other $a_{0}$ values is.. By definition $\alpha(R_{0},z_{0})=0^{\circ}$. We further set $z_{0}=0$ kpc (i.e. the symmetry plane of $\alpha$ is set to be the Galactic midplane). Moreover, $a_{1}=a_{4}=0$, since for most realistic models the tilt angle does not vary at the midplane. By symmetry arguments the coefficients of all even powers of $z$ (including $a_{5}$) must be zero, since $\alpha$ is expected to be either antisymmetric with respect to the midplane or zero. Since we have found that at $R\sim 6$ kpc the tilt angles are consistent with spherical alignment for all $z$ probed (see left panels of Fig. 3), we additionally set $a_{2}=0$ such that at $R=R_{0}$: $\alpha(R_{0},z)=0^{\circ}$. With these choices:

[TABLE]

We thus fit this functional form to the data to derive values for $R_{0}$ and $a_{3}$ such that the $\chi^{2}$-statistic defined as:

[TABLE]

is minimised. Here $j$ runs over the number of bins $N_{\mathrm{bins}}$ where a measurement is made, in other words where $N>100$ stars.

For most bins at $|z|\leq 2.0$ kpc and $5\leq R\leq 12$ kpc the inferred statistical errors on the tilt angles are very small (e.g. see the dashed contours in Fig. 2). In that case systematic errors need to be considered. One such source of systematic errors are substructures. We performed tests to estimate the effect of substructures in velocity space on the tilt angle. To this end we inserted $N_{\textrm{sub}}=1$, $4$, $9$, $16$, $25$, or $36$ substructures on smooth non-tilted velocity distributions with velocity dispersions of $20$ km/s and $35$ km/s in $v_{z}$ and $v_{R}$ (i.e. values representative of the thin disk near $R\sim R_{\odot}$), respectively. Each substructure was assigned a random number of stars such that the total fraction of stars in substructures is $f_{\textrm{sub}}=5\%$, $10\%$, $15\%$, or $20\%$. We randomly assigned velocity dispersions to the substructures, drawn uniformly from $1$ km/s to $5$ km/s in both directions. For each combination of $(N_{\textrm{sub}}$, $f_{\textrm{sub}})$ we considered $100$ realisations. The median (absolute) tilt angle found from these experiments is $\sim 1$ degree, implying that this value is representative of the error introduced by neglecting the presence of substructures in a velocity distribution. This result is independent of the total number of stars $N$ for $N\gtrsim 10,000$ (a value that is representative of the number of stars in the bins with $\epsilon[\alpha(R_{j},z_{j})]<1^{\circ}$). Thus, when minimising the $\chi^{2}$ we consider a floor for the statistical error $\epsilon[\alpha(R_{j},z_{j})]$ in each bin of $1^{\circ}$.

We fit to find $R_{0}=(6.16\pm 0.16)$ kpc and $a_{3}=(0.72\pm 0.04)^{\circ}/{\mathrm{kpc}}^{2}$ resulting in a reduced $\chi^{2}$ of $1.65$. The cyan line in Fig. 3 shows the tilt angles predicted by this fit, which reproduces relatively well the trends observed in the data. The model goes through the $1\sigma$-error bars for approximately $60\%$ of all spatial bins, while for $98\%$ of bins the model matches the data within 3$\times$ the estimated uncertainty. This indicates that our simple model provides a fair description of the behaviour of the tilt of the velocity ellipsoid across the Galactic volume probed by our dataset.

The fact that the total reduced $\chi^{2}$-value is greater than unity indicates that the tilt angles for some bins are not fitted very well by the model. For example at $R\sim 10$ kpc the tilt angles as inferred from the data are asymmetric with respect to the $z=0$ plane: at $z>0$ kpc they more or less attain a constant value of $\sim 2.0^{\circ}$, whereas below the midplane the tilt angles become steeper with $z$ (e.g. $-15^{\circ}$ at $z=-3.0$ kpc). The fits at such radii are therefore relatively poor. For the bins between $R=11$ kpc and $R=12$ kpc, we notice that the observed tilt angles seem to have a small positive offset from zero near $z=0$. These offsets are small (of order $2$ degrees), although they do affect the goodness of fit measure.

5 Discussion

5.1 The impact of (parallax) measurement errors on the recovered tilt angles.

Gaia Collaboration et al. (2018a) have reported the presence of a systematic error on the Gaia DR2 parallaxes in the form of a zero-point offset of a few 10 of $\mu$as (in the sense that Gaia parallaxes are too small) and whose exact amplitude depends on location on the sky. Such systematic zero-point offset affects the tangential velocities of the stars, which are determined from both distances and proper motions. The overall systematic parallax offset in Gaia DR2 was determined using distant quasars by Lindegren et al. (2018) to be approximately $-29\mu$as, with a large RMS of $\sim 43\mu$as. Arenou et al. (2018) using different samples of objects (RR Lyrae stars, Magellanic Clouds, open clusters, dwarf spheroidal galaxies, etc.) report important variations in the zero-point offsets, highlighting the complexity of the offset. Nonetheless all values are consistent given the large estimated RMS.

Around the time this paper was submitted, Schönrich et al. (2019) reported a new estimate of the parallax zero-point offset based on the distance estimation method used in Schönrich & Aumer (2017) (also see Schönrich et al. 2012). These authors argue for a much larger zero-point for the parallaxes in the RVS subset of Gaia DR2, namely of magnitude $-54\pm 6\mu$as. Zinn et al. (2018) and Khan et al. (2019) applied asteroseismology to determine distances to Red Giant Branch (RGB) and Red Clump (RC) stars with Gaia $G$-band magnitudes similar to those present in the RVS subset of Gaia DR2 and determined an offset close to $-50\mu$as, while Sahlholdt & Silva Aguirre (2018), using asteroseismology information on dwarfs, report that the offset could be $\sim-35\pm 16\mu$as. More recently Hall et al. (2019), using RC stars with asteroseismology, estimate the mean offset to be $-41\pm 10\mu$as. These comparisons suggest that the offset could well be larger for the brighter stars of the Gaia RVS sample but that its amplitude is quite uncertain.

5.1.1 Quantification of the impact of a zero-point offset

We first quantify how the tilt angles are affected if parallaxes are underestimated. For illustration purposes, we estimate the impact on the recovered tilt angles induced by a systematic error (with mean $-29\mu$as) while also including the effects of random errors333In Appendix B we analytically compute how the $v_{R}$- and $v_{Z}$-velocities (and thus their moments and tilt angles) are affected by the parallax zero-point offset alone.. Their effect is examined by using the Gaia Universe Model Snapshot (GUMS, Robin et al. 2012), which is based on the Besançon Galaxy Model (Robin et al. 2003).

We mimic the Gaia DR2 subsample with full phase-space information, by selecting stars in GUMS that have $G<13$ mag, as this is roughly the magnitude limit for radial velocities in Gaia’s current data release. We generate $100$ data realisations by convolving the (error-free) GUMS sample with a Gaussian with Gaia DR2-like random and systematic errors for the parallaxes (Lindegren et al. 2018). The systematic parallax offsets for the stars are drawn from a Gaussian with mean $-29\mu$as and standard deviation of $30\mu$as444as estimated in Gaia Collaboration et al. (2018b).. To obtain a distance estimate we invert the parallaxes and consider only those stars that satisfy $\varpi/\epsilon(\varpi)>5$ and $\varpi\gtrsim 200\mu$as. Here $\varpi$ is the observed parallax and $\epsilon(\varpi)$ the random parallax error and thus the same quality criteria are applied as to the real data (see Sect. 2).

For each spatial bin the median (over all realisations) of the distribution of tilt angles is compared to the tilt angles from the error-free model, on the meridional plane. The error-free GUMS model has close to cylindrically aligned velocity ellipses ($\gamma_{\textrm{GUMS}}\sim 0^{\circ}$). The impact of the random and systematic parallax uncertainties on the tilt angles depends on location as can be seen in the left panel of Fig. 8. At $R\lesssim 7$ kpc the orientations of the velocity ellipses change towards the direction of spherical alignment ($\Delta\gamma_{\textrm{GUMS}}>0$ for $z>0$ and $\Delta\gamma_{\textrm{GUMS}}<0$ for $z<0$), while for $R\gtrsim 9$ kpc the change is in the opposite sense.

The middle panel of Fig. 8 shows the spread in tilt angles over all realisations, and reveals that the errors result in a spread with a typical amplitude of $\lesssim 4^{\circ}$, except for the outermost bins, where it can be twice as large, and hence comparable to $\Delta\gamma_{\textrm{GUMS}}$. The right panel shows at which locations the median change in tilt angles, caused by parallax errors, is larger than the RMS from realisation to realisation. For bins located at distances of $\sim 2$ kpc a change in tilt angle due to parallax errors is thus likely to occur in a preferential direction, with the amplitude of this change varying from realisation to realisation.

These findings imply that, if parallaxes are underestimated, the tilt angles inferred may appear steeper than they really are in the inner Galaxy, while the opposite happens in the outer Galaxy, thus the tilt angles become shallower there. If we take the results from GUMS at face value, $|\Delta\gamma_{\textrm{GUMS}}|\approx 6^{\circ}$ at $(R,|z|)\sim(5,3)$ kpc, which means that an unaccounted for zero-point offset of magnitude $29\mu$as in the parallaxes affects the inferred tilt angles such that they appear steeper by $\sim 6^{\circ}$. This does not radically change the type of alignment at this location (where spherical alignment would imply $\gamma\sim 30^{\circ}$). For $(R,|z|)\sim(11,2)$ kpc we find that $|\Delta\gamma_{\textrm{GUMS}}|$ can attain values close to $5^{\circ}$, which is of similar amplitude as the misalignment seen in Fig. 2. Although the GUMS tilt angles intrinsically have $\gamma_{\textrm{GUMS}}\sim 0^{\circ}$, we find similar amplitudes for the cases explored in Appendix B, where we start from both intrinsically spherically and cylindrically aligned ellipsoids.

In the analysis presented in previous sections, we have effectively corrected for the parallax offset by using the McMillan (2018) distances. If the assumed parallax zero-point is too small, the results presented in this section indicate that, especially towards the outer Galaxy, the zero-point offset could produce tilts that are less steep than what they are intrinsically. We explore such a larger offset next.

5.1.2 A zero-point offset as large as -54$\mu$as

Everall et al. (2019) have derived tilt angles using the Schönrich et al. (2019) Bayesian distance estimates (with parallax zero-point of $-54\mu$as). These authors showed that the tilt angles appear to be much more consistent with spherical alignment when using those distances.

Since the method used in Schönrich et al. (2019) assumes spherical alignment, we preferred not to directly use their distances while testing for the effect of a large $-54\mu$as offset. Therefore we here also explore how the tilt angles change if the parallax offset would be as large as $-54\mu$as, by comparing them to the case in which the offset is $-29\mu$as. For both cases we take the extended sample and invert the parallaxes after correcting for the zero-point offset (as in Sect. 4.3), such that the changes due to the differences in parallax offset can be easily compared.

In Fig. 9 we show the results. The blue squares have been calculated after correcting for a parallax zero-point offset of $-29\mu$as, whereas for the orange diamonds a value of $-54\mu$as is assumed. For the outer Galaxy ($R=9$ kpc and $R=11$ kpc) such a larger parallax zero-point can modify the tilt angles such that they are more consistent with spherical alignment, in agreement with our analysis of the previous section.

A direct comparison of the tilt angles obtained using McMillan (2018) Bayesian distances (who assumes a zero-point of $-29\mu$as, green crosses) with the results obtained from inverting the parallaxes after correcting for a zero-point of $-29\mu$as (blue squares), shows good agreement except for $R=5$ kpc. At this location, it would seem as if the choice of the distance estimator would play a role in the determination of the tilt angle. The Bayesian distances result in tilt angles that are just slightly steeper than expected for spherical alignment, while inverting the parallaxes results in much shallower tilt angles (the larger the offset assumed the shallower the tilt angles). On the other hand, comparing the tilt angles obtained using Schönrich et al. (2019) Bayesian distances (who find a zero-point of $-54\mu$as, red starred symbols), with the results obtained from inverting the parallaxes after correcting for a zero-point of $-54\mu$as (orange diamonds), shows rather similar trends at $R=5$ kpc. At the other radii shown, these Bayesian distances also result in tilt angles that are in good agreement with inverting the parallaxes.

The analysis presented in the last two sections shows that the amplitude of the systematic error in the parallax, in the form of a zero-point offset, plays a role in the determination of the tilt angles for the outer Galaxy ($R>9$ kpc). Since the offset is known to vary with celestial position, magnitude and colour, it is difficult at this point to properly correct for it, and this impairs a very accurate determination of the tilt angle throughout the range of distances probed. However, recall that the range of zero-point offsets is bracketed by the values explored (i.e. from $-54\mu$as to $-29\mu$as), so the analysis presented here gives us a handle on the possible outcomes.

5.2 Constraints to models of the Milky Way

Several models of the Milky Way have been proposed by matching a variety of constraints (e.g. McMillan 2011; Piffl et al. 2014; McMillan 2017). Particularly useful for the interpretation of the findings reported in this paper are Stäckel models (e.g. de Zeeuw 1985; Dejonghe & de Zeeuw 1988). Axisymmetric models with a potential of Stäckel form have the property that the equations of motion are separable in their spheroidal coordinates. Therefore the principal axes of the velocity ellipsoids are always aligned with these coordinates (also see: Eddington 1915). The foci of such a coordinate system then determine the alignment at each position. For a composite model to be of a Stäckel form, the locations of the foci must be identical for all components.

Famaey & Dejonghe (2003), for example, have extended the two-Stäckel component work of Batsleer & Dejonghe (1994) by adding a third component, such that the model could allow for a thin and thick disk, in addition to a halo component. The authors use constraints such as the (flat) rotation curve, circular velocity at the position of the Sun, the Oort constants, and the local total mass density in the disk to search for a set of consistent parameters for their Stäckel models. Here we take the set of prolate spheroidal coordinates, $(\lambda,\phi,\nu)$, from Famaey & Dejonghe (2003, mass model II). The foci of this oblate mass model are located at ($R$, $z$) = ([math], $\pm 0.88$) kpc. At $R\sim 0$ and $|z|\lesssim 0.88$ kpc such spheroidal coordinates align with the cylindrical coordinate system (see Fig. 10). Outside of these foci and with increasing distance from the Galactic centre the spheroidal coordinates approach the spherical coordinate system. In general, any (composite) Stäckel model predicts a change in the tilt of the velocity ellipse from cylindrical to spherical alignment. The transition radius depends on the location of the foci.

Since the observed tilt angles at $R\sim 4$ kpc already show near spherical alignment, this implies foci at $|z|\lesssim 4$ kpc. Their exact position would depend on whether the innermost region of the Galaxy, not probed by our dataset, is cylindrically aligned or not, and if so at what distance the transition occurs. However, the tilt angles in the outer Galaxy ($9\lesssim R\lesssim 12$ kpc) derived using the McMillan (2018) distances are not consistent with Stäckel models that have foci at $|z|\lesssim 4$ kpc, and would require a larger focal distance. We have numerically checked these statements by comparing the predicted tilt angles of both oblate and prolate Stäckel models (for a large range of different focal distances) to the observed tilt angles while taking into account their errors.

There are of course many more models with bulge, disk and halo components, for example spherical bulge, exponential disk, Navarro-Frenk-White (NFW, Navarro et al. 1996) halo, or Miyamoto & Nagai (1975) models. The separable models are in that sense a subset but have the advantage that for them the tilt of the velocity ellipsoid is dictated by the coordinate system in which the equations of motion (Hamilton-Jacobi equation to be more precise) separate.

Piffl et al. (2014) have applied a five component mass model (gas disk, thin and thick disk, flattened bulge and dark halo) to RAVE DR4 stars. Using their best-fitting parameters we computed the relevant velocity moments from the distribution function for a similar range in $R$ and $z$ as probed in our dataset. The tilt angles for this model are spherically aligned for $R\gtrsim 7$ kpc and are, as in the separable models discussed above, changing towards cylindrical alignment with decreasing $R$.

In Fig. 11 we show the tilt angles for both the Stäckel model (purple line) of Famaey & Dejonghe (2003) and the Piffl et al. (2014) model (orange line), for radii $R=6$ kpc and $R=10$ kpc. The green diamonds indicate the tilt angles as found by Method 1. Since this Stäckel model has focii at $|z|\lesssim 0.88$, which is very close to the Galactic centre with respect to the innermost radius probed in our dataset, the Stäckel model is almost indistinguishable from spherical alignment for all positions probed. The Piffl et al. (2014) model has tilt angles that are shallower at $R=6$ kpc, but also approaches spherical alignment with increasing Galactic radius. At $R=10$ kpc, for example, the tilt angles from the Piffl et al. (2014) model are seen to nearly coincide with the expectation for spherical alignment.

We note that if the parallax zero-point is larger than assumed here the tilt angles do become more consistent with spherical alignment for large radii (see 5.1.2). This is in line with predictions for both composite Stäckel models as well as for the Piffl et al. (2014) model. In addition, it would be interesting to know whether the tilt angles become shallower towards the central regions of the Galaxy (at $R\lesssim 4$ kpc). In principal it would then be possible to solve for the focal distance. However, the effects of both the type of distance estimator and the assumed parallax zero-point are too large to make firm statements in this region. Future data releases will for sure enable to probe regions closer to the Galactic centre more robustly.

6 Conclusions

We have studied the trends in the tilt angle of the velocity ellipsoids in the meridional plane for a high-quality sample of more than $5$ million stars located across a large portion of the Galaxy, from $R\sim 4$ kpc to $R\sim 13$ kpc, and reaching a maximum distance from the plane of $\sim 3.5$ kpc.

We find that the tilt angles are somewhat dependent on the offset of the Gaia DR2 parallaxes, and that the effects are particularly important for the outer Galaxy. When using the McMillan (2018) Bayesian distances, derived assuming an offset of $-29\mu$as, we find that the tilt angles are consistent with (near) spherical alignment at $R\lesssim 7$ kpc for all heights probed ($|z|\lesssim 3$ kpc). Beyond $R\gtrsim 9$ kpc the tilt angles clearly become more shallower than expected for spherical alignment. These trends remain when the stars are separated into ‘populations’ according to their metallicity (as given by LAMOST DR4). We provide a simple analytic function for the tilt angle in spherical coordinates $\alpha(R,z)/[{\rm deg}]\approx 0.72(R-6.16)z$, that fits well the trend observed as a function of Galactic radius and height, after projecting the stars onto the $(R,z)$-plane.

We find that if the amplitude of the zero-point offset in the parallax is underestimated, the angles tend to appear shallower than they intrinsically are in the outer Galaxy (i.e. changing into the direction of cylindrical alignment if the ellipsoid is intrinsically spherically aligned). We quantify the impact on the tilt angles when assuming a parallax zero-point as large as $-54\mu$as, as estimated in Schönrich et al. (2019) (also see Everall et al. 2019). Such a large offset (the upper limit of estimates reported in the literature by other authors) does indeed lead to tilt angles that are more consistent with spherical alignment than obtained when using the McMillan (2018) distances. Therefore it will be particularly important to pin-down, in future Gaia data releases, the amplitude of the parallax zero-point as well as its local variations as these affect our ability to constrain the mass distribution in our Galaxy.

Acknowledgements.

The authors thank the anonymous referee whose insightful comments helped improving the quality of the manuscript. JH thanks Helmer Koppelman for helping in creating the dataset used in this work (Koppelman et al. 2018). AH acknowledges financial support from a VICI grant from the Netherlands Organisation for Scientific Research, N.W.O. TdZ is grateful to the Kapteyn Astronomical Institute for the hospitality during his Blaauw Professorship. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. We have also made use of data from: (1) the APOGEE survey, which is part of Sloan Digital Sky Survey IV. SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration (http://www.sdss.org). (2) the RAVE survey (http://www.rave-survey.org), whose funding has been provided by institutions of the RAVE participants and by their national funding agencies. (3) the LAMOST survey (www.lamost.org), funded by the National Development and Reform Commission. LAMOST is operated and managed by the National Astronomical Observatories, Chinese Academy of Sciences. For the analysis, the following software packages have been used: vaex (Breddels & Veljanoski 2018), NumPy (Oliphant 2015), matplotlib (Hunter 2007), Jupyter Notebook (Kluyver et al. 2016), TOPCAT and STILTS (Taylor 2005, 2006).

Appendix A Standard errors of sample (co)variances.

To estimate the error on the inferred tilt angles from Method 1 of Sect. 3.2 we propagate the errors of the relevant velocity moments from Eq. 1.

The error on a sample variance, $s^{2}$, can be estimated (e.g. Rao 1973; Mood et al. 1974) by using

[TABLE]

for $N$ stars. Here, $\mu_{4}$ denotes the intrinsic $4^{\mathrm{th}}$ central moment and $s^{2}=\frac{1}{N-1}\sum^{N}_{i=1}\left(v_{i}-\langle v\rangle\right)^{2}$, for which $v_{i}$ is the relevant velocity component, either $v_{R}$ or $v_{z}$, of star $i$ and $\langle v\rangle$ its mean taken over all stars in the bin considered. The intrinsic velocity moments are estimated by their observed values, which is a good approximation given the relatively small errors in the data for the bins explored.

The error on a sample covariance $S_{xy}$ of $x$ and $y$ can be estimated (see Stuart & Ord (1987) or Rose & Smith (2002) for using mathStatica) by

[TABLE]

where $\mu_{22}=E[\{x-E(x)\}^{2}\{y-E(y)\}^{2}]$ for $E$ denoting the expectation value. We have defined $S_{xy}=\frac{1}{N-1}\sum^{N}_{i=1}(x_{i}-\langle x\rangle)\,(y_{i}-\langle y\rangle)$. In our application $x$ is replaced for $v_{R}$ and $y$ for $v_{z}$. The intrinsic moments are again estimated by taking the equivalent moments directly from the observed velocity distribution.

As an example for Sect. 3.2 we show in Fig. 12 the error distributions for the bin at $R=11.5$ kpc and $z=1.5$ kpc. This bin is near the edge of the volume investigated, but still contains $2,016$ stars. The vertical grey dashed lines indicate the values of the velocity moments that would be derived by using the data directly (i.e. not accounting for the errors). The medians of the error distributions are indicated by the vertical grey dotted lines. The recovered intrinsic velocity moments from Method 1 are visualised by the vertical black solid lines (as here, these usually coincide with the vertical grey dashed lines). Thus, even for this outer bin, the effects of measurement errors are relatively small.

Appendix B The impact of a systematic parallax offset on the recovered tilt angles.

Here we explain how a systematic parallax offset can affect the inferred tilt angles. For this purpose, we now only consider the $(x,z)$-plane and we assume that all parallaxes are shifted by the same offset $\Delta\varpi=-0.029$ mas.

For Galactic longitude $l$ and latitude $b$ the $(U,V,W)$-velocities in km/s can be computed the usual way (Johnson & Soderblom 1987; Bovy 2011):

[TABLE]

Here, $\mu_{l^{\star}}=\mu_{l}\cos(b)$ and $\mu_{b}$ denote the proper motions in mas/yr in the direction of $l$ and $b$, respectively, $\varpi$ is the parallax in mas, and $k=4.74047\frac{\textrm{km/s}}{\textrm{kpc mas/yr}}$ (assuming a Julian year).

When only considering an error in the parallaxes the ‘observed’ velocities are affected as:

[TABLE]

Subscript 0 denotes the true position and velocities, subscript 1 the ‘observed’ quantities. Furthermore:

[TABLE]

and:

[TABLE]

Let us now define the tilt angle $\delta$ as:

[TABLE]

In a steady state axisymmetric system $\langle v_{R}\rangle=\langle v_{z}\rangle=0$. Therefore, at the $(x,z)$-plane $\langle U\rangle=\langle W\rangle=0$, and thus $\textrm{var}(U)=\langle U^{2}\rangle$, $\textrm{var}(W)=\langle W^{2}\rangle$, and $\textrm{cov}(U,W)=\langle UW\rangle$. For $l=0^{\circ}$ and $l=180^{\circ}$ we also notice that $U=-v_{R}$ and $W=v_{z}$ , and therefore that $\delta=-\gamma$. In the remainder of this Appendix we refer to $\delta$ when we use ‘tilt angle’ (unless stated otherwise).

Plugging Eq. 13 up to first order in $\frac{\Delta\varpi}{\varpi_{0}}$ into Eq. 16 we get:

[TABLE]

in which:

[TABLE]

where $\pm$ holds for $l\in\begin{Bmatrix}0^{\circ}\\ 180^{\circ}\end{Bmatrix}$.

To further explore the effect of a shift in the parallaxes we now investigate what would happen to the tilt angles in two different cases of alignment: spherical alignment and cylindrical alignment.

We start by rewriting Eq. 17 into the form of

[TABLE]

We then get:

[TABLE]

in which:

[TABLE]

Then, under the assumptions that $\left|\epsilon_{C}\right|\ll 1$ and $\left|\epsilon_{D}\right|\ll 1$, we get:

[TABLE]

We highlight the effects for four different latitudes:

[TABLE]

Since the velocity ellipse is mostly non-tilted ($\delta_{0}=0^{\circ}$) at the Galactic midplane the inferred tilt angles at $b=0^{\circ}$ are not affected by an error in the parallax. Geometrically this is not surprising since, at $b=0^{\circ}$, the $U$-component of the velocities are not affected. The $W$-velocities are only inflated and do not change the tilt angle. However, if $\delta_{0}\neq 0^{\circ}$, then the term between the square brackets becomes larger than one, since for typical values of the velocity moments at the midplane $\sigma(v_{R})>\sigma(v_{z})$ (see e.g. Gaia Collaboration et al. 2018c). The inferred tilt angle is therefore steeper (more positive if $\delta_{0}>0^{\circ}$ and more negative if $\delta_{0}<0^{\circ}$). At $|b|=90^{\circ}$, the effect is reversed and the tilt angle becomes shallower (less positive if $\delta_{0}>0^{\circ}$ and less negative if $\delta_{0}<0^{\circ}$) due to the parallax offset. For the case of spherical alignment the relation $\tan(2\delta_{0})=\tan(2\theta)$ can be applied.

The approximations used so far fail for $\langle U_{0}W_{0}\rangle\simeq 0$, since then $|\epsilon_{\mathrm{C}}|\not\ll 1$, and for $\langle U^{2}_{0}\rangle\simeq\langle W^{2}_{0}\rangle$, since then $|\epsilon_{\mathrm{D}}|\not\ll 1$, and thus $|x|\not\ll 1$. In the case of cylindrical alignment ($\delta_{0}=\langle U_{0}W_{0}\rangle=0$) and for $\left|\epsilon_{\mathrm{D}}\right|\ll 1$ we get555If, hypothetically, both $\langle U^{2}_{0}\rangle=\langle W^{2}_{0}\rangle$ and $\langle U_{0}W_{0}\rangle=0$, then $\tan(2\delta_{1})=\pm\tan(2b)$.:

[TABLE]

where we used that $\tan(2\delta_{1})\simeq 2\delta_{1}$ for small deviations around $\delta_{1}=0^{\circ}$. This means that at $l=0$ ($l=180^{\circ}$) and for $\sigma(v_{R})>\sigma(v_{z})$ the tilt angles appear to be negative (positive) for $b>0^{\circ}$, and positive (negative) for $b<0^{\circ}$.

We have inserted the relevant Galactic velocity dispersions as a function of $R$ and $z$ and set the covariance term such that there is either spherical or cylindrical alignment throughout the extent of the dataset. We find that the tilt angles are affected very similarly. This is visualised in Fig. 13 (recall that $\gamma=-\delta$ since we here consider $l=0^{\circ}$ and $l=180^{\circ}$ only). We therefore think that our test performed in Sect. 5.1 is realistic, even though the intrinsic tilt angles of the GUMS catalogue are more or less cylindrically aligned.

Besides the fact that the orientation of the velocity ellipse changes due to the parallax offset, obviously the stars under consideration also move in position. Thus, in fact a sample of stars with tilt angle $\delta_{0}$ at parallax $\varpi_{0}$ gets ‘observed’ at $\varpi_{1}$ with tilt angle $\delta_{1}$. We have not taken this effect into account in the analytic description from this Appendix.