DETECTION OF A DEARTH OF STARS WITH ZERO ANGULAR MOMENTUM IN THE SOLAR NEIGHBORHOOD

We are entering a golden age for Milky Way astronomy. The European Space Agency's (ESA) Gaia mission (Gaia Collaboration et al. 2016b), which launched 2013 December 19, has recently published its first data release (Gaia Collaboration et al. 2016a), giving us a new window into our Galaxy and, in particular, the solar neighborhood. The primary astrometric catalog in Gaia-DR1 is the Tycho-Gaia Astrometric Solution (TGAS; Michalik et al. 2015; Lindegren et al. 2016), which uses data from the Tycho-2 catalog (Høg et al. 2000), to provide a baseline of approximately 30 years upon which to calculate astrometric values for stars in common between Tycho-2 and Gaia. There is also significant overlap between stars in the TGAS catalog and stars observed by the Radial Velocity Experiment (RAVE; e.g., Steinmetz et al. 2006), enabling the full six-dimensional phase-space information to be known for over 200,000 stars in the solar neighborhood. This enables us to explore local dynamics in unprecedented detail.

Many aspects of the structure and dynamics of the Milky Way are difficult to measure, owing to our position within the Galaxy, and complex observational selection effects such as dust extinction. Thus, many of the fundamental parameters of the Milky Way carry significant uncertainty. One such parameter is the velocity at which the Sun rotates around the Galaxy, for which plausible values span the range from $\approx 240$ to $260\,\mathrm{km}\ {{\rm{s}}}^{-1}$ (e.g., Bovy et al. 2012; Reid et al. 2014) and are dependent on the data set and technique used.

Carlberg & Innanen (1987) proposed that the rotation velocity of the Sun may be measured by searching for a lack of stars exhibiting zero angular momentum.⁴ Stars with zero angular momenta are expected to plunge into the Galactic nucleus and subsequently experience scattering onto chaotic orbits with a high scale height, henceforth spending the majority of their time in the stellar halo (Martinet 1974). If stars with very low angular momentum are indeed not present in the solar neighborhood, the tail of the tangential velocity distribution will exhibit a dip centered at the solar reflex value. This method for measuring the solar velocity is attractive because it should depend only on the existence of low angular momentum orbits within the Milky Way's disk. With six-dimensional phase-space measurements for nearby stars it becomes possible to calculate the tangential velocity distribution with respect to the Sun, and thus we can test this prediction by searching for a dearth of stars around the assumed value of the negative of the solar motion.

This Letter is constructed as follows. In Section 2, we discuss our treatment of the data and the feature observed in the resulting velocity distribution. In Section 3, we present simulated models that can explain the feature and make predictions for the size and shape of the observed feature. In Section 4, we fit our model to the data and present our measurement of ${v}_{\odot ,\mathrm{reflex}}$. Finally, in Section 5, we discuss the implications of the detection and look forward to future measurements.

We cross match the TGAS catalog (Lindegren et al. 2016) and the RAVE DR5 data (Kunder et al. 2016) to add RAVE line of sight velocities to the TGAS astrometric data. The resulting sample consists of 216,201 stars with six-dimensional phase-space measurements. Where available, we employ the RAVE spectrophotometric distance estimates because estimating distances from the TGAS parallax is non-trivial (e.g., Bailer-Jones 2015). However, where no RAVE distance is available, we naively invert the TGAS parallax π to obtain distance estimates for the remaining stars, removing any star with ${\sigma }_{\pi }/\pi \gt 0.1$ to avoid large distance uncertainties. Then, we convert the velocities from equatorial coordinates to standard Galactic Cartesian coordinates centered on the Sun, $(X,Y,Z,{v}_{X},{v}_{Y},{v}_{Z})$, with v_X positive in the direction of the Galactic center (toward Galactic longitude l = 0) and v_Y positive in the direction of Galactic rotation (toward $l=90^\circ$), both measured with respect to the Sun.

