A Map of the Local Velocity Substructure in the Milky Way Disk

Recent results have shown that the stellar disk of the Milky Way exhibits coherent substructure in both spatial density and velocity. Using data from the Sloan Digital Sky Survey (SDSS; York et al. 2000) and the companion Sloan Extension for Galactic Understanding and Exploration (SEGUE; Yanny et al. 2009), both Widrow et al. (2012) and Yanny & Gardner (2013) found that at the position of the Sun there are more stars 1 kpc above the Galactic plane than 1 kpc below it, but fewer stars 400 pc above the Galactic plane than 400 pc below it. This points to a "bending" mode of the disk ("breathing" modes are symmetric in density about the Galactic plane). SDSS/SEGUE also uncovered vertical velocity substructure (Widrow et al. 2012), and oscillations of the stellar phase-space density (Gómez et al. 2012) that were thought to be induced by minor merger events in the Milky Way. Data from the RAdial Velocity Experiment (RAVE; Steinmetz et al. 2006) allowed Williams et al. (2013) to discover that the disk exhibits wave-like substructure in the V_Z, V_R, and V_θ directions, with amplitudes as high as 17 km s⁻¹. These disk substructures are also present in data from the Chinese Large Area Multi-Object Spectroscopic Telescope (LAMOST, also named the Guo Shou Jing Telescope after a Chinese astronomer and engineer who is famous for his instrumentation; Su & Cui 2004; Cui et al. 2012). Carlin et al. (2013) used LAMOST radial velocities to show that in the vicinity of the Sun the stars above the Galactic plane exhibited net flow outward from the Galactic center and downward toward the Galactic plane. Stars below the disk showed bulk flow in the opposite direction. Although all of the surveys showed similar types of structure, they used different tracers and probed different portions of the Galactic disk. It has been difficult to combine the results of each publication into a global view of disk substructure.

Xu et al. (2015) found that the midplane of the disk oscillates vertically, as one looks toward the Galactic anticenter from the position of the Sun, indicative of a bending mode of the disk. This was later supported by Morganson et al. (2016), using data from the Panoramic Survey Telescope and Rapid Response System 1 Survey (PS1; Kaiser et al. 2010) $3\pi$ data set. These ripples in the disk roughly follow the pattern of our Galaxy's spiral arms, as traced by H ii regions, giant molecular clouds, and methanol masers (Hou & Han 2014). The Monoceros Ring (Ibata et al. 2003; Yanny et al. 2003) at 20 kpc from the Galactic center, and the more distant TriAnd Cloud (Rocha-Pinto et al. 2004) might also be explained as ripples in a stellar disk that extends at least 25 kpc from the Galactic center. This latter suggestion is supported by the observation that the TriAnd and TriAnd2 substructures have metallicities indicative of association with a disk population rather than a dwarf satellite (Price-Whelan et al. 2015).

Although there might be other mechanisms for creating density and velocity substructure in the disk (Debattista 2014; Faure et al. 2014; Monari et al. 2015, 2016), the observed oscillations are roughly consistent with simulations of the disk's response to a large ($\gt {10}^{9}\,{M}_{\odot }$) satellite or dark subhalo passing near, or through, the disk (Kazantzidis et al. 2008; Younger et al. 2008; Minchev et al. 2009; Purcell et al. 2011; Gómez et al. 2013; Widrow et al. 2014; D'Onghia et al. 2016; Laporte et al. 2016). Corrugations have also been observed in the Milky Way's H i gas and young stellar populations in the Milky Way (Spicker & Feitzinger 1986) and in external galaxies (Matthews & Uson 2008a). It has been proposed that these oscillations were induced by dwarf galaxy satellites (Chakrabarti et al. 2011). There is the potential that disk substructure could be used to probe the properties of Milky Way satellites, even if the satellites are composed of dark matter and are not directly visible (Feldmann & Spolyar 2015), a possibility called galactoseismology by Widrow et al. (2012). Chakrabarti et al. (2015) used ripples observed in the Milky Way gas to infer the presence of a previously undetected dwarf satellite of the Milky Way, though the existence of this dwarf galaxy is disputed (Pietrukowicz et al. 2015).

