Vertical dynamics of disk galaxies in MOND

TL;DR

This study compares the vertical dynamics of disk galaxies under MOND and Newtonian gravity with dark matter, finding that differences in vertical velocity dispersion are often too subtle for current observations to distinguish.

Contribution

It provides a detailed analysis of vertical kinematics in disk galaxies under MOND and Newtonian models, highlighting the challenges in observationally discriminating between the two.

Findings

Deep-MOND disks can have much larger vertical velocity dispersions than Newtonian models.

High-surface density disks show reduced differences, complicating observational discrimination.

Current observations may not be sufficient to distinguish between MOND and dark matter models.

Abstract

We investigate the possibility of discriminating between Modified Newtonian Dynamics (MOND) and Newtonian gravity with dark matter, by studying the vertical dynamics of disk galaxies. We consider models with the same circular velocity in the equatorial plane (purely baryonic disks in MOND and the same disks in Newtonian gravity embedded in spherical dark matter haloes), and we construct their intrinsic and projected kinematical fields by solving the Jeans equations under the assumption of a two-integral distribution function. We found that the vertical velocity dispersion of deep-MOND disks can be much larger than in the equivalent spherical Newtonian models. However, in the more realistic case of high-surface density disks this effect is significantly reduced, casting doubts on the possibility of discriminating between MOND and Newtonian gravity with dark matter by using current observations.

Figures

Equations

$${\nabla }^2{\varphi }^{\mathrm{N}}=4\pi G\rho ,$$

$$\nabla \cdot \left[\mu \left(\frac{\Vert \nabla \varphi \Vert }{a_0}\right)\nabla \varphi \right]=4\pi G\rho .$$

$$\nabla \cdot \left(\Vert \nabla \varphi \Vert \nabla \varphi \right)=4\pi G{a}_0\rho .$$

$$\mu (g/a_0)\mathbf{g}={\mathbf{g}}^{\mathrm{N}}+\mathbf{S},$$

$$v_{\mathrm{c}}^2(R)=R\frac{\partial {\varphi }_{\mathrm{T}}(R,0)}{\partial R},$$

$${\rho }_{\mathrm{DM}}(r)=\frac{1}{4\pi G{r}^2}\frac{d}{dr}\left[r{v}_{\mathrm{c}}^2(r)-r{v}_{\mathrm{cN}}^2(r)\right],$$

$$\frac{\partial \rho {\sigma }^2}{\partial z}=-\rho \frac{\partial {\varphi }_{\mathrm{T}}}{\partial z},$$

$$\frac{\partial \rho {\sigma }^2}{\partial R}+\frac{\rho ({\sigma }^2-\overline{{\mathrm{v}}_{\phi }^2})}{R}=-\rho \frac{\partial {\varphi }_{\mathrm{T}}}{\partial R},$$

$$\rho {\sigma }^2={\int }_z^{\infty }\rho \frac{\partial {\varphi }_{\mathrm{T}}}{\partial z^{\prime }}d{z}^{\prime }.$$

$$\rho {\sigma }^2={\int }_z^{\infty }\rho \frac{\partial {\varphi }^{\mathrm{N}}}{\partial z^{\prime }}d{z}^{\prime }+{\int }_z^{\infty }\rho \frac{d{\varphi }_{\mathrm{DM}}}{dr}\frac{z^{\prime }}{r}d{z}^{\prime },$$

$$u_{\phi }^2=k^2(\overline{{\mathrm{v}}_{\phi }^2}-{\sigma }^2),$$

$${\sigma }_{\phi }^2={\sigma }^2+(1-k^2)(\overline{{\mathrm{v}}_{\phi }^2}-{\sigma }^2),$$

$$\Sigma v_{\mathrm{p}}\equiv {\int }_{-\infty }^{\infty }\rho u_{\mathrm{n}}dl,$$

$$\Sigma V_{\mathrm{p}}^2\equiv {\int }_{-\infty }^{\infty }\rho u_{\mathrm{n}}^{2}dl,$$

$$\Sigma {\sigma }_{\mathrm{p}}^2\equiv {\int }_{-\infty }^{\infty }\rho {\sigma }_{\mathrm{n}}^{2}dl,$$

