ON THE STABILITY OF A GALACTIC DISK IN MODIFIED GRAVITY

Disk galaxies encompass non-axisymmetric features such as spiral arms and bars, which have significant influences on the galactic evolution (Kormendy & Kennicutt 2004; Sellwood 2014). These spiral arms trigger star formation in the outer part of the disk galaxy, drive secular changes in the stellar orbits and change the mass distribution within the disk. Understanding the origin and evolution of these spiral structures is one of the hardest problems in astrophysics. However, it is widely believed that spiral patterns are gravitationally driven density waves (e.g., those produced by a spontaneous disturbance or, in some cases, a companion system) and propagating in the surface of the galaxies. Therefore the density wave theory is an indispensable tool in this subject. Assuming that density waves are small perturbations on the density of the galactic disk, one can use linear analysis in order to find the linear oscillatory modes of the disk. However, analyzing the behavior of a density wave in a galactic disk and calculating the relevant oscillatory modes is a difficult task. The main reason for this difficulty is that gravity plays the dominant role in the evolution of the density waves. On the other hand, gravity is a long range force. Correspondingly, perturbations in all parts of the system are coupled. This coupling makes the analysis very difficult and consequently the oscillatory modes can be found analytically only for Kalnajs disks (Kalnajs 1972). Also numerical mode calculations have been performed for some simple models including the isochrone disk, the Kuzmin disk, and power-law disks (Zang 1976; Kalnajs 1978; Vauterin & Dejonghe 1996; Pichon & Cannon 1997; Evans & Read 1998; Jalali & Hunter 2005).

In the early 1960s it was understood that for tightly wound spiral density waves (waves whose wavelength is much less than the radius), the long-range coupling of gravity is negligible, and the response of the background disk to these tightly wound spiral density waves can be determined locally (Lin & Shu 1964). In fact, in this approximation, known as Wentzel–Kramers–Brillouin (WKB) approximation, the relevant analysis is completely analytic and independent of the given model for the background disk. In other words, using this approximation, one can easily determine the stability of the given disk against axisymmetric perturbations (a density wave). At an arbitrary location in the stellar disk, stability requires that

$$\begin{eqnarray}&&{{Q}_{s}}=\frac{\kappa {{\sigma }_{r}}}{3.36G{{{\Sigma }}_{0}}}\gt 1\end{eqnarray} \tag{ 1 }$$

where κ is the epicyclic frequency, ${{\sigma }_{r}}$ is the radial velocity dispersion, ${{{\Sigma }}_{0}}$ is the disk unperturbed density, and G is the gravitational constant. The requirement (1) is known as Toomre's local stability criterion and has a simple physical interpretation. In fact this criterion represents a competition between the stabilizing influences of pressure (${{\sigma }_{r}}$) and the angular momentum (κ) against the destabilizing effect of gravity (${{{\Sigma }}_{0}}$). In the case of a fluid (gaseous) disk, this criterion takes the following form

$$\begin{eqnarray}&&{{Q}_{g}}=\frac{\kappa {{v}_{s}}}{\pi G{{{\Sigma }}_{0}}}\gt 1\end{eqnarray} \tag{ 2 }$$

where v_s is the sound speed. It is worth mentioning that it was first shown by Safronov (1960) that a thin fluid rotating disk can be unstable against local axisymmetric disturbances under the effect of its own gravity. These criteria, i.e., (1) and (2), have been modified to take into account various effects. For example the effect of the thickness of the disk has been investigated by Toomre (1964) and Vandervoort (1970). Two-component (fluid + stars) disks of finite thickness have been considered by Kato (1972), Bertin & Romeo (1988), Romeo (1992), Wang & Silk (1994) and Romeo & Wiegert (2011). See also Shadmehri & Khajenabi (2012) for a two-component study where the gaseous component is turbulent. Rafikov (2001) studied the general case of multiple components. The stability of a two-fluid disk was investigated by Jog & Solomon (1984), Elmegreen (1995) and Jog (1996). Effects of magnetic field on linear gravitational instabilities in thin galactic disks have been studied by Elmegreen (1987, 1994), Gammie (1996), Fan & Lou (1997), Kim & Ostriker (2001). Gammie (1996) has also considered the effect of viscosity on the local stability. Also the local stability in a two-component galactic disk with gas dissipation has been studied by Elmegreen (2011).

The purpose of this paper is to introduce a modified stability criterion in the context of Modified Gravity (MOG) for a single-component galactic disk (fluid or stellar). This theory is a covariant modification of Einstein's general relativity (GR) and has been introduced for addressing the dark matter problem (see Moffat 2006). In the weak field limit, the Poisson equation in this theory is different from that of Newtonian gravity. Therefore, one might naturally expect that the local gravitational stability criterion is different from (1) and (2). Finding the generalized version of Toomre's criterion in modified theories of gravity is not a new idea. For example, Milgrom (1989) has found the local stability criterion in modified Newtonian dynamics (MOND). Also a Toomre-like criterion in the context of metric f(R) gravity has been calculated for both gaseous and stellar disks by Roshan & Abbassi (2014b). Furthermore, this issue has been investigated in the context of some modified gravity theories which introduce a Yukawa like term in the gravitational force by Habibi et al. (2014). In that paper (Habibi et al. 2014) the local stability of disks in the context of MOG was also studied briefly as a special example for the above-mentioned theories. Also some approximate stability criteria have been derived. Here, we give a more careful treatment and find the exact form of the stability criteria and use the relevant observational data to find a detectable difference between Newtonian gravity and MOG.

In this section we briefly review the current status of MOG among other extended theories of gravity. Also, we briefly review its weak field limit, which is necessary to achieve the main goal of this paper.

As we mentioned in the introduction, MOG is a fully relativistic and covariant generalization of GR. Moffat (2006) introduced this theory as an alternative to dark matter theories. MOG is much more complicated than GR in the sense that its associated gravitational fields are more than GR. In fact, MOG is a scalar–tensor–vector theory of gravity, while GR is a tensor theory. In other words, in GR the components of the metric tensor are the gravitational fields of the theory. On the other hand, in MOG in addition to the metric tensor, there are three scalar fields ($\omega ({{x}^{\beta }}),\mu ({{x}^{\beta }})$, and $G({{x}^{\beta }})$) and also a massive Proca vector field ${{\phi }^{\beta }}$ (in the current literature of MOG, the scalar field ω is assumed to be constant). Therefore, MOG has been equipped with some degrees of freedom in order to handle the dark matter problem without invoking dark matter particles. This theory has been successfully applied to explain the rotation curves of spiral galaxies and the mass discrepancy in galaxy clusters (Brownstein & Moffat 2006, 2007; Moffat 2006; Moffat & Toth 2008, 2009, 2013; Moffat & Rahvar 2013, 2014).

