Milky Way Tomography with K and M Dwarf Stars: the Vertical Structure of the Galactic Disk

TL;DR

This study uses SDSS data of K and M dwarf stars to map the Milky Way's disk structure, revealing vertical asymmetries, wave-like features, and variations in scale height related to stellar populations and Sun's position.

Contribution

It provides a detailed analysis of the Milky Way's vertical structure using a large stellar sample, highlighting asymmetries and population differences with implications for galactic modeling.

Findings

Detection of wave-like vertical asymmetries in star counts.

Variation in disk scale height depending on star sample and location.

Potential indication of stellar population differences across latitudes.

Abstract

We use the number density distributions of K and M dwarf stars with vertical height from the Galactic disk, determined using observations from the Sloan Digital Sky Survey, to probe the structure of the Milky Way disk across the survey's footprint. Using photometric parallax as a distance estimator we analyze a sample of several million disk stars in matching footprints above and below the Galactic plane, and we determine the location and extent of vertical asymmetries in the number counts in a variety of thin- and thick-disk subsamples in regions of some 200 square degrees within 2 kpc in vertical distance from the Galactic disk. These disk asymmetries present wave-like features as previously observed on other scales and at other distances from the Sun. We additionally explore the scale height of the disk and the implied offset of the Sun from the Galactic plane at different locations, noting that the scale height of the disk can differ significantly when measured using stars only above or only below the plane. Moreover, we compare the shape of the number density distribution in the north for different latitude ranges with a fixed range in longitude and find the shape to be sensitive to the selected latitude window. We explain why this may be indicative of a change in stellar populations in the latitude regions compared, possibly allowing access to the systematic metallicity difference between thin- and thick-disk populations through photometry.

Figures

Tables

Table 1. Table 1: Best-fit Results to the North only (n), South only (s), and North-South Combined (c) Selections of the Uniform Latitude Sample.

$l$ range (deg)		$N_0$ $(\times {10}^6)$	$H_1$ (kpc)	$H_2$ (kpc)	$f$	${\chi }^2/\mathrm{dof}$
70 75	s	6.11 $\pm$ 0.33	0.211 $\pm$ 0.007	0.688 $\pm$ 0.028	0.112 $\pm$ 0.008	1.18
	c	4.83 $\pm$ 0.17	0.239 $\pm$ 0.006	0.773 $\pm$ 0.034	0.108 $\pm$ 0.008	1.19
	n	3.86 $\pm$ 0.16	0.278 $\pm$ 0.009	1.004 $\pm$ 0.118	0.083 $\pm$ 0.013	0.96
75 85	s	6.01 $\pm$ 0.24	0.210 $\pm$ 0.005	0.662 $\pm$ 0.019	0.116 $\pm$ .006	1.07
	c	5.07 $\pm$ 0.13	0.234 $\pm$ 0.004	0.742 $\pm$ 0.022	0.104 $\pm$ 0.006	1.12
	n	4.35 $\pm$ 0.13	0.262 $\pm$ 0.006	0.899 $\pm$ 0.059	0.082 $\pm$ 0.008	0.95
85 95	s	5.71 $\pm$ 0.22	0.213 $\pm$ 0.005	0.662 $\pm$ 0.020	0.113 $\pm$ 0.007	1.23
	c	4.84 $\pm$ 0.12	0.239 $\pm$ 0.004	0.735 $\pm$ 0.024	0.100 $\pm$ 0.006	1.32
	n	4.19 $\pm$ 0.12	0.271 $\pm$ 0.006	0.933 $\pm$ 0.077	0.069 $\pm$ 0.009	1.14
95 105	s	5.76 $\pm$ 0.23	0.211 $\pm$ 0.005	0.654 $\pm$ 0.021	0.108 $\pm$ 0.007	0.83
	c	4.60 $\pm$ 0.11	0.242 $\pm$ 0.004	0.734 $\pm$ 0.026	0.098 $\pm$ 0.007	1.20
	n	3.86 $\pm$ 0.10	0.278 $\pm$ 0.006	0.948 $\pm$ 0.088	0.066 $\pm$ 0.009	1.18
105 120	s	6.17 $\pm$ 0.24	0.196 $\pm$ 0.005	0.595 $\pm$ 0.013	0.119 $\pm$ 0.005	1.10
	c	4.74 $\pm$ 0.10	0.235 $\pm$ 0.004	0.695 $\pm$ 0.019	0.099 $\pm$ 0.006	1.48
	n	4.01 $\pm$ 0.09	0.276 $\pm$ 0.005	1.028 $\pm$ 0.092	0.052 $\pm$ 0.006	1.13
120 135	s	5.69 $\pm$ 0.20	0.206 $\pm$ 0.005	0.626 $\pm$ 0.017	0.105 $\pm$ 0.006	1.09
	c	4.72 $\pm$ 0.10	0.236 $\pm$ 0.004	0.718 $\pm$ 0.022	0.085 $\pm$ 0.005	1.38
	n	4.18 $\pm$ 0.10	0.268 $\pm$ 0.005	0.952 $\pm$ 0.074	0.051 $\pm$ 0.006	1.10
135 150*	s*	0.67 $\pm$ 0.04	0.588 $\pm$ 0.015	0.193 $\pm$ 0.005	9.580 $\pm$ 0.471	1.20
	c	4.76 $\pm$ 0.10	0.232 $\pm$ 0.004	0.686 $\pm$ 0.020	0.088 $\pm$ 0.005	1.63
	n	4.00 $\pm$ 0.09	0.273 $\pm$ 0.005	1.006 $\pm$ 0.095	0.045 $\pm$ 0.006	1.21
150 165	s	6.20 $\pm$ 0.27	0.190 $\pm$ 0.005	0.574 $\pm$ 0.012	0.115 $\pm$ 0.005	1.12
	c	4.78 $\pm$ 0.11	0.229 $\pm$ 0.004	0.666 $\pm$ 0.019	0.092 $\pm$ 0.006	1.65
	n	4.14 $\pm$ 0.10	0.269 $\pm$ 0.004	1.015 $\pm$ 0.097	0.041 $\pm$ 0.005	1.15

Table 2. Table 2: Quantifying the Effects of Latitude Selection in the North-only Analysis

$l$ range (deg)	$b$ range (deg)	$⟨z_{\mathrm{raw}}⟩$ (kpc)	$z$ range (kpc)	$\Delta$[Fe/H] (dex)	$\Delta M_r$ (mmag)
35 40	30 50	1.57	0.35 0.5	0.01	-10
	50 90	2.09	0.8 1.3	-0.08	80
			1.7 2.0	-0.49	450
40 45	30 50	1.55	0.35 0.5	-0.04	40
	50 90	2.10	0.8 1.3	-0.1	100
			1.7 2.0	-0.6	500
45 50	30 50	1.54	0.35 0.5	-0.03	30
	50 90	2.08	0.8 1.3	-0.05	50
			1.7 2.0	-0.7	600
50 55	30 50	1.53	0.35 0.5	0.05	-50
	50 90	2.07	0.8 1.3	-0.06	60
			1.7 2.0	-0.53	480
55 60	30 50	1.51	0.35 0.5	0.05	-50
	50 90	2.07	0.8 1.3	-0.06	60
			1.7 2.0	-0.62	550
165 180	30 50	1.29	0.35 0.5	-0.3	280
	50 90	1.95	0.8 1.3	-0.48	440
			1.7 2.0	-1.32	1010
180 195	30 50	1.29	0.35 0.5	-0.44	410
	50 90	1.95	0.8 1.3	-0.46	420
			1.7 2.0	-1.4	1100
195 203	30 50	1.30	0.35 0.5	-0.5	500
	50 90	1.95	0.8 1.3	-0.34	320
			1.7 2.0	-1.26	977
203 210	30 50	1.28	0.35 0.5	-0.29	280
	50 90	1.92	0.8 1.3	-0.38	360
			1.7 2.0	-1.3	1000

Equations

$$M_r$$

$$d=d[r_0,M_r((r-i{)}_0,[Fe/H])]+0.1415(r-i{)}_0+0.0436,$$

$$d=d[r_0,M_r((g-i{)}_0,[Fe/H])]+0.08982(g-i{)}_0-0.0726.$$

$$\mathcal{V}(z)=\frac{1}{2}\delta (l_2-l_1)z^2\left(\frac{1}{{\sin }^2{b}_1}-\frac{1}{{\sin }^2{b}_2}\right),$$