$$\Sigma {\sigma }_{\mathrm{los}}^2\equiv {\int }_{-\infty }^{\infty }\rho \overline{{\left(⟨\mathbf{v},\mathbf{n}⟩-v_{\mathrm{p}}\right)}^2}dl=\Sigma \times \left({\sigma }_{\mathrm{p}}^2+V_{\mathrm{p}}^2-v_{\mathrm{p}}^2\right).$$

$$\Sigma {\sigma }_{\mathrm{p}}^2=2{\int }_0^{\infty }\rho {\sigma }^{2}dz=2{\int }_0^{\infty }\rho \frac{\partial {\varphi }_{\mathrm{T}}}{\partial z}zdz,$$

$$u_{\mathrm{n}}=-u_{\phi }\sin \phi ,$$

$${\sigma }_{\mathrm{n}}^2={\sigma }^2+(1-k^2)(\overline{{\mathrm{v}}_{\phi }^2}-{\sigma }^2){\sin }^2\phi ,$$

$$\rho (R,z)=\frac{M}{4\pi b^3}\frac{s{\tilde{R}}^2+(s+3\zeta )(s+\zeta {)}^2}{[{\tilde{R}}^2+(s+\zeta {)}^2{]}^{5/2}{\zeta }^3},$$

$$\varphi (R,z)=-\frac{GM}{b}\frac{1}{\sqrt{{\tilde{R}}^2+(s+\zeta {)}^2}}.$$

$$v_{\mathrm{cN}}^2(R)=\frac{GM}{b}\frac{{\tilde{R}}^2}{[{\tilde{R}}^2+(s+1{)}^2{]}^{3/2}}.$$

$$v_{\mathrm{c}}^2(R)=\frac{\sqrt{GM{a}_0}{\tilde{R}}^{3/2}}{[{\tilde{R}}^2+(s+1{)}^2{]}^{3/4}},$$

$${\rho }_{\mathrm{DM}}(r)\simeq \frac{M}{4\pi b^3}\sqrt{\frac{a_0{b}^2}{GM}}\frac{2{\tilde{r}}^2+5(s+1{)}^2}{2\sqrt{\tilde{r}}[{\tilde{r}}^2+(s+1{)}^2{]}^{7/4}},$$

$$g_{z,\mathrm{disk}}=\sqrt{\frac{a_0}{||{\mathbf{g}}^{\mathrm{N}}||}}{g}_z^{\mathrm{N}},$$

$$g_{z,\mathrm{DM}}^{\mathrm{N}}=-\frac{\partial {\varphi }_{\mathrm{DM}}(r)}{\partial z}=-\frac{z{v}_{\mathrm{c}}^2(r)}{r^2}.$$

$$\frac{g_{z,\mathrm{DM}}^{\mathrm{N}}}{g_{z,\mathrm{disk}}}\sim \frac{1}{1+s}+\left[\frac{s}{2(1+s{)}^2}+\frac{s(5+s)}{4(1+s){\tilde{R}}^2}\right]{\tilde{z}}^2,$$

$$\rho (R,z)=\frac{M}{4\pi R_0^2{z}_0}\exp \left(-\frac{R}{R_0}\right){\mathrm{sech}}^2\left(\frac{z}{z_0}\right),$$

$$\mu (t)=\frac{t}{1+t}.$$

Full text

††pubyear: 2007

1 Introduction

The flatness of the rotation curves of disk galaxies in their external regions has been the primary focus in the lasting on dark matter (DM) and Modified Newtonian Dynamics (MOND). As a consequence of the formalisation of MOND in fundamental physics (Bekenstein 2004; Zlosnik, Ferreira & Starkman 2007), combined with some difficulties of the cold dark matter scenario on galaxy scales (e.g. Binney 2004, and references therein), the debate extended to other astrophysical contexts. For example, MOND models have been recently tested against the Cosmic Background Radiation (Skordis et al. 2006, Skordis 2006), gravitational lensing (Chen & Zhao 2006; Zhao et al. 2006; Angus et al. 2007; Takahashi & Chiba 2007), star cluster and galaxy dynamics (Zhao & Tian 2006; Sanchez-Salcedo, Reyes-Iturbide & Hernandez 2006; Haghi, Rahvar & Hasani-Zonooz 2006; Nipoti, Londrillo & Ciotti 2007; Scarpa et al. 2007; Tiret & Combes 2007), and the solar system (Bekenstein & Magueijo 2006, Sanders 2006, Sereno & Jetzer 2006). Due to the surprising ability of MOND to reproduce the kinematics of different systems (e.g., Milgrom 2002; Sanders & McGaugh 2002; Bekenstein 2006), it is of obvious interest to look for other tests to discriminate between DM and MOND.

