THE MILKY WAY'S DARK MATTER HALO APPEARS TO BE LOPSIDED

Asymmetries are common in disk galaxies and are seen in lopsidedness (Baldwin et al. 1980; Rix & Zaritsky 1995; Bournaud et al. 2005; Saha et al. 2007), in warps (Garcia-Ruiz et al. 2002; Sanchez-Saavedra et al. 2003; Saha & Jog 2006), in the rotation curves of receding side and approaching side (Manabe & Miyamoto 1975; Swaters 1999), and in the distribution of the neutral hydrogen gas in galaxies (Richter & Sancisi 1994). These asymmetries in the dynamical phenomena can provide valuable information about the nature of the underlying dark matter potential. In fact, they have often been considered a reflection of the asymmetry in the dark matter distribution (Weinberg 1994; Jog 1997). However, one needs to be careful and systematically discard other possibilities such as dynamical instabilities as the origin of these asymmetries (e.g., Saha & Jog 2006). Our motivation in this paper is to check whether a nonaxisymmetric dark matter halo can explain the recently measured asymmetry in the atomic hydrogen gas distribution in our Galaxy.

In order to study the shape of the dark matter halo, it is best to look for tracers which lie mostly outside the main baryonic disk component. Examples of such tracers are the motion of satellites, neutral hydrogen (H i) gas in the outer region of the galaxy etc. In particular, over the last two decades it has become well known that the H i gas layer flares significantly beyond the optical disk of our Galaxy (Kulkarni et al. 1982; Knapp 1987; Wouterloot et al. 1990; Diplas & Savage 1991; Merrifield 1992; Nakanishi & Sofue 2003) and the recent Leiden/Argentine/Bonn (LAB) survey (Kalberla et al. 2005) of Galactic H i reveals the most comprehensive, uniformly sampled flaring map of the H i gas extended out to a very large radius from the Galactic center. Many external edge-on galaxies seem to show flaring in the thickness of the neutral hydrogen gas (Brinks & Burton 1984; Olling 1996a; Matthews & Wood 2003). The most likely reason for the flaring is that the total gravitational force acting perpendicular to the disk plane decreases with radius while the velocity dispersion of H i is observed to be nearly constant (Lewis 1984). The contribution to the total perpendicular gravitational force comes mainly from the stellar disk, gas, and the dark matter halo. Since the midplane density of the stars falls off rapidly compared to that of dark matter halo at large galactocentric radii (typically beyond the optical disk), the dark matter halo is expected to take over the major role in determining the vertical distribution of H i gas in the outer region. This makes the H i layer in the Galactic outskirts extremely sensitive to the distribution of dark matter and presently available neutral hydrogen gas (from the LAB survey) extended out to a very large radius from the Galactic center provides us an unique opportunity to examine the detailed nature of the dark matter distribution in the Milky Way.

This approach has been used to investigate the nature of dark matter halos by studying the thickness of the neutral hydrogen gas in our Galaxy (Olling & Merrifield 1998, 2000, 2001; Narayan et al. 2005; Kalberla et al. 2007) and for M31 (Banerjee & Jog 2008). Of these, recent work by Narayan et al. (2005) based on the Wouterloot et al. (1990) H i data shows that the observed H i flaring of Milky Way is best fit by a spherical dark matter halo with density falling faster than the isothermal halo. On the other hand, the most recent and very detailed work by Kalberla et al. (2007), based on the LAB survey data, shows that the dark matter distribution in the Galaxy is rather complicated. They needed a massive extended dark matter halo, a self-gravitating dark matter disk, and a dark matter ring to explain the flaring in the H i gas in our Galaxy. Now, one of the most striking feature in the LAB survey is that the neutral hydrogen gas thickness map shows systematic north–south (hereafter N–S) asymmetry, where Galactic north refers to 0° < φ < 180° and Galactic south by 180° < φ < 360°, in the Galctocentric cylindrical polar coordinate system (R, φ, z). This asymmetry was previously seen in previous H i maps of the Milky Way (Henderson et al. 1982), but rarely remarked on. So it may be worth investigating if the Milky Way H i flaring can be explained by a nonaxisymmetric dark matter halo.

Such nonaxisymmetries in the collisionless dark matter halo are probably not uncommon. In the current cosmological paradigm, the Λ cold dark matter (ΛCDM) halos are formed by dissipationless gravitational collapse of the material associated with the peaks of the primordial density fluctuation field and then grow via mergers and accretion in a highly nonlinear fashion. The resulting dynamical structures of these halos can be highly asymmetric. In fact, the ΛCDM halos formed in the recent Millennium Simulation (Springel et al. 2005) show asymmetry in their mass distribution (Gao & White 2006).

In the present study, we derive, numerically, the flaring in the thickness of the neutral hydrogen gas using self-consistent model for the Galaxy (Narayan & Jog 2002; Narayan et al. 2005) including a nonaxisymmetric dark matter halo and disk. We take into account the self-gravity of the gas. Our analysis produces a nonaxisymmetric flaring curve for the gas and in this respect our study is probably different from all the previous studies which tried to derive H i flaring in spiral galaxies. Based on our analysis, we show that an elliptically perturbed (readers are referred to Section 2.2 for details) lopsided dark matter halo with density falling faster than that of an isothermal halo can explain the observed N–S asymmetry in the H i thickness map in the Galaxy.

The paper is organized in the following order. Section 2 describes the formulation of the problem and models of the dark matter halo. In Section 3, we present the definition of the thickness of the H i gas and asymmetry measurements. Method and input parameters are discussed in Section 4, while Section 5 describes the results and various possible models of dark matter halo. Comparison of different models is done in Section 6. Section 7 describes a comparison with previous works. Discussion and conclusions are in Sections 8 and 9, respectively.

We studied the dynamics of an H i gas disk under vertical hydrostatic equilibrium in a generalized nonaxisymmetric dark matter halo potential. The H i disk is gravitationally coupled to the stellar counterpart. Basically, we are going to study the equilibrium vertical structure of a two-component star–gas system in an asymmetric potential. We used the cylindrical polar coordinate system (R, φ, z) suitable for the disk geometry. Presuming that the H i gas is in hydrostatic equilibrium, the vertical dynamics of each component under the force field due to the surrounding dark matter halo can be described by coupling the Poisson equation and the equation of hydrostatic equilibrium along the normal to the midplane for each component.

The Poisson equation for the system can be written as

$$\begin{equation} \frac{1}{R}\frac{\partial }{\partial R}\left(R\frac{\partial {\Phi _{t}}}{\partial R}\right) + \frac{1}{R^2} \frac{\partial ^2 \Phi _{t}}{\partial {\varphi ^2}} + \frac{\partial ^2 \Phi _{t}}{\partial z^2} = 4\pi G \left(\sum _{i}^{2} \rho _{i} + \rho _{h} \right), \end{equation} \tag{ 1 }$$

where Φ_t is the total potential due to the stars, H i gas, and dark matter halo, ρ_i with i = 1 to 2 denotes the mass density for the stellar and H i components, and ρ_h is the density of dark matter halo.

Next, we made a comparative estimate of the Poisson equation's radial and azimuthal terms (which we denote as T_R and T_φ, respectively). We assumed that the potential of the Galaxy is asymmetric by a small quantity _p and the dominant asymmetry is in the form of a lopsidedness (i.e., the isopotential contours are distributed according to cos φ, corresponding to m = 1 azimuthal wavenumber). We then wrote the total potential, Φ_t, in the simplistic form (Rix & Zaritsky 1995; Jog 1997):

$$\begin{equation} \Phi _{t}(R,\varphi,z)=\Phi _{0}(R,z)(1 + \epsilon _{p}\cos (\varphi -\varphi _{p})), \end{equation} \tag{ 2 }$$

where Φ₀(R, z) is the axisymmetric part of the total potential and φ_p is the constant phase factor. Then, the azimuthal term can be written as

$$\begin{equation} T_{\varphi }= \frac{\Phi _{0}(R,z)}{R^{2}} \epsilon _{p}\cos (\varphi - \varphi _{p}), \end{equation} \tag{ 3 }$$

and the radial term as,

$$\begin{equation} T_{R} \equiv \frac{1}{R}\frac{\partial V_{c}^2}{\partial R}(1+ \epsilon _{p}\cos (\varphi - \varphi _{p})), \end{equation} \tag{ 4 }$$

