Stellar Bar Evolution in the Absence of Dark Matter Halo

Mahmood Roshan

doi:10.3847/1538-4357/aaaaad

1. Introduction

It is well-known that bulgeless numerical galactic models that are initially in a rotational equilibrium state are violently unstable to a pressure-dominated system; for primary simulations, see Miller et al. (1970) and Hohl (1971), for example. This fact is known in the literature as the bar instability. More specifically, the disk shape changes into a rotating bar over a timescale that is small compared to the rotational period of the outermost particles. However, it was understood that the existence of a rigid spherical component can suppress the bar instability (Ostriker & Peebles 1973). After Ostriker & Peebles' influential work, the role of dark matter halo on the global stability of galaxies has been extensively investigated; see, for example, Sellwood (1981, 2014b), Efstathiou et al. (1982), Athanassoula & Sellwood (1986), and Athanassoula (2002) for a recent review of the subject.

Although, historically, the dark matter halo was introduced to prevent bar formation, observations show that one-third of disk galaxies are strongly barred. On the other hand, one-third of them are weakly barred, and the rest do not contain a bar. In fact, one of the main unsolved problems with these galaxies is the frequency of bars. There have been several attempts to explain the presence or absence of a bar; for a very short review on the subject, see Berrier & Sellwood (2016). In this paper, we do not attempt to present a new explanation for the observed bar frequency. We are interested in bar evolution in models where there is no dark matter halo and where the gravitational force is different from that in the Newtonian case.

It is important to mention that the bar instability is directly linked to the dark matter problem. Conversely, modified gravity is another path to address this problem. Therefore, the global stability of disk galaxies is an important issue that might help distinguish between dark matter and modified gravity. In fact, dark matter particles have not yet been observed although there are several laboratories throughout the world searching for these particles using completely different techniques; see Bertone et al. (2010). This means that modified gravity is still an important approach to address this problem. These theories can be divided into two main categories: dark energy models and alternative theories that dispense with the dark matter paradigm. For a review on dark energy models, we refer the reader to Capozziello & De Laurentis (2011). On the other hand, Modified Newtonian Dynamics (MOND) is one of the well-studied theories modified for dark matter (Milgrom 1983; Famaey & McGaugh 2012). The scalar–tensor–vector theory of gravity known as MOG in the literature is another theory that has been presented to address the dark matter problem (Moffat 2006).

Before explaining our aim in this paper, let us briefly mention some well-established results regarding the stability of disk galaxies in both viewpoints, i.e., in particle dark matter and modified gravity. It is well understood that disks in live dark matter halos form bars more readily than those in rigid halos (Athanassoula 2002; Sellwood 2016). In fact, the angular momentum transfer between the bar and the halo, in principle, can increase the growth rate. This means that the so-called disk stability criteria derived using rigid halos (Ostriker & Peebles 1973; Efstathiou et al. 1982) cannot be used for real galaxies (Athanassoula 2008). The rotation of the halo can also influence the growth rate. For example, Saha & Naab (2013) reported that if the halo rotates in the same direction as the disk, bar growth will increase. On the other hand, counterrotating halos may suppress the instability.

Another important result of stellar N-body simulations is that the bar pattern speed decreases with time (Debattista & Sellwood 2000; Algorry et al. 2017). However, observations do not confirm this point (Aguerri et al. 2015). Therefore, it is still a challenge for standard dark matter halo models.

On the other hand, few investigations have been done on the stability of disk galaxies in modified gravity theories. The local stability of disk galaxies has been investigated in Milgrom (1989). The global stability of disk galaxies in MOND has been investigated using N-body simulations in Christodoulou (1991), Brada & Milgrom (1999), and Tiret & Combes (2007). Using a low-resolution simulation, it has been shown in Christodoulou (1991) that disk galaxies are more stable in MOND than in equivalent models in Newtonian gravity. A similar result has been reported in Brada & Milgrom (1999). However, advanced simulations with a larger number of particles show that the bar growth rate in MOND is larger than in Newtonian models (Tiret & Combes 2007). This means that galactic disks without a dark matter halo in MOND develop a bar sooner than in Newtonian models with a dark matter halo. However, in MOND, the bar starts to weaken, while in the dark matter model it does not stop growing. In this case, the final magnitude of the bar at the end of the simulation in MOND is smaller than the dark matter halo models.

Roshan & Abbassi (2015) investigated the local stability of galactic disks in MOG. In other words, the generalized version of the Toomre criterion has been found and applied to a sample of spiral galaxies. The result shows that there is no significant difference between MOG and the standard picture in the local stability issue. More specifically, the effects of MOG appear at large distances and, naturally, does not affect the dynamics of local perturbations. Roshan et al. (2016) investigated the global stability of the Maclaurin disk in MOG using a semi-analytic method. On the other hand, the global stability of exponential and Mestel disk galaxies in MOG has been investigated in Ghafourian & Roshan (2017) using low-resolution N-body simulations. Although there is no meaningful difference in the local stability of disk galaxies in MOG and the dark matter model, the global stability is substantially different. More specifically, it was found in Ghafourian & Roshan (2017) that disks are more stable in MOG, and the growth rate is smaller.

However, the maximum number of particles in the simulations of Ghafourian & Roshan (2017) is N = 2 × 10⁴. Therefore, the results may suffer from shot noise and numerical artifacts. Although some tests have been done to ensure the reliability of the results, it still seems necessary to perform high-resolution simulations. Furthermore, galactic models in Ghafourian & Roshan (2017) do not include dark matter halos. On the other hand, the dark matter halo is one of the main components of a spiral galaxy in the standard picture. Therefore, without adding a dark matter halo, it is not possible to compare appropriately the galactic dynamics in MOG with those in the standard case. Therefore, we generalize the simulations of Ghafourian & Roshan (2017) by using another N-body code which allows a large number of particles. In fact, in this paper, we will use 2 × 10⁶–10⁷ particles for each live component. On the other hand, we add the dark matter component and make more realistic galactic models.

The outline of this paper is as follows. In Section 2, we briefly introduce MOG and discuss its weak field limit. In Section 3, we describe the numerical method for setting the initial conditions and evolving the point particles. In fact, we make three models that have identical initial conditions, one model in MOG without a dark matter halo and two models in the standard case including a dark matter halo. In Section 4, we compare the bar evolution in these three models and discuss their significant differences and their relevance to current observations. Finally, we discuss future developments.