Figure 1 shows the distribution of v_Y for stars within 700 pc as a histogram (top panel) and with Kernel Density Estimation (e.g., Wasserman 2006; bottom panel) using a Gaussian kernel with bandwidth 9. We chose 700 pc as a balance between quantity of stars and quality of data. There is a clear underdensity of stars in the approximate region of $-210\,\mathrm{km}\ {{\rm{s}}}^{-1}$ to $-270\,\mathrm{km}\ {{\rm{s}}}^{-1}$, marked with an arrow (see also the top panel of Figure 4). Although there are only $\approx 10$ stars per bin at these velocities, this dip in the distribution is clear across $\approx 8$ bins. We defer a discussion of the significance until Section 4 where we fit a model for the dip obtained from simulations. We also confirmed that the stars in this low-velocity tail are not systematically metal-poor compared with the overall sample, for which the RAVE distances are expected to be of lower accuracy.

Zoom In Zoom Out Reset image size

Furthermore, in the two-dimensional distribution of $({v}_{X},{v}_{Y})$, we find a lack of stars with both positive v_X and large negative v_Y. That is, there are very few stars observed on disk orbits ($| {v}_{X}| \lt 150\ \mathrm{km}\ {{\rm{s}}}^{-1}$ and $| {v}_{Z}| \lt 150\,\mathrm{km}\ {{\rm{s}}}^{-1}$) plunging toward the Galactic center. Quantitatively, there are 24 stars with ${v}_{X}\gt 0$ in the range $-290\,\mathrm{km}\ {{\rm{s}}}^{-1}\lt {v}_{Y}\lt -230\,\mathrm{km}\ {{\rm{s}}}^{-1}$ and 34 stars with ${v}_{X}\lt 0$. These rates are inconsistent with being drawn from the same distribution at $\gtrsim 1.7\sigma$.

Our analysis does not explicitly take into account the uncertainties on the distance or velocity estimates. Propagating the uncertainties through the coordinate transformation results in an uncertainty of approximately 10% on the measured Cartesian velocities. At the v_Y range of the dip, this corresponds to uncertainties of approximately $24\,\mathrm{km}\ {{\rm{s}}}^{-1}$. Because the width of the dip feature appears to be approximately $60\,\mathrm{km}\ {{\rm{s}}}^{-1}$ we can neglect the uncertainties in this initial investigation.

A likely explanation for these missing stars, as mentioned in Section 1, is that they have been scattered onto chaotic orbits with larger scale heights by interaction with the Galactic nucleus. This is expected for disk stars with approximately zero angular momentum as discussed in Carlberg & Innanen (1987). Thus, as these stars would then spend the majority of their orbits far from the Galactic plane, it is very unlikely that they would be observed in the solar neighborhood at any one given time. In Galactocentric coordinates, such stars have tangential velocities ${v}_{T}\approx 0$, corresponding to heliocentric ${v}_{Y}\approx {v}_{\odot ,\mathrm{reflex}}$, minus the solar tangential velocity measured in the Galactocentric frame.

Carlberg & Innanen (1987) performed an analysis of this effect and matched models of the dip constructed within an analytic potential to data from local stellar catalogs complete to 25 pc. They found the dip to be centered at $250\,\mathrm{km}\ {{\rm{s}}}^{-1}$, with a depth greater than 80%. However, they are careful to note that there are only 18 stars with ${v}_{Y}\lt -140\,\mathrm{km}\ {{\rm{s}}}^{-1}$ and that a larger sample will improve the measurement. The TGAS+RAVE sample used here has 374 stars with $-310\,\mathrm{km}\ {{\rm{s}}}^{-1}\leqslant {v}_{Y}\lt -150\,\mathrm{km}\ {{\rm{s}}}^{-1}$, allowing far greater confidence in our subsequent analysis of the feature.

First, we make a fresh prediction of the feature that we expect to observe in the TGAS data, similar to Carlberg & Innanen (1987), but with an updated potential, and drawing the distribution of initial positions, radial, and vertical velocities from the observed data to tailor the prediction to the current data set. It is challenging to construct an N-body model with sufficient resolution to predict the high-velocity tail in the local neighborhood. Thus, we integrate test particles in a Milky Way–like potential and observe the resulting orbits.

For our Milky Way potential we use MWPotential2014 from galpy (Bovy 2015). MWPotential2014 consists of a power-law spherical bulge potential with an exponential cutoff, a Miyamoto–Nagai disk potential, and an NFW halo potential. The parameters of this potential have been fit to a wide variety of dynamical data in the Milky Way; the full parameters are given in Bovy (2015). To model the hard Galactic nucleus that is not included in MWPotential2014, we also include a Plummer potential ${\rm{\Phi }}(R,z)=-M/\sqrt{{R}^{2}+{z}^{2}+{b}^{2}}$, where $M=2\times {10}^{9}$ M_⊙ and b = 250 pc. Note that we do not include a non-axisymmetric bar potential in this initial work, which may affect the feature slightly. This should be considered when applying this technique to future Gaia data releases.

