UPDATED KINEMATIC CONSTRAINTS ON A DARK DISK

Since the original study by Oort (1932, 1960), the question of disk dark matter has been a subject of controversy. Over the years, several authors have suggested the idea of a dark disk to explain various phenomena. Kalberla et al. (2007) proposed a thick dark disk as a way to explain the flaring of the interstellar gas layer. Read et al. (2008) showed, using cosmological simulations, that a thick dark disk is formed naturally in a ΛCDM cosmology as a consequence of satellite mergers. A phenomenologically very different idea is that of a thin dark disk. Fan et al. (2013) put forward a model for dark matter, coined Double Disk Dark Matter (DDDM), where a small fraction of the total dark matter is self-interacting and dissipative, forming a thin disk. Following the 2013 DDDM paper, Randall & Reece (2014) showed that a dark matter disk of surface density ∼10 $\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$ and scale height ∼10 pc could potentially explain periodicity of comet impacts on Earth, giving a target surface density and scale height. If such a dark disk were found to exist, we would want to know what are its surface mass density, Σ_D, and its scale height, h_D. Its density in the plane would then be approximately ${\rho }_{D}(0)\simeq {{\rm{\Sigma }}}_{D}/4{h}_{D}$.¹ A dark disk model allows us to test the viability of these benchmark values. We do not assume this model is favored, but since the literature clearly supports no dark disk as a possibility, we ask whether it also allows for a disk with this or greater density.

A possible concern in a thin dark-disk model is if the system can emerge dynamically when instabilities and fragmentation are accounted for. These concerns may be resolved with simulations and more careful theoretical considerations, or possibly additional model building. For the purposes of this paper, however, we consider the dark disk from a purely phenomenological perspective: we ask only what can the data tell us about the presence of a thin dark disk, and how to better determine this in the future.

The existing literature suggests that such a model is highly constrained. In principle, there are many ways in which one might aim to constrain a dark disk kinematically. All of these rely on the Poisson–Jeans theory (cf. Section 2). Based on our survey of the literature, we categorize the attempted constraints into three main categories:

• 1.  
  Using stellar kinematics to fit a model for known matter containing interstellar gas, Galactic disk, and dark halo. The parameter being fit here is the surface density of the Galactic disk. Once this is found, one has a measure of the total surface density of the galactic disk. The authors then compare this result to the inventories of known matter to try to constrain the dark disk surface density.
• 2.  
  Using stellar kinematics to directly measure the vertical gravitational acceleration K_z as a function of height z above the Galactic plane. By the Poisson equation, this is proportional to the total surface density Σ_z integrated to height z above the plane. Again an attempt is made to compare this result to the inventories of known matter.
• 3.  
  Using an isothermal component model (Bahcall 1984a, 1984b) for known matter, and fitting the remaining dark matter mass self-consistently using stellar kinematics.

Table 1 contains a list of such constraints and which category they fall into. In the face of all these constraints, one may wonder what hope is left for a dark disk. We claim, however, that, at least in their present state, the majority, if not all, of the constraints in Table 1 do not apply.

Table 1.  Bounds on the Surface Density and Local Density of Total Matter and Visible Matter in the Galaxy

Authors	Year	Bound (M⊙ pc−n)	Category
Oort	1932	Σtot(100 pc) = 31	2
Oort	1960	Σtot(100 pc) = 29 ± 10%	2
Bahcall	1984a	${{\rm{\Sigma }}}_{D,{\rm{thin}}}\leqslant 17\qquad {\rho }_{{\rm{tot}}}(0)\leqslant 0.24$	3
Bahcall	1984b	${{\rm{\Sigma }}}_{{\rm{tot}}}=55\mbox{--}83\qquad {\rho }_{{\rm{tot}}}(0)\quad =$ 0.17–0.25	3
Bienayme et al.	1987	ρDM(0) ≤ 0.03 for thick dark disk	3
Kuijken & Gilmore	1991	Σtot(1.1 kpc) = 71 ± 6	1
Bahcall et al.	1992	${{\rm{\Sigma }}}_{{\rm{tot}}}={70}_{-16}^{+24}$	3
Flynn & Fuchs	1994	Σtot = 52 ± 13	3a
Pham	1997	ρtot(0) = 0.11 ± 0.01	NA
Creze et al.	1998	ρtot(0) = 0.076 ± 0.015 (assumed constant density)	1a
Holmberg & Flynn	2000	ρtot(0) = 0.102 ± 0.010 ρvis = 0.095	3a
Korchagin et al.	2003	Σtot(350 pc) = 42 ± 6	2
Siebert et al.	2003	${{\rm{\Sigma }}}_{{\rm{tot}}}(800\;\mathrm{pc})={76}_{-12}^{+25}$	1
Holmberg & Flynn	2004	Σtot(1.1 kpc) = 74 ± 6	3
Bienaymé et al.	2006	Σtot(800 pc) = 57–66	1
Garbari et al.	2011	ρhalo = 0.003–0.033	3
Moni Bidin et al.	2012b	Σtot(1.5 kpc) = 55.6 ± 4.7	2
Bovy & Tremaine	2012	ρhalo = 0.008 ± 0.003	2
Zhang et al.	2013	Σtot(1 kpc) = 67 ± 6
		ρhalo(0) = 0.0065 ± 0.0023	1
Bovy & Rix	2013	Σ1100 = 68 ± 4	1
Bienaymé et al.	2014	Σtot(1.1 kpc) = 68.5 ± 1
		${{\rm{\Sigma }}}_{{\rm{tot}}}(350\;\mathrm{pc})={44.2}_{-2.9}^{+2.3}$	1

Note.

^aDenotes bounds derived using HF technique.

Download table as:  ASCIITypeset image

Method 1 assumes that there is no dark disk and fits a model using known matter to stellar kinematics to show that a dark disk is not necessary. For example, Kuijken & Gilmore (1989a, 1989b, 1991) showed that the density and velocity distributions of K-dwarfs in the Milky Way were consistent with little or no disk dark matter (besides the usual dark halo). Although it is correct that the stellar kinematics do not require a dark disk, the converse is not correct: they do not exclude a dark disk. One might then try to exclude a dark disk by comparing the total Galactic disk surface density found in this way to the inventories of known matter, which at face value seems to be correct. For example, using the distribution function techniques developed by Binney (Binney & Tremaine 2008), Bovy & Rix (2013) dynamically constrained the amount of matter in the Galactic disk, claiming that their results, which are consistent with only standard Galactic components, "leave little room for a dark disk component" in the Milky Way.

However, there is an important detail that has so far been overlooked: the dark disk can actually make room for itself. This is because the distributions of known components in the Galactic disk are extrapolated from observations near the midplane. However, depending on the gravitational pull in the disk, this extrapolation can produce thinner scale heights for these components. For example, suppose we know the density, ρ_i, of certain mass components near the plane. The surface density contributions of these components will be Σ_i ∼ ρ_ih_i, where h_i is the thickness of components i. The dark disk will affect this value by "pinching" the matter distributions and reducing the thickness, h_i, of the components, resulting in a lower surface density determination. The surface density estimation for known matter is therefore lower in the presence of a dark disk than without one. Qualitatively, we can write this "triangle inequality" as

$$\begin{eqnarray}{\rm{\Sigma }}({\rm{visible}}\;{\rm{matter}}\;{\rm{+}}\;{\rm{dark}}\;{\rm{disk}}) & \leqslant & {\rm{\Sigma }}({\rm{visible}}\;{\rm{matter}}\;{\rm{alone}})\\ & & +{\rm{\Sigma }}({\rm{dark}}\;{\rm{disk}}).\end{eqnarray} \tag{ 1 }$$

We therefore need to keep this in mind when comparing kinematic determinations of the total surface density to known mass inventories. Figure 1 compares the various bounds in the literature to models with and without a dark disk, as well as to a dark disk model that does not take into account the pinching effect back-reaction of the dark disk.

Zoom In Zoom Out Reset image size

An additional problem with this method is that the specific models that are fit to the kinematics do not allow for a dark disk. Most of the studies using Method 1 (Kuijken & Gilmore 1989a, 1989b, 1991; Bienaymé et al. 2006, 2014; Zhang et al. 2013) fit models containing only a halo, but no thin component. Other studies (Bovy & Rix 2013; Zhang et al. 2013) considered a thin gas component as well but kept its mass fixed and did not allow for any additional mass in a thin component. The study of Creze et al. (1998), in fact, did not consider a stellar disk component at all and assumed a constant density potential. Their result is therefore very difficult to compare to the literature and in particular to a thin dark disk model.

A similar argument would apply to Method 2. Indeed, the results of Bovy & Tremaine (2012) are consistent with a thin dark disk scenario for precisely the same reason. However, it should be noted that the results of Bovy & Tremaine were based on kinematic data high above the Galactic midplane (1.5–4.5 kpc). Uncertainties in asymmetric drift at these heights could possibly allow total surface densities that are much lower. In fact, analyzing the same data, Moni Bidin et al. (2012a) find that the data do not even allow a dark halo, let alone a dark disk. In Moni Bidin et al. (2015a), however, they point out that this method is effective only in constraining the mass above z = 1.5 kpc, i.e., ${{\rm{\Sigma }}}_{{\rm{tot}}}(z)-{{\rm{\Sigma }}}_{{\rm{tot}}}(1.5\;\mathrm{kpc})$. (This must be the case, as their 2012 value of Σ_tot(1.5 kpc) is significantly lower than all literature measurements of Σ_tot(1.1 kpc).) This argument applies to the results of Binney & Tremaine (2008) as well. These results therefore do not apply in constraining a thin dark disk.

Another study measuring Σ_tot(z) directly is that of Korchagin et al. (2003). Their result ${{\rm{\Sigma }}}_{{\rm{tot}}}(10\;\mathrm{pc})=10\pm 1\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$ is highly dependent on their choice for the form of ρ_tot(z). By assuming that ${\rho }_{{\rm{tot}}}(z)\sim {\rm{sech}}(z/2h)$ or ${\mathrm{sech}}^{2}(z/2h)$, they are actually not allowing a thin massive component from the outset. Their only robust result is therefore ${{\rm{\Sigma }}}_{{\rm{tot}}}(350\;\mathrm{pc})=42\pm 6\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$ since it does not depend on the type of extrapolation to z = 0 used. Even at 350 pc, these results seem to be in disagreement with a dark disk of the size we are considering, as can be seen in Figure 1. However, it was pointed out to us (T. Girard 2013, private communication) that the extinction corrections were applied with the wrong sign. Intuitively, without reddening corrections, we expect preferential reddening near the Galactic plane to increase the apparent number of red giants near the plane, giving a cuspy profile such as we would expect to arise in the presence of a thin dark disk. Although reddening corrections should remove this potentially large bias, applying the reddening corrections with the wrong sign should make it twice as large. It is therefore unclear why Korchagin et al. did not obtain a tighter bound on Σ_tot(50 pc). At any rate, we cannot take their present results at face value.

On the other hand, the older studies of Oort (1932, 1960) seem to particularly favor a dark disk model, as can be seen in Figure 1. (Here, we assume the robust result to be Σ_tot(100 pc) and ignore the endpoint result Σ_tot(50 pc).) However, Read (2014) argues that these results were based on "poorly calibrated photometric distances," stars that were too young to be in equilibrium, and other questionable assumptions.

An important lesson here is that, as with particle physics experimental studies, it helps to have a definite model in order to find a self-consistent bound. Previous analyses may not apply because of changes in measured parameters, but also because the constraint itself depends on the model, which determines the distribution of other components. On top of this, most previous studies did not use the dark disk thickness as an independent parameter. With a model it becomes clear that it makes sense to analyze the data with such a parameter included. Given the large amount of data becoming available, it makes sense to view the data as a way of measuring standard parameters, but also constraining—or hopefully discovering—new ones. It is also true that, were one to analyze a different model aside from the dark disk model we have in mind, one may have to do the analysis differently to account for any different parameters that we (and previous authors) did not include.

Method 3, on the other hand, can include a dark disk in a manner that is gravitationally self-consistent. For example, Bahcall (1984a) used this method to show that a thin dark disk similar in size to the interstellar gas disk had to be lighter than ${{\rm{\Sigma }}}_{D}=17\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$ in the solar region. Bienayme et al. (1987) also used this technique to show that a thick dark disk had to be lighter than ${\rho }_{D}(0)\leqslant 0.03\;{M}_{\odot }\;{{\rm{pc}}}^{-3}$, as did Flynn & Fuchs (1994) to show that a model with no dark halo but a very light thick dark disk, with a total galactic surface density Σ_tot = 52 ± 13, was a good match to the kinematics of K giants in the solar region. Most importantly, this method was used by Holmberg & Flynn (2000) to show that the kinematics of A and F stars in the solar region were consistent with visible matter distributions alone, without the need for a dark disk. By adding dark matter to known components, they were able to compute a midplane density for the total matter of $\rho (0)=0.100\pm 0.006\;{M}_{\odot }\;{{\rm{pc}}}^{-3}$. Since they considered dark matter as thin as the molecular hydrogen mass component, their result effectively sets a bound on a dark disk of scale height of 40 pc.

