Stellar Inventory of the Solar Neighborhood using Gaia DR1

TL;DR

This study uses Gaia DR1 data to map the stellar density, luminosity, and mass functions in the solar neighborhood, providing detailed profiles and confirming the feasibility of future comprehensive galactic mapping.

Contribution

It presents the first detailed inventory of local stellar populations using Gaia DR1, including density profiles, luminosity functions, and mass functions, with analysis of Gaia's selection effects.

Findings

Vertical stellar density profiles follow sech^2 distributions.

Total mid-plane stellar density is approximately 0.040 M_sun/pc^3.

High-mass end of the present-day mass function follows a power law with index -4.7.

Abstract

The absolute number and the density profiles of different types of stars in the solar neighborhood are a fundamental anchor for studies of the initial mass function, stellar evolution, and galactic structure. Using data from the Gaia DR1 Tycho-Gaia Astrometric Solution, we reconstruct Gaia's selection function and we determine Gaia's volume completeness, the local number density, and the vertical profiles of different spectral types along the main sequence from early A stars to late K stars as well as along the giant branch. We clearly detect the expected flattening of the stellar density profile near the mid-plane for all stellar types: All vertical profiles are well represented by sech^2 profiles, with scale heights ranging from ~50 pc for A stars to ~150 pc for G and K dwarfs and giants. We determine the luminosity function along the main sequence for M_V < 7 (M >~ $0.72 M_\odot$) and along the giant branch for M_J >~ -2.5. Converting this to a mass function, we find that the high-mass (M > $1\,M_\odot$) present-day mass function along the main sequence is d n / d M = 0.016 $(M/M_\odot)^{-4.7}$ stars/pc^3/$M_\odot$. Extrapolating below M = $0.72\,M_\odot$, we find a total mid-plane stellar density of 0.040+/-0.002 $M_\odot$/pc^3. Giants contribute 0.00039+/-0.00001 stars/pc^3 or about 0.00046+/-0.00005 $M_\odot$/pc^3. The star-formation rate surface density is \Sigma(t) = 7+/-1 exp(-t/[7+/-1 Gyr]) $M_\odot$/pc^2/Gyr. Overall, we find that Gaia DR1's selection biases are manageable and allow a detailed new inventory of the solar neighborhood to be made that agrees with and extends previous studies. This bodes well for mapping the Milky Way with the full Gaia data set.

Tables

Table 1. Table 1: Properties, luminosity function, and mass function of stars along the main-sequence.

Spectral type	$J-K_s$	$M_V$	$M_J$	Mass	$n$		$dn/d{M}_V$		$dn/dM$		$\rho \left(>M\right)$		$z_d$
				$M_{\odot }$	$\left({10}^{-4}{\mathrm{pc}}^{-3}\right)$		$\left({10}^{-4}{\mathrm{pc}}^{-3}{\mathrm{mag}}^{-1}\right)$		$\left({10}^{-4}{\mathrm{pc}}^{-3}{M}_{\odot }^{-1}\right)$		$\left({10}^{-4}{M}_{\odot }{\mathrm{pc}}^{-3}\right)$		(pc)
A0V	​​​​​$-0.004$	1.34	1.25	2.15	0.38$\pm$	0.04	1.4$\pm$	0.2	2.2$\pm$	0.3	0.8$\pm$	0.1	49$\pm$	6
A1V	0.007	1.48	1.32	2.05	0.38$\pm$	0.04	2.0$\pm$	0.2	3.0$\pm$	0.3	1.6$\pm$	0.2	51$\pm$	6
A2V	0.024	1.55	1.35	2.00	0.41$\pm$	0.04	3.9$\pm$	0.3	5.5$\pm$	0.5	2.4$\pm$	0.2	51$\pm$	6
A3V	0.032	1.76	1.47	1.90	0.69$\pm$	0.04	4.9$\pm$	0.3	9.1$\pm$	0.6	3.7$\pm$	0.3	46$\pm$	3
A4V	0.059	1.84	1.51	1.85	0.75$\pm$	0.05	5.2$\pm$	0.3	10.1$\pm$	0.6	5.1$\pm$	0.4	52$\pm$	3
A5V	0.069	1.89	1.54	1.83	0.35$\pm$	0.04	5.4$\pm$	0.6	10.1$\pm$	1.1	5.8$\pm$	0.5	37$\pm$	6
A6V	0.074	2.07	1.64	1.76	0.60$\pm$	0.04	5.2$\pm$	0.4	13.4$\pm$	0.9	6.8$\pm$	0.6	53$\pm$	4
A7V	0.095	2.29	1.78	1.67	1.01$\pm$	0.05	5.0$\pm$	0.3	12.6$\pm$	0.7	8.5$\pm$	0.6	55$\pm$	3
A8V	0.117	2.30	1.78	1.67	0.63$\pm$	0.04	5.5$\pm$	0.3	14.0$\pm$	0.8	9.6$\pm$	0.7	52$\pm$	4
A9V	0.121	2.51	1.90	1.59	0.65$\pm$	0.03	5.9$\pm$	0.3	16.1$\pm$	0.9	10.6$\pm$	0.8	56$\pm$	4
F0V	0.143	2.79	2.13	1.50	1.19$\pm$	0.06	4.8$\pm$	0.2	14.0$\pm$	0.7	12.4$\pm$	0.8	64$\pm$	3
F1V	0.166	2.99	2.26	1.44	1.41$\pm$	0.06	5.9$\pm$	0.2	18.9$\pm$	0.8	14.4$\pm$	0.9	67$\pm$	3
F2V	0.190	3.08	2.32	1.43	1.23$\pm$	0.05	8.5$\pm$	0.4	35.1$\pm$	1.5	16.2$\pm$	1.0	64$\pm$	7
F3V	0.198	3.23	2.42	1.39	0.88$\pm$	0.04	7.4$\pm$	0.3	35.3$\pm$	1.7	17.4$\pm$	1.1	75$\pm$	4
F4V	0.211	3.40	2.55	1.33	1.48$\pm$	0.05	9.3$\pm$	0.3	29.6$\pm$	1.1	19.4$\pm$	1.1	62$\pm$	7
F5V	0.227	3.70	2.78	1.25	2.54$\pm$	0.07	10.8$\pm$	0.3	36.3$\pm$	1.0	22.5$\pm$	1.2	80$\pm$	3
F6V	0.256	3.87	2.90	1.21	3.19$\pm$	0.08	13.6$\pm$	0.3	53.2$\pm$	1.3	26.4$\pm$	1.3	89$\pm$	4
F7V	0.273	4.01	3.01	1.18	2.41$\pm$	0.07	15.5$\pm$	0.4	68.8$\pm$	1.9	29.2$\pm$	1.4	97$\pm$	5
F8V	0.286	4.15	3.11	1.14	2.54$\pm$	0.06	18.1$\pm$	0.5	72.5$\pm$	1.8	32.1$\pm$	1.5	104$\pm$	6
F9V	0.300	4.45	3.34	1.08	4.46$\pm$	0.09	20.3$\pm$	0.4	89.3$\pm$	1.8	36.9$\pm$	1.6	118$\pm$	8
G0V	0.329	4.50	3.38	1.07	4.10$\pm$	0.11	23.4$\pm$	0.6	117.0$\pm$	3.1	41.3$\pm$	1.7	108$\pm$	11
G1V	0.335	4.79	3.57	1.02	4.21$\pm$	0.09	24.7$\pm$	0.5	140.2$\pm$	3.1	45.6$\pm$	1.8	177$\pm$	17
G2V	0.366	4.86	3.64	1.00	4.76$\pm$	0.12	26.4$\pm$	0.7	135.9$\pm$	3.5	50.4$\pm$	1.9	123$\pm$	20
G3V	0.373	4.94	3.70	0.99	2.05$\pm$	0.07	27.4$\pm$	0.9	136.9$\pm$	4.6	52.4$\pm$	2.0	133$\pm$	39
G4V	0.382	4.98	3.74	0.98	1.49$\pm$	0.07	24.8$\pm$	1.1	148.9$\pm$	6.7	53.9$\pm$	2.0	—
G5V	0.386	5.13	3.84	0.97	2.51$\pm$	0.07	26.4$\pm$	0.8	251.2$\pm$	7.3	56.3$\pm$	2.1	—
G6V	0.403	5.18	3.88	0.96	2.59$\pm$	0.08	25.9$\pm$	0.8	259.1$\pm$	7.6	58.8$\pm$	2.2	—
G7V	0.409	5.32	3.97	0.94	2.39$\pm$	0.05	25.1$\pm$	0.5	159.0$\pm$	3.5	61.0$\pm$	2.2	—
G8V	0.424	5.55	4.14	0.90	4.67$\pm$	0.14	25.2$\pm$	0.8	155.5$\pm$	4.7	65.2$\pm$	2.4	—
G9V	0.452	5.76	4.29	0.87	5.22$\pm$	0.14	23.7$\pm$	0.6	149.1$\pm$	4.1	69.8$\pm$	2.5	—
K0V	0.478	5.91	4.47	0.85	4.01$\pm$	0.12	22.3$\pm$	0.7	160.4$\pm$	4.7	73.2$\pm$	2.6	—
K1V	0.499	6.19	4.57	0.82	4.65$\pm$	0.12	21.6$\pm$	0.6	185.9$\pm$	4.8	77.0$\pm$	2.7	—
K2V	0.531	6.57	4.76	0.78	8.89$\pm$	0.20	26.9$\pm$	0.6	254.1$\pm$	5.7	83.9$\pm$	2.8	—
K3V	0.600	7.04	4.98	0.73	13.72$\pm$	0.29	32.3$\pm$	0.7	304.9$\pm$	6.4	93.9$\pm$	3.0	—
K4V	0.669	7.25	5.11	0.72	10.87$\pm$	0.35	32.0$\pm$	1.0	362.4$\pm$	11.7	101.8$\pm$	3.3	—

Table 2. Table 2: Luminosity function of giants.

$M_J$	$n$		$dn/d{M}_J$
	$\left({10}^{-4}{\mathrm{pc}}^{-3}\right)$		$\left({10}^{-4}{\mathrm{pc}}^{-3}{\mathrm{mag}}^{-1}\right)$		$\left({10}^{-4}{\mathrm{pc}}^{-3}\right)$
$-2.375$	0.05$\pm$	0.01	0.21$\pm$	0.05	0.05$\pm$	0.01
$-2.125$	0.05$\pm$	0.01	0.19$\pm$	0.03	0.10$\pm$	0.01
$-1.875$	0.06$\pm$	0.01	0.24$\pm$	0.05	0.16$\pm$	0.02
$-1.625$	0.05$\pm$	0.02	0.22$\pm$	0.07	0.21$\pm$	0.03
$-1.375$	0.12$\pm$	0.01	0.48$\pm$	0.04	0.33$\pm$	0.03
$-1.125$	0.30$\pm$	0.02	1.21$\pm$	0.07	0.64$\pm$	0.03
$-0.875$	0.44$\pm$	0.02	1.77$\pm$	0.07	1.08$\pm$	0.04
$-0.625$	0.26$\pm$	0.02	1.03$\pm$	0.07	1.34$\pm$	0.04
$-0.375$	0.17$\pm$	0.01	0.70$\pm$	0.05	1.51$\pm$	0.04
$-0.125$	0.11$\pm$	0.01	0.44$\pm$	0.03	1.62$\pm$	0.04
$0.125$	0.10$\pm$	0.01	0.40$\pm$	0.03	1.72$\pm$	0.04
$0.375$	0.13$\pm$	0.01	0.51$\pm$	0.04	1.85$\pm$	0.04
$0.625$	0.15$\pm$	0.01	0.61$\pm$	0.05	2.00$\pm$	0.05
$0.875$	0.16$\pm$	0.02	0.63$\pm$	0.07	2.16$\pm$	0.05
$1.125$	0.22$\pm$	0.01	0.89$\pm$	0.05	2.38$\pm$	0.05
$1.375$	0.22$\pm$	0.01	0.90$\pm$	0.04	2.61$\pm$	0.05
$1.625$	0.30$\pm$	0.01	1.19$\pm$	0.05	2.90$\pm$	0.05
$1.875$	0.30$\pm$	0.01	1.20$\pm$	0.06	3.20$\pm$	0.06
$2.125$	0.31$\pm$	0.01	1.25$\pm$	0.06	3.51$\pm$	0.06
$2.375$	0.27$\pm$	0.02	1.09$\pm$	0.07	3.78$\pm$	0.06

Equations

$$\lambda (O|\theta )$$

$$={\nu }_{\ast }(X,Y,Z|\theta )D^2\cos \delta$$

$${\rho }_{\text{mbox}CMD}(M_J,[J-K_s{]}_0|X,Y,Z)S(J,J-K_s,\alpha ,\delta ),$$

$$\ln \mathcal{L}(\theta )$$

$$=\sum _i\ln \lambda (O_i|\theta )-\int \mathrm{d}O\lambda (O|\theta ),$$

$$=\sum _i\ln {\nu }_{\ast }(X_i,Y_i,Z_i|\theta )-\int \mathrm{d}D{D}^2\mathrm{d}\alpha \mathrm{d}\delta \cos \delta {\nu }_{\ast }(X,Y,Z|\theta )$$

$$\int \mathrm{d}J\mathrm{d}(J-K_s){\rho }_{\text{mbox}CMD}(M_J,[J-K_s{]}_0|X,Y,Z)S(J,J-K_s,\alpha ,\delta ),$$

