The Stellar Velocity Distribution Function in the Milky Way Galaxy

Understanding galaxy formation is a major goal of current astrophysical research. Detailed studies of the stellar components of the Milky Way galaxy play a fundamental role in building this understanding because they provide direct and robust tests of and constraints on theories of galaxy formation. The latter are presently formulated within the framework of hierarchical formation in a ΛCDM universe (e.g., White & Frenk 1991; Mo et al. 1998), and the past few decades have seen a myriad of results affirming the role that mergers have had in the evolution of the Milky Way, no more vividly than the countless extant substructures that have been discovered and mapped throughout the Galactic stellar halo, whose origin appears to be dominated by the debris of these mergers.

Meanwhile, explorations of the evolution of the Galactic disk have a long and rich history that reaches back to the kinematical studies of Eddington (1915), Jeans (1915), Corlin (1920), Kapteyn (1922), Oort (1926), and Strömberg (1927), who explored the local stellar velocity distribution to gain information about large-scale Galactic structure. The field significantly advanced when Roman (1950) and Parenago (1950) discovered a strong correlation between spectral type and UV excess (i.e., metallicity) and the kinematical properties of stars near the Sun. One of the explanations of this phenomenon, which led to the seminal paper by Eggen et al. (1962), is that the oldest and least chemically enriched stars we observe today were born kinematically hot, and these turbulent birth velocities declined steadily as the gas in the Galaxy dissipated energy and settled into a dynamically cold disk—i.e., older populations were dynamically hotter ab initio (e.g., Brook et al. 2004; Minchev et al. 2013). However, it was later suggested that within some fraction of the thin disk, stars are born with random velocities that are small at birth but that increase with time (e.g., Spitzer & Schwarzschild 1951; Fuchs & Wielen 1987; Binney et al. 2000; Nordström et al. 2004; Aumer & Binney 2009).

The velocity distribution function (DF) of stars in the Galaxy can be a valuable aid in uncovering the relationships between kinematics, metallicity, and age for disk and halo stars, and lends insights into the dynamical history of stellar populations (Bienaymé 1999; de Simone et al. 2004; Famaey et al. 2005; Bovy & Hogg 2010). In the Galactic stellar halo, where the timescales for dynamical mixing are long, the history of minor mergers imprints long-lived, dynamically cold substructures that are quite legible fossils from which Galactic archeology can readily synthesize a sequence of events (e.g., Ibata et al. 1994; Odenkirchen et al. 2001; Majewski et al. 2003; Belokurov et al. 2006; Grillmair & Dionatos 2006).

However, a multitude of processes can affect the velocity DF in the Galactic disk, including the potential mixture of the aforementioned ab initio "born hot" and"born cold" components, and with the latter affected by numerous disk-heating mechanisms. The sum of these effects that increase velocity dispersions in disk stars complicates the simplest model of a spiral galaxy, where stars in the disk librate about circular orbits in the Galactic plane. From state-of-the-art magnetohydrodynamical cosmological simulations, Grand et al. (2016) found that in the secular evolution of the disk, bar instabilities are a dominant heating mechanism. The authors also reported that heating by spiral arms (e.g., de Simone et al. 2004; Minchev & Quillen 2006), radial migration (e.g., Sellwood & Binney 2002), and adiabatic heating from midplane density growth (e.g., Jenkins 1992) are all subdominant to bar heating. Another fundamental astrophysical process that scatters stars onto more eccentric orbits or onto orbits that are inclined to the disk's equatorial plane is the accretion of lower-mass systems (see, e.g., Kazantzidis et al. 2008; Martig et al. 2014, and references therein). This is particularly relevant given that recent studies suggest that the Milky Way likely underwent such a merger ∼8 Gyr ago (e.g., Nissen & Schuster 2010; Nidever et al. 2014; Linden et al. 2017; Belokurov et al. 2018; Kruijssen et al. 2019). Moreover, Laporte et al. (2018), using N-body simulations of the interaction of the Milky Way with a Sagittarius-like dSph galaxy, showed that the patterns in the large-scale velocity field reported from observations, like vertical waves (e.g., Widrow et al. 2012), can be described by tightly wound spirals and vertical corrugations excited by Sagittarius (Sgr) impacts.

Studies of the density distribution of stars above and below the Galactic plane have shown that the Galactic disk is in fact composed of two distinct components (Yoshii 1982; Gilmore & Reid 1983): a thin disk with a vertical scale height of 300 ± 50 pc and a thick component with a scale height of 900 ± 180 pc at the distance of the Sun from the Galactic Center (GC; see Bland-Hawthorn & Gerhard 2016 for a review of the spread in these quantities found by various authors using different methods). On the other hand, detailed stellar-abundance patterns for disk stars show a clear bimodal distribution in the [α/Fe] versus [Fe/H] plane. This bimodality is associated with—and firm evidence for—distinctly separated thin and thick disks (Edvardsson et al. 1993; Bensby et al. 2003; Fuhrmann 2011; Adibekyan et al. 2013). The high-[α/Fe] sequence, corresponding to the thick disk, exists over a large radial and vertical range of the Galactic disk (e.g., Nidever et al. 2014; Hayden et al. 2015; Duong et al. 2018).

The existence of two chemically distinguished disk components points to a different chemical evolution and, hence, distinct disk-formation mechanism and epoch (e.g.; Chiappini et al. 1997; Masseron & Gilmore 2015; Clarke et al. 2019) for the two disk components. Kinematically, the velocity dispersion for the chemically selected thick-disk component is, on average, larger than the dispersion reported for the low-[α/Fe] thin disk (Soubiran et al. 2008; Lee et al. 2011; Anguiano et al. 2018). However, there is substantial overlap in the thin- and thick-disk velocity DFs, which makes kinematics a less-reliable diagnostic for discriminating populations. From the standpoint of isolating and studying the chemodynamical properties of these populations free of cross-contamination, this is a bit unfortunate, because, despite many large, dedicated spectroscopic surveys of Galactic stars from which detailed chemical-abundance patterns can be measured for millions of stars (see below), there are now several orders of magnitude more stars with kinematical data available, thanks to ESA's Gaia mission.