Gómez et al. (2016, 2017) examined fully cosmological simulations of Milky Way-like galaxies, and concluded that satellite-induced vertical oscillations of the disk are common. In addition, these oscillations can explain density substructures like the ones observed in the Monoceros Ring. More impressively, they show that as a satellite falls into the Milky Way, it can perturb the density substructure of the halo, which can in turn produce perturbations in the disk (see also Weinberg 1998; Vesperini & Weinberg 2000). This means that previous simulations of the effects of the Sagittarius dwarf galaxy infall on the disk were missing an important contribution—the effect of a redistribution of halo mass on the disk. This might explain why previous simulations had difficulty producing the magnitude of the disk oscillations observed in the Milky Way.

The similarity between the observed ripple pattern, the Milky Way's spiral arms, and the spiral structure produced in the simulations, leads one to wonder whether the spiral structure in galaxies is induced and sustained by orbiting dwarf galaxies (and possibly dark subhalos). Although it has been known for some time that galaxy interactions can produce spiral arms (Toomre & Toomre 1972), the source of spiral structure, like that of the Milky Way, has traditionally been explained under the assumption of quasi-stationary spiral structure (QSSS; Shu 2016). The possibility that dwarf galaxies could be exciting the disk substructure observed in the Milky Way, and the knowledge that the dwarf galaxy interactions also produce spiral structure, could mean that the QSSS hypothesis, and thus spiral density waves, may not be required to sustain spiral structure in Milky Way-like galaxies (Newberg & Xu 2017). In support of this, note that de la Vega et al. (2015) showed that disk substructure generated from dwarf galaxy interactions is the result of the phase wrapping of epicyclic perturbations, without the self-gravity required by wave propagation in breathing or bending modes. While some effort has gone into understanding the effects of the Galactic bar and existing spiral arms on disk substructure (e.g.,—Debattista 2014; Faure et al. 2014; Monari et al. 2016; Vera-Ciro & D'Onghia 2016), the connection between satellites and the origin of normal spiral structure has not been explored.

In this paper, we reproduce the velocity substructure within ∼1.5 kpc of the Sun, following a technique and data set very similar to that of Carlin et al. (2013). However, we have more data, greatly improved distance estimates, improved proper motions, and we propagate errors in our measurements to produce a map of the velocity substructure in a portion of the local disk. Our maps confirm the compression of the disk reported outside of the solar circle by Williams et al. (2013), though we observe the compression in the third quadrant and not in the second. The maps also confirm the nearer bulk motion of the disk found by Carlin et al. (2013); the more distant parts of the Carlin et al. data had larger velocity errors and were not used in our study. We observe an increase in disk rotation speed at greater radii outside the solar circle, which is comparable to the trend found in López-Corredoira (2014), but we map the trend in three dimensions.

We present a three-dimensional map, with errors, of the disk velocity substructure in a region approximately $1.5\times 1.5\,\times 1.5\,\mathrm{kpc}$ in size, centered on the Galactic plane and extending out from the Sun toward the anticenter. While the disk appears to have zero Galactocentric radial velocity in the third quadrant, the stars in the second quadrant, just half a kiloparsec away, exhibit a bulk motion toward the Galactic center of about 15 km s⁻¹ and also a compression in the Z-direction toward a point about 0.5 kpc above the Galactic plane. We observe that the magnitude of the circulation velocity decreases as a function of distance from the Galactic plane, as has been known for many years. The circulation velocity of stars within 0.2 kpc of the Galactic plane is 205 km s⁻¹ near the Sun's position, and increases to 225 km s⁻¹ 1.5 kpc farther from the Galactic center. The velocity substructure map we present is the best to date because it shows the statistical significance of the measured bulk velocity with a granularity of 0.2 kpc in a local small region.