We will see that, using Method 3, the bounds on a dark disk become even tighter as the scale height is decreased. However, there are at least two questions to be asked on the interpretation of this result: (i) Can the combined errors in visible components and kinematics allow for a dark disk on their own? (ii) Do the assumptions of their kinematic methods (e.g., of thermal equilibrium) hold?

To answer the first question we note that a non-zero dark disk mass could certainly have been hiding in the previously assumed error bars on Σ_vis, including the errors on the interstellar gas, which had until now been around 50%. Bovy & Rix (2013) reported an uncertainty on Σ_tot of $\pm 4\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$, and an uncertainty on ${{\rm{\Sigma }}}_{\mathrm{stars}+\mathrm{remnants}}$ also of $\pm 4\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$. If we combine this with the traditionally assumed $\pm 7\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$ error on the interstellar gas, we can already allow a dark disk with mass $\sqrt{{4}^{2}+{4}^{2}+{7}^{2}}\simeq 9\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$.

The baryonic matter components of the Galaxy have, however, been more carefully measured in recent years. Revised values for the interstellar gas parameters are the subject of Kramer & Randall (2016). We also include new values for the stellar components, revised by McKee et al. (2015). We repeated the study of HF2000, including radial velocities (which were not available at the time), and using the recent three-dimensional reddening map of Schlafly et al. (2014). With these updates, we show in this paper that the traditional HF2000 method would not allow a dark disk heavier than $\sim 4\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$.

Regarding the second question, however, we note that one of the major assumptions in the kinematic method of HF2000 is that the populations under study are in statistical-mechanical equilibrium. We claim (and show in Appendix D) that the definition of equilibrium needed for these analyses is one that precludes the possibility of oscillating solutions, in which the tracer population is in a distribution whose center oscillates above and below the Galactic plane. Such behavior, if present, would smooth out the effects of any thin Galactic component, in which case incorrectly assuming a static distribution would result in an upper bound on a dark disk that is too low.

We show that the tracer populations studied by HF2000, indeed, do possess a net vertical velocity, as well as a center that is vertically displaced from the Galactic midplane, which are evidence for such deviations from equilibrium. The HF2000 analysis, which assumes that the tracer populations are in equilibrium, produces an overly strong bound. We show here how a consistent analysis can be performed without any need for the assumption of equilibrium. Analyzed in this way, we show that the kinematics currently allow for a thin dark disk of up to $14\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$, with a weaker bound for thicker disks. (These considerations apply equally to the studies of Flynn & Fuchs 1994, and of Creze et al. 1998. The latter sample, in fact, does show deviations from equilibrium including a negative net vertical displacement from the midplane, as reported by the authors, and a net vertical velocity, as can be inferred from their numbers. This may, along with the reasons described above, account for the unusually low values of ρ_tot(0) found by these authors.)

The Poisson–Jeans theory will be important both for constructing a model of the Galactic potential, and for gaining a qualitative understanding of the effect of a thin dark disk. We will now explain the Poisson–Jeans theory, following which we will use the P–J theory to solve a simple toy model. We will then explain the theory behind the HF2000 study, and derive the Holmberg & Flynn relation.

2.1. Poisson–Jeans Theory

Consider the phase-space distribution f_i(${\boldsymbol{x}},{\boldsymbol{v}}$) for a stellar population i, satisfying

$$\begin{eqnarray}&&\int {d}^{3}v\;{f}_{i}({\boldsymbol{x}},{\boldsymbol{v}})={\rho }_{i}({\boldsymbol{x}})\end{eqnarray} \tag{ 2 }$$

and, for a more general function, g,

$$\begin{eqnarray}&&\int {d}^{3}v\;{f}_{i}({\boldsymbol{x}},{\boldsymbol{v}})g({\boldsymbol{v}})={\rho }_{i}({\boldsymbol{x}})\langle g\rangle \end{eqnarray} \tag{ 3 }$$

where ρ_i(${\boldsymbol{x}}$) is the population density. The f_i must also satisfy the collisionless Boltzmann equation (Liouville's theorem)

$$\begin{eqnarray}&&\displaystyle \frac{{{Df}}_{i}}{{Dt}}\equiv \displaystyle \frac{\partial {f}_{i}}{\partial t}+\dot{{\boldsymbol{x}}}\cdot \displaystyle \frac{\partial {f}_{i}}{\partial {\boldsymbol{x}}}+\dot{{\boldsymbol{v}}}\cdot \displaystyle \frac{\partial {f}_{i}}{\partial {\boldsymbol{v}}}=0.\end{eqnarray} \tag{ 4 }$$

If we assume our stellar population is in equilibrium then the first term on the right vanishes. Also, we can replace $\dot{{\boldsymbol{x}}}$ and $\dot{{\boldsymbol{v}}}$ by ${\boldsymbol{v}}$ and $-\partial {\rm{\Phi }}/\partial {\boldsymbol{x}}$, respectively, where Φ is the gravitational potential. From the Boltzmann equation we can derive the vertical Jeans equation in the standard way, found e.g., in Binney & Tremaine (2008). The vertical Jeans equation then reads:

$$\begin{eqnarray}&&\displaystyle \frac{1}{R}\displaystyle \frac{\partial }{\partial R}(R{\rho }_{i}{\sigma }_{i,{Rz}})+\displaystyle \frac{\partial }{\partial z}({\rho }_{i}{\sigma }_{i,z}^{2})+{\rho }_{i}\displaystyle \frac{\partial {\rm{\Phi }}}{\partial z}=0\end{eqnarray} \tag{ 5 }$$

The first term in Equation (5) describes the rate of change of the correlation between v_z and v_R, so it is commonly referred to as the "tilt" term. The tilt term is expected to be very small near the plane, since, being asymmetric in z, it must vanish at z = 0 (Moni Bidin et al. 2012b). If we restrict our analysis to heights that are small compared to the tracer population's scale height, we can safely neglect this term. We thus have (dropping the subscript on ${\sigma }_{i,z}$)

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial z}({\rho }_{i}{\sigma }_{i}^{2})+{\rho }_{i}\displaystyle \frac{\partial {\rm{\Phi }}}{\partial z}=0\end{eqnarray} \tag{ 6 }$$

which we can rewrite as

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\rho }_{i}{\sigma }_{i}^{2}}\displaystyle \frac{\partial }{\partial z}({\rho }_{i}{\sigma }_{i}^{2})+\displaystyle \frac{1}{{\sigma }_{i}^{2}}\displaystyle \frac{\partial {\rm{\Phi }}}{\partial z}=0,\end{eqnarray} \tag{ 7 }$$

from which we can immediately see that the solution is

$$\begin{eqnarray}&&{\rho }_{i}\;{\sigma }_{i}^{2}\propto \mathrm{exp}\left(-\int {dz}\;\displaystyle \frac{1}{{\sigma }_{i}^{2}}\displaystyle \frac{\partial {\rm{\Phi }}}{\partial z}\right).\end{eqnarray} \tag{ 8 }$$

Additionally, if we take our stellar population to be isothermal, that is σ_i(z) = constant, then the solution reduces further to

$$\begin{eqnarray}&&{\rho }_{i}(R,z)={\rho }_{i}(R,0)\;{e}^{-({\rm{\Phi }}(R,z)-{\rm{\Phi }}(R,0))/{\sigma }_{i}^{2}}.\end{eqnarray} \tag{ 9 }$$

We will now consider only R = R_⊙, where we will fix Φ(R_⊙, 0) = 0. This gives

$$\begin{eqnarray}&&{\rho }_{i}(z)={\rho }_{i}(0)\;{e}^{-{\rm{\Phi }}(z)/{\sigma }_{i}^{2}}.\end{eqnarray} \tag{ 10 }$$

Indeed, once we have assumed equilibrium, vanishing tilt, and isothermality, what we are left with is simply a gas at temperature per unit mass ${{kT}}_{i}/2M={\sigma }_{i}^{2}/2$, and Equation (10) is merely a Boltzmann factor. A possible objection to the method is that it was found by Garbari et al. (2011) that neglecting the tilt term in Equation (5) in general leads to a biased determination of ρ_dm. However, this was only found to be a problem at heights of z ≳ 0.5 kpc or greater. Garbari et al. describe the HF2000 sample, which is relatively close to the plane, as "unlikely to be biased." With this solution to the Jeans equation we can now relate the potential Φ to the mass distribution consisting of the different components ρ_i. By the Poisson equation,

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}{\rm{\Phi }}=4\pi G\displaystyle \sum _{i}{\rho }_{i}.\end{eqnarray} \tag{ 11 }$$

We can split up the Laplacian in the standard way, following Kuijken & Gilmore (1989b):

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}{\rm{\Phi }}=\displaystyle \frac{{\partial }^{2}{\rm{\Phi }}}{\partial {z}^{2}}+\displaystyle \frac{1}{R}\displaystyle \frac{\partial }{\partial R}\left(R\displaystyle \frac{\partial {\rm{\Phi }}}{\partial R}\right)\end{eqnarray} \tag{ 12a }$$

$$\begin{eqnarray}&&=\;\displaystyle \frac{{\partial }^{2}{\rm{\Phi }}}{\partial {z}^{2}}+\displaystyle \frac{1}{R}\displaystyle \frac{\partial }{\partial R}{V}_{c}^{2}\end{eqnarray} \tag{ 12b }$$

$$\begin{eqnarray}&&=\;\displaystyle \frac{{\partial }^{2}{\rm{\Phi }}}{\partial {z}^{2}}+2\displaystyle \frac{{V}_{c}}{R}\displaystyle \frac{\partial {V}_{c}}{\partial R}\end{eqnarray} \tag{ 12c }$$

$$\begin{eqnarray}&&=\;\displaystyle \frac{{\partial }^{2}{\rm{\Phi }}}{\partial {z}^{2}}+2({B}^{2}-{A}^{2})\end{eqnarray} \tag{ 12d }$$

where V_c is the local circular velocity, and A and B are the "Oort constants"

$$\begin{eqnarray}&&A\equiv \displaystyle \frac{1}{2}\;\left(-\displaystyle \frac{\partial {V}_{c}}{\partial R}+\displaystyle \frac{{V}_{c}}{R}\right)\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&B\equiv -\displaystyle \frac{1}{2}\;\left(\displaystyle \frac{\partial {V}_{c}}{\partial R}+\displaystyle \frac{{V}_{c}}{R}\right).\end{eqnarray} \tag{ 14 }$$

We have assumed azimuthal symmetry in Equation (12a). Although in reality the azimuthal symmetry of the Galaxy is spoiled by spiral arms and other structures, we will assume that these are not relevant over the timescales needed for our tracer populations to relax to their current distributions. At the Sun's position the second term in Equation (12a) is very small, so we shall rename it 4πGδρ. From Catena & Ullio (2010), we have A − B = 29.45 ± 0.15 $\mathrm{km}\;{{\rm{s}}}^{-1}\;{\mathrm{kpc}}^{-1}$ and A + B = 0.18 ± 0.47 $\mathrm{km}\;{{\rm{s}}}^{-1}\;{\mathrm{kpc}}^{-1}$, giving

$$\begin{eqnarray}&&\delta \rho =(-2\pm 5)\times {10}^{-4}\;{M}_{\odot }\;{\mathrm{pc}}^{-3}\end{eqnarray} \tag{ 15 }$$

near the z = 0 plane. This is comparable in magnitude to the density of red giant stars and to the stellar halo density (Holmberg & Flynn 2000), but is about 2–3 orders of magnitude smaller than the total density and can safely be neglected:

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}{\rm{\Phi }}}{\partial {z}^{2}}=4\pi G\displaystyle \sum _{i}\;{\rho }_{i}.\end{eqnarray} \tag{ 16 }$$

Combining Equation (16) with the Jeans Equation (10), we have the Poisson–Jeans equation for the potential Φ:

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}{\rm{\Phi }}}{\partial {z}^{2}}=4\pi G\displaystyle \sum _{i}{\rho }_{i}(0){e}^{-{\rm{\Phi }}/{\sigma }_{i}^{2}}.\end{eqnarray} \tag{ 17 }$$

where we have dropped the R coordinate label. This can also be cast in integral form (assuming z-reflection symmetry)

$$\begin{eqnarray}&&\displaystyle \frac{{\rho }_{i}(z)}{{\rho }_{i}(0)}=\mathrm{exp}\left(-\displaystyle \frac{4\pi G}{{\sigma }_{i}^{2}}\displaystyle \sum _{k}{\int }_{0}^{z}{dz}^{\prime} \;{\int }_{0}^{z^{\prime} }{{dz}}^{\prime\prime }\;{\rho }_{k}({z}^{\prime\prime })\right)\end{eqnarray} \tag{ 18 }$$

which is the form used in our Poisson–Jeans solver.

2.2. A Toy Model

We now study the Poisson–Jeans theory in some simple models following Spitzer (1942). This will be useful in understand the qualitative effects of a thin dark disk. We first study a self-gravitating, single-component galaxy. In this case, Equation (17) reads

$$\begin{eqnarray}&&\displaystyle \frac{1}{4\pi G}\displaystyle \frac{{d}^{2}{\rm{\Phi }}}{{{dz}}^{2}}={\rho }_{0}{e}^{-{\rm{\Phi }}(z)/{\sigma }^{2}}\end{eqnarray} \tag{ 19 }$$

which has solution