where the rotation velocity, $V_c =\sqrt{R \partial {\Phi _{0}(R,z)}/{\partial R}}$ is determined at the midplane (z = 0). To proceed further with the comparison, we need the form of the combined axisymmetric potential Φ₀(R, z). We assumed that the rotation curve in the outer region of the stellar disk is mainly dominated by the dark matter halo.

Let V_c ∼ C₀R^−α, where the index α determines the shape of the rotation curve and C₀ is some constant of proportionality. Note that α < 0, i.e., a negative α would produce rising rotation curves usually found in dwarf galaxies. In normal spiral galaxies α ⩾ 0. For example, α = 0 produces a flat rotation curve that corresponds to a screened isothermal dark matter halo (p = 1 in our notation, see Equation (12) in Section 2.2). In the other extreme, when α = 0.5 the rotation curve falls in a Keplerian fashion. The dark matter halos that produce an asymptotically Keplerian rotation curve are with indices 1.5 < p ⩽ 2 (for the present work, we restrict ourselves to p = 2). Therefore, we can say that the range of the indices in the rotation curves, 0 ⩽ α ⩽ 1, roughly maps to the range 1 ⩽ p ⩽ 2 in the density distribution of the dark matter halos. Now clearly, for an asymptotically flat rotation curve the radial term is exactly equal to zero at the disk midplane (z = 0), whereas the azimuthal term is given by T_φ ∼ log(R)/R²_pcos(φ − φ_p). For a very slowly falling rotation curve, the azimuthal term dominates over the radial term and with little algebraic manipulations it can be shown:

$$\begin{equation} T_{R} + T_{\varphi } \sim \left[ 1 + \frac{\epsilon _{p}}{\alpha ^2} \cos (\varphi -\varphi _{p})\right]T_{R}. \end{equation} \tag{ 5 }$$

Then, the Poisson equation for a disk embedded in a nonaxisymmetric dark matter halo with a slowly falling rotation curve (α ≠ 0 and α ≪ 1) can be written as

$$\begin{equation} \frac{\partial ^2 \Phi _{t}}{\partial z^2} \simeq 4\pi G \left(\sum _{i}^{2} \rho _{i} + \rho _{h} \right) - \left[ 1 + \frac{\epsilon _{p}}{\alpha ^2} \cos (\varphi -\varphi _{p})\right]T_{R}. \end{equation} \tag{ 6 }$$

Since T_R is a negative quantity, the combined effect of the radial and azimuthal terms in the Poisson equation (Equation (6)) is either to increase or decrease the vertical oscillation frequency of the disk depending upon the orientation of the dark matter halo and the values of α and _p. In other words, this term may either try to confine the gas more toward the disk midplane or help flaring. However, the contribution from the combined radial and azimuthal terms to the thickness of H i in the outer region is not significant compared to that due to the first term in Equation (6). The second term on the right-hand side (rhs) of Equation (6) contributes to ∼10% in the H i thickness in the outer region for a very slowly falling rotation curve. So in the zeroth-order approximation, we can safely say that the thickness of the disk components is largely determined by the vertical pull due to the first term on the rhs of the Poisson Equation (6).

The vertical hydrostatic equilibrium equation under the total potential for each component can be written as (Rohlfs 1977; Binney & Tremaine 1987)

$$\begin{equation} \sigma _{zi}^2 \frac{\partial \ln \rho _{i}}{\partial z} = -\frac{\partial \Phi _{t}}{\partial z}. \end{equation} \tag{ 7 }$$

In the above equation, σ_zi denotes the vertical velocity dispersion of the ith disk component in the problem. On combining the above Equations (6) and (7), the vertical equilibrium of each component in the disk under the dark matter halo potential is

$$\begin{eqnarray} \frac{\partial }{\partial z}\left(\sigma _{zi}^2 \frac{\partial \ln \rho _{i}}{\partial z} \right) \!\!\!\!&=&\!\!\!\! {-} 4\pi G \left(\sum _{i}^{2} \rho _{i} + \rho _{h} \right) \nonumber \\ &&\!\!\!\! + \left[ 1 + \frac{\epsilon _{p}}{\alpha ^2} \cos (\varphi -\varphi _{p})\right]T_{R}. \end{eqnarray} \tag{ 8 }$$

2.1. An Analytic Approximation to the Thickness

An analytic approximation for the thickness can be obtained near the disk midplane by direct integration of the vertical equilibrium equation. This is a useful guide to the numerical integration of Equation (8).

On integrating Equation (7) of the vertical equilibrium equation and using the boundary condition at z = 0

$$\begin{eqnarray*} &&\rho _{i}(R,\varphi,z) = \rho _{\rm mid}^{i}(R,\varphi,0), \end{eqnarray*}$$

where ρⁱ_mid(R, φ, 0) is the midplane volume density of each component, we get

$$\begin{equation} \sigma _{zi}^2 \ln { \left[\frac{\rho _{i}(R,\varphi,z)}{\rho _{\rm mid}^{i}(R,\varphi,0)} \right]}+ \Phi _{t}(R,\varphi,z) - \Phi _{t}(R,\varphi,0) = 0. \end{equation} \tag{ 9 }$$

Near the disk midplane. the potential along the vertical direction will not be very different than that of the disk midplane. Using Taylor's expansion along the vertical direction and the vertical equilibrium of each component, it can be shown that

$$\begin{equation} \rho _{i}(R,\varphi,z) = \rho _{\rm mid}^{i}(R,\varphi,0)\times e^{-\frac{z^{2}}{2 H_{i}^2(R,\varphi)}}. \end{equation} \tag{ 10 }$$

Hence, the density distribution under vertical equilibrium follows a Gaussian near the disk midplane with nonaxisymmetric thickness. Here, H_i(R, φ) = σ_zi/ν(R, φ) and the vertical frequency ν is given by

$$\begin{equation} \nu (R,\varphi) = \sqrt{ 4 \pi G \left(\sum _{i}^{2} \rho _{\rm mid}^{i}(R,\varphi) + \rho _{h}(R,\varphi,0) \right)}. \end{equation} \tag{ 11 }$$

In deriving the above form, we used Equation (6) with T_R = 0, since the thickness of the gas is largely determined by the vertical pull of the matter toward the disk midplane. The reader should bear in mind that the above analytical formula for the thickness (H_i) is not valid for high z.

2.2. Models of Nonaxisymmetric Dark Matter Halo

We consider here for the first time a seven-parameter dark matter halo model for describing the neutral hydrogen gas thickness in the Milky Way. The density profile of the nonaxisymmetric dark matter halo is given by

$$\begin{equation} \rho _{h}(R,\varphi,z) = \frac{\rho _{0}(q)}{\big[ 1 + \frac{\Lambda _{\varphi }^2}{R_{c}^2(q)}\big]^p}, \end{equation} \tag{ 12 }$$

where ρ₀ is the central mass density of the halo, R_c is the core radius, q determines the oblateness of the halo, and p is the index determining the nature of the density profile. The index p = 1 denotes a softened isothermal dark matter halo with density falling as R⁻² at large radii. The name "softened isothermal halo" is derived from the singular isothermal halo by adding a suitable core radius; these softened halos are also called screened or pseudoisothermal halos. Since the mass M(R) within a radius R is proportional to R, the p = 1 halo produces an asymptotically flat rotation curve. For p = 1.5, the density ρ_h ∝ R⁻³ at large radii resembling the Navarro–Frenk–White (NFW) halo profile (Navarro et al. 1996). The mass M(R) within a radius R of such a p = 1.5 halo is proportional to log R and goes to infinity as R goes to infinity but much more gradually than the p = 1 halo. Whereas p = 2 denotes a perfect ellipsoid dark matter halo with density falling like R⁻⁴ at large radii leading to essentially a finite mass halo. So the asymptotic rotation curve would be like a Keplerian one. The mass profile of the p = 2 axisymmetric dark matter halo can be written as

$$\begin{equation} M_{2}(R) \sim 2\pi \rho _{0} R_{c}^3 q \left[\tan ^{-1}(R/R_{c}) - \frac{R/R_c}{1 + R^2/R_{c}^2}\right]. \end{equation} \tag{ 13 }$$

In the above Equation (12)

$$\begin{eqnarray} && \Lambda _{\varphi }^2 = \big(R_{\varphi }^2 + z^2/q^2 \big); \nonumber \\&&[-5pt]\\[-5pt] && R_{\varphi } = R \big(1 - \epsilon _{h}^{1}\cos {(\varphi - \varphi _{h})} + \epsilon _{h}^{2}\cos \,{2(\varphi - \varphi _{h})} \big), \nonumber \end{eqnarray} \tag{ 14 }$$