LAMOST (Cui et al. 2012) is one of the National Major Scientific Projects undertaken by the Chinese Academy of Sciences. It is a 4 m class, downward-facing Schmidt reflector with 4000 optical fibers in the focal plane, that was designed for spectroscopic surveys. The spectra have a resolution $R\,=\lambda /{\rm{\Delta }}\lambda =1800$, which is similar to the resolution of SDSS/SEGUE. The telescope is located at the Xinglong Observing Station ($40^\circ 23^{\prime} 39^{\prime\prime} \,{\rm{N}}$ $117^\circ 34^{\prime} 30^{\prime\prime}$ E), one hundred miles northeast of Beijing. Due to the fixed second mirror (called "mirror B" because it is the second reflecting surface but corresponds to the primary mirror of a Schmidt telescope), the optics favor targets closer to the zenith. Local weather limits the observing season to the months of September to May. As a result of the fixed mirror B and restricted observing season, survey observations are restricted to decl. $-10^\circ \lt \delta \lt 60^\circ$, and R.A. near the Galactic anticenter are preferred.

The optical properties of the telescope system and the site location limit the magnitudes of spectroscopic targets to $V\sim 19.5$ (Luo et al. 2012), though most of the observed objects are close to $V\sim 16$. The telescope has been used to carry out a spectroscopic survey (Zhao et al. 2012) that primarily yields Milky Way stars, and most of those stars are in the Galactic disk. For example, in the 2014/2015 observing season, the survey netted 1,558,924 spectra, of which 1.4M are classified as Galactic stars, 23,000 are classified as galaxies, and 7,000 are classified as QSOs (He et al. 2016). Details of the LAMOST spectral survey can be found in Zhao et al. (2012). The science plan for the LEGUE portion of the survey can be found in Deng et al. (2012). The pilot survey began in 2011 October, and the general survey started in 2012 September and is still in progress. The survey has publicly released the pilot survey (Luo et al. 2012), Data Release 1 (DR1; Luo et al. 2015), and Data Release 2 (DR2; Luo et al. 2016). LAMOST DR2 includes measured parameters for 3.9M spectra obtained through spring 2014.

The data used in this work include stellar parameters measured by the LAMOST pipeline for all stars in LAMOST Data Release 3 (DR3, which includes data from the first 3 years of survey operations plus the Pilot Survey), and data from the first 2 quarters of the fourth year of survey operations (data observed through 2016 February). This encompasses a total of ∼3.9 million stellar spectra of ∼3.1 million unique stars. Extinction corrections throughout are based on the maps of Schlegel et al. (1998). In addition to the measured stellar parameters, distances derived via the method of Carlin et al. (2015) are included in our analysis. For stars that have been observed multiple times, the stellar parameters and distances represent weighted mean values of all measurements, weighted by spectroscopic signal-to-noise ratio (S/N).

Positions and Proper Motions "Extra Large" (PPMXL; Roeser et al. 2010) is a full-sky catalog containing the positions and proper motions for about 900 million objects. However, it is known (Wu et al. 2011; Grabowski et al. 2015; Vickers et al. 2016) that the proper motions calculated by PPMXL contain significant systematic error. We again confirm this systematic error by examining the PPMXL proper motion measurements of quasars and galaxies identified by LAMOST. Since these objects are far too distant to contain any true proper motion, the proper motion measurements result entirely from a combination of statistical and systematic error.

Any systematic error in proper motion measurements within a particular region of the sky causes a systematic shift in the perpendicular component of velocity that increases linearly with distance. This can become very significant, even for observations within the local disk region (1 kpc × 1 mas ${\mathrm{yr}}^{-1}=4.74$ km s⁻¹). Therefore, it is extremely important to remove as much of this systematic error as possible before calculating the three-dimensional stellar velocities.