$$x=-8+d\cos b\cos l;y=d\cos b\sin l;z=d\sin b$$

$$n(z)=n_0\left({\mathrm{sech}}^2\left(\frac{(z+z_{\odot })}{2H_1}\right)+f{\mathrm{sech}}^2\left(\frac{(z+z_{\odot })}{2H_2}\right)\right),$$

$${\chi }^2=\sum _{i=1}^{N_{\mathrm{bin}}}{\left(\frac{\left(N_i-n(z_i)\right)}{{\sigma }_i}\right)}^2,$$

$$\exp ((R-R_0)/(2.7\mathrm{kpc})),$$

$$n_{\mathrm{eff}}(z)=\sum _i{n}_{\mathrm{raw}}(z,R_i)\exp ((R_i-R_0)/(2.7\mathrm{kpc})),$$

$${\chi }_b^2=\sum _{i=1}^{N_{\mathrm{bin}}}\frac{{\left(n_i^{>}-n_i^{<}\right)}^2}{n_i^{<}},$$

Full text

galaxies: structure – galaxies: evolution

††slugcomment: FERMILAB-PUB-17-187-AE, UCI-TR-2017-05

1 Introduction

The study of the structure of the Milky Way galaxy emerges as a key tool in the analysis of the nature of its dark halo and the most massive dwarf galaxies within it, through the manner in which they interact with the visible gas, dust, and stars of the Galactic disk. Astrometric observations at large scales support a cosmology of cold dark matter (CDM) and dark energy, the so-called $\Lambda$CDM model (Planck Collaboration et al., 2016), and an “inside-out” formation history of the cosmos (White & Rees, 1978; Frenk et al., 1983) and indeed of the Milky Way. At small scales this model predicts many more satellite galaxies in the Milky Way than have been observed (Kauffmann et al., 1993; Klypin et al., 1999). Discoveries of faint Milky Way satellites continue to be made (see, e.g., Belokurov et al. (2007)), however, and the technical limitations of current observations explain at least some of the mismatch (Tollerud et al., 2008; Walsh et al., 2009; Bullock et al., 2010) in the numbers of observed and expected satellite galaxies, or subhalos. Nevertheless, questions persist concerning the number, evolution, and mass distribution of Milky Way subhalos (see Kravtsov (2010) for a review), as well as their stellar content (Boylan-Kolchin et al., 2011b). Yet more dwarf galaxies could potentially be detected through a search for their tidal imprint on outer gaseous disks (Chakrabarti et al., 2011), and low-mass subhalos can be destroyed through interactions with the disk (D’Onghia et al., 2010): these disjoint ideas argue for the importance of the observational study we effect here, namely, of the vertical structure of the Galactic disk within 2 kpc of the Sun in height and in-plane distance.

The hierarchical nature of structure formation, i.e., the ongoing merging of clumps with time, suggests that dark matter could have surviving phase-space structure, which could also be imprinted on the stars. This could occur, e.g., through the formation of stellar streams in the halo. In this paper we are interested in identifying possible traces of halo-disk interactions and thus wish to search for the appearance of stellar number count distributions that break the spatial symmetries expected from the integral of motions embedded in virialized probability distribution functions. To do this we focus on an observational study of that possibility through use of a sample of some 3.6 million K and M dwarfs from the Sloan Digital Sky Survey (SDSS).

The appearance of asymmetric structures in the Milky Way has been noted through HI gas, dust, and stellar tracers, particularly in association with the warping of the disk (Binney, 1992). Various scenarios have been proposed for the appearance of warps, though analytic (Nelson & Tremaine, 1995) and numerical (Shen & Sellwood, 2006) studies suggest a dynamical origin: warps can appear and disappear through interactions of the disk with the halo and/or the satellites it contains over time scales short with respect to the age of the universe. The stellar disk is more complicated and thus exhibits additional features. Rings (Newberg et al., 2002; Morganson et al., 2016), as well as ripples (Price-Whelan et al., 2015; Xu et al., 2015), have been noted, with the latter at distances in excess of 10 kpc from the Galactic center, where the disk is thinning out. In the solar neighborhood vertical asymmetries have been noted in the number counts (Widrow et al., 2012; Yanny & Gardner, 2013), as well as in the radial velocities (Widrow et al., 2012), which have been confirmed by other studies (Carlin et al., 2013; Williams et al., 2013). These asymmetries could be temporally recent as well and be generated by non-axisymmetric features in the disk, such as the spiral arms or galactic bar (Debattista, 2014; Faure et al., 2014; Monari et al., 2015, 2016), or by interactions of the disk with the halo and its embedded satellite galaxies (Gómez et al., 2013; Widrow et al., 2014; Gómez et al., 2016, 2017; Laporte et al., 2016). It is thought that the asymmetries in the number counts are more suggestive of a dynamical origin (Laporte et al., 2016), although, on the other hand, it appears that the spiral arms possess out-of-plane structure as well (Camargo et al., 2015c).

Only two of the Milky Way’s satellites are known to be possibly massive enough and close enough to be able to perturb the structure of the dark halo and Galactic disk: the Sagittarius dwarf spheroidal (dSph) galaxy at 20 kpc and the Large Magellanic Cloud at 50 kpc from the Sun (Weinberg, 1998; Jiang & Binney, 1999; Bailin, 2003; Purcell et al., 2011; Gómez et al., 2013). Such massive objects are believed to have formed at late cosmic times (Boylan-Kolchin et al., 2011a). Numerical simulations of the tidal interactions of the satellites reveal significant vertical perturbations of the disk (Gómez et al., 2013; Widrow et al., 2014; Gómez et al., 2016, 2017; Laporte et al., 2016). However, simulations of such effects do not yield the precise vertical asymmetries observed in the existing data (Laporte et al., 2016).

In this paper we revisit and expand the earlier studies of Widrow et al. (2012) and Yanny & Gardner (2013) to scrutinize the vertical structure of the Galactic disk and its variation across the Galactic plane. We use a larger stellar photometric sample to subdivide the disk of the Milky Way within $|z|<2$ kpc in vertical distance from the Galactic plane into bins of longitude $l$ and latitude $b$ to explore the extent to which the large-scale asymmetries persist to smaller scales. We show that these asymmetries are more pronounced toward the Galactic anticenter, and we find similarities to the ripples seen at further distances from the Sun by Xu et al. (2015) and Price-Whelan et al. (2015). We provide quantitative locations and amplitudes of these asymmetries so that those who model Milky Way structure can compare with them; associated with these features are changes in the disk scale heights, north and south — we report these as well. Interestingly, numerical simulations of the tidal interactions of the disk-satellite-halo system reveal an induced spiral arm and barred structure in the Galactic disk (Purcell et al., 2011; Laporte et al., 2016). This also motivates our observational study, and we compare the variations we do observe with the known spiral arm structure of the Galactic disk (Camargo et al., 2015c).

Our studies are also pertinent to the long-standing problem of the determination of the matter density in the vicinity of the Sun (Kapteyn, 1922), as inferred from the measured kinematics of the local stars (Oort, 1932). In this so-called Oort problem, an assumed gravitationally relaxed population of stars is used to trace the local gravitational potential and to infer the local matter and dark matter densities, where we note Binney & Tremaine (2008) and Read (2014) for reviews. The appearance of vertical oscillations in the stellar number counts and velocity distributions of the Milky Way (Widrow et al., 2012; Yanny & Gardner, 2013) suggests that the local stars may not be sufficiently gravitationally relaxed, incurring additional uncertainty in the assessment of the local dark matter density (Banik et al., 2017), a parameter key to current efforts to detect dark matter directly (Peter, 2010). We believe that further observational study and analysis, as we help realize in this paper, will be able to resolve the origin of the vertical asymmetries and ultimately that of their impact on the assessment of the local dark matter density. Further studies with Gaia (Gaia Collaboration, 2016a; Gaia Collaboration et al., 2016b) will also be key (Banik et al., 2017).