Gaia is an all-sky astrometric satellite from which we can obtain accurate sky position, parallax, and proper motion, along with the estimated uncertainties and correlations of these properties, for ∼1.3 billion sources. Moreover, the existence of a Gaia subsample containing line-of-sight velocities provides unprecedented accuracy for individual space velocities for more than 7 million stellar objects (Gaia Collaboration et al. 2018). Unfortunately, there are no α-element abundance measurements for the vast majority of this Gaia subsample. This prevents an unbiased and comprehensive study of the velocity DF for the thin disk, thick disk, and halo treated separately.

Fortunately, the astronomical community has invested great effort into building and developing massive multiobject spectroscopic surveys, where the measurement of abundances beyond the simple overall metallicity is possible. Projects like SEGUE (Yanny et al. 2009), RAVE (Zwitter et al. 2008), LAMOST (Zhao et al. 2012), Gaia-ESO (Gilmore et al. 2012), GALAH (de Silva et al. 2015), and the Apache Point Observatory Galactic Evolution Experiment (APOGEE) (Majewski et al. 2017) are transforming our understanding of the Milky Way galaxy through their generation of vast chemical-abundance databases on Milky Way stars. The future is even more promising, with even larger spectroscopic surveys planned, such as Sloan Digital Sky Survey (SDSS) V (Kollmeier et al. 2017), WEAVE (Dalton et al. 2012), MOONS (Cirasuolo & MOONS Consortium 2016), DESI (DESI Collaboration et al. 2016), 4MOST (Feltzing et al. 2018), and MSE (Zhang et al. 2016), all aiming to provide individual abundances for millions more Galactic stars in both hemispheres, together with line-of-sight velocities, v_los, with a precision of a few hundred m s⁻¹. While the sum total of observed stars in these spectroscopic surveys will span only ∼1% of even the present Gaia subsample containing line-of-sight velocities, the chemical data in the previous surveys can be used to begin characterizing and interpreting the vastly larger number of sources having kinematical data in the Gaia database.

One of the main goals of the present study is to perform, for the first time, a detailed and unbiased study of the Galactic velocity DFs—derived from Gaia data—for individual, chemically separated stellar populations and to explore how these distributions change for different Galactocentric radii and distances from the Galactic midplane. For this study, we use the individual stellar abundances from the APOGEE survey, specifically [Mg/Fe] and [Fe/H], to relatively cleanly discriminate the thin and thick disks, associated with the low and high α sequences, respectively, in the [α/Fe]–[Fe/H] plane, as well as the halo stars, which predominantly inhabit other regions of the same plane. Using the kinematical properties of these chemically defined subsamples, we build a data-driven kinematical model, which we then apply to the full Gaia database to ascertain the contribution of the different Galactic structural components to the velocity-space DF as a function of Galactic cylindrical coordinates, R and z. We also create two-dimensional maps in the R–z plane, where we explore the behavior of the thick-to-thin-disk density normalization and the halo-to-disk density normalization.

This paper is organized as follows. In Section 2, we describe the APOGEE and Gaia data sets we employ in the analysis. Section 3 examines the velocity DF of the Galactic disk and halo, and describes the building of the data-driven model. We discuss the thick-disk normalization and the halo-to-disk density in Section 4, and the most relevant results of the present study are summarized and discussed in Section 5.

Our study makes use of the data products from Data Release 16 (DR16) of the APOGEE (Majewski et al. 2017). Part of both SDSS-III (Eisenstein et al. 2011) and SDSS-IV (Blanton et al. 2017) via APOGEE and APOGEE-2, respectively, the combined APOGEE enterprise has been in operation for nearly a decade, and, through the installation of spectrographs on both the Sloan (Gunn et al. 2006) and du Pont 2.5 m telescopes (Bowen & Vaughan 1973; Wilson et al. 2019), has procured high-resolution, H-band spectra for more than half a million stars across both hemispheres. The survey provides ${v}_{\mathrm{los}}$, stellar atmospheric parameters, and individual abundances for on the order of fifteen chemical species (Holtzman et al. 2015; Zamora et al. 2015; García Pérez et al. 2016). A description of the latest APOGEE data products from DR14 and DR16 can be found in Holtzman et al. (2018) and H. Jönsson et al. (2020, in preparation), respectively.

In this study, we also use data from the second data release (DR2) of the Gaia mission (Gaia Collaboration et al. 2018). This catalog provides full six-dimensional space coordinates for 7,224,631 stars: positions (α, δ), parallaxes (ϖ), proper motions (${\mu }_{\alpha }^{* }$, μ_δ), and radial line-of-sight velocities (${v}_{\mathrm{los}}$) for stars as faint as G = 13 (Cropper et al. 2018). The stars are distributed across the full celestial sphere. Gaia DR2 contains ${v}_{\mathrm{los}}$ for stars with effective temperatures in the range ∼3550–6900 K. The median uncertainty for bright sources (G < 14) is 0.03 mas for the parallax and 0.07 mas yr⁻¹ for the proper motions (Arenou et al. 2018; Lindegren et al. 2018). The precision of the Gaia DR2 ${v}_{\mathrm{los}}$ at the bright end is on the order of 0.2–0.3 km s⁻¹, while at the faint end, it is on the order of 1.4 km s⁻¹ for T_eff = 5000 K stars and ∼3.7 km s⁻¹ for T_eff = 6500 K. For further details about the Gaia DR2 subsample containing ${v}_{\mathrm{los}}$ measurements, we refer the reader to Sartoretti et al. (2018) and Katz et al. (2019). We follow the recommendations from Boubert et al. (2019) for studies in Galactic dynamics using Gaia to remove stars where the color photometry is suspect, as well as stars where the ${v}_{\mathrm{los}}$ measurement is based on fewer than four transits of the instrument.