In this paper we explore the possibility of differentiating DM and MOND by studying the vertical dynamics of disk galaxies. While it is often taken for granted that—as far as disk kinematics is concerned—the DM and MOND interpretations are nearly degenerate, in fact this is not true. An important issue is the choice of the Newtonian system used for comparison with the MOND results. In the present context, the proper comparison is between a baryonic disk in MOND, and the same disk, in Newtonian gravity, immersed in a DM halo (which for simplicity we assume spherically symmetric), such that the circular velocity in the equatorial plane is the same as in the MOND model111Other cases are known in which MOND and Newtonian systems are equivalent with respect to some dynamical properties, but different with respect to others (e.g. Ciotti & Binney 2004).. In particular, we compare the kinematical fields of equivalent models under the assumption of two-integral distribution function. A work in some respect complementary to ours, though with important differences, was done by Milgrom (2001), who studied, for a given baryonic disk, the shape of the DM halo of the Newtonian model with exactly the same dynamics (not only the same rotation curves) as the disk in MOND gravity. Recently, a test of MOND based on the comparison of disk circular velocity and vertical velocity dispersion has been proposed also by Stubbs & Garg (2005)222The presented test is affected by the unjustified assumption that the vertical velocity dispersion depends only on the local gravitational field..

The paper is organised as follows. In Section 2 we describe the general method adopted, and we perform a preliminary investigation of the problem. In Section 3 we study in detail the projected kinematical fields of Miyamoto-Nagai and thick exponential disks in deep MOND regime, and we compare them with their equivalent Newtonian models with DM, for different values of flattening/thickness. We also address the problem of MOND dynamics in Milky-Way like galaxy models. Our results are summarised in Section 4.

2 The method

In the present work we consider MOND in Bekenstein & Milgrom’s (1984) formulation, in which the Poisson equation

[TABLE]

where $\phi^{\rm N}$ is the Newtonian gravitational potential generated by the density distribution $\rho$, is substituted by the non-relativistic field equation

[TABLE]

In equation above $\|...\|$ is the standard Euclidean norm, $\phi$ is the MOND gravitational potential produced by $\rho$, and for finite mass systems $\nabla\phi\to 0$ for $\|{\bf x}\|\to\infty$. The interpolating function $\mu(t)$ is not constrained by theory except that it must run smoothly from $\mu(t)\sim t$ at $t\ll 1$ to $\mu(t)\sim 1$ at $t\gg 1$, with a dividing acceleration scale ${a_{0}}\simeq 1.2\times 10^{-10}{\rm m}\,{\rm s}^{-2}$. In the so-called ‘deep MOND regime’ (hereafter dMOND), describing low-acceleration systems ($\|\nabla\phi\|\ll{a_{0}}$), $\mu(t)=t$ and equation (2) reduces to

[TABLE]

As well known, equation (2) can be combined with equation (1) so that the MOND ${\bf g}=-\nabla\phi$ and Newtonian ${\bf g}^{\rm N}=-\nabla\phi^{\rm N}$ gravitational fields are linked as

[TABLE]

where ${\bf S}$ is a solenoidal field dependent on the specific $\rho$ considered. It can be proved that the equation above reduces to Milgrom’s (1983) empirical relation (i.e., ${\bf S}=0$) in case of one-dimensional symmetries and in the special case of the razor-thin Kuzmin (1956) disk density distribution (Brada & Milgrom 1995), but in general one cannot impose ${\bf S}=0$. The presented results are derived by solving the field equation (2) with the numerical potential solver presented in Ciotti, Londrillo & Nipoti (2006).

2.1 Two-integral dynamics

In the following we compare the dynamics of two disk galaxy models with the same circular velocity curve in the equatorial plane

[TABLE]

where $\phi_{\rm T}(R,z)$ is the total potential of the model, and $(R,z,\varphi)$ are the standard cylindrical coordinates. One model is built in MOND for a baryonic disk galaxy described by the density distribution $\rho(R,z)$. The associated Newtonian model, which we call equivalent, consists of the same baryonic disk plus a spherical DM halo, whose radial density profile is determined by the condition of matching the rotation curve obtained in the MOND model. In practise, the density distribution of the DM halo of the equivalent model is given by