2. Weak Field Limit of MOG

MOG is a relativistic generalization of general relativity (GR), which, besides the metric tensor, has two scalar fields, i.e., G and μ, and one massive Proca vector field ϕ^α. In other words, MOG is a scalar–tensor–vector theory of gravity. In particle physics language, MOG possesses a spin two graviton, i.e., the metric tensor g_μν; a spin one massive graviton, namely the vector field ϕ^α; and a spin zero massless graviton, i.e., the scalar field G.

The existence of these fields enables the theory to address the dark matter problem without invoking dark matter particles. More specifically, in the weak field limit, MOG explains the flat rotation curve of spiral galaxies (Brownstein & Moffat 2006a), as well as the mass discrepancy in galaxy clusters (Brownstein & Moffat 2006b). This theory has also been applied to the Bullet Cluster and Train Wreck Cluster in order to explain the unusual lensing data; see Israel & Moffat (2016). On the other hand, the cosmological aspects of the theory has been widely investigated; for example, see Moffat & Toth (2009), Moffat (2015a), Roshan (2015), Jamali & Roshan (2016), Shojai et al. (2017), and Jamali et al. (2017). It is important to mention that MOG has a viable sequence of cosmological epochs. On the other hand, in its original form, it cannot be considered as a dark energy model. In other words, this theory is designed to address the dark matter problem and uses a cosmological constant Λ to be consistent with dark energy observations.

In order to study the dynamics of a disk galaxy, the weak field limit of the theory is required. The weak field limit of MOG has been studied in Moffat & Rahvar (2013) and Roshan & Abbassi (2014). In the following, we briefly review the weak field limit. The field equations of the theory, which are obtained by varying the generic action of the theory with respect to the fields, i.e., the metric tensor g_μν, G, ϕ_α, and μ, are given by

$\begin{eqnarray}&&{G}_{\mu \nu }+{g}_{\mu \nu }\left(\displaystyle \frac{2}{{G}^{2}}{{\rm{\nabla }}}_{\alpha }G{{\rm{\nabla }}}^{\alpha }G-\displaystyle \frac{\square G}{G}+G\omega {\mu }^{2}{\phi }_{\alpha }{\phi }^{\alpha }\right.\\ &&\quad -\,\left.4\pi \displaystyle \frac{{{\rm{\nabla }}}_{\alpha }\mu {{\rm{\nabla }}}^{\alpha }\mu }{{\mu }^{2}}\right)-\displaystyle \frac{2{{\rm{\nabla }}}_{\mu }G{{\rm{\nabla }}}_{\nu }G}{{G}^{2}}+\displaystyle \frac{{{\rm{\nabla }}}_{\mu }{{\rm{\nabla }}}_{\nu }G}{G}\\ &&\quad -\,8\pi \left(\displaystyle \frac{1}{4\pi }G\omega {\mu }^{2}{\phi }_{\alpha }{\phi }^{\alpha }-\displaystyle \frac{{{\rm{\nabla }}}_{\mu }G{{\rm{\nabla }}}_{\nu }G}{{G}^{2}}-\displaystyle \frac{{{\rm{\nabla }}}_{\mu }\mu {{\rm{\nabla }}}_{\nu }\mu }{{\mu }^{2}}\right)\\ &&\quad +\,\displaystyle \frac{G\omega }{4\pi }\left({B}_{\,\mu }^{\alpha }{B}_{\nu \alpha }+\displaystyle \frac{1}{4}{B}^{\alpha \beta }{B}_{\alpha \beta }{g}_{\mu \nu }\right)=-8\pi {{GT}}_{\mu \nu },\end{eqnarray} \tag{ 1 }$

$\begin{eqnarray}&&\square G+\displaystyle \frac{G}{2}\left(\displaystyle \frac{{{\rm{\nabla }}}^{\alpha }\mu {{\rm{\nabla }}}_{\alpha }\mu }{{\mu }^{2}}+\displaystyle \frac{R}{8\pi }\right)-\displaystyle \frac{3}{2G}{{\rm{\nabla }}}^{\alpha }G{{\rm{\nabla }}}_{\alpha }G=0,\end{eqnarray} \tag{ 2 }$

$\begin{eqnarray}&&\omega {{\rm{\nabla }}}_{\mu }{B}^{\mu \nu }+\omega {\mu }^{2}{\phi }^{\nu }=4\pi {J}^{\alpha },\end{eqnarray} \tag{ 3 }$

$\begin{eqnarray}&&\square \mu -\displaystyle \frac{{{\rm{\nabla }}}_{\alpha }\mu {{\rm{\nabla }}}^{\alpha }\mu }{\mu }-\displaystyle \frac{{{\rm{\nabla }}}_{\alpha }\mu {{\rm{\nabla }}}^{\alpha }G}{G}+\displaystyle \frac{1}{4\pi }G\omega {\mu }^{3}{\phi }^{\alpha }{\phi }_{\alpha }=0,\end{eqnarray} \tag{ 4 }$

where $\square ={{\rm{\nabla }}}_{\mu }{{\rm{\nabla }}}^{\mu }$ , B_μν = ∇_μϕ_ν–∇_νϕ_μ, G_μν is the Einstein tensor, R is the Ricci scalar, T_μν is the energy–momentum tensor, ω and κ are coupling constants, and J^α is a "matter current" obtained from the variation of the matter action with respect to the vector field and its time component is written as J⁰ = κρ. In this section only, we use units in which the speed of light is unity. In order to find the weak field limit, let us perturb the Minkowski metric η_μν as follows

$\begin{eqnarray}{g}_{\mu \nu }={\eta }_{\mu \nu }+{h}_{\mu \nu }=\left(\begin{array}{cc}1+2{\rm{\Phi }} & 0\\ 0 & (-1+2{\rm{\Psi }}){\delta }_{{ij}}\end{array}\right),\end{eqnarray} \tag{ 5 }$

where $| {h}_{\mu \nu }| \ll 1$ . On the other hand, using the field equations, it is easy to show that the background values of J^α and ϕ^α are zero. Therefore, the linear perturbation of other fields and the matter current can be written as