Individual space velocities in a Cartesian Galactic system were obtained by following the equations in Johnson & Soderblom (1987). That is, from Gaia DR2 line-of-sight velocities, proper motions, and parallaxes, we derive the space velocity components (U, V, W). In the case of the APOGEE targets, we use the APOGEE-2 DR16 catalog ${v}_{\mathrm{los}}$ (Ahumada et al. 2019), where the internal precision is better than 0.1 km s⁻¹ (Nidever et al. 2015). When the relative uncertainty in parallaxes becomes larger, the inverse of a measured parallax is a biased estimate of the distance to a star (Strömberg 1927; Luri et al. 2018). For this reason, we select stars with positive parallaxes and a relative parallax uncertainty smaller than 20% (σ_ϖ/ϖ ≤ 0.2). We also set the Gaia flag astrometric-excess-noise = 0 to drop stars with poor astrometric fits. In addition, the flag ${rv}-{nb}-\mathrm{transits}\gt 4$ is set to secure enough Gaia transits that robust ${v}_{\mathrm{los}}$ measurements are in hand. Following Boubert et al. (2019), the selected Gaia stars for this study must have reported G_BP and G_RP magnitudes. That leaves 4,774,723 Gaia targets for this exercise. The remaining targets have high-precision Gaia parallaxes (the vast majority of the stars have σ_ϖ/ϖ ≤ 0.05), so that their distances can be determined by simple parallax inversion without significantly biasing the error derivation (e.g., Bailer-Jones et al. 2018 and references therein). We adopt a right-handed Galactic system, where U is pointing toward the GC, V in the direction of rotation, and W toward the North Galactic Pole. For the peculiar motion of the Sun, we adopt the values U_⊙ = 11.1 km s⁻¹, V_⊙ = 12.2 km s⁻¹, and W_⊙ = 7.2 km s⁻¹ (Schönrich et al. 2010). We also transform the velocities from Cartesian to a cylindrical Galactic system: υ_R, υ_ϕ, υ_z. The Sun is assumed to be located at X = 8.34 ± 0.16 kpc and the circular rotation speed at the location of the Sun to be 240 ± 8 km s⁻¹ (Reid et al. 2014). We define $R={({X}^{2}+{Y}^{2})}^{1/2}$ as the distance from the GC, projected onto the Galactic plane.

In the end, the typical uncertainty in the velocities used here is Δυ ∼ 1.5 km s⁻¹ per dimension. To remove outlier velocities that can yield unrealistic velocity dispersions, we select stars where (${\upsilon }_{R}^{2}+{({\upsilon }_{\phi }-240)}^{2}+{\upsilon }_{z}^{2}{)}^{1/2}$ < 600 km s⁻¹, which removes a total of 503 stars. In the next section, we explore the kinematical properties of the Milky Way using the data set just described.

Stars of different mass synthesize and, upon death, eject into the interstellar medium different chemical elements, on different timescales. The overall metallicity measured in a star's atmosphere, for the most part, represents an integral over star formation and chemical enrichment prior to that star's birth, while the abundances of individual elements can be used to track the ratio of recent to past star formation rates.

Thus, stellar chemical-abundance patterns offer a means to discriminate stellar populations that have experienced differing star formation and chemical-enrichment histories. To procure unbiased velocity DFs for individual stellar populations in the nearby disk, we exploit the [Fe/H] and [Mg/Fe] abundances measured from APOGEE spectra, the combination of which has been show to give very good discrimination of the thin disk, thick disk, and halo (e.g., Bensby et al. 2003; Hayes et al. 2018; Weinberg et al. 2019). We select from the APOGEE database those stars for which the signal-to-noise ratio (S/N) > 70 and no aspcapbad flag is set (Holtzman et al. 2015). The APOGEE survey has a number of focused science programs that target specific objects like the Sgr dSph galaxy (Majewski et al. 2013; Hasselquist et al. 2017), the Large and Small Magellanic Clouds (Nidever et al. 2020), numerous star clusters (Mészáros et al. 2020), and many other stellar astrophysics programs (Zasowski et al. 2013, 2017) that are not germane to this study of normal field stars. To remove these specialized targets from our database, we identified fields associated with the special programs and removed all targets within those fields.

We summarize our chemistry-based discrimination and selection of our three primary stellar populations in Figure 1. The higher [Fe/H] population with low-α abundances (−0.7 < [Fe/H] < +0.5, −0.1 < [Mg/Fe] < +0.17) is associated with the Galactic thin disk, and our selection contains a total of 211,820 of these stars. The [Fe/H] < 0 stars with high-α abundances (+0.18 < [Mg/Fe] < +0.4 within the metallicity range −1.0 < [Fe/H] < 0.0) we associate with the Galactic thick disk, and the number of stars we have in that sample is 52,709. Finally, the Galactic halo population is identified with stars having [Fe/H] < −1.0 (e.g., Hayes et al. 2018). The number of APOGEE stars in this population is 5795. The magnitude-limited Gaia data set used in this study contains different stellar populations associated with different structures in the Galaxy. This is evident in Figure 2, where we show the distribution of υ_R, υ_ϕ, and υ_z for the population-integrated Gaia sample and for the population-separated APOGEE-2 data. From a comparison of the upper and lower panels, it becomes obvious that an assignment of population membership to a particular star based exclusively on kinematical properties is not generally possible. Moreover, it is not possible on the basis of kinematical data alone to reliably determine even the relative contributions of the different populations to the net velocity DF on a statistical basis. Figure 2 shows that the velocity DF of the different Galactic components clearly overlap, but also that individual abundances from high-resolution spectroscopy surveys are a useful tool for apportioning stars to their relative stellar populations (Navarro et al. 2011; Nidever et al. 2014).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

While severe overlap is expected in the ${\upsilon }_{R}$ and υ_z dimensions, where all three stellar populations share a mean value around 0 km s⁻¹, we expect more separation of populations in υ_ϕ due to variations in asymmetric drift. Nevertheless, stars with υ_ϕ >70 km s⁻¹ dominate the distribution and yet have contributions from all three populations, though, naturally, most strongly from the thin and thick disks. We also observe that the thin-disk population shows a very small number of slow-rotating stars; these could represent a small contaminating portion of stars from the accreted halo, which can have metallicities as high as [Fe/H] ∼ −0.5 (Hayes et al. 2018).

Interestingly, for the population with [Fe/H] < −1.0, we observe an extended velocity DF with a peak around υ_ϕ ∼ 0 km s⁻¹, together with a second peak around 120 km s⁻¹ (see the middle panel in Figure 2). We discuss this metal-poor population in more detail in Section 3.1.3.