where Λ_φ represents the surface of the concentric ellipsoid with lopsided (m = 1) distribution superposed with second harmonics (m = 2). The degree of lopsidedness in the dark matter distribution is ¹_h and the degree of second harmonic is determined by ²_h. φ_h denotes the phase of the asymmetric halo with respect to the Galactic axes. It is very hard to say without a detailed stability analysis which harmonics (m = 1 or m = 2) will dominate in the dark matter halo.

Since the higher harmonics (m > 2 or 3) perturbations will be associated with smaller spatial scales compared to the first or second harmonics, in a collisionless dark matter halo such small-scale higher harmonics are most likely to be Landau damped and large-scale perturbations would be weakly damped (Weinberg 1994). Our expectation is that the dark matter halos are likely to be dominated by the first two harmonics: a lopsided (m = 1) perturbation and an elliptical (m = 2) perturbation. The effects of a lopsided halo perturbation and m = 2 component in the halo perturbation onto the disk dynamics were studied by Jog (1999, 2000). Note that with ²_h = 0, we end up with a purely lopsided dark matter halo. With all _h = 0, the dark matter halo becomes the usual four-parameter axisymmetric halo used in previous studies (Narayan et al. 2005; de Zeeuw & Pfenniger 1988; Becquaert & Combes 1997). The parameter φ_h contains important information as to how the dark matter halo is oriented with respect to the Galactic axes. It is important to note that the first two harmonics in the dark matter halo are out of phase. Writing the density distribution of the halo in the above form results in freedom to investigate a wide variety of dark matter halo potentials, while simpler to use than the triaxial ones.

By solving the coupled Poisson equation and the vertical hydrostatic equilibrium equation (Equation (8)), we obtain the volume density of the atomic hydrogen gas ρ₂(R, φ, z). To solve Equation (8), we used the midplane volume density and its dispersion of each component as inputs (see Section 4). We considered an exponential surface density distribution for the stellar disk with a mild lopsidedness into it.

$$\begin{equation} \Sigma (R, \varphi) = \Sigma _{0} e^{-R/R_d}(1 + \epsilon _{s} \cos (\varphi - \varphi _{h})), \end{equation} \tag{ 15 }$$

where Σ₀ and R_d are the central surface density and scale length of the stellar mass distribution, respectively. The value of _s, in principle, can be determined self-consistently by calculating the response of the stellar disk to the nonaxisymmetric dark matter halo considered in Section 2.2 but the task is beyond the scope of this paper and not necessary for obtaining a correct zeroth-order model. Instead, we use reasonable numerical values of _s and calculate its effect on the H i thickness in the region of our interest (R ⩾ 16 kpc). We have checked that the gas thickness remains almost unchanged as the value of _s is reduced from 0.15 to 0.0 (corresponding to an axisymmetric stellar disk). We have also checked if _s alone can account for the observed systematic N–S asymmetry (prominent beyond ∼16 kpc) in the H i thickness map. We found that even unrealistically high values of _s alone cannot account for such high asymmetry in the H i flaring. Note that the contribution of the axisymmetric stellar disk itself almost drops to zero beyond ∼16 kpc compared to that due to the dark matter halo. So in the present calculation, we considered only the axisymmetric stellar disk in deriving the H i thickness beyond ∼16 kpc. One important assumption in solving Equation (8) is that the velocity dispersion remains constant (isothermal approximation) along the vertical direction.

The thickness d(R, φ) of the gas is determined by using the second moment of the volume density distribution and is given by the following relation:

$$\begin{equation} d^{2}(R,\varphi) = \frac{\int _{-\infty }^{\infty }{z^{2}\rho _{2}(R,\varphi,z)\, dz}}{\int _{-\infty }^{\infty }{\rho _{2}(R,\varphi,z)\, dz}}. \end{equation} \tag{ 16 }$$

d(R, φ) gives the radial variation of the thickness along a particular azimuthal direction (φ) in the disk. Using this thickness map, we can examine the degree of asymmetry in the gas thickness distribution.

We note that the thickness of the gas layer in this case is different from the method used by Levine et al. (2006a) but gives qualitatively similar result. Kalberla et al. (2007) calculate the gas scale height yet another way, and get results which differ both from the second moment of the distribution (Equation (16)) and from Levine et al. (2006a). All this says is that because of the exclusion of certain areas (especially 90° ⩽ φ ⩽ 110°), and the possible effects of optical depth, different methods of calculating the scale height can give rise to different numerical values. The important point is that all methods show values of scale height that are ∼2–2.5 times larger in the northern part of the Galaxy (l = 0°–180°) than the south (l = 180°–360°) and it is this difference that we trying to model and explain.

Let A(φ_k) be the area under the thickness curve along a particular azimuthal angle (φ_k):

$$\begin{equation} A(\varphi _{k}) = \int _{R_{\min }}^{R_{\max }}{d(R,\varphi _{k})\, dR}. \end{equation} \tag{ 17 }$$

In the above equation, R_min and R_max represent the initial and the final radius, respectively in the available observation of the gas thickness. The degree of asymmetry in the thickness distribution, presuming that we are considering a smooth curve, can then be described by

$$\begin{equation} \eta (\varphi _{k}) = \frac{\left|{A(\varphi _{k}) - A(\varphi _{k} +\pi)}\right|}{ A(\varphi _{k}) + A(\varphi _{k} +\pi)}. \end{equation} \tag{ 18 }$$

The range of η is 0 ⩽ η ⩽ 1, with η = 0 denoting a symmetric distribution. On the other hand, η = 1 denotes a highly asymmetric gas/star distribution leading to a one-sided flaring in the galaxy. It is unlikely that η ≈ 1, because this would mean that the velocity dispersion of the gas in one-half of the galaxy is almost zero (cold component) compared to the other half. This may lead to an instability in the disk. Considering φ_k = 90°, the above relation would produce the numerical value of the observed N–S asymmetry in the thickness map.

The present study focuses on deriving the thickness map of the neutral hydrogen gas (H i) in the very outer region (R ⩾ 16 kpc) of the Galactic disk. Since the various baryonic components of the stars and H i dominate the mass in this region, we neglect the effect of molecular hydrogen gas (H₂) on the thickness distribution of H i. Equation (8), which describes the zeroth-order vertical equilibrium under a nonaxisymmetric dark matter distribution, represents two coupled differential equations for the two disk components: H i and stars. The vertical density distribution for each component, responding to the total potential due to the disk and the dark matter halo, was solved numerically as an initial value problem using the fourth-order Runge–Kutta method of integration (Press et al. 1994, Chapter 6). The details of this method are presented in Narayan & Jog (2002) and Narayan et al. (2005). Because of the underlying nonaxisymmetry, Equation (8) is solved for each azimuthal direction along the Galactocentric radius. Along a particular azimuth, at a particular radius, we calculate the second moment of the vertical density distribution for each component according to Equation (16) and call it the thickness of that component. Repeating this procedure for regular intervals in the azimuth along the galactocentric radius gives us the thickness map of the neutral hydrogen gas layer in the Galaxy.

The primary input parameters needed for the method are the midplane volume density and the vertical velocity dispersion for each disk component. The H i midplane volume density is obtained from the LAB survey (Kalberla et al. 2005; Levine et al. 2006a). For the stellar disk, we used its surface density according to Equation (15) for which we need to know the central surface density and scale length. Note that Equation (15) represents a lopsided stellar disk and as discussed in Section 3, the value of _s does not effect much the H i thickness in the very outer region of the Galaxy. So we considered effectively an axisymmetric stellar disk in the region of our interest. Using the following measured/inferred quantities: the stellar surface density at the solar region Σ_⊙, the disk scale length R_d, and the distance of Sun from the Galactic center R_☉, we can derive the central surface density Σ₀ and in fact the surface density at any radius. We use Σ_☉ = 45 M_☉ pc⁻² which is consistent with 48 ± 9 M_☉ pc⁻² obtained by Kuijken & Gilmore (1991) and 52 ± 13 M_☉ pc⁻² obtained by Flynn & Fuchs (1994) for the total surface density, after the gas density is subtracted. We use the IAU-recommended value for R_☉ (=8.5 kpc) and this was also used to determine the observed thickness map for H i gas (Levine et al. 2006a). The scale length R_d was set equal to 3.2 kpc (Mera et al. 1998) in accordance with the recent determinations of smaller disk scale length for our Galaxy.

4.1. Stellar and H i Velocity Dispersion

4.2. Model Rotation Curves