$\begin{eqnarray} & & {\phi }^{\mu }={\phi }_{1}^{\mu }({\boldsymbol{r}},t),\,\,\,\,\,G={G}_{0}+{G}_{1}({\boldsymbol{r}},t)\\ & & \mu ={\mu }_{0}+{\mu }_{1}({\boldsymbol{r}},t)\,\,\,\,\,\,\,{J}^{\mu }={J}_{1}^{\mu }({\boldsymbol{r}},t),\end{eqnarray} \tag{ 6 }$

where the subscripts "0" and "1" stand for the background and perturbed quantities, respectively. Now, assuming a perfect fluid energy–momentum tensor and substituting Equations (5) and (6) into the field Equations (1)–(4), we find the linearized field equations

$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}{\rm{\Psi }}=4\pi {G}_{0}\left(\displaystyle \frac{16\pi -2}{16\pi -3}\right)\rho ,\end{eqnarray} \tag{ 7 }$

$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}{\rm{\Phi }}=4\pi {G}_{0}\left(\displaystyle \frac{16\pi }{16\pi -3}\right)\rho ,\end{eqnarray} \tag{ 8 }$

$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}{\phi }_{1}^{\mu }-{\mu }_{0}{\phi }_{1}^{\mu }=-\displaystyle \frac{4\pi }{\omega }{J}_{1}^{\mu },\end{eqnarray} \tag{ 9 }$

$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}{G}_{1}=\displaystyle \frac{8\pi {G}_{0}^{2}}{16\pi -3}\rho ,\end{eqnarray} \tag{ 10 }$

$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}{\mu }_{1}=0.\end{eqnarray} \tag{ 11 }$

Now, in order to find the modified version of the Poisson equation, we need to linearize the equations of motion of a test particle, or equivalently, the geodesic equation, in MOG. The geodesic equation in MOG differs from GR and is given by

$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}{x}^{i}}{d{\tau }^{2}}+{{\rm{\Gamma }}}_{\alpha \beta }^{i}\displaystyle \frac{{{dx}}^{\alpha }}{d\tau }\displaystyle \frac{{{dx}}^{\beta }}{d\tau }=\kappa {B}_{\,\alpha }^{i}\displaystyle \frac{{{dx}}^{\alpha }}{d\tau },\end{eqnarray} \tag{ 12 }$

in which the Γ's are the Christoffel symbols. It is easy to show that in the weak field limit, this equation takes the following form,

$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}{\boldsymbol{r}}}{{{dt}}^{2}}\simeq -{\rm{\nabla }}{{\rm{\Phi }}}_{\mathrm{eff}},\end{eqnarray} \tag{ 13 }$

where the effective potential Φ_eff is defined as

$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{eff}}={\rm{\Phi }}+\kappa {\phi }_{1}^{0}.\end{eqnarray} \tag{ 14 }$

On the other hand, Equations (8) and (9) can be straightforwardly solved as

$\begin{eqnarray}&&{\rm{\Phi }}=-\displaystyle \frac{16\pi {G}_{0}}{16\pi -3}\int \displaystyle \frac{\rho ({\boldsymbol{r}}^{\prime} )}{| {\boldsymbol{r}}-{\boldsymbol{r}}^{\prime} | }{d}^{3}x^{\prime} ,\end{eqnarray} \tag{ 15 }$

$\begin{eqnarray}&&{\phi }_{1}^{0}=\displaystyle \frac{\kappa }{\omega }\int \displaystyle \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{ 16 }$

It is clear from Yukawa Equation (9) that the scalar field μ appears as a mass for the vector field. Using these solutions, we rewrite Equation (14) as

$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{eff}}({\boldsymbol{r}})=-{G}_{N}\int \displaystyle \frac{\rho ({\boldsymbol{r}}^{\prime} )}{| {\boldsymbol{r}}-{\boldsymbol{r}}^{\prime} | }(1+\alpha -\alpha {e}^{-{\mu }_{0}| {\boldsymbol{r}}-{\boldsymbol{r}}^{\prime} | }){d}^{3}x^{\prime} ,\end{eqnarray} \tag{ 17 }$

where α is a dimensionless parameter defined as $\alpha =\tfrac{\kappa }{\omega {G}_{N}}$ , and G_N is Newton's gravitational constant. Also, in order to recover the Newtonian gravity in small distances, it is necessary to set 16πG₀/(16π − 3) − κ²/ω = G_N. Now, using Equations (8), (9), (14) and (17), one can readily verify that the Poisson equation in Newtonian gravity is replaced by

$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}{{\rm{\Phi }}}_{\mathrm{eff}}({\boldsymbol{r}})=4\pi {G}_{N}\rho +\alpha {\mu }_{0}^{2}{G}_{N}\int \displaystyle \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{ 18 }$

Hereafter, for simplicity, we drop the subscript "N " from Newton's constant G_N. It is clear that we have two free parameters, α and μ, which should be fixed using relevant dark matter observations. From the rotation curve data of spiral galaxies, we have $\alpha =8.89\,\pm \,0.34$ and μ₀ = 0.042 ± 0.004 kpc⁻¹ (Moffat & Rahvar 2013).

Now, let us find the gravitational potential of a test particle M located at r' = 0. In this case, the density is given by ρ(r') = M δ(r'), where δ is the Dirac delta function. By substituting this density into Equation (17), we arrive at

$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{eff}}(r)=-\displaystyle \frac{{GM}}{r}[1+\alpha -\alpha {e}^{-{\mu }_{0}r}].\end{eqnarray} \tag{ 19 }$

Consequently, the gravitational force experienced by another unit mass located at a distance r from the point mass M is given by

$\begin{eqnarray}&&{\boldsymbol{F}}=-\displaystyle \frac{{GM}}{{r}^{2}}[1+\alpha -\alpha (1+{\mu }_{0}r){e}^{-{\mu }_{0}r}]\displaystyle \frac{{\boldsymbol{r}}}{r}.\end{eqnarray} \tag{ 20 }$

At small radii, the Newtonian force is recovered. On the other hand, at larger scales, the Yukawa term is suppressed, and the gravitational force in MOG is enhanced by a factor (1 + α) compared to the Newtonian force.