The SDSS (York et al., 2000) is a mature observational platform that gives us access to a very large photometric sample of stars, observed with the same telescope under carefully calibrated conditions. The uniformity of the photometric calibration is ensured through multiple observations of the same stars, rendering it nearly free of astrophysical assumptions (Padmanabhan et al., 2008). The SDSS is the first survey to provide significant coverage in the south, allowing us to study and compare regions in both the north and the south. We study red, main-sequence stars, particularly K and M dwarfs, because main-sequence stars vary little in their intrinsic luminosity for a given color and metallicity, and the population of faint red stars is less likely to be infiltrated by non-main-sequence stars, such as giants. Our largest analysis sample contains 3.6 million of such stars. By breaking the sky into regions, we study the vertical structure of the Galaxy and search for changes across the footprint. We use Galactic longitude and latitude with the Sun at the center such that the longitude is measured counterclockwise from the line toward the Galactic center and latitude is measured from the plane of the Galactic equator such that positive latitude is north, noting the standard transformation from equatorial to Galactic coordinates (Binney & Merrifield, 1998). As in Widrow et al. (2012) and Yanny & Gardner (2013), we study the vertical symmetry of the Galactic disk by comparing the number of stars observed north and south. For continuity and simplicity we follow Widrow et al. (2012) and Yanny & Gardner (2013) in our modelling of the vertical structure of the disk. That is, we modify the solution of Spitzer (1942) for the vertical structure of an axially symmetric thin disk to give a parameterization of the stellar density that has both a thin disk and a thick disk, as inferred from observations (Gilmore & Reid, 1983). Apparently, the two disks have distinct metallicities and scale lengths (Bensby et al., 2014; Bensby, 2014), though we refer the reader to Bovy & Rix (2013) for further discussion; we do not include these refinements, however, in our analysis. We probe the local structure of the stellar disk by determining the parameters of our model across the footprint, or, specifically, as we change the chosen range of longitude for fixed latitude.

In this paper we not only compare the structure, north and south, with longitude, across the footprint, but we also exploit the complete northern coverage of the SDSS to determine (i) the rate at which the mass in the stars changes with distance from the Galactic center, to compare with galactic models, and (ii) how changing the selected latitude interval for a fixed longitude interval impacts the shape of the vertical stellar number count distribution. If there were no change in vertical structure across the Galactic plane and if the vertical metallicity distribution were uniform, then the shape of the vertical distribution in number counts would not change when the latitude window is changed, once corrections are made for the geometric acceptance. We do, however, observe definite changes in shape; consequently, we have a photometric proxy for changes in in-plane structure and/or metallicity. We believe that such studies can be refined and sharpened with further observations from the Gaia mission (Gaia Collaboration, 2016a; Gaia Collaboration et al., 2016b), and not only through more photometric observations, because many more spectra will also become available. Spectral information makes it possible to measure the vertical gradient in metallicity (Hayden et al., 2014), and further studies with better resolution across the Galactic plane should help reveal the specific origin of the shape differences.

We begin by discussing our data selection and consider the various systematic effects that must be understood to determine the vertical distribution of stellar counts. We then turn to a discussion of our fitting procedure before presenting the results of our north-south combined analysis and north-only analysis.

2 Data Selection and Photometric Distance Assessment

This study uses data from the SDSS DR9 (Ahn et al., 2012), the Ninth Data Release of the SDSS, which was also used in Widrow et al. (2012) and Yanny & Gardner (2013). As in our earlier papers, the errors are predominantly systematic and derive from the use of a photometric “parallax” relation to relate the color of a main-sequence star to its intrinsic luminosity and thus to determine its distance. Particular sources of error include those associated with stellar identification, as well as possible inadequacies in the application of corrections for reddening and absorption due to dust. The first source of error includes not only the possible confusion of dwarfs with giants, or the appearance of unresolved binaries, but also the possibility of admixtures of stellar populations of varying age and metallicity. We refer to Yanny & Gardner (2013) for an extended discussion and for tests of our corrections for dust effects. We follow the procedures described in Yanny & Gardner (2013) and now summarize them briefly.

To minimize dust effects, we analyze data only above $30^{\circ}$ in absolute latitude, and we correct the data for the reddening and absorption due to dust by using the maps of Schlegel et al. (1998). In selecting the stars for our analysis, we require $15<r_{0}<21.5$, where $r_{0}$ is the apparent brightness of the star, and $1.8<(g-i)_{0}<2.4$. Note that $g-i$ is the color defined as the ratio of the intensities of the $g$-band and the $i$-band, which are parts of a five-band survey where each band is associated with a range of frequencies, whereas $(g-i)_{0}$ is the color once the effect of reddening from dust has been removed. For our selection of $|b|>30^{\circ}$, the color excess $E(B-V)$ is typically less than 0.03 mag. This translates to a color correction in $g-r$ of about 0.03 mag, and an error on that correction of much less than that amount. Recall that a color error of 0.01 mag corresponds to a distance error of $<5\%$ and that the typical SDSS photometric color errors are about 0.02 mag for fainter stars.

The majority of stars in our data selection are K- and M-type main-sequence stars, which are the reddest and coldest. They make up the majority of the main sequence, with M-type stars being the most prevalent. By using the reddest stars, we reduce the possibility of pollution from non-main-sequence stars. However, it is still possible for giants to infiltrate the sample. Figure 1 shows the results of a photometric test for the presence of giants. This plot analyzes stars with $14.9<r_{0}<15.4$ in a region of the sky where giants are expected to occur; we thus expect our test sample to have a much greater proportion of giants than our analysis sample. We have found the number of giants in our photometric test sample to be a very small fraction of the total number of stars and thus conclude that the giant admixture in our data set is trivially small.

In order to determine the distance to each of our selected stars, we use the photometric parallax relations devised by Jurić et al. (2008) for $(r-i)_{0}$ color and by Ivezić et al. (2008) for $(g-i)_{0}$ color, along with the refinements of Yanny & Gardner (2013) from fits to globular clusters for the red stars we consider here — we refer the reader to Yanny & Gardner (2013) for a detailed discussion and comparison with Widrow et al. (2012). The photometric parallax method determines the distance to a star from its apparent brightness and its intrinsic brightness — the latter is inferred from its color and metallicity. The relation we have used in $(r-i)_{0}$ color for the absolute magnitude $M_{r}$ is that of Jurić et al. (2008) but also includes a metallicity correction $\Delta M_{r}$:

[TABLE]

with $\Delta M_{r}([\hbox{Fe/H}])=-1.11[\hbox{Fe/H}]-0.18[\hbox{Fe/H}]^{2}$ (Ivezić et al., 2008). We have taken $[\hbox{Fe/H}]$ to be $-$0.3 universally because we are unable to use the photometric metallicity assessment of Ivezić et al. (2008) due to the absence of sufficiently precise $u$-band information for the faint, red stars we consider in this paper. Moreover, for calibration purposes, only the stellar clusters M67 and NGC 2420 are sufficiently red and out of the Galactic plane (Yanny & Gardner, 2013), and their metallicities are much larger than expected for thick-disk stars. (In Section 3.3 we report evidence for decreasing metallicity in the K/M dwarfs as we sample them well above the Galactic plane.) A similar color-magnitude equation, adapted from Ivezić et al. (2008), has been used in $(g-i)_{0}$ color:

[TABLE]

with the same metallicity correction. Using the intrinsic brightness $M_{r}$ and the apparent brightness $r_{0}$ in magnitudes, we calculate the distance $d[r_{0},M_{r}]$ in kiloparsecs using the relation $r_{0}-M_{r}=-5+5\log_{10}d[r_{0},M_{r}]$. We refine these distances using the analysis of Yanny & Gardner (2013), so that, for $(r-i)_{0}$ color, the distance $d$ is finally determined to be

[TABLE]

whereas in $(g-i)_{0}$ color, it is

[TABLE]

For the faint, red stars we analyze in this paper the color-color diagram, as shown in Figure 14 of Yanny & Gardner (2013), in $(r-i)_{0}$ versus $(g-i)_{0}$ shows a tight correlation. We thus expect the two photometric parallax relations to work comparably well, and this is borne out by both our findings and those of previous studies (Widrow et al., 2012; Yanny & Gardner, 2013).

The color precision of the SDSS photometry ($\pm 0.02$ mag) leads to typically $\pm 0.2$ mag precision in $M_{r}$ and thus $\pm 10\%$ distance errors (Jurić et al., 2008). This assumes that the assigned metallicity is correct. If an incorrect metallicity were used, then the distance errors would be larger. Supposing, as usual, that the stars close to the plane (thin disk) have an [Fe/H] metallicity of $-0.3$, but the stars far from the plane (thick disk) have an [Fe/H] metallicity of $-0.8$, then one has an additional error of up to $+0.4$ mag (-20% in distance) for a full mismatch of $\Delta{\rm[Fe/H]}=-0.5$ (Ivezić et al., 2008). We thus expect the reported distances to become gradually too large as the population shifts from mostly thin-disk stars at $|z|<0.5$ kpc to that of mostly thick-disk stars at $|z|>1.5$ kpc because the distances have been overestimated by 10%-20% as a consequence of our metallicity assumption.