We model the rotation curve for the Galaxy using a bulge, an exponential disk for the stars and gas and a dark matter halo. Modeling the observed rotation curve is a nontrivial task given that there are various kinds of uncertainties and turns out to be a nonlinear regression problem in a multidimensional parameter space. Naturally, such modeling would suffer from uniqueness problem. So we aim here to reproduce the main features in the observed rotation curve and try to generate such rotation curves for which the circular speed lies in the range determined by the relation Θ° = (27 ± 2.5)R° km s⁻¹ due to Kerr & Lynden-Bell (1986). Since we use the IAU-recommended value for R° = 8.5 kpc, the above range implies V_c = 230 ± 21 km s⁻¹ at the solar radius. Keeping these constraints in mind, we proceed to derive the model rotation curves in the following way. We assume that the disk and the bulge are aligned with the symmetry axis of the dark matter halo. We also assume the virial equilibrium in the Galaxy. Then, the square of the total circular velocity in the disk midplane (z = 0) can be written as

$$\begin{equation} V_{c}^2 = V_{\rm bulge}^2 + V_{\rm stars}^2 + V_{\rm gas}^2 + V_{dmh}^2, \end{equation} \tag{ 19 }$$

where we adopt a Plummer–Kuzmin bulge model to derive the bulge contribution to the rotation curve. The density profile of the bulge is given by the following formula (Binney & Tremaine 1987):

$$\begin{equation} \rho _{b}(R)=\frac{3M_{b}}{4\pi R_{b}^{3}}\left(1+\frac{R^{2}}{R_{b}^{2}}\right)^{-5/2}, \end{equation} \tag{ 20 }$$

where R_b is the bulge scale length and M_b is the total bulge mass. We have used same values of these two parameters in all of our model rotation curves, e.g., R_b = 2.5 kpc and M_b = 2.8 × 10¹⁰ M_☉ (see Blum 1995). The inclusion of this simple spherical bulge reproduces a reasonable looking rotation curve for the Galaxy. However, the contribution of the bulge to the H i thickness in the region of interest (R ⩾ 16 kpc) is almost negligible, the bulge contributes to ∼5% to the overall gas thickness.

The disk contribution to the circular speed is derived using the potential for the exponential mass distribution ($\Sigma (R)=\Sigma _{0}e^{-R/R_d}$ where values of the parameters are mentioned in Section 4) and similarly for the various dark matter halo models for which the parameters are mentioned in the appropriate places below. The gas contribution to the total circular speed is derived from the following exponential distribution for the H i beyond 14 kpc (Levine et al. 2006a):

$$\begin{eqnarray*} &&\Sigma _{{\rm H}\,\hbox{\scriptsize {\sc i}}}(R)=4.5\times e^{[-(R - 14\,{\rm kpc})/4.3\,{\rm kpc}]} \end{eqnarray*}$$

Below 14 kpc, we use a constant surface density for the H i gas (as observed) to derive its contribution to the rotation curve. Although there is a noteable difference of our model rotation curves with the observed one (Brand & Blitz 1993), the slowly falling rotation curves due to the p = 2 halo model follow the trend found in the Milky Way's rotation curve in the very recent analysis of Xue et al. (2008) based on the SDSS data.

There is a clear N–S asymmetry in the thickness map of the H i gas in our Galaxy. It is not obvious what would have caused such asymmetry in the H i gas distribution. What is the underlying nature of this asymmetry? It remains to be shown whether this is purely a gas dynamical effect or reflects some gravitational effect. Below we describe a step by step analysis of the cause of this asymmetry, which gradually reveals the nature of dark matter halo in our Galaxy. We used averaged flaring data over the north (0° ⩽ϕ⩽ 180°) and south (180° ⩽ϕ⩽ 360°), respectively, excluding about 15° region about the Sun–Galactic center line. The vertical volume density distribution of the gas is derived using the rotation curve due to Brand & Blitz (1993). The thickness of the gas is then derived by taking the second moment of the density distribution after the thickness filter has been applied. In this respect, the thickness measurement of the gas is different from Levine et al. (2006a) who use the tail integration method and it is also different from Kalberla et al. (2007) who use half-width at half-maximum (HWHM) for the half-thickness of the gas later. In all our analyses below, we use the same thickness map for the gas, although it is understood that a different rotation curve would indeed produce different thickness map for the gas. However, on doing an error analysis we find that the error in calculating the distance due to an error in the rotation velocity is small (∼10% at large galactocentric distances) to produce an appreciable change in the observed thickness of the gas. The readers are referred to Section 8 for a discussion on the dependence of the thickness data on the rotation curve. The observed N–S asymmetry denotes the average asymmetry in the thickness map of the H i gas and according to Equation (18), it is given by η_obs = 0.262.

5.1. Models of Axisymmetric Dark Matter Halo

We first considered a simple axisymmetric dark matter halo model (with all _h = 0 in Equation (14)) for the Galaxy and used different velocity dispersions for the H i in the two halves to see if the observed N–S asymmetry could be reproduced.

5.1.1. p = 1, Softened Isothermal Halo

The softened isothermal dark matter halo produces naturally the asymptotically flat rotation curve in a spiral galaxy. It has been shown previously based on the Wouterloot et al. (1990) data that a p = 1 softened isothermal halo of any shape (oblate or prolate) cannot explain the observation (Narayan et al. 2005). In the present study, we found the same trend, so we considered a nearly spherical halo to begin with for further investigation. We adopt the parameters of the p = 1 isothermal dark matter halo from the mass model of our Galaxy based on microlensing observation by Mera et al. (1998). We found that the p = 1 halo with core radius R_c = 5 kpc, central density ρ₀ = 0.035 M_☉ pc⁻³, which gives rise to an asymptotically flat rotation curve with terminal velocity of 220 km s⁻¹ (Brand & Blitz 1993; see bottom right of Figure 1), and with $\sigma _{{\rm H}\,\hbox{\scriptsize {\sc i}}}=9$ km s⁻¹ cannot explain the present observation (LAB data). The fact that p = 1 softened isothermal halo cannot explain the flaring in the gas thickness distribution has already been verified by Narayan et al. (2005) in the case of Wouterloot et al. (1990) data and more recently by Kalberla et al. (2007) based on the LAB survey data. Now, there are two aspects of the H i thickness map derived from the LAB survey data: the radial variation of the gas thickness in the north and the south and the prominent asymmetry between the north and thesouth as mentioned already. An axisymmetric p = 1 halo with the above-mentioned parameters (giving rise to V_c = 220 km s⁻¹) and $\sigma _{{\rm H}\,\hbox{\scriptsize {\sc i}}} = 9$ km s⁻¹ appears to be quite massive and hence difficult to reproduce the present observation. By inspection of Figure 1 (top left), it is easy to check that there has to be an order of magnitude decrease in the dark matter midplane density to reproduce the observed gas thickness at 30 kpc. Moreover, the striking difference between the two slopes makes it clear that even a linearly increasing velocity dispersion with radial distance (although unphysical) cannot fit the observation. Obviously, an axisymmetric p = 1 halo alone cannot explain the other important aspect of the thickness map, namely, the observed N–S asymmetry. This fact leads us to explore another possibility of a model of an axisymmetric p = 1 halo accompanied by a nonaxisymmetric distribution of H i velocity dispersion. It can be explained in a simple way that a different H i velocity dispersion on both sides of the hemispheres is not going to improve the situation either. In a naive theory, the thickness of H i gas can be written as $h_{{\rm H}\,\hbox{\scriptsize {\sc i}}} \propto \sigma _{{\rm H}\,\hbox{\scriptsize {\sc i}}}/\sqrt{G \rho _{\rm mid}}$. Using this simple formula, it is easy to see that in order to increase the thickness by a factor of 2 at a particular radius, one needs to raise the dispersion by a factor of ∼2 (because ρ_mid remains the same). At R = 30 kpc, the thickness of H i gas in the northern hemisphere is roughly a factor of 2 more than that in the southern hemisphere and to explain this asymmetry based on purely gas dynamical effect and the axisymmetric p = 1 isothermal halo (giving rise to a flat rotation curve), one needs to have $\sigma _{{\rm H}\,\hbox{\scriptsize {\sc i}}} \sim 20$ km s⁻¹ which is unlikely according to the standard models of interstellar medium (ISM) of our Galaxy.

Zoom In Zoom Out Reset image size

5.1.2. p = 1.5, an NFW Type Halo