As in other alternative theories of gravity for dark matter, the gravitational force is stronger than the Newtonian force. It is interesting to mention that in the strong field limit, MOG predicts new effects that are absent in GR. For example, MOG leads to gravito-electrical and gravito-magnetic effects that are directly related to the existence of the vector field; for more details, see Armengol & Romero (2017) and Moffat (2015b). However, we do not expect this effect in the galactic scale where almost all of the dynamics can be explained in the Newtonian limit.

3. Numerical Method

We use the N-body code developed by Sellwood (2014a). For galaxy models with a live halo, we use the hybrid N-body code described by Sellwood (2003, Appendix B) in which two grids are used to compute the gravitational field: a 3D cylindrical polar grid for the disk component and a spherical grid for the live halo component. One may find the full description of the code in the online manual (Sellwood 2014a).

For galaxy models in MOG, we have modified the code in order to take into account the modified gravity corrections to the gravitational force. In our galaxy model in MOG, as expected, there is no dark matter halo. The energy conservation has been tested, and the error in the total energy is always smaller than 1%.

3.1. Initial Conditions

All of the models we present in this paper start with an exponential disk given by the following surface density,

$\begin{eqnarray}&&{\rm{\Sigma }}(R)=\displaystyle \frac{{M}_{d}}{2\pi {R}_{d}^{2}}\,{e}^{-R/{R}_{d}},\end{eqnarray} \tag{ 21 }$

where M_d is the disk mass and R_d is its length scale. The extent of the disk is limited by tapering the surface density with a function that varies as a cubic polynomial from unity at R = 3.2 R_d to zero at R = 4 R_d. Using the epicyclic approximation, in which particles move on nearly circular orbits, we set the radial velocity dispersion σ_R in such a way to ensure the local stability of the disk. In fact, the local stability criterion is written as Q > 1, where Toomre's stability parameter Q (Toomre 1964) is given by

$\begin{eqnarray}&&Q(R)=\displaystyle \frac{{\sigma }_{R}}{{\sigma }_{R,{\rm{c}}{\rm{r}}{\rm{i}}{\rm{t}}}},{\rm{w}}{\rm{h}}{\rm{e}}{\rm{r}}{\rm{e}}\,\,\,\,{\sigma }_{R,{\rm{c}}{\rm{r}}{\rm{i}}{\rm{t}}}=\displaystyle \frac{3.36G{\rm{\Sigma }}}{\kappa },\end{eqnarray} \tag{ 22 }$

and κ(R) is the epicyclic frequency. Using the epicycle relation, we also find the initial azimuthal velocity dispersion σ_ϕ =κσ_R/2Ω, where Ω(R) = v_c/R, and v_c is the circular velocity. On the other hand, the disk has a Gaussian density profile in the vertical direction with the scale height z₀ = 0.05R_d. Furthermore, the vertical velocities are determined by integrating the vertical one-dimensional Jeans equation as

$\begin{eqnarray}&&{\sigma }_{z}^{2}(R,z)=\displaystyle \frac{1}{\rho (R,z)}{\int }_{z}^{\infty }\rho (R,z^{\prime} )\displaystyle \frac{\partial {\rm{\Phi }}}{\partial z^{\prime} }{dz}^{\prime} .\end{eqnarray} \tag{ 23 }$

For the halo component, we have chosen the isotropic Plummer sphere, which has the following density profile

$\begin{eqnarray}&&\rho (r)=\displaystyle \frac{3{M}_{h}}{4\pi {b}^{3}}{\left[1+{\left(\displaystyle \frac{r}{b}\right)}^{2}\right]}^{-5/2},\end{eqnarray} \tag{ 24 }$

where M_h is the halo mass and b is the characteristic radius. The values of these parameters are selected so that the Newtonian model has the same rotation curve as the disk in modified gravity (which has no spherical component). The halo is truncated at r = 3 b. It is important to mention that this mass distribution is consistent with the rising rotation curves in late-type galaxies of lower luminosity (e.g., de Blok et al. 2008).

It is necessary to mention that we have two main reasons for choosing this halo. In fact, with this choice, it is easy to match the rotation curves of the dark matter and modified gravity models. On the other hand, this profile has also been used in Combes (2014) in order to compare bulge formation in MOND and in standard Newtonian dynamics. Also, it has been used in Tiret & Combes (2007) to study the stellar bar evolution in MOND. Therefore, this choice of halo enables us to make a direct comparison with Tiret & Combes (2007), and consequently, a direct comparison between MOND and MOG.

Dejonghe (1987) gives expressions for the distribution function (DF) of the Plummer sphere, characterized by a parameter q. We set q = 0, which corresponds to an isotropic sphere. Therefore, we start with a sphere with a known DF when it is in isolation. However, when a disk component is added, one needs to find the equilibrium DF of the spherical component in the composite halo and disk model. To do so, we use a halo compression algorithm described by Young (1980) and Sellwood & McGaugh (2005). This algorithm uses the fact that both radial action and the total angular momentum are conserved during compression. Finally, the particles will be selected from the new DF function as described in Debattista & Sellwood (2000). For more details, see Sellwood & McGaugh (2005) in which the iterative solution for a single spheroid with an added disk was studied using the above-mentioned procedure.

We emphasize that no bulge is included in our models. More specifically, we study two single-component models, "EN" and "EM," with an exponential disk. EN is the exponential disk model in Newtonian gravity, and EM is the corresponding model in MOG. On the other hand, "EPR" and "EPL" are two-component models including an exponential disk and a Plummer halo. In EPR, the halo is rigid, and in EPL, it is live and responsive. The properties of these models are summarized in Table 1. For each responsive component, we employ two million particles. For example, in EPL, we employed two million particles for the disk and two million particles for the live halo. Other models have two million particles. In order to check that the results do not depend on the number of particles, we also increase the number of particles by a factor of 10 in all of the models.

Table 1. Model Information

	Disk	Disk	Disk	Initial	Halo	Halo	Halo	Halo	Time	Cylindrical		Spherical	Softening
