THE FUTILE SEARCH FOR GALACTIC DISK DARK MATTER

The mass distribution of our Galaxy can be derived from a knowledge of forces both parallel and perpendicular to its midplane. Our knowledge of the parallel forces comes by inference from measurements of the generally flat rotation profiles of other spiral galaxies and by direct measurements of Milky Way stellar motions that are used to estimate galactic rotation. The conventional explanation for the observations is that dark matter (DM) in a halo somehow "conspires" with visible matter to keep the rotation profiles flat.

Studies of forces perpendicular to the galactic midplane have also led to DM hypotheses, but this DM is confined to the disk where it purportedly has a higher density than the halo DM. The search for disk DM, started 95 years ago by Öpik (1915), has a controversial history. It has lasted even longer than the 75 year futile search for Vulcan (Baum & Sheehan 1997).

Kapteyn (1922) was the first to apply the kinetic theory of gases to interpret star counts and estimate the gravitational field of the disk in the vicinity of the Sun. Oort (1932) refined the theory and wrote the seminal paper on the topic. The gravitational potential can be estimated from observations of the midplane or column densities and the velocity distributions of all observable stars and gases, and it can independently be estimated from observations of limited sets of tracer stars. A discrepancy between the two estimates may be indicative of DM. Numerous papers were published following Oort's lead, culminating with the disk model of Bahcall et al. (1992, hereafter BFG). The data they used to determine the gravitational potential comprised 13 stellar components (excluding a very small spheroid component) plus four gas components. They found that their data supported the hypothesis that disk DM exists. Yet the significance of this finding is debatable because it has not been confirmed on using K giant data from the Hipparcos catalog (Korchagin et al. 2003; Holmberg & Flynn 2004). Specifically, Holmberg & Flynn (2004, hereafter HF) updated the BFG model and proffered the most comprehensive attempt yet at modeling the galactic disk mass density in the solar vicinity. The HF model has 11 stellar components and the same 4 gas components as the BFG model, plus a thick disk constituent, but no disk DM.

Another approach is the study of the thickness of the neutral hydrogen component of the Milky Way. Kalberla (2003) and Kalberla et al. (2007) analyzed flaring data in a Galaxy model with six gaseous and four stellar components (including a central bulge and bar). To explain the observations, they required DM in a thick disk and in a ring—plus the usual DM halo. Using a computer simulation of the evolution of galaxies, Read et al. (2008) concluded that our Galaxy includes disk DM. Conversely, analyses of gravitational microlens observations toward the galactic center (Popowski et al. 2005; Hamadache et al. 2006) found no evidence for DM in the disk; and, alternatively, the semblance of a similar, yet fictitious, DM disk is predicted by MOND (Milgrom 2001). Disk DM, like halo DM, is a broadly accepted working hypothesis, but it remains just a hypothesis that demands further investigation. Independently being able to confirm or deny the existence of a DM disk would be important to fundamental physics and cosmology, so we have devised a technique to analyze the z and v_z distributions of Milky Way stars and gases by applying the virial theorem (Section 2.3) to the HF disk model (Section 3.1).

2.1. The Virial Theorem

We adopt the BFG toy model for which the mass density of a uniformly flat Milky Way disk is ρ(z) = ρ(−z). In the vicinity of the Sun, at distance r = r_☉ from the Galaxy rotation axis, the disk mass surface density (also called the column density) is

$$\begin{equation*} \Sigma =\int _{-\infty }^{\infty }\rho (z)\,{\it dz}; \end{equation*}$$

its kinetic energy surface density is

$$\begin{equation} T=\Sigma \, \sigma ^2/2 \end{equation} \tag{ 1 }$$

(which, in effect, defines the variance σ² of the vertical velocity v_z); and its gravitational energy surface density is

$$\begin{equation} V=-\int _{-\infty }^{\infty }\rho \Psi ^{\prime }\,z\,{\it dz} = \int _{-\infty }^{\infty }\rho gz\,{\it dz}, \end{equation} \tag{ 2 }$$