$$\begin{eqnarray}&&\displaystyle \frac{d{\rm{\Phi }}(z)}{{dz}}=\displaystyle \frac{{\sigma }^{2}}{h}\mathrm{tanh}(z/2h)\end{eqnarray} \tag{ 20 }$$

or

$$\begin{eqnarray}&&\rho (z)={\rho }_{0}\;{{\rm{sech}}}^{2}(z/2h)\end{eqnarray} \tag{ 21 }$$

where

$$\begin{eqnarray}&&h=\displaystyle \frac{\sigma }{\sqrt{8\pi G{\rho }_{0}}}.\end{eqnarray} \tag{ 22 }$$

Although this was a toy model, it illustrates two very important features of gravitating disks: (1) The scale height, h, of the disk grows with its vertical dispersion, σ, approximately as h ∼ σ, and (2) The scale height decreases with the total midplane density. These features will remain qualitatively true even when more components are included.

We can now ask what happens if we include an additional very thin delta-function component (as a toy approximation to a dark disk). In this case, we have

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}{\rm{\Phi }}=4\pi G{\rho }_{s}+4\pi G{{\rm{\Sigma }}}_{D}\delta (z)\end{eqnarray} \tag{ 23 }$$

and the Poisson–Jeans equation takes the form

$$\begin{eqnarray}&&\displaystyle \frac{1}{4\pi G}\displaystyle \frac{{d}^{2}{\rm{\Phi }}(z)}{{{dz}}^{2}}-{{\rm{\Sigma }}}_{D}\;\delta (z)={\rho }_{s0}\mathrm{exp}\left(-\displaystyle \frac{{\rm{\Phi }}(z)}{{\sigma }^{2}}\right)\end{eqnarray} \tag{ 24 }$$

where ρ_s stands for the density of stars and Σ_D is the dark disk surface density. Defining $Q\equiv {{\rm{\Sigma }}}_{D}/4{\rho }_{s0}{h}_{s}$, (with h_s defined as in Equation (22)), we can write down the exact solution, which is

$$\begin{eqnarray}&&{\rho }_{s}(z)={\rho }_{s0}(1+{Q}^{2})\;{{\rm{sech}}}^{2}\left(\displaystyle \frac{\sqrt{1+{Q}^{2}}}{2{h}_{s}}(| z| +{z}_{0})\right)\end{eqnarray} \tag{ 25 }$$

with

$$\begin{eqnarray}&&{z}_{0}\equiv \displaystyle \frac{2{h}_{s}}{\sqrt{1+{Q}^{2}}}{\rm{arctanh}}\;\left(\displaystyle \frac{Q}{\sqrt{1+{Q}^{2}}}\right).\end{eqnarray} \tag{ 26 }$$

Effectively, the dark disk has the effect of "pinching" the density distribution of the other components, as we can see in Figure 2. In particular, it reduces their scale heights, and for a fixed midplane density, ρ_i(0), it implies that their surface density, Σ_i, is less than what it would be without the disk. Actually, this is what would happen were we to include any other additional mass components.

Zoom In Zoom Out Reset image size

An important consequence of this for the stellar disk, especially in the context of our analysis, is that for a fixed midplane density the effect of the dark disk is to decrease the stellar surface density. Integrating Equation (25), the stellar surface density find is

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{s}(Q)={{\rm{\Sigma }}}_{s}(0)(\sqrt{1+{Q}^{2}}-Q)\end{eqnarray} \tag{ 27a }$$

$$\begin{eqnarray}&&=\;\sqrt{{{\rm{\Sigma }}}_{s}{(0)}^{2}+{{\rm{\Sigma }}}_{D}^{2}\;}-{{\rm{\Sigma }}}_{D}\end{eqnarray} \tag{ 27b }$$

where ${{\rm{\Sigma }}}_{s}(0)\equiv 4{\rho }_{s0}{h}_{s}$ is what the surface density would have been without the dark disk. We can see that the stellar surface density, Σ_s(Q), is monotonically decreasing with Q. We can then write the total surface density as

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{\rm{tot}}}={{\rm{\Sigma }}}_{s}(Q)+{{\rm{\Sigma }}}_{D}\end{eqnarray} \tag{ 28a }$$

$$\begin{eqnarray}&&=\;\sqrt{{{\rm{\Sigma }}}_{s}{(0)}^{2}+{{\rm{\Sigma }}}_{D}^{2}}\end{eqnarray} \tag{ 28b }$$

verifying the triangle inequality (Equation (1)) for this simple case.

2.3. The HF Study

We now explain the analysis of Holmberg & Flynn (2000), and derive the HF relation. This analysis, developed under different forms by Kuijken & Gilmore (1989a, 1989b), Fuchs & Wielen (1992), Flynn & Fuchs (1994), and Holmberg & Flynn (2000), is related to the P–J theory. Near R = R_⊙, we can write

$$\begin{eqnarray}&&f({\boldsymbol{x}},{\boldsymbol{v}})\propto {f}_{z}(z,{v}_{z}).\end{eqnarray} \tag{ 29 }$$

where f(z, v_z) represents the one-dimensional phase space density in z and v_z. We therefore work solely with f_z(z, v_z) and drop the subscript z. This function satisfies

$$\begin{eqnarray}&&\int {{dv}}_{z}\;f(z,{v}_{z})\;g({v}_{z})=\rho (z)\langle g\rangle .\end{eqnarray} \tag{ 30 }$$

The Boltzmann equation also tells us that

$$\begin{eqnarray}&&{v}_{z}\displaystyle \frac{\partial f}{\partial z}-\displaystyle \frac{\partial {\rm{\Phi }}}{\partial z}\displaystyle \frac{\partial f}{\partial {v}_{z}}=0\end{eqnarray} \tag{ 31 }$$

which has solution

$$\begin{eqnarray}&&f(z,{v}_{z})=F\left(\displaystyle \frac{1}{2}{v}_{z}^{2}+{\rm{\Phi }}(z)\right).\end{eqnarray} \tag{ 32 }$$

That is, since this function is constant it must be a function solely of the "vertical energy" $\tfrac{1}{2}{v}_{z}^{2}+{\rm{\Phi }}(z)$. Using Equations (31) and (32) we can write

$$\begin{eqnarray}&&{\rho }_{i}(z)={\int }_{-\infty }^{\infty }{dw}\;{f}_{i}(z,w)\end{eqnarray} \tag{ 33a }$$

$$\begin{eqnarray}&&=\;{\int }_{-\infty }^{\infty }{dw}\;{f}_{i}(0,\sqrt{{w}^{2}+2{\rm{\Phi }}(z)})\end{eqnarray} \tag{ 33b }$$

and since we can write

$$\begin{eqnarray}&&{f}_{i}(z,w)={\rho }_{i}(z){f}_{z;i}(w)\end{eqnarray} \tag{ 34 }$$

where f_z(w) is the vertical velocity probability distribution at height z, normalized so that $\int \;{{dw}{f}}_{z}(w)=1$, we have

$$\begin{eqnarray}&&\displaystyle \frac{{\rho }_{{f}_{i}}(z)}{{\rho }_{i}(0)}={\int }_{-\infty }^{\infty }{dw}\;{f}_{z=0;i}(\sqrt{{w}^{2}+2{\rm{\Phi }}(z)}).\end{eqnarray} \tag{ 35 }$$

By extracting the in-the-plane velocity distribution ${f}_{z=0}(w)$ from tracer populations of A and F-stars, and integrating to find ${\rho }_{{\rm{A,F}}}(z)$ for specific potentials Φ(z), Holmberg & Flynn were able to test models for Φ by comparing the ${\rho }_{{\rm{A,F}}}(z)$ obtained in this fashion to the observed tracer densities. The model they considered was a modified version of the "Bahcall models" (Bahcall 1984a, 1984b, 1984c). The Bahcall models consist of splitting the visible matter into a series of isothermal components, each with a distinct vertical dispersion, σ_i. By assuming dispersions and densities in the plane, ρ_i(0), for all components, the Poisson–Jeans Equation (17) then gives a unique solution for Φ(z). The Bahcall model used by HF2000 is shown in Table 1 of the latter. An updated version of this model, featuring slight modifications, was used in Flynn et al. (2006). Holmberg & Flynn found that a mass model with little or no dark matter in the disk was in good agreement with the data, as Kuijken & Gilmore (1989a, 1989b) did previously using a similar method. By adding or subtracting invisible mass to the various components in the model, HF then obtained a range on the acceptable mass models, which gave a range of acceptable densities as a function of height. Figure 3 (top) compares the tracer densities with those predicted using the HF technique with an updated mass model containing no dark disk. We also show in Figure 3 (bottom) the same result with a dark disk of surface density ${{\rm{\Sigma }}}_{D}=10\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$ and h_D = 10 pc. Clearly it is of interest to know what the maximum value of Σ_D is compatible with data for a given h_D, within statistical errors.

Zoom In Zoom Out Reset image size

We now describe the numerical procedures used in the selection of our sample, as well as our methods for performing extinction corrections and for extracting the tracers' velocity distribution.

3.1. Sample Selection

We worked with the new reduction of the Hipparcos data (van Leeuwen 2008). For our sample we used mostly the same stars as did HF2000, which contained A and F-stars. We performed the same cuts as HF2000. That is, for A-type stars, $-0.2\lt B-V\lt 0.6$, and 0.0 < M_V < 1.0, and the completeness limit $V\leqslant 7.9+1.1\mathrm{sin}\;| b|$. By the completeness limit we have, using the definition of absolute magnitude,

$$\begin{eqnarray}&&{M}_{V}+5\;{\mathrm{log}}_{10}\left(\displaystyle \frac{d}{10\;{\rm{pc}}}\right)\leqslant 7.9+1.1\mathrm{sin}\;| b| \end{eqnarray} \tag{ 36 }$$

and so

$$\begin{eqnarray}&&\displaystyle \frac{d}{10\;{\rm{pc}}}\leqslant {10}^{(7.9+1.1\mathrm{sin}| b| -{M}_{V})/5}\end{eqnarray} \tag{ 37 }$$

and since $| b| \geqslant 0$ and M_V < 1.0, we have for our A-star sample, d ≲ 240 pc as our completeness limit. We therefore limited our sample to a vertical cylinder centered at the Sun with a radius and half-height of r_c = h_c = 170 pc. This ensures that the diagonals of the cylinder will be shorter than 240 pc. Note that HF2000 used a cylinder with radius 200 pc. The part of their A-star sample higher than 130 pc was therefore not complete. (By the same reasoning, their F-star sample was only complete to a height of 66 pc.) Using this method, our data set contained approximately 1500 A-stars. The precise numbers depend on the reddening corrections, as will be explained below.

For the F-star sample, the cuts in color were also $-0.2\lt B-V\lt 0.6$, but here we used 1.0 < M_V < 2.5. Here the maximum distance was found to be 120 pc. We therefore restricted our analysis to a cylinder of radius 50 pc and half-height of 109 pc. Additionally, to avoid any bias from the Coma Star Cluster (Trumpler 1938), we further cut off this sample from above at +40 pc. The final sample contained only about 500 stars, with the exact number depending on the extinction/reddening corrections and the value of the height of the Sun above the Galactic plane, Z_⊙.

3.2. Extinction Corrections

3.3. The Velocity Distributions

For the vertical velocity of their stars, HF2000 used the approximation

$$\begin{eqnarray}&&w=\displaystyle \frac{\kappa {\mu }_{b}}{\tilde{\pi }\;\mathrm{cos}b}+u\mathrm{cos}l\mathrm{tan}b+v\mathrm{sin}l\mathrm{tan}b.\end{eqnarray} \tag{ 38 }$$

Where $\tilde{\pi }$ is the parallax and κ = 4.74047 km s⁻¹ mas (mas yr)⁻¹. This approximation is valid for stars with $\mathrm{sin}b\ll 1$. However, we make use of the exact equation

$$\begin{eqnarray}&&w-{w}_{0}=\displaystyle \frac{\kappa {\mu }_{b}}{\tilde{\pi }}\mathrm{cos}b+{V}_{R}\mathrm{sin}b\end{eqnarray} \tag{ 39 }$$

where w₀ = 7.25 km s⁻¹ is the vertical velocity of the Sun (Schönrich et al. 2010) relative to the LSR. Since radial velocities are not known, we simply take an average:

$$\begin{eqnarray}&&\langle w\rangle ={w}_{0}+\displaystyle \frac{\kappa {\mu }_{b}}{\tilde{\pi }}\mathrm{cos}b+\langle {V}_{R}\rangle \mathrm{sin}b,\end{eqnarray} \tag{ 40 }$$

where we take

$$\begin{eqnarray}\langle {V}_{R}\rangle & = & {-{\boldsymbol{V}}}_{{\rm{Sun}}}\cdot \hat{{\boldsymbol{r}}}=-{u}_{0}\mathrm{cos}l\mathrm{cos}b\\ & & -{v}_{0}\mathrm{sin}l\mathrm{cos}b-{w}_{0}\mathrm{sin}b\end{eqnarray} \tag{ 41 }$$