[TABLE]

where $v_{\rm c}$ and $v_{\rm cN}(R)$ are the disk circular velocity in MOND and Newtonian gravity, respectively. Note that the positivity of $\rho_{\rm DM}$ is not guaranteed for a generic disk density, so we always check the positivity of the equivalent DM halo density distribution.

For a two-integral distribution function, i.e., $f=f(E,J_{z})$, the Jeans equations for the disk are

[TABLE]

[TABLE]

where $\sigma^{2}\equiv\sigma_{R}^{2}=\sigma_{z}^{2}$, $\overline{\rm v_{\varphi}^{2}}=u_{\varphi}^{2}+\sigma_{\varphi}^{2}$, $u_{\varphi}=\overline{\rm v}_{\varphi}$, and a bar over a symbol indicates its phase-space mean over the velocity space (e.g., see Binney & Tremaine 1987, hereafter BT). Note that in a MOND model $\phi_{\rm T}$ is the solution of equation (2) for the disk density distribution, while in the equivalent Newtonian model $\phi_{\rm T}=\phi^{\rm N}(R,z)+\phi_{\rm DM}(r)$, with $\phi_{\rm DM}$ determined from the density profile in equation (6).

The velocity dispersion $\sigma^{2}(R,z)$ is obtained from integration of equation (7) with boundary condition $\rho\sigma^{2}=0$ for $z\to\infty$, i.e.

[TABLE]

In the equivalent model

[TABLE]

where $r=\sqrt{R^{2}+z^{\prime 2}}$ and $d\phi_{\rm DM}/dr=(v_{\rm c}^{2}-v_{\rm cN}^{2})/r$.

To split $\overline{\rm v_{\varphi}^{2}}$ into streaming motion $u_{\varphi}$ (that for simplicity we assume nowhere negative), and azimuthal dispersion $\sigma_{\varphi}^{2}$, we adopt the Satoh (1980) $k$–decomposition

[TABLE]

and

[TABLE]

with $0\leq k\leq 1$; this procedure can be applied only when $\overline{{\rm v}_{\varphi}^{2}}-\sigma^{2}\geq 0$ everywhere. For $k=0$ no ordered motions are present, and the velocity dispersion tensor is maximally tangentially anisotropic, while in the isotropic rotator ($k=1$) the galaxy flattening is due to azimuthal streaming velocity. In principle, by allowing for $k=k(R,z)$, even more rotationally supported models can be constructed, up to the maximum rotation case considered in Ciotti & Pellegrini (1996), where $\sigma_{\varphi}^{2}=0$ everywhere.

The explicit projection formulae for the kinematical fields of axisymmetric two-integral systems viewed along a generic direction ${\bf n}$ of the line of sight (los) can be found elsewhere (e.g., see Lanzoni & Ciotti 2003, Ciotti & Bertin 2005, Riciputi et al. 2005): here we just recall that at any given place in the galaxy the *los *component of the streaming velocity field is $u_{\rm n}\equiv\overline{\langle{\bf v},{\bf n}\rangle}=u_{i}n_{i}$, where ${\langle,\rangle}$ is the standard inner product, the repeated index convention is applied, and Cartesian coordinates are assumed. The analogous quantity associated with the velocity dispersion tensor is $\sigma_{\rm n}^{2}\equiv\overline{\langle{\bf v}-{\bf u}\rangle^{2}}=\sigma_{ij}n_{i}n_{j}$, and the corresponding (mass–weighted) kinematical projected fields are obtained by integration along the *los *as

[TABLE]

[TABLE]

and

[TABLE]

where $\Sigma=\int_{-\infty}^{\infty}\rho\,dl$ is the disk density projected along ${\bf n}$. In general $\sigma_{\rm p}$ is not the velocity dispersion measured in observations: in the presence of a non–zero projected velocity field $v_{\rm p}$, this quantity is

[TABLE]

For simplicity, in this paper we focus on the limit cases of disk observed face-on and edge-on. The face-on projected velocity dispersion is given by

[TABLE]

where $\Sigma(R)=2\int_{0}^{\infty}\rho(R,z)dz$ is the face-on disk density. In the edge-on case, assuming the *los *directed along the $x$ axis, the *los *direction in the natural coordinate system is given by333The *los *vector points towards the observer, and so positive velocities correspond to a blue–shift. ${\bf n}=(1,0,0)$ and the projection plane is $(y,z)$, so that