At large distances (R ≫ R_c) from the center, p = 1.5 halo resembles an NFW profile (Navarro et al. 1996) for the dark matter halo. In our model of vertical equilibrium, this seems to be preferable compared to the p = 1 softened isothermal dark matter halo. Because of the density falling faster than the p = 1 isothermal halo, the resulting rotation curve also starts falling beyond about 10 kpc. We use Equation (2.96b) of Binney & Tremaine (1987) to generate the rotation curve of the p = 1.5 halo. The axisymmetric p = 1.5 halo with R_c = 8 kpc and ρ₀ = 0.025 M_☉ pc⁻³ produces a reasonable rotation curve (see bottom right of Figure 1). At R_☉, the rotation velocity is 223 km s⁻¹ and at 2 R_☉, V_c = 216 km s⁻¹. At 25 kpc, the difference in rotation velocities between the p = 1 and p = 1.5 halo is approximately 20 km s⁻¹. In the top right panel of Figure 1, we show the half-thickness of H i gas due to the p = 1.5 halo model considered here. The axisymmetric halo model with $\sigma _{{\rm H}\,\hbox{\scriptsize {\sc i}}}=9.2$ km s⁻¹ fits well the observation in the southern part beyond about 15 kpc. However, the same model does not fit the thickness curve at all in the northern halves. The solid line shows the curve with H i velocity dispersion ∼12.2 km s⁻¹ and yet does not give a good fit. An increase in the H i velocity dispersion beyond 12 km s⁻¹ would force the model curve to intersect the observed one only at a single point in the north. This again demonstrates that the N–S asymmetry in the thickness map of the H i is probably not due to a gas dynamical effect, rather it arises due to some gravitational effect.

5.1.3. p = 2, a Perfect Ellipsoidal Halo

Previous studies by Narayan et al. (2005) have investigated the axisymmetric p = 2 halo in considerable detail to explain the H i thickness in the Galaxy. Their p = 2 halo provides a good fit to the Wouterloot et al. (1990) flaring data. In the present study, we also considered an axisymmetric p = 2 halo and see if different H i velocity dispersions can explain the N–S asymmetry. We consider a core radius R_c = 9.4 kpc and central density ρ₀ = 0.035 M_☉ pc⁻³ which produces a reasonable rotation curve (see Figure 1) within the uncertainties in the observation out to 30 kpc where the difference between the p = 1 (flat rotation curve) and p = 2 rotation curves is ∼25 km s⁻¹. Note that the rotation curves for both p = 1.5 and p = 2 halos are slowly falling with radial distance. Recent analysis of the kinematics of a large number of blue horizontal-branch halo stars from the SDSS database by Xue et al. (2008) shows that the rotation curve of our Milky Way is actually falling slowly with the Galactocentric radius. So this direct observational analysis supports the falling trend in the rotation curve of Milky Way. In fact, the choice of the halo parameters for the p = 2 case makes the rotation curve more realistic as observed. At 25 kpc, the circular speed due to our p = 2 halo differs from the flat rotation curve only by ∼15 km s⁻¹ which is well within the observed error bars. Beyond about 15 kpc, with linearly decreasing H i velocity dispersion from 9.2 (at 10 kpc) to 8.0 (at 30 kpc) km s⁻¹, the p = 2 halo considered here gives a good fit to the observed data in the southern halves. However, even with $\sigma _{{\rm H}\,\hbox{\scriptsize {\sc i}}}=12.2$ or 14.2 km s⁻¹, the p = 2 halo does not fit well the observation in the north. With 14.2 km s⁻¹, the model curve overestimates the observation between 13 and 23 kpc and underestimates beyond 24 kpc. On the other hand with 12.2 km s⁻¹, the model underestimates the observation beyond 20 kpc. The solid curve in Figure 1 (bottom left panel) is with a constant $\sigma _{{\rm H}\,\hbox{\scriptsize {\sc i}}}=13.2$ km s⁻¹ in the northern halves. Compared to p = 1 and p = 1.5, an axisymmetric p = 2 halo with different H i velocity dispersion comes closer to the observation. We have also tried to use a linearly decreasing H i velocity dispersion from 13.2 or 12.2 km s⁻¹, but no good fit found. It becomes quite clear that the observed H i flaring in the north cannot be explained with any physically meaningful variation of the H i velocity dispersion.

5.2. Nonaxisymmetric Dark Matter Distribution

Based on the above three axisymmetric cases we have investigated so far, we conclude that an asymmetric state of ISM combined with axisymmetric dark matter potential is not able to explain the observed asymmetry in the thickness map indicating that the observed asymmetry is probably not due to a gas dynamical effect. The large-scale asymmetry in the N–S thickness map is instead likely to be gravitational in origin. Beyond 16 kpc, the stellar contribution is not significant compared to the dark matter, which implies that the N–S asymmetry is not likely due to the asymmetric stellar disk and also the observed asymmetry in the stellar disk is not very high in the Galaxy. This gives us a room to explore whether the N–S asymmetry originates due to an asymmetric dark matter halo. As explained in Section 1, such asymmetric dark matter halos are not uncommon even in the cosmological N-body simulations. By inspecting the observed H i thickness map, we guess that the N–S asymmetry is primarily lopsided (m = 1) in nature. A lopsided dark matter halo has been used to provide an explanation for the observed lopsidedness in the underlying stellar disk and for the asymmetry in the rotation curves of disk galaxies (Jog 1997, 2002). Below we test our hypothesis that the Milky Way's dark matter halo is lopsided. We show next that the observed nature of the H i distribution demands that the lopsided dark matter halo to be oriented with phase angle ϕ_h = 270°. In the present case, ϕ_h is the angle between the direction in which the dark matter distribution is elongated more on one side compared to the other and the Galactocentric coordinate axes. So the dark matter density maximum of the lopsided halo is along the southern direction (ϕ = 270°) and the density minimum is along the northern direction (ϕ = 90°). We test p = 1.5 and p = 2 halos as the likely candidates in our further investigation; p = 1 is not used because it fails to produce the observed flaring even in the south.

5.2.1. p = 1.5 Lopsided Halo

Having oriented the dark matter halo with phase ϕ_h = 270°, we found that an axisymmetric (_h = 0) p = 1.5 halo with H i velocity dispersion 9.2 km s⁻¹ fits quite well the observation in the southern halves (see Figure 1). However, a purely lopsided configuration (²_h = 0 in Equation (14)) with reasonable parameters cannot explain the observed asymmetry in the thickness map of H i gas. The parameter ranges we have tried are $\sigma _{{\rm H}\,\hbox{\scriptsize {\sc i}}}= 7\hbox{--}10$ km s⁻¹ and ¹_h = 0.02–0.4. We found that no combination of these parameters could provide a good fit to the observed data. Further, we tried to fit the combined model of the dark matter halo (given in Equation (12)) with lopsidedness as the dominant component. Again, for the above-mentioned parameter range no good fit to the observed data was found. Figure 2 shows a configuration of the dark matter halo in which both the first (m = 1) and second (m = 2) harmonics are equal in strength and out of phase with each other. Of course, such configuration of p = 1.5 halo with reasonable H i velocity dispersion (9.2 km s⁻¹) provides better fit to the observation than the axisymmetric configuration of p = 1.5 halo with high $\sigma _{{\rm H}\,\hbox{\scriptsize {\sc i}}}$.

Zoom In Zoom Out Reset image size

5.2.2. p = 2 Lopsided Halo

Since the axisymmetric p = 2 halo already seemed promising compared to others, we carried out a detailed analysis of the p = 2 lopsided halo here. We used the basic parameters like core radius and central density of the axisymmetric p = 2 configuration (see Section 5.1.3). In order to find a best-fit model of the Galaxy under a purely lopsided halo, we make a two-dimensional grid of two independent and free parameters $\big(\sigma _{{\rm H}\,\hbox{\scriptsize {\sc i}}}$, ¹_h) spanning a large dynamical range, i.e., $\big[\sigma _{{\rm H}\,\hbox{\scriptsize {\sc i}}}, \epsilon _{h}^{1} \big] = [(7$–12) km s⁻¹, (0.05–0.40)] for this halo model. We found that with $\sigma _{{\rm H}\,\hbox{\scriptsize {\sc i}}} \le 9.0$ km s⁻¹, a purely lopsided p = 2 halo is not sufficient to explain the observation. With some more exploration, we found that a purely lopsided p = 2 halo with ¹_h = 0.17 and $\sigma _{{\rm H}\,\hbox{\scriptsize {\sc i}}} = 11$ km s⁻¹ does fit the observation quite well (see Figure 3). We call this configuration model p2L. The rotation curve for this model is shown in the bottom panel of Figure 3. The difference between the rotation velocities in the north and the south is ∼20 km s⁻¹ for this model. The density contours of the dark matter halo are shown in the right panel of Figure 3. Note that the H i velocity dispersion for this model p2L is fairly high. So we further explored the parameter space consisting of ¹_h, ²_h, and $\sigma _{{\rm H}\,\hbox{\scriptsize {\sc i}}}$ for the p = 2 dark matter halo to find the best possible model to explain the observed radial variation of the thickness map and its N–S asymmetry.