$$\text{textswab}S(\alpha ,\delta ,D)$$

$$\equiv \int \mathrm{d}J\mathrm{d}(J-K_s){\rho }_{\text{mbox}CMD}(M_J,[J-K_s{]}_0|X,Y,Z)S(J,J-K_s,\alpha ,\delta ),$$

$$\ln \mathcal{L}(\theta )=$$

$$-\int \mathrm{d}D{D}^2\mathrm{d}\alpha \mathrm{d}\delta \cos \delta {\nu }_{\ast }(X,Y,Z|\theta )\text{textswab}S(\alpha ,\delta ,D).$$

$${\nu }_{\ast }(X,Y,Z|\theta )=\sum _k{n}_k{\Pi }_k(X,Y,Z),$$

$$\ln \mathcal{L}(\{n_k{\}}_k)$$

$$=\sum _i\ln \sum _k{n}_k{\Pi }_k(X_i,Y_i,Z_i)$$

$$-\int \mathrm{d}D{D}^2\mathrm{d}\alpha \mathrm{d}\delta \cos \delta \sum _k{n}_k{\Pi }_k(X,Y,Z)\text{textswab}S(\alpha ,\delta ,D),$$

$$\ln \mathcal{L}(\{n_k{\}}_k)$$

$$\displaystyle\ \ =\sum_{k}\Big{[}N_{k}\ln n_{k}$$

$$\displaystyle\qquad-n_{k}\int\mathrm{d}D\,D^{2}\mathrm{d}\alpha\mathrm{d}\delta\,\cos\delta\,\Pi_{k}(X,Y,Z)\,\textswab{S}(\alpha,\delta,D)\Big{]}\,,$$

$${\hat{n}}_k=\frac{N_k}{\int \mathrm{d}D{D}^2\mathrm{d}\alpha \mathrm{d}\delta \cos \delta {\Pi }_k(X,Y,Z)\text{textswab}S(\alpha ,\delta ,D)}.$$

$${\hat{n}}_k=\frac{N_k}{{\int }_{{\Pi }_k}{\mathrm{d}}^{3}x\text{textswab}S(\alpha ,\delta ,D)}.$$

$$\Xi ({\Pi }_k)=\frac{{\int }_{{\Pi }_k}{\mathrm{d}}^{3}x\text{textswab}S(\alpha ,\delta ,D)}{{\int }_{{\Pi }_k}{\mathrm{d}}^{3}x},$$

$${\hat{n}}_k=\frac{N_k}{\Xi ({\Pi }_k)V({\Pi }_k)}=\frac{1}{\Xi ({\Pi }_k)}\frac{N_k}{V({\Pi }_k)}.$$

$${\sigma }_{{\hat{n}}_k}=\frac{{\hat{n}}_k}{\sqrt{N_k}}.$$

$$\text{textswab}S(\alpha ,\delta ,D)$$

$$\equiv \int \mathrm{d}J\mathrm{d}(J-K_s){\rho }_{\text{mbox}CMD}(M_J,[J-K_s{]}_0|X,Y,Z)S(J,J-K_s,\alpha ,\delta ),$$

$$\approx \sum _{j}S(M_{J,j}+\mu +A_J,[J-K_s{]}_{0,j}+E(J-K_s),\alpha ,\delta ),$$

$$\text{textswab}S(\alpha ,\delta ,D)$$

$${\nu }_{\ast }(Z)=n\left[(1-\alpha )\mathrm{sech}{\left(\frac{Z+Z_{\odot }}{2z_d}\right)}^2+\alpha \mathrm{sech}{\left(\frac{Z+Z_{\odot }}{2z_{d,2}}\right)}^2\right],$$

$${\log }_{10}\tau (M)=1.015-3.491{\log }_{10}M+0.8157({\log }_{10}M{)}^2.$$

$$\frac{\mathrm{d}n}{\mathrm{d}M}=0.016{\left(\frac{M}{{\mathrm{M}}_{\odot }}\right)}^{-4.7}{\mathrm{M}}_{\odot }^{-1}{\mathrm{pc}}^{-3},M>1{\mathrm{M}}_{\odot }.$$

$${\rho }_{\ast }^{\mathrm{MS}}=0.040\pm 0.002{\mathrm{M}}_{\odot }{\mathrm{pc}}^{-3},$$

$$\frac{\mathrm{d}N}{\mathrm{d}M}=9.3{\left(\frac{M}{{\mathrm{M}}_{\odot }}\right)}^{-6.5}{\mathrm{M}}_{\odot }^{-1}\mathrm{pc}{}^{-2},M>1{\mathrm{M}}_{\odot }.$$

$${\Sigma }_{\mathrm{SFR}}(t)=7.2\pm 1.0\exp (-t/7\pm 1\mathrm{Gyr})M_{\odot }{\mathrm{pc}}^{-2}{\mathrm{Gyr}}^{-1},$$

$$\frac{{\Sigma }_{\ast }(\tau )}{{\Sigma }_{\mathrm{form}}}$$

$$\mathrm{with}x={\log }_{10}\left(\tau /\mathrm{Gyr}\right)(\mathrm{solar}\mathrm{metallicity}).$$

$$\frac{{\Sigma }_{\mathrm{WD}}(\tau )}{{\Sigma }_{\mathrm{form}}}$$

$$\mathrm{with}x={\log }_{10}\left(\tau /\mathrm{Gyr}\right)(\mathrm{solar}\mathrm{metallicity}).$$

$$\frac{\mathrm{d}n}{\mathrm{d}{M}_J}=10^{-4.35+0.21M_J}\times \left[1+3.20\mathcal{N}(-0.92,0.27)\right],$$

Full text

††pagerange: Stellar Inventory of the Solar Neighborhood using Gaia DR1–LABEL:LastPage††pubyear: 2017

$\mathsection$ 1 Introduction

The Milky Way is a cornerstone in our understanding of the structure and evolution of galaxies. And within the Milky Way, the solar neighborhood provides a fundamental basis for studies of Galactic structure, star formation and stellar evolution, and Galactic dynamics. Within a few hundred parsec from the Sun, we can see stars that span the extremes from the most luminous main-sequence and giant stars to the faintest M dwarfs. This provides an essential basis for understanding the baryonic content of the Milky Way and external galaxies as it allows the initial-mass function (IMF; e.g., Gould et al. 1996; Kroupa 2001; Chabrier 2001), the mass-to-light ratio (e.g., Flynn et al., 2006), and the local star-formation history (e.g., Binney et al., 2000) to be determined directly. A complete baryonic census of the solar neighborhood is also important for comparing with dynamical determinations of the local mass distribution, which are only sensitive to the combined mass in baryons and dark matter (e.g., Holmberg & Flynn, 2000; Bovy & Rix, 2013).

The most precise censuses of stars in the solar neighborhood are largely based on small, volume-complete surveys of local stars (for example, within 25$\,\mathrm{pc}$, Reid et al. 2002; within 50$\,\mathrm{pc}$, Jahreiss et al. 1998), which contain only a handful of the most luminous stars. While this does not matter for determining the IMF from the number density of long-lived stars, for the purpose of, e.g., measuring the local stellar density distribution and comparing it to dynamical estimates (McKee et al., 2015), small volumes centered on the Sun are potentially biased, because the Sun is likely offset from the mid-plane of the Galaxy by 15 to 25$\,\mathrm{pc}$ (e.g., Binney et al. 1997; Chen et al. 2001; Jurić et al. 2008). Extending the local stellar census to a few hundred parsecs would allow a much better determination of the bright end of the stellar luminosity function and of the present-day mass function.

The ESA Gaia mission has been designed to investigate the luminous and dark-matter distribution within the Milky Way to constrain its formation and evolution (Gaia Collaboration et al., 2016a). While Gaia will observe more than 1 billion stars, these constitute only about 1 % of all of the stars in the Milky Way and the volume completeness for different types of stars will be complex even in the final data release due to variations in stellar colors, interstellar extinction, stellar density, and the observation pattern in different parts of the sky. Gaia released its first data in the Fall of 2016 (Gaia Collaboration et al., 2016b), which consists of the primary Tycho-Gaia Astrometric Solution (TGAS) containing positions, parallaxes, and proper motions for a subset of the Tycho-2 catalog (H$\o$g et al., 2000) and the secondary data set with approximate positions for stars brighter than $G\approx 20.7$ in the wide Gaia photometric bandpass $G$ (Lindegren et al., 2016). Based on 14 months of data and a preliminary astrometric processing, the Gaia DR1 TGAS catalog contains 2,057,050 stars with parallaxes with typical uncertainties of $\approx 0.3\,\mathrm{mas}$. As such, TGAS is by far the largest catalog of trigonometric, high-precision parallaxes for a uniformly selected sample of stars. However, no estimate of the completeness of TGAS is provided in the data products. This makes it difficult to use the TGAS data to investigate the spatial and dynamical structure of the solar neighborhood and many of the first uses of the TGAS data have therefore focused on purely kinematic studies that do not require knowledge of the selection function (e.g., Allende Prieto et al., 2016; Hunt et al., 2016; Bovy, 2017; Helmi et al., 2017).

This paper has two main parts. The first part consists of an extensive discussion of the volume completeness of the TGAS catalog for different stellar types. This may be useful to other studies employing TGAS data to study various aspects of the Milky Way’s structure. It can also be straightforwardly extended to the Gaia DR2 data when they appear. The second part uses the volume completeness combined with star counts for different classes of stars to perform a new stellar census of the solar neighborhood, covering main-sequence stars from early A-type to late K-type dwarfs and giant stars from the subgiant branch to the upper red-giant branch. We also determine the vertical density profiles for these stellar types up to a maximum of 400$\,\mathrm{pc}$ above the Galactic mid-plane. The resulting census is consistent with previous work and constitutes the most precise determination of the mass distribution in the solar neighborhood within the mass range observed in TGAS.

The structure of this paper is as follows. We introduce our approach to determining stellar densities and their spatial dependence from incomplete surveys in $\mathsection$ 2. In particular, we introduce the effective volume completeness, which for a given small volume in the Galaxy and a given stellar type represents the fraction of stars observed by the survey. In $\mathsection$ 3 we apply this formalism to the Gaia DR1 catalog and determine the completeness of the TGAS catalog for different stellar types: spectral types along the main sequence and giants in various luminosity ranges. The results from $\mathsection$ 3 depend on a detailed understanding of the “raw” TGAS selection function: the fraction of true stars contained in the catalog as a function of sky position, color, apparent magnitude, etc. We determine this raw selection function directly from the TGAS catalog by comparing it to the 2MASS catalog (Skrutskie et al., 2006); this is discussed in detail in Appendix A.

In $\mathsection$ 4, we apply the formalism from $\mathsection$ 2 to measure the mid-plane densities and the vertical density profile of different stellar types along the main sequence and fit parametric density laws. From this we determine the luminosity function, present-day mass function, total mid-plane number and mass densities, and the star-formation history of the solar neighborhood. In $\mathsection$ 5, we similarly measure the mid-plane densities and the vertical density profile for different types of giant stars and determine the luminosity function along the giant branch and the total mid-plane number and mass density of giants. One of the parameters that we fit for each stellar type is the Sun’s offset from the mid-plane as defined by that stellar type and we discuss the resulting offsets in detail in $\mathsection$ 6. We discuss some additional aspects of our results in $\mathsection$ 7 and conclude in $\mathsection$ 8. The basic data that we use throughout this paper comes from the TGAS catalog (Gaia Collaboration et al., 2016b) matched to photometry from the 2MASS catalog using a $4^{\prime\prime}$ nearest-neighbor search. Various cuts on this basic data set for different applications are described in the text below. We express the volume completeness in rectangular Galactic coordinates $(X,Y,Z)$, centered on the Sun with $X$ toward the Galactic center, $Y$ in the direction of Galactic rotation, and $Z$ directed toward the North Galactic Pole. We do not make any use of Galactocentric coordinates or any sort of kinematics and our results therefore do not depend on any assumed value of the standard Galactic constants (e.g., the distance to the Galactic center or the Sun’s velocity).

$\mathsection$ 2 Stellar Number Density Laws from Incomplete Surveys

Our goal is to determine the intrinsic stellar number density distribution $\nu_{*}(X,Y,Z)$, where $(X,Y,Z)$ is a set of heliocentric cartesian coordinates, of different stellar populations from TGAS observations of their three-dimensional position ($\alpha,\delta,D$) [right ascension, declination, and distance]. If TGAS observed all stars in the Galaxy, we would simply count the stars in small volumes and divide this number by the volume to get the density in terms of stars$\,\,\mathrm{pc}^{-3}$. However, because TGAS only observes stars in a finite magnitude range and because even at peak completeness it does not contain all stars, the situation is not quite this simple.

The problem of determining stellar number density laws $\nu_{*}(X,Y,Z)$ of stellar populations in the Milky Way has a long history (e.g., Bok, 1937; Bahcall, 1986; Jurić et al., 2008) and has been discussed recently by Bovy et al. (2012b) and Bovy et al. (2016). The discussion there focused primarily on spectroscopic surveys and on the problem of fitting a parametric density law to star count data. The focus here is different, in that we aim to determine the density $\nu_{*}(X,Y,Z)$ non-parametrically in bins in $(X,Y,Z)$ from data that cover a substantial fraction of the sky. Much of the discussion, especially about the effect of dust extinction, in this section follows that in Bovy et al. (2016).