Figure 2 displays the orbits of three stars integrated for 2 Gyr. The top panel shows the orbit of a star with ${v}_{T}=0\,\mathrm{km}\ {{\rm{s}}}^{-1}$, which exhibits chaotic behavior upon interaction with the galactic nucleus. The middle and lower panels show the orbits of stars with ${v}_{T}=10\,\mathrm{km}\ {{\rm{s}}}^{-1}$ and ${v}_{T}=-10\,\mathrm{km}\ {{\rm{s}}}^{-1}$, respectively, which exhibit well-behaved disk orbits. The star with no angular momentum spends little of its orbital period near the galactic plane, whereas the two stars that do not approach the nucleus remain within a few hundred parsecs of the plane.

Zoom In Zoom Out Reset image size

Without an ab initio model for the disk, it is difficult to predict how many stars are expected to be missing near zero angular momentum and what the exact profile of the dip should be. Therefore, we use the fraction of the orbital period of the test-particle stars that they spend near the mid-plane of the disk as a proxy for whether it has been scattered to a much higher scale height. We integrate many orbits with positions, radial, and vertical velocities drawn from the TGAS+RAVE sample and from an initial, uniform distribution in v_Y covering the low-velocity tail. We then re-weight the stars using the fraction of time they spend near the plane and construct the dip profile by dividing this weighted, final v_Y distribution by the initial, uniform distribution.

Figure 3 shows the computed dip profiles for velocities within $80\,\mathrm{km}\ {{\rm{s}}}^{-1}$ of v_T = 0. The top panel shows the results of seven simulations with reflex solar motion, ${v}_{\odot ,\mathrm{reflex}}$, of −220, −230, −235, −240, −245, −250, and $-260\,\mathrm{km}\ {{\rm{s}}}^{-1}$, assuming $| z| \lt 300$ pc as the criterion "close to the disk plane." It is clear that the shape and depth of the overlaid distributions are very similar and thus not dependent on the value of the solar motion. The dip profile is not entirely symmetric around zero ${v}_{Y}-{v}_{\odot ,\mathrm{reflex}}$, because the spatial distribution of the TGAS+RAVE sample is asymmetric around the Sun. The middle panel of Figure 3 displays the same as the upper panel, for ${v}_{\odot ,\mathrm{reflex}}=-240\,\mathrm{km}\ {{\rm{s}}}^{-1}$, but assuming $| z| \lt 1\,\mathrm{kpc}$ to weight the orbits. The profile of the dip is notably shallower, with a depth of ∼0.85, compared to ∼0.7 for the models in the top panel, but the shape of the dip profile is similar. The bottom panel of Figure 3 shows same as the middle panel, but assuming $| z| \lt 50$ pc to weight the orbits. The depth is consistent with the models in the top panel.

Zoom In Zoom Out Reset image size

We create a function $D({v}_{Y}-{v}_{\odot ,\mathrm{reflex}})$ of the dip by smoothly interpolating the model with ${v}_{\odot ,\mathrm{reflex}}=-240\,\mathrm{km}\ {{\rm{s}}}^{-1}$ and $| z| \lt 300$ pc displayed in the top panel of Figure 3. Because the dip profile computed above is an approximation, we give the model more freedom by allowing the dip's amplitude to vary. We model the full distribution $f({v}_{Y})$ of v_Y as an exponential multiplied by the dip:

$$\begin{eqnarray}&&f({v}_{Y})\\ &&\quad \propto \exp (m\,{v}_{Y}+b)\times \left(1-\alpha \left[\displaystyle \frac{1-D({v}_{Y}-{v}_{\odot ,\mathrm{reflex}})}{1-D(0)}\right]\right).\end{eqnarray} \tag{ 1 }$$

Written in this manner, the depth is parameterized by a parameter α, normalized such that $\alpha =0$ corresponds to no dip and $\alpha =1$ corresponds to a dip reaching all the way to zero, i.e., a complete absence of stars.