with u₀ = 11.1 km s⁻¹ and v₀ = 12.24 km s⁻¹, also given by Schönrich et al. (2010). This follows from the fact that from the LSR frame we expect $\langle {V}_{R}^{{\rm{(LSR)}}}\rangle =0$ for any axisymmetric stellar population. Since we are also working with stars close to the plane ($\mathrm{sin}b\ll 1$), approximating V_R by its average still gives a good approximation, and this approximation will still be better than Equation (38). This gives

$$\begin{eqnarray}\langle w\rangle & = & \displaystyle \frac{\kappa {\mu }_{b}}{\tilde{\pi }}\mathrm{cos}b+{w}_{0}{\mathrm{cos}}^{2}b-{u}_{0}\mathrm{cos}l\mathrm{cos}b\mathrm{sin}b\\ & & -{v}_{0}\mathrm{sin}l\mathrm{cos}b\mathrm{sin}b.\end{eqnarray} \tag{ 42 }$$

We can do even better since the radial velocities have been measured for many stars in the solar region since 2000. Therefore, rather than simply use Equation (42) for all our stars, we supplemented the Hipparcos data with radial velocities from Barbier-Brossat & Figon (2000). We could thereby use measured radial velocities for 52% of our stars within $| b| \lt 12^\circ$. It was reported in Binney et al. (1997) that stars with measured radial velocities constitute a kinematically biased sample. However, Korchagin et al. (2003) showed, for their red giant sample with measured V_R, the result only appears biased when judging from their proper motions; once these are converted to Cartesian x, y, and z-velocities, the kinematic bias is no longer observed. Moreover, we did not restrict our analysis to stars with measured radial velocities, since we completed the sample using the average $\langle {V}_{R}\rangle$. In light of these factors, along with the fact that we restrict our analysis to low b, we expect the radial velocities to introduce very little bias to our results.

3.4. The Height of the Sun

As previously stated, one of the main differences between our analysis and that of HF2000 is that we update the mass model for the Galactic disk in light of more recent observations. One important feature of this update is the use of midplane densities for the interstellar gas as inputs to the Poisson–Jeans equation. These were not available at the time of the HF2000 study. The gas surface densities can also be used to place further constraints on the dark disk by demanding self-consistency of the model. This is explained in a separate paper (Kramer & Randall 2016) where we show how the combination of midplane and surface densities, which depend on the density profile in the vertical direction, can also set a bound on a dark matter disk. For now, we will discuss Holmberg & Flynn's model and how we modify it.

The Galactic disk model we update is Flynn et al.'s slightly updated model from their 2006 study, which includes fainter stars than HF2000 and allows for the presence of a thick stellar disk. It is shown in Table 2. In the far right column we show the updated values, which we discuss in Sections 4.1 and 4.2.

Table 2.  The Bahcall Model Used by Flynn et al. (2006)

i	Description	ρi(0)	σi	Σi	${\rho }_{i}{(0)}_{{\rm{new}}}$
		(M⊙ pc−3)	(km s−2)	(M⊙ pc−2)	(M⊙ pc−3)
1	H2	0.021	4.0	3.0	0.014a
2	HI(1)	0.016	7.0	4.1	0.015a
3	HI(2)	0.012	9.0	4.1	0.005a
4	warm gas	0.0009	40.0	2.0	0.0011a
5	giants	0.0006	20.0	0.4	0.0006a
6	MV < 2.5	0.0031	7.5	0.9	0.0018
7	$2.5\lt {M}_{V}\lt 3.0$	0.0015	10.5	0.6	0a
8	3.0 < MV < 4.0	0.0020	14.0	1.1	0.0018a
9	4.0 < MV < 5.0	0.0022	18.0	1.7	0.0029
10	5.0 < MV < 8.0	0.007	18.5	5.7	0.0072
11	MV > 8.0	0.0135	18.5	10.9	0.0216
12	white dwarfs	0.006	20.0	5.4	0.0056
13	brown dwarfs	0.002	20.0	1.8	0.0015
14	thick disk	0.0035	37.0	7.0	0.0035
15	stellar halo	0.0001	100.0	0.6	0.0001

Note. The Σ_i were calculated by Flynn et al. from the solution to the Poisson–Jeans equation, except in the case of the interstellar gas, where they were held fixed by HF and the midplane densities were chosen to give the correct Σ_i. Note that revised values of ${\rho }_{i}{(0)}_{{\rm{new}}}$ give revised values of Σ_i and that these are dependent on Σ_D.

^aMarks components whose dispersions σ_i have been revised.

Download table as:  ASCIITypeset image

4.1. Interstellar Gas

4.2. Other Components

4.3. Poisson–Jeans Solver

We implement our Poisson–Jeans solver with the densities and dispersions as inputs, as did HF2000. Our solver is implemented by using Equation (18) in the following way. An initial distribution is assumed for each component ${\rho }_{i}^{(0)}(z)$. The potential Φ(z) due to these components is then calculated via

$$\begin{eqnarray}&&{{\rm{\Phi }}}^{(0)}(z)=4\pi G\displaystyle \sum _{i}{\int }_{0}^{z}{dz}^{\prime} \;{\int }_{0}^{{z}^{\prime }}{dz}^{\prime\prime} \;{\rho }_{i}^{(0)}(z^{\prime\prime} )\end{eqnarray} \tag{ 43 }$$

which then gives the next iteration of ρ_i(z):

$$\begin{eqnarray*}&&{\rho }_{i}^{(N+1)}={\rho }_{i}(0)\mathrm{exp}(-{{\rm{\Phi }}}^{(N)}(z)/{\sigma }_{i}^{2})\end{eqnarray*}$$

and this process is repeated until the solution converges. Remarkably, in only 5 iterations the solution for a single component disk converges to the Spitzer (1942) solution (Equation (21)) to better than one part in 10⁷! This is affected only very slightly by the initial distributions ${\rho }_{i}^{(0)}$ assumed.

5.1. Traditional Method—Statistics

In order to compare the stellar kinematics to a given dark disk model, we define a χ²-type statistic X that measures the distance between the predicted and observed densities:

$$\begin{eqnarray}&&X[{\rm{\Phi }}]\equiv {\int }_{{z}_{{\rm{min}}}}^{{z}_{{\rm{max}}}}{dz}\;\displaystyle \frac{| {\rho }_{{\rm{obs}}}(z)-{\rho }_{f,{\rm{obs}}}[{\rm{\Phi }}(z)]{| }^{2}}{{{\rm{\Delta }}}_{\rho }^{2}(z)}\end{eqnarray} \tag{ 44 }$$

where

$$\begin{eqnarray}&&{\rho }_{f,{\rm{obs}}}[{\rm{\Phi }}(z)]\equiv {\rho }_{f}(z)=\rho (0){\int }_{-\infty }^{\infty }{dw}\;{f}_{z=0}(\sqrt{{w}^{2}+2{\rm{\Phi }}(z)})\end{eqnarray} \tag{ 45 }$$

as defined in (35), Δ_ρ(z) represents the uncertainty in ρ(z) at the position z, and where z_min and z_max are given by the completeness limits, Z_⊙ ± 170 pc for A stars and ${Z}_{\odot }-92\;{\rm{pc}}$, ${Z}_{\odot }+40\;{\rm{pc}}$ for F stars.

As explained in Appendix A, we can then bootstrap from the data to obtain the probability distribution P(X) and obtain total probabilities on the A and F-star data $p({X}_{{\rm{A}}},{X}_{{\rm{F}}}| {\rm{\Phi }})$ given a dark matter disk model contained in the gravitational potential Φ.

5.2. A Cross-check: Measuring ${{\rm{\Phi }}}_{{\rm{obs}}}(z)$ Directly

It is worth pointing out that besides the traditional HF analysis, which compares measured and model densities, it is also possible to use the HF equation to measure the Galactic potential ${{\rm{\Phi }}}_{{\rm{obs}}}(z)$ directly from the data. This has the advantage that it treats the errors in density and velocity on more equal footing. This can be done by interpreting Equation (35) as an equation for ${\rho }_{{\rm{A,F}}}({\rm{\Phi }})$ (as in Kuijken & Gilmore 1989b, Equation (20)):

$$\begin{eqnarray}&&\displaystyle \frac{{\rho }_{f}({\rm{\Phi }})}{\rho (0)}={\int }_{-\infty }^{\infty }{{dw}}_{0}\;{f}_{z=0}(\sqrt{{w}_{0}^{2}+2{\rm{\Phi }}})\end{eqnarray} \tag{ 46 }$$

and then inverting this to obtain Φ_f(ρ). The "observed" gravitational potential is then given by combining this with the observed density ${\rho }_{{\rm{A,F}}}(z)$:

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{{\rm{obs}}}(z)\equiv {{\rm{\Phi }}}_{f}({\rho }_{{\rm{obs}}}(z))\end{eqnarray} \tag{ 47 }$$

Note that Equation (46) always gives ρ(Φ)/ρ(0) ≤ 1 for real w. For values of the observed density ρ_obs(z) > ρ(0), we analytically continued ${f}_{z=0}(w)$ to imaginary w (negative Φ) by fitting ${f}_{z=0}(w)$ to a sum of three analytic Gaussians. Once we obtain ${{\rm{\Phi }}}_{{\rm{obs}}}(z)$, we can then perform a cross-check on our analysis by comparing the measured ${{\rm{\Phi }}}_{{\rm{obs}}}(z)$ with our model Φ(z) directly. The statistics can be accounted for similarly as in Appendix A, by sampling ${{\rm{\Phi }}}_{{\rm{obs}}}(z)$ and computing their distribution.

5.3. Non-equilibrium Method

As was mentioned in Section 1, our tracer populations show evidence for deviations from equilibrium behavior. This evidence is a non-zero mean velocity and a displacement from the assumed position of the Galactic plane. Figure 4 shows the velocity and density distributions of the A stars. They feature a net velocity of 1.3 ± 0.3 km s⁻¹ and a net displacement from the Galactic plane of 19 ± 5 pc. A bulk velocity and vertical displacement from the plane can clearly be seen.

Zoom In Zoom Out Reset image size

Because of these features, the assumption ∂f/∂t = 0 in Equation (4) will not be satisfied. Consequently, Equations (31) and (32) no longer hold, and the HF relation (35) that, in principle, constrains the potential will no longer be useful for constraining the Galactic model. On the other hand, we might expect that, if the tracer distribution is oscillating about the plane in the Galactic potential, the long-time average of ∂f/∂t will vanish. We therefore expect that the HF relation will hold for long-time averages. However, the only way to compute the long-time average of the star dynamics is to assume some form for the potential Φ(z). We show in Appendix D that, in this case, the HF relation (Equation (35)) is trivially satisfied for the potential Φ(z):

$$\begin{eqnarray}&&{\bar{\rho {(0)}^{-1}}}_{[{\rm{\Phi }}]}{\overline{\rho (z)}}_{[{\rm{\Phi }}]}=\int {dw}\;{\bar{{f}_{z=0}}}_{[{\rm{\Phi }}]}(\sqrt{{w}^{2}+2{\rm{\Phi }}(z)})\end{eqnarray} \tag{ 48 }$$

where ${\bar{\cdot }}_{[{\rm{\Phi }}]}$ represents the time average under evolution in the potential Φ(z). Since Equation (48), the correct statement of the HF relation, is trivially satisfied for any potential it therefore cannot be used to constrain a dark matter model. Although we do indeed expect a dark disk to "pinch" the tracer distributions as explained in Section 2.2, this pinch cannot be measured with the HF method if the sample is oscillating.

The solution we propose is to observe over time the shape of the distribution relative to the position of its center z₀(t), i.e., $\rho (z-{z}_{0}(t),t)$, and to predict its time average, $\overline{\rho (z-{z}_{0})}$. Figure 5 shows this time average for a model with Σ_D = 0 as well as with ${{\rm{\Sigma }}}_{D}=20\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$. Comparing this to the Hipparcos "snapshot," which, in our case, is used as initial data, gives a method of constraining Σ_D by assuming that the instantaneous distribution should not deviate strongly from its time average. We therefore define an appropriate χ² parameter:

$$\begin{eqnarray}&&{\chi }^{2}=\int {dz}\;{\left|\displaystyle \frac{\rho (z-{z}_{0}(0),0)-\overline{\rho (z-{z}_{0})}}{{{\rm{\Delta }}}_{\rho }(z)}\right|}^{2}.\end{eqnarray} \tag{ 49 }$$

Where ${{\rm{\Delta }}}_{\rho }^{2}(z)$ represents the expected variance in $\rho (z-{z}_{0})$. For this, we use the variance ${{\rm{\Delta }}}_{\rho }^{2}(z)$ computed in Section 5.1. This assumes that the distribution $\rho (z-{z}_{0})$ is static in time, and that at any time, including the present t = 0, we expect ${\rho }_{{\rm{obs}}}(z-{z}_{0})$ to be equal to its average up to sampling error Δ_ρ(z). Although a more careful analysis would include the aforementioned sampling variance ${{\rm{\Delta }}}_{\rho }^{2}$ as well as the variance due to the time dependence of $\rho (z-{z}_{0})$, and will therefore be larger than ${{\rm{\Delta }}}_{\rho }^{2}$, for the purposes of this analysis we neglected this contribution and considered only the sampling error ${{\rm{\Delta }}}_{\rho }^{2}$. More careful determinations may be important in the future.

Zoom In Zoom Out Reset image size