In this paper we consider a sample of stars several times larger than those analyzed previously, to the end of discerning and interpreting variations in the stellar number distributions across the footprint. As a result, in the current analysis, systematic errors associated with the possible inadequacy of the dust corrections for reddening and extinction across the sky, as well as with any nonuniformity in the photometric calibration itself, become potentially pertinent and need to be considered. We now address each of these effects in turn.

We have used the dust maps of Schlegel et al. (1998) to assess the effects of dust; they are constructed from observations of dust emission in the far infrared after correcting for its temperature. Photometric (Schlafly et al., 2010) and spectroscopic (Schlafly & Finkbeiner, 2011) tests have shown these maps to be accurate for the higher-latitude ($|b|>30^{\circ}$) sightlines we employ in our study, though this work has also revealed, working in $(g-r)_{0}$ color and averaging over large angular scales, that the stars in the south are redder than those in the north. Quantitatively the two assessments are not the same (note Table 5 in Schlafly & Finkbeiner (2011)), with the photometric assessment using blue-tip stars being larger, yielding a result of 21.8 mmag redder in the south rather than $8.8\pm 1.5$ mmag. We believe that the common reddening difference is better attributed to an inadequacy in the calibration of the $g_{0}$-band magnitude in the south, rather than to one in the dust reddening correction per se (Betoule et al., 2013; Yanny & Gardner, 2013), though metallicity variations in the photometric sample can appear as well. Indeed Schlafly & Finkbeiner (2011) suggests that the blue-tip population could be fundamentally different in the north, due to age and metallicity differences. Subsequent work on different fronts supports this view.

Photometric studies with Pan-STARRS1 of stellar reddening under assumptions of variation in stellar metallicity and population across the Galactic disk have recently been employed to construct an independent dust map, based on dust absorption (Schlafly et al., 2014; Green et al., 2015); these maps agree closely with the older emission-based maps if compared out of the Galactic plane (Schlafly et al., 2014). using Pan-STARRS1 data, and its claimed stability of $(20,10,10,10,20)$ mmag in $ugriz$ bands has been confirmed in the northern hemisphere by Finkbeiner et al. (2016). Additionally, they note small differences in the median magnitudes of Pan-STARRS1 and SDSS data, north and south, of less than $\pm 10\,{\rm mmag}$ in the rms for $|b|>20^{\circ}$ and $griz$ colors. Were these shifts uniform over smaller angular scales, then they would be of little relevance. However, their maps of the differences between Pan-STARRS1 and SDSS data reveal that they can differ up to $\pm 40$ mmag over $15^{\prime}$ pixels. Nevertheless, the shifts in $(r-i)$ color are smaller than those in $(g-i)$ color, even in localized regions (Finkbeiner et al., 2016). Thus we think that if our conclusions are insensitive to the choice of the photometric parallax scheme then they should also be robust with respect to a refinement of the photometric calibration. Nevertheless, in later sections we consider how the errors we discuss compare to the size of the effects we discover.

3 Functional Fits

A geometric selection function is used in order to correct for geometric effects and to calculate a stellar density. The selection function determines the effective volume for each bin of observed stars in height $z$ from the Galactic plane. Each bin is divided by the associated selection function to determine the number density at that $z$. The selection function is

[TABLE]

where $\delta$ is the width of the bin in kpc and $z$ is in the middle of the bin and is given in kpc. Note that $l_{1}<l<l_{2}$, $b_{1}<|b|<b_{2}$, and both quantities are in radians. The coordinates

[TABLE]

are all measured from the Galactic center with the convention that the Sun is located at (-8, 0, 0) kpc in this coordinate system. We define distances from the Sun in the Galactic plane as $\Delta x$ and $\Delta y$. Figure 2 shows a plot of $z$ versus $\Delta x$ or $\Delta y$, demonstrating how different slices in latitude and $z$ yield a different projected area on the $z-x$ or $z-y$ plane and impact the volume ${\cal V}(z)$, Slices of different width in longitude also affect ${\cal V}(z)$, but the specific longitudes themselves do not.

The stellar densities $n(z)\equiv n_{\rm raw}(z)/{\cal V}(z)$ are then plotted and fit to the thin- and thick-disk model of Widrow et al. (2012), a modification of the model of Spitzer (1942) that adds a second disk:

[TABLE]

where $n_{0}$ is a normalization factor, $H_{1}$ is the scale height of the thin disk, $H_{2}$ is the scale height of the thick disk, and $f$ is the fraction of the stellar population that belongs to the thick disk. The parameter $z_{\odot}$ is the height of the Sun above the Galactic plane, though we emphasize that $z_{\odot}$ does not correspond to a point-to-point distance because the Galactic plane is a locally determined quantity: the Galactic disk need not be perfectly flat. Our model assumes a two-disk structure that is north-south-symmetric with respect to the Galactic plane. Due to brightness saturation effects, we analyze only those stars that possess a vertical height $|z|>0.35\,{\rm kpc}$ from the Galactic plane. Thus we believe that we are considering a mix of thin- and thick-disk stars across the SDSS footprint. It is pertinent to compare our choice with other recent models of the structure of the Galactic disk. The recent work of Robin et al. (2014) and Bovy (2017), for example, also supports the use of a ${\rm sech}^{2}$ vertical distribution for thick-disk stars in the Milky Way, though the thin-disk model of Czekaj et al. (2014) employs Einasto ellipsoids (generalized exponentials). In addition, the study of Schwarzkopf & Dettmar (2000) suggests, rather, that the vertical distribution of stars in galaxies with a merger history might be better described by a ${\rm sech}^{1}$ distribution. Since our particular purpose is to study the appearance of north-south (N/S) symmetry breaking across the Galactic plane, we believe that our conclusions do not rely on the particular N/S symmetric parameterization of the vertical distributions that we use. If the stellar disk is flat over small variations in $x$ and $y$, rather than merely flat on average, then one can expect $z_{\odot}$ to be universal for all lines of sight. However, if the stellar disk has local warps or ripples, then the determined $z_{\odot}$ need not be universal, and can serve as a proxy for such local variability. Thus variations in the determination of $z_{\odot}$ in the context of a fit to a distribution that is N/S symmetric probe how well one can meaningfully define the Galactic “plane.”

To fit Eq. (7) to the vertical distribution of star counts, we optimize $\chi^{2}$, namely,

[TABLE]

over $N_{\rm bin}$ bins, where $N_{i}$ stars are in bin $i$ centered on a height $z_{i}$ and $n(z_{i})$ is the theoretical model of Eq. (7). The weights $\sigma_{i}$ are given by $\sqrt{N_{i}}$ of $N_{i}$ stars. We use a combination of C and Python and the fitting routines of ROOT 111https://root.cern.ch, a powerful data analysis framework. At low $z$, the selection function becomes ineffective if the stars cannot be seen because they are too bright. A minimum $z$ cut is employed to remove this effect. This cut was determined for each data selection by fitting the model with trial cuts and minimizing the $\chi^{2}/{\rm dof}$ of the fit; we have determined that a minimum $z$ cut of $\lvert z\rvert>0.35\,{\rm kpc}$ can be used for all our fits.

Thus far we have addressed the vertical structure of the disk exclusively. However, in certain stellar data sets, the variation in in-plane radial distance $R$, noting $R\equiv\sqrt{x^{2}+y^{2}}$, can exceed 1.5 kpc. To investigate the consequences of such variations, we also employ fits in which we bin our data set in $R$ as well and scale out the $R$ variation, noting a similar procedure in Banik et al. (2017), by modifying the selection function of Eq. (5) to include the factor

[TABLE]

where $R_{0}=8\,{\rm kpc}$ and $2.7\,{\rm kpc}$ is the central value of the radial scale length (with an uncertainty of 4%) found by Bovy & Rix (2013). That is, after binning the data in $R$ and $z$, where $R_{i}$ is the coordinate at the center of a bin in $R$, we construct

[TABLE]