While there is no precedent for the accuracy and precision of the Gaia astrometry, and for such an enormous number of objects, there is still no detailed chemistry available for the vast majority of stars in the Gaia data set. Thus, one cannot yet leverage the huge statistical power of Gaia to explore the kinematical properties of individual stellar populations in an unbiased manner, nor use these kinematics to sort stars reliably into populations to help define other gross properties of the populations. However, we will now show how, with a data-driven model trained with the velocity DFs defined for the combined APOGEE+Gaia data set for relatively nearby stars, we can harness the power of the greater Gaia database without chemical data to sort stars into populations based on their kinematics and use these statistical memberships to ascertain other bulk properties of the populations—e.g., to assess the relative densities of stars in these populations over a broad range of the Galactic locations. We describe the steps toward these goals in the next sections.

3.1. APOGEE Data-driven Model

To build a simple, kinematical data-driven model, we do not assume a triaxial Gaussian DF (Schwarzschild 1907; Bienaymé 1999); instead, we use the velocity DF for the chemically selected thin-disk, thick-disk and halo Galactic components in the APOGEE data. The DF, $f({\boldsymbol{\upsilon }}$), is defined such that $f({\boldsymbol{\upsilon }}$) $d{\boldsymbol{\upsilon }}$ is the number of stars per unit volume with velocity in the range $[{\boldsymbol{\upsilon }}$, ${\boldsymbol{\upsilon }}$ + $d{\boldsymbol{\upsilon }}$]. We now explore the characteristics of the DF for each primary Galactic stellar population.

3.1.1. Chemically Selected Thin Disk

Figure 3 shows the velocity DF for the υ_R, υ_ϕ, and υ_z of the low-α sequence population selected in the [Fe/H]–[Mg/Fe] plane (see Figure 1). We also show the best Gaussian fit for each of the three components of velocity (red lines in Figure 3).

Zoom In Zoom Out Reset image size

By and large, these Gaussian fits are reasonable descriptors of the DFs. Even for the υ_ϕ component, where we expect a skew in the observed DF due to asymmetric-drift effects, a normal distribution nevertheless reproduces the distribution reasonably well. We find the largest discrepancies between the observed velocity DF and a Gaussian distribution in all three cases to be mainly at the very peaks and in the wings, and, for the latter, especially in the cases of υ_ϕ and υ_z. As a demonstration of the utility of Gaussians as descriptors of the DFs, we find that for the chemically selected thin disk, the DF skewness(υ_R)¹⁰  = 0.04 and the kurtosis(υ_R) = 0.61. For the azimuthal velocity, we have skewness(υ_ϕ) = −0.44 and kurtosis(υ_ϕ) = 0.98. Furthermore, for the vertical component of the velocity, we find skewness(υ_z) = −0.06 and a kurtosis(υ_z) = 1.33. These skewness and kurtosis values lie within the allowable range for normal univariate distributions; e.g., Darren & Mallery (2010) argue that values for asymmetry and kurtosis between −2 and +2 are indicative of normal univariate distributions, while even by the more conservative limits of −1.5 and +1.5 advocated by Pilkington et al. (2012), the Figure 3 DFs are found to be well-described as Gaussians.

For the low-α sequence population, we find the following velocity dispersion values, (σ_R, σ_ϕ, σ_z) = (36.81 ± 0.07, 24.35 ± 0.04, 18.03 ± 0.03) km s⁻¹. Based on these values, the shape of the velocity ellipsoid for the thin disk is found to be (σ_R:σ_ϕ:σ_z) = (1.00:0.66:0.49). We also find the vertex deviation for this population to be α_Rϕ = −401 ± 009 and the tilt of the velocity ellipsoid to be α_Rz = +141 ± 002; these quantities are defined so that α_ij corresponds to the angle between the i-axis and the major axis of the ellipse formed by projecting the three-dimensional velocity ellipsoid onto the ij plane, where i and j are any of the stellar velocities. We can also calculate the orbital anisotropy parameter (Binney 1980) in spherical polar coordinates; for the thin disk, we have β = +0.57 ± 0.01. These findings are summarized in Tables 1 and 2.

Table 1.  Summary of the Kinematical Properties for the Different Population Subsamples

	$\overline{{\upsilon }_{\phi }}$	σR	σϕ	σz	σϕ/σR	σz/σR
	(km s−1)	(km s−1)	(km s−1)	(km s−1)
Chemically selected thin disk	+229.43 ± 0.54	37.61 ± 0.07	25.01 ± 0.04	18.53 ± 0.03	0.66 ± 0.01	0.49 ± 0.01
Chemically selected thick disk	+191.82 ± 0.24	64.68 ± 0.20	50.82 ± 0.15	43.60 ± 0.13	0.78 ± 0.01	0.67 ± 0.01
Halo ([Fe/H] < −1)	+35.53 ± 1.28	150.57 ± 1.58	115.67 ± 1.21	86.67 ± 0.91	0.77 ± 0.02	0.57 ± 0.01
Disk-like rotation	+186.10 ± 1.66	59.23 ± 1.39	50.00 ± 1.17	47.55 ± 0.85	0.84 ± 0.03	0.80 ± 0.02
Nonrotation	−2.35 ± 1.57	165.51 ± 1.94	95.12 ± 1.11	94.10 ± 1.16	0.57 ± 0.01	0.57 ± 0.01

Download table as:  ASCIITypeset image

Table 2.  Tilt of the Velocity Ellipsoid, Vertex Deviation, and Orbital Anisotropy Parameter for the Different Populations

	${\alpha }_{\mathrm{Rz}}$	αRϕ	β
	(degrees)	(degrees)
Chemically selected thin disk	+1.41 ± 0.02	−4.01 ± 0.09	+0.57 ± 0.01
Chemically selected thick disk	+5.16 ± 0.15	−6.01 ± 0.47	+0.32 ± 0.01
Halo ([Fe/H] < −1)	+11.43 ± 0.35	−2.30 ± 0.15	+0.20 ± 0.02
Disk-like rotation	+0.86 ± 0.38	+1.18 ± 0.59	+0.56 ± 0.02
Nonrotation	+0.15 ± 0.02	+0.16 ± 0.02	+0.42 ± 0.02

Download table as:  ASCIITypeset image