Published proper motion corrections are based on a variety of fitting methods of extragalactic objects with varying ranges of precision, from constant shifts (Wu et al. 2011) to correction tables (Grabowski et al. 2015). Other corrections implement analytic fitting methods (Carlin et al. 2013; Vickers et al. 2016), which are effective in smoothing out large-scale systematics, but do not necessarily work well in all areas of the sky, especially in regions with low Galactic latitude. In this work, we analyze the correction from Vickers et al. (2016), which is magnitude-dependent and utilizes spherical harmonics fits, since it is the previous best global correction for the PPMXL database.

The template-matching algorithm in the LAMOST pipeline categorizes objects based on their best template match as "STAR," "GALAXY," "QSO," or "UNKNOWN," where the UNKNOWN category is a catch-all for spectra that do not match any of the templates well. To derive proper motion corrections, we select all extragalactic objects classified as GALAXY or QSO in the data set up through DR4Q2. These should constitute a fixed reference frame for proper motion zero-points. As with the stellar catalog, we match the extragalactic objects to PPMXL with a $1\buildrel{\prime\prime}\over{.} 5$ matching radius. We remove all duplicates and select objects with a proper motion less than $| \mu {| }_{\max }=30$ mas yr⁻¹ (Figure 1; the determination of this value is explained further in Section 2.3) in order to remove outliers and possibly nearby stars with large proper motion that were misidentified. After all selections, we obtained a sample of 100,367 LAMOST spectra of quasars and galaxies, along with their corresponding proper motions from the PPMXL catalog.

Zoom In Zoom Out Reset image size

Since the true proper motion of these objects must be zero, we expect to see an approximately uniform Gaussian random distribution around zero along each direction to account for statistical error. From a quick glance at a plot of proper motions over the sky (Figure 2), it is very clear that there is systematic clumping of proper motions based on their celestial coordinates, in addition to an overall negative net motion in each component. In comparison, the Vickers et al. (2016) correction significantly decreases the error in most parts of the sky. However, we show that there are still many parts of the sky where this error is actually increased (Figure 3). This occurs very often in the surrounding area of regions with large systematic error, indicating that the systematics of these objects deviate over a scale too small to be handled well by an analytic function. We show in the next section that systematic error of the PPMXL proper motions is greatly reduced if only information in a small region of the sky around each object is used to determine the systematic offsets (Figure 4).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Due to the apparently random small angle variations in the systematic error (see Section 2.6), fitting the proper motions of extragalactic objects to any kind of analytic function would need to be of extremely high order to ensure an improvement in all regions of the sky. It is also unclear, using this method, which regions of the sky can be confidently used without a function to give an estimate of remaining systematic error. Giving a quantifiable estimation for this remaining systematic error is absolutely necessary for measuring the statistical significance of stellar velocities. An alternative to this method, which we choose to use in this work, is to construct a table of systematic offsets over many bins across the sky, independently calculating the correction needed at each position to produce a result somewhat similar to Grabowski et al. (2015).

The concept of our correction is to construct a large table of corrections over R.A. α and decl. δ, where at each correction point (α,δ), we return correction values based on the systematics of extragalactic objects near this position, but not necessarily only those within the region constrained by the bin. From a given correction point, we search within a surrounding search radius for nearby quasars and galaxies from our LAMOST data. The opposite of the mean of the proper motions of these objects is calculated as the systematic correction value for this point, in each component. This is the value that we add to any object in the 025 by 025 bin centered at this coordinate (see Section 2.5).

Using this method, there are a few subjective parameters that we define here: since our method is fairly vulnerable to outliers and the LAMOST data set contains objects with proper motion upwards of several hundred mas yr⁻¹, we initially eliminate all objects with a proper motion above an upper bound of $| \mu {| }_{\max }=30$ mas yr⁻¹ during our data selection (Section 2.2). We also define a search radius ${r}_{\mathrm{search}}=2\buildrel{\circ}\over{.} 5$ (Figure 5), which is the maximum angular separation allowed in order to be considered a nearby object and used in the correction. Objects further apart than this are not used in the correction. In order to eliminate any rare instances of two or three objects with large and similar random errors mimicking a precise and accurate correction, we set a minimum number of objects, ${n}_{\min }=4$, that must be found in order to return a correction. Lastly, in order to cut down on excessive computation time, we stop searching after the nearest ${n}_{\max }=100$ objects are found.