We fit the distribution of v_Y of disk stars ($| {v}_{X}| \lt 150\,\mathrm{km}\ {{\rm{s}}}^{-1}$ and $| {v}_{Z}| \lt 150\,\mathrm{km}\ {{\rm{s}}}^{-1}$) in the range $-310\,\mathrm{km}\ {{\rm{s}}}^{-1}\leqslant {v}_{Y}\ \lt -150\,\mathrm{km}\ {{\rm{s}}}^{-1}$ using the model in Equation (1). The top panel of Figure 4 displays the best-fit $\alpha =0$ model overlaid on the zoomed-in tail of the histogram of the distribution of v_Y. This clearly demonstrates the dearth of stars between approximately ${v}_{Y}=-210\,\mathrm{km}\ {{\rm{s}}}^{-1}$ and $-270\,\mathrm{km}\ {{\rm{s}}}^{-1}$ when compared with the exponential model.

Zoom In Zoom Out Reset image size

The best-fitting dip model is shown in the bottom panel of Figure 4. The data prefer a dip: the log likelihood of the best-fit model is $\mathrm{ln}{ \mathcal L }=282.0$, compared with $\mathrm{ln}{ \mathcal L }=276.7$ for the best-fit exponential model with two fewer parameters ($\alpha =0$). For comparison, we also fit a model where we allow the width to be controlled by a new parameter w that stretches the profile in Figure 3 along the x axis. For this model we find $\mathrm{ln}{ \mathcal L }=282.7$ and w = 1.8. Although the likelihood is marginally higher for this model, the significance of the detection is lower owing to the extra free parameter.

Using the Akaike information criterion (AIC; Akaike 1974), the difference in ln likelihood corresponds to a $2.1\sigma$ detection for the w = 1 model and a $1.8\sigma$ detection for the w free case. When combined with the $1.7\sigma$ detection from the v_X to v_Y analysis using a Bonferroni correction to account for the multiple tests (e.g., Dunn 1959) this gives a $2.7\sigma$ detection for w = 1 and a $2.5\sigma$ detection when w is left free. We explore the allowed range of models for w = 1 using a Markov Chain Monte Carlo (MCMC) analysis with emcee (Foreman-Mackey et al. 2013). We find that $\alpha =0.52\pm 0.1$ and ${v}_{\odot ,\mathrm{reflex}}=-239\pm 9\,\mathrm{km}\ {{\rm{s}}}^{-1}$ at $1\sigma$.

Even though we do not have an ab initio prediction of the dip profile, we find that the best-fitting α is consistent within its uncertainty to that obtained using the models (which have $\alpha \approx 0.3$). Thus, systematics due to our approximate model for the dip in a simple axisymmetric Milky Way potential are likely quite small.

We also repeat the analysis on a larger sample of stars, including radial velocities from both RAVE and the Large sky Area Multi-Object fiber Spectroscopic Telescope (LAMOST; e.g., Zhao et al. 2012) DR2. There is a small overlap of these catalogs and for sources in both catalogs we use the RAVE radial velocity measurements that are typically more accurate. We use the distance estimates from Astraatmadja & Bailer-Jones (2016) using the "Milky Way" prior because the LAMOST stars do not come with spectrophotometric distance estimates. With this sample, we obtain a measurement of ${v}_{\odot ,\mathrm{reflex}}=-236\pm 9\,\mathrm{km}\ {{\rm{s}}}^{-1}$ at $1\sigma$, with a detection significance of $2.5\sigma$ ($3\sigma$ when combined with the v_X to v_Y analysis). Note that we do not present this as our main result because the distances derived in Astraatmadja & Bailer-Jones (2016) are heavily dependent on the choice of prior, and thus at this stage we expect the RAVE spectrophotometric distances to be more reliable. However, we find it extremely encouraging that the feature is visible for two independent distance estimates and that the values measured are similar.