3.1.2. Chemically Selected Thick Disk

In Figure 4, the velocity DF for the three individual space velocities is shown for the high-α sequence population defined in Figure 1. For this population, we find a mean rotational velocity of $\overline{{\upsilon }_{\phi }}$ = 191.82 ± 0.24 km s⁻¹ and velocity dispersion of (σ_R, σ_ϕ, σ_z) = (62.44 ± 0.21, 44.95 ± 0.15, 41.45 ± 0.15) km s⁻¹. As in the case for the thin disk, we find that normal distributions (red lines in Figure 4) give a good match for the velocity components υ_R and υ_z. For example, we calculate skewness = 0.03 and kurtosis = 0.62 for the radial velocity component, while for the vertical velocity, we have 0.02 and 0.86, respectively.

Zoom In Zoom Out Reset image size

However, for the rotational velocity (υ_ϕ) of the thick-disk population, the spread in asymmetric drift is more prominent than for the thin disk, so that the υ_ϕ DF is skewed to low velocities. For this velocity component, we have skewness = −1.38 and kurtosis = −4.46. Nevertheless, remarkably, it is possible to account for this skewness adequately with only a single other normal distribution (gray line) fit to the υ_ϕ DF (see middle panel in Figure 4). While the simplicity of a two-Gaussian fit is very convenient, it is not clear whether it has a physical meaning representing two distinct subpopulations, or that the mathematical contrivance just happens to work well as an explanation for what could be a more complex combination of subpopulations. If it really has a physical meaning as an apparent second subpopulation in υ_ϕ, this less-dominant subpopulation has Gaussian parameters μ ∼ 92 km s⁻¹ and σ ∼ 115 km s⁻¹, and could represent the oldest population in the thick disk, which is, on average, lagging some 100 km s⁻¹ behind a more rapidly rotating and significantly dynamically colder and presumably younger thick-disk population (green distribution in Figure 4).

We find differences in the shape of the velocity ellipsoid for the thick disk with respect to the thin disk (Table 1), where (σ_R:σ_ϕ:σ_z) = (1.00:0.78:0.67). For purposes of this calculation, we treat the total thick-disk sample as one population; this is certain to result in some increase in the measured dispersion in the υ_ϕ direction.

Theoretical studies in the formation of the thick disk predict a wide range of values for σ_z/σ_R. For example, Villalobos et al. (2010) found a range from ∼0.4 to 0.9 for a model with formation of the thick disk through heating due to accretion events. Assuming that a merger led to the dynamical heating of a preexisting precursor disk to the thick disk (e.g., Haywood et al. 2018; Di Matteo et al. 2019; Gallart et al. 2019 and references therein), our findings may be suggestive of an encounter with a satellite on a low/intermediate orbital inclination.

Interestingly, we find the tilt of the velocity ellipsoid for the thick disk to be larger than that found for the thin disk; however, our tilt α_Rz = +516 ± 015 is lower than previous values reported for the thick disk. For example, Casetti-Dinescu et al. (2011) found α_Rz = +86 ± 18, where the authors selected the thick-disk population using stellar density laws and the vertical height with respect to the Galactic plane. With a sample of ∼1200 red giants, Moni Bidin et al. (2012) found an even larger value, α_Rz = +100 ± 05. They created the thick-disk sample by selecting stars with ∣z∣ > 1.3 kpc, and they used the individual velocities of each star to remove outliers they associated with the halo. However, a direct comparison of our results with previous measurements from the literature is difficult, because the tilt angle varies with z (Moni Bidin et al. 2012; Hayden et al. 2020).

We also find the orbital anisotropy parameter for the Galactic thick disk, β = +0.32 ± 0.01, to be lower than the value found for the thin disk, which suggests that the orbital anisotropy is mildly radial for this population.

3.1.3. Halo

The stellar density of halo stars, simply defined here by stars having [Fe/H] < −1.0, is significantly lower than that for the disk-like populations in our volume of study; hence, the selection effects driven by the APOGEE survey design (Zasowski et al. 2017) may have a larger impact on this population. The black lines in Figure 5 show the velocity DF for the halo population. For the distribution in υ_ϕ, we use a mixture of two Gaussian components. The metal-poor stars with disk-like kinematics (gray Gaussian in the middle panel of Figure 5) may be associated with the thick disk, while the components with a nonrotation average motion may be associated with the inner halo (see Chiba & Beers 2000; Carollo et al. 2010; Hayes et al. 2018; Mackereth et al. 2019 and references therein for more details). The stars in the entire halo sample ([Fe/H] < −1.0) are characterized by a radially elongated velocity ellipsoid, where we have (σ_R, σ_ϕ, σ_z) = (150.57 ± 1.52, 115.67 ± 1.08, 86.67 ± 0.93) km s⁻¹. For the entire halo population defined in this study (red lines in Figure 5), we find a small mean prograde rotation of 35 km s⁻¹. These results are in good agreement with the halo properties reported in Chiba & Beers (2000) and Bond et al. (2010), who also defined the halo by simple metallicity cuts. Finally, for the entire [Fe/H] < −1.0 halo sample, we find a large angle for the tilt of the velocity ellipsoid, α_Rz = +1143 ± 035, while β is found to be nearly isotropic for these stars.

Zoom In Zoom Out Reset image size

Meanwhile, the "halo" subpopulation with disk-like kinematics in our sample shows a mean rotational velocity of $\overline{{\upsilon }_{\phi }}$ = 186.10 ± 1.36 km s⁻¹ (gray Gaussian in Figure 5); this velocity is very similar to the mean rotational velocity we found for the Galactic thick disk (see the middle panel in Figure 4). For the velocity dispersion of this subpopulation, we find (σ_R, σ_ϕ, σ_z) = (59.25 ± 1.12, 50.00 ± 0.91, 47.52 ± 0.94) km s⁻¹. These results are consistent with the values reported for the chemically selected thick disk (see Section 3.1.2), supporting previous results showing that the thick disk might exhibit an extended metal-poor tail—more metal poor than [Fe/H] = −1.0 and reaching values of [Fe/H] ∼ −1.5 or even lower (i.e., Norris 1986; Morrison et al. 1990; Beers et al. 2002; Nissen & Schuster 2010; Robin et al. 2014; Fernández-Alvar et al. 2018; Hayes et al. 2018; An & Beers 2020).