Zoom In Zoom Out Reset image size

We tested many different values for $| \mu {| }_{\max }$ and ${r}_{\mathrm{search}}$ before arriving at these values, in order to optimize the correction. They were tested by attempting to correct each object using only its nearby objects, and not itself, to ensure self-correction did not contribute to the success of the entire correction. For testing purposes, the correction points were defined to be the positions of each extragalactic object with which we tested the algorithm; for the correction table, correction points were selected to be at the centers of each $0\buildrel{\circ}\over{.} 25\times 0\buildrel{\circ}\over{.} 25$ bin. Using the resulting correction, we optimized two main criteria. First, a lower value for the average magnitude of proper motion would indicate high precision in the values output by the correction; however, most of this is due to random error, which is impossible to correct. Second and more importantly, a normal distribution of positive and negative values centered around zero in each proper motion component, independent of position, would indicate a lack of systematic error.

Both criteria were tested partially by eye through attempting to eliminate large values and obvious systematics in Figure 2 versus Figure 3 versus Figure 4. Additionally, we measured the second criterion via the trend of the average magnitude of proper motion at various sampling rates (Figure 6). This test gradually increases the subsampling interval, starting at 0°, which means the magnitude of the proper motion of each individual object is calculated and averaged. The interval was tested up to 45°, at which each $45^\circ \times 45^\circ$ bin first averages the proper motions of all objects within it and should be close to zero due to cancellation of many positive and negative values, resulting in a near zero average of the magnitude of each bin if there is no net systematic shift. Since we are comparing the same set of objects in exactly the same bins in Figure 6, a trend that approaches mean magnitude of proper motion equaling zero more quickly, indicates a smaller amount of systematic error. In other words, the first value at a sampling rate of zero is predominantly indicative of the magnitude of random error in proper motion. At larger sampling rates, the value is dominated by systematic error, if present.

Zoom In Zoom Out Reset image size

Using our correction method, we are also able to produce an estimation of the accuracy of the correction. Since we subtract the mean proper motion of a set of data to create the correction, the uncertainty is simply the standard error of the mean of this data set. This uncertainty in the correction is how we calculate our estimation of remaining systematic error, since any statistical error in the correction will propagate through every star in its particular region of the sky. Therefore, systematic error (in the mean) of each proper motion component is calculated as

$$\begin{eqnarray}&&{\mathrm{SE}}_{\mu }=\displaystyle \frac{{\sigma }_{\mu }}{\sqrt{n}},\end{eqnarray} \tag{ 1 }$$

where n is the number of nearby quasars and ${\sigma }_{\mu }$ is the standard deviation of their PPMXL proper motion measurements in the respective component. This value is included along with every correction value in order to determine our confidence in proper motion measurements, and therefore estimate the remaining systematic error present, in each region of the sky.

We constructed a table of bins over very short intervals, spanning the entire region of the sky covered by LAMOST, in which the center position of each bin was assigned a correction value via the method described in Section 2.3. This resulted in a very long proper motion correction table that is a high-resolution representation of the sky between a decl. of $-15^\circ$ and 75^◦, where each row corresponds to a bin of size $0\buildrel{\circ}\over{.} 25$ by $0\buildrel{\circ}\over{.} 25$. Entries in our table with a value of 100.0 in the systematic error columns represent a location lacking sufficient data ($n\lt 4$) to justify any attempt at correction, and therefore contain a value of zero in their proper motion shift columns. Note that there is only sufficient data in the range of $-10\buildrel{\circ}\over{.} 5\lt \delta \lt 61\buildrel{\circ}\over{.} 75$; however, we defined our range of decl. using round numbers for more convenient implementation of the correction. We present our correction table in Table 1, in order to show the format of each column, and also in Figure 7, which is an image representation of the entire correction and where it is reliable.