We have shown that a dip in the low-velocity tail of stars in the solar neighborhood is present at $2.7\sigma$ significance. This feature can be plausibly explained by the absence of stars on orbits with approximately zero angular momentum, due to such stars being scattered onto halo orbits after interaction with the Galactic nucleus as predicted by Carlberg & Innanen (1987). Using orbit integration in a realistic model of the Galactic potential, we have demonstrated that this dip should be centered on minus the solar rotational velocity ${v}_{\odot ,\mathrm{reflex}}=-{v}_{\odot }$. In this way, we were able to measure ${v}_{\odot }=239\pm 9\,\mathrm{km}\ {{\rm{s}}}^{-1}$. Our modeling does not include the Galactic bar or spiral arms, which affect the orbits of stars. However, because the dip is due to zero-angular-momentum stars moving on orbits that are highly eccentric that are being lost from the solar neighborhood over many dynamical times, the gravitational perturbations and current position of the relatively diffuse stellar bar and of the spiral structure should not greatly affect the profile of the dip and we expect their effect to be minor on our derived value of ${v}_{\odot ,\mathrm{reflex}}$.

Observations of the Galactic center have precisely measured the proper motion μ_SgrA^∗ of Sgr A* (Reid & Brunthaler 2004). Assuming that Sgr A* is at rest with respect to the dynamical center of the Milky Way, this apparent proper motion is due to the reflex motion of the Sun. The ratio of the linear and angular measurements of the reflex motion of the Sun ${v}_{\odot ,\mathrm{reflex}}/$ μ_SgrA^∗then becomes a measurement of the distance R₀ to the Galactic center. For μ_SgrA^∗ = 30.24 ± 0.11 km s⁻¹ kpc⁻¹, our measurement of ${v}_{\odot ,\mathrm{reflex}}$ implies ${R}_{0}=7.9\pm 0.3\,\mathrm{kpc}$. This is an essentially dynamical measurement of R₀. If the lack of zero-angular-momentum stars is confirmed in future Gaia data releases and systematics can be controlled, this feature can in principle saturate the precision implied by the uncertainty in μ_SgrA^∗, which is $\approx 30\,\mathrm{pc}$.

Future Gaia data releases will have many orders of magnitude more stars, with substantially higher accuracy on the measured parameters. Additionally, future Gaia data releases will contain line of sight velocities for stars in the solar neighborhood, removing the asymmetry in the sample introduce by the necessity to cross match with RAVE. Thus, with future Gaia data releases we hope that this technique may provide a very precise measurement of the solar reflex motion with respect to the Galactic center. We can make a rough prediction for the strength of the feature expected in future Gaia data releases by using the Gaia Object Generator (GOG; Luri & Palmer et al. 2014) based on the Gaia universe Model Snapshot (GUMS; Robin et al. 2012). GOG predicts that approximately 17.5 million stars will be observed with radial velocities within 700 pc, and thus, by weighting that number by the fraction of stars in the tail compared to the full sample matched with RAVE, e.g., 374/216,201 = 0.0017, we can predict that Gaia will observe roughly 30,000 stars in the tail within 700 pc. This should allow us to confirm the existence of this feature, and if present allow us to increase the precision of the measurement of ${v}_{\odot ,\mathrm{reflex}}$ to within approximately $1\,\mathrm{km}\ {{\rm{s}}}^{-1}$, assuming the errors decrease on the order of $\sqrt{N}$. In turn, the error on R₀ would become dominated by the error in the proper motion measurement. At this level, the systematics will then become extremely important, and in future analyses we must take into consideration other factors such as the bar, spiral structure, and local giant molecular clouds, although the tangential deflections are expected to be small.

Further work is also needed to fully determine the origin of this feature, and to make more detailed predictions on the shape and depth that we expect to observe. The exact dimensions of the feature are likely slightly dependent on the potential of the inner Galaxy, and thus matching the observed shape in future Gaia data releases to dip functions constructed from various models for the potential may give us valuable insight into the potential of the inner Galaxy.

J.A.S.H. is supported by a Dunlap Fellowship at the Dunlap Institute for Astronomy & Astrophysics, funded through an endowment established by the Dunlap family and the University of Toronto. J.B. and R.G.C. received partial support from the Natural Sciences and Engineering Research Council of Canada. J.B. also received partial support from an Alfred P. Sloan Fellowship. The MCMC analyses in this work were run using emcee (Foreman-Mackey et al. 2013).

This work has made use of data from the European Space Agency (ESA) mission Gaia (http://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC; http://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular, the institutions participating in the Gaia Multilateral Agreement. Funding for RAVE (http://www.rave-survey.org) has been provided by institutions of the RAVE participants and by their national funding agencies.