For the other, "more traditional" halo component with a nonrotating average motion of $\overline{{\upsilon }_{\phi }}=-2.31\pm 1.59$ km s⁻¹ (green Gaussian in Figure 5), we find that the velocity dispersion is (σ_R, σ_ϕ, σ_z) = (165.52 ± 1.97, 95.10 ± 1.12, 94.14 ± 1.14) km s⁻¹. This halo population has been extensively discussed in the literature (e.g., Searle & Zinn 1978; Ratnatunga & Freeman 1985; Gilmore & Wyse 1998; Gratton et al. 2003; Di Matteo et al. 2019, and references therein).

These properties for the hotter halo component are consistent with current descriptions of the halo as an accreted population of the Milky Way. For example, Nissen & Schuster (2010) used ∼100 halo stars in the solar neighborhood to identify a "low-α" population, for which the kinematics suggest that it may have been accreted from dwarf galaxies (Zolotov et al. 2009; Purcell et al. 2010), some specifically originating from the ω Cen progenitor galaxy (Bekki & Freeman 2003; Meza et al. 2005; Majewski et al. 2012). Meanwhile, Belokurov et al. (2018) studied the orbital anisotropy for the local stellar halo, and, by comparing the observational results with cosmological simulations of halo formation, concluded that the inner halo was deposited in a major accretion event by a satellite with M_vir > 10¹⁰ M_⊙, inconsistent with a continuous accretion of dwarf satellites. Belokurov et al. (2018) also highlighted the non-asymmetric structure of the remains of the merger in the velocity distribution. Their findings echo discussions presented in simulations including a similar accretion event by Meza et al. (2005), as well as the exploration of the Venn et al. (2004) data set by Navarro et al. (2011), where the locally sampled velocity distribution of a high-energy accretion event appears to produce a mixture of kinematical populations.

In addition to the Belokurov et al. (2018) study, Helmi et al. (2018) selected the retrograde halo population using APOGEE DR14 (Abolfathi et al. 2018) and Gaia DR2 to conclude that the inner halo is dominated by debris from an object that, at infall, was slightly more massive than the Small Magellanic Cloud. They also argued that the merger must have led to the dynamical heating of the precursor of the Galactic thick disk approximately 10 Gyr ago. The dwarf galaxy progenitor of this debris—now variously called "Gaia-Enceladus," the "Gaia Sausage," or collectively the "Gaia-Enceladus-Sausage"—is thought to have fallen into the Milky Way on a highly eccentric orbit (e ∼ 0.85), accounting for the predominance of stars with such radial orbits in the inner Galaxy (Belokurov et al. 2018; Haywood et al. 2018; Helmi et al. 2018; Fattahi et al. 2019).

3.2. The 3D Velocity DF as a Function of R and z

In a larger context, our chemical-selection approach also allows us to explore the 3D velocity DF as a function of R and z for the different subpopulations. Figure 6 compares the regions of the sky probed by the APOGEE (gray line) and Gaia (black line) samples. The samples probe a region between 4 < R < 13 kpc, but the majority of stars are located at distances between 7 and 10.5 kpc. The spatial distribution for the vertical height shows that most of the selected stars have $| z|$ < 2.0 kpc (right panel in Figure 6). The spatial distributions for the two samples are remarkably similar, despite the fact that the APOGEE survey works at brighter magnitudes but is predominantly targeted at giant stars, whereas the larger Gaia sample studied here probes fainter magnitudes, is parallax-error limited, and is dominated by dwarf stars. Differences in the spatial distributions mainly reflect the fact that Gaia is an all-sky survey while the APOGEE DR16 subsample is a pencil-beam survey dominated by northern hemisphere observations (where APOGEE has been surveying longer), which explains why we see an asymmetry at larger R and an excess for z > 0 with respect to Gaia (see right panel in Figure 6).

Zoom In Zoom Out Reset image size

Figure 7 shows the velocity DF for three different ranges in Galactocentric radius and vertical height for the chemically selected thin disk (top panels), thick disk (middle panels), and halo (bottom panels). We also show the cumulative fraction for each distribution. The different colors show the velocity DFs for three ranges in Galactocentric radius and vertical height.

Zoom In Zoom Out Reset image size

In most cases, we do not find large discrepancies between the velocity DF for the three individual components across different ranges of R and z, especially for the thin disk and thick disk. The negligible spatial variations in the velocity DF for the different, chemically selected populations justifies and motivates our strategy to use the velocity distributions discussed in Sections 3.1.1–3.1.3 to build a data-driven model of Galactic stellar populations.

However, it is notable that we observe a small asymmetric-drift variation with z for the thick disk (right figure in the middle panels in Figure 7). We also find that the innermost part of the thin disk tends to lag with respect to the rest of the population at larger Galactocentric radii (left figure in the upper panels in Figure 7). On the other hand, for the [Fe/H] < −1 group of stars, we find strong variations in the contributions of the two components in υ_ϕ described in Section 3.1.3, with the thick-disk-like population more prominent in the inner region (R < 7.3 kpc) and closer to the Galactic plane ($| z|$ < 0.5 kpc), and the "inner halo" population with a nonrotating average motion dominating the distribution for the outer regions in R and z (see bottom panels in Figure 7).

3.3. The Stellar Disk and Halo Contribution in the Gaia Sample

To carry out an unbiased study of the spatial distributions of the Galactic components for Gaia stars, the bulk of which lack abundance information, we use the results presented above to model the contribution from the thin-disk, the thick-disk, and the halo stars in the much larger Gaia sample. The Gaia velocity DF observed in the three components of space velocity (see Figure 2) is a combination of the different components of the stellar disk and the Galactic halo. Using the APOGEE velocity DF measured for the chemically distinguished disk and halo components, we can estimate the fraction of stars that belong to each of these components in the much larger Gaia data set.

In the previous section, we verified that the APOGEE subsample is a legitimate proxy for the larger Gaia sample. In particular, we checked which regions of the Galaxy are probed by the Gaia sample employed in this exercise as well as the APOGEE sample and find (Figure 6) that they survey comparable volumes.