Zoom In Zoom Out Reset image size

The columns of our correction table are as follows: "ra" and "dec" give the center values of α and δ of each bin, in degrees. For a star to be assigned to a bin, it must be within a range of ±0.125 from the given center values in both α and δ. "pmra_shift" and "pmde_shift" give the value that needs to be added to the uncorrected proper motion from PPMXL in each component, in mas yr⁻¹. "pmra_se" and "pmde_se" give the expected remaining systematic error for this location in the sky for their respective components, as calculated in Section 2.4. Note that proper motion in R.A. is ${\mu }_{\alpha }\cos \delta$, while proper motion in decl. is simply ${\mu }_{\delta }$.

In the vast majority of the sky, the systematic error in proper motion was significantly decreased by our correction, even compared to the current standard correction to the PPMXL database (Vickers et al. 2016). This can be seen in both Figures 4 and 6. The difficulty in creating a large area correction to PPMXL seems to be due to the nature of how the systematic error arises. The error is clumped randomly throughout the sky (Figure 2), on the order of only a few degrees, which may suggest that the systematic error is dependent on the photographic plate on which it was measured. This makes fitting the systematic error to any analytic function extremely difficult, and justifies the use of our correction table (Table 1, Figure 7).

The area of greatest concern with this, as well as many other proper corrections, is unfortunately near the Galactic plane. It appears that our correction still does slightly better at minimizing the systematic error in this region, but it is difficult to compare due to the small amounts of data. However, our correction may still be more useful in this region because it includes an estimation of remaining systematic error, which is inversely correlated with the number of nearby extragalactic objects from LAMOST, and directly with the dispersion of nearby extragalactic proper motion measurements from PPMXL. This allows for one to easily identify which parts of the data are meaningful, which is extremely important for discovering actual velocity substructure of disk stars and suppressing deceiving trends caused by systematic error.

In this paper, we adopt the same Galactic Cartesian $(X,Y,Z)$ and right-handed Galactic cylindrical $(R,\theta ,Z)$ coordinate systems used in Carlin et al. (2013) to represent positions. Note that both coordinate systems are Galactocentric, and θ increases in the same direction as Galactic longitude, with the Sun centered at $\theta =0$. We use $({V}_{X},{V}_{Y},{V}_{Z})$ and $({V}_{R},{V}_{\theta },{V}_{Z})$ to represent the Galactic Standard of Rest velocities in these coordinate systems. Note that the disk rotates opposite to the θ direction, with a local velocity of approximately ${V}_{\theta }=-210$ km s⁻¹. In these coordinates, we adopt a solar position and velocity of ${(X,Y,Z)}_{\odot }=(-8,0,0)$ kpc and ${({V}_{X},{V}_{Y},{V}_{Z})}_{\odot }\,=(10.1,224.0,6.7)$ km s⁻¹, respectively (Hogg et al. 2005).

For this work we use cataloged radial velocities and stellar parameters from the LAMOST survey, including all objects available internally to collaborators through quarter 2 of the fourth observing season (DR4Q2; including the Pilot Survey, Data Release 1: DR1, DR2, and DR3). Our initial stellar sample includes all objects classified as "STAR" by the LAMOST stellar parameters pipeline. Many stars in the LAMOST catalogs have been observed multiple times over the ∼5 years the telescope has operated; for these objects, we combine the velocities, stellar parameters, and distances via a weighted average, with weights determined by the spectroscopic S/N. The final catalog after removal (and combination of) duplicate measurements includes ∼3.1 million unique stars, of which roughly 3.04 million have r-band S/N larger than 10 (median ${\rm{S}}/{{\rm{N}}}_{r}\sim 54$). The median uncertainties in measured/derived spectroscopic quantities for the ${\rm{S}}/{{\rm{N}}}_{r}\gt 10$ sample are ${\sigma }_{\mathrm{RV}}=13.5$ km s⁻¹, ${\sigma }_{{T}_{\mathrm{eff}}}=114$ K, ${\sigma }_{[\mathrm{Fe}/{\rm{H}}]}=0.15$, ${\sigma }_{\mathrm{log}g}\,=0.40$, and ${\sigma }_{\mathrm{dist}}=0.25$, where ${\sigma }_{\mathrm{dist}}$ is the median of the fractional distance errors, $\tfrac{| \delta d| }{d}$. This catalog was matched positionally (with a maximum matching radius of $1\buildrel{\prime\prime}\over{.} 5$) to the PPMXL (Roeser et al. 2010) catalog of proper motions.