after summing over all over bins in $R$ for fixed $z$, and then compute $n(z)=n_{\rm eff}(z)/{\cal V}(z)$. As the result of this procedure we end up with larger relative errors in $n(z)$, though we find no indications that such an analysis is warranted. We find that using $n_{\rm eff}(z)$ in place of $n(z)$ does not affect our extracted parameters in a significant way and thus does not affect our conclusions. Thus in what follows, we do not employ a simultaneous analysis in $R$ and $\phi$ outright, but rather bin our data in the vertical coordinate using the selection function of Eq. (5) to determine $n(z)$.

3.1 North and South Combined Analysis

We begin by studying regions selected such that the north and south have completely matched coverage. The first selection has been designed for uniformity by selecting stars with Galactic coordinates in $70^{\circ}<l<165^{\circ}$ and $51^{\circ}<|b|<65^{\circ}$, which are then analyzed in longitude increments of $5^{\circ}$, $10^{\circ}$, or $15^{\circ}$ so that each selection has roughly the same number of stars (Figure 3(a)). This uniform latitude cut creates a sample that can be directly compared as a function of longitude without the added complication of changing latitude. The second selection has been made to maximize the number of stars included in the analysis, by dividing the sample into regions that cover $5^{\circ}$, $10^{\circ}$, or $15^{\circ}$ of longitude, and choosing the maximum range in latitude possible in each wedge while maintaining complete coverage in the north and south (Figure 3(b)). The $R$ and $\phi$ distributions of the $(l,b)$ wedges in Figure 3(b), for the particular range of colors $1.8<(g-i)_{0}<2.4$ that we analyze, are shown in Figure 4.

Using the previously discussed cuts and selections, we apply the selection function in Eq.(5) and employ the fitting function of Eq. (7) with the photometric parallax method in $(g-i)_{0}$ color as described in Section 2 to determine the distance to each star. We require $0.35\,{\rm kpc}\leq z\leq 2.0\,{\rm kpc}$, where the maximum value of $z$ is determined by color completeness (Yanny & Gardner, 2013). We first analyzed the region shown in Figure 3(a) to determine what variations occur when changing the selected longitude range and then repeated our analysis with the expanded latitude sample of Figure 3(b). There are 1.05 and 3.59 million stars in total in the selections of Figures 3(a) and  3(b), respectively. The changes in Galactic parameters for the two sets are sufficiently small that we have shown the stellar density profiles for the expanded latitude selections in Figure 3(b). Figure 5 shows the results, displaying the stellar number density as a function of vertical distance. We find the quality of fit ($\chi^{2}/{\rm dof}$) to range from 1.4 to 2 for the fits shown. Although the vertical stellar distribution is grossly smooth and symmetric, the plots do show a definite north-south asymmetry of ${\cal O}(10\%)$ in the stellar number counts, with the effect becoming more marked at larger $l$, as noticed in Yanny & Gardner (2013), though the expanded latitude data set is four times larger. The wave-like nature of the departures from the north-south-symmetric fit argues against either metallicity variations or a failure of the north-south photometric calibration in explaining this result; indeed the possible sizes of such effects are too small to have an impact. We discuss the residuals, north and south, in Section 3.2. We note that the quality of fit is poorest in regions where the observed wave-like nature of the N/S asymmetry is most prominent. In the remainder of this section we consider the changes in the Galactic parameters across the regions described in Figures 3(a) and (b), to provide a context for our analysis of the N/S variations.

Figures 6, 7(a), and 7(b) show the fit parameters $z_{\odot}$, thin-disk thickness, and thick-disk thickness, respectively, as functions of longitude. Our analysis shows only modest variation in $z_{\odot}$ across the footprint, both for the selections with a uniform latitude-cut (which we henceforth refer to as “uniform”) and for the expanded latitude selections (which we henceforth refer to as “expanded”). The average $z_{\odot}$ values obtained for the uniform and expanded analyses, respectively, are $z_{\odot}=14.9\pm 0.5$ pc and $z_{\odot}=15.3\pm 0.4$ pc. The average for the uniform analysis is consistent with the value of $z_{\odot}=14.3\pm 0.6$ pc found by Yanny & Gardner (2013). Fitting the expanded results with a straight line reveals a $\chi^{2}/{\rm dof}$ of 3.53, so that there appears to be a genuine variation with longitude. However, there is little difference in the average if the uniform or expanded latitude selection is employed. Our values for $z_{\odot}$ are about 2$\sigma$ lower than the recent literature using similar stars, noting $z_{\odot}\sim 25$ pc (Jurić et al., 2008) and $z_{\odot}=27.5\pm 6\rm pc$ (Chen et al., 1999, 2001). Other determinations of $z_{\odot}$, such as the value of $26\pm 3$ pc found in a study of Cepheid variables (Majaess et al., 2009), or that of Joshi (2007), using younger population tracers, with values ranging from 6 to 20 pc, differ from these results and support the existence of the “environmental” sensitivity we have found (see also Bovy (2017)). Figure 7 shows the analogous results for the disk thicknesses. The differences between the averages for uniform and expanded thickness are small, particularly for the thin disk, which is most pertinent to the north-south asymmetry seen in Figure 5, as will become clear in Section 3.2. Recall, too, that the errors in photometric distances are at least $\pm 10\%$. The evidence for variations with $l$ is also weaker than in the case of $z_{\odot}$, though a change in the vertical scale height farther from the Galactic center — and hence at larger $l$ — is expected (Kent et al., 1991; Narayan & Jog, 2002). Sensitivity to such variations is evidently limited by our statistics, so that we will return to this point in our north-only analysis, for which we have a vastly larger number of stars.

Although other analyses have revealed a range in the value of $z_{\odot}$, determined from stars of differing spectral class in differing regions of the sky, our analysis is the first to reveal variations in its value across the Galactic plane with stars of the same spectral class. Barring the existence of significant in-plane metallicity gradients, which would be at odds with the results of existing observational studies (Hayden et al., 2014), we believe the variations we have found speak to a Galactic disk that possesses ripples and thus is not globally flat. We have confirmed that this feature, and indeed all the features we have found thus far, also appear in our analyses using photometric parallax with $(r-i)_{0}$ color. We note, moreover, that a positive warp in the north, toward $l\sim 90^{\circ}$, has also been observed (Russeil, 2003), though this does not explain the variations in $z_{\odot}$ that we have observed. The ripples could potentially arise from tidal effects on the disk, which appear in numerical simulations of disk-satellite encounters (Gómez et al., 2013; Widrow et al., 2014; Gómez et al., 2016, 2017; Laporte et al., 2016). However, it is also important to compare our results with known in-plane features that vary across the SDSS footprint, and thus we consider the Milky Way’s spiral arm structure.

The dust obscuring the Milky Way’s disk has limited our ability to study its spiral arms. Questions persist as to their nature and origin (Baba et al., 2009; Sellwood, 2011; Sellwood & Carlberg, 2014), as well as to their precise structure (Georgelin & Georgelin, 1976; Levine et al., 2006; Hou et al., 2009; Lépine et al., 2011; Francis & Anderson, 2012; Camargo et al., 2013, 2015c; Bobylev & Bajkova, 2014; Griv et al., 2014; Pettitt et al., 2014; Vallée, 2014). The absence of bright, massive 13CO clouds in the interarm space (Roman-Duval et al., 2009) promotes the association of these apparently short-lived features with the spiral arms themselves. As a result, young, red embedded clusters of stars should act as tracers of the spiral arms (Camargo et al., 2013, 2015a, 2015b, 2015c). Camargo et al. (2015a) have challenged the notion of spiral arms as in-plane structures with north-south-symmetric vertical extent (Cox & Gómez, 2002; Monari et al., 2016), revealing embedded clusters with locations extending both above and below the plane. These authors also note that cluster formation can be associated with the halo as well (Camargo et al., 2015b, 2016). Nevertheless, it is possible that the variations in $z_{\odot}$ are related to an out-of-plane structure of the spiral arms. Figure 4 shows the spiral arm structure detailed in Camargo et al. (2015c) overlaid on the analysis sample of Figure 3(b) plotted in Galactocentric coordinates $R,\phi$, and the analyzed wedges in longitude have been marked with a $+$ in regions where the value of $z_{\odot}$ is greater than the global average and with a $-$ in regions where it is less than the global average. While the work of Camargo et al. (2015a, b, c) shows evidence for some parts of the Milky Way’s stellar spiral arms being located a few hundred parsecs above or below the Galactic plane, we do not see any correlation between these out-of-plane stellar spiral arms and variations in $z_{\odot}$ with $l$ or any correlation with regions of the Galactic plane where significant north-south asymmetries are apparent at $|z|>0.35$ kpc.