We discuss how we determine the TGAS completeness as a function of $(J,J-K_{s},\alpha,\delta)$ by comparing the number counts in TGAS to those in 2MASS in Appendix A. Thus, we have a function $S(J,J-K_{s},\alpha,\delta)$ that gives the fraction of stars at a given $(J,J-K_{s},\alpha,\delta)$ contained in the TGAS catalog. As discussed in Appendix A, we only determine the completeness in the 48 % of the sky with ‘good’ TGAS observations and within this region, the completeness is independent of $(\alpha,\delta)$; therefore essentially $S(J,J-K_{s},\alpha,\delta)=S(J,J-K_{s})\,\Theta([\alpha,\delta]\ \mbox{in\ `good'\ part\ of\ the\ sky})$, where $\Theta(a)$ is the function that is one when $a$ holds true and zero otherwise. However, for generality, we will write all expressions in terms of the full $S(J,J-K_{s},\alpha,\delta)$.

2.1 Generalities

We determine $\nu_{*}(X,Y,Z)$ from a data set of stars that is not complete in any (simple) geometric sense and thus we need to take the selection function $S(J,J-K_{s},\alpha,\delta)$ into account. The connection between the selection function $S(J,J-K_{s},\alpha,\delta)$ and the three-dimensional position $(X,Y,Z)$ is made through a color-magnitude density (CMD) $\rho_{\mbox{\tiny CMD}}(M_{J},[J-K_{s}]_{0}|X,Y,Z)$ that gives the distribution in (absolute magnitude, unreddened color), potentially a function of position, that allows us to relate the observed $(J,J-K_{s},\alpha,\delta)$ to $(X,Y,Z)$ through the distance and three-dimensional extinction map. While these could (and should) in principle be inferred simultaneously with the stellar density, we assume in what follows that the CMD $\rho_{\mbox{\tiny CMD}}(M_{J},[J-K_{s}]_{0}|X,Y,Z)$ and the three-dimensional extinction map $(A_{J},E(J-K_{s}))[X,Y,Z]$ are known a priori.

To determine $\nu_{*}(X,Y,Z)$ we use a likelihood approach that models the full rate function $\lambda(O|\theta)$ that gives the number of stars as a function of all observables $O$ of interest for a set of model parameters $\theta$. These observables $O$ are in this case $(\alpha,\delta,D,J,J-K_{s})$—we will use $(X,Y,Z)$ and $(\alpha,\delta,D)$ interchangeably because they are related by coordinate transformation, but will keep track of the Jacobian—and we can write

[TABLE]

where we have assumed that the model parameters only affect $\nu_{*}$. In this decomposition, the factor $D^{2}\,\cos\delta$ comes from the Jacobian of the transformation between $(\alpha,\delta,D)$ and $(X,Y,Z)$.

An observed set of stars indexed by $i$ is a draw from a Poisson process with rate function $\lambda(O|\theta)$ with the likelihood $\mathcal{L}(\theta)$ of the parameters $\theta$ describing the density law given by

[TABLE]

where in the second equation we have dropped terms that do not depend on $\theta$. As in Bovy et al. (2016) we simplify this expression by defining the effective selection function $\textswab{S}(\alpha,\delta,D)$ defined by

[TABLE]

where $5\log_{10}\left(D/10\,\mathrm{pc}\right)=J-M_{J}-A_{J}[X,Y,Z]$ and $[J-K_{s}]_{0}=J-K_{s}-E(J-K_{s})[X,Y,Z]$. The ln likelihood then becomes

[TABLE]

Unlike the selection function $S(J,J-K_{s},\alpha,\delta)$ which depends on the survey’s operations only (which parts of the sky were observed, for how long, …), the effective selection function $\textswab{S}(\alpha,\delta,D)$ is a function of both the survey operations and the stellar population under investigation. Its usefulness derives from the fact that it encapsulates all observational effects due to selection and dust obscuration and turns the inference problem into a purely geometric problem. It directly gives the fraction of stars in a given stellar population that are observed by the survey in a given direction and at a given distance. Like the survey selection function $S(J,J-K_{s},\alpha,\delta)$, $\textswab{S}(\alpha,\delta,D)$ takes values between zero (fully incomplete) and one (fully complete). Under the assumption that the CMD $\rho_{\mbox{\tiny CMD}}(M_{J},[J-K_{s}]_{0}|X,Y,Z)$ and the three-dimensional extinction map $(A_{J},E(J-K_{s}))[X,Y,Z]$ are known (or at least fixed in the analysis), the effective selection function can be computed once for a given (survey, stellar population) pair.

To determine the best-fit parameters $\hat{\theta}$ of a parameterized density law $\nu_{*}(X,Y,Z|\theta)$ one has to optimize the ln likelihood given above. This ln likelihood can be marginalized analytically over the overall amplitude of the density (the local normalization if you will); this is discussed in Bovy et al. (2016) and similar expressions would apply here if such a marginalization is desired.

2.2 Non-parametric binned density laws

Now suppose that one wants to determine the density $\nu_{*}(X,Y,Z)$ of a stellar population in a set of non-overlapping bins in $(X,Y,Z)$. The bins are given by a set $\{\Pi_{k}\}_{k}$ of rectangular functions that evaluate to unity within the domain of the bin and zero outside of it. The domain can have an arbitrary shape, but typically this would be an interval in each of $X$, $Y$, and $Z$ or perhaps in $R_{xy}$ and $Z$, where $R_{x}y=\sqrt{X^{2}+Y^{2}}$. We can then write the density as

[TABLE]

where $\theta\equiv\{n_{k}\}_{k}$ is a set of numbers that give the density in each bin and that therefore parameterizes the density law.

The ln likelihood then becomes

[TABLE]

which, because the $\{\Pi_{k}\}_{k}$ are a set of non-overlapping bins, can be simplified to

[TABLE]

where $N_{k}$ is the number of points $i$ in the observed set that fall within bin $k$. We can maximize this likelihood for each $n_{k}$ analytically and find best-fit $\hat{n}_{k}$

[TABLE]

The denominator in this expression is known as the effective volume. Using the same symbol $\Pi_{k}$ to denote the three-dimensional integration region and using $x=(X,Y,Z)$ and $(\alpha,\delta,D)$ interchangeably because they are related through coordinate transformation, this can be written as the following simple expression

[TABLE]

Thus, the effective volume corresponding to a given spatial region $\Pi_{k}$ is the spatial integral of the effective selection function over $\Pi_{k}$. We can then define the effective volume completeness $\Xi(\Pi_{k})$ of the spatial region $\Pi_{k}$ as

[TABLE]

where $\int_{\Pi_{k}}\mathrm{d}^{3}x=V(\Pi_{k})$ is simply the actual geometric volume of $\Pi_{k}$. Because $\textswab{S}(\alpha,\delta,D)$ is a function bounded by zero and one, the effective volume completeness $\Xi(\Pi_{k})$ also takes values between zero and one. In terms of $\Xi(\Pi_{k})$, Equation (9) becomes

[TABLE]

This expression makes sense, because for a complete sample $\textswab{S}(\alpha,\delta,D)=1$, such that $\Xi(\Pi_{k})=1$ and this expression simplifies to the number divided by the volume of the bin, the standard way to compute a number density.

From the second derivative of the ln likelihood, we find the uncertainty on the $\hat{n}_{k}$

[TABLE]

$\mathsection$ 3 The TGAS completeness for different stellar populations

From the discussion above, it is clear that for a given survey and a given stellar population, the completeness functions of interest are (i) the effective selection function $\textswab{S}(\alpha,\delta,D)$ and (ii) the effective volume completeness $\Xi(\Pi_{k})$. In this section, we compute these functions for different stellar types along the main sequence and along the giant branch using the TGAS selection function determined in Appendix A. These functions will give a sense of how TGAS samples the extended solar neighborhood for a given stellar type.

3.1 Definitions of stellar types

The near-infrared CMD for stars with well-determined parallaxes $\varpi$ in TGAS is displayed in Figure 1. “Well-determined parallaxes” for the purpose of this figure means (i) $\varpi/\sigma_{\varpi}>10$ if $M_{J}<0$, (ii) $\varpi/\sigma_{\varpi}>20$ for $M_{J}>5$, and (iii) $\varpi/\sigma_{\varpi}>20-2\,(M_{J}-5)$ for $0<M_{J}<5$, with $M_{J}$ computed based on the TGAS $\varpi$ without correcting for extinction. These cuts are chosen to have a relatively well-populated upper main sequence and giant branch, because simple cuts on $\varpi/\sigma_{\varpi}$ tend to select few intrinsically-bright, and thus typically distant, stars. We use stars contained within the darkly-shaded region as a sampling $(M_{J},[J-K_{s}]_{0})_{j}$ of the intrinsic CMD $\rho_{\mbox{\tiny CMD}}(M_{J},[J-K_{s}]_{0})$, which we assume to be independent of position. This is a reasonable assumption for the $D\lesssim 1\,\mathrm{kpc}$ probed by TGAS. The darkly-shaded region is defined by a shifted and stretched version of the mean dwarf stellar locus from Pecaut & Mamajek (2013) and of the giant locus described below.

Along the main sequence, we define stellar types everywhere in this paper using the mean dwarf stellar locus from Pecaut & Mamajek (2013)111Specifically, we use version 2016.08.21 downloaded from http://www.pas.rochester.edu/~emamajek/EEM_dwarf_UBVIJHK_colors_Teff.txt ., displayed as a gray curve in Figure 1. Thus, we define A dwarfs as those stars along the main sequence with $J-K_{s}$ colors between the color of A0V and A9V stars according to the stellar locus and similarly for F, G, and K stars. For M stars we only go as red as the color of M5V dwarfs ($[J-K_{s}]_{0}=0.89$). In $\mathsection$ 4 we further subdivide each stellar type into subtypes A0V, A1V, etc. For those cases, the $J-K_{s}$ boundaries are located halfway between the stellar subtype in question and the adjoining subtypes. For example, for A1V the boundary goes from halfway between A0V and A1V to halfway between A1V and A2V. The sampling $(M_{J},[J-K_{s}]_{0})_{j}$ for a given stellar type then consists of those stars in the correct color range that fall within the darkly-shaded region of Figure 1. We always limit this sample to 1,000 stars for computational reasons.

We proceed similarly to define different types of giants. As a fiducial giant locus, we use a solar metallicity ($Z=0.017$) PARSEC isochrone with an age of $10^{0.8}\,\mathrm{Gyr}$ (Bressan et al., 2012). For simplicity, we define different types of giants using a simple cut on $M_{J}$: (i) “subgiants” for $1<M_{J}<4$, (ii) “lower red-giant branch” (RGB) for $-0.5<M_{J}<1$, (iii) “red-clump” (RC) for $-1.5<M_{J}<-0.5$, and (iv) “upper RGB” for $-4<M_{J}<-1.5$. In $\mathsection$ 5 below, we will further subdivide these types into $\Delta M_{J}=0.25\,\mbox{mag}$ bins. Similar to the main-sequence stars above, a sampling $(M_{J},[J-K_{s}]_{0})_{j}$ for a given giant type then consists of those stars in the correct $M_{J}$ range that fall within the darkly-shaded region to the right of the dashed line in Figure 1. This selection of giants largely avoids the asymptotic giant branch.

3.2 Effective selection function for different stellar types in TGAS

With the definitions of different stellar types from the previous subsection, we compute the effective selection function $\textswab{S}(\alpha,\delta,D)$ for different stellar types. In particular, we use the sampling $(M_{J},[J-K_{s}]_{0})_{j}$ from $\rho_{\mbox{\tiny CMD}}(M_{J},[J-K_{s}]_{0})$ for each stellar type to approximate Equation (3) as a Monte Carlo integration

[TABLE]

where we substitute $(J,J-K_{s})=(M_{J,j}+\mu+A_{J},[J-K_{s}]_{0,j}+E(J-K_{s}))$ in $S(J,J-K_{s},\alpha,\delta)$. In this equation, $\mu=5\log_{10}\left(D/10\,\mathrm{pc}\right)$ is the distance modulus and we have suppressed the dependence of the extinction map $(A_{J},E(J-K_{s}))$ on $(\alpha,\delta,D)$. In the absence of extinction, the effective selection function becomes

[TABLE]

Because our TGAS selection function $S(J,J-K_{s},\alpha,\delta)$ does not depend on $(\alpha,\delta)$ within the region of the sky over which it is defined, this is a one-dimensional function of distance.

The effective selection function in the absence of extinction is shown for different stellar types in Figure 2. Along the main sequence, the overall trend is that later types of stars can be seen out to a smaller distance than earlier types of stars. This is because of two reasons: not only are later-type dwarfs intrinsically fainter, but the TGAS selection function also has a brighter faint-end cut-off for redder, late-type stars (see Figure 29 in the appendix). On the other hand, the TGAS selection function has a slightly higher plateau value at intermediate magnitudes for redder stars, which has the effect that the M-dwarf effective selection function is slightly higher than for earlier types at its peak. The lower distance limit results from the bright cut-off of the TGAS selection function stemming from the exclusion of the brightest ($G\lesssim 6$) stars in TGAS. For A, F, and G dwarfs this causes a hole around the location of the Sun for these stars extending up to about $100\,\mathrm{pc}$ for A stars.

The effective selection function for different types of giants in Figure 2 behaves similarly to that of the dwarfs. Intrinsically-brighter giants can be seen out to larger distances, but not much larger distances, because the intrinsically-brighter (in the near infrared) giants on the upper RGB are also redder, where the TGAS selection function cuts off at brighter magnitudes. For both early-type dwarfs and luminous giants the effective selection function becomes small around $2\,\mathrm{kpc}$. The bright cut-off of the TGAS selection function again causes a hole around the location of the Sun of $\approx 50\,\mathrm{pc}$ to $\approx 200\,\mathrm{pc}$ for giants depending on their luminosity.

The curves in Figure 2 are for the case of zero extinction. For the intrinsically-brighter stellar types that extend out to $\gtrsim 1\,\mathrm{kpc}$, extinction has important effects. Because the extinction varies as a function of three-dimensional position, the effective selection function computed using Equation (13) becomes a function of $(\alpha,\delta,D)$. We compute the extinction using the combined extinction model from Bovy et al. (2016), which merges the three-dimensional extinction maps from Marshall et al. (2006), Green et al. (2015), and Drimmel et al. (2003). These maps are combined in that order in the case of overlap to create a full-sky three-dimensional extinction map (see the appendix of Bovy et al. 2016 for full details). In what follows, we use this map at its full resolution, which varies between $\approx 4^{\prime}$ to $\approx 20^{\prime}$.

The effect of extinction on the zero-extinction curves in Figure 2 is displayed in Figure 3, focusing on F dwarfs. The effect of extinction is mainly to lower the distance out to which a stellar population can be observed. A secondary effect is that the peak of the effective selection function can become higher, because of the higher completeness to redder stars in TGAS, but this is a minor effect. The narrow dark-gray region contains $68\,\%$ of the part of the sky for which we have determined the TGAS selection function; the light-gray region contains 95 %. To give another sense of this, the colored lines are lines of sight (in $N_{\mathrm{side}}=32$ HEALPix pixels) at the 10 through 90th quantile of the distribution of the mean extinction to $400\,\mathrm{pc}$. Thus, the effect of extinction is overall small, except for the about 10 % of the sky with high extinction, located near the Galactic plane.

3.3 Effective volume completeness maps for TGAS

Similar to the effective selection function $\textswab{S}(\alpha,\delta,D)$, we can compute the effective volume completeness $\Xi(\Pi_{k})$ for different kinds of volumes $\Pi_{k}$ and for different types of stars. The effective volume completeness in cylindrical regions with height $50\,\mathrm{pc}$ centered on $(X,Y)=(0,0)$ with the cylinder’s axis parallel to the $Z$ axis of the heliocentric coordinate frame is shown in Figure 4 for different stellar types along the main sequence. The top panel considers a narrow cylinder with a radius of $R_{xy}=50\,\mathrm{pc}$, the bottom panel a broader cylinder with a radius of $R_{xy}=250\,\mathrm{pc}$. Because the effective selection function of giants in Figure 2 closely resembles that of earlier-type main-sequence stars, the effective volumes for giants are similar to those of A, F, and G dwarfs and are not shown in Figure 4.

Because of the bright cut-off of TGAS, A, F, and G dwarfs have a hole in their effective volume centered on the Sun (see above) for the narrow cylinder. This hole is absent for the broader cylinder, because it gets filled in by stars at $|Z|\approx 0$, but with $(X,Y)=\mathcal{O}(100,100)\,\mathrm{pc}$. Because we only determined the TGAS selection function over $48\,\%$ of the sky, the broader cylinder has a complicated effective volume as a function of $Z$ for A and F stars. The effective volume for K dwarfs is small beyond a few $100\,\mathrm{pc}$ and M dwarfs are only sampled substantially by TGAS out to about $50\,\mathrm{pc}$.

The dashed curves in Figure 4 represent the effective volume when taking the full three-dimensional distribution of interstellar extinction into account, while the solid curves assume that there is no extinction. It is clear that extinction only has a very minor, percent-level effect on the effective volume in vertical slices near the Sun.

To get a sense of how TGAS samples stars of different types over the full three-dimensional volume around the Sun, we compute the effective volume completeness $\Xi(\Pi_{k})$ for cubic volumes in $(X,Y,Z)$ with a volume of $(100\,\mathrm{pc})^{3}$. This effective volume completeness is displayed for different main-sequence stellar types in vertical slices in Figure 5 and in slices in $X$ in Figure 6. The same for different types along the giant branch is shown in Figures 7 and 8. These effective volume completenesses take the full three-dimensional extinction map into account.

These maps of the TGAS volume completeness give a direct sense of the volume probed by different stellar types, with early main-sequence stars sampling a substantial fraction out to $\approx 2\,\mathrm{kpc}$ out of the plane and $\approx 1\,\mathrm{kpc}$ in the plane and late-type dwarfs only extending out to $\approx 50\,\mathrm{pc}$. The intrinsically-bright giants ($M_{J}\approx-2.5$) extend about as far as A dwarfs.

It is clear from these completeness maps that the TGAS completeness is complex, because of the patchy nature of the ‘good’ TGAS footprint over which we determined the selection function and because of the three-dimensional distribution of dust. Even if the remaining 52 % of the sky were included in the selection function, the volume completeness would remain complex, because the raw TGAS completeness is lower and complicated in this 52 % of the sky. These completeness maps need to be taken into account in any investigation using TGAS for which the stellar density or the manner in which a stellar population samples the local volume matters. Tools to compute the effective completeness and the effective volume completeness for a given stellar type are available in the gaia_tools.select.tgasEffectiveSelect class in the gaia_tools package available at

https://github.com/jobovy/gaia_tools .

$\mathsection$ 4 Stellar density laws for main-sequence stars

We now combine the method for determining non-parametric, binned stellar densities from $\mathsection$ 2.2 with the determination of the effective volume completeness from $\mathsection$ 3 to determine the stellar density and its vertical dependence for different stellar types along the main sequence. In $\mathsection$ 4.1 we test the effective volume completeness by directly comparing simple models for the stellar density of various types of stars combined with the effective volume completeness with the observed number counts in TGAS. In $\mathsection$ 4.2 we determine the stellar density and its vertical dependence in narrow $Z$ bins and fit these with simple analytic profiles. $\mathsection$ 4.3 focuses on the measurement of the mid-plane stellar densities and its implications for the present-day mass function, the star-formation history, and the density of stellar remnants in the solar neighborhood. We postpone a discussion of the Sun’s offset from the mid-plane defined by different types of stars to $\mathsection$ 6, where we discuss the measured offset for both main-sequence stars and giants.

In this and the subsequent section, we count as stars those TGAS sources that fall within the lightly-shaded region shown in Figure 1, including the overlap with the darkly-shaded region. This light-shaded region is a shifted and stretched version of the mean dwarf stellar locus from Pecaut & Mamajek (2013) and of the giant locus described above, similar to the dark-shaded region discussed in $\mathsection$ 3.1. It is designed, by hand, to encompass the vast majority of stars that are plausibly part of the dwarf and giant sequences, without including outliers. Because there are only very few stars outside of the light-shaded area and all we use are the number counts of stars, none of our results are significantly affected by changes to this area. We estimate distances simply as inverse parallaxes and only consider stars out to distances of $444\,\mathrm{pc}=0.2/[0.45\,\mathrm{mas}]$. Because the parallax uncertainties in the part of the sky that we consider are typically $<0.45\,\mathrm{mas}$ (see Appendix A), these stars typically have relative parallax uncertainties better than 20 %, small enough that the distance estimate from the inverse parallax is not strongly biased (Bailer-Jones, 2015). We only count stars up to $|Z|=412.5\,\mathrm{pc}$ and to a maximum $\sqrt{X^{2}+Y^{2}}<250\,\mathrm{pc}$, so this distance cut only affects these very furthest bins for the most luminous stars.

Because the effects of three-dimensional extinction are small (see discussion in \mathsection$$\mathsection 3.2 and 3.3), yet computationally laborious to compute, we ignore extinction in this section. We have performed the analysis described in this section taking into account the three-dimensional dependence of extinction for the main spectral classes (all A dwarfs, all F dwarfs, etc.) and find differences of only 1 to $2\,\%$ in the stellar density distributions and inferred parameters of the vertical density profile compared to assuming no extinction. These are a factor of a few or more smaller than our statistical uncertainties.

4.1 Forward modeling of the observed star counts

Before discussing our determination of the underlying stellar density profiles of main sequence stars, we test how well the observed star counts for different stellar types are represented by the effective volume completeness that we compute based on the raw TGAS selection function from Appendix A combined with the model for the CMD of different stellar types from $\mathsection$ 3.

The observed number counts of late K dwarfs (K5V through K9V) in a $100\,\mathrm{pc}$ wide cylinder centered on the Sun in $25\,\mathrm{pc}$ wide bins in $Z$ are displayed in Figure 9. From the discussion above in $\mathsection$ 3.3, it is clear that late K dwarfs can only be observed at relatively high completeness out to $|Z|\lesssim 100\,\mathrm{pc}$. From our prior understanding of the vertical density profile of late K dwarfs, we know that the scale height of these typically old stars is $\gg 100\,\mathrm{pc}$ (e.g., Jurić et al. 2008) and the vertical density profile should thus be close to constant within the observed TGAS volume. Therefore, we expect the observed number counts to largely reflect the completeness of the survey rather than the underlying stellar density profile. If the underlying density is constant, then the counts in equal-volume bins are proportional to the effective volume completeness: $N_{k}\propto\Xi(\Pi_{k})$ (Equation 11). The blue curve in Figure 9 shows the vertical dependence of the effective volume completeness, re-scaled to the median of the three central points (to fix the proportionality constant). This simple ‘model’ matches the observed number counts well.

To compare the observed number counts to those predicted by the effective volume completeness for stars that are more luminous or younger than late K dwarfs, we need to model the underlying density profiles. Below, we fit simple $\mathrm{sech}^{2}$ profiles to the binned stellar density laws of different stellar types. Figure 10 compares the observed number counts for early F dwarfs (F0V to F4V) and for RC stars (discussed in more detailed in $\mathsection$ 5 below) with underlying stellar densities $\hat{n}_{k}$ fit as $\mathrm{sech}^{2}$ profiles and multiplied by the effective volume ($N_{k}=\hat{n}_{k}\times[\Xi(\Pi_{k})\,V(\Pi_{k})]$ from Equation 11). It is clear that these models fit the observed number counts well: the vertical dependence and the dips near $Z=50\,\mathrm{pc}$ are the same in the model and the observed number counts. The dips are caused by the geometry of the ‘good’ part of the sky.

We conclude that we can successfully model the observed number counts of stars in TGAS using the effective volume completeness determined in $\mathsection$ 3.

4.2 Binned stellar densities along the main sequence

The binned vertical stellar density of stellar subtypes from A0V to K9V are displayed in Figure 11. To account for the fact that luminous A dwarfs can be seen to much larger distances than faint K dwarfs, these vertical densities are determined in $25\,\mathrm{pc}$ wide ranges in $Z$ in a cylinder that extends out to $R_{xy}=250\,\mathrm{pc}$ for A dwarfs, $200\,\mathrm{pc}$ for F dwarfs, $150\,\mathrm{pc}$ for G dwarfs, and $100\,\mathrm{pc}$ for K dwarfs. We do not consider M dwarfs, because the volume over which TGAS has relatively high completeness to M dwarfs is very small. While we determine absolute stellar densities, we shift the profiles in the $y$ direction in this figure to better display the vertical dependence. For each stellar type, we only show bins with (a) more than four stars, (b) effective volume completeness larger than $3\,\%$ of the maximum effective volume completeness for the stellar type in question, and (c) effective volume completeness larger than 1 part in $10^{5}$, except for the $Z=0$ bin. These cuts are designed to weed out bins where our stellar densities are highly noisy.

It is clear that the vertical profiles in Figure 11 broaden significantly when going from the earliest types of A dwarfs to later G and K dwarfs. By G5V, TGAS does not extend far enough to see a substantial decrease in the stellar density with vertical distance from the Sun. By K5V, the volume over which TGAS has significant completeness has become so small that it is difficult to measure the density at all.

We fit each stellar type’s number density profile by minimizing the $\chi^{2}$ deviation between the observed $\hat{n}_{k}\pm\sigma_{\hat{n}_{k}}$ and a model composed of two $\mathrm{sech}^{2}$ profiles

[TABLE]

where $n$ is the mid-plane density, $z_{d}$ is the scale height of the main $\mathrm{sech}^{2}$ component, $Z_{\odot}$ is the Sun’s offset from the mid-plane, $\alpha$ is the fraction of the mid-plane density that is part of the second $\mathrm{sech}^{2}$, and $z_{d,2}$ is the scale height of the second $\mathrm{sech}^{2}$ component, constrained to be larger than $z_{d}$. This $z_{d,2}\gg 1\,\mathrm{kpc}$ in almost all cases, such that the second $\mathrm{sech}^{2}$ profile is essentially a constant density. The second component typically has a small amplitude $\alpha$, especially for the earlier types (for later types, it becomes degenerate with the first component, because these types’ densities are close to constant). Thus, the measured number density profiles are well-represented by $\mathrm{sech}^{2}$ profiles within $Z\lesssim 4\,z_{d}$ or $Z<400\,\mathrm{pc}$ (if $4\,z_{d}>400\,\mathrm{pc}$). The fit are shown as blue curves in Figure 11; the gray band displays the $68\,\%$ confidence region obtained by MCMC sampling of the likelihood defined as $\exp(-\chi^{2}/2)$.

The parameters for each stellar type together with their uncertainties are given in Table 1. We list the mean color, absolute magnitude in $V$ and $J$, and the mass of each type determined from the dwarf locus from Pecaut & Mamajek (2013). We also give the derived quantities $\mathrm{d}n/\mathrm{d}M_{V}$ (the luminosity function), $\mathrm{d}n/\mathrm{d}M$ (the mass function), and the mid-plane mass density $\rho(>M)$ of stars with masses down to that of the stellar type in question. For types later than K4V, we are not able to determine reliable parameters because of the lack of data (see Figure 11); for types later than G3V, we are unable to determine the scale height because the distance out to which these stars can be seen in TGAS is too small compared to their scale height.

The luminosity function $\mathrm{d}n/\mathrm{d}M_{V}$ is shown in Figure 12. We compare the luminosity function to that determined by Reid et al. (2002) from a volume complete sample within $25\,\mathrm{pc}$ from the Sun. The agreement between these two determinations is good; importantly, they are essentially based on non-overlapping volumes, because for G dwarfs and earlier types, TGAS only starts observing stars at distances greater than about $25\,\mathrm{pc}$. This good agreement is partly due to the fact that we find that the Sun is at the mid-plane of the density of A and F stars (the measured Solar offsets using different stellar types are displayed in Figure 19 below and discussed in more detail together with the equivalent measurements for giants in $\mathsection$ 6). If this were not the case, the $25\,\mathrm{pc}$ volume of Reid et al. (2002) would be offset from the mid-plane by half a scale height. The new determination based on TGAS data has much higher resolution in absolute magnitude than Reid et al. (2002), which has $\Delta M_{V}=1$ compared to typical $\Delta M_{V}\approx 0.15$ here (of course, Reid et al. 2002 extend down to $M_{V}=20$, while we are limited to $M_{V}\lesssim 7.5$).

The scale height for different stellar types is shown in Figure 13 as a function of their main-sequence lifetime $\tau$, computed as in Reid et al. (2002)

[TABLE]

Stars of a given stellar type will largely be younger than this age. The scale height rises gradually from $\approx 50\,\mathrm{pc}$ for A dwarfs to $\approx 100\,\mathrm{pc}$ for late F dwarfs and larger for older G dwarfs. The scale height of young populations of stars is less than half of that of the atomic and molecular hydrogen (McKee et al., 2015), but similar to that of young open clusters ($\approx 50$ to $75\,\mathrm{pc}$, e.g., Bonatto et al., 2006; Buckner & Froebrich, 2014; Joshi et al., 2016) and of OB stars ($\approx 50\,\mathrm{pc}$, Reed, 2000).

4.3 The mass function and star-formation history of the solar neighborhood

Using a relation for the mass as a function of stellar type along the main sequence, we can translate the luminosity function of Figure 12 into a mass function $\mathrm{d}n/\mathrm{d}M$. We use masses for stellar types derived from the dwarf mean stellar locus from Pecaut & Mamajek (2013). The resulting mass function is shown in Figure 14. This function gives the mass density at the mid-plane in units of numbers of stars per solar mass and cubic pc. This is the present-day mass function. For long-lived stars this should reflect the IMF, while for short-lived stars the present-day mass function results from the combination of the IMF, stellar evolution, and kinematic heating. Our latest-type stars for which we can measure reliable mid-plane densities (K2V through K4V) are long-lived enough to trace all mass formed over the history of the disk (e.g., Bressan et al., 2012) and thus reflect the IMF for these types. Therefore, we anchor IMF models to the observed $\mathrm{d}n/\mathrm{d}M$ of these populations. In particular, we employ the lognormal and exponential IMF models from Chabrier (2001) and the broken-power law IMF model from Kroupa (2001). The range spanned by these three IMFs when anchored to the K dwarfs (using the median proportionality constant) is shown as the gray band in Figure 14.

We fit the high-mass end of the mass function using a power-law and find

[TABLE]

This can be compared to the mass function obtained by Scalo (1986), who finds $\mathrm{d}n/\mathrm{d}M=0.019\,(M/\,\mathrm{M}_{\odot})^{-5.4}\,\,\mathrm{M}_{\odot}^{-1}\,\mathrm{pc}^{-3}$; these differ by less than $25\,\%$ within the relevant mass range. Combining the power-law high-mass mass function at $M>1\,\mathrm{M}_{\odot}$ with the IMF fit to the K dwarfs at $M\leq 1\,\mathrm{M}_{\odot}$, we can compute the total mid-plane density of main-sequence stars. This is

[TABLE]

where the uncertainty includes a contribution from the uncertainty in the raw TGAS selection function for K dwarfs. As shown in Table 1, only $0.01\,\mathrm{M}_{\odot}\,\mathrm{pc}^{-3}$ of this is directly determined from our stellar density measurements and the rest is an extrapolation down to $M=0.08\,\mathrm{M}_{\odot}$ (assumed here to be the smallest mass above which core-hydrogen-burning occurs). This is in good agreement with previous determinations, e.g., the measurement of $\rho^{\mathrm{MS}}_{*}=0.036\pm 0.004\,\mathrm{M}_{\odot}\,\mathrm{pc}^{-3}$ from McKee et al. (2015) based in large part on the measurements from Reid et al. (2002).

Above $1\,\mathrm{M}_{\odot}$, the difference between the present-day mass function and the IMF is $\approx 0.014\,\mathrm{M}_{\odot}\,\mathrm{pc}^{-3}$ between $1<M/\,\mathrm{M}_{\odot}<3$. If we extrapolate to $M=8\,\mathrm{M}_{\odot}$, this difference is $\approx 0.027\,\mathrm{M}_{\odot}\,\mathrm{pc}^{-3}$ between $1<M/\,\mathrm{M}_{\odot}<8$. Assuming that all of the stars in this initial mass range turn into white dwarfs, the corresponding mass density of white dwarfs would be $0.0035\,\mathrm{M}_{\odot}\,\mathrm{pc}^{-3}$ using the initial-mass to final-mass relation of Kalirai et al. (2008) (as we will see below, the mass contained in giant populations is negligible in this context). This is a lower limit, because this estimate does not take into account that the older stars have more extended vertical profiles. Accounting for a typical factor of two difference in $z_{d}$ between old and young populations, the predicted density of white dwarfs becomes $0.0065\,\mathrm{M}_{\odot}\,\mathrm{pc}^{-3}$, in good agreement with the value of $0.0056\pm 0.0010\,\mathrm{M}_{\odot}\,\mathrm{pc}^{-3}$ from McKee et al. (2015).

When corrected for the effect of vertical heating, we can translate the present-day mass function into a cumulative star-formation history. To do this, we integrate the observed number volume density corresponding to the dominant $\mathrm{sech}^{2}$ component over $Z$ and obtain the observed number surface density for all stellar types. Similar to the volume density above, we can represent the high-mass end as a power-law:

[TABLE]

We turn the observed number surface density into an estimate of the total stellar mass formed going backwards in time from the present day up to the main-sequence lifetime of the stellar type, by running an IMF through the surface density measurement of each stellar type and integrating over all masses. That is, the current surface density $\mathrm{d}N_{i}/\mathrm{d}M(\tau_{i})$ for a stellar type $i$ with main-sequence lifetime $\tau_{i}$ corresponds to a total stellar mass $\Sigma_{*}(\tau_{i})$ up to $\tau_{i}$ in the past if we anchor an IMF to $\mathrm{d}N_{i}/\mathrm{d}M(\tau_{i})$ and integrate over all masses. The (reverse) cumulative star-formation history of the solar neighborhood thus determined is displayed in Figure 15, where the uncertainties are due to the uncertainty in the IMF (marginalizing over the three IMF models that we consider). We see that about $2\,\mathrm{M}_{\odot}\,\mathrm{pc}^{-2}$ was formed in the last Gyr and that the increase towards higher ages is faster than linear, indicating that the star-formation rate has decreased in time.

We fit a model for an exponentially-declining (or increasing) star-formation rate to these measurements and find that the star-formation rate declines in time and is given by

[TABLE]

where $t$ is time measured from $10\,\mathrm{Gyr}$ ago. The total amount of stellar mass formed in the last 10 Gyr is $\Sigma_{\mathrm{form}}=38.5\pm 2.5\,\mathrm{M}_{\odot}\,\mathrm{pc}^{-2}$. The declining star-formation rate found here is consistent with the analysis of Aumer & Binney (2009), who find the star-formation rate to be $2--7$ times lower now than at the start of star formation in the disk.

From this measurement of the total mass formed, we can estimate the mass that is presently still contained in visible stars. Based on maximum masses as a function of age from PARSEC (Bressan et al., 2012) isochrones and a Kroupa (2001) IMF, we find that the ratio of current-to-formed stellar mass is

[TABLE]

At approximately half-solar metallicity, this fraction is only about a percent lower; for other IMFs (e.g., those from Chabrier 2001), this fraction can be different up to $\approx 5\,\%$. For the exponentially-declining star-formation history, this gives $\Sigma_{*}=23.0\pm 1.5\,\mathrm{M}_{\odot}\,\mathrm{pc}^{-2}$ or about $60\,\%$ of the formed mass; this agrees, largely by construction, with directly counting all of the mass above $1\,\mathrm{M}_{\odot}$ using Equation (19) and combining it with mass below $1\,\mathrm{M}_{\odot}$ obtained in a similar manner as for $\rho^{\mathrm{MS}}_{*}$ above: $\Sigma_{*}=22.8\pm 1.0\,\mathrm{M}_{\odot}\,\mathrm{pc}^{-2}$. For a flat star-formation history, the ratio would be $\approx 63\,\%$; for a steeply-declining star-formation history with $e$-folding time of $1\,\mathrm{Gyr}$, the ratio would be $55\,\%$. Considering that we only include stellar mass contained in simple $\mathrm{sech}^{2}$ components within $400\,\mathrm{pc}$ from the mid-plane (what would traditionally be called the “thin” disk), this estimate is in good agreement with previous determinations of the surface-density of the thinner component of the disk (Flynn et al., 2006; Bovy et al., 2012a; McKee et al., 2015).

Similarly, we determine the ratio of the current mass in white dwarfs versus the total stellar mass formed using the initial–final mass relation from Kalirai et al. (2008) and find that

[TABLE]

At lower metallicity, this ratio is about $0.5\,\%$ higher. For the exponentially-declining star-formation history, this gives $\Sigma_{\mathrm{WD}}=3.6\pm 0.5\,\mathrm{M}_{\odot}\,\mathrm{pc}^{-2}$ or about $10\,\%$ of the formed mass (and about $16\,\%$ of the current stellar mass). For a flat star-formation history, the ratio would be $\approx 8.4\,\%$ (and about $13\,\%$ of the current stellar mass); for a star-formation history with $e$-folding time of $1\,\mathrm{Gyr}$, the ratio would be $11.5\,\%$ (and about $21\,\%$ of the current stellar mass). Accounting for the about $7\,\mathrm{M}_{\odot}\,\mathrm{pc}^{-2}$ currently in the thicker component of the disk (Bovy et al., 2012a) that we have ignored and that likely has a star-formation history more sharply peaked in the past, we estimate a total disk $\Sigma_{\mathrm{WD}}=5\pm 1\,\mathrm{M}_{\odot}\,\mathrm{pc}^{-2}$, again in good agreement with McKee et al. (2015)’s $\Sigma_{\mathrm{WD}}=4.9\pm 0.8\,\mathrm{M}_{\odot}\,\mathrm{pc}^{-2}$.

All in all, the mass distribution of the solar neighborhood and more broadly the solar cylinder determined from the TGAS star counts is in excellent agreement with previous studies and in many instances substantially more precise.

$\mathsection$ 5 Stellar density laws for giant stars

In the previous section, we discussed the stellar density profiles of different stellar types along the main sequence in detail. We do the same in this section for stars along the giant branch. Similar to the case of the main sequence above, we count as stars those objects that fall within the lightly-shaded region indicated in Figure 1, including the darkly-shaded region. We again ignore dust extinction, because its effect on the effective volume completeness is negligible within the volume that we consider for the purpose of determining vertical density profiles.

5.1 Binned stellar densities along the giant branch

Figure 16 shows the vertical stellar density profiles for giants in $\Delta M_{J}=0.25$ bins, determined from the number counts in a cylinder with a radius of $250\,\mathrm{pc}$ in $\Delta Z=25\,\mathrm{pc}$ bins. Comparing these to each other and to the density profiles of the later stellar types along the main sequence in Figure 11, we find that these profiles are all very similar. The giant profiles are significantly more noisy than the dwarf profiles, because of the relative paucity of giants in the TGAS sample (each $\Delta M_{J}=0.25$ bin typically only has about 1,000 stars). The measured profiles are well represented by single $\mathrm{sech}^{2}$ profiles: we find that the amplitude $\alpha$ of a second component when included in the fit is always so small that it is practically insignificant.

To obtain density profiles with smaller uncertainties, we also determine the vertical density of all stars belong to each major type of giant, because from Figure 16 it appears that they all have the same vertical profile. This is shown for the subgiants, lower RGB, and RC in Figure 17; there are too few stars on the upper RGB to obtain a high precision measurement of the vertical profile. These high-precision profiles demonstrate that the $\mathrm{sech}^{2}$ fit is indeed a good representation of the giants’ density profiles: the profiles clearly display the exponential decline at $|Z|>z_{d}$ and the flattening near $Z=0$. The scale heights of these $\mathrm{sech}^{2}$ fits are $\approx 150\,\mathrm{pc}$ (subgiants: $130\pm 20\,\mathrm{pc}$, lower RGB: $100\pm 33\,\mathrm{pc}$, RC: $154\pm 10\,\mathrm{pc}$).

5.2 The luminosity function of giants

From the $\mathrm{sech}^{2}$ fits to the density profiles in Figure 16 we determine the mid-plane number densities of giants. The luminosity function obtained from these is displayed in Figure 18 and tabulated in Table 2. This luminosity function has an overall exponential decline toward the upper RGB, with a bump at the location of the RC. The luminosity function is well fit by an exponential function of $M_{J}$ (or, equivalently, a power-law of luminosity), with a Gaussian bump superimposed

[TABLE]

in number $\mathrm{pc}^{-3}\,\mathrm{mag}^{-1}$ and where $\mathcal{N}(\mu,\sigma)$ is the (normalized) Gaussian probability distribution with mean $\mu$ and standard deviation $\sigma$.

The orange curve in Figure 18 shows the predicted luminosity function from the PARSEC (Bressan et al., 2012) isochrones weighted using a lognormal Chabrier (2001) IMF, a solar-neighborhood-like metallicity distribution function from Casagrande et al. (2011), and a uniform age distribution, using the same cuts to define the giant branch as used for the data. The amplitude of this predicted luminosity function is set such that the total stellar mass represented by $M>0.72\,\mathrm{M}_{\odot}$ stars in the model is $0.01\,\mathrm{M}_{\odot}\,\mathrm{pc}^{-3}$, the total mid-plane density that we directly measured from number counts along the main sequence in $\mathsection$ 4.3 above. The amplitude is thus not fit to the observed luminosity function of giants. It is clear that the agreement between the predicted and the observed luminosity functions for giants is excellent. The overall amplitude, overall decline toward more luminous giants, and the RC bump are all in good agreement.

The total mid-plane density of giants from integrating over the giant luminosity function between $-4<M_{J}<2.5$ is

[TABLE]

Jahreiss et al. (1997) find $n=0.00049\pm 0.00009\,\mathrm{giants}\,\mathrm{pc}^{-3}$ from a complete sample of giants within $25\,\mathrm{pc}$. This agrees well with the more precise value found here.

To turn this measurement of the mid-plane number density of giants into an estimate of the mid-plane mass density of giants requires the average stellar mass of giants of different luminosities. We determine the stellar mass of giants of a given luminosity $M_{J}$ using PARSEC (Bressan et al., 2012) isochrones, selecting giants as stars with surface gravities $\log g<3.75$ and $J-K_{s}>0.4$ and marginalizing over a flat age distribution and the solar-neighborhood-like metallicity distribution function from Casagrande et al. (2011). The mass function thus obtained is well fit by

[TABLE]

in $\,\mathrm{M}_{\odot}\,\mathrm{pc}^{-3}\,\mathrm{mag}^{-1}$. The total mid-plane mass density in giants is

[TABLE]

where the uncertainty is dominated by an (estimated) uncertainty of stellar mass along the giant branch. This is close to the value of $\rho^{\mathrm{giants}}_{*}=0.00060\pm 0.00012\,\mathrm{M}_{\odot}\,\mathrm{pc}^{-2}$ from Flynn et al. (2006).

$\mathsection$ 6 The Sun’s height above the mid-plane

The value of the Sun’s offset from the mid-plane defined by each spectral subtype from A0V to G3V is displayed in Figure 19. For later-type dwarfs we are unable to determine $Z_{\odot}$, because the stellar density for these stars is almost entirely flat within the observed $Z$ range (see Figure 11). Remarkably, the Sun is consistent with being at the mid-plane defined by each spectral type, with a combined measurement of $Z_{\odot}=-0.9\pm 0.9\,\mathrm{pc}$. For A and F dwarfs, $Z_{\odot}$ is determined from the vertical distribution over $\approx 4$ scale heights in both directions away from the mid-plane, leading to the high precision in the measured $Z_{\odot}$.

We measure the Sun’s offset from the mid-plane defined by the three major, different types of giants (subgiants, lower RGB, and RC) using the density profiles shown in Figure 17. The measured $Z_{\odot}$ for these three populations are: $Z_{\odot}=42\pm 7\,\mathrm{pc}$ (subgiants), $Z_{\odot}>45\,\mathrm{pc}\,(1\sigma)$ (lower RGB), and $Z_{\odot}=22\pm 6\,\mathrm{pc}$ (RC). For the lower RGB we only measure a lower limit with an a priori upper limit of $Z_{\odot}<50\,\mathrm{pc}$. This inferred value is strongly affected by the density profile at $|Z|>200\,\mathrm{pc}$ for lower RGB stars, which is somewhat asymmetric. Only considering lower RGB stars below $200\,\mathrm{pc}$, we find $Z_{\odot}=25\pm 13\,\mathrm{pc}$. Combining this value with the measurements using subgiants and RC stars, we find $Z_{\odot}=29\pm 4\,\mathrm{pc}$ with respect to the mid-plane defined by giants. If we only include giants at $|Z|<200\,\mathrm{pc}$ for all three types, we find $Z_{\odot}=35\pm 7\,\mathrm{pc}$.

We have investigated whether the inferred $Z_{\odot}$ depends on the assumed extinction model. The basic fits shown in Figure 19 ignore extinction in the calculation of the effective volume completeness, because it has only a minor effect. We have also determined the vertical stellar profiles of entire spectral types—all stars of type AV, all stars of type FV, etc.—taking into account the three-dimensional extinction model discussed in $\mathsection$ 3.2. The inferred $Z_{\odot}$ only differ by $\approx 0.1\,\mathrm{pc}$ from those determined without taking extinction into account for A and F dwarfs and by only $1\,\mathrm{pc}$ for G dwarfs.

As discussed in Appendix A, our model for the raw TGAS selection function is not perfect at the bright end. This may affect the results derived in this section, especially for the intrinsically-bright A and F dwarfs and for red giants, many of which are located near the bright-end cut-off of TGAS. Therefore, we have re-determined the stellar densities, removing all stars with $J_{G}<8$ (where $J_{G}$ is defined in Appendix A), which is where our model for the raw selection function performs poorly. The resulting density profiles are very similar to those derived from the full data set, although they are much noisier for the earliest A dwarfs. The combined value of the Sun’s height above the mid-plane is $Z_{\odot}=-0.3\pm 1.1\,\mathrm{pc}$, consistent with the value determined from the full data set. The scale height and mid-plane densities are also the same. For giants, we measure $Z_{\odot}=36\pm 5\,\mathrm{pc}$ when removing the bright stars and the bright-end behavior of the TGAS selection can therefore not explain the difference that we find between the $Z_{\odot}$ determined from A/F dwarfs and from giants.

We also investigate whether the wide cylinders (radii of $250\,\mathrm{pc}$ and $200\,\mathrm{pc}$) used to determine the vertical density profiles of A and F dwarfs affect the Solar offset from the mid-plane defined by these populations. We determine the density profiles for all types of main-sequence stars for cylinders with width $150\,\mathrm{pc}$. We find density profiles and $\mathrm{sech}^{2}$ fits that are almost the same as those based on the wider cylinders for A and F dwarfs. In particular, the Sun’s height above the mid-plane defined by A and F dwarfs is $Z_{\odot}=-0.7\pm 1.3\,\mathrm{pc}$, fully consistent with the value determined from the wider cylinders.

Finally, we determine the Solar offset only considering density measurements at $|Z|>100\,\mathrm{pc}$. We do this because the effective volume completeness for A and F stars has a mildly-peaked structure centered on $Z=0$, due to the bright-end TGAS cut-off and the sky cut (see Figure 4). If the effective volume completeness is even slightly incorrect, this could induce an apparent peak in the stellar density profile, which is shallow near the plane. At $|Z|>100\,\mathrm{pc}$, the effective volume for A and F stars is almost constant with $Z$. In this case we find the Sun’s height as defined by A and F dwarfs to be $Z_{\odot}=-2.5\pm 1.2\,\mathrm{pc}$, consistent with the result from the full data set.

Therefore, the TGAS data strongly prefer a value of $Z_{\odot}$ when measured using A and F dwarfs that is consistent with zero within a few pc. However, while we have performed a number of tests, described in the previous paragraphs, of the impact of the TGAS selection function on the measurement of $Z_{\odot}$, our TGAS selection function is only an empirical model derived using a number of assumptions, including that the selection function does not vary with position on the sky (within the ‘good’ part of the sky). These assumptions are unlikely to be correct in detail and may affect the $Z_{\odot}$ measurement in ways that are not covered by the tests discussed above. The next Gaia data release will not depend on the Tycho-2 catalog and will likely have a simpler selection function than TGAS. This will allow for a more robust determination of $Z_{\odot}$ using various stellar types.

Our results for $Z_{\odot}$ with respect to A/F dwarfs and with respect to giants are at odds with each other. The measured $Z_{\odot}$ with respect to A/F dwarfs also appears to be at odds with previous determinations, which find a canonical value of $Z_{\odot}=15$ to $25\,\mathrm{pc}$ (e.g., Binney et al., 1997; Chen et al., 2001; Jurić et al., 2008) and agree with the Solar offset we determine using giants. However, these literature determinations are based on old stellar populations, rather than the young populations traced by A and F stars here. A value of $Z_{\odot}=25\,\mathrm{pc}$ corresponds to one half of a scale height for the A dwarfs and can be clearly ruled out by the TGAS star counts. For younger tracers (young open clusters and OB stars), $Z_{\odot}$ is typically found to be smaller, with $Z_{\odot}=6$ to $15\,\mathrm{pc}$ (e.g., Bonatto et al., 2006; Joshi, 2007; Joshi et al., 2016), closer to what we determine here. Moreover, measurements of the vertical distribution of atomic (Malhotra, 1995) and molecular gas (Nakanishi & Sofue, 2006) in the Milky Way find typical fluctuations in the location of their mid-plane of tens of pc. These are most likely inherited by newly-formed stars and could persist for a few dynamical times for the younger stars traced by A and F dwarfs.

$\mathsection$ 7 Discussion

7.1 Computational complexity

In this paper, we have determined stellar densities using the $N=\mathcal{O}(10^{6})$ stars in the TGAS catalog. The full Gaia data set will eventually consist of $\mathcal{O}(10^{9})$ stars and it is therefore urgent to ask how the method used in this paper scales to this much larger number of stars. When determining stellar densities using the formalism in $\mathsection$ 2.2, the main computational complexity arises from the calculation of the effective volume in the denominator in Equation (9). The numerator in that expression can be efficiently computed for all volume bins simultaneously (e.g., using a histogram function) using $\mathcal{O}(N)$ operations and is fast even for $\mathcal{O}(10^{9})$ stars. The calculation of the effective volume does not scale directly with $N$, because it does not involve the data, only the selection function. However, because future data releases will cover a larger region of the Galaxy, the numerical integration of the effective volume will extend over a larger area and take longer to compute.

In the limit of vanishing dust extinction, computing the effective volume is easy, because the distance integral can be factored out from the integral on the sky. This is the case as long as the raw selection function does not depend on sky position, except for a purely geometric sky cut, like in our model for the TGAS selection function. While the selection function will certainly depend on sky position in small parts of the sky with high stellar density (where Gaia cannot catalog all sources due to telemetry limitations), over most of the sky future Gaia selection functions will likely not depend on position on the sky. However, extinction cannot be avoided when determining stellar densities close to the mid-plane and this complicates the effective-volume integral. Figures 4, 5, 6, 7, and 8 take the full three-dimensional dependence of dust extinction into account and took altogether a few days to compute using dozens of processors. Because the effective volume does not directly scale with the number of data points, as long as the extinction model is kept fixed (or only a small number of alternative models are explored), this computational time will not increase much in future data releases.

7.2 Other uses of the effective completeness for TGAS

We have used the effective completeness and the effective volume completeness for TGAS from $\mathsection$ 3 only for measuring the vertical density profiles of different stellar types. However, the effective completeness should be taken into account in any use of the TGAS data that depends on how the underlying stellar density distribution is catalogued in TGAS and this, in principle, includes essentially all use of the TGAS parallaxes as distance indicators. Because Gaia measures parallaxes rather than distances, any distance inferred from the TGAS parallaxes requires a distance prior (Bailer-Jones, 2015). This prior should take the effective completeness of TGAS into account. For example, to get the best inferred distances using the method of Astraatmadja & Bailer-Jones (2016a), their priors should be multiplied by the effective selection function, which they currently assume is uniform down to a magnitude limit (and they fully ignore the bright cut-off in TGAS; Astraatmadja & Bailer-Jones 2016b).

Similarly, because the effective completeness of TGAS is a strong function of stellar type (see Figure 2), any representation of the CMD from TGAS data (e.g., in Gaia Collaboration et al., 2016b) has strong selection biases and densities within the color–magnitude plane do not reflect true, underlying densities, but rather the density of objects contained in TGAS. This has strong implications for any analysis that models the observed CMD and uses it to, for instance, improve distances inferred from the TGAS parallaxes without attempting to correct for TGAS’s selection biases (e.g., Leistedt & Hogg 2017, L. Anderson et al. 2017, in prep.). Such models of the CMD cannot be applied to data other than the TGAS data themselves, because other data will have a different effective completeness and thus have different biases in the color–magnitude plane.

More broadly, many analyses of the TGAS data that use more than a small volume can be affected by the selection function, even if they do not directly depend on the number counts or underlying stellar densities. For example, analyses of the kinematics—typically assumed to be free of selection biases—of different stellar populations based on TGAS may depend on the selection function if they average over large volumes and gradients across this volume are important. For example, suppose one examines the kinematics of stars along the main sequence for a volume with $\sqrt{X^{2}+Y^{2}}<500\,\mathrm{pc}$ and $-100\,\mathrm{pc}<Z<100\,\mathrm{pc}$. Then the most luminous main-sequence stars will sample this volume almost uniformly, while the faintest main-sequence stars will only cover the nearest few tens of pc and the measured kinematic properties will effectively cover widely different volumes. If Galactic gradients are important, the effective completeness needs to be taken into account when comparing different stellar populations. The completeness maps displayed in Figures 5, 6, 7, and 8 can act as a guide to determine how important selection effects may be.

7.3 Astrophysical implications

We have obtained a number of new constraints on the stellar-populations structure of the local Milky-Way disk from the analysis of the vertical densities of main-sequence and giant stars in TGAS. One of the most directly measured quantities is the distribution of stellar mass in the mid-plane along the main sequence covering $0.7\,\mathrm{M}_{\odot}\lesssim M\lesssim 2.2\,\mathrm{M}_{\odot}$. While the TGAS data represent only a small fraction of the full, final Gaia catalog, this range extends both to high enough $M$ such that there is little mass in stellar populations with higher $M$ and to low enough $M$ that we reach the part of the main sequence where stars of any age have not evolved into giants yet. Thus, assuming an IMF model for $M\lesssim 0.7\,\mathrm{M}_{\odot}$, we obtain a census of all stellar mass in the solar neighborhood, an important ingredient in the baryonic mass budget in the solar neighborhood (McKee et al., 2015). Because—beside filling in the $G\lesssim 6$ part of the magnitude distribution and the part of the sky that we have ignored in the selection function—improvements in the sampling of the Gaia catalog will primarily come from going to fainter magnitudes and thus larger distances, the census described in this paper is practically the definitive census of $M\gtrsim 1\,\mathrm{M}_{\odot}$ stars in the solar neighborhood. However, future data releases will allow the low-mass IMF to be determined from the Gaia data and thus improve the census by directly including the contribution of low-mass stars.

We have used the measurements of the intrinsic stellar number counts at $M>1\,\mathrm{M}_{\odot}$, after translating them to surface densities to account for vertical heating, to determine the star-formation history of the solar neighborhood. Our measurement in Figure 15 is a direct measurement of the reverse-cumulative star-formation history and thus presents an important new constraint on models for the evolution of the Milky Way disk. However, to make this measurement we had to assume that the slope of the high-mass IMF is constant over the history of the disk and that it is known. Our results are degenerate with the slope of the high-mass IMF. For example, if we assume that the high-mass IMF tracks the present-day mass function that we measure, the points in Figure 15 would scatter around a constant value and we would infer that stars only formed during the last few hundred Myr. Using the kinematics of stars along the main sequence in addition to their number densities can help break this degeneracy (Binney et al., 2000) and using the number densities along the giant branch would provide additional constraints (see Figure 18).

The star-formation history that we measure is corrected for vertical heating. However, stars may also heat and migrate in the radial direction (Sellwood & Binney, 2002) and the locally-measured star-formation history in that case is essentially that of annuli in the disk convolved with the migration history (e.g., Schönrich & Binney, 2009). If migration is extremely efficient and the disk is fully mixed, the locally-measured star-formation history would simply be the star-formation history of the entire disk. Conversely, if migration is inefficient, the locally-measured star-formation history reflects the local evolution. The truth most probably lies somewhere in the middle and it is thus difficult to directly relate the locally-measured star-formation history to the evolution of the solar circle.

7.4 Future work

Future Gaia data releases will allow the stellar census performed here to be extended in many ways. Star counts down to $G\approx 20.7$ will allow all stars down to $M=0.08\,\mathrm{M}_{\odot}$ to be mapped within $\approx 25\,\mathrm{pc}$ and in rapidly-increasing volumes for higher masses (e.g., $\approx 190\,\mathrm{pc}$ for $M=0.15\,\mathrm{M}_{\odot}$, $\approx 500\,\mathrm{pc}$ for $M=0.3\,\mathrm{M}_{\odot}$). Thus the completeness maps shown for TGAS in Figures 5, 6, 7, and 8 will cover much larger volumes of the Galaxy beyond the closest kpc. For all but the faintest stars, the vertical and radial density law within the Galaxy can be determined to high precision, which will provide an unprecedented view of the distribution of mass and light in a large disk galaxy.

To get the most out of the Gaia data, various improvements to the analysis presented in the paper will have to be made. First of all, because no all-sky, complete stellar catalog currently exists down to $G\approx 20.7$, we will not be able to determine the catalog completeness by comparing to an external, complete catalog, as we did for TGAS using the 2MASS catalog in Appendix A. Ultimately, the completeness of the Gaia catalog should be determined from the Gaia data and instrument themselves, by considering how likely it is that sources with a given set of properties (sky location, brightness, color, …) appear in the final catalog. From the way in which Gaia scans the sky it is expected that, but for the most crowded regions of the sky, the final Gaia catalog will be 100 % complete down to some faint limit. For making the most complete stellar census based on Gaia data, determining the completeness around the faint limit will be of great importance, as small magnitude differences at the faint end represent large numbers of stars and large volumes. For example, by determining the completeness of TGAS using 2MASS rather than the Tycho-2 catalog, we were able to extend the magnitude range over which the completeness is understood by $\approx 1.5\,\mbox{mag}$ or a factor of two in distance. For Gaia DR2, we can still use 2MASS down to $J\approx 15$ to determine the completeness of the catalog, an improvement of $\approx 3\,\mbox{mag}$, or a factor of four in distance, over TGAS.

A second necessary improvement concerns the methodology used to infer the underlying stellar density. We have presented the statistical framework for inferring the underlying density from an incomplete survey in $\mathsection$ 2. In the current application of this general methodology, we have assumed that only the underlying stellar density is unknown, while the survey selection function, the density in the CMD $\rho_{\mbox{\tiny CMD}}$, and the three-dimensional extinction map are presumed to be known. While the survey selection function will hopefully be well known enough from the data processing itself (see previous paragraph), the density in the CMD and the three-dimensional extinction map are far from well known. Determining them is among the other main scientific goals of the Gaia mission (Gaia Collaboration et al., 2016a). Because we have focused in this paper on the nearest few hundred pc and have used near-infrared photometry, the three-dimensional extinction has only a minor effect on the observed number counts. For the CMD, we were able to approximate the underlying density by focusing on the high-precision parallaxes.

Ultimately, to gain a full, empirical understanding of the stellar distribution, CMD, and three-dimensional extinction throughout the Milky Way, all three should be obtained through a single analysis using the methodology described in $\mathsection$ 2. The likelihood given in the first line of Equation (2.1) is general and could include parameters of flexible models for the CMD and the extinction map (we have only assumed that these are known in the second line of Equation 2.1). So far, studies of the CMD or the three-dimensional extinction map have assumed that the other ingredients are perfectly known. For example, studies of the CMD using TGAS data assume that the density distribution and extinction are known (e.g., Leistedt & Hogg 2017, L. Anderson et al. 2017, in prep.). Similarly, the sophisticated determination of the three-dimensional extinction by Green et al. (2015) assumes that the stellar locus (the CMD) and the density distribution of the Milky Way are known (Green et al., 2014). With more Gaia data, especially with the BP/RP photometry, it should be possible to constrain all three ingredients simultaneously.

$\mathsection$ 8 Conclusion

We have conducted the first detailed stellar inventory of the solar neighborhood using the Gaia DR1 TGAS data. To do this, we have performed a detailed analysis of the raw TGAS selection function and have successfully derived the TGAS completeness over 48 % of the sky with ‘good’ TGAS observations. Using this raw completeness, we have determined the effective completeness in distance and volume of different stellar populations in TGAS, taking into account their intrinsic distribution in color and absolute magnitude as well as the three-dimensional dependence of dust extinction. Maps of the completeness of TGAS for different stellar tracers are given in Figures 5, 6, 7, and 8 and these should provide a useful guide for many studies making use of TGAS data.

Using our determination of the completeness of TGAS, we have measured the intrinsic stellar density distribution for different stellar types along the main sequence and along the giant branch. This results in a detailed new inventory of the stellar mass distribution in the solar neighborhood given in Tables 1 and 2. We have determined the luminosity function along the main sequence for $7.25>M_{V}>1.34$ (Figure 12), the present-day mass function for stars with masses $M\gtrsim 0.72\,\mathrm{M}_{\odot}$ (Figure 14), and total mid-plane density in stars (Equation 18). From these, we have determined the implications for the mass in white dwarfs and we have measured the star-formation history of the solar neighborhood (Figure 15 and Equation 20). We have further determined the luminosity function of stars along the giant branch (Figure 18) and the total number and mass density of giants in the mid-plane (Equations 24 and 26).

We have also measured the vertical density profiles of different types of dwarfs and giants. For all stellar types, we clearly see that the vertical density profiles flatten at $|Z|\lesssim 100\,\mathrm{pc}$ and all profiles are well represented as $\mathrm{sech}^{2}$ profiles within a few scale heights ($Z\lesssim 4\,z_{d}$). The scale height of these profiles increases smoothly when going from the earliest A-type dwarfs ($z_{d}\approx 50\,\mathrm{pc}$) to late K-type dwarfs ($z_{d}\approx 150\,\mathrm{pc}$); giants have similar profiles as late dwarfs. Surprisingly, we find that the Sun is at the mid-plane defined by A-type and F-type dwarfs ($Z_{\odot}=-0.9\pm 0.9\,\mathrm{pc}$), in tension with previous measurements. However, we are unable to identify any systematic in our analysis that would produce this result. With respect to older stars on the main sequence and on the giant branch, we find that the Sun is offset from the mid-plane by $Z_{\odot}=29\pm 4\,\mathrm{pc}$, in good agreement with previous measurements.

The new stellar inventory made possible by Gaia DR1 is in good agreement with previous studies, but substantially more precise for the stars that we directly observe ($M\gtrsim 0.72\,\mathrm{M}_{\odot}$). The detailed determination of the completeness of TGAS opens up many avenues of investigation using TGAS data that depend on how TGAS samples the underlying stellar distribution. The methodology for determining the completeness can also be straightforwardly extended three magnitudes fainter for Gaia DR2 and the tools developed here will thus remain useful. As described in Appendix A and $\mathsection$ 3, we have made code available that allows the raw and effective selection function for TGAS to be evaluated. All of the code used to perform the analysis presented in this paper is available at

https://github.com/jobovy/tgas-completeness ,

and can serve as an example of how to use the selection function.

Acknowledgments I thank Dustin Lang for providing the TGAS-matched 2MASS data and Chris Flynn and Wilma Trick for helpful comments. JB received support from the Natural Sciences and Engineering Research Council of Canada. JB also received partial support from an Alfred P. Sloan Fellowship and from the Simons Foundation.

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. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation.

Some of the results in this paper have been derived using the HEALPix (Górski et al., 2005), astropy (Astropy Collaboration et al., 2013), and emcee (Foreman-Mackey et al., 2013) software packages.

Appendix A The completeness of TGAS in color–magnitude–sky-position

Gaia DR1 consists of two astrometric solutions (Gaia Collaboration et al., 2016b; Lindegren et al., 2016): the primary TGAS data set containing positions, parallaxes, and proper motions for a subset of the Tycho-2 catalog (H$\o$g et al., 2000) and the secondary data set with approximate positions for stars brighter than $G\approx 20.7$. The primary solution consists of 2,057,050 stars out of 2,539,913 stars in Tycho-2. The Tycho-2 catalog is 99 % complete down to $V\approx 11$ and its completeness drops quickly at fainter magnitudes. However, TGAS does not share this simple completeness limit, as many of the $\approx 20\,\%$ of missing stars are at much brighter magnitudes and the very brightest stars ($G\lesssim 6$) are missing, because the specialized observing mode that they require is not yet sufficiently calibrated to produce reliable results. Due to the inhomogeneity of the scanning law, the completeness also strongly varies over the sky. In this Appendix, we investigate the completeness of TGAS.

A.1 Overall completeness

Almost all of the Tycho-2 stars were considered in the primary astrometric solution and the main reasons that stars failed to be included in TGAS are: (i) they are too bright ($G\lesssim 6$) or (ii) the quality of their astrometric solution is too low (as evidenced by a parallax uncertainty larger than 1 mas or a position uncertainty larger than 20 mas; Lindegren et al. 2016). While the hope is that the astrometric uncertainties will eventually be dominated by photon noise and, therefore, apparent magnitude, the quality of the astrometry in Gaia DR1 is primarily set by observational limitations and systematic uncertainties: the lack of high numbers of observations in certain parts of the sky, limitations in the current attitude model, and the simplistic modeling of the point-spread function. The completeness is therefore also in large part a function of broad observational properties (such as the number of astrometric transits) and color, but to a lesser degree of apparent magnitude (above the faint limit of Tycho-2).

Because the data processing for the secondary solution, which contains sky positions and broad band $G$ magnitudes down to $G\approx 20.7$, is still in a preliminary state, we cannot assume that the secondary solution is a complete sample that we can use to assess the completeness of TGAS. We could use Tycho-2 as the reference catalog, but this would limit any use of the TGAS selection function to $V\lesssim 11$, the 99 % completeness limit of Tycho-2 (without modeling the completeness of Tycho-2 itself). As we will see below, TGAS is 50 % complete in at least half of the sky down to $V\approx 12$ and about 20 % complete down to $V\approx 12.5$, allowing us to extend TGAS coverage about a factor of two in distance and eight in volume beyond the nominal Tycho-2 limit.

We therefore use the Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006) Point Source Catalog, which is $>99\,\%$ complete down to $J=15.8$ and $K=14.3$ over almost the entire sky. As we will see below, this is at least two magnitudes fainter than the TGAS completeness limit even in the best parts of the sky and therefore suffices for our purposes. We select unique ($\texttt{use\_src}=^{\prime\prime}\negmedspace\negmedspace 1^{\prime\prime}$), reliable point sources in 2MASS that either have signal-to-noise ratio greater than 10 in $J$ ($\texttt{ph\_qual}=^{\prime\prime}\negmedspace\negmedspace A^{\prime\prime}$) or are brighter than the detector saturation limit in either the “Read 2 $-$ Read 1” or in the “Read 1” exposures and similarly in $K_{s}$.

We assess the overall completeness by comparing the number of point sources with $6<J<10$ and $0<J-K_{s}<0.8$ in TGAS and 2MASS. The number counts on the sky in the two catalogs are shown in Figure 20. Here and in what follows we bin the sky using HEALPix222See http://healpix.sourceforge.net . level $N_{\mathrm{side}}=2^{5}=32$. At this level, the sky is divided into 12,288 equal-area pixels that have an approximate size of $1.8^{\circ}$ and an approximate area of $3.36\deg^{2}$. This size was chosen as a compromise between large pixels for good Poisson statistics on number counts in TGAS and 2MASS and small pixels that resolve the small-scale structure in the Gaia DR1 scanning law. These pixels are small enough to clearly display the structure in the scanning law (see Figure 21 below). By comparing the 2MASS counts to the TGAS counts in this relatively bright magnitude range, it is clear that while overall the counts are similar, the TGAS counts have significant features that are absent in 2MASS.

We extract some of the main properties of the Gaia DR1 scanning law and data processing directly from the TGAS catalog. Figure 21 displays the mean number of astrometric transits, the spread in this number, the mean spread in the scan directions (the direction along which stars transit the CCDs, using the catalog entry scan_direction_strength_k4), and the mean parallax uncertainty as a function of position on the sky. The mean number of observations clearly shows the imprint of the scanning law, with very few transits near the ecliptic, many near the ecliptic poles, and narrow ridges of abundant transits in ellipsoidal regions around the ecliptic. The spread in the number of transits is small, except in the ellipsoidal regions, which in reality are narrower than our pixel size. Below, we will remove the small part of the sky with high spread in the number of observations, because our sky pixelization is inadequate there. Over the majority of the sky our pixelization captures the properties of the observations well.

The lower-left panel in Figure 21 displays a measure of how well distributed the scan directions are in different parts of the sky. Because Gaia essentially performs one-dimensional scans, the astrometric accuracy is higher when different transits scan through a given field at a large variety of angles. As expected, the spread in fields near the ecliptic is small, because these fields have had few transits and thus cannot have a large spread in scan directions. However, the regions near the ecliptic poles also have a relatively small spread in scan directions, even though they have had many observations. As is clear from the lower-right panel with the mean parallax uncertainty, the regions near the ecliptic have relatively large astrometric uncertainties because of this small spread in scan directions.

Comparing Figures 20 and 21 one can see that many of the artifacts in Figure 20 are aligned with features in the scanning law. In many cases, these features coincide with the narrow ridges of the ellipsoidal rings around the ecliptic, which have large numbers of transits and one would, thus, naively expect to have high completeness.

We compute the overall completeness as the ratio of the TGAS counts in the range $6<J<10$ and $0<J-K_{s}<0.8$ to those in 2MASS. These ranges were chosen to contain enough stars in TGAS to allow a precise measurement (that is, not too hobbled by Poisson noise) of the completeness as a function of position on the sky using the above pixelization. The overall completeness is shown as a function of the properties of the scanning law and of the mean $J$ and $J-K_{s}$ color in Figure 22. It is clear that a low number of transits is most strongly associated with low overall completeness (top-left panel). The overall completeness quickly drops for regions with less than 10 transits. The overall completeness is also low when the spread in the scan directions is small (top-right panel). As discussed above, some of the fields with a large number of transits, but also a large spread in the number of transits have low completeness.

Based on these considerations of the Gaia DR1 scanning law, we define a ‘good’, or well observed, part of the sky as those pixels that satisfy the following cuts:

$\bullet$ Mean number of AL observations $\geq 8.5$;

$\bullet$ Spread in number of AL observations $\leq 10$;

$\bullet$ Spread in scan directions $\leq 0.7$;

$\bullet$ $|\mathrm{ecliptic\ latitude}|\geq 20^{\circ}$.

The last cut serves only to remove about $11\,\%$ of the sky that would otherwise remain as small, isolated islands around the ecliptic. These cuts select $48\,\%$ of the sky. Of the remaining $52\,\%$, $7\,\%$ is removed because it has a low number of transits (but otherwise good observations), $5\,\%$ because it has a small spread in scan directions, and $28\,\%$ because it has both a low number of transits and a small spread in scan directions. About $1\,\%$ of the sky is excluded because it has a large spread in the number of transits (the ellipsoidal ridges in Figure 21).

The overall completeness in the ‘good’ part of the sky is displayed in Figure 23. It is clear that the overall completeness is largely isotropic and does not contain sharp features. The overall completeness is slightly lower near the Galactic plane (which runs through much of the light-green parts of this map), which is because the mean $J$ is somewhat fainter near the plane and TGAS is not complete down to $J=10$ at all colors (see below). The sky pixels selected by the above cuts are shown as blue dots in Figure 22. The cut to the ‘good’ part of the sky selects regions that by and large have the same color and apparent magnitude distribution as the rest of the sky. The mean parallax uncertainty in the ‘good’ part of the sky is typically $\sigma_{\varpi}<0.45\,\mathrm{mas}$.

A.2 Completeness as a function of color and magnitude

The overall completeness in the ‘good’ part of the sky selected in the previous subsection does not appear to have any significant residual dependence on position on the sky (see Figure 22). Therefore, we continue under the assumption that the completeness in the ‘good’ region of the sky does not depend on sky position. We can then use the large number of stars in this region of the sky to determine the dependence of the TGAS completeness on color and magnitude in detail.

Figure 24 displays number counts in TGAS and 2MASS in the ‘good’ part of the TGAS sky in the range $-0.05<J-K_{s}<1$ and $J<14$ as a function of $(J-K_{s},J)$. The main features of the distribution at brighter magnitudes are the same in both catalogs, but there is both a lack of stars in TGAS at the bright and faint end. At the faint end in particular, there is a steep cut-off in the TGAS number counts that is absent in 2MASS and therefore needs to be because of the TGAS selection. This cut-off has a strong dependence on color, the approximate shape of which is indicated by the dashed line. The overall number counts (summed over $J-K_{s}$ in the $-0.05<J-K_{s}<1$ range) are shown in Figure 27. These number counts demonstrate the same bright- and faint-end cut-offs. The faint-end cut-off appears broad in this representation because it is summed over color. It is also clear that the TGAS counts nowhere quite reach the 2MASS counts, they fall short at every magnitude.

The TGAS completeness computed as the number counts in TGAS divided by those in 2MASS are displayed in Figure 28. This clearly shows the sharp drop in the completeness at faint magnitudes. To remove the strong color dependence of the completeness cut-off, we adjust the $J$ magnitude to a new $J_{G}=J+(J-K_{s})^{2}+2.5\,(J-K_{s})$ that runs approximately parallel to the completeness cut-off. This line for $J_{G}=12$ is shown in Figure 28. The $J-K_{s}$ color dependence is essentially caused by our use of a near-infrared color and magnitude for a survey whose completeness is more appropriately a function of an optical magnitude. The $J\rightarrow J_{G}$ relation is close to the relation that maps $J\rightarrow G$ along the stellar locus.

The number counts in TGAS and 2MASS as a function of $J_{G}$ are shown in Figure 27; the TGAS number counts drop more steeply in $J_{G}$ than they do in $J$. The completeness (ratio of TGAS to 2MASS number counts) as a function of $(J-K_{s},J_{G})$ is displayed in Figure 28. It is clear that the definition of $J_{G}$ has succeeded in removing the color dependence of the faint-end cut-off. A slight color dependence in both the level of the completeness around $J_{G}=10$ and in the bright cut-off remains. We ignore the latter, because for star counts there are very few stars at these bright magnitudes that contribute to the stellar density. To deal with the former, we divide the color range in three equal-sized bins with $\Delta J-K_{s}=0.35$ in the range $-0.05<J-K_{s}<1$ and approximate the selection function as being a function of $J_{G}$ only. The completeness in these bins as a function of $J_{G}$ is shown in Figure 29 (for easier interpretability, we have translated $J_{G}$ to $J$ at the center of each bin). The gray curve is a smooth spline fit.

The final model for the TGAS selection function $S(J,J-K_{s},\alpha,\delta)$ is therefore a the function that is (i) zero for ($\alpha,\delta$) outside of the ‘good’ part of the sky and (ii) given by the smooth spline model for the color bin in which $J-K_{s}$ is located, evaluated at $J_{G}$. We do not model stars bluer than $J-K_{s}=-0.05$ (mainly O and B stars) or redder than $J-K_{s}=1$. This model for the raw TGAS selection function is available in the gaia_tools.select.tgasSelect class in the gaia_tools package available at

https://github.com/jobovy/gaia_tools .

A.3 Comparison with Tycho-2

We test our model for the TGAS selection function by determining the completeness of TGAS with respect to its parent catalog Tycho-2 at apparent magnitudes brighter than the 99 % completeness limit of Tycho-2 ($V\approx 11$). Similar to how we determine the TGAS selection function above by comparing TGAS number counts to number counts in 2MASS, we compute the number counts as a function of color $B_{T}-V_{T}$ and magnitude $V_{T}$ in the ‘good’ part of the sky (where our model for the selection function is applicable). These number counts for TGAS and Tycho-2 are displayed in Figure 30. The ratio of these number counts gives the completeness of TGAS with respect to Tycho-2 and this ratio is shown in the left panel of Figure 31.

Our model for the TGAS selection function is a function of $(J,J-K_{s})$ and to compare it to the Tycho-2 number counts we need to translate the model to $(V_{T},B_{T}-V_{T})$. We do this using the following color–color transformations

[TABLE]

Evaluating our model for the selection function as a function of $(V_{T},B_{T}-V_{T})$ using these relations, we obtain the model displayed in the right panel of Figure 31. Because we use deterministic color–color relations that do not account for the scatter in this transformation, our model evaluated as a function of $B_{T}-V_{T}$ is constant in three $B_{T}-V_{T}$ ranges (owing to the fact that this is true for the model as a function of $J-K_{s}$). Comparing the measured TGAS/Tycho-2 completeness to the model in Figure 31, we find that the overall agreement of our model with the measurement is good. The model captures the shape and amplitude of the completeness from the bright end to the faint end and from the blue end to the red end.

A more detailed comparison is shown in Figure 32. This figure compares the completeness as a function of $V_{T}$ of TGAS versus Tycho-2 derived from the TGAS and Tycho-2 number counts to the model in eight bins in $B_{T}-V_{T}$. The model agrees well with the data in almost all cases, except (a) at the bright end in almost all color bins and (b) in its overall amplitude around $B_{T}-V_{T}\approx 1.1$. These discrepancies are expected based on how we determined the selection function above. To produce a simple model that accurately describes the faint end of the completeness, we performed a simple transformation $(J,J-K_{s})\rightarrow J_{G}$ in which the faint end cut-off is independent of color. However, the completeness at the bright end does depend on color in $J_{G}$ and our model does not capture this dependence. This causes the discrepancy at the bright end of the completeness in $(V_{T},B_{T}-V_{T})$. Clearly, the completeness has a simpler color dependency in $B_{T}-V_{T}$ than in $J-K_{s}$, so if we had a complete optical survey we could determine a more accurate selection function. For the purpose of this paper, a small discrepancy at bright magnitudes does not influence our results much, because very few stars have such bright magnitudes.

The discrepancy near $B_{T}-V_{T}\approx 1.1$ is caused by our choice of binning in $J-K_{s}$. For ease of interfacing with the large 2MASS database, the $J-K_{s}$ color bins were chosen to be equal size and our binning does not perfectly capture the higher completeness at the red end. This shows up most prominently around $J-K_{s}\approx 0.6$ or $B_{T}-V_{T}\approx 1.1$, which falls in the middle color bin of the model, but in reality is closer to the higher completeness of the reddest color bin of the model. The discrepancy amounts to only a few percent underestimation of the completeness in a narrow color strip.

Appendix B Overall completeness of the low-quality portion of the sky

We present the overall completeness (see $\mathsection$ A.1) of the badly-observed part of the TGAS sky in Figure 33. Comparing to the overall completeness of the ‘good’ part of the sky in Figure 23, the completeness in this badly-observed part is significantly lower and has artifacts due to the scanning law. The overall completeness of this part of the sky is also shown in Figure 34 split into the four main reasons why a part of the sky is not included in the ‘good’ part.