Run	Type	Edge	Mass	Q	Type	Mass	Scale	${r}_{\max }$	Step	Polar Grid	δz	Grid	Length
EN	Exp	4	1	1.5	none	⋯	⋯	⋯	0.01	390 × 407 × 135	0.08	⋯	0.16
EM	Exp	4	1	1.5	none	⋯	⋯	⋯	0.01	193 × 224 × 125	0.08	⋯	0.16
EPR	Exp	4	1	1.5	Plum	7.5	5.5	⋯	0.01	99 × 128 × 125	0.08	⋯	0.16
EPL	Exp	4	1	1.5	Plum	8	10	30	0.01	100 × 128 × 125	0.08	1001	0.16

Note. Column 1: letter identification for the simulation. Column 2: type of initial disk; "Exp" is exponential disk. Column 3: initial outer radius of the disk in units of R_d. Column 4: initial disk mass M_d. Column 5: initial Toomre's Q parameter. Column 6: type of halo, "Plum" is Plummer isotropic sphere. Column 7: mass of the halo component, M_b. Column 8: radial scale of the halo component b in units of R_d. Column 9: outer edge of the halo, ${r}_{\max }$ , in units of R_d. Column 10: basic time step in units of τ₀. Column 11: number of rings, spokes, and planes in the cylindrical polar grid. Column 12: vertical distance between grid planes. Column 13: number of shells in the spherical grid. This grid extends to R = 180R_d. Column 14: gravity softening length in units of R_d.

Download table as: ASCII Typeset image

3.2. Units

We use M_d as the unit of mass and R_d as the unit of length throughout this paper. In other words, we use units such that R_d = M_d = 1. Also, we scale the Newtonian gravitational constant as G = 1. In this case, the unit of velocity is V₀ = (GM_d/R_d)^1/2, and the time unit is τ₀ = R_d/V₀ =(R_d³/GM)^1/2. For example, a suitable choice for our models can be R_d = 2.6 kpc and τ₀ = 10 Myr, which yields M_d ≃4 × 10¹⁰M_⊙ and V₀≃254 km s⁻¹. For all models, we compute the evolution for 500τ₀. We use a time step Δt = 0.01τ₀, and the time step is increased by successive factors of 2 in three radial zones. The cubic spline softening length (Monaghan 1992) is set to 0.16 R_d in all models. Also, except in the model EPR, which has a rigid component, the grid is recentered every 16 time steps. This dynamic recentering helps avoid numerical artifacts.

3.3. Rotation Curves

Using the above-mentioned units, the rotation curves of our four models have been illustrated in Figure 1. It is necessary to mention that for an exponential disk, the rotation curve in MOG is computed for a given α and μ₀ first, i.e., the model EM, and then equivalent Newtonian systems are built with a dark matter halo added, in order to obtain the same rotation curve. More specifically, by changing the halo parameters, i.e., the mass and the characteristic length, in both the rigid and live cases, we construct the models EPR and EPL, whose rotation curve is similar to that of the model EM. It should be emphasized that in order to compare the stellar bar evolution in MOG and Newtonian gravity, it is necessary to start with the same initial rotation curves. In this case, the dispersion velocities are also almost the same. One should note that in model EM, we have set α = 8.9 and μ₀ = 0.015R_d⁻¹, where α is close to the observed value. However, we are constructing toy galactic models and naturally we do not have to use the observational values. On the other hand, we know that these parameters change from galaxy to galaxy.

3.4. Bar Amplitude

As we mentioned before, the main aim of this paper is to study and compare the bar growth in Newtonian gravity and MOG. Therefore, the main quantity to be measured in our simulations is the bar amplitude. It is important to mention that the code determines the gravitational forces on the cylindrical polar grid using sectoral harmonics $0\leqslant m\leqslant 8$ and surface harmonics 0 ≤ l ≤ 4 on the spherical grid. For more details, see the online manual available at http://www.physics.rutgers.edu/~sellwood/manual.pdf.

The amplitude of the non-axisymmetric features can be measured by computing

$\begin{eqnarray}&&{A}_{m}(t)=\displaystyle \sum _{j}{\mu }_{j}{e}^{{im}{\phi }_{j}(t)},\end{eqnarray} \tag{ 25 }$

where μ_j is the mass of the jth particle particle, and ϕ_j is its cylindrical polar angle at time t. One should note that this summation is over the disk particles only. So, the bar amplitude, for which m = 2, would be obtained from the ratio A₂/A₀.

4. Results

4.1. Stabilizing Effects of MOG

In the recent paper Ghafourian & Roshan (2017), we compared two galaxy models which are similar to EM and EN. In fact, without adding a spherical halo, it has been shown that MOG has stabilizing effects. More specifically, increasing the MOG free parameters supports the global stability of the exponential disk. This point has also been checked for the Mestel disk. In fact, the procedure in Ghafourian & Roshan (2017) is similar to that presented in Ostriker & Peebles (1973), in which it has been proven that a rigid halo has stabilizing effects. The maximum number of particles in Ghafourian & Roshan (2017) is 2 × 10⁴. Therefore, the results may suffer from numerical artifacts caused by employing a small number of particles. Here we have increased the particle number to 2 × 10⁷ and used a code that employs a completely different technique to compute the gravitational field. In model EN, the disk rapidly extends to larger radii. Therefore, we have to enlarge the grid in order to avoid particles from escaping from the grid. In the following, let us briefly show that the main results of Ghafourian & Roshan (2017) do not change by increasing the number of particles. Then, we will discuss more realistic models by adding dark matter halos to the Newtonian models.

In the model EM, we have chosen $\alpha =8.9$ and μ₀ =0.015R_d⁻¹. The corresponding rotation curve is shown in Figure 1. As expected, the circular velocity at large distances from the center of the disk is higher than the Newtonian model EN. The evolution of the bar amplitude is shown in Figure 2. In both top and bottom panels in this figure, the black dashed curve belongs to the model EN and the red solid curve belongs to the model EM. We have used a logarithmic scale so that the growth rate can be estimated from the slope of the curve in the period of linear growth. The small boxes inside each panel show the evolution of the bar amplitude at early times τ < 100. It is clear from the top panel that the exponential instability growth starts t around τ ≃ 25 in both models EN and EM. However, in the model EM, the growth rate is smaller than that in the Newtonian case. This point is clearer from the bottom panel, where we have employed a larger number of particles. In both cases, after reaching a maximum, the bar decays slowly at late times. The slope of this decrease is smaller in the bottom panel. However, the qualitative behavior is almost the same. It is clear from the top panel that the final magnitude of the bar in EM is smaller than that in EN. In the bottom panel, this will happen at later times. This fact has also been reported in Ghafourian & Roshan (2017).