[TABLE]

and

[TABLE]

where $\sin\varphi=y/R=y/\sqrt{y^{2}+x^{2}}$. Note that, after projection, the $y$ coordinate can been identified with $R$.

2.2 A preliminary analysis

As an introductory exercise we compare, as a function of disk flattening, the vertical force field near the equatorial plane of a Miyamoto & Nagai (1975, hereafter MN) disk in dMOND regime and its equivalent Newtonian model. In order to carry out the calculations analytically we assume ${\bf S}=0$, while in Section 3 the exact MOND gravitational field with the appropriate ${\bf S}$ field is used. The MN density distribution with scale length $b$ and flattening parameter $s\equiv a/b$ is

[TABLE]

where $M$ is the total disk mass, $\zeta\equiv\sqrt{1+\tilde{z}^{2}}$, $\tilde{z}\equiv z/b$, and $\tilde{R}\equiv R/b$. The corresponding Newtonian potential is (BT)

[TABLE]

For $a=0$ the MN density distribution is a Plummer (1911) sphere with scale radius $b$, while for $b=0$ one obtains the Kuzmin (1956) disk. The Newtonian disk circular velocity is given by equations (5) and (21) as

[TABLE]

In the ${\bf S}=0$ approximation the dMOND circular velocity is obtained from equation (4) as

[TABLE]

and according to equation (6), $\rho_{\rm DM}\propto A(\tilde{r})-\sqrt{GM/{a_{0}}b^{2}}B(\tilde{r})$, where $A$ and $B$ are two dimensionless functions, and $\tilde{r}\equiv r/b$. In the dMOND limit $\sqrt{GM/{a_{0}}b^{2}}=0$ and so

[TABLE]

which is positive everywhere and has the expected $r^{-2}$ behaviour for $r\to\infty$. Note that equation (24) implies that in the equivalent Newtonian model the MN disk is just a tracer of the DM gravitational field.

The obtained formulae are useful because they allow to evaluate the vertical gravitational field near the disk plane in the two models. In the ${\bf S}=0$ approximation the dMOND vertical field is

[TABLE]

where $g_{z}^{\rm N}$ is the vertical Newtonian gravitational field of the MN disk, and the vertical gravitational field of the equivalent Newtonian DM halo is

[TABLE]

Simple algebra then shows that for $R\to\infty$ and $z\to 0$,

[TABLE]

i.e., the vertical field is stronger in dMOND than in DM models, and the discrepancy between the two cases is larger for flatter models (i.e. for larger values of $s$). Thus, we expect $\sigma^{2}$ in the disk to be different in the two models. In turn, equations (8) and (11) imply that also the streaming velocity field in the disk equatorial plane is different in the two models, even though the circular velocity is the same by construction.

3 Results

We now obtain the kinematical fields (for simplicity in the fully isotropic case) of a few disk galaxy models by computing the exact MOND and Newtonian potentials with the numerical code presented in Ciotti et al. (2006). The potential solver is based on a spherical grid of coordinates ($r$, $\vartheta$, $\varphi$), with $N_{r}\times N_{\vartheta}\times N_{\varphi}$ points, uniformly spaced in ($\arctan{r}$, $\vartheta$, $\varphi$). To compute the intrinsic and projected velocity dispersions, we interpolate the potential from the spherical grid to a cylindrical grid with $N_{R}\times N_{z}$ points, uniform in ($\arctan{R}$, $\arctan{z}$), on which we evaluate the integrals (9) and (13-15). The numerical integration routines are verified by deriving the intrinsic and projected velocity dispersions of MN disks in Newtonian gravity, and comparing the results with the analytical expressions reported in Ciotti & Pellegrini (1996). For instance, with $N_{r}=N_{\vartheta}=N_{R}=N_{z}=128$ the analytical results are reproduced with relative errors $\sim 0.1$ %.

3.1 Miyamoto-Nagai disks