We calculated the proper motion of each star in each component (${\mu }_{\alpha }\cos \delta$ and ${\mu }_{\delta }$) by adding its uncorrected proper motion from PPMXL to the correction from the row in Table 1 corresponding to the star's position. We also assigned a value of systematic error ${\mathrm{SE}}_{\mu }$ to each component of proper motion, as determined by the corresponding row from Table 1 for each star, in addition to the given statistical random error from the PPMXL catalog. We derive distances to stars, and their uncertainties, based on the method of Carlin et al. (2015), which uses a Bayesian method to estimate the absolute magnitude of each star via a comparison of its measured stellar parameters (${T}_{\mathrm{eff}}$, [Fe/H], and $\mathrm{log}g$) to a grid of isochrones. We use the distance d from an error-weighted mean of the distance calculations made from the the estimated absolute magnitude in the 2MASS K-band and the SDSS r-band, whose extinction-corrected (Schlegel et al. 1998) apparent magnitudes are obtained from LAMOST. We use radial velocity v_r, and its statistical random error, as measured by LAMOST.

We transform the equatorial positions $(d,\alpha ,\delta )$ and velocities $({v}_{r},d{\mu }_{\alpha }\cos \delta ,d{\mu }_{\delta })$ of every star into Galactic cylindric coordinates $(R,\theta ,Z)$ and $({V}_{R},{V}_{\theta },{V}_{Z})$. We also transform and propagate all statistical errors σ and systematic errors SE into the same Galactic cylindric coordinates.

We reduced our data set down to 335,599 stars via many selection cuts (a full list of our cuts is shown in Table 2). In order to conserve direct comparability with the results from Carlin et al. (2013), we use only F-type stars by selecting stars that have been classified by LAMOST with a subclass beginning with "F." Our R cut confines the data to a 2 kpc ring beyond the solar radius, while the θ and Z cuts are intended to be as inclusive as possible; any data outside of this region were far too sparse to be used. Our velocity cuts are intended to throw away outliers that are most likely not disk stars, but rather nearby high velocity halo stars. We cut random error and S/N in order to remove potential outliers, while systematic error cuts remove some especially bad regions of the sky where the systematic error may be very large, due to little or no extragalactic data. All position and velocity cuts were performed after the coordinate transformation described in Section 3.3. Note that these cuts only define the volume of study, remove statistical outliers and error code values, and remove data with large systematic or random errors. A heat map of the Galactic coordinates of our stellar sample is shown in Figure 8. We also present an HR Diagram of our stellar sample in Figure 9.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2.  Stellar Sample Selection Criteria

Classified as F-type star by LAMOST

${\rm{S}}/{\rm{N}}\gt 5$ in g, r, and i

$8\lt R\lt 10$ kpc

$-2\lt Z\lt 2$ kpc

$9^\circ \gt \theta \gt -12^\circ$

$-150\lt {V}_{R}\lt 150$ km s−1

$-400\lt {V}_{\theta }\lt -100$ km s−1

$-150\lt {V}_{Z}\lt 150$ km s−1

${\mathrm{SE}}_{{V}_{R}}\lt 8$ km s−1 (Systematic Error)

${\mathrm{SE}}_{{V}_{Z}}\lt 10$ km s−1

${\mathrm{SE}}_{{V}_{\theta }}\lt 10$ km s−1