To assess further the ability of the derived data-driven model to quantify the number of objects in the Gaia sample that are attributable to the thin disk, thick disk, and halo, we apply an Anderson–Darling-style statistical test sensitive to discrepancies at low and high values of ${\boldsymbol{\upsilon }}$ (Press et al. 1992), adopted as follows:

$$\begin{eqnarray*}&&\langle {S}_{G}({\boldsymbol{\upsilon }})-{S}_{A}({\boldsymbol{\upsilon }}){\rangle }^{2}/{S}_{A}({\boldsymbol{\upsilon }})(1-{S}_{A}({\boldsymbol{\upsilon }})).\end{eqnarray*}$$

The Anderson–Darling statistic is based on the empirical cumulative DF. The cumulative distribution for the Gaia sample (S${}_{G}({\boldsymbol{\upsilon }})$) is calculated directly from the individual velocities in the Gaia data set. For the cumulative distribution calculated using the velocity DF from the APOGEE data set (S ${}_{A}({\boldsymbol{\upsilon }})$), we have a mixture of different distributions for the thin disk, thick disk, and halo that we draw directly from the APOGEE data (black distributions in Figures 3–5). The sum of the three stellar velocity distributions associated with the three Galactic components allows us to derive estimates, guided by Anderson–Darling statistical tests, of the fraction of objects that are part of the thin disk, thick disk, and halo in the Gaia data set. We can write $f{({\boldsymbol{\upsilon }})}_{\mathrm{gaia}}={\sum }_{i=1}^{3}{\epsilon }_{i}f{({\boldsymbol{\upsilon }})}_{\mathrm{apogee}}$, where i represents each velocity component, and _i = (_thin, _thick, _halo) is the fraction of stars associated with that Galactic component. Following the results in Section 3.2, we assume that the shape of the velocity distribution for each component does not change in the local volume of study; what changes is the fraction of stars in each Galactic component.

Figure 8 shows the velocity DF for the three velocity components. The thick black line is the Gaia DR2 data set employed in this study; red is the data-driven model built from the chemically selected thin disk, thick disk, and halo. That the thick black Gaia data lines in Figure 8 are almost completely obscured by the red model lines is a testament to the veracity of the model to explain the data. For the υ_R component, we find the fraction of Galactic components (thin disk, thick disk, and halo, respectively)  = (0.777, 0.208, 0.015), while for υ_ϕ, we find  = (0.881, 0.103, 0.016). Finally, for the vertical component, υ_z, we obtain  = (0.799, 0.186, 0.015). Figure 8 also shows the decomposition into the thin-disk, thick-disk, and halo velocity DFs (gray lines) from the data-driven model showing the estimation of the fraction of objects that are part of the different Galactic structures. The radial and vertical velocity DFs give similar fractions for the three components; however, the υ_ϕ velocity yields a larger number of thin-disk stars compared to the thick disk. It is likely that this difference has to do with the fact that while the distribution of υ_R and υ_z can be reproduced with a Gaussian function (see Sections 3.1.1 and 3.1.2), the distribution of υ_ϕ is strongly non-Gaussian. The latter is highly skewed to low velocities due to the asymmetric drift, especially for the Galactic thick disk (as seen in Figures 4 and 5), and this non-asymmetric behavior is more challenging for our model to describe.

Zoom In Zoom Out Reset image size

Nevertheless, by combining the results calculated using the three individual velocities and taking the average value, we estimate that 81.9% ± 3.1% of the objects in the selected Gaia data set are thin-disk stars, 16.6% ± 3.2% are thick-disk stars, and 1.5% ± 0.1% belong to the Milky Way halo.

In this section, we determine the local density normalization (f_ρ = ρ_T/ρ_t) of the thick disk compared to the thin disk using the velocity DF derived above. Many previous derivations of f_ρ were based on star-count data derived from photometric parallaxes (e.g., Gilmore & Reid 1983; Jurić et al. 2008), with the estimates ranging from 1% to 12%. Bland-Hawthorn & Gerhard (2016) analyzed the results from 25 different photometric surveys conducted since the discovery of the thick disk (Yoshii 1982) and concluded that f_ρ = 4% ± 2% at the solar circle. However, star-count approaches to this problem are subject to the degeneracy between the derived scale lengths and scale heights for each of the thin-disk and thick-disk components, which drives a large uncertainty in f_ρ (Chen et al. 2001; Siegel et al. 2002). Our novel approach, using the observed velocity DFs for the thin disk, thick disk, and halo from APOGEE as a data-driven model applied to the Gaia data set, is not affected by this degeneracy. Moreover, we are able to analyze the behavior of f_ρ in the R–z plane, averaged over the Galactocentric polar angle, ϕ.

To build a given R–z pixel, we create intervals of 0.1 kpc in position, and use only pixels where the number of stars for a given interval is N(R_i, z_i) ≥ 50. Figure 9 shows the stellar number density as a function of Galactic cylindrical coordinates R and z created through this exercise. The largest number of stars per pixel in the magnitude-limited Gaia sample is within −1 < z < 1 kpc and 6 < R < 10 kpc. Using the data-driven model described in Section 3.1, we create different velocity DFs, where our free parameters are the fraction of thin-disk, thick-disk, and halo stars in intervals of 1% for each population. For a given R–z pixel, we use the observed velocity DF in the area of study, and we perform an Anderson–Darling-style statistical test (as was done in Section 3.3). The minimum value in the statistical test (i.e., the maximum deviation between the cumulative distribution of the data and the data-driven model function) is used to estimate the fraction, and hence the density normalization, of each population. We perform this analysis for the three individual velocity components separately. Figure 10 shows the results from the statistical Anderson–Darling test for the three velocities (top panel). We find that the distributions are very similar regardless of the velocity components used. The bottom panel of Figure 10 also shows the Anderson–Darling test values for the three individual velocities added in quadrature in the R–z plane. R–z pixels with lower numbers of stars tend to have larger Anderson–Darling test values, i.e., a larger deviation between the data and the model.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Once we know the fraction of stars that belong to the thin disk, thick disk, and halo, calculation of the thick-disk normalization, f_ρ, in the R–z plane is straightforward. For the final value, we average between the f_ρ given for each velocity component. We present the results in Figure 11, where the R–z plane is color-coded by f_ρ. These figures, for the first time, reveal in detail how f_ρ varies across the Galaxy.