Zoom In Zoom Out Reset image size

Based on the different values of H i velocity dispersion and different values of the ¹_h and ²_h, we construct three models (model A, model B, and model C) which give a very good fit to the averaged observation in both halves of the Galaxy and thereby explain the observed asymmetry.

Model A. We assumed the velocity dispersion of H i, $\sigma _{{\rm H}\,\hbox{\scriptsize {\sc i}}} = 8.5$ km s⁻¹, to be constant out to 30 kpc for simplicity. The value of the velocity dispersion is more reasonable one and it is very closed to what has been used by Kalberla et al. (2007) in their best-fit model. Once $\sigma _{{\rm H}\,\hbox{\scriptsize {\sc i}}}$ is fixed, we are essentially left with again a two-dimensional grid of free parameters ¹_h and ²_h for the dark matter halo. The parameter range explored is [¹_h, ²_h] = [(0.05–0.40), (0.05–0.40)]. We found that the inclusion of the second harmonic component in the dark matter halo dramatically improved the quality of the fit compared to all the previously explored cases. The best-fit model parameters, for the above-mentioned $\sigma _{{\rm H}\,\hbox{\scriptsize {\sc i}}}$, are ¹_h = 0.20 and ²_h = 0.18. Such configuration of the dark matter halo which is indeed lopsided, reproduced the N–S asymmetry quite well in the Galaxy (see Figure 4). The rotation curves for the model in the two halves of the Galaxy are shown in the right panel of Figure 4. The circular speeds at the solar radius are quite comparable to the p = 1 flat rotation curve within the uncertainties. At 25 kpc, the north and south rotation curves differ from each other by ∼19 km⁻¹, the value in the north being ∼209 km⁻¹. The density contours of the model dark matter halo are shown in the bottom panel of Figure 4. Because of the presence of the second harmonics (m = 2), the contours are elongated in the east–west direction of the Galaxy too. It is quite clear now that a simple model that consists only of a lopsided dark matter halo does not reproduce the observed asymmetry in the thickness map of the atomic hydrogen gas in our Milky Way. The lopsided dark matter halo with some amount of second harmonic (m = 2) superposed onto it appears to give the best fit to the data compared to the purely lopsided dark matter halo. We call the emerging picture of dark matter halo as elliptically perturbed lopsided dark matter halo.

Zoom In Zoom Out Reset image size

Model B. Here, we consider the velocity dispersion of H i, $\sigma _{{\rm H}\,\hbox{\scriptsize {\sc i}}} = 9$ km s⁻¹ slightly higher compared to model A and again flat out to 30 kpc. The increased velocity dispersion improves the fit by a very small amount below 16 kpc. With the same halo parameters as in model A, the fit is not good beyond 16 kpc. So we reduced the strength of the second harmonic (²_h) and when the second harmonic is decreased to ∼2/3 of the first harmonic (¹_h) we again recovered a good fit to the data. The resulting model halo with parameters ¹_h = 0.2 and ²_h = 0.14 reproduced the N–S asymmetry in the Galaxy quite well (see Figure 5). The rotation curves for the model in the two halves of the Galaxy are shown in the right panel of Figure 5. At 25 kpc, the rotation speed in the north is ∼191 km⁻¹ and in the south it is ∼212 km⁻¹. It appears that the rotation curve in the south is more closer to the flat rotation curve (with V_c = 220 km⁻¹). The density contours of the model dark matter halo are shown in the bottom panel of Figure 5. Due to the presence of a small component of second harmonic, the contours are elongated by a small amount compared to model A in the east–west direction of the Galaxy. This dark matter halo can be considered as a dominantly lopsided halo. This model again confirmed our earlier findings from model A.

Zoom In Zoom Out Reset image size

Model C. We constructed this model primarily to see if we can explain the observation also in the region 10 < R < 16 kpc. We found that the midplane density distributions of H i and stars are almost same in the radial range 10 ⩽ R ⩽ 16 kpc (see Figure 8) and interestingly enough the half-thickness of H i in this radial range is almost the same as that of stars. Another important fact to be noted is that the H i half-thickness is roughly constant in this radial range (10 ⩽ R ⩽ 16 kpc). Now, we know that the thickness of H i ∼$\sigma _{{\rm H}\,\hbox{\scriptsize {\sc i}}}/\sqrt{(}G \rho _{\rm mid})$, where ρ_mid is the total midplane volume density. So we vary the H i velocity dispersion in the following manner in this radial range:

$$\begin{eqnarray} \sigma _{{\rm H}\,\hbox{\scriptsize {\sc i}}} & = & 55.6 e^{-R/{2 R_{\sigma }}}\;{\rm km\;s^{-1}}\;{\rm for}\; 10 \le R \le 16\;{\rm kpc} \nonumber \\ & = & 8.5\;{\rm km\;s^{-1}}\;R \ge 16\;{\rm kpc}. \end{eqnarray} \tag{ 21 }$$

From 16 kpc (∼5 R_d) onward $\sigma _{{\rm H}\,\hbox{\scriptsize {\sc i}}}=8.5$ km s⁻¹ and stays flat out to 30 kpc. We consider R_σ = 4370 pc because a similar value in the case of stellar disk can produce a flat thickness (Lewis & Freeman 1989) in the Galaxy.

With the H i velocity dispersion varied in the above fashion (Equation (21)) and ¹_h = 0.2 and ²_h = 0.18 (same as model A), we got an excellent fit to the observation (see Figure 6). Model C provides the best fit to the observed data and reproduces the observed asymmetry in the thickness map of H i. In terms of halo asymmetry parameters, model C is exactly equal to the model A and so are the rotation curves. The density contours are the same as in Figure 4. So model C confirms that our best-fit dark matter halo model and again it is a lopsided dark matter halo with some elliptical perturbation.

Zoom In Zoom Out Reset image size

6.1. Measurement of Asymmetry

Given the radial variation of the H i thickness in different azimuths (ϕ), we can quantify the degree of observed asymmetry by η (see Equation (18)) in the two halves of the Galaxy. We used this parameter η as a measure of the underlying asymmetry in the thickness distribution and to differentiate between the above models which seem to reproduce the observation quite well. The value of η_obs is computed by assuming φ_k = 90° as discussed at the end of Section 3. Just by comparing the observed asymmetry (η_obs) and the model predicted asymmetry (η_model), we find that model p2L is the best match to the observation (see Table 1) but the velocity dispersion used in this model is quite high. Our next best model is model C based on the comparison of asymmetry and this model also has a very high velocity dispersion between 10 and 16 kpc. On the other hand,model A is probably more reasonable because the H i velocity dispersion is not very high. The H i dispersion of model A is almost the same as that used recently by Kalberla et al. (2007); their solution is close to a single component model like we have considered with an effective constant velocity dispersion of 8.3 km s⁻¹. In all the four models of p = 2 dark matter halo one thing is common and it is the parameter ¹_h ∼ 0.2, the degree of lopsidedness. So the different models are essentially builtup around the same basic configuration, a lopsided halo. To make the picture more clear, these models, being a function of two apparently independent parameters ($\sigma _{{\rm H}\,\hbox{\scriptsize {\sc i}}}$ and ²_h), can be thought of as a family of models of lopsided dark matter halos. From Table 1, amongst the first three models, it is clear that as ²_h increases, $\sigma _{{\rm H}\,\hbox{\scriptsize {\sc i}}}$ decreases to make the fit better. One interesting fact about these families of models is that a better fit to the observed data demands that the second harmonic be out of phase with the first one.

Table 1. Asymmetry and the Models of the p = 2 Lopsided Dark Matter Halos