${\sigma }_{{V}_{R,Z,\theta }}\lt 40$ km s−1 (Random Error)

${\sigma }_{R}\lt 0.2$ kpc

${\sigma }_{Z}\lt 0.2$ kpc

${\sigma }_{\theta }\lt 1\buildrel{\circ}\over{.} 3$

Note. Selection criteria applied to LAMOST data. These criteria limit the dataset to one type of star, in a region of the galaxy in which there is sufficient data to study disk substructure. Other criteria eliminate data with large errors in position or velocity.

Download table as:  ASCIITypeset image

We binned our stellar data into approximately $0.2\times 0.2\,\times 0.2\,\mathrm{kpc}$ large spatial bins in R,Z, and θ (we use $1\buildrel{\circ}\over{.} 3$ bins in θ, which have a length between 0.18 and 0.23 kpc within $8\lt R\lt 10$ kpc). These bin sizes were selected because the distance error for each star is about 20% (Carlin et al. 2015), so the positions of stars are typically smoothed out along our line of sight by 200 pc. Because the coordinate system in which we are observing the bulk motions is Galactocentric and the direction of smoothing is heliocentric, there are significant errors in our knowledge of the spatial coordinates of stars in all three axes ($R,\theta ,Z$). We did not observe variations in the mean velocities of stars at scales smaller than 0.2 kpc, but that does not mean that smaller-scale substructure does not exist; it would have been blurred out by errors in our knowledge of the location of each star.

In each bin, we calculated the number of stars N, the mean velocity in each component $(\langle {V}_{R}\rangle ,\langle {V}_{\theta }\rangle ,\langle {V}_{Z}\rangle )$, the mean systematic error value of each velocity component $(\langle {\mathrm{SE}}_{{V}_{R}}\rangle ,\langle {\mathrm{SE}}_{{V}_{\theta }}\rangle ,\langle {\mathrm{SE}}_{{V}_{Z}}\rangle )$, and the propagated error in each velocity component (${\sigma }_{{V}_{R}},{\sigma }_{{V}_{\theta }},{\sigma }_{{V}_{Z}}$). In our analysis, we measured the uncertainty of each velocity component in each bin by summing in quadrature $\langle {\mathrm{SE}}_{V}\rangle$ and ${\sigma }_{V}$.

We present these three-dimensional mean velocities, along with their uncertainties, in Figure 10 by averaging over the bins in one dimension per panel in order to get the best overall representation of the data from each perspective. Since the bins that we combined are separated by large distances ($\gt 3^\circ$) in the sky, the systematic errors should be uncorrelated (note that in Section 2.6 and from Figure 6 we deduced that any remaining systematic error is only correlated to a few degrees in the sky), so we did not average all systematic error values this time. Instead, we simply propagate the uncertainties of each bin ($\langle {\mathrm{SE}}_{V}\rangle$ and ${\sigma }_{V}$ summed in quadrature), thereby substantially reducing the uncertainty with the number of bins.

Zoom In Zoom Out Reset image size

In order to show exactly where our data lie, we created Figure 11 with the individual $0.2\times 0.2\times 0.2\,\mathrm{kpc}$ bins by stepping through $1\buildrel{\circ}\over{.} 3$ slices of θ in each panel. It is important to compare apparent trends from Figure 10 to the individual bins in this plot to be sure that the trends are not misleading due to regions of missing data. Along with Figure 11, we present the data file of individual bins containing each component of $\langle V\rangle$, ${\sigma }_{V}$, and $\langle {\mathrm{SE}}_{V}\rangle$.

Zoom In Zoom Out Reset image size

Although our data set only encompasses a small portion of the Milky Way disk, it currently contains the most precise bulk velocity calculations to date, as well as estimates of their precision. This type of data could be very useful for fitting to simulations of the Milky Way disk, and would be appropriate to use with a fitting algorithm that assigns higher weights to data with smaller uncertainties. We also hope to increase the scope of this data in the future by incorporating observations from additional facilities in order to cover a larger region of the sky.