Zoom In Zoom Out Reset image size

From inspection of Figure 11, one can clearly see where the thin-disk population dominates (dark blue), where the f_ρ values range from nearly zero to ∼15%. We also see a transition region, where we have a more even mix between thin-disk and thick-disk populations (green), indicated by where the maps have f_ρ ∼ 50%. Finally, we also observe the region where the population that belongs to the Galactic thick disk dominates (red; see the right panel in Figure 11), shown by where f_ρ > 85%. Interestingly, we observe that the density ratio in the transition region shows differences between the north and south. The area around R ∼ 10 kpc does not show the same transition in density at z ∼ −1.5 kpc as at z ∼ 1.5 kpc. This could be an effect of the Galactic warp (e.g., Levine et al. 2006; Romero-Gómez et al. 2019 and references therein), but this could also be related to excitation of wave-like structures creating a wobbly galaxy (Widrow et al. 2012; Williams et al. 2013). Furthermore, Figure 11 shows a radial trend, in the sense that at increasing R the thin disk dominates at larger z (i.e., the blue region in Figure 11 gets thicker at larger Galactocentric radius). This phenomenon could be interpreted as thin-disk flaring (e.g., Feast et al. 2014; Thomas et al. 2019), but also can be related to the different scale lengths of the thin disk with respect to the thick disk, where the scale length for the thick disk is shortened with respect to the one for the thin disk. The latter suggestion is in line with the recent literature on this debated topic. For example, Bensby et al. (2011) and Bovy et al. (2012) suggested a shorter scale length for the chemically distinguished high-α disk compared to the low-α disk.

Moreover, in the right panel of Figure 11, we have the same Gaia data R–z plane color-coded by f_ρ, but in this case, the figure is restricted from 0% to 15%. The visualization of this new range in f_ρ allows us to see smaller differences in the disk normalization for a given position. We suggest that vertical oscillations in the disk (e.g., Chequers et al. 2018) may be responsible for these small fluctuations in f_ρ for the thin Galactic disk by introducing small changes in the shape of the velocity DF in such a dynamically active disk.

In the end, using our novel approach based on velocities, rather than star counts, we find the local thick-to-thin-disk density normalization to be ρ_T(R_⊙)/ρ_t(R_⊙) = 2.1% ± 0.2%.

Figure 12 shows the halo-to-disk density normalization. We find that the halo is most dominant in the Gaia data set in the region R < 8 kpc and ∣z∣ > 1.3 kpc. We find the local halo-to-disk density normalization to be ρ_H(R_⊙)/(ρ_T(R_⊙) + ρ_t(R_⊙)) = 1.2% ± 0.6%.

Zoom In Zoom Out Reset image size

Combining the precise stellar abundances from the APOGEE survey with the astrometry from Gaia, we study the velocity DF for chemically selected low-α (thin-disk), high-α (thick-disk), and halo ([Fe/H] < −1) stars. Using the kinematical properties of these subsamples, we built a data-driven model and used it to dissect a 20% parallax-error-limited Gaia sample to understand the contribution of the different Galactic structural components to the velocity-space DF as a function of Galactic cylindrical coordinates R and z. We find that 81.9% ± 3.1% of the objects in the selected Gaia data set are thin-disk stars, 16.6% ± 3.2% are thick-disk stars, and 1.5% ± 0.1% belong to the Milky Way halo.

The local fraction of the Milky Way thick disk, ρ_T(R_⊙)/ρ_t(R_⊙), is still under debate, as evidenced by the large spread in derived values—ranging from 2% (Gilmore & Reid 1983) to 12% (Jurić et al. 2008)—from star-counting methods, which are notoriously fraught with model degeneracies. Our analysis, based on the velocity characteristics of chemically selected populations, helps to break the aforementioned degeneracies and favors lower values for the normalization (e.g., Bovy et al. 2012). We find the local thick-to-thin-disk density normalization to be ρ_T(R_⊙)/ρ_t(R_⊙) = 2.1% ± 0.2%, a result consistent with the lower end of values derived using star-count/density analyses, but determined in a completely different way.

Using the same methodology, the local halo-to-disk density normalization is found to be ρ_H(R_⊙)/(ρ_T(R_⊙) + ρ_t(R_⊙)) = 1.2% ± 0.6%.

This is several times larger than other values found recently using kinematically selected samples. For example, Amarante et al. (2020), using high tangential velocity stars (υ_t > 200 km s⁻¹) in Gaia DR2, found the local halo-to-disk density normalization to be 0.47%. This is in agreement with the value (0.45%) found in Posti et al. (2018), where the local stellar halo was selected through action angle distributions using Gaia DR1 and RAVE. The differences between these kinematics-based halo samples and our chemically selected halo ([Fe/H] < −1.0) might be explained by our halo sample, including a non-negligible fraction of metal-poor, but kinematically colder stars that may be related to the thick disk. Their inclusion may inordinately elevate our derived halo-to-disk density normalization.

While overall a chemical discrimination of stellar populations, as undertaken here, does a remarkably good job and provides many benefits over other population-separation methods, ultimately our approach is limited by the degree to which the metal-weak thick-disk and halo populations overlap in chemical-abundance spaces like that shown in Figure 1. In the future, this chemical overlap may be overcome by either the use of a larger number of chemical dimensions or the combined use of kinematics and chemistry.

B.A. acknowledges helpful conversations with B. K. Gibson, K. C. Freeman, I. Minchev, and A. C. Robin, as well as comments from the anonymous referee that improved the manuscript. B.A., S.R.M., C.R.H., and X.C. appreciate funding from NSF grant AST-1909497. T.C.B. acknowledges partial support from grant PHY 14–30152 (Physics Frontier Center/JINA-CEE), awarded by the U.S. National Science Foundation. Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the US Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS website is www.sdss.org.

SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU)/University of Tokyo, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatório Nacional/MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

This publication made use of NASA's Astrophysics Data System. 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. This research has made use of the Spanish Virtual Observatory (http://svo.cab.inta-csic.es), supported from the Spanish MINECO/FEDER through grant AyA2014-55216.