where Ψ is the gravitational potential and g is the gravitational specific force (acceleration). For (z/r_☉)² ≪ 1, the gravitational energy between horizontal layers parallel to the midplane is effectively independent of the gravitational energy of the orbit about the Galaxy rotation axis. Then if the disk does not secularly disperse, the virial theorem (Goldstein 1959; Collins 1978) applies:

$$\begin{equation} \overline{2T+V}=0, \end{equation} \tag{ 3 }$$

where the overbar is the time averaging operator. If the vertical mass column is in a state of equilibrium [$\dot{\rho }(z)=0$], then the instantaneous value of 2T + V is 0. If instead the mass column is not in a state of equilibrium, then there is an interchange of energy density back and forth between T and V, and consequently 2T + V oscillates about 0.

2.2. The Poisson–Boltzmann Equation

For a disk vertical mass column that is in an equilibrium state, the Poisson equation applies:

$$\begin{equation} \Psi ^{\prime \prime }(z)=4\pi G\rho (z), \end{equation} \tag{ 4 }$$

where G is the gravitational constant. The boundary conditions are Ψ(0) = Ψ'(0) = 0. If an assembly of disk stars is isothermal, then σ² is independent of z, and if it is not isothermal, then it is composed of multiple sub-assemblies that are isothermal. The density distribution of all disk stars and gas clouds, ultimately partitioned among N isothermal components, is the sum of N Boltzmann distributions,

$$\begin{equation} \rho (z)=\sum _{n=1}^{N}\rho _{n}(z)=\sum _{n=1}^{N}\rho _n(0) \, \exp \big [-\Psi (z)/\sigma _{n}^{2}\big]. \end{equation} \tag{ 5 }$$

Combining the Poisson (4) and Boltzmann (5) equations gives the Poisson–Boltzmann equation,

$$\begin{equation} \Psi ^{\prime \prime }=4\pi G \sum _{n=1}^{N} \rho _{n}\exp \big [-\Psi /\sigma _{n}^{2}\big]. \end{equation} \tag{ 6 }$$

Spitzer derived an analytic solution for the Poisson–Boltzmann equation with a single (n = 1) isothermal component,

$$\begin{equation} \Psi ^{\prime \prime }=4\pi G\rho (0) \exp \,[-\Psi /\sigma ^{2}]=2k^2\sigma ^2\exp \,[-\Psi /\sigma ^{2}], \end{equation} \tag{ 7 }$$

where

$$\begin{equation} k^2\sigma ^2=2\pi G\rho (0). \end{equation} \tag{ 8 }$$

Multiply both sides of Equation (7) by Ψ', and then integrate,

$$\begin{equation} (\Psi ^{\prime })^{2}=4k^2\sigma ^4 [1-\exp \, (-\Psi /\sigma ^{2})]. \end{equation} \tag{ 9 }$$

Define

$$\begin{equation*} u^{-2}=\rho (z)/\rho (0)=\exp (-\Psi /\sigma ^{2}\,); \end{equation*}$$

then

$$\begin{equation*} \Psi =2\sigma ^2\,\ln u, \end{equation*}$$

so the left-hand side of Equation (9) is

$$\begin{equation} (\Psi ^{\prime })^{2}=4\sigma ^4 (u^{\prime }/u)^2, \end{equation} \tag{ 10 }$$

and the right-hand side of Equation (9) is

$$\begin{equation} 4k^2\sigma ^4 [1-\exp \, (-\Psi /\sigma ^{2})]=4k^2\sigma ^4[1-u^{-2}]. \end{equation} \tag{ 11 }$$

Equating Equations (10) and (11),

$$\frac{du}{\sqrt{u^{2}-1}}= k {\it dz}.$$

Integrating this equation from u(0) = 1 to u(z) gives u(z) = cosh kz, and then the Spitzer solution,