Model	$\sigma _{{\rm H}\,\hbox{\scriptsize {\sc i}}}$ (km s−1)	ρ0 (M☉ pc−3)	Rc (kpc)	q	1h	2h	Σ1.1 (M☉ pc−2)	ηobs	ηmodel
p2L	11.0	0.035	9.4	0.95	0.17	0.0	72.0	0.262	0.261
A	8.5	0.035	9.4	0.95	0.20	0.18	80.0	0.262	0.265
B	9.0	0.035	9.4	0.95	0.20	0.14	78.0	0.262	0.271
C	Equation (21)	0.035	9.4	0.95	0.20	0.18	81.0	0.262	0.260

Download table as:  ASCIITypeset image

In Table 1, η_obs is obtained from the observed average thickness of H i on both the halves. η_model is from the fitted curves.

6.2. Constraints from the Surface Mass Density at Solar Radius

Apart from the constraint on the rotation curve due to Kerr & Lynden-Bell (1986), we provide here our estimations of total surface density at the solar radius for various models and compare them with other measurements from the literature. The value of the total surface density at the solar radius put a strong and important constraint on any mass model of our Galaxy. Based on K dwarfs, Kuijken & Gilmore (1991) provide us the total surface density within 1.1 kpc from the Galactic midplane Σ_1.1 = 71 ± 6 M_☉ pc⁻². By analyzing HIPPARCOS K giants, Holmberg & Flynn (2004) measure Σ_1.1 = 74 ± 6 M_☉ pc⁻². By modeling the H i thickness in the Galaxy, Kalberla et al. (2007) estimate Σ_1.1 = 79 ± 2 M_☉ pc⁻². Our best-fit models are in good agreement with all the previous measurements. For example, for the model, p2L Σ_1.1 = 72 M_☉ pc⁻²; for the model A, Σ_1.1 = 80 M_☉ pc⁻²; for the model B, Σ_1.1 = 78 M_☉ pc⁻², and for model C, Σ_1.1 = 81 M_☉ pc⁻². We used the total baryon mass surface density at the solar radius to be Σ_star+gas = 50 M_☉ pc⁻² which is within the limit (48 ± 8 M_☉ pc⁻² at R_☉) provided by Kuijken & Gilmore (1989). Within 800 pc from the Galactic plane, Holmberg & Flynn (2004) report a surface density of Σ_0.8 = 65 ± 6 M_☉ pc⁻² and from Kalberla et al. (2007) Σ_0.8 = 66.9 ± 2 M_☉ pc⁻². In this context, our model p2L gives Σ_0.8 = 65 M_☉ pc⁻² and from model A, we get Σ_0.8 = 71 M_☉ pc⁻². In Figure 7, we plot the total surface mass density (derived from modeling of the vertical hydrostatic equilibrium) of the Galaxy (including stars, gas, and dark matter) within 1.1 kpc from the Galactic midplane as function of radius. It shows clearly that the surface density is higher in the southern part of the Galaxy. In any case, our results on the surface density measurements in the solar neighborhood are in close agreement with previously quoted values or within the quoted error bars.

Zoom In Zoom Out Reset image size

6.3. Total Mass of the Galaxy

The most recent work which has extensively used the LAB survey H i data are Levine et al (2006a) and Kalberla et al (2007). Of these two, Kalberla et al. (2007) have studied in considerable detail the behavior of the H i thickness in the Galaxy and found that in order to explain the observation they needed beside a massive dark matter halo, a dark matter disk and a dark matter ring. Such configurations of the dark matter for our Galaxy are troublesome for the CDM paradigm.

Here, we would like to point out few similarities and differences of our method with the previous ones. First of all, one of our best-fit models, namely model A, uses a constant H i velocity dispersion of 8.5 km s⁻¹ which is close to what Kalberla et al. (2007) used for the their model namely 8.3 km s⁻¹. To make it more clear, although we do not use explicitly the different phases of ISM (e.g., cold neutral medium and warm neutral medium) to model the H i thickness map, the state of the ISM is roughly the same in our model as it is in Kalberla et al. (2007). In contrast to Kalberla et al. (2007), we consider nonaxisymmetric models of dark matter halo and also a mild lopsided stellar disk. From Figure 8, it is clear that the gas self-gravity is comparable to the stellar one and we take into account the full self-gravity of the H i gas in deriving the thickness map of the gas.

Zoom In Zoom Out Reset image size

There is a marked difference in the averaged flaring curves in the north and the south derived in our paper with that in Kalberla et al. (2007) which excludes a region 90° < φ < 110° in the north showing very high flaring in the H i thickness map. On the other hand, if this region is included in the averaging process, the averaged observed gas thickness would increase by ∼20% and if the H i velocity dispersion remains constant, one would expect a less massive dark matter halo from a simple minded calculation. Similarly, if we exclude this region the thickness in the north reduces to ∼3.4 kpc.

Another important fact about the LAB survey data is the pronounced N–S asymmetry in the gas thickness distribution. Certainly, an axisymmetric model of the Galaxy cannot reproduce such asymmetry in the data. In this context, recent work by Sanchez-Salcedo et al. (2008) has also shown that MOdified Newtonian Dynamics (MOND) provides a reasonably good fit to the azimuthaly averaged flaring in the H i gas using the same LAB survey data. Now for an axisymmetric baryon distribution, the MOND potential would also be axisymmetric, because the differential operator acting on the potential in the Poisson equation is rotationally invariant. Hence, it would probably be very hard to reproduce the observed N–S asymmetry in the H i thickness map. Whereas our family of lopsided dark matter halos naturally explain the observed N–S asymmetry.

The density distribution of the dark matter halo (namely the p = 2) in our model is similar to that obtained by Narayan et al. (2005; namely, a p = 2 halo as the best-fit model) based on the Wouterloot et al. (1990) H i flaring data. The axisymmetric version of our model is similar to that used by Narayan et al. (2005).

• 1.  
  Self-gravity of the atomic hydrogen gas. In the radial range from 10 to 16 kpc, we find that azimuthally averaged midplane volume density of the atomic hydrogen gas is comparable to that of the stars (see Figure 8) and beyond this radial range gas is in fact dominating over the stars. This suggests that the self-gravity of the H i gas is quite important in determining its thickness and thereby the nature of dark matter halo. By self-gravity of the gas, we mean that the gas is held by its own gravity. In other words, we include the gas density in the Poisson equation to derive its contribution to the total potential of the Galaxy. Without any self-gravity, the gas would move like a test particle under the potential of the Galaxy. So it would be interesting to check the change brought about by excluding the H i self-gravity on its vertical scale height. The difference is not negligible, it is seen to be about 10%–15% within the optical disk (R < 4R_d). This could be because these regions are dominated by the stellar disk and by the dark matter, respectively. For the region 4 < R/R_d < 6, the difference is substantial (∼20%–30%) suggesting that in this range, the gas gravity is very important in negotiating the hydrostatic equilibrium for the H i layer. Thus, neglecting it may lead to a serious overestimate of the H i scale height in general at all radii in the outer Galaxy and to explain the observed gas scale height one may need to invoke a heavier dark matter halo.
• 2.  
  Molecular hydrogen gas. In the present work, we are mostly concerned with the very outer region of the Galactic disk especially beyond 16 kpc. The particular reason for this is that the thickness curve of the atomic hydrogen gas in the northern hemisphere of the Galaxy deviates noticeably from the southern one beyond this region. The molecular hydrogen gas extends upto about 17 kpc in the Galaxy (Wouterloot et al. 1990) and beyond around this region there is little data. Since our calculation of gas thickness is local (in the sense that we do not consider the gas gravity in global sense), we have neglected the effect of molecular hydrogen on the H i thickness throughout our work. In fact, the absence of molecular hydrogen gas beyond about 16 kpc makes it more convenient in disentangling the effect of dark matter halo on the H i gas.
• 3.  
  Effect of Galactic constants. We have built the stellar disk model based on the IAU-recommended values for the galactic constants—R° = 8.5 kpc and Θ° = 220 km s⁻¹. It would definitely be interesting and also worthwhile to know how the results for the halo density profile vary with the assumed galactic constants. For example, Olling & Merrifield (2001) find the effect of varying these constants on the inferred axis ratio of the halo. Unfortunately, all the observational inputs for our model, such as the H i and H₂ surface densities, H i scale height, and the stellar velocity dispersion, are based on the IAU-recommended galactic constants and rescaling them for different values of the constants is beyond the scope of this paper.