3.2 North and South Comparative Analysis

We now turn to a detailed analysis of the north-south differences shown in Figure 5. Working in a coordinate system with the Galactic plane as its origin (by shifting $z$ so that $z\rightarrow z+z_{\odot}$), we show residuals calculated as (data-model)/model as a function of vertical distance in Figure 8. We overlay the residual in the north and that in the south in each region of Figure 5 using the parameters of its fit. This visualizes the north-south asymmetry of Figure 5 and confirms the earlier results of Widrow et al. (2012) and Yanny & Gardner (2013). The residuals have wave-like features that grow in amplitude and vertical extent with increasing longitude, refining and extending what was observed by Yanny & Gardner (2013). Interestingly, at low $z$ the south always has a positive residual, whereas the north always has a negative residual. Note that although the figures go out to 2.5 kpc, the fits themselves only extend to 2 kpc. Repeating this analysis on regions with uniform latitude reveals the same pattern; the changes due to latitude are small compared to the effects seen in the figure. The data set with the expanded latitude selection is used here because it allows for a larger range in longitude and better statistics. We find that similar results emerge when we repeat our analysis using distances computed using $(r-i)_{0}$ color.

Given the visual differences between the north and the south shown in Figure 8, a quantitative comparison has been made by performing fits on the north and south independently. We employ the same thin- and thick-disk model but for one modification. The original model had difficulties fitting stars that are only in the north or only in the south; this is solved by fixing the value of $z_{\odot}$ — the other parameter values are determined by fitting. We use a $z_{\odot}$ of $14.3\,{\rm pc}$, the average value obtained by Yanny & Gardner (2013), though our results are not sensitive to that particular choice. The results of these fits are shown in Figure 9, which shows the thicknesses of the thin and thick disks as functions of longitude. In these fits we found it pertinent to employ the uniform sample of Figure 3(a) in order to focus on the possibility of longitudinal variations in a crisp way.

Separate fits to the northern and southern samples reveal a much greater thickness in the north for both the thin and thick disks, though it is also the case that the fraction $f$ of the thick disk relative to the thin disk is also smaller in the north. We show an explicit illustration of this correlation in Table 1, which lists the fit parameters for the two-disk component fits (note Eq. (7)) to the matching north and south data samples in bins of $5^{\circ}-15^{\circ}$ bins in longitude (see Figure 3(a)), including the normalization $N_{0}$, the thin- and thick-disk scale heights $H_{1}$ and $H_{2}$, respectively, the relative fraction $f$ of the thick- and thin-disk populations, and the $\chi^{2}/{\rm dof}$ for each fit. We show the fit results for the sample in the north fitted by itself, the sample in the south fitted by itself, and for the combined north-south sample. In the combined north-south fits, $f$ is approximately 10%, and the thin- and thick-disk scale heights are 0.24 and 0.7 kpc, respectively, with a relatively good $\chi^{2}/{\rm dof}$. When the fits are done separately for the north and south, however, the scale heights for both the thin and thick disks are always larger in the north, by about 20%-30%, with a corresponding decrease in $f$, by about 10%-50%. This strong correlation between $f$ and the scale heights has been noted several times previously in the literature, and it has been speculated that it is a consequence of a fitting degeneracy, rather than an indication of a variation in the scale heights and $f$ with position or hemisphere. We refer to Figure 9 of Chen et al. (2001), as well as Figure 21 of Jurić et al. (2008) for concrete examples. The combined values of the scale heights and thick-disk fraction are in good agreement with Jurić et al. (2008), which analyzed observations from both the northern and southern hemispheres, finding $H_{1}=0.24$ kpc and $H_{2}=0.8$ kpc with $f$ around 0.1; they, too, noted a strong correlation between $f$ and the disk scale heights. Given the good $\chi^{2}/{\rm dof}$ of the separate north and south fits shown in Table 1 (to our knowledge such separate fits have not been done before), and the poorer combined fits, we suggest, rather, that the difference in scale heights and fractions with hemisphere is a real effect. Its existence supports that of a large north-south asymmetry in the vertical stellar distribution in the vicinity of the Sun. Figure 10 shows this is not due to calibration errors.

The N/S differences shown in Figures 8–10 are depicted in an $(l,b,z)$, Sun-centered coordinate system. One may also present the same star-count data in an $(R,\phi,z)$ coordinate system based on the Galactic center, and this is shown in Figures 11 and 12. Panels (a) and (b) of Figure 11 shows two azimuthal bins of star-count data from the footprint of Figure 3(a) at radial distances up to 0.5 kpc from the sun, and up to 3 kpc from the Galactic plane. The pattern of asymmetries seen in several panels of Figure 8 is clearly present in panel (a), which combines several $(l,b)$ bins into one nearby $(R,\phi)$ bin. Panels (c)-(e) offer $(R,\phi)$ presentations of the asymmetries of the data in the expanded footprint of Figure 3(b). While the selection function of the data points with $|z|$ has not been applied, though it has been in Figure 8, the data samples are in all cases taken from matching north and south regions, plotting $|z_{\rm GP}|=|z+z_{\odot}|$ with $z_{\odot}$ = 14.9 pc. Note that the N/S thickness difference emerges as a gross feature close to the Galactic plane in both the uniform and expanded latitude samples.

Figure 12 offers the broadest look at the data arranged in Galactic $(R,\phi)$ coordinates. Here raw star counts, uncorrected by the selection function, are compared in N/S matched samples. In this figure we relax the $|z|<$ 2 kpc constraint for better illustration. The circular inset at left portrays the plane of the Galaxy marked with the azimuthal coordinate $\phi$ and radial coordinate $R$, where ($R_{\odot},\phi_{\odot}$) = (8 kpc, 180∘). We have broken our coverage in the ($R$, $\phi$) plane into smaller regions indicated by the small solid dots, and we have picked out the dots (regions) with the most star counts, as indicated by the red box, for explicit illustration. The faint large gray arcs indicate curves of constant $R$ in steps of 0.5 kpc. Each row of overlapping inset panels selects a different $\phi$ range, as marked. Each panel shows in black (red) the northern (southern) Galactic hemisphere counts as a function of $|z_{\rm GP}|$ from 0-6 kpc, where $|z_{\rm GP}|$ has been adjusted to account for a universal offset of the Sun above the plane of $z_{\odot}$ = 14.9 pc. The fact that the north and south counts converge at large $|z_{\rm GP}|$ suggests that this offset is appropriate. While the fine asymmetrical structure in $(l,b)$ is not visible, one still may clearly see, especially in the panel for $8<R<8.5$ kpc, $175^{\circ}<\phi<180^{\circ}$, that there are vertical oscillations in the stellar density near the Sun. Several other panels show an excess of counts in the south at vertical distances up to 1 kpc from the plane. We see explicitly that (i) there are more stars in the south over most of the footprint and (ii) the feature of Figure 11(c), which drives the wave-like N/S asymmetry we have found, as well as the N/S variation in Galactic parameters, also breaks axial symmetry. It also has no analogue at slightly larger $R$. We encourage further studies of the combined stellar and dark matter mass distributions in the solar neighborhood to delve into its origin.

It should be noted that a partial explanation of these differences could come from a global calibration difference. The photometric calibration of the SDSS survey is done by studying adjacent regions across the footprint and calibrating overlapping regions together (Padmanabhan et al., 2008). Since the entire footprint in the north is connected, it can be calibrated in a uniform way. However, the northern and southern regions of the footprint do not overlap, so that they cannot be calibrated in this manner. Studies of the blue-tip stars suggest that the stars are systematically redder in the south (Schlafly et al., 2010; Yanny & Gardner, 2013). To study the implications of this, a downward shift of -16 mmag for $(g-i)_{0}$ color was implemented in the south, which would redress the typical differences in reddening found by Yanny & Gardner (2013). By repeating the fits with $(g-i)_{0}$ color as well as repeating the fits using distances determined with $(r-i)_{0}$ color, Figure 10 is obtained. The shift in color calibration has a negligible effect on the thicknesses, as we would expect from our discussion in Section 2. Although using the $(r-i)_{0}$ color relation reveals slightly different thicknesses, the north-south offsets still exist, encouraging the stance that there is a real difference in the stellar densities between north and south.

3.3 North-only Analysis: Scale Height Changes Across the Footprint