An alternate proposal for how the present distribution $\rho (z-{z}_{0})$ fits the model is to measure the stability of the distribution over time. Were we to observe that the initial (possibly oscillating) distribution $\rho (z-{z}_{0})$ decayed to a different (possibly oscillating) distribution, we would conclude that there was a dynamical mismatch between the distribution and the potential. On the other hand, were we to observe that the initial distribution retained its behavior over time, we would conclude that the distribution and potential were appropriately matched. Such an analysis would also have to account for spiral arm crossing and is beyond the scope of this manuscript. In any case, if evidence persists of nonstatic distributions, more careful analyses using nonstatic distributions will be required.

Static Method

We assigned probabilities to models with various Σ_D and h_D in the manner described in Section 5.1. The scale height, h_D, was defined so that the value of the density of the dark disk at z = h_D is ρ(z = 0) sech²(1/2).

The probabilities and probability density derived without applying reddening corrections are shown in Figure 6. For low scale heights, h_D, the probability has a 95% upper bound at ${{\rm{\Sigma }}}_{D}\simeq 2\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$ and a 95% lower bound of $1\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$. For larger scale heights the bound is slightly weaker, with the 95% upper bound growing beyond ${{\rm{\Sigma }}}_{D}\simeq 5\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$ at h_D = 80 pc. However, without reddening corrections, the entire parameter space appears to be disfavored at 94%. Moreover, without reddening corrections, the value Σ_D = 0 appears to be excluded at greater than 95% confidence.

Zoom In Zoom Out Reset image size

Figure 7, on the other hand, shows, on the left, the results using the standard HF method used in Figure 6 but with reddening corrections applied. This plot shows the results determined using the old values for the mass parameters from Flynn et al. (2006). We see that the upper bound is roughly the same with as without reddening corrections. Now, however, the entire parameter space is only disfavored at 77%. On the right, we have the bounds using the new values for the gas parameters (Kramer & Randall 2016) for comparison. The 95% upper bound in this case (for low scale height) is around $4\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$. Figure 8, on the left, shows the values determined by using both the updated gas values and the updated values for the stellar components (McKee et al. 2015). Here, the bound is $3\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$, which is remarkably similar to the bound using the 2006 values. This is because changes the changes in the various mass parameters tend to compensate each other on average, as would be expected statistically.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In terms of midplane densities, the dark matter densities are less than $0.02\;{M}_{\odot }\;{{\rm{pc}}}^{-3}$ for thick scale heights (h_D ≳ 75 pc) but are unbounded for low scale heights, varying as ${h}_{D}^{-1}$. Note that, according to this analysis, the value Σ_D = 0 is still ruled out at more than 70% confidence, so one might still question the robustness of this method as applied to this data set and mass parameters.

In order to know whether the uncertainty in the interstellar gas mass parameters have an important effect on the bound, we compute the bounds using gas densities one standard deviation below their average values. The results are shown in Figure 8 on the right. Here, the 95% upper bound at low scale height (h_D ∼ 10 pc) is $3.5\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$, which is slightly higher than the bound using the mean gas values. The value Σ_D = 0 is still ruled out at more than 70% confidence.

The plot in Figure 9 shows the results of the cross-check mentioned in Section 5.2 where the HF equation was interpreted as an equation for ${\rm{\Phi }}(\rho )$. The results of this method are consistent with the results from the standard analysis, although the bounds are weaker. This is because Φ(ρ) blows up as $\rho \to {0}^{+}$, causing large fluctuations in Φ(z) at high z and thus reducing the sensitivity of this method relative to the standard method.

Zoom In Zoom Out Reset image size

Non-equilibrium Constraint

Figure 10 shows the results obtained using the non-equilibrium HF method described in Section 5.3 and Appendix D, applied to the A star sample of Section 3.1. Here, we assume the distribution moves in the gravitational potential determined by the Galactic model and ask that the shape of the measured distribution not be far from that of the average distribution. The top plot shows the results computed assuming the low value of Z_⊙ = 7 pc (cf. Section 3.4), i.e., assuming our distribution is centered at the Galactic plane. The bottom plot shows the results computed assuming the more accepted value of Z_⊙ = 26 pc. Note that the upper subplots in both figures show relative probability density, which, for Z_⊙ = 26 pc, very slightly seems to favor ${{\rm{\Sigma }}}_{D}\simeq 2\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$. In any case, the absolute probabilities clearly show that any values between Σ_D = 0 and quite high density are allowed.

Zoom In Zoom Out Reset image size

The bound here is significantly weaker that of the static method. One reason is that the oscillations reduce the amount of time that the tracers spend in the dark disk, thus reducing their sensitivity to it. Applying the static method to such an oscillating population would not yield a bound on the parameter space since, in this case, as explained in Section 5.3, the static HF relation is trivially satisfied. Another reason the bound here is weaker simply depends on the amount of usable data in this method. This is reduced because the stars in the data are oscillating vertically, but, in order to take a time average of the stars' distribution, we can include only values of z that are covered by the data at all times. In other words, the oscillation of the stars means that the z-cutoffs on the data are oscillating as well. The z-cutoffs on the time average are therefore the minimum values of these oscillating z-cutoffs. It is for this reason that the non-equilibrium analysis was performed only on the A stars: the amount of F star data available (c.f. Section 3.1) after taking into account the oscillating z-cutoffs was not sufficient for analysis.

For low scale heights (h_D = 10 pc), in the non-equilibrium method, we find that the 95% upper bound is between ${{\rm{\Sigma }}}_{D}\simeq 10\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$ and ${{\rm{\Sigma }}}_{D}\simeq 14\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$, and that the bound grows with h_D, as for the static case.²

Galactic Disk Parameters

Another important result of our analysis is the value for the surface density of the Galactic disk, which we compare against existing measurements. Figure 11 shows the values of Σ_ISM, the surface density of the interstellar medium; Σ_*, the surface density of the stellar disk (not including brown dwarfs and stellar remnants); and Σ_1.1, the total surface density to 1.1 kpc, including all visible and dark components (including e.g., brown dwarfs, dark halo, and dark disk). These were computed for a dark disk model with scale height 10 pc, and were obtained by integrating the self-consistent Poisson–Jeans solutions of each isothermal component. For the value of ${{\rm{\Sigma }}}_{D}=10\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$ at this scale height (black dashed line in Figure 11), we find ${{\rm{\Sigma }}}_{{\rm{ISM}}}\simeq 10\pm 2\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$, the uncertainty being attributed to uncertainty in gas midplane densities, consistent with the value ${{\rm{\Sigma }}}_{{\rm{ISM}}}=11.0\pm 0.8\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$ of Kramer & Randall (2016) and with the value ${{\rm{\Sigma }}}_{{\rm{ISM}}}=12.8\pm 1.5\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$ of McKee et al. (2015). At the Σ_D upper bound obtained from the non-equilibrium method (gray dashed line), we have ${{\rm{\Sigma }}}_{{\rm{ISM}}}=9.7\pm 0.7\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$. Values of Σ_D higher than about ${{\rm{\Sigma }}}_{D}=20\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$ give lower values of Σ_ISM that are at odds with the literature. In Kramer & Randall (2016), we perform a more detailed analysis of the self-consistency of the interstellar gas parameters and derive an independent bound on the DDDM parameter space.

Zoom In Zoom Out Reset image size

We also find ${{\rm{\Sigma }}}_{* }=30\pm 2\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$ for the stellar surface density (not including brown dwarfs and stellar remnants) for ${{\rm{\Sigma }}}_{D}=10\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$, h_D = 10 pc. This agrees with the value ${{\rm{\Sigma }}}_{* }=30\pm 1\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$ measured by Bovy et al. (2012). (On the other hand, the latter measurement was made assuming exponential profiles for the stellar disk components. Subleading corrections to this value obtained by assuming more realistic disk profiles near the plane give anywhere between 26 and $28\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$ depending on the specific model.) Note that our value was computed using the recent values of McKee et al. (2015), and is somewhat higher than the values inferred using Flynn et al.'s (2006) model. The latter's parameters yield ${{\rm{\Sigma }}}_{* }=30\pm 2\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$ without the dark disk, and ${{\rm{\Sigma }}}_{* }=26\pm 2\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$ when it is included. The uncertainties on our value of Σ_* were estimated from those on the individual components given in McKee et al. At the upper bound, we find ${{\rm{\Sigma }}}_{* }=28.5\pm 2.0\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$. Note that Σ_D = 0 gives ${{\rm{\Sigma }}}_{* }=34\pm 2\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$, which is too large. The Σ_* values therefore seem to favor ${{\rm{\Sigma }}}_{D}\gtrsim 4\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$.

For the total surface density to 1.1 kpc (still with h_D = 10 pc), we find Σ_1.1 grows between $73.5\pm 6.0\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$ for Σ_D = 0 and ${{\rm{\Sigma }}}_{D}=76.5\pm 6.0\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$ for ${{\rm{\Sigma }}}_{D}=10\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$. These are slightly higher than the values in Bovy & Rix (2013) who measured $68\pm 4\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$, but agree within their combined uncertainty. The major source of uncertainty in our measurement of this quantity is that of the dark halo density, which we set as described in Section 4.2. Here, too, we find that ${{\rm{\Sigma }}}_{D}\gtrsim 20\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$ is at odds with the literature. For thicker scale heights (h_D ≳ 50 pc), we find that even for ${{\rm{\Sigma }}}_{D}\gtrsim 13\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$, Σ_1.1 is too large.

It is of interest to use existing and future kinematical data to ascertain the possible existence of a dark disk component in the Milky Way disk. Previous analyses argued that a dark disk is not necessary to match the data but, as has been often demonstrated, that is a far cry from ruling it out. Inspired by the Holmberg & Flynn (2000) study, we have rederived the kinematic constraints on a dark disk by considering a dark disk in a self-consistent manner. Our analysis features updated kinematics and extinction corrections, an updated model for the interstellar gas, and careful statistics.

We find that for a dark disk of sech² scale height of 10 pc, the static method rules out surface densities greater than $3\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$ at 95% confidence. In terms of midplane density, the favored value is ${\rho }_{D}={0.0}^{+0.1}\;{M}_{\odot }\;{{\rm{pc}}}^{-3}$, giving a total matter density in the plane of $\rho ={0.1}^{+0.1}\;{M}_{\odot }\;{{\rm{pc}}}^{-3}$. These bounds increase with scale height. We note, however, that in the static method, even a model with Σ_D = 0 is disfavored at around 70%, pointing to the inadequacy of the method for the current data set. Using values of gas midplane densities a standard deviation lower than their average moves this bound closer to $4\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$. We also find that the updated values of the mass model increase the bound relative to the old Flynn et al. (2006) values by about 50% from $2\pm 5\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$ to $3\pm 1\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$.

These results were derived by arbitrarily removing the non-equilibrium features from our sample, which consisted of a bulk vertical velocity and a net vertical displacement from the Galactic midplane. However, if these features are instead taken into account using a non-equilibrium version of the HF analysis, the bound increases to ${{\rm{\Sigma }}}_{D}\leqslant 14\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$, assuming Z_⊙ = 26 pc. This is partly because, in this method, it becomes more difficult to constrain a dark disk for the same amount of data.

We also showed that the static HF method can be modified to directly measure the Galactic potential.

For thin dark disk models, the total surface density of the Galactic disk to 1.1 kpc, Σ_1.1, as well as the surface densities of visible matter, Σ_* and Σ_ISM, were found to be consistent with literature values both at the target values of Σ_D and h_D, and at their 95% upper bound. According to our Poisson–Jeans solver, literature Σ_* values appear to favor a dark disk model.

Since the non-equilibrium analysis allows surface densities as low as zero and up to $14\;{M}_{\odot }\;{{\rm{pc}}}^{-2}$, it follows that a dark disk may account for the comet periodicity that was extrapolated from the crater measurements. The results using the Gaia data promise to be more constraining, as more data close to the Galactic midplane will be available over a much larger area. However, it will be important to verify whether the sample used is in equilibrium. The correct method to use may indeed have to take account of the motion of the star populations, in which case non-equilibrium methods will be more appropriate.

We would like to thank Chris Flynn for all his comments and suggestions, and for reviewing our results. We would also like to thank Johann Holmberg for his comments. We would like to thank Jo Bovy for sharing his insights into the DDDM model and its constraints, and for pointing us toward the Holmberg & Flynn technique for setting a bound on it, as well as for reviewing our results. We would also like to thank Chris McKee for sharing his most recent results with us and for reviewing our work. We thank our referee for useful suggestions and comments. We would like to thank Matt Reece and Doug Finkbeiner for discussions of the statistical analysis. Thanks also to Alexander Tielens, Katia Ferriere, and Matt Walker for help on the ISM parameters, and to Chris Stubbs for his interest and comments. E.D.K. was supported by N.S.F. grants of L.R., and by Harvard F.A.S., Harvard Department of Physics, and Center for the Fundamental Laws of Nature. L.R. was supported by N.S.F grants PHY-0855591 and PHY-1216270. Calculations were performed using MATLAB 2015a.

In order to compare the stellar kinematics to a given dark disk model, we define a χ²-type statistic X that measures the distance between the predicted and observed densities:

$$\begin{eqnarray}&&X[{\rm{\Phi }}]\equiv {\int }_{{z}_{{\rm{min}}}}^{{z}_{{\rm{max}}}}{dz}\;\displaystyle \frac{| {\rho }_{{\rm{obs}}}(z)-{\rho }_{f,{\rm{obs}}}[{\rm{\Phi }}(z)]{| }^{2}}{{{\rm{\Delta }}}_{\rho }^{2}(z)}\end{eqnarray} \tag{ 50 }$$

where

$$\begin{eqnarray}&&{\rho }_{f,{\rm{obs}}}[{\rm{\Phi }}(z)]\equiv {\rho }_{f}(z)=\rho (0){\int }_{-\infty }^{\infty }{dw}\;{f}_{z=0}(\sqrt{{w}^{2}+2{\rm{\Phi }}(z)})\end{eqnarray} \tag{ 51 }$$

as defined in (35), Δ_ρ(z) represents the uncertainty in ρ(z) at the position z, and where where z_min, z_max are given by the completeness limits, Z_⊙ ± 170 pc for A stars and ${Z}_{\odot }-92\;{\rm{pc}}$, ${Z}_{\odot }+40\;{\rm{pc}}$ for F stars. As explained in Section 3.1, for the F stars, we do not include any data higher than 40 pc above the Sun.

The observed densities ${\rho }_{{\rm{obs}}}(z)\equiv {\rho }_{{\rm{A,F}}}(z)$ were constructed using the kernel histogram technique. That is, we represented each star as a Gaussian with unit area in position and velocity space, and then obtained total distributions by summing these Gaussians. The uncertainties in position and velocity of the individual stars, Δz and Δw, do not, however, fully account for the error. The chief source of error on this sparse distribution is the likelihood that stars will actually fill in the purported distribution (i.e., Poisson error).

We estimate the latter in the following way. If, in the case of ρ(z), we assume the data are arranged in "bins" of half-width, Δz, then the density ρ(z) should be

$$\begin{eqnarray}&&\rho (z)=\displaystyle \frac{N(z)}{2A{\rm{\Delta }}z}\end{eqnarray} \tag{ 52 }$$

where N(z) is the number of stars in each "bin" and A is the cross-sectional area of the bin at height z. Both Δz and N(z) are unknown. However, we know that the fluctuations in ρ(z) should be given, according to Poisson statistics, by

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\rho }(z)=\displaystyle \frac{\sqrt{N(z)}}{2A{\rm{\Delta }}z}\end{eqnarray} \tag{ 53 }$$

Eliminating N(z), we have

$$\begin{eqnarray}&&{\rm{\Delta }}z=\displaystyle \frac{\rho }{2A{{\rm{\Delta }}}_{\rho }^{2}}.\end{eqnarray} \tag{ 54 }$$

Since the determination of ρ(z) and Δ_ρ(z) themselves depend on the value of Δz, Equation (54) should be used recursively. Also, since the derivation of Equation (54) was heuristic, we estimate this computation of Δz to be correct only to within a factor of two or so. In this way, we obtain values of Δz in the range 6–12 pc. To this we must add the Δz arising from measurement error. This error, which grows with z², is derived by propagating the error in parallax available in the Hipparcos catalog, and grows with z². Summing these two contributions in quadrature gives values of Δz between Δz ≃ 7 pc near z = 0 and Δz ≃ 20 pc near the extremities in the case of A stars, and Δz ≃ 7–13 pc for the F stars. On the other hand, standard kernel density estimation techniques (which assume a constant kernel width and unimodal distribution) suggest widths closer to (Silverman 1986):

$$\begin{eqnarray}&&{\rm{\Delta }}z\simeq \displaystyle \frac{{\rm{stdev}}\{{z}_{i}\}}{{N}^{1/5}}\simeq 18\;{\rm{pc}}\end{eqnarray} \tag{ 55 }$$

and

$$\begin{eqnarray}&&{\rm{\Delta }}w\simeq \displaystyle \frac{{\rm{stdev}}\{{w}_{i}\}}{{N}^{1/5}}\simeq 2.1\;{\rm{pc}},\end{eqnarray} \tag{ 56 }$$

for A stars. We find similar values for the F stars. We find that the final result is roughly independent of this width for the range Δz ≃ 6–30 pc. The weakest bound is obtained for Δz ≃ 18 pc, before the width of the kernels becomes so large that it biases the distribution and begins to remove the signature of any potential structure near z = 0.

In order to construct the relevant uncertainty Δ_ρ(z) we considered the density distributions, ρ_obs(z), obtained by repeatedly sampling a fraction q < 1 of the stars in the respective samples. We did this by either including or not including each star in the original data set with probability q. The variance between these distributions thus gives the spread ${{\rm{\Delta }}}_{\rho }^{(q)}(z)$ at every point z. In terms of ${{\rm{\Delta }}}_{\rho }^{(q=1/2)}(z)$ obtained from repeatedly sampling half the stars in the sample, we can obtain the uncertainty Δ_ρ(z) on the original sample as (Equation (75a))

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\rho ,{\rm{obs}}}(z)=2\;{{\rm{\Delta }}}_{\rho }^{(q=1/2)}(z)\end{eqnarray} \tag{ 57 }$$

which can be derived using the binomial distribution with probability q (see Appendix B). In numerical simulations, we find this factor of 2 is closer to 1.97. Although we do not know the source of this discrepancy, we note that its effect is small compared to the remaining errors in the analysis. The errors ${{\rm{\Delta }}}_{{\rho }_{f,{\rm{obs}}}}$ were constructed in a similar way by sampling velocities from the data set, giving (Equation (81))

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{{\rho }_{f},{\rm{obs}}}(z)={{\rm{\Delta }}}_{{\rho }_{f}}^{(q=1/2)}(z).\end{eqnarray} \tag{ 58 }$$

The errors were then added in quadrature:

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\rho }^{2}(z)={{\rm{\Delta }}}_{\rho ,{\rm{obs}}}^{2}(z)+{{\rm{\Delta }}}_{{\rho }_{f,{\rm{obs}}}}^{2}(z).\end{eqnarray} \tag{ 59 }$$

Uncertainties in the Galactic potential, obtained from Sections 4.1 in the case of interstellar gas and from McKee et al. in the case of stellar components, were were found to be negligible compared to the errors above.

In order to associate a value of X[Φ] with a probability, we constructed a probability distribution for the values of X[Φ] obtained from fluctuations of the density and velocity distributions. To do this, we note that, given the true potential Φ_true of the Galaxy, the value

$$\begin{eqnarray}&&X[{{\rm{\Phi }}}_{{\rm{true}}}]=\int {dz}\displaystyle \frac{{| {\rho }_{{\rm{obs}}}-{\rho }_{f,{\rm{obs}}}[{{\rm{\Phi }}}_{{\rm{true}}}]| }^{2}}{{{\rm{\Delta }}}_{\rho }^{2}}\end{eqnarray} \tag{ 60 }$$

is itself a fluctuation with respect to its equilibrium value

$$\begin{eqnarray}&&{X}_{{\rm{eq}}}[{{\rm{\Phi }}}_{{\rm{true}}}]=\int {dz}\displaystyle \frac{| {\rho }_{{\rm{eq}}}-{\rho }_{f,{\rm{eq}}}[{{\rm{\Phi }}}_{{\rm{true}}}]{| }^{2}}{{{\rm{\Delta }}}_{\rho }^{2}}=0,\end{eqnarray} \tag{ 61 }$$

as ρ_obs and f_obs are assumed to be fluctuations of their equilibrium values ρ_eq and f_eq, respectively. We can similarly define Φ_obs so that

$$\begin{eqnarray}&&{\rho }_{f,{\rm{obs}}}[{{\rm{\Phi }}}_{{\rm{obs}}}]={\rho }_{{\rm{obs}}}\end{eqnarray} \tag{ 62 }$$

and X[Φ_obs] = 0. We then compute fluctuations in X[Φ_obs] by sampling fluctuations in ρ_obs and f_obs. That is, each time, k, that we sample a set ${\{z^{\prime} ,w^{\prime} \}}_{k}$ from the parent populations {z, w}, this gives an observed density ${\rho }_{{\rm{obs}},k}(z)$ and ${f}_{{\rm{obs}},k}(w)$. We therefore have

$$\begin{eqnarray}&&{X}_{k}[{{\rm{\Phi }}}_{{\rm{obs}}}]=\int {dz}\displaystyle \frac{{| {\rho }_{{\rm{obs}},k}-{\rho }_{f,{\rm{obs}},k}[{{\rm{\Phi }}}_{{\rm{obs}}}]| }^{2}}{{{\rm{\Delta }}}_{\rho }^{2}}.\end{eqnarray} \tag{ 63 }$$

Thus, by repeatedly sampling values of X_k, we obtain the distribution P(X). We can then use the distribution P(X) to associate probabilities with every point in parameter space by computing the cumulative X distribution:

$$\begin{eqnarray}&&C(X)={\int }_{0}^{X}{{dX}}^{\prime }\;P({X}^{\prime }).\end{eqnarray} \tag{ 64 }$$

C(X[Φ]) represents the probability that the fluctuations in ρ_obs, ${\rho }_{f,{\rm{obs}}}[{\rm{\Phi }}]$ with respect to the equilibrium distributions, ρ_eq and f_eq, assuming that Φ_true = Φ could result in a value X ≤ X[Φ]. The probability $1-C(X[{\rm{\Phi }}])$ therefore represents the the probability that X ≥ X[Φ]. In Bayesian terms, this is the probability $p(X| {\rm{\Phi }})$. Combining the independent results from the A and F stars, we have the probability

$$\begin{eqnarray}&&p({X}_{{\rm{A}}},{X}_{{\rm{F}}}| {\rm{\Phi }})=(1-C({X}_{{\rm{A}}}[{\rm{\Phi }}]))\times (1-C({X}_{{\rm{F}}}[{\rm{\Phi }}])).\end{eqnarray} \tag{ 65 }$$

We reject any model for which $p({X}_{{\rm{A}}},{X}_{{\rm{F}}}| {\rm{\Phi }})\lt 0.05$. This will be our criterion for "exclusion at 95% confidence." Note that this does not represent the probability of the model itself, but rather the probability of the results given the model.

As an aside, probability densities in model spaces Σ_D and h_D can also be computed by fitting the P(X) distribution to a χ² distribution, as explained in Appendix C. However, since this fit involves two slightly degenerate parameters, as explained therein, these probability densities are only approximate and are included in our results for illustration purposes only. We will see that, even so, they match qualitatively what is suggested by the absolute probabilities.

In Appendix A, we claimed that we could construct the errors on the total distributions Δ_ρ(z) or Δ_f(w) by repeatedly sampling a fraction q of the stars. We now proceed to derive the relationship between these values and the values ${{\rm{\Delta }}}_{\rho }^{(q)}(z)$, ${{\rm{\Delta }}}_{f}^{(q)}(w)$ obtained by sampling.

We model the positions and velocities {z_i, w_i} of the stars as being drawn from some "parent" probability distribution, f(z, w). If we sample a very small fraction q ≪ 1 of the stars, we expect the randomness in the resulting set $\{z{^{\prime} }_{{i}^{\prime }},w{^{\prime} }_{{i}^{\prime }}\}$ relative to the set z_i, w_i to mimic the randomness of the set {z_i, w_i} relative to the parent probability distribution f(z, w), with the the errors reduced by a Poisson factor $\sqrt{q}$. We can therefore estimate the randomness in the parent set {z_i, w_i} by repeatedly sampling a small fraction q from this set and dividing the standard deviation ${{\rm{\Delta }}}^{(q)}$ in the resulting set by $\sqrt{q}$:

$$\begin{eqnarray}&&{\rm{\Delta }}=\underset{q\to 0}{\mathrm{lim}}\displaystyle \frac{{{\rm{\Delta }}}^{(q)}}{\sqrt{q}}.\end{eqnarray} \tag{ 66 }$$

However, since the data set $\{{z}_{i},{w}_{i}\}$ is of finite size, for q ≪ 1 we cannot generate a large enough child set $\{z{^{\prime} }_{{i}^{\prime }},w{^{\prime} }_{{i}^{\prime }}\}$ to be able to obtain reliable statistics. On the other hand, if we use a larger sampling fraction q, Equation (66) will no longer hold because different samplings i' will repeatedly contain the same values. The trick will therefore be to relate ${\mathrm{lim}}_{q\to 0}{{\rm{\Delta }}}^{(q)}/\sqrt{q}$ to a reference value ${{\rm{\Delta }}}^{({q}_{0})}$ for some reference ${q}_{0}={ \mathcal O }(1)$.

In order to derive the dependence of ${{\rm{\Delta }}}^{(q)}$ on q, consider binning the data {z_i, w_i} into bins in z and w space. For concreteness, let us consider only z space for now. Let N_k be the number of stars in bin k. Each star in bin k has a probability q of being sampled, and a probability $1-q$ of not being sampled. The distribution of outcomes is therefore binomial:

$$\begin{eqnarray}&&1={(q+1-q)}^{{N}_{k}}\end{eqnarray} \tag{ 67a }$$

$$\begin{eqnarray}&&=\;\displaystyle \sum _{r=0}^{{N}_{k}}\displaystyle \left(\displaystyle \genfrac{}{}{0em}{}{{N}_{k}}{r}\right){q}^{r}{(1-q)}^{{N}_{k}-r},\end{eqnarray} \tag{ 67b }$$