For the details of deriving the weak field limit of MOG, we refer the reader to Moffat & Rahvar (2013) and Roshan & Abbassi (2014a). In the weak field limit of this theory, a particle's equation of motion can be written as

$$\begin{eqnarray}&&\frac{{{d}^{2}}{\boldsymbol{r}} }{d{{t}^{2}}}=-\nabla {\Phi }\end{eqnarray} \tag{ 3 }$$

where Φ is an effective gravitational potential and is given by

$$\begin{eqnarray}&&{\Phi }={\Psi }+\chi {{\phi }^{0}}\end{eqnarray} \tag{ 4 }$$

where χ is a coupling constant and ${{\phi }^{0}}$ is the zeroth component (or time component) of the vector field. Furthermore, Ψ and ${{\phi }^{0}}$ satisfy the following equations

$$\begin{eqnarray}&&{{\nabla }^{2}}{\Psi }=4\pi (1+\alpha )G\rho \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&({{\nabla }^{2}}-\mu _{0}^{2})\chi {{\phi }^{0}}=-4\pi \alpha G\rho \end{eqnarray} \tag{ 6 }$$

where ${{\mu }_{0}}$ and α are the free parameters of the theory in the weak field limit, G and ρ are the gravitational constant and the matter density, respectively. In fact, ${{\mu }_{0}}$ is the background value of the scalar field μ and α is related to the coupling constant χ as $\alpha ={{\chi }^{2}}/{{\omega }_{0}}G$. Note that ${{\omega }_{0}}$ is the background value of the scalar field ω. For moe details see Roshan & Abbassi (2014a). The observational values of the free parameters α and ${{\mu }_{0}}$ are known from rotation curve data of spiral galaxies. It has been shown by Moffat & Rahvar (2013) that the best values for these parameters are $\alpha =8.89\pm 0.34$ and ${{\mu }_{0}}=0.042\;\pm 0.004\;{\rm kp}{{{\rm c}}^{-1}}$. It is worth mentioning that Equations (5) and (6) can be combined to write the generalized Poisson equation as

$$\begin{eqnarray}&&{{\nabla }^{2}}{\Phi }=4\pi G\rho +\alpha \mu _{0}^{2}G\int \frac{{{e}^{-{{\mu }_{0}}|{\boldsymbol{r}} -{\boldsymbol{r}} ^{\prime} |}}}{|{\boldsymbol{r}} -{\boldsymbol{r}} ^{\prime} |}\rho ({\boldsymbol{r}} ^{\prime} ){{d}^{3}}x^{\prime} .\end{eqnarray} \tag{ 7 }$$

For considering the dynamics of a self-gravitating fluid system in the context of MOG, in addition to the generalized Poisson Equation (7), we need the generalized version of the continuity and Euler equations. It can be shown that the mathematical forms of these equations are the same as in Newtonian gravity (Roshan & Abbassi 2014a). The continuity and Euler equations are:

$$\begin{eqnarray}&&\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho {\boldsymbol{v}} )=0\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\frac{\partial {\boldsymbol{v}} }{\partial t}+({\boldsymbol{v}} \cdot \nabla ){\boldsymbol{v}} =-\frac{\nabla p}{p}-\nabla {\Phi }\end{eqnarray} \tag{ 9 }$$

where p is the pressure and ${\boldsymbol{v}}$ is the velocity of the fluid. Equations (7)–(9) combined with the equation of state of the fluid make a complete set of equations for describing the dynamics of a fluid system in the framework of MOG. Equivalently, one can use Equations (5) and (6) instead of (7).

On the other hand, in order to investigate the dynamics of a stellar system we need the Boltzmann equation and the generalized Poisson Equation (7). The Boltzmann equation in MOG has been investigated in Roshan & Abbassi (2014a) and Moffat & Rahvar (2014). Again, one can verify that the mathematical form of the Boltzmann equation does not change. The main difference is that the effective potential Φ appears in the Boltzmann equation instead of the Newtonian gravitational potential ${{{\Phi }}_{N}}$. Therefore, the collisionless Boltzmann equation in MOG is

$$\begin{eqnarray}&&\frac{\partial f}{\partial t}+{\boldsymbol{v}} \cdot \nabla f-\nabla {\Phi }\cdot \frac{\partial f}{\partial {\boldsymbol{v}} }=0\end{eqnarray} \tag{ 10 }$$

where f is the phase-space distribution function. In Section 4, we will use (10) for studying the local stability of a stellar disk.

In this section we study the behavior of a tightly wound density wave propagating in the surface of a fluid disk in the context of MOG. First, we use the modified Poisson Equation (7) in order to calculate the gravitational potential of a tightly wound surface density. Then, by linearizing the field equations, we find the relevant dispersion relation. We assume that the background disk is razor-thin and axisymmetric and its plane corresponds to the $x-y$ plane. Also, we assume that the disk has a barotropic equation of state as $p=K{{{\Sigma }}^{\gamma }}$, where K and γ are constant real parameters. The model assumes a non-rotating cylindrical system $(r,\varphi ,z)$, such that the z axis coincide with the rotation axis of the disk and the angle $\varphi$ increases in the direction of rotation. In cylindrical coordinates, by assuming that the fluid is axisymmetric, continuity and two components of the Euler equation can be written as

$$\begin{eqnarray}&&\frac{\partial {\Sigma }}{\partial t}+\frac{1}{r}\frac{\partial }{\partial r}\left( {\Sigma }r{{v}_{r}} \right)+\frac{1}{r}\frac{\partial }{\partial \varphi }\left( {\Sigma }{{v}_{\varphi }} \right)=0\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&\frac{\partial {{v}_{r}}}{\partial t}+{{v}_{r}}\frac{\partial {{v}_{r}}}{\partial r}+\frac{{{v}_{\varphi }}}{r}\frac{\partial {{v}_{r}}}{\partial \varphi }-\frac{v_{\varphi }^{2}}{r}=-\frac{\partial }{\partial r}\left( {\Phi }+h \right)\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&\frac{\partial {{v}_{\varphi }}}{\partial t}+{{v}_{r}}\frac{\partial {{v}_{\varphi }}}{\partial r}+\frac{{{v}_{\varphi }}}{r}\frac{\partial {{v}_{\varphi }}}{\partial \varphi }+\frac{{{v}_{\varphi }}{{v}_{r}}}{r}=-\frac{1}{r}\frac{\partial }{\partial \varphi }\left( {\Phi }+h \right)\end{eqnarray} \tag{ 13 }$$