$$\begin{equation} \rho (z)=\rho (0) u^{-2}=\rho (0)\;\mbox{sech}^{2} kz \,, \end{equation} \tag{ 12 }$$

and

$$\begin{equation} \Psi (z)=2\sigma ^{2}\ln u=2\sigma ^{2}\ln \cosh kz \,. \end{equation} \tag{ 13 }$$

This solution conforms with the steady state virial theorem (2T + V = 0) because, by Equations (12) and (1),

$$\begin{equation} \Sigma =\rho (0)\int _{-\infty }^{\infty }\,\mbox{sech}^{2} kz\,{\it dz}=2\rho (0)/k, \end{equation} \tag{ 14 }$$

$$\begin{equation} 2T=2\rho (0)\sigma ^2/k; \end{equation} \tag{ 15 }$$

and, by Equations (12), (13), (2), and (15),

$$\begin{eqnarray} V&=&-2\rho (0)\sigma ^2\int _{-\infty }^{\infty }\mbox{sech}^{2} kz\;\frac{d(\ln \cosh kz)}{{{\it dz}}}\,z\,{\it dz} \nonumber \\ & = &-2\rho (0)\sigma ^2\int _{-\infty }^{\infty }\mbox{sech}^{2} kz\,\tanh kz\,kz\,{\it dz}\nonumber\\ & =& -2\rho (0)\sigma ^2/k = -2T. \end{eqnarray} \tag{ 16 }$$

We are aware of no analytic solution to Equation (6) for N>1, but it can be integrated numerically to find a unique, steady state solution. (This is not a trivial point because no such solutions can exist for three-dimensional Poisson–Boltzmann integrations such as modeling the structure of an isothermal star (Pestaña & Eckhardt 2007).) We are accordingly poised to examine the problem by using the virial theorem. The virial theorem approach has an advantage over Poisson–Boltzmann integration because it merely requires that $\overline{2T+V}=0$, so an oscillatory solution is feasible, whereas the Poisson–Boltzmann integration requires that 2T + V = 0 and so only a time invariant solution is allowed.

2.3. The Vertical Structure of the Disk

An infinitesimal layer of matter, ρ(z) dz, is the source of a specific force

$$\begin{equation*} g=2\pi G\rho (z)\,{\it dz} \end{equation*}$$

in the direction of the layer. Thus the gravitational energy surface density between infinitesimal layers of matter at z and z+a is

$$\begin{equation*} {-}{\it ga}\,\rho (z+a)=-2\pi Ga\,\rho (z)\,\rho (z+a)\,{\it dz}, \end{equation*}$$

and the gravitational energy surface density due to all layers that are separated by distance a is

$$\begin{eqnarray*} {-}2\pi Ga\!\!\int _{-\infty }^{\infty }\!\!\rho (z)\rho (z+a)\,{\it dz} &=& -2\pi G\Sigma ^2a\!\!\int _{-\infty }^{\infty }\!\!p(z)p(z+a)\,{\it dz}\nonumber\\ &=& -2\pi G\Sigma ^2 q(a)a, \end{eqnarray*}$$

where p(z) is the frequency function

$$\begin{equation} p(z)=\rho (z)/\Sigma, \end{equation} \tag{ 17 }$$

and q(a) is the frequency function given by the convolution

$$\begin{equation} q(a)= \int _{-\infty }^{\infty }p(z)\,p(z+a)\,{\it dz}. \end{equation} \tag{ 18 }$$

The total gravitational energy surface density is then

$$\begin{equation} V=-2\pi G\Sigma ^2\int _{0}^{\infty }q(a)a\,da. \end{equation} \tag{ 19 }$$

The lower limit of this integral is 0 because a ⩾ 0 to avoid double counting.

Combining Equations (1), (3) (not time averaged), (17), and (19) gives

$$\begin{equation} \frac{\rho (0)\sigma ^2}{\pi G\Sigma ^2}=2p(0)\int _0^{\infty }q(a)a\,da. \end{equation} \tag{ 20 }$$