We consider two MN galaxy models with flattening parameter $s=1$ (i.e., an almost spherical system), and $s=5$. Their exact dMOND fields are computed by solving numerically equation (3). The density distributions of the spherical DM haloes of the equivalent Newtonian models are then obtained from equation (6) with $v_{\rm cN}=0$, as appropriate for the dMOND case, in which the disk is just a tracer of the halo potential. As expected, the resulting halo profiles (Fig. 1, left panel; solid and dashed lines) decrease as $r^{-2}$ in the outer regions, matching the analytical formula (24), obtained in the ${\bf S}=0$ approximation (dotted lines). However, for r{\mathrel{\mathchoice{\lower 2.9pt\vbox{\halign{\m@th\displaystyle\hfil#\hfil\cr<\crcr\sim\crcr}}}{\lower 2.9pt\vbox{\halign{\m@th\textstyle\hfil#\hfil\cr<\crcr\sim\crcr}}}{\lower 2.9pt\vbox{\halign{\m@th\scriptstyle\hfil#\hfil\cr<\crcr\sim\crcr}}}{\lower 2.9pt\vbox{\halign{\m@th\scriptscriptstyle\hfil#\hfil\cr<\crcr\sim\crcr}}}}}b the exact halo density profiles are considerably flatter than those predicted by equation (24), having inner logarithmic slope d\log\rho_{\rm DM}/d\log r{\mathrel{\mathchoice{\lower 2.9pt\vbox{\halign{\m@th\displaystyle\hfil#\hfil\cr>\crcr\sim\crcr}}}{\lower 2.9pt\vbox{\halign{\m@th\textstyle\hfil#\hfil\cr>\crcr\sim\crcr}}}{\lower 2.9pt\vbox{\halign{\m@th\scriptstyle\hfil#\hfil\cr>\crcr\sim\crcr}}}{\lower 2.9pt\vbox{\halign{\m@th\scriptscriptstyle\hfil#\hfil\cr>\crcr\sim\crcr}}}}}-1/2. The difference (which is larger in the flatter disk) indicates that the contribution of the field ${\bf S}$ is non-negligible in the central regions, as also found by Ciotti et al. (2006) for triaxial and axisymmetric Hernquist (1990) models in dMOND regime. However, in the outer regions the agreement is perfect.

In Fig. 2 (left panel), we plot the ratio $g^{\rm N}_{z,{\rm DM}}/g_{z,{\rm disk}}$ as a function of $z$ at the representative radius $R=b$. Again, the results nicely confirm the asymptotic estimate of equation (27) when $z=0$: the effect is bigger in the flatter disk, and just above the disk plane $g_{z,{\rm disk}}$ is significantly stronger than $g^{\rm N}_{z,{\rm DM}}$. We also note that for $z\gg b$ the vertical Newtonian field exceeds the dMOND one, though slightly. As a consequence, also the intrinsic and projected velocity dispersions are different in the equivalent Newtonian and dMOND MN disks, while the circular velocities in the equatorial plane are the same by construction. The midplane vertical velocity dispersion $\sigma(R,0)$ and the projected face-on velocity dispersion $\sigma_{\rm los}$ are larger in the dMOND systems than in the equivalent Newtonian systems at all radii $R<5b$, with the larger discrepancies (up to a factor of 2) in the flatter disk (Fig. 3, left column). The situation is more complicate in the edge-on projection: the projected streaming velocity $v_{\rm p}$ in the Newtonian models are larger than in the dMOND cases, while the edge-on projected velocity dispersion $\sigma_{\rm los}$ in the galactic plane is higher in dMOND than in Newtonian models.

3.2 Thick exponential disks

We now investigate the more realistic thick exponential disk

[TABLE]

where $M$ is the total disk mass, and $R_{0}$ and $z_{0}$ are characteristic scale-lengths. As for the MN case, we consider here two dMOND models (and so in the equivalent Newtonian model the disk does not contribute to the potential): a thicker disk with $R_{0}=3z_{0}$, and a thinner disk with $R_{0}=15z_{0}$.

As in MN disks, the spherical DM haloes density profiles of the equivalent Newtonian models are rather flat in the central regions and $\propto r^{-2}$ at large radii (Fig. 1, right panel). Interestingly, these DM haloes “predicted” by MOND for low-surface brightness disks are not characterised by the steep central cusps expected in the context of cold dark matter. The dMOND vertical force near the disk is significantly stronger than in the equivalent Newtonian models (Fig. 2, right panel): at $R=R_{0}$ in the thinner model the dMOND vertical force ratio is $\sim 10$ just above the plane. This difference reflects in the velocity fields (Fig. 3, right panels), which behave qualitatively as the MN models, but the discrepancies are larger: for example, the ratio of intrinsic velocity dispersion is as high as $\sim 4$ in the thinner model.

3.3 Milky-Way like galaxies

The previous analysis restricted to models in the dMOND regime, in order to estimate the largest possible differences in the kinematical fields of equivalent disk galaxies. We now focus on a Milky-Way like galaxy model, with the purpose of quantifying the MOND kinematical effects in realistic high surface-brightness galaxies. We will not enter the challenging problem of finding the best-fitting MOND model of the Milky Way (which is beyond the aim of the present paper; see Famaey & Binney 2005); instead, we apply our method to one of the currently available mass models of the Milky Way that have been derived in the context of DM. In particular, we consider Model 1 of Dehnen & Binney (1998; hereafter model DB), which consists of a truncated oblate spheroidal power-law bulge, a stellar disk made of a thin and a thick component, a gaseous exponential disk with a central hole, and finally an oblate spheroidal DM distribution. Model DB represents a good Newtonian dynamical model of the solar neighbourhood and the Galaxy, and can be classified as a disk dominated model, because in Newtonian gravity the disk provides the main contribution to the rotation curve in the central regions.

As we solve now equation (2), we need to specify the interpolating function $\mu$. Following Famaey & Binney (2005), we fix ${a_{0}}\simeq 1.2\times 10^{-10}{\rm m}\,{\rm s}^{-2}$ and

[TABLE]

This function gives a better fit to the terminal velocity curve of the Milky Way than the commonly adopted function $\mu(t)=t/\sqrt{1+t^{2}}$ and fits extremely well the rotation curve of NGC3198, also being consistent with Bekenstein’s (2004) TeVeS theory (Zhao & Famaey 2006; see also Famaey et al. 2007). We consider two models with the same baryonic distribution as model DB: model DBMOND (a purely baryonic MOND model), and model DBEQ (the equivalent Newtonian model with spherical DM halo). Clearly, the DM halo density profile in model DBEQ is not identical to that of model DB: Fig. 4 (top panel) shows that in the equatorial plane the two distributions are very similar in the radial range $R\sim 5-15\,{\rm kpc}$, while the halo profile of DBEQ is steeper than that of DB at small radii. As a consequence, the rotation curve of models DBMOND and DBEQ (which is the same by construction) differs from the rotation curve of DB, being systematically higher at small radii and lower at large radii (Fig. 4, bottom). However, the differences in the radial range $0-15\,{\rm kpc}$ are within $\sim 5\%$, and $v_{\rm c}(8\,{\rm kpc})=226\,{\rm km\,s}^{-1}$ in model DBMOND (and DBEQ), consistent with the observational constraint at the solar neighbourhood $220\pm 15\,{\rm km\,s}^{-1}$ (solid symbol and vertical bar in the diagram; e.g., BT).

In contrast with the dMOND cases discussed in Sections 3.1 and 3.2, now the contribution of the baryonic component to the gravitational field is important also in the Newtonian model, and this results in smaller differences between the MOND and Newtonian equivalent models. For instance, considering the vertical force as a function of $z$ at $R=8\,{\rm kpc}$ (Fig. 5), we find that the MOND force is stronger than the Newtonian one, with a difference of $\sim 30\%$ at $z\sim$1 kpc. The solid dot marks an observational estimate (and associated uncertainty) of the vertical force 1.1 kpc above the plane at the solar radius (Kuijken & Gilmore 1989, 1991). Remarkably, the difference between the vertical forces at $z=1.1\,{\rm kpc}$ predicted by the MOND and Newtonian models is larger than the observational error. So, taken at face value, this result would favour DM rather than MOND models, though we stress again that we did not attempt to build the best-fit MOND model of the Milky Way, and we cannot exclude the possibility of finding a MOND model of the Milky Way satisfying this and other observational constraints (see also Famaey & Binney 2005). In any case, this result shows that in principle it could be possible to try to falsify MOND by using the observational constraints on the vertical force above the plane.

Following the treatment of Sections 3.1 and 3.2 we computed the kinematical fields of models DBMOND and DBEQ, again under the same assumptions of two-integral distribution function444Because of this assumption our results on the velocity fields cannot be intended, strictly speaking, to represent the specific case of the Milky Way, whose distribution function is not two-integral, being $\sigma_{R}^{2}>\sigma_{z}^{2}$ in the solar neighbourhood. and full isotropy. For both models we computed the intrinsic and projected velocity fields of the thin-disk (thin curves) and thick-disk (thick curves) stars (Fig. 6, left panels), finding that the discrepancy between the two models is up to 20-25% for intrinsic and face-on velocity dispersion, while just few per cent in case of edge-on projection (Fig. 6, right panels). These results show that the differences in the vertical kinematical fields in a realistic high-surface density galaxy are quite small, and it is not obvious that one can discriminate between MOND and Newtonian gravity (plus DM) by using the information currently derived by the observations.

4 Discussion and conclusions

In this paper we quantified the differences in the vertical force and kinematical fields between MOND disk galaxy models and equivalent (i.e., having identical circular velocity in the equatorial plane) Newtonian models with spherical DM haloes. We showed that, in principle, MOND and Newtonian gravity with DM can be differentiated with accurate measurements of the vertical force near the disk and of the projected kinematical fields. In particular, MOND models have stronger vertical force near the plane and higher vertical velocity dispersion than equivalent Newtonian models with DM. Our results are not inconsistent with the finding of Read & Moore (2005) that satellite orbits in DM and MOND behave very similarly, because far from the disk plane the force fields predicted by the two theories are very similar.

From the observational point of view, the robust discrimination between MOND and Newtonian gravity with DM is still challenging due to several factors: (i) the strongest effects are expected in low-acceleration (low-surface-brightness) systems, in which measurements of the velocity fields are difficult; (ii) there are often several stellar populations with different scale-lengths and heights, and the mass-to-light of the stellar component is tunable; (iii) most disks have gaseous components whose density has non-smooth features which could make differences in any local volume of kpc scale; (iv) the models should be compared to data with noise so our comparisons of MOND and DM models of exactly same circular velocity curves remains a theoretical exercise; (v) DM haloes are not necessarily spherically symmetric.

With regard to point (v), we briefly discuss how our results would change by relaxing the assumption of an equivalent spherical DM halo. Following Milgrom (2001), we start by constructing, for a given disk distribution $\rho$, the “fully equivalent” DM density distribution as $\rho_{\rm DM}=\nabla^{2}\phi/4\pi G-\rho$, where $\phi$ is the MOND potential; as well known, $\rho_{\rm DM}$ can be negative in some region of space (Milgrom 1986; Ciotti et al. 2006). By construction, the fully equivalent Newtonian model has gravitational field identical to the MOND disk field in all the space, not just the same midplane circular velocity (as the equivalent model with spherical halo). Interestingly, we found that the density distributions $\rho_{\rm DM}$ of the fully equivalent DM haloes of the dMOND models presented in Sections 3.1 and 3.2 are nowhere negative (in these cases $\rho_{\rm DM}=\nabla^{2}\phi/4\pi G$, because the baryonic disk is just a tracer in dMOND-equivalent models). Figure 7 shows that $\rho_{\rm DM}$ presents a significant disk structure (more pronounced in the case of more flattened baryonic disks). This result clearly indicates that it is not possible to obtain a Newtonian model with spheroidal DM distribution strictly equivalent to a MOND disk model (see also Milgrom 2001). However, we found that the discrepancy in the vertical dynamics in the galactic plane between MOND and DM models can be significantly reduced considering spheroidal DM haloes: for instance, a dMOND MN disk with $s=5$ and an oblate spheroidal halo with axis ratio $\sim 0.25-0.3$ (and the same asymptotic circular velocity as the baryonic dMOND disk) would produce very similar vertical velocity dispersion and rotational velocity in the galactic plane (excluding the very central regions).

For all these reasons, it appears difficult to differentiate between MOND and DM using the currently available observational measures of stellar kinematics of external disk galaxies (e.g., Bottema 1993; Kregel & van der Kruit 2005, and references therein). The most promising application of our method seems to be the case of the Milky Way: in particular, for the specific baryonic mass distribution considered in this paper the difference between the MOND and DM vertical forces above the plane is larger than the uncertainty on the observational estimates. The estimates of both the baryonic mass distribution and the vertical force above the plane are expected to be greatly improved in the near future thanks to data from the GAIA mission.

Acknowledgments

We thank Renzo Sancisi for helpful discussions, and the Referee, Robert Sanders, for his useful suggestions. L.C. and P.L. were partially supported by a MIUR grant Cofin 2004.