**Figure 2.** The top panel shows the evolution of the bar amplitude for our models when each responsive component has 2 × 10⁶ particles. The bottom panel shows the bar amplitude when there are 2 × 10⁷ particles in each responsive component. For example, the total number of particles in the model EPL is 4 × 10⁷. The small boxes in each panel show the evolution of the bar at early times, i.e., τ < 100.
Download figure:
Standard image High-resolution image

The bar pattern speed has been shown in Figure 3. Similar to the results presented in Ghafourian & Roshan (2017), the pattern speed in the model EM is larger than that in the EN model.

4.2. Bar Growth

As mentioned before, it is necessary to compare MOG with Newtonian models that include a dark matter halo. In this section, we compare the model EM with the two-component models EPR and EPL. EPR has a rigid Plummer halo and EPL includes a live Plummer halo. In fact, we have chosen the halo's physical properties such that the rotation curve coincides with that of the model EM. As is clear from Table 1, the mass and the characteristic length scale in the rigid halo are different from those of the live halo. In the case of live halo, if we use the same parameters used in the rigid halo, the rotation curves will not match. In fact, in model EPL, the halo is compressed by the growth of the initial disk and leads to a different rotation curve than that in model EPR. For example, see models "GR" and "GL" in Berrier & Sellwood (2016).

It is clear from Figure 2 that the bar grows rapidly in model ERL. On the other hand, the growth rate in model EPR is substantially lower than that in both EM and EPL models. Also, there are rapid oscillations in the bra amplitude of the EPR model. However, it does not reveal a physical process. In fact, for models with rigid components, we do not use moving grids. Therefore,the center of mass of the system rotates around the grid center and gives rise to these unreal oscillations in the bar growth.

It is a well-known fact that disks in live halos form bars more readily than those in rigid halos (Athanassoula 2002). It has been shown in Sellwood (2016) that the increased growth rate of the bar instability in systems with a live halo results from the angular momentum exchange between the halo and the disk. We see from the top panel of Figure 2 that the exponential growth in both models EPL and EM starts at around τ ≃ 25 and reaches its maximum at τ ≃ 100. However, the instability growth rate is smaller in the model EM. More specifically, we found that the exponential growth rate, i.e., e^ωτ, in the EPL and EM models are ω ≃ 0.067 and ω ≃ 0.054, respectively.

It is important to mention that there is a significant difference between the bar growth in MOG and MOND. It is shown in Tiret & Combes (2007) that the growth rate in MOND is substantially higher than in the dark matter halo case. It is not the case in MOG, and MOG postpones the occurrence of the bar instability. Although the halo density in Tiret & Combes (2007) is the same as ours, the disk density in Tiret & Combes (2007) is different and given by a Miyamoto–Nagai disk. Therefore, for a more careful comparison between MOG and MOND, it is necessary to construct equivalent galactic models with the same surface densities. We leave this work to a future study.

On the other hand, it is interesting to mention that the maximum value of the bar magnitude in the model EM is smaller than that in the model EPL. This means that MOG leads to weaker bars during bar instability. This behavior is completely contrary to MOND, which leads to stronger bars during the instability. It is clear in Figure 2 that the final bar magnitude in MOG is smaller than that in EPL. This is also the case in MOND. In other words, although the maximum magnitude of the bar is larger in MOND compared with the dark matter halo model, its final value is smaller (Tiret & Combes 2007).

The projected positions of the particles on the x–y, x–z, and y–z planes for the EM and EPL models are shown in Figure 4. The two left columns belong to the model EM, and the two right columns belong to the model EPR. In both models, a strong two-arm spiral develops at τ ≃ 100 and rapidly disappears and gives way to a bar. The bar in model EPL starts to grow slowly. However, in the EM model, the bar magnitude decreases and oscillates at the same time.

This oscillation can be related to the existence of some beating between different modes. This oscillatory behavior in model EM can also be seen in the pattern speed evolution. In fact, one should note that A₂ is not only a bar strength but also a spiral strength. In principle, these pattern speeds of waves can be slightly different, and this can give rise to the above-mentioned oscillations.

In Figure 5, we have plotted the relative overdensities associated with the m = 2 mode at different times. The four left columns belong to the model EPL, and the four right columns to the model EM. In the model EPL, after τ ≃ 250, the bar mode retains its dominance. The twofold symmetry is clear until the end of the simulations. On the other hand, in the EM model, the evolution is completely different. As seen in the right panel of Figure 5, the twofold symmetry fades and appears frequently. This is completely in agreement with the period of oscillations observed in Figure 2.

4.3. Power Spectra

The clearest way to support the existence of beating between different modes in the model EM is to find the power spectrum. As we already mentioned, we compute the azimuthal Fourier coefficients of the mass distribution at each grid radius and save it at regular time intervals. In order to find the power in terms of frequency and radius for each sectoral harmonic m, we compute the Fourier transformation in time of the coefficients.

In Figure 6, we have shown the power spectra for three sectoral harmonics m = 2, 3, and 4 for the model EM with 2 × 10⁷ particles. The top row is for the first part of the evolution, 22 < τ < 260, and the bottom row is for the second part of the evolution. In each panel, the solid (red) curve shows mΩ_c, and the dashed (blue) curves mΩ_c ± κ, where Ω_c(R) is the angular frequency for the circular motion and κ(R) is the Linblad epicycle frequency. Since each horizontal ridge indicates a coherent density wave, it is clear from the bottom-left panel that there are two m = 2 waves in the second half of the simulation with pattern speeds mΩ₁ ≃ 0.52 and mΩ₂ ≃ 0.3. In this case, the beat angular frequency is Ω_b = m(Ω₁ − Ω₂) = 0.22. Consequently, the beat period is τ_b ≃ 28.6, which is completely in agreement with the period of oscillations in the bar magnitude shown in the bottom panel of Figure 2. Therefore, the power spectra confirm the existence of beating between modes.