• 4.  
  Rotation curve and gas thickness data. Different rotation curves produce different degrees of flaring because of the dependence of R on Θ through Equation (1) of Levine et al. (2006a). The rotation curves used in the present analysis differ by a small amount <10% from the flat rotation curve, Θ(R) = 220 km s⁻¹, used by Levine et al. (2006a). An error analysis shows that R changes approximately by ∼Θ/Θ_☉, where Θ_☉ is the circular speed at the solar position. Thus, dR/dΘ, the sensitivity of an error in the distance to an error in the circular speed produces an error in the distance of a parcel of gas from the center of about 10% at 30 kpc. Because all values of Galactic latitude are small at large R, this translates to an error in the thickness of the gas layer by no more that 10%, and then only at the largest distances.Although the differences in the rotation curves used by us and by Levine et al. (2006a) imply that the two analyses are not exactly commensurate, we are trying to explain a factor of 2 increase in the thickness of the gas layer of the Milky Way from one hemisphere to the other, an order of magnitude larger than the maximum 10% effect caused by differences in the rotation curves. Note that the factor of 2 change in the scale height occurs at all radii beyond about R = 17 kpc and the effect on ΔR will be smaller at smaller radii because ΔΘ is also smaller. Therefore, while there is a 10% effect at the largest distances, the most it will do is to have a 10% effect on the scaling of the model.It is worth mentioning at this point that we do not solve the vertical hydrostatic equilibrium considering epicylic orbit correction for the gas. For this, one needs to solve the Jeans equation and Poisson equation self-consistently, which is considerably more complex than what we do here, and in any event probably results in only a small correction. Our primary aim is to understand the nature of the observed asymmetry and build a first-order model to explain it.
• 5.  
  Effect of Galactic warp. The vertical structure of the Milky Way's disk is fairly complicated. The warp in the H i disk is very asymmetric; in the northern hemisphere, the disk midplane rises to a height of about 4 kpc, while in the south it goes down to about 1 kpc and then again comes back to the undisturbed midplane (Kerr 1957; Burton 1988; Levine et al. 2006a). Since the warp is a global feature in the Galaxy, it will have its own self-gravity to affect the vertical oscillations in the disk. The effective vertical oscillation frequency in the disk can be written as $\nu _{\rm eff} = \sqrt{\nu _{d}^2 + \nu _{\rm warp}^2}$, where ν_d is the vertical frequency of the unperturbed disk and ν_warp is the extra vertical frequency due to the enhanced self-gravity of the warp. Because the warp is global in nature, ν_warp is an integral quantity (see Equation (4b) in Saha & Jog 2006). However, an analytic form for ν_warp can be written in a local sense via the WKB analysis and it is ∼$\sqrt{2\pi G \Sigma _{d}|k|}$ (Binney & Tremaine 1987), where |k| is the wavenumber for the warp. Since we are solving the vertical hydrostatic equilibrium locally, the contribution from the global warp ($|k| \longrightarrow 0$) becomes insignificant. In order to account for the warp in the gas thickness, one has to formulate the vertical hydrostatic equilibrium as an integral problem which we plan for the future.
• 6.  
  Galactic spiral structure. The recent work by Levine et al. (2006b), based on the 21 cm LAB Galactic H i survey, reveals a multiarmed spiral structure of our Galaxy. Their study shows a good correlation between the positions of the spiral arms and the H i thickness. This is an expected behavior because the presence of the spiral arms would cause enhanced gravity and scattering of the H i clouds would not change the vertical velocity dispersion appreciably, leading to a decrement in the H i thickness. However, we do not expect the thickness of H i to change appreciably because the strength of the perturbed surface density does not vary strongly as a function of the Galactocentric radius along an arm or even from arm to arm (Levine et al. 2006b). In the present paper, under the zeroth-order approximation, we worked with the average data in the northern and southern halves of the Galaxy, respectively. This is certainly an incomplete modeling of the data and in future we expect to construct a self-consistent formulation of the problem to include the spiral structure of the Galaxy.
• 7.  
  Uniqueness of the lopsided halo model. In astronomy, it is quite a hard job to prove directly the uniqueness of a proposed model which fits a given set of observational data well. One obvious way of approach is to compare different theoretical models against the given data set. Given the various uncertainties in the observation and our incomplete understanding of the physics of our Galaxy, such a comparison may even lead to degeneracies and Bayesian analysis would probably be the appropriate way out to lifting such degeneracies or discard other obvious models based on physical grounds. Before we go into discussing other possible theoretical models, we would like to remind the reader some basic facts just for a recap. We are dealing with the outer region of the Galaxy and especially beyond R > 16 kpc which is ∼5 disk scale length and this is also the typical size of the stellar disk in a disk galaxy. So the influence of the stellar disk on the vertical distribution of the gas is almost negligible compared to that due to the dark matter halo.

Certainly, the internal disk instabilities alone (being normally weak, except the bar instability) are unlikely to produce the observed asymmetry in the gas thickness distribution.

One of the possible alternative models is an off-centered dark matter halo with respect to the Milky Way's disk. Such a model of an off-centered axisymmetric halo has been used previously by Levine & Sparke (1998) to generate lopsidedness in the disk. As pointed out by these authors that this method is effective when the disk lies within the core radius (almost constant density region of the halo) of the halo and as result is more efficient in dwarf galaxies rather than in luminous galaxies. In this context, one of our halo model, namely, p = 2 halo whose core radius is ∼9 kpc as compared to the disk scale length ∼3 kpc could have been a possible case for investigation. However, note that the N–S asymmetry in the gas distribution begins beyond 16 kpc which is roughly ∼2 × R_c of the halo and at this radius, the disk is no more in a nearly constant density region of the halo, making it harder to maintain the lopsidedness. In any case, this is an interesting possibility to be investigated in a future problem.

Other possible model is an asymmetric gas accretion onto the disk as proposed by Bournaud et al. (2005) to reproduce lopsidedness observed in the galactic disk. To produce strong lopsidedness as observed, one needs the accretion of the cold gas through the cosmological filaments to be highly asymmetric and in reality it is not clear if the gas accretion is asymmetric to such a degree. On the other hand, it makes sense to think in this direction because the H i distribution looks more disturbed in the north rather than in the south. At this point, it is worthwhile to mention that the gas actually contributes very little to the rotation curve (mostly dominated by the dark matter); so an axisymmetric p = 1 halo with an asymmetric gas accretion again may not be the good candidate for the present observation. However, an axisymmetric p = 2 halo with an asymmetric gas accretion could have been a possible candidate, but it is beyond the scope of the present paper to examine such a model in considerable detail.

On the other hand, our lopsided halo models are more natural to occur in cosmological scenario. Tidal interactions or large-scale tidal harassment or major mergers of neighboring dark matter halos are more likely to produce large-scale perturbations in the dark matter halo. And since the CDM halos are collisionless object, the survival of such global perturbations in the halo is less of a problem.

We have analyzed both axisymmetric and nonaxisymmetric configurations of dark matter halos with a consistent picture of the ISM to explain the observed nature of the N–S asymmetry in the thickness distribution of the H i gas. Below we summarize the main results.

With a model of the ISM that has reasonable values of the gas velocity dispersion, an isothermal dark matter halo producing a flat rotation curve with V_c = 220 km s⁻¹ cannot produce the observed flaring in the H i gas thickness.

We show that the nature of the systematic N–S asymmetry in the H i thickness map is gravitational in nature. An axisymmetric dark matter halo with different values of the H i velocity dispersion in the two halves of the Galaxy cannot reproduce this asymmetry. The observed asymmetry in the thickness map of neutral hydrogen gas is apparently not the result of purely gas dynamical effects.

We show that a purely p = 1.5 or p = 2 lopsided dark matter halo also cannot explain the observed N–S asymmetry. For a purely p = 2 lopsided halo, the H i velocity dispersion has to be unreasonably large to come close to the observation and even then, the fit is not very good.

Finally, we come up with a configuration of the dark matter halo in which some amount of second harmonic (m = 2) is superposed out of phase onto a purely p = 2 lopsided halo. For the best-fit models (A and C), the values of lopsidedness and elliptical perturbation are ¹_h = 0.2 and ²_h = 0.18, respectively. We call such a halo an elliptically perturbed lopsided dark matter halo which can explain the observed N–S asymmetry. Basically, the emerging picture of the dark matter halo of the Milky Way is dominantly lopsided in nature. In such a halo, the density falls off faster than the p = 1 isothermal dark matter halo. The emerging model of the asymmetric dark matter halo is supported by the halos formed in the recent cosmological N-body simulation.

We thank the referee for the critical and constructive comments which have improved the contents of the paper substantially. K.S. thanks the Raman Research Institute for supporting a visitorship during which an initial draft of the paper started. L.B. and E.L. like to acknowledge support of NSF grant AST-0540567 to the University of California.