in the sense that the probability of sampling r stars in the bin k is given by

$$\begin{eqnarray}&&{P}_{k,r}=\displaystyle \left(\displaystyle \genfrac{}{}{0em}{}{{N}_{k}}{r}\right){q}^{r}{(1-q)}^{{N}_{k}-r}\end{eqnarray} \tag{ 68a }$$

$$\begin{eqnarray}&&=\;{(1-q)}^{{N}_{k}}\displaystyle \left(\displaystyle \genfrac{}{}{0em}{}{{N}_{k}}{r}\right){\left(\displaystyle \frac{q}{1-q}\right)}^{r}\end{eqnarray} \tag{ 68b }$$

$$\begin{eqnarray}&&\equiv {(1-q)}^{{N}_{k}}\displaystyle \left(\displaystyle \genfrac{}{}{0em}{}{{N}_{k}}{r}\right)\tilde{q}\;{}^{r}\end{eqnarray} \tag{ 68c }$$

where we have defined $\tilde{q}\equiv q/(1-q)$. The average number of stars sampled in bin k will therefore be

$$\begin{eqnarray}&&\langle {r}_{k}\rangle =\displaystyle \sum _{r=1}^{{N}_{k}}{P}_{k,r}r\end{eqnarray} \tag{ 69a }$$

$$\begin{eqnarray}&&=\;{(1-q)}^{{N}_{k}}\displaystyle \sum _{r=1}^{{N}_{k}}\displaystyle \left(\displaystyle \genfrac{}{}{0em}{}{{N}_{k}}{r}\right)\tilde{q}\;{}^{r}r\end{eqnarray} \tag{ 69b }$$

$$\begin{eqnarray}&&=\;{(1-q)}^{{N}_{k}}\;\tilde{q}\displaystyle \frac{d}{d\tilde{q}}\;\displaystyle \sum _{r=1}^{{N}_{k}}\displaystyle \left(\displaystyle \genfrac{}{}{0em}{}{{N}_{k}}{r}\right)\tilde{q}\;{}^{r}\end{eqnarray} \tag{ 69c }$$

$$\begin{eqnarray}&&=\;{(1-q)}^{{N}_{k}}\;\tilde{q}\displaystyle \frac{d}{d\tilde{q}}\;{(1+\tilde{q})}^{{N}_{k}}\end{eqnarray} \tag{ 69d }$$

$$\begin{eqnarray}&&=\;{(1+\tilde{q})}^{-{N}_{k}}\;\tilde{q}\displaystyle \frac{d}{d\tilde{q}}\;{(1+\tilde{q})}^{{N}_{k}}\end{eqnarray} \tag{ 69e }$$

$$\begin{eqnarray}&&=\;{(1+\tilde{q})}^{-{N}_{k}}\;\tilde{q}{N}_{k}{(1+\tilde{q})}^{{N}_{k}-1}\end{eqnarray} \tag{ 69f }$$

$$\begin{eqnarray}&&\ =\;\displaystyle \frac{\tilde{q}}{1+\tilde{q}}\;{N}_{k}\end{eqnarray} \tag{ 69g }$$

$$\begin{eqnarray}&&=\;q\;{N}_{k}\end{eqnarray} \tag{ 69h }$$

which just says that if we sample a fraction q of the total number of stars we expect to sample that same fraction q of the stars in each bin. We can similarly calculate

$$\begin{eqnarray}&&\langle {r}_{k}^{2}\rangle ={(1+\tilde{q})}^{-{N}_{k}}\;{\left(\tilde{q}\displaystyle \frac{d}{d\tilde{q}}\right)}^{2}{(1+\tilde{q})}^{{N}_{k}-1}\end{eqnarray} \tag{ 70a }$$

$$\begin{eqnarray}&&=\;{(1+\tilde{q})}^{-{N}_{k}}\;\tilde{q}\;\displaystyle \frac{d}{d\tilde{q}}\;\tilde{q}{N}_{k}{(1+\tilde{q})}^{{N}_{k}-1}\end{eqnarray} \tag{ 70b }$$

$$\begin{eqnarray}&&=\;{{qN}}_{k}+{q}^{2}{N}_{k}({N}_{k}-1)\end{eqnarray} \tag{ 70c }$$

which gives

$$\begin{eqnarray}&&{\rm{\Delta }}{r}_{k}^{2}\equiv \langle {r}_{k}^{2}\rangle -\langle {r}_{k}{\rangle }^{2}=q(1-q){N}_{k}.\end{eqnarray} \tag{ 71a }$$

Since we expect ${{\rm{\Delta }}}^{(q)}\propto {\rm{\Delta }}r$, we thus find that

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\rho }^{(q)}(z)=C(z)\sqrt{q(1-q)}\end{eqnarray} \tag{ 72 }$$

for some function C(z). We therefore find, in the limit $q\to 0$,

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\rho }(z)=\underset{q\to 0}{\mathrm{lim}}\displaystyle \frac{{{\rm{\Delta }}}^{(q)}}{\sqrt{q}}\end{eqnarray} \tag{ 72a }$$

$$\begin{eqnarray}&&=\;\underset{q\to 0}{\mathrm{lim}}C(z)\sqrt{1-q}\end{eqnarray} \tag{ 72b }$$

$$\begin{eqnarray}&&=\;C(z).\end{eqnarray} \tag{ 72c }$$

We therefore have

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\rho }^{(q)}(z)={{\rm{\Delta }}}_{\rho }(z)\sqrt{q(1-q)}.\end{eqnarray} \tag{ 73a }$$

Which we can use to solve for Δ_ρ(z) by numerically calculating ${{\rm{\Delta }}}_{\rho }^{({q}_{0})}(z)$ at some reference value q₀:

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\rho }(z)=\displaystyle \frac{{{\rm{\Delta }}}_{\rho }^{({q}_{0})}(z)}{\sqrt{{q}_{0}(1-{q}_{0})}}.\end{eqnarray} \tag{ 74a }$$

For the reference value q₀ = 1/2 (sampling half the stars), we therefore have

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\rho }(z)=2\;{{\rm{\Delta }}}_{\rho }^{(1/2)}(z).\end{eqnarray} \tag{ 75a }$$

When simulating this sampling technique with a Hipparcos sample containing more stars than our data set, and extrapolating to q = 0, we found that Δ_ρ(z) was more correctly given by $1.97{{\rm{\Delta }}}_{\rho }^{(1/2)}(z)$. We are unsure of the reason for this discrepancy. At any rate, the effect of using the factor of 1.97 instead of 2 is small compared to the remaining errors in the analysis.

For f(w), the story is slightly different because of the normalization $\int {dw}\;f(w)=1$. Although Δr_k (where here k represents the kth bin in w-space) is still given by $\sqrt{q(1-q){N}_{k}}$, the quantity corresponding to f(w) here is ${N}_{k}/{\sum }_{k^{\prime} }{N}_{k^{\prime} }$ in the parent set or ${r}_{k}/{\sum }_{k^{\prime} }{r}_{k^{\prime} }$ in the sample set. Since for a large number of bins,

$$\begin{eqnarray}&&\displaystyle \sum _{k}{r}_{k}\simeq \displaystyle \sum _{k}\langle {r}_{k}\rangle =q\displaystyle \sum _{k^{\prime} }{N}_{k^{\prime} },\end{eqnarray} \tag{ 76 }$$

we have

$$\begin{eqnarray}&&{{\rm{\Delta }}}^{(q)}{f}_{k}=\displaystyle \frac{\sqrt{q(1-q){N}_{k}}}{q{\displaystyle \sum }_{k^{\prime} }{N}_{k^{\prime} }}.\end{eqnarray} \tag{ 77 }$$

Similarly to the case of Δ(z), we expect, for small q, that

$$\begin{eqnarray}&&\underset{q\to 0}{\mathrm{lim}}{{\rm{\Delta }}}^{(q)}{f}_{k}=\underset{q\to 0}{\mathrm{lim}}\displaystyle \frac{\sqrt{{{qN}}_{k}}}{q{\displaystyle \sum }_{k^{\prime} }{N}_{k^{\prime} }}.\end{eqnarray} \tag{ 78 }$$

We therefore find that

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{f}(w)=\sqrt{\displaystyle \frac{q}{1-q}}{{\rm{\Delta }}}_{f}^{(q)}(w).\end{eqnarray} \tag{ 79 }$$

Since by Equation (35), ${\rho }_{f}(z)\sim \int f$ scales linearly with f(w), we therefore have the same relationship for ρ_f (z):

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{{\rho }_{f}}(z)=\sqrt{\displaystyle \frac{q}{1-q}}{{\rm{\Delta }}}_{{\rho }_{f}}^{(q)}(z)\end{eqnarray} \tag{ 80 }$$

and, for q = 1/2,

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{{\rho }_{f}}(z)={{\rm{\Delta }}}_{{\rho }_{f}}^{(1/2)}(z).\end{eqnarray} \tag{ 81 }$$

To compute approximate probability densities, we can define a distance between expected and observed density

$$\begin{eqnarray}&&{\chi }_{{\rm{model}}}^{2}\equiv {\int }_{0}^{{z}_{{\rm{max}}}}\displaystyle \frac{{dz}}{{\rm{\Delta }}z}{\left|\displaystyle \frac{{\rho }_{{\rm{model}}}(z)-{\rho }_{{\rm{obs}}}(z)}{{{\rm{\Delta }}}_{\rho }(z)}\right|}^{2}\equiv \displaystyle \frac{{X}_{{\rm{model}}}}{{\rm{\Delta }}z},\end{eqnarray} \tag{ 82 }$$

where χ² is defined in the usual way as

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{n=1}^{k}\displaystyle \frac{{({x}_{n}-\langle {x}_{n}\rangle )}^{2}}{{{\rm{\Delta }}}_{n}^{2}}\end{eqnarray} \tag{ 83 }$$

with k representing the number of relevant degrees of freedom, and with x_n representing the degrees of freedom. The probability density in model space would then be proportional to

$$\begin{eqnarray}&&p({\chi }^{2})\sim \mathrm{exp}(-{\chi }^{2}/2).\end{eqnarray} \tag{ 84 }$$

However, we do not know the proportionality factor Δz since we do not know the number of relevant degrees of freedom for Φ(z). Although we know the number of points we are using, the correlation between nearby points reduces the number of relevant degrees of freedom. This factor is important because when in the exponential in Equation (84), it will affect the sharpness of the probability curve. A simple way to determine the degree of freedom length Δz as well as the number k of relevant degrees of freedom is to fit the probability distribution of P(X/Δz) to a χ² distribution:

$$\begin{eqnarray}&&P({\chi }^{2})=\displaystyle \frac{1}{{2}^{\tfrac{k}{2}}{\rm{\Gamma }}\;\left(\tfrac{k}{2}\right)}{({\chi }^{2})}^{\tfrac{k}{2}-1}{e}^{-{\chi }^{2}/2}.\end{eqnarray} \tag{ 85 }$$

The reason this method is only approximate is that there is degeneracy between the two parameters we are fitting, k and 1/Δz. To see this, note that the position of the peak of the χ² distribution grows linearly with k. Clearly, the position of the peak of P(X/Δz) also grows linearly with 1/Δz. Although other aspects of P(χ²) also depend on k, this degeneracy still persists and makes it difficult to fit the distributions. While it does not solve this degeneracy problem, it will prove slightly simpler to fit the cumulative distribution $C(X)={\int }_{0}^{X}{{dX}}^{\prime }P({X}^{\prime })$, given by the incomplete Gamma function:

$$\begin{eqnarray}&&C({\chi }^{2})=\displaystyle \frac{{\rm{\Gamma }}\left(\tfrac{k}{2},{\chi }^{2}\right)}{{\rm{\Gamma }}\left(\tfrac{k}{2}\right)},\end{eqnarray} \tag{ 86 }$$

Another factor that complicates this is that Φ(z) is functional depending on all the values of f(w). What this implies is that although the errors on f(w) and ρ(z) are approximately Gaussian, the errors on Φ(z) will not be. Figure 12 below shows the X distributions for A and F stars with the best fit Δz and k = 12 and 8 (respectively) degrees of freedom. However, because of the degeneracy between k and 1/Δz, the best fit values are very approximate and the derived probability density, p(Σ_D), should therefore only be regarded as qualitative.

Zoom In Zoom Out Reset image size

As explained in Section 5.3, we expect that the HF relation will hold for long-time averages. We will demonstrate that

$$\begin{eqnarray*}&&{\overline{\rho (z)}}_{[{\rm{\Phi }}]}{\bar{\rho {(0)}^{-1}}}_{[{\rm{\Phi }}]}=\int {dw}\;{\bar{{f}_{z=0}}}_{[{\rm{\Phi }}]}(\sqrt{{w}^{2}+2{\rm{\Phi }}(z)})\end{eqnarray*}$$

where the ${\bar{\cdot }}_{[{\rm{\Phi }}]}$ represents the time average under evolution in the potential Φ(z). However, if this potential Φ is the same potential as under the square root $\sqrt{{w}^{2}+{\rm{\Phi }}(z)}$, then (87) will be satisfied trivially for any initial conditions. This can be demonstrated by taking the time-dependent density and in-plane velocity distributions to be:

$$\begin{eqnarray}&&\rho (z,t)=\displaystyle \frac{1}{{A}_{0}}{\displaystyle \sum }_{i}\delta (z-{z}_{i}(t))\end{eqnarray} \tag{ 87a }$$

$$\begin{eqnarray}&&{f}_{z=0}(w,t)=\displaystyle \frac{{\displaystyle \sum }_{i}\delta (w-{w}_{i}(t))\theta (| {z}_{i}(t)| \lt \epsilon /2)}{{\displaystyle \sum }_{j}\theta (| {z}_{j}(t)| \lt \epsilon /2)}\end{eqnarray} \tag{ 87b }$$

where the sums are over all stars in the tracer population, A₀ is the cross-sectional area of our sample, and $\theta (| {z}_{i}(t)| \lt \epsilon /2)$ assures that z_i(t) is within some appropriate small distance /2 of the plane. We can now compute the long-time average:

$$\begin{eqnarray}&&\overline{\rho (z)}=\underset{t\to \infty }{\mathrm{lim}}\displaystyle \frac{1}{t}{\int }_{0}^{t}{dt}^{\prime} \;\rho (z,t^{\prime} )\end{eqnarray} \tag{ 88a }$$

$$\begin{eqnarray}&&=\;\displaystyle \frac{1}{{A}_{0}}\displaystyle \sum _{i}\underset{t\to \infty }{\mathrm{lim}}\displaystyle \frac{1}{t}{\int }_{0}^{t}{dt}^{\prime} \;\delta (z-{z}_{i}(t^{\prime} ))\end{eqnarray} \tag{ 88b }$$

$$\begin{eqnarray}&&=\;\displaystyle \frac{1}{{A}_{0}}\displaystyle \sum _{i}\underset{{N}_{i}\to \infty }{\mathrm{lim}}\displaystyle \frac{1}{{N}_{i}{T}_{i}/2}{N}_{i}{\int }_{0}^{{T}_{i}/2}{dt}\;\delta (z-{z}_{i}(t))\end{eqnarray} \tag{ 88c }$$

$$\begin{eqnarray}&&=\;\displaystyle \frac{1}{{A}_{0}}\displaystyle \sum _{i}\displaystyle \frac{2}{{T}_{i}}{\int }_{-{z}_{{\rm{max}}}}^{{z}_{{\rm{max}}}}\displaystyle \frac{{{dz}}_{i}}{{w}_{i}({z}_{i})}\;\delta (z-{z}_{i})\end{eqnarray} \tag{ 88d }$$

$$\begin{eqnarray}&&=\;\displaystyle \frac{1}{{A}_{0}}\displaystyle \sum _{i}\displaystyle \frac{2}{{T}_{i}| {w}_{i}(z)| }\end{eqnarray} \tag{ 88e }$$

where we have used the fact that the trajectories of the stars are periodic with individual periods, T_i, and where $| {w}_{i}({z}_{i})|$ is star i's speed (fixed by energy conservation) at height z_i. We can proceed similarly for the velocity distribution. For a large number of stars, we can write the denominator as

$$\begin{eqnarray}&&\displaystyle \sum _{i}\theta (| {z}_{i}(t)| \lt \epsilon /2)=\epsilon \;\rho (0,t)\;{A}_{0}.\end{eqnarray} \tag{ 89 }$$

We now proceed to write

$$\begin{eqnarray}&&\overline{{f}_{z=0}(w)}=\underset{t\to \infty }{\mathrm{lim}}\displaystyle \frac{1}{t}{\int }_{0}^{t}{dt}^{\prime} \;{f}_{z=0}(w,t^{\prime} )\end{eqnarray} \tag{ 90a }$$

$$\begin{eqnarray}&&=\;\displaystyle \sum _{i}\underset{t\to \infty }{\mathrm{lim}}\displaystyle \frac{1}{t}{\int }_{0}^{t}{dt}^{\prime} \;\displaystyle \frac{\delta (w-{w}_{i}(t^{\prime} ))\theta (| {z}_{i}(t^{\prime} )| \lt \epsilon /2)}{\epsilon \;\rho (0,t){A}_{0}}\end{eqnarray} \tag{ 90b }$$

Now, when star i is at z = 0, its velocity will be maximum and will be equal to ${w}_{i}=\sqrt{2{E}_{i}}$. The amount of time that star i spends between $z=\mp \epsilon /2$ will therefore be given by $\epsilon /\sqrt{2{E}_{i}}$. In each period, the integral will therefore receive a contribution $\epsilon /\sqrt{2{E}_{i}}\;\delta (w-\sqrt{2{E}_{i}})$ on the way up and $\epsilon /\sqrt{2{E}_{i}}\;\delta (w+\sqrt{2{E}_{i}})$ on the way down. We can achieve this by replacing

$$\begin{eqnarray}&&\theta (| {z}_{i}(t^{\prime} )| \lt \epsilon )=\displaystyle \frac{\epsilon }{\sqrt{2{E}_{i}}}\displaystyle \sum _{{n}_{i}=0}^{{N}_{i}}\delta (t-{t}_{0i}-{n}_{i}{T}_{i}/2)\end{eqnarray} \tag{ 91a }$$

$$\begin{eqnarray}&&=\;\displaystyle \frac{\epsilon }{\sqrt{2{E}_{i}}}\displaystyle \sum _{{n}_{i}=0}^{{N}_{i}/2}(\delta (t-{t}_{0i}-{n}_{i}{T}_{i})\\ &&\quad +\delta (t-{t}_{0i}-({n}_{i}+1/2){T}_{i}))\end{eqnarray} \tag{ 91b }$$

where the first crossing of the plane for star i happens at time t_0i. We thus have

$$\begin{eqnarray}&&\overline{{f}_{z=0}(w)}=\displaystyle \frac{1}{\epsilon {A}_{0}}\displaystyle \sum _{i}\displaystyle \frac{\epsilon }{\sqrt{2{E}_{i}}}\underset{{N}_{i}\to \infty }{\mathrm{lim}}\displaystyle \frac{1}{{N}_{i}{T}_{i}/2}{\int }_{0}^{{N}_{i}{T}_{i}/2}{dt}\;\displaystyle \frac{1}{\rho (0,t)}...\end{eqnarray} \tag{ 92a }$$

$$\begin{eqnarray}&&\times \;\left[\delta (w-{\rm{sign}}({w}_{i}({t}_{0i}))\sqrt{2{E}_{i}})\displaystyle \sum _{{n}_{i}=0}^{{N}_{i}/2}\delta (t-{t}_{0i}-{n}_{i}{T}_{i})\right....\end{eqnarray} \tag{ 92b }$$

$$\begin{eqnarray}&&+\delta (w+{\rm{sign}}({w}_{i}({t}_{0i}))\sqrt{2{E}_{i}})\\ &&\quad \left.\times \quad \displaystyle \sum _{{n}_{i}=0}^{{N}_{i}/2}\delta (t-{t}_{0i}-({n}_{i}+1/2){T}_{i})\right]\end{eqnarray} \tag{ 92c }$$

$$\begin{eqnarray}&&=\;\displaystyle \frac{1}{{A}_{0}}\displaystyle \sum _{i}\displaystyle \frac{1}{\sqrt{2{E}_{i}}}\underset{{N}_{i}\to \infty }{\mathrm{lim}}\displaystyle \frac{1}{{N}_{i}{T}_{i}/2}...\end{eqnarray} \tag{ 92d }$$

$$\begin{eqnarray}&&\times \;\left[\delta (w-{\rm{sign}}({w}_{i}({t}_{0i}))\sqrt{2{E}_{i}})\displaystyle \sum _{{n}_{i}=0}^{{N}_{i}/2}\displaystyle \frac{1}{\rho (0,{t}_{0i}+{n}_{i}{T}_{i})}\right....\end{eqnarray} \tag{ 92e }$$

$$\begin{eqnarray}&&+\;\left.\delta (w+{\rm{sign}}({w}_{i}({t}_{0i}))\sqrt{2{E}_{i}})\displaystyle \sum _{{n}_{i}=0}^{{N}_{i}/2}\displaystyle \frac{1}{\rho (0,{t}_{0i}+({n}_{i}+1/2){T}_{i})}\right].\end{eqnarray} \tag{ 92f }$$

Now, unless the periods of ρ(0, t) and of the trajectories of the individual stars divide each other, ${\sum }_{{n}_{i}=0}^{{N}_{i}/2}1/\rho (0,{t}_{0i}+{n}_{i}{T}_{i})$ will just be, in the large N_i limit, N_i/2 times a time average of 1/ρ(0, t). Since the stars with periods dividing that of ρ(0, t) is a set of measure zero, we have

$$\begin{eqnarray}&&=\;\displaystyle \frac{1}{{A}_{0}}\displaystyle \sum _{i}\displaystyle \frac{1}{\sqrt{2{E}_{i}}}\underset{{N}_{i}\to \infty }{\mathrm{lim}}\displaystyle \frac{1}{{N}_{i}{T}_{i}/2}...\end{eqnarray} \tag{ 93a }$$

$$\begin{eqnarray}&&\times \;\left[\delta (w-{\rm{sign}}({w}_{i}({t}_{0i}))\sqrt{2{E}_{i}})\displaystyle \frac{{N}_{i}}{2}\;\bar{\rho {(0)}^{-1}}\right....\end{eqnarray} \tag{ 93b }$$

$$\begin{eqnarray}&&+\;\left.\delta (w+{\rm{sign}}({w}_{i}({t}_{0i}))\sqrt{2{E}_{i}})\displaystyle \frac{{N}_{i}}{2}\;\bar{\rho {(0)}^{-1}}\right]\end{eqnarray} \tag{ 93c }$$

$$\begin{eqnarray}&&=\;\displaystyle \frac{1}{{A}_{0}}\displaystyle \sum _{i}\displaystyle \frac{1}{\sqrt{2{E}_{i}}}\displaystyle \frac{1}{{T}_{i}}\;\bar{\rho {(0)}^{-1}}[\delta (w-\sqrt{2{E}_{i}})+\delta (w+\sqrt{2{E}_{i}})]\end{eqnarray} \tag{ 93d }$$

$$\begin{eqnarray}&&=\;\displaystyle \frac{1}{{A}_{0}}\displaystyle \sum _{i}\displaystyle \frac{1}{\sqrt{2{E}_{i}}}\displaystyle \frac{1}{{T}_{i}}\;\bar{\rho {(0)}^{-1}}\;2\sqrt{2{E}_{i}}\;\delta ({w}^{2}-2{E}_{i})\end{eqnarray} \tag{ 93e }$$

$$\begin{eqnarray}&&=\;\displaystyle \frac{\bar{\rho {(0)}^{-1}}}{{A}_{0}}\displaystyle \sum _{i}\displaystyle \frac{2}{{T}_{i}}\;\delta ({w}^{2}-2{E}_{i}).\end{eqnarray} \tag{ 93f }$$

We can now evaluate ${\bar{{f}_{z=0}(w^{\prime} )}| }_{w^{\prime} =\sqrt{{w}^{2}+2{\rm{\Phi }}(z)}}$, using the fact that ${w}_{i}{(z)}^{2}/2+{\rm{\Phi }}(z)={E}_{i}$:

$$\begin{eqnarray}&&{\bar{{f}_{z=0}(w^{\prime} )}| }_{w^{\prime} =\sqrt{{w}^{2}+2{\rm{\Phi }}(z)}}\\ &&\quad =\displaystyle \frac{\bar{\rho {(0)}^{-1}}}{{A}_{0}}\displaystyle \sum _{i}\displaystyle \frac{2}{{T}_{i}}\;\delta ({w}^{2}-2({E}_{i}-{\rm{\Phi }}(z)))\end{eqnarray} \tag{ 94a }$$

$$\begin{eqnarray}&&=\;\displaystyle \frac{\bar{\rho {(0)}^{-1}}}{{A}_{0}}\displaystyle \sum _{i}\displaystyle \frac{2}{{T}_{i}}\;\delta ({w}^{2}-{w}_{i}{(z)}^{2})\end{eqnarray} \tag{ 94b }$$

where $| {w}_{i}(z)|$ is the speed of star i at height z. This gives

$$\begin{eqnarray}&&\displaystyle \int {dw}\;{\overline{{f}_{z=0}({w}^{\prime })}| }_{{w}^{\prime }=\sqrt{{w}^{2}+2{\rm{\Phi }}(z)}}\\ &&\quad =\displaystyle \int {dw}\;\displaystyle \frac{\bar{\rho {(0)}^{-1}}}{{A}_{0}}\displaystyle \sum _{i}\displaystyle \frac{2}{{T}_{i}}\;\delta ({w}^{2}-{w}_{i}{(z)}^{2})\end{eqnarray} \tag{ 95a }$$

$$\begin{eqnarray}&&=\;\displaystyle \frac{\bar{\rho {(0)}^{-1}}}{{A}_{0}}\displaystyle \sum _{i}\displaystyle \frac{2}{{T}_{i}| {w}_{i}(z)| }\end{eqnarray} \tag{ 95b }$$

which, by Equation (88e), is equal to $\overline{\rho (z)}$.