By plotting the magnitude of the m = 3 and m = 4 modes, it turns out that they also have oscillatory behavior. Similar to the twofold symmetric mode m = 2, it is clear from the bottom row of Figure 6 that beating also occurs for m = 3 and m = 4 modes. For these modes, there are at least four different beating waves with different frequencies.

4.4. Vertical Structure

Particles can resonate both in the plane of the disk and perpendicular to it. In this case, the vertical oscillations of the particles can be amplified, and a peanut shape forms (Combes & Sanders 1981). One may infer from the top panel of Figure 2 that the bar strength drops substantially in both models EPL and EM at τ ≃ 150. On the other hand, at this time, the vertical thickness of the disk grows significantly. To see this more clearly, we have plotted the root mean square (rms) thickness, i.e., $\sqrt{\langle {z}^{2}\rangle }$ , in Figure 7 for all models for the same instants presented in Figure 4. Therefore, the drop in the bar amplitude coincides with the peanut resonance in both models. The vertical resonances thicken the disk and consequently stabilize the disk against bar formation. In this case, it is natural to expect that the bar weakens. On the other hand, in the EPR model, there is no drop in the evolution of the bar strength, and consequently, the thickness does not grow, and the peanut shape does not appear. This confirms that the drop in the bar magnitude is due to the peanut occurrence.

It is seen in Figure 7 that the position of the peanut lobe in both models EM and EPL is in the interval 1.5 R_d < R < 2 R_d and does not shift significantly during the disk evolution. Therefore, let us measure the rms thickness at R = 1.5 R_d with respect to time. The result has been shown in Figure 8. The bar instability heats the disk and thickens it continuously. The steps correspond to the occurrence of the peanut, when the particles suddenly leave the plane of the disk. The buckling occurs around τ ≃ 80 in model EPL, and τ ≃ 100 in the EM model. This is consistent with the above-mentioned claim that the weakening of the bar coincides with the appearance of the peanut.

**Figure 8.** Evolution of the rms thickness of the disk with respect to time. The bar instability continuously heats the disk and consequently thickens it. The obvious steps correspond to the occurrence of peanut resonance.
Download figure:
Standard image High-resolution image

In both models, buckling occurs in a short time interval compared to the simulation time. In model EPL, during the buckling, the rms thickness increases from 0.05 R_d to 0.3 R_d. On the other hand, in the model EM, the increase rate is slower, and the rms thickness increases from 0.05 R_d to 0.2 R_d. Therefore, this simulation shows that the dark matter stellar model produces stronger peanuts compared with the model in MOG. This is also the case in stellar models in MOND (Tiret & Combes 2007).

4.5. Pattern Speed Ω_p

It is instructive to calculate the pattern speed Ω_p(t). The result has been illustrated in Figure 3. There are 2 × 10⁶ particles in the each live component in the top panel. On the other hand, for comparison, there are 3 × 10⁶ particles in the bottom panel. As reported in Ghafourian & Roshan (2017), the pattern speed in model EM is larger than that in the EN model. On the other hand, Ω_p(t) oscillates in model EM. This is not the case in the other models.

There are some important features in Figure 3. The pattern speed starts with almost the same value for all models. However, although all models begin with the same initial conditions, the pattern speed's magnitude and evolution are different when the disks evolve. More specifically, Ω_p(t) oscillates in EM but its mean magnitude is almost constant. However, in the model EPL, when the bar is formed at τ ≃ 100, the pattern speed starts to decrease at τ ≃ 125. In fact the particles of the halo slow down the bar through dynamical friction. It is clear that there is no such decrease in the model EPR, where the halo is rigid. So, one can be assured that the substantial decrease in the model EPL is due to the dynamical friction induced by the live halo particles. As expected, dynamical friction does not appear in the MOG model, i.e., in model EM. This is also the case in MOND; see Tiret & Combes (2007) for more details. The absence of dynamical friction in modified gravity can be a crucial point in distinguishing between dark matter halo and modified gravity from an observational point of view. For example, the inefficiency of dynamical friction can substantially reduce the merger frequency of the galaxies. As a consequence, bulge formation and its morphology and kinematics, in principle, can be different from the standard view (Combes 2016).

Let us emphasize an important observational fact that is relevant to the pattern speed of disk galaxies. Observations show that independent of the Hubble type of disk galaxies, bars have been formed and then evolve with time as fast rotators; see Aguerri et al. (2015) and references therein. On the other hand, these observations are not in complete agreement with the numerical simulations. For example, a recent ΛCDM cosmological hydrodynamical simulation by Algorry et al. (2017) verifies that bars slow down as they grow. Therefore, the existence of fast bars can be considered a challenge to ΛCDM models (Debattista & Sellwood 2000).

In order to compare the pattern speeds in models EM and EPL, let us use the parameter ${ \mathcal R }$ defined as

$\begin{eqnarray}&&{ \mathcal R }=\displaystyle \frac{{D}_{L}}{{a}_{B}},\end{eqnarray} \tag{ 26 }$

where D_L is the corotation radius and a_B is the bar's semimajor axis. Bars are fast if ${ \mathcal R }\simeq 1$ , and slow if ${ \mathcal R }\gg 1$ . Finding the corotation radius in the simulation is not hard in the sense that we have the bar pattern speed and the angular frequency Ω(R) with respect to time. Therefore, one can simply estimate the corotation radius. On the other hand, in order to estimate the bar's semimajor axis at a given time τ, we assume a rectangle with width ΔL and length 5 around the line y = tan ϕ(τ)x, where ϕ(τ) is the angular displacement of the bar. Then, we divide this rectangle into small elements and choose the bar length to be the length at which the density of the particles is less than b per cent of the central element. We vary ΔL between 0.2 and 1, and choose two values for b, i.e., 10 and 20. Naturally, these choices lead to an error bar in the final value of ${ \mathcal R }$ each time.