Since both sides of Equation (20) are dimensionless, they need not have the same units of length. Then for computational convenience, we implicitly define the unit of length on the right-hand side (and equations leading to it) by setting

$$\begin{equation} p(0)=1/2, \end{equation} \tag{ 21 }$$

so Equation (20) becomes the virial theorem criterion for a state of equilibrium,

$$\begin{equation} R=Q, \end{equation} \tag{ 22 }$$

where

$$\begin{equation} R=\frac{\rho (0)\sigma ^2}{\pi G\Sigma ^2}\,, \end{equation} \tag{ 23 }$$

and

$$\begin{equation} Q=\int _0^{\infty }q(a)a\,da. \end{equation} \tag{ 24 }$$

In Equation (23), Σ/ρ(0) and σ²/GΣ retain the dimension pc. Using Equations (8) and (14), this ratio for the Spitzer solution is

$$\begin{equation} R_{S}=1/2. \end{equation} \tag{ 25 }$$

We have examined various trial frequency functions, p(z), using them in Equations (18) and (24) to evaluate Q. These functions are not arbitrary. Each must be bell shaped: an even function of z (that is, p(z) = p(−z)), with a single peak at p(0) = 1/2 (that is, zp'(z) ⩽ 0) and, of course, ∫^∞_−∞p(z) dz = 1. Moreover, we assume that p(z) and all its derivatives are continuous, and that p(z) is not a pathological frequency function, meaning that its nth moment ∫^∞_−∞zⁿp(z) dz exists. Therefore, for large z, p(z) ∼ exp(−αz^2m), where α is a positive constant and m ⩾ 1 is an integer. With these provisos, it turns out that for all plausible frequency functions,

$$\begin{equation} Q\approx 1/2. \end{equation} \tag{ 26 }$$

We offer two simple examples (for which m = 1).

• 1.  
  For the Spitzer function (see Equations (12), (17), and (21), with k = 1),

  $$\begin{equation} p_S(z)=\frac{1}{2}\,\mathrm{sech}^2 z, \end{equation} \tag{ 27 }$$

  which applies exactly for a single component model, but still might be a suitable approximation for a multiple component model. For the Spitzer function,

  $$\begin{equation} Q_S= 1/2, \end{equation} \tag{ 28 }$$

  so, as expected (see Equation (25)), R_S = Q_S.
• 2.  
  For the Gaussian frequency function,

  $$\begin{equation} p_G(z)=\frac{1}{2}\exp \left[-\pi \left(\frac{z}{2}\right)^2 \right], \end{equation} \tag{ 29 }$$

  for which

  $$\begin{equation} Q_G=\sqrt{2}/\pi =0.45. \end{equation} \tag{ 30 }$$

These are merely examples. The exact value for Q is not important for this study. Equation (26) is entirely adequate to restate the condition for equilibrium, Equation (22), as

$$\begin{equation} R=Q\approx 1/2. \end{equation} \tag{ 31 }$$

3.1. The HF Model

The 15 components of the HF disk mass model were itemized by Flynn et al. (2006). From these, we calculate

$$\begin{equation} \Sigma =\sum _1^{15}\Sigma _i=49.3\,{\mathrm{pc}^{-2}}\, M_{\odot }, \end{equation} \tag{ 32 }$$

$$\begin{eqnarray} \pi G\Sigma ^2&=&\pi \times 49.3^2\,{\mathrm{pc}^{-4}}\; (GM_{\odot })\,M_{\odot } \nonumber\\ &=&34.4\times 10^{-27}\,{\mathrm{pc}^{-1}{\mbox{s}}^{-2}}\, M_{\odot }, \end{eqnarray} \tag{ 33 }$$

$$\begin{equation} \sigma ^2=\sum _1^{14}(\Sigma _i/\Sigma)\sigma _i^2=587\, \mathrm{km^2\,s^{-2} =614\times 10^{-27}\,pc^2\,s^{-2}}, \end{equation} \tag{ 34 }$$

$$\begin{equation} \rho (0)=\sum _1^{15}\rho _i(0)=0.0914\, {\mathrm{pc}^{-3}}\,M_{\odot }, \end{equation} \tag{ 35 }$$

$$\begin{equation} \rho (0)\sigma ^2=56.1\times 10^{-27}\,{\mathrm{pc}^{-1}\mbox{s}^{-2}}\,M_{\odot }, \end{equation} \tag{ 36 }$$

and (see Equation (23)),

$$\begin{equation} R_{{\rm HF}}= 1.63, \end{equation} \tag{ 37 }$$

which is too large compared with Q ≈ 1/2 for an equilibrium solution to exist. Hence, a Poisson–Boltzmann integration using the HF model, or any other such model, is an inappropriate and futile exercise. Nevertheless, HF performed this integration without disk DM and found the result to be compatible with vertical density profiles of K giant tracers. We look at the same data and conclude that the HF model can be compatible with a Poisson–Boltzmann steady state solution only if disk DM exists. For example, suppose that DM has the same vertical profile and v_z variance as visible matter, but that visible matter accounts for only 30% of the total surface mass density. Then ρ(0) and Σ increase by the factor 10/3, and thus R_HF is modified to 1.63 × 3/10 = 0.49, which is a plausible value for Q. However, the revised model no longer satisfies the Poisson–Boltzmann equation and the change increases the vertical accelerations of the tracer stars by the factor 10/3 and decreases their scale height by 70%, and such a model conflicts drastically with K giant vertical profiles. Introducing DM solves one problem but creates another. The overall problem can only be resolved by abandoning the requirement that the disk z profile is in a steady state or that gravity is Newtonian.

3.2. The BFG Models

A Poisson–Boltzmann solution is feasible only if 2T + V = 0, but this requirement is incompatible with observations; compelling a steady state solution is tantamount to compelling the existence of DM or non-Newtonian gravity. Suppose instead that a disk column's vertical profile varies periodically with time. Then 2T + V oscillates about $\overline{2T+V}=0$, so that at any epoch, 2T + V could be substantially greater or less than zero. In other words, R oscillates about Q. We use this broadened aspect to interpret the HF data and uncover interesting possible ramifications.

The analysis so far has been quantitative, but from here on it will be descriptive and speculative—and accordingly less defensible. Using dimensional analysis, we reckon that the z oscillation period is of the order [Gρ(0)]^−1/2 = 50 million years. The estimated period is roughly a quarter of a galactic year, but this estimate is inexact because it is based on dimensional considerations only. It would have been larger if we had used the mean midplane density $\overline{\rho (0)}$ over one cycle instead of the current midplane density (which we believe to be near the maximum value), and there undoubtedly are dimensionless factors that would arise in a precise analysis, decreasing the estimated period by as much as a factor of 10, but increasing it no more than fourfold (each disk star orbiting the galactic center has a retrograde precession caused by the other disk stars, so the z oscillation period must be less than the orbital period). What we find puzzling about the stellar orbits is the cause of the symmetric clustering of their ascending and descending nodes.

Disk column oscillations are exchanges between gravitational, V ∝ −Σ³/ρ(0), and kinetic, T ∝ Σσ², energy surface densities. Since Σ does not vary with time, oscillations in ρ(0) entail oscillations in σ²; and the ρ(0) and σ² variations are in phase, so that the peaks in midplane mass and kinetic energy densities coincide. Equation (37) tells us that at the current oscillation phase ρ(0)σ² is about three times higher than its mean value, but it does not tell us how close ρ(0)σ² is to its peak value. At the peak, there would be no net flow of stars toward or away from the midplane. Unaware of any evidence for such a flow, we suppose that ρ(0)σ² is close to its peak and corroborate this with data from Table A.1 of Famaey et al. (2005), selecting single giants (* tag in Column 24) in the smooth background (g = 3 in Column 34) that are in the two zones defined by 100 ⩽ |z| ⩽ 150 pc. There are 467 stars in the northern zone and 224 stars in the southern zone. We calculate that the two zones are approaching each other at 1.3 ± 1.4 km s⁻¹. The virial theorem requires an oscillation with substantial energy flow, but these numbers tell us that there is negligible net mass flow now (at least for the giants), so the present epoch appears to be one of the maximum ρ(0)σ². With little net mass flow toward or away from the midplane, the v_z velocity distributions may appear to be isothermal, but at epochs of maximum flow, the net flow velocities would be of the order Σ[G/ρ(0)]^1/2 = 10 km s^-1. Considering its imprecision, this is not much different from σ_HF = 24 km s^-1. A net outward flow would be obvious in about 10 million years.

We extend our z coordinate system to a standard cylindrical coordinate (r, ϑ, z) system. Then, we expect oscillations elsewhere on the disk with periods proportional to $\overline{\rho (r,\vartheta,0)}^{\;-1/2}$. Because these mean midplane densities vary, so do the z oscillation periods. Just small differences in $\overline{\rho (r,\vartheta,0)}$ between neighboring disk columns can effect large changes in their relative z oscillation phases and, therefore, large differences in their instantaneous midplane densities ρ(r, ϑ, 0).

The increased midplane mass and kinetic energy densities caused by oscillations are conducive to the formation of stars, including the massive O and B stars whose luminosity lifetimes would be relatively ephemeral compared with the disk oscillation periods. We therefore conjecture that the disk spirals are where the midplane mass and kinetic energy densities are now high, so stars are being formed there at relatively high rates. Vertical oscillations also occur in the gaps between the spirals, but their midplane mass and kinetic energy densities are now low (and their column gravitational energies are high) because of phase differences, so brighter stars that were created when their midplane mass and kinetic energy densities had been high have since waned. A computer simulation might tell us the shape that p(z) takes in the gaps, for we do not even know whether the frequency function remains unimodal. But neither do we know whether our conjecture concerning the formation of spirals is correct; nor, indeed, can we be sure that the disk actually oscillates, for conceivably the HF data can be explained with a non-Newtonian gravitational theory. However, we maintain that the conventional Poisson–Boltzmann interpretation of the data is demonstrably wrong, and we offer the disk oscillation hypothesis as one plausible alternative.

We have used the virial theorem to study the problem of determining the density of matter near the Sun because it is less restrictive than the Poisson–Boltzmann approach. In fact, we have shown that all the comparisons of tracer star densities with Poisson–Boltzmann models or with equilibrium models in general, for example, Kapteyn (1922), Oort (1932), Bahcall et al. (1992), Kalberla (2003), and Holmberg & Flynn (2000, 2004), have been wrong. The reason is that the disk is not in a state of equilibrium (2T + V ≠ 0). Furthermore, any model made compatible with 2T + V = 0 by adding DM will then be incompatible with Poisson–Boltzmann integration models and with observed tracer star densities. No evidence can be found for DM by using this approach. Explaining the evidence (e.g., the BFG, Kalberla (2003), and HF disk models, plus tracer star densities) requires disk oscillations or, perhaps, non-Newtonian gravity. A cursory inquiry into the possibility of disk oscillations without DM leads to intriguing hints concerning the structure and dynamics of a galactic disk.

We contend that any astrophysical argument predicated on the existence of a DM disk—for example, Stothers (1998), Sánchez-Salcedo (1999), de Boer et al. (2005), Kalberla et al. (2007), Revaz et al. (2009), and Gates & Gyuk (2004)—does not warrant consideration. We reject the recently proposed model for the formation of disk galaxies within the ΛCDM framework because it would produce a DM disk (Read et al. 2008; Purcell et al. 2009). Finally, there is no fictitious semblance of a DM disk as predicted by MOND (Milgrom 2001; Bienaymé et al. 2009), so our analysis argues against MOND.

We thank the referee for positive comments and constructive suggestions, and especially for drawing out attention to a key numerical error and encouraging us to resolve its effect.