where v_r, ${{v}_{\varphi }}$, and ${\Sigma }$ are the velocity components in the radial and azimuthal directions and the surface density, respectively. Furthermore, h is the specific enthalpy defined as $h=\int dp/{\Sigma }$. Now, we linearize Equations (11)–(13) by assuming that ${\Sigma }={{{\Sigma }}_{0}}+{{{\Sigma }}_{1}}$, ${{v}_{r}}={{v}_{r0}}+{{v}_{r1}}={{v}_{r1}}$, ${{v}_{\varphi }}={{v}_{\varphi 0}}+{{v}_{\varphi 1}}$, ${\Phi }={{{\Phi }}_{0}}+{{{\Phi }}_{1}}$ and $h={{h}_{0}}+{{h}_{1}}$. The subscript "0" refers to the background value of the given quantity and "1" refers to the corresponding perturbed quantity. After linearizing, the fluid Equations (11)–(13) can be written as (Binney & Tremaine 2008)

$$\begin{eqnarray}&&\frac{\partial {{{\Sigma }}_{1}}}{\partial t}+\frac{1}{r}\frac{\partial }{\partial r}\left( {{{\Sigma }}_{0}}r{{v}_{r1}} \right)+{\Omega }\frac{\partial {{{\Sigma }}_{1}}}{\partial \varphi }+\frac{{{{\Sigma }}_{0}}}{r}\frac{\partial {{v}_{\varphi 1}}}{\partial \varphi }=0\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\frac{\partial {{v}_{r1}}}{\partial t}+{\Omega }\frac{\partial {{v}_{r1}}}{\partial \varphi }-2{\Omega }{{v}_{\varphi 1}}=-\frac{\partial }{\partial r}\left( {{{\Phi }}_{1}}+{{h}_{1}} \right)\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&\frac{\partial {{v}_{\varphi 1}}}{\partial t}+{\Omega }\frac{\partial {{v}_{\varphi 1}}}{\partial \varphi }+\frac{{{\kappa }^{2}}}{2{\Omega }}{{v}_{r1}}=-\frac{1}{r}\frac{\partial }{\partial \varphi }\left( {{{\Phi }}_{1}}+{{h}_{1}} \right)\end{eqnarray} \tag{ 16 }$$

where ${\Omega }(r)$ is the angular rotation rate and the epicyclic frequency κ is defined as

$$\begin{eqnarray}&&\kappa (r)=\sqrt{r\frac{d{{{\Omega }}^{2}}}{dr}+4{{{\Omega }}^{2}}}.\end{eqnarray} \tag{ 17 }$$

Also, Equations (5) and (6) can be linearized as

$$\begin{eqnarray}&&{{\nabla }^{2}}{{{\Psi }}_{1}}=4\pi (1+\alpha )G{{{\Sigma }}_{1}}\delta (z)\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&({{\nabla }^{2}}-\mu _{0}^{2})\chi \phi _{1}^{0}=-4\pi \alpha G{{{\Sigma }}_{1}}\delta (z)\end{eqnarray} \tag{ 19 }$$

where $\delta (z)$ is the Dirac delta function. We choose infinitesimal disturbances of the form ${{Q}_{1}}={{Q}_{a}}{{e}^{i(kr+m\varphi -\omega t)}}$ where ω is the oscillation frequency and $k=\frac{2\pi }{\lambda }$ is the wavenumber. For the local stability analysis we use the WKB approximation (or tight winding approximation) which requires that $\frac{kr}{m}\gg 1$ (Binney & Tremaine 2008). This approximation allows us to neglect terms proportional to $1/r$ compared with the terms proportional to k. Employing this approximation, Equations (14)–(16) can be combined to give

$$\begin{eqnarray}&&(m{\Omega }-\omega ){{{\Sigma }}_{a}}+k{{{\Sigma }}_{0}}{{v}_{{\rm ra}}}=0\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{{v}_{{\rm ra}}}=\frac{(m{\Omega }-\omega )k({{{\Phi }}_{a}}+{{h}_{a}})}{{\Delta }}\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{{v}_{\varphi a}}=\frac{2iB}{\omega -m{\Omega }}{{v}_{{\rm ra}}}\end{eqnarray} \tag{ 22 }$$