The result has been shown in Figure 9 for both models EM and EPL. In this figure, we computed ${ \mathcal R }$ for models with 2 × 10⁷ particles in each live component, at four different times. In fact, we realize that the corotation radius in model EPL increases with time. More specifically, at τ = 500, the corotation radius is D_L ≃ 3.6, while at τ = 800, it reaches to D_L ≃ 5.3. On the other hand, D_L stays constant in model EM around 2.7. Furthermore, the bar length for both models is almost constant for τ > 500. Therefore, we expect that ${ \mathcal R }$ increases with time in model EPL and stays constant in model EM. This behavior can be seen in Figure 9. It is important to mention that the bar is faster in model EM compared to the dark matter halo model. In fact, in reality, most observed bars appear to lie in the range $0.9\lesssim { \mathcal R }\lesssim 1.4$ , and consequently, the observed bars are fast. From this perspective, it is clear from Figure 9 that the model EM is closer to the observed interval, while both models lie outside the interval.

**Figure 9.** The ratio ${ \mathcal R }={D}_{L}/{a}_{B}$ for models EM and EPL with 2 × 10⁷ particles in each live component, calculated at four different times. The error bars are calculated from the averages introduced for measuring the bar length. The dashed line indicates ${ \mathcal R }=1.4$ . It is clear that bars are faster in the modified gravity model than in the standard dark matter halo case.
Download figure:
Standard image High-resolution image

However, more careful studies are needed to compare MOG with the relevant observations. More specifically, cosmological hydrodynamical simulations would be extremely helpful in performing a statistical study for comparing MOG and observations. To the best of our knowledge, cosmological simulations in the context of alternative theories, which dispense with the need for dark matter particles, have not been developed.

5. Discussion and Conclusion

In this paper, we have presented N-body simulations of isolated galaxies in both approaches, in Newtonian gravity with a dark matter halo and in MOG. The main aim is to compare bar growth in the presence and absence of a dark matter halo. More specifically, in the absence of a dark matter halo, we used the modified gravitational force introduced in MOG. We construct three galactic models with almost the same initial conditions, i.e., EM, EPR, and EPL. The model EM belongs to MOG and has no spherical component. On the other hand, the model EPR is a galactic model with an exponential disk and a rigid Plummer spherical component. The model EPL is more realistic, and its halo is a live Plummer sphere. In fact, we have a fourth toy model named EN. This model is an exponential disk in Newtonian gravity without a dark matter halo. Comparing this model with EM, we recovered the main results already reported in Ghafourian & Roshan (2017) using low-resolution simulations.

Our results show that the bar growth in model EM is slower than in model EPL, and faster than in the EPR model. In fact, MOG behaves neither as a live dark matter halo nor a rigid dark matter halo. However, when compared with the EN model, it is clear that MOG has stabilizing effects and suppresses the bar instability for a while. Also, we found that the final strength of the bar in MOG is smaller than that in both the EPR and EPL models. Moreover, the maximum magnitude of the bar in MOG is smaller than that in the Newtonian models. As we have already mentioned, this is in contrast with MOND, where stellar simulations lead to stronger bars. However, although the bar's maximum strength is unchanged for a relatively long time in MOND, it retains its maximum value for a short time in MOG.

We also found that the evolution of the pattern speed in these models is different. More specifically, the pattern speed in the model EM oscillates around an almost constant value and does not change substantially. However, the pattern speed in the model EPL decreases significantly due to the dynamical friction induced by the dark matter particles. In other words, not only is the morphology of the bar in our modified gravity theory different from that in the standard cases, but its dynamics is also different. We explicitly showed that the model EM leads to faster bars compared to the dark matter model.

Furthermore, the rms thickness of the disk in the model EM at the end of the simulation is substantially smaller than that in the model EPL. In other words, the peanut formation is more effective in model EPL than in model EM.

More investigations are still required for galactic dynamics in MOG. For example, it would be interesting to study the bar growth in the presence of a massive spherical bulge and a galactic gas component. In fact, the gas component has crucial effects on the system. For example, bars exert torque on the gas component, which inflows toward the center, and consequently the gas angular momentum is given to the bar. In this case, the bar is weakened, and finally, the bar will be faster and shorter (Bournaud & Combes 2002).

Also, the stabilizing effects of MOG should be compared with other well-known dark matter halo profiles. Furthermore, careful study is required to compare the bar evolution in MOG and MOND. In this case, the initial conditions, i.e., the disk density profile and the velocities, should be the same in both theories. Our results here show that MOG has some features that are completely different from MOND. We leave these issues as subjects of study for future works.

To summarize the results, let us make a quick comparison between the bar evolution in MOG and that in dark matter halo models. The main differences are: (1) the bar growth rate is smaller in MOG, (2) the final magnitude of the bar is smaller in MOG, (3) the pattern speed is oscillatory in MOG, and it also does not decrease with time; this is the main difference between MOG and the dark matter models studied in this paper. (4) There is no dynamical friction experienced by the bar in MOG. (5) The maximum value of the bar strength is smaller in MOG. We reiterate that the bar magnitude in MOG stays in its maximum value for a relatively short time interval. In other words, MOG predicts that the strongly barred galaxies are not that frequent today. This fact is in agreement with the relevant observations, which show that less than 20% of barred galaxies are strongly barred (Whyte et al. 2002; Díaz-García et al. 2016).

As our final remark, we stress that in order to make a precise comparison between observations and MOG predictions for the abundance, size, and pattern speeds of bars, as well as for their evolution with redshift, it is necessary to perform high-resolution cosmological simulations. For instance, for a recent work in the standard cosmological model, see Algorry et al. (2017). However, there are no similar works in modified theories of gravity designed to dispense with the need for dark matter particles, even in MOND and its relativistic versions. Therefore, large numerical effort is still required to be able to compare modified gravity predictions for galactic properties with relevant observations.

I would like to appreciate Francoise Combes for insightful and constructive comments. I am particularly grateful to Jerry Sellwood for guidance in running the GALAXY code. The main part of the calculations have been done on the ARGO cluster in the International Center for Theoretical Physics (ICTP). I would like to thank ICTP for very kind hospitality, during which the main part of this work has been done. This work is supported by the Ferdowsi University of Mashhad under grant No. 2/44972 (18/07/1396).

Stellar Bar Evolution in the Absence of Dark Matter Halo

Article metrics

Permissions

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction

2. Weak Field Limit of MOG