Since more stars have been observed in the north we use these observations to study variations in Galactic parameters across the footprint. Once again, we have made two different types of selections. Both divide the footprint into slices of $5^{\circ}$, $10^{\circ}$, and $15^{\circ}$ slices in longitude, but each uses different latitudes. The first analysis uses a fixed latitude range of $50^{\circ}<b<90^{\circ}$ and longitude range $35^{\circ}<l<210^{\circ}$ (Figure 13(a)), in which we have 1.67 million stars in total. The second analysis uses a fixed latitude range of $30^{\circ}<b<90^{\circ}$ but only looks at the longitudes $35^{\circ}<l<60^{\circ}$ and $165^{\circ}<l<210^{\circ}$ (Figure 13(b)). These regions are selected in order to have the largest latitude coverage possible, and we analyze 2.22 million stars in total. Since we use $b>30^{\circ}$ in order to avoid dust effects, this is the maximum range of latitude allowed by our analysis.

Figures 14 and 15 show the thicknesses of the thin and thick disks as functions of longitude and the in-plane radial distance $R$ from the Galactic center for both of the regions depicted in Figure 13. The difference in thickness for regions with the same longitudes but different latitudes reveals that the scale heights are sensitive to the value of the lower latitude cut. We have checked that this result persists once we repeat our fits to a data sample in which the $R$ dependence has been scaled out, as per the modified selection function of Eq. (9). It can also be observed, particularly for the analysis of Figure 13(b) shown in Figures 14 and 15, that there is an increase in thickness, for both disks, with an increase in $R$. The value of $R$ used is the average value of $\sqrt{x^{2}+y^{2}}$ for the stars in that wedge that are included in the fit. The red points, corresponding to the inclusion of stars with latitudes down to $30^{\circ}$, sample distances further from the Galactic center at large $l$ and those closer to the Galactic center at small $l$. In a single-disk model, in particular, if $\rho(z)=\rho_{0}{\rm sech}^{2}(z/(2z_{0}))$, then the surface mass density $\Sigma$ has the form $\Sigma=4\rho_{0}z_{0}$ (Spitzer, 1942; Binney & Tremaine, 2008). Figure 14 shows an increase in scale height of about 10% over a change in $R$ from 7 to 9 kpc, but from our fits we determine the parameter $z_{0}$ to be more nearly constant and, moreover, we infer that $\Sigma$ is decreasing. A change in the vertical scale height with $R$ has been noticed previously by Kent et al. (1991) (and in galaxies other than the Milky Way by de Grijs & Peletier (1997)), with an estimated scale height of $\sim 247$ pc at $R=8$ kpc, in good agreement with our result for the thin-disk scale height. Its origin has been discussed by Narayan & Jog (2002). It is important to note that this analysis covers the entire range of latitude from $30^{\circ}<b<90^{\circ}$; this provides us with more statistics than our analyses limited to higher latitudes.

In all of our results, the thin- and thick-disk thicknesses appear to depend on the lower limit of the latitudes included in the analysis. To show this effect explicitly, Figure 16 shows how the thickness depends on te minimum latitude chosen, for a region with $55^{\circ}<l<60^{\circ}$ and a maximum latitude of $90^{\circ}$. We thus see that a higher minimum latitude cut leads to a larger thickness, for both the thin and thick disks; the change with minimum latitude can be as large as 20%. Since we employ a selection function to eliminate geometric effects and the sizes of the systematic errors that could account for such a change are ruled out (see Section 2), we expect that this apparent change arises from one or more physical effects.

In order to understand this apparent change in thickness with latitude cut, regions of a low minimum latitude cut ($30^{\circ}<b<50^{\circ}$) and a high minimum latitude cut ($50^{\circ}<b<90^{\circ}$) are compared in Figure 17. This has been repeated for all the regions shown in Figure 13(b). Although the shapes do not change grossly, it can nevertheless be seen at high $z$ that the higher-latitude data points have a greater stellar density than the lower-latitude data points. This illustrates how the measured thickness could increase with increasing minimum latitude cut. The apparent increase in thickness could be the result of vertically changing stellar populations or of a change in in-plane structure. The manner in which a change in the selected latitude window changes the sampling of the Galactic plane is illustrated in Figure 2. The different latitude selections also sample the stellar populations at different average heights above the Galactic plane as shown in Table 2. Note that $\langle z_{\rm raw}\rangle$ determines the average in $z$ of the raw stellar number counts — no selection function has been applied. We observe that the approximate difference in $\langle z_{\rm raw}\rangle$ in the two latitude samples is about $0.5$ kpc, with the difference increasing slightly, to about $0.6$ kpc, at large $l$. Since we expect the metallicity of a selected star to be smaller as its value of $z$ grows larger, noting that this has been established in a recent spectroscopic study of red giants (Hayden et al., 2014), we expect that the vertical changes in the metallicity of the stellar populations should play a role to some extent.

We now quantify the latitude-dependent difference in shape in each panel of Figure 17. To do this, we have determined the metallicity shift required to minimize the shape differences between the high- and low-latitude samples. Applying the selection function to each latitude sample, we multiplied the number of counts in each bin in the high-latitude sample by an overall scaling factor such that the total counts for each latitude sample are the same. This has been done to permit a comparison of the vertical distribution of the stars regardless of their net number. We then examined the three different regions in $z$, in which we have seen shape differences, separately. In each one we executed a $\chi^{2}$ test to determine whether the shapes of the two latitude samples are the same, computing, namely,

[TABLE]

where $n_{i}^{>}$ and $n_{i}^{<}$ are the numbers of stars in the (rescaled) high-latitude and low-latitude samples in bin $i$ of $N_{\rm bin}$ equal-width bins. We then modified the metallicity of the higher-latitude sample by steps of -0.01 until $\chi_{b}^{2}$ is minimized. We repeated this for each region in $z$ and then repeated the entire process for each slice of longitude. The results of this analysis are shown in Table 2 as a metallicity shift $\Delta$[Fe/H] in dex and as a shift in the intrinsic magnitude $\Delta M_{r}$ in mmag as per Eq. (2). We can see that the shifts in $\Delta M_{r}$ are typically much larger, particularly at large $z$, than the photometric calibration shifts of less than 10 mmag discussed by Finkbeiner et al. (2016) and thus speak to a physical effect. (We have also checked that similar features emerge when we repeat our analysis in $(r-i)_{0}$ color.) These shifts are largely, but not completely, explained by the common paradigm of two Galactic disks — a thin disk with a [Fe/H] of -0.3 and a thick disk with a [Fe/H] of -0.8 with different scale heights. Different $(l,b)$ bins sample different ratios of thin- and thick-disk stars, and thus our approximation that all stars have a [Fe/H] of -0.3 is more incorrect as $|z|$ moves further above 0.5 kpc. Interestingly, too, the manner in which the metallicity changes, that it decreases substantially at large $z$, is also consistent with the trend found in spectroscopic studies (Hayden et al., 2014), though the change is also larger in the large $l$ samples at large $z$ than the typical shifts observed in those studies (Hayden et al., 2014). It is intriguing that a simple thick- and thin-disk paradigm is not enough to explain the very large shifts seen for $1.7<|z|<2.0$ kpc, toward the anticenter ($165^{\circ}<l<210^{\circ}$), and here a significantly different stellar population, with lower metallicity than the thick disk may be present in significant amounts. This is worth exploring further. It is also possible that the vertical structure changes as one moves across the Galactic plane. We anticipate that the rich spectroscopic data sets to emerge from the Gaia mission (Gaia Collaboration, 2016a; Gaia Collaboration et al., 2016b) will help to resolve the origin of the effect we have found regarding the change in the vertical distributions of stars with latitude. We nevertheless view photometric studies of the sort we have pioneered here as being of continuing utility, because they serve as an efficient way of discerning what the most interesting regions of the sky might be for detailed spectroscopic studies.

4 Summary and Future Prospects

We have used a photometric sample of up to 3.6 million K and M dwarf stars from the SDSS to study the structure of the Galactic disk in the vicinity of the Sun. By selecting different regions of Galactic latitude and longitude, we have been able to study changes in structure as a function of the stars’ location within and above the Galactic plane. As shown in Figure 5, the stellar number count distributions are remarkably smooth across the footprint, though we have discovered significant changes in their distribution nonetheless. Our sample of red, main-sequence stars are the longest-lived stars, so that we would expect its one-body distribution function to be the solution of a collisionless Boltzmann equation, which describes the interaction of a star with a mean-field mass distribution. We have presumed this solution to be both axially and mirror (north-south) symmetric in adopting a ${\rm sech}^{2}(z/2H_{1})$ form for its spatial projection in $z$, $n(z)$, as we have used in Eq. (7). With this expectation, we have analyzed the stars and their distribution through matched observations in the north and the south, as well as through observations made in the north only. The departures we have found break these supposed symmetries, though we cannot definitely say whether this is the result of recent dynamical interactions of the disk stars with external agents or because of other internal effects, such as the presence of the spiral arms or the Galactic bar, so that they do not hold precisely.