where Δ and B (Oort's constant of rotation) are

$$\begin{eqnarray}&&{\Delta }={{\kappa }^{2}}-{{(m{\Omega }-\omega )}^{2}}\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&B(r)=-\frac{1}{2}\left( {\Omega }+\frac{d({\Omega }r)}{dr} \right).\end{eqnarray} \tag{ 24 }$$

In order to complete the local stability analysis, let us find the gravitational potential of a tightly wound spiral perturbation, ${{{\Sigma }}_{1}}={{{\Sigma }}_{a}}{{e}^{i(kr+m\varphi -\omega t)}}$ in the neighborhood of an arbitrary point $({{r}_{0}},{{\varphi }_{0}})$. Locally, the WKB density perturbation can be considered as a plane wave with wavenumber ${\boldsymbol{k}} =k{{{\boldsymbol{e}} }_{r}}$, where ${{{\boldsymbol{e}} }_{r}}$ is the unit vector in the radial direction. Therefore, finding the potential of a WKB perturbation reduces to finding the gravitational potential of a plane density wave. In order to find the gravitational potential of this perturbation, without loss of generality, we choose the x axis to be parallel to ${\boldsymbol{k}} ({{r}_{0}})$. Then we use Equations (18) and (19) to find the gravitational potential. To do so, we guess the solution to (19) as

$$\begin{eqnarray}&&\phi _{1}^{0}(x,y,z,t)=\phi _{a}^{0}\;{{e}^{i(kx-\omega t)-|\zeta z|}}\end{eqnarray} \tag{ 25 }$$

where $\phi _{a}^{0}$ and ζ are arbitrary constants. The disk is razor-thin and so there is no matter outside ($z\ne 0$) the disk. Therefore, at $z\ne 0$ Equation (18) can be written as ${{\nabla }^{2}}\phi _{1}^{0}=\mu _{0}^{2}\phi _{1}^{0}$. Substituting the potential (25) into this equation, one can easily show that $\zeta =\pm \sqrt{{{k}^{2}}+\mu _{0}^{2}}$. On the other hand, since matter is located at z = 0, the derivative of $\phi _{1}^{0}$ with respect to z is not continuous on the disk. In order to fix the parameter ζ, we integrate Equation (18) with respect to z in the interval $z=-\xi$ to $z=\;+\xi$, where ξ is a positive constant, and then let $\xi \to 0$. The result is

$$\begin{eqnarray}&&\phi _{1}^{0}=\frac{2\pi G\alpha }{\chi \sqrt{{{k}^{2}}+\mu _{0}^{2}}}{{{\Sigma }}_{a}}.\end{eqnarray} \tag{ 26 }$$

A similar procedure for Equation (19) gives

$$\begin{eqnarray}&&{{\varphi }_{a}}=-\frac{2\pi G(1+\alpha )}{|k|}{{{\Sigma }}_{a}}.\end{eqnarray} \tag{ 27 }$$

Finally, the gravitational potential of a plane wave in the context of MOG is given by

$$\begin{eqnarray}&&{{{\Phi }}_{a}}=-\frac{2\pi G}{|k|}{{{\Sigma }}_{a}}\left( 1+\alpha -\frac{|k|\alpha }{\sqrt{{{k}^{2}}+\mu _{0}^{2}}} \right).\end{eqnarray} \tag{ 28 }$$

Equations (20)–(22) and (28) combined with ${{h}_{a}}=v_{s}^{2}{{{\Sigma }}_{a}}/{{{\Sigma }}_{0}}$ form a closed set of linear equations. We solve these equations to obtain the dispersion relation. Using Equations (21) and (28) we can find a purely algebraic expression between v_ra and ${{{\Sigma }}_{a}}$ as

$$\begin{eqnarray}&&\begin{array}{l} {{v}_{{\rm ra}}}=\frac{(m{\Omega }-\omega )k}{{\Delta }}(\frac{v_{s}^{2}}{{{{\Sigma }}_{0}}}-\frac{2\pi G}{|k|}[1+\alpha -\frac{|k|\alpha }{\sqrt{{{k}^{2}}+\mu _{0}^{2}}}]){{{\Sigma }}_{a}}. \\ \end{array}\end{eqnarray} \tag{ 29 }$$

Substituting this equation into (20), and restricting ourselves to the axisymmetric (m = 0) perturbations, we obtain the following dispersion relation:

$$\begin{eqnarray}&&{{\omega }^{2}}={{\kappa }^{2}}+{{k}^{2}}v_{s}^{2}-2\pi G{{{\Sigma }}_{0}}|k|\left( 1+\alpha -\frac{|k|\alpha }{\sqrt{{{k}^{2}}+\mu _{0}^{2}}} \right).\end{eqnarray} \tag{ 30 }$$

We derive the local gravitational stability criterion by using this dispersion relation. As expected, in the limit ${{\mu }_{0}}\to 0$ or $\alpha \to 0,$ the dispersion relation for WKB perturbations in Newtonian gravity is recovered. Since all the quantities in the right-hand side (rhs) of (30) are real quantities, then the disk is stable against local disturbances if ${{\omega }^{2}}\gt 0$ and unstable if ${{\omega }^{2}}\lt 0$.

In order to find the local stability criterion, we write the dispersion relation (30) in a dimensionless form. To do so, we define the dimensionless wavenumber X, and the parameter β as

$$\begin{eqnarray}&&X=\frac{k}{{{\mu }_{0}}},\;\;\;\;\beta =\frac{{{\mu }_{0}}{{v}_{s}}}{\kappa },\end{eqnarray} \tag{ 31 }$$

therefore the dispersion relation can be rewritten in the dimensionless form

$$\begin{eqnarray}&&\frac{{{\omega }^{2}}}{{{\kappa }^{2}}}=1+{{\beta }^{2}}{{X}^{2}}-\frac{2\beta |X|}{{{Q}_{g}}}(1+\alpha )+\frac{2\beta }{{{Q}_{g}}}\frac{\alpha {{X}^{2}}}{\sqrt{1+{{X}^{2}}}}\end{eqnarray} \tag{ 32 }$$

where Q_g is the dimensionless Toomre parameter given by Equation (2). Therefore, the stability criterion ${{\omega }^{2}}\gt 0$ takes the following form:

$$\begin{eqnarray}&&{{Q}_{g}}\gt \frac{2\beta }{1+{{\beta }^{2}}{{X}^{2}}}\left( (1+\alpha )|X|-\frac{\alpha {{X}^{2}}}{\sqrt{1+{{X}^{2}}}} \right).\end{eqnarray} \tag{ 33 }$$

This is the exact criterion for the stability of the fluid disk against a perturbation mode with given wavenumber X. Once we have the observational values of α and ${{\mu }_{0}}$, and also the magnitude of the wavenumber X, then we can explicitly determine the required value of Toomre's parameter for providing the stability at the given point. More specifically, we should know the magnitude of the ratio ${{v}_{s}}/\kappa$ at the given location. It is natural that if we find the maximum value of the rhs of (33) for given α and β, then we will get a criterion that is independent of the value of the wavenumber. Since this cannot be done analytically, we have listed the local stability criterion for various values of α and β in Table 1. It is worth remembering that, by fitting MOG to observed rotation curves of galaxies, the free parameters α and ${{\mu }_{0}}$ have been determined to have the values: $\alpha =8.89\pm 0.34$ and ${{\mu }_{0}}=0.042\;\pm 0.004\;{\rm kp}{{{\rm c}}^{-1}}$. So in Table 1, we have chosen α as 8.89 ± 0.34. On the other hand, in order to choose rational values for β, we have estimated it in the solar neighborhood. In the solar neighborhood the epicycle frequency is $\kappa \sim 37\;{\rm km}\;{{{\rm s}}^{-1}}\;{\rm kp}{{{\rm c}}^{-1}}$, and the radial velocity dispersion is ${{\sigma }_{r}}\sim 38\;{\rm km}\ {{{\rm s}}^{-1}}$. Therefore $\beta \sim 0.04$. It is obvious from Table 1 that in the solar neighborhood there is no significant difference between the local stability criterion of MOG and Newtonian gravity. Also, it is clear that the difference between these theories considering the local stability will be significant in the locations where the ratio of the sound speed to the epicycle frequency is large. In fact, ${{v}_{s}}/\kappa$ is an important length scale for comparing MOG with Newtonian gravity in the disk galaxies. If this length scale satisfies the inequality ${{v}_{s}}/\kappa \gtrsim 2.4\;{\rm kpc}$, the deviation between the local stability criterion of MOG and Newtonian gravity will be large. At the end of this section, we will discuss the relevant observational data for this length scale.

It is interesting to plot the dispersion relation for the growing (unstable) modes and discuss the growing rates. Our purpose is to compare the growth rates in MOG and Newtonian gravity. To do so, we rewrite the dispersion relation (30) as follows:

$$\begin{eqnarray}&&s{{^{\prime} }^{2}}=\frac{2|q|}{{{Q}_{g}}}\left( 1+\alpha -\frac{\alpha |q|}{\sqrt{{{q}^{2}}+{{\beta }^{2}}}} \right)-(1+{{q}^{2}})\end{eqnarray} \tag{ 34 }$$

where the dimensionless wavenumber q is defined as $q=k{{v}_{s}}/\kappa ,$ and $s^{\prime} =i\omega /\kappa$. Figure 1 shows the fluid dispersion relation for various values of Q_g and β. The dotted curves are for Newtonian gravity. The maximum growth rate in MOG is higher than that of Newtonian gravity, and also occurs at smaller wavenumber. As in Newtonian gravity (dotted curves), the growth rates increase with decreasing Q_g. It is also clear from Figure 1 that in MOG the range of instability is stretched to slightly larger wavenumber. Furthermore, by increasing the parameter β the range of instability and also the maximum growth rate increase.

Zoom In Zoom Out Reset image size

It is worth mentioning that in metric f(R) gravity, by increasing the corresponding parameter β, the growth rate decreases (Roshan & Abbassi 2014b). In this theory the parameter β can be defined as $\beta ={{m}_{0}}{{v}_{s}}/\kappa$, where m₀ is one of the free parameters. It is surprising that, although both MOG and f(R) gravity lead to a stronger gravitational force relative to Newtonian gravity, MOG increases the growth rate while f(R) gravity decreases it.

It is also instructive to illustrate the neutral stability curves, i.e., $\omega =0$, for axisymmetric tightly wound spiral density waves in a fluid disk. These curves determine the boundary of the stable and unstable perturbations. Let us define a critical wavelength ${{\lambda }_{{\rm crit}}}=2\pi /{{k}_{{\rm crit}}}$ where ${{k}_{{\rm crit}}}={{\kappa }^{2}}/2\pi G{{{\Sigma }}_{0}}$. In Newtonian gravity ${{\lambda }_{{\rm crit}}}$ is the largest unstable wavelength for a fluid disk with zero sound speed. Using this definition and the dispersion relation (30), the line separating stable and unstable perturbations is the solution of the following equation

$$\begin{eqnarray}&&{{Q}_{g}}(y)=\sqrt{4y[1+\alpha -\frac{\alpha }{\sqrt{1+{{\left( \frac{2y\beta }{{{Q}_{g}}} \right)}^{2}}}}]-4{{y}^{2}}}\end{eqnarray} \tag{ 35 }$$

where $y=\lambda /{{\lambda }_{{\rm crit}}}$ and $\lambda =2\pi /k$. We have shown the neutral curves for various values of β in Figure 2. Furthermore, in this figure, we have assumed $\alpha =9.23$. The dotted curve corresponds to Newtonian gravity. Also, the dashed curves show the neutral curves of a fluid disk and solid lines are the corresponding lines for a stellar disk. It is clear that, at small wavelengths, the neutral curves of MOG and Newtonian gravity coincide. This coincidence is almost true for the interval $0\lt y\lt 0.5$, i.e., $\lambda \lt 0.5{{\lambda }_{{\rm crit}}}$. However, for $y\gt 0.5$ deviations between MOG and Newtonian gravity appear. More specifically, the range of instability is stretched to longer wavelength. It is surprising that the longest unstable wavelength in MOG is more than 10 times the Newtonian case. In fact one can easily show that the longer unstable mode in MOG occurs at $\lambda =(1+\alpha ){{\lambda }_{{\rm crit}}}$. This is quite a big difference between these theories. However, one should note that at long wavelengths, there is a serious doubt about the validity of the WKB approximation. For example, for the solar neighborhood ${{\lambda }_{{\rm crit}}}\sim 10\;{\rm kpc}$. This wavelength is longer than the distance from the sun to the center of the Galaxy, $r\sim 8.5\;{\rm kpc}$. Therefore ${{k}_{{\rm crit}}}r\sim 5.3$ and for the largest unstable wavelength in MOG we have $kr\sim 0.5$. Thus, obviously, the WKB approximation is not satisfied. Therefore, the width of the range of instability in MOG compared with Newtonian gravity cannot be considered as a significant and detectable difference.

Zoom In Zoom Out Reset image size

On the other hand, it is clear in Figure 2 that increasing the parameter β leads to a higher Q_g for stability of all wavelengths. In other words, in MOG, a higher Q_g parameter compared with Newtonian gravity is needed to overcome the local gravitational instability. This result is also consistent with Table 1. As we discussed, this difference between MOG and Newtonian gravity could lead to a detectable deviation.

The origin of this difference can be simply understood. Using the modified Poisson Equation (7), the gravitational force between two point masses m₁ and m₂ that are located at ${{{\boldsymbol{r}} }_{{\bf 1}}}$ and ${{{\boldsymbol{r}} }_{2}}$ respectively, is given by

$$\begin{eqnarray}&&{\boldsymbol{F}} ={{{\boldsymbol{F}} }_{N}}+\alpha {{{\boldsymbol{F}} }_{N}}(1-{{e}^{-{{\mu }_{0}}r}}(1+{{\mu }_{0}}r))\end{eqnarray} \tag{ 36 }$$

where ${{{\boldsymbol{F}} }_{N}}=-G{{m}_{1}}{{m}_{2}}/{{r}^{2}}{{{\boldsymbol{e}} }_{r}}$ is the Newtonian gravitational force and ${\boldsymbol{r}} ={{{\boldsymbol{r}} }_{{\bf 1}}}-{{{\boldsymbol{r}} }_{{\bf 2}}}$ and ${{{\boldsymbol{e}} }_{r}}$ is the unit vector in the direction of ${\boldsymbol{r}}$. Taking into account the current values of α and ${{\mu }_{0}}$, one can easily verify that ${\boldsymbol{F}}$ is stronger than Newton's gravitational force ${{{\boldsymbol{F}} }_{N}}$. Now, keeping in mind the physical meaning of the gravitational Jeans instability, it is natural to expect more gas pressure to prevent the gravitational collapse. In other words, it seems like that, when the gravitational force is stronger than Newtonian gravity, then we need a higher Q_g parameter. However, this interpretation is not always true. For example, in the case of metric f(R) gravity, where the corresponding gravitational force in the weak field limit is stronger than Newtonian gravity, the required value for Q_g is smaller than the Newtonian one.

The stability of stellar disks was first studied by Toomre (1964), Kalnajs (1965), and Lin & Shu (1966) who pointed out that in certain regimes the stellar disk is stable against disturbances. When dealing with stellar disk stability, we need to use the collisionless Boltzmann Equation (10), as well as the modified Poisson Equation (7). Since we assume that the stellar disk is collisionless, the stars cannot create any pressure and consequently the sound speed for this system is zero. However, because of the epicycle motion at any given patch in the disk there are stars from different parts of the perturbed disk. This causes an important cancellation effect (Toomre 1964). In fact, it turns out that the stellar disk can be described by the Jeans equations; their mathematical form is rather similar to the hydrodynamic Equations (11) and (9). However, although Equation (20) does not change, Equation (21) is replaced with (see Roshan & Abbassi 2014b for more detail)

$$\begin{eqnarray}&&{{v}_{{\rm ra}}}=\frac{m{\Omega }-\omega }{{\Delta }}k{{{\Phi }}_{a}}\mathcal{F}\end{eqnarray} \tag{ 37 }$$

where $\mathcal{F}$ is the reduction factor given by (Binney & Tremaine 2008)

$$\begin{eqnarray}&&\mathcal{F}(s,{{q}^{2}})=\frac{1-{{s}^{2}}}{{\rm sin} \pi s}\int _{0}^{\pi }{{e}^{-{{q}^{2}}(1+{\rm cos} \tau )}}{\rm sin} s\tau {\rm sin} \tau d\tau \end{eqnarray} \tag{ 38 }$$

where $s=-is^{\prime} =(\omega -m{\Omega })/\kappa$, and $q=k{{\sigma }_{r}}/\kappa$ and ${{\sigma }_{r}}$ are the dimensionless wavenumber and the radial velocity dispersion, respectively. Since the reduction factor does not depend on the form of Poisson equation (Roshan & Abbassi 2014b), in MOG it possesses the same general form and exhibits similar dependencies on dimensionless frequency and wavenumber as it has for Newtonian gravity. Now, substituting (37) into Equation (20), we get the following dispersion relation

$$\begin{eqnarray}\begin{array}{rcl} {{(\omega -m{\Omega })}^{2}} & = & {{\kappa }^{2}}-2\pi G{{{\Sigma }}_{0}}|k|\left[ 1+\alpha \right. \\ {} & {} & \left. -\frac{\alpha |k|}{\sqrt{{{k}^{2}}+\mu _{0}^{2}}} \right]\mathcal{F}(s,{{q}^{2}}). \\ \end{array}\end{eqnarray} \tag{ 39 }$$

At the limit $\alpha \to 0$ or ${{\mu }_{0}}\to 0$, the standard dispersion relation is recovered. Again, we restrict ourselves to axisymmetric perturbations. The system is stable against all axisymmetric modes if Equation (39) does not have a solution with negative ${{\omega }^{2}}$. In order to check this expectation, we rewrite dispersion relation (39) as follows:

$$\begin{eqnarray}&&1=\frac{2\pi |q|}{3.36{{Q}_{s}}}\left( 1+\alpha -\frac{\alpha |q|}{\sqrt{{{q}^{2}}+{{\beta }^{2}}}} \right)\frac{\mathcal{F}(-is^{\prime} ,{{q}^{2}})}{1+s{{^{\prime} }^{2}}}.\end{eqnarray} \tag{ 40 }$$

In this case $\beta ={{\mu }_{0}}{{\sigma }_{r}}/\kappa$. It is straightforward to verify that the maximum value of the rhs of (40) occurs at $s^{\prime} =0$. Therefore, if the rhs of Equation (40) for $s^{\prime} =0$ is smaller than 1, then there is no perturbation with negative ${{\omega }^{2}}$ that satisfies the dispersion relation (40). Consequently, the disk will be stable against all perturbations. Therefore, the stability criterion reads

$$\begin{eqnarray}&&{{Q}_{s}}\gt \frac{2\pi |q|}{3.36}\left( 1+\alpha -\frac{\alpha |q|}{\sqrt{{{q}^{2}}+{{\beta }^{2}}}} \right)\mathcal{F}(0,{{q}^{2}}).\end{eqnarray} \tag{ 41 }$$

This criterion can be compared with the stability criterion for a fluid disk, i.e., Equation (33). The criterion (41) can be replaced with a more simplified form:

$$\begin{eqnarray}&&{{Q}_{s}}\gt \frac{2\pi }{3.36|q|}\left( 1+\alpha -\frac{\alpha |q|}{\sqrt{{{q}^{2}}+{{\beta }^{2}}}} \right)\left( 1-{{e}^{-{{q}^{2}}}}{{I}_{0}}({{q}^{2}}) \right)\end{eqnarray} \tag{ 42 }$$

where I₀ is the modified Bessel function of order 0. In deriving (42), we have used the following identity

$$\begin{eqnarray}&&\mathcal{F}(0,{{q}^{2}})=\frac{1}{{{q}^{2}}}\left( 1-{{e}^{-{{q}^{2}}}}{{I}_{0}}({{q}^{2}}) \right).\end{eqnarray} \tag{ 43 }$$

Equation (42) is the main result of this section. For a given perturbation with wavenumber q, if we know the parameters α and β then we will easily calculate the required value of Q_s for stability of that perturbation. On the other hand, if we could find the maximum value of the rhs of Equation (42), then we should obtain a criterion that guarantees stability against all wavelengths at every location on the disk. In Table 1, the maximum value of the rhs has been calculated for various values of α and β. The tabulated data show that the required values of Q_s for stabilizing the disk are larger than those for the Newtonian case. As we mentioned before, this fact is also the case for a fluid disk.

In order to compare the growth rates in fluid and stellar disks, we have plotted the stellar disk dispersion relation (40) in Figure 3. This figure can be compared with Figure 1. The solid curves correspond to Q_s = 0.3 and $\beta =0.04$–0.2. The dashed curves correspond to Q_s = 0.5 and $\beta =0.04$–0.2. Furthermore, the dotted lines are the corresponding Newtonian dispersion relations for the above-mentioned values of Q_s. As in the case of a fluid disk, the growth rates increase with decreasing parameter Q_s. Also larger β leads to larger growth rate for long wavelengths. In fact, from Figure 3 it is evident that for short wavelengths there is no difference between the growth rates in MOG and Newtonian gravity. This point is also completely the same as for the fluid disk.

Zoom In Zoom Out Reset image size

It is also instructive to plot the boundary between the stable and unstable modes. As we mentioned before, this curve can be specified by the requirement ${{\omega }^{2}}=0$. In this case, the dispersion relation (40) can be rewritten as

$$\begin{eqnarray}&&\begin{array}{l} \frac{f({{Q}_{s}},y)y}{1-{{e}^{f({{Q}_{s}},y)}}{{I}_{0}}(f({{Q}_{s}},y))}=\left[ 1+\alpha \right.\left. -\frac{\alpha }{\sqrt{1+{{\beta }^{2}}/f({{Q}_{s}},y)}} \right] \\ \end{array}\end{eqnarray} \tag{ 44 }$$

where $f({{Q}_{s}},y)={{\left( \frac{3.36{{Q}_{s}}}{2\pi y} \right)}^{2}}$. The curves corresponding to (44) have been plotted in Figure 2 with solid lines for different values of β (0.02–0.2). It is obvious that, similar to the fluid case, the range of instability is stretched to longer wavelength. The longest unstable wavelength is $\lambda =10.23{{\lambda }_{{\rm crit}}}$ (note that we have chosen $\alpha =9.23$). However as we discussed in the previous section, this obvious deviation from the standard case does not necessarily lead to a significant difference between MOG and Newtonian gravity.

As we show in Sections 3 and 4, the deviation between the stability criteria of Newtonian gravity and MOG is highly controlled by the magnitude of parameter β. Therefore, in order to compare these theories, it is essential to find the magnitude of this parameter with respect to radius from observational data in some real spiral galaxies. To do so, we use a sample of six galaxies of THINGS catalog of spiral galaxies (Leroy et al. 2008). We have chosen this sub-sample of galaxies just for illustration, and our results in this section are not restricted to these galaxies. For every galaxy in this catalog, the rotation curve v(r), epicycle frequency $\kappa (r)$ and the velocity dispersion ${{\sigma }_{r}}(r)$ and v_s are known. So we can easily calculate β with respect to r for both the stellar and gaseous components of the galaxies.

Before computing β from observational data, it would be useful to make an order of magnitude estimation. It is worth mentioning that, at large radii in spiral galaxies, the rotation curve is flat, i.e., $v(r)\sim {{v}_{{\rm flat}}}$ where ${{v}_{{\rm flat}}}$ is constant. Therefore, ${\Omega }(r)\sim \frac{{{v}_{{\rm flat}}}}{r}$, and consequently $\kappa (r)=\sqrt{2}\frac{{{v}_{{\rm flat}}}}{r}$. For the gaseous component of the galaxies, v_s is almost constant with respect to r and is ${{v}_{s}}\sim 11\;{\rm km}\ {{{\rm s}}^{-1}}$. In fact we have assumed that the sound speed is equal to the gas velocity dispersion ${{\sigma }_{{\rm gas}}}$, see Leroy et al. (2008). Thus at large radii β rises linearly with radius as

$$\begin{eqnarray}&&\beta (r)=\frac{{{\mu }_{0}}{{v}_{s}}}{\sqrt{2}{{v}_{{\rm flat}}}}r.\end{eqnarray} \tag{ 45 }$$

This means that the magnitude of β is greater near the edge of the disk galaxy and consequently higher values for Q_g are required to overcome the gravitational instability there. An immediate and important consequence follows from this equation. Near the edge of the galactic disk the difference between Newtonian gravity and MOG might be significant. Therefore, the outer disks of spiral galaxies are exposed to instability more than the inner disks. Consequently, there might be some regions near the edge of the disk which can be stable in Newtonian gravity but unstable in the context of MOG. Such a deviation between these theories at large radii is expected because the free parameter $\mu _{0}^{-1}\sim 24\;{\rm kpc}$ is a large length scale compared to the characteristic size of the galaxies. As a consequence, taking into account the presence of the Yukawa term in this theory, we know that the effects of MOG appear at large radii.

From Table 1, it is evident that the stability criterion is significantly different from the standard case if $\beta \gtrsim 0.2$. Using Equation (45) the requirement $\beta \gtrsim 0.2$ reads

$$\begin{eqnarray}&&r\gtrsim 6.7\left( \frac{{{v}_{{\rm flat}}}}{{{v}_{s}}} \right){\rm kpc}.\end{eqnarray} \tag{ 46 }$$

If we set an average value $\frac{{{v}_{{\rm flat}}}}{{{v}_{s}}}\sim 10$, then deviation between Newtonian dynamics and MOG will appear at $r\gtrsim 67\;{\rm kpc}$. On the other hand, almost for all spiral galaxies, the characteristic optical radius is smaller than $67\;{\rm kpc}$. Therefore, this quick estimation shows that it is unlikely to find a galaxy where the predictions of MOG are substantially different from those of Newtonian gravity.

However, let us do a more careful analysis by finding function $\beta (r)$ for our sample. The rotation curve v(r) for the spiral galaxies of THINGS catalog can be written as (Boissier et al. 2003)

$$\begin{eqnarray}&&v(r)={{v}_{{\rm flat}}}\left[ 1-{\rm exp} \left( \frac{-r}{{{l}_{{\rm flat}}}} \right) \right]\end{eqnarray} \tag{ 47 }$$

where ${{v}_{{\rm flat}}}$ and ${{l}_{{\rm flat}}}$ are free parameters that refer to the velocity at which the rotation curve is flat and the length scale over which it approaches this rotation speed, respectively. These parameters are not universal and are different from galaxy to galaxy. For these parameters we use Table 4 of Leroy et al. (2008). Also, using the rotation speed v(r) the epicycle frequency can be written as

$$\begin{eqnarray}&&\kappa (r)=1.41\frac{v(r)}{r}\sqrt{1+b}\end{eqnarray} \tag{ 48 }$$

where $b=d{\rm ln} v(r)/d{\rm ln} r$. Furthermore, the radial velocity dispersion is given by (van der Kruit 1988; Leroy et al. 2008)

$$\begin{eqnarray}&&{{\sigma }_{r}}(r)=1.55\sqrt{G{{l}_{*}}{{{\Sigma }}_{*}}(r)}\end{eqnarray} \tag{ 49 }$$

where ${{{\Sigma }}_{*}}$ is the stellar surface density and l_* is the stellar scale length. We use Table 4 in Leroy et al. (2008) for l_*. For THINGS galaxies, ${{{\Sigma }}_{*}}$ has been measured with respect to radius in Leroy et al. (2008), for example see Table 7. Finally, using Equations (47)–(49), we have calculated β for both the gaseous and stellar components of the sample. The result has been illustrated in Figure 4. The red circles show $\beta (r)$ for the gaseous component and green circles correspond to the stellar component. As we have already expected, for the gaseous component β at large radii increases linearly with r. However, as illustrated in Figure 4, although β increases with radius, it remains small even at the outer disks. This happens also for the stellar component. In other words, β remains small for the stellar component. Of course, the behavior of β with respect to r for the stellar component is completely different from that of the gaseous component. In fact, MOG predicts that in outer disks the gaseous component must be more unstable against local perturbations than the stellar disk. However, as we mentioned above, this difference is not too big to be significant and detectable.

Zoom In Zoom Out Reset image size

In order to make the comparison between Newtonian gravity and MOG more clear, let us rewrite the stability criterion as follows:

$$\begin{eqnarray}&&{\Sigma }(r)\lt {{{\Sigma }}_{{\rm crit}}}(r).\end{eqnarray} \tag{ 50 }$$

In other words, a local axisymmetric perturbation at r is stable against gravitational collapse if the background matter density at that point is smaller than the critical density ${{{\Sigma }}_{{\rm crit}}}(r)$. Here, ${{{\Sigma }}_{{\rm crit}}}(r)$ in Newtonian gravity for gaseous and stellar disks is respectively

$$\begin{eqnarray}&&{{({{{\Sigma }}_{{\rm crit}}})}_{N}}=\frac{\kappa {{v}_{s}}}{\pi G},\;\;{{({{{\Sigma }}_{{\rm crit}}})}_{N}}=\frac{\kappa {{\sigma }_{r}}}{3.36G}.\end{eqnarray} \tag{ 51 }$$

On the other hand using criteria (41) and (33), one can write the critical surface density in MOG as follows:

$$\begin{eqnarray}&&{{({{{\Sigma }}_{{\rm crit}}})}_{{\rm MOG}}}=\frac{{{({{{\Sigma }}_{{\rm crit}}})}_{N}}}{H(\alpha ,\beta )}\end{eqnarray} \tag{ 52 }$$

where $H\gt 1$ and is a function of α and β. Note that for our purpose here, the exact form of the function H is not needed. Therefore, we can conclude that the critical surface density in the context of MOG is smaller than the Newtonian gravity. This makes a disk more unstable than it would be in the standard case. In order to see the difference between ${{({{{\Sigma }}_{{\rm crit}}})}_{{\rm MOG}}}$ and ${{({{{\Sigma }}_{{\rm crit}}})}_{N}}$, we have plotted their ratio in Figure 5. Although for the gaseous component the ratio ${{({{{\Sigma }}_{{\rm crit}}})}_{N}}/{{({{{\Sigma }}_{{\rm crit}}})}_{{\rm MOG}}}$ increases linearly, it does not get significantly larger than 1. In fact this ratio can be approximated as

$$\begin{eqnarray}&&\frac{{{({{{\Sigma }}_{{\rm crit}}})}_{N}}}{{{({{{\Sigma }}_{{\rm crit}}})}_{{\rm MOG}}}}\sim 1+\beta .\end{eqnarray} \tag{ 53 }$$

Since β remains small even at the outer disks, then one can conclude that the threshold matter density for gravitational collapse in Newtonian gravity and MOG is almost the same. Therefore, from this perspective, it seems that there is no difference between these theories. It should be stressed that we have done the same analysis for all galaxies in the THINGS catalog. And the result of this section is applicable to all of them.

Zoom In Zoom Out Reset image size

In this paper, the local gravitational stability of galactic disks has been investigated in the context of MOG. The dispersion relation in the WKB approximation has been derived for both gaseous and stellar disks in the context of MOG. Using the relevant dispersion relations, we have found the local stability criteria for gaseous and stellar disks.

Finally, we have used a sample of spiral galaxies of the THINGS catalog in order to compare these theories using the current observational data. Our final results have been summarized in the Figures 4 and 5. In Figure 5, we have compared the critical surface density over which the local perturbations collapse, in Newtonian gravity and MOG. For all galaxies of the sample, although the ratio of ${{({{{\Sigma }}_{{\rm crit}}})}_{N}}/{{({{{\Sigma }}_{{\rm crit}}})}_{{\rm MOG}}}$ increases with radius, it does not differ substantially from unity. Therefore, at least for all the spirals in the THINGS catalog, there is no significant difference between Newtonian gravity and MOG. It is also important to note that it could be more precise to find the stability criterion for a two-fluid disk (gas + stars). However, taking into account the main results of this paper, it is unlikely that a two-fluid analysis changes the results.

Of course, our results in this paper do not mean that the dynamics of spiral galaxies are the same in Newtonian gravity and MOG. Only the fate of local perturbations is the same in these theories. It is worth remembering that MOG can explain the rotation curves of the spiral galaxies without using the dark matter halos. However, this cannot be achieved in Newtonian gravity without invoking dark matter particles. Also, one should note that the global stability analysis of the disk galaxies in MOG, in principle, can be different from the standard case, for example see Brandao & de Araujo (2010) for a numerical study of disk galaxies in MOG.

As a final remark, let us mention that there is a well-known instability in the dynamics of spiral galaxies called bar instability. In fact, in the early 1970s, the N-body simulations of self-gravitating disks showed that a cold rotationally supported disk is globally unstable to the formation of a non-axisymmetric system. In other words, a rotation-dominated disk of particles will undergo a rapid transition to a bar-like pressure-dominated system (for example see Hohl 1971). This result was puzzling because real galaxies like the Milky Way do not show such a behavior. However, it is known that the presence of a dark matter halo can help to stabilize the disk, see Ostriker & Peebles (1973). Furthermore, the presence of a dense bulge-like mass component near the center of the disk can also stabilize a cold self-gravitating disk, see Sellwood (2014) for a comprehensive review of the subject. Another approach to this problem can be provided via modified gravity theories. In alternative theories of dark matter there is no dark matter halo and the bar instability problem should be addressed without invoking dark matter. For example see Brada & Milgrom (1999) and Tiret & Combes (2007). In these papers the global stability of the disk galaxies has been investigated in MOND. As we mentioned before, MOG is an alternative theory of gravity presented for addressing the dark matter problem. Therefore, it would be interesting to study the bar instability problem in this theory. The result of the current paper can be used in numerical study of disk galaxies in MOG. More specifically, in order to set the initial conditions of the particles in an N-body simulation for studying the global stability of the disk, the local stability criterion on the disk is required to be satisfied.

We would like to thank the referee for useful comments. This work is supported by Ferdowsi University of Mashhad under Grant No. 100835 (25/05/1393).