Our combined north-south analysis reveals a difference in the stellar densities north and south, providing further evidence for a vertical wave in number counts found by Widrow et al. (2012) and Yanny & Gardner (2013), as well as circumstantial evidence for dynamical symmetry breaking of the presumed integrals of motion. We have presented, in Figure 8 and Table 1, quantitative measures of the differences in star counts in matched areas above and below the plane of the Milky Way as a function of distance in $(l,b)$ bins of ${\cal O}(10^{\circ})$. We find significant offsets from north-south symmetry in individual bins of up to 20% in star counts, even after correcting for the location of the Sun at some 14.9 pc above the plane. These differences cannot be explained solely by misestimations of the metallicities of the stars, because the differences are too large and oscillate in $|z|$ with a vertical period of about $400$ pc.

We have also studied the evolution of the north-south differences with longitude and find that they become more pronounced as $l$ approaches $l\approx 180^{\circ}$. Panels (c)-(d) in Figure 11 (and Figure 12) indicate variations in the azimuthal, in addition to the vertical, direction. Although we have observed vertical wave-like features (see Figure 8), such are not apparent in the breaking of axial symmetry we observe. A close look at Figure 12 reveals, rather, an impulse localized within 5 degrees of the Galactic anticenter. This is intriguingly reminiscent of the numerical study of the impact of the Sagittarius dSph galaxy with the Galactic disk by Purcell et al. (2011). A comparison of the disk thicknesses, north and south, also reveals that the thicknesses are larger in the north, though this is also correlated with a smaller thick-disk fraction in the north. These results remain after correcting for possible color calibration offsets and repeating the analysis with $(r-i)_{0}$ color.

In addition, we have used our comparative analysis to determine the inferred location of the Sun above the Galactic plane across the footpoint: we can say that this distance, $z_{\odot}$, changes across the footprint, speaking possibly to ripples in the stellar density across the plane. No significant correlation is seen between sightlines toward spiral arms where the stellar spiral arms of the Galaxy appear out of the plane by up to several hundred parsecs and the wave-like overdensities we have seen at distances of $>500$ pc — but the data are very limited: such a correlation cannot be ruled out.

Studying the observations in the north exclusively admits the possibility of studying the effects of changing latitude, allowing the analysis of observations up to $b=90^{\circ}$ with improved statistics. We have found that the determined thicknesses, for both thin and thick disks, do depend on the value of the minimum latitude analyzed. An increase in thickness was also found with increasing lower latitude cut, which is likely representative of vertically changing stellar populations or of changes in the in-plane structure. Finally we have compared the shape in the stellar distributions with vertical height for different latitude windows, finding definite changes in shape that we think are reflective of either or both stellar population changes and/or changes in the in-plane structure.

A variety of systematic effects have been considered in this study. The effects of dust are minimized by employing dust corrections and restricting our analysis to out-of-plane sightlines with $|b|>30^{\circ}$. We have also performed a photometric test for giants and have found that their infiltration into our analysis sample is completely negligible. The possibility of a color calibration effect that could be different in the north and south has also been studied explicitly and determined to be insignificant.

Our error in converting from stellar colors to distance is a combination of that in the photometric calibration in magnitude and color and of metallicity misestimation, noting that we have used [Fe/H] of the thin disk for a mixed population of thin- and thick-disk stars in the absolute magnitude calculation. The distance errors from calibration and color error have been estimated to be approximately 10% rms, or 0.2 in absolute magnitude (Jurić et al., 2008), and we have refined our distance assessment using red globular clusters. As for the metallicity error, we reiterate that we fix [Fe/H] to be -0.3, i.e., we assume a 100% thin-disk population for our photometric metallicity correction. An additional error would appear if we were, rather, sampling a 100% thick-disk population with [Fe/H] of -0.8. The empirical effects of sampling different ratios of thin- and thick-disk populations stars in different $(l,b)$ bins are tabulated in Table 2. We leave the exercise of “inverting” the results of this table to deduce the thin-disk and thick-disk fractions in each bin to a future work. We note that the results are broadly consistent with a switchover from a thin-disk to a thick-disk population as one moves from heights of $<0.5$ kpc to over 1.5 kpc above or below the plane, with a few places where halo or other low-metallicity stars, with [Fe/H] of -1.6, may represent a significant fraction of the mix toward the Galactic anticenter.

We ask the question of whether the larger differences of up to 20% in stellar density between the north and south in various longitude bins seen in Figure 8, could be due to a significant difference in the metallicity of the population of stars between the north and the south, rather than a difference in overall star density at a given distance. For instance, following the lower right panel in Figure 8, which has $165^{\circ}<l<180^{\circ}$, $46.6^{\circ}<|b|<64.3^{\circ}$, there is a 20% excess in density of star counts at $|z|=0.7$ kpc in the north versus the south. To explain this difference as a difference in population metallicity and not as an overall stellar density wave, one would need, for instance, a very large difference in the metallicity as a function of $|z|$ in the ratio of thin- to thick-disk stars between north and south, and that large difference in metals is not compatible with the relatively smooth change in metals content of stars seen around the Galaxy (Ivezić et al., 2008; Hayden et al., 2014).

Plotting the matched N/S stellar count data in Galactic $(R,\phi,z)$ coordinates shows a 20% excess of counts in the south just beyond (0-0.5 kpc) the solar radius toward the anticenter at a height of 0.4 kpc below the Galactic plane. This excess becomes a 10% deficit (or excess in the north) once one reaches a height of 0.8 kpc above the plane (see Figure 12). It is interesting to note that at lower vertical heights, the Orion spur is in this same part of the sky. In other directions, as in we see in Figure 8, at similar (sub-kpc) distances from the Sun, the asymmetry mostly manifests itself as an excess of counts in the south at about 0.4 kpc below the plane. In all subsamples of the matched north-south star-count data, we report here that the scale height of both the thin- and thick-disk stellar populations appears to be systematically larger in the north than the south, when the scale heights are fit separately. This may indicate a slight displacement of the center plane of the thick disk above the thin disk near the Sun or other configurations suggesting that the disks are not fully in stationary equilibrium.

These quantitative locations and amplitudes of over- and underdensities should be useful in serving as a constraint on dynamical models of the Galaxy, which would describe its dark matter distribution, its satellites, and its stellar disk as they interact and evolve over time. Insights into the past history and distribution of matter can be inferred by a model that reproduces these disk asymmetries. We note, too, that comparing the vertical distribution of number counts as a function of selected latitude may prove an efficient way of locating possible stellar populations, or streams, of ultralow metallicity for subsequent spectroscopic study. The $\Lambda$CDM model speaks to dark matter with small-scale phase-space structure and motivates the search for stellar streams. Follow-up studies at yet higher resolution with our methods could look for variations and asymmetries with greater sensitivity, though they would likely require improved photometric calibrations and reddening corrections. These improvements already largely exist (Green et al., 2015; Finkbeiner et al., 2016), so that we look forward to data from the Gaia era (Gaia Collaboration, 2016a; Gaia Collaboration et al., 2016b) with confidence in the ability to further and refine the studies pioneered here.

This manuscript evolved from D.F.’s 2016 May senior honors thesis at the University of Kentucky, and she acknowledges the support of The University of Kentucky Singletary Scholarship as well as an Undergraduate Summer Research Fellowship from the University of Kentucky during its completion. D.F. and S.G. acknowledge partial support from the U.S. Department of Energy under contract DE-FG02-96ER40989. S.G. thanks Jonathan Feng and his colleagues at the University of California, Irvine for gracious hospitality during the completion of this work. We thank the referee for several important suggestions that improved the presentation of this paper, including presentation of the data in Galactic coordinates. Funding for the Sloan Digital Sky Survey has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is http://www.sdss.org$\,.$ SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatório Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

∗Present Address: School of Physics, Georgia Institute of Technology, Atlanta, GA 30332, USA