Does Modified Gravity Predict Fast Stellar Bars in Spiral Galaxies?

As mentioned earlier, the main purpose of this work is to investigate the model-independency of previously reported results on the difference in disk evolution in modified gravity and DM. Because previous studies focused on exponential (EXP) disks and Plummer halos (Roshan 2018; Roshan & Rahvar 2019), to ensure the independence of the results for the adopted disk/halo models, we first changed the disk's density distribution from EXP to Kuzmin–Tommre (KT) while we keep the live Plummer halo model unchanged (LPH model). In the next step, we change both the halo and disk profiles. In this case, we employ a KT disk and a live isothermal halo (LIH model). To investigate further, we also consider another halo model with the Hernquist profile (LHH model). This halo model has a cusp in the inner radii, which increases the initial rotational velocity in the inner radii and influences the evolution of the disk. As we will discuss comprehensively, the dynamical friction plays a central role in discriminating between DM and modified gravity. On the other hand, it is natural to expect that this mechanism works differently in the cored and cuspy halo models. Therefore, to draw a complete and consistent picture, it seems necessary to compare both the cored and cuspy types with modified gravity models. For this purpose, we include the LHH model in the simulations.

In this study, the KT disk, model 1 of Toomre (1963), has a surface density of the form

$$\begin{eqnarray}&&{\rm{\Sigma }}(R)=\displaystyle \frac{{M}_{d}}{2\pi {R}_{d}^{2}}{\left[{\left(\displaystyle \frac{R}{{R}_{d}}\right)}^{2}+1\right]}^{-3/2},\end{eqnarray} \tag{ 2 }$$

where R is the radial coordinate in the cylindrical coordinate system (R, ϕ, z). The Plummer halo, a polytrope model of index 5, has the density

$$\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{ 3 }$$

while for the simple cored Isothermal halo model (Evans 1993), the density is given by

$$\begin{eqnarray}&&\rho (r)=\displaystyle \frac{{V}_{h}^{2}}{4\pi {{Ga}}^{2}}\displaystyle \frac{3+{x}^{2}}{{\left(1+{x}^{2}\right)}^{2}}.\end{eqnarray} \tag{ 4 }$$

The Hernquist halo model (Hernquist 1990) also has density of the form

$$\begin{eqnarray}&&\rho (r)=\displaystyle \frac{{M}_{h}c}{2\pi r{\left(r+c\right)}^{3}}.\end{eqnarray} \tag{ 5 }$$

In the above equations, M_d and M_h are the disk and halo mass, respectively; R_d is the radial disk scale length; a, b, and c are characteristic radii of the halos; and x = r/a. Furthermore, V_h in the isothermal sphere is ${V}_{h}=\sqrt{{{GM}}_{h}/a}$. Note also that the system of units is considered so that M_d = R_d = G = 1, where G is the gravitational constant. Therefore, the units of velocity and time can be written as ${V}_{0}={({{GM}}_{d}/{R}_{d})}^{1/2}$ and ${\tau }_{0}={R}_{d}/{V}_{0}={({R}_{d}^{3}/{GM})}^{1/2}$. As a suitable example, selecting ${R}_{d}=2.6\,\mathrm{kpc}$ and τ₀ = 10 Myr yields M_d ≃ 4 × 10¹⁰M_⊙ and V₀ ≃ 254 km s⁻¹.

To make a reasonable comparison, we constrain our models to start their evolution from a similar state. Therefore, the models are constructed such that the initial rotational velocities of the four models are as close as possible. Furthermore, the initial velocity dispersions of the models have a very similar trend. To do this, the employed free parameters of the DM halo density and modified gravity should be varied until the same initial state for the baryonic matter is achieved in all models. Notice that we use the same parameters for the disk density in all models.

The truncation radius of the KT disk is chosen to be $r=6{R}_{d}$, where the density starts to taper smoothly from $r=5{R}_{d}$ by using a cubic function. For our Plummer model, the halo mass is chosen as M_h = 3.6M_d and the halo scale length as b = 8.5R_d, where the truncation radius is at r = 3.4b. In the Isothermal model, the halo mass is ${M}_{h}=1.7{M}_{d}$ and its scale length is $a=6{R}_{d}$ with truncation radius $r=10a$. In the Hernquist model, the adopted halo mass is M_h = 14M_d, the scale length is c = 19R_d, and the truncation radius is $r=8c$. As mentioned earlier, the cusp in the Hernquist halo results in a relatively higher rotational velocity in the inner radii, but, with this choice of parameters, the rotation curve in the outer radii matches the other models very well. The LHH model seems too massive compared to the other models. However, regarding Gauss's theorem, at least at the beginning of the simulation, only the DM mass inside a sphere with radius 6R_d, that is, $M(6{R}_{d})$, contributes to the gravitational force. It is easy to show that $M(6{R}_{d})\simeq 0.8$, 0.7, and 0.85 for the LHH, LPH, and LIH models, respectively. In other words, in all of the models, the DM mass fraction inside $r=6{R}_{d}$ is around 0.4. Therefore, one may ensure that all of the models start with suitable initial conditions.

For the NLG model, the values selected for the free parameters are α = 10.9 and $\mu =0.0045{R}_{d}^{-1}$. Also, for the MOG model we employ α = 2.8 and $\mu =0.031{R}_{d}^{-1}$. The resulting rotational velocities of our models are illustrated in Figure 1.

Zoom In Zoom Out Reset image size

For all of the models, the Toomre parameter Q (Toomre 1964) is considered to be 1.5. This value is high enough for our models to protect the disk from local instability and fragmentation.

It should also be mentioned that for every active component (DM halo/disk) in our main simulations, we have adopted N = 2 × 10⁶ particles. However, to make sure that the results are independent of the particle number, we have performed the simulations for a higher number of particles (Section 4.6).

Recall that the adopted mesh for the disk and DM halo is cylindrical polar (P3D) and spherical (S3D), respectively. The grid size selected for our S3D mesh is 1001 × 2501, where the first number is the number of radial grid shells and the second is the radius of the grid outer boundary. For the P3D mesh, we have 193 × 224 × 125, where the numbers demonstrate the number of mesh points in the radial, azimuthal, and vertical directions, respectively. Also, "lscale," which is the number of mesh points on the length unit, is chosen to be 12.5. Furthermore, the softening length in the P3D mesh is  = 0.16R_d. This choice of mesh points is enough to certify that the majority of particles do not leave the mesh during the simulations. This property is checked for all of the models, and less than 2% of the particles escape the mesh before the end of the simulation.

It is necessary to mention that the computation of the potential is achieved by using sectoral harmonics $0\leqslant m\leqslant 8$ in the P3D mesh and surface harmonics $0\leqslant l\leqslant 8$ in the S3D. More details can be found in the online manual.³

The duration of our simulations is selected to be τ = 800τ₀ and the time step is Δτ = 0.01τ₀. It should also be mentioned that three different time zones with factors of 1, 2, and 4 Δτ are selected to account for the evolution of particles in regions of different density; shorter time steps are required for particles in denser regions. However, to ensure that the results are not affected by the choice of time step, we also checked the results for lower values of Δτ (Section 4.6). Moreover, the grid is recentered every 16 steps to avoid numerical artifacts.

According to the selected options, the code computes the initial positions and velocities of the particles to form the initial equilibrium state of the system. Note that if the model consists of more than one component, it is necessary to find the distribution function (DF) in the composite system and then choose the particles from the resulting DF. Accounting for this matter in our DM model, we used the Compress algorithm (Young 1980; Sellwood & McGaugh 2005) implemented in the code, which finds the potential and changes the halo density so that the equilibrium condition in the presence of the disk is met. The initial positions and velocities of the system are then obtained.

In order to investigate the formation and evolution of the bar in our models, we use the Fourier expansion of the mass distribution of the particles in the disk plane. It is straightforward to show that the mth Fourier coefficient is written as

$$\begin{eqnarray}&&{A}_{m}(\tau )=\left|\sum _{j}{\mu }_{j}{e}^{{im}{\phi }_{j}}\right|,\end{eqnarray} \tag{ 6 }$$

where μ_j is the mass and ϕ_j is the cylindrical polar angle of the jth particle at time τ. Note that we have used the same mass for all of the particles in our simulation, so every particle has an equal contribution to the total mass of the system. Also, recall that the above summation is only over the disk's particles and naturally does not take the DM particles into account. Considering the above equation, one can find the bar/spiral amplitude by calculating the ratio of ${A}_{2}/{A}_{0}$. Although m = 2 includes both the spiral and bar-like features, because the spirals are short-lived patterns, it would be reasonable to use this quantity to account for the effect of the bar mode.

Figure 2 shows the evolution of the bar for our models. From this figure, it is apparent that the bar amplitude in both modified gravity theories starts with an initial peak, and then after a drop, it continues to oscillate around an almost constant value. Compared to MOG, the formation of the bar is faster in NLG and the bar is stronger. Also, the bar amplitude reaches its constant value at τ ≃ 150, while in the MOG model it happens later at about τ = 200. In comparison to NLG, the bar formation in DM models starts later and at a lower rate. It also has an initial peak that appears at τ ≃ 75 for the Isothermal (orange) and at τ ≃  100 for the Plummer model, but after a drop at τ ≃ 170, it starts to grow again at an almost constant rate. Both models experience a decreasing phase around τ ≃ 400, and then it continues to grow until the end of the simulation. However, note that although the bar strengths in these two DM models have an overall similar trend and reach almost the same value at final states, the initial behavior in the formation and growth of the bar is not independent of the adopted halo profile. One can infer that, in comparison to the Plummer model, the bar starts to form in the Isothermal model faster and at an earlier time, which matches MOG in early stages. However, a clear discrepancy between the evolution of the disks in the modified gravity and cored DM models is visible after about τ ≃ 200, regardless of the selected mass model.

Zoom In Zoom Out Reset image size

The behavior of the LHH model, which includes a cuspy DM halo at the beginning of the simulation, shows some differences from the other DM models. From Figure 2, it is clear that the growth rate of bar formation in the LHH model is higher. Moreover, the value of bar strength remains around its initial peak for a longer duration in comparison to the other halo models. Then, at about τ ≃ 250, another phase in the evolution of the bar amplitude in the LHH model starts, in which the bar weakens. This period is also long-lasting in comparison to the other DM models. At about τ ≃ 500, the bar amplitude of the LHH model starts to grow again. However, at τ ≃ 800, the bar in LHH is clearly weaker than in the LPH and LIH models. Since the evolution of the LHH model does not seem to be completed at this stage, we continued its evolution up to τ = 1600 in the subplot of Figure 2. As one can see, the growing phase of the bar amplitude in this model continues to the end of the simulation. According to this figure, it could be concluded that although the details in the evolution of a cuspy DM model are different from a cored DM model and it shows some delays, its overall behavior does not change and is still very different from the models of modified gravity.

According to this behavior, one may conclude that the comparison of bar evolution in our modified gravity/DM models shows two different phases. In the first phase, which is common in all five models, an initial strong nonaxisymmetric m = 2 mode forms in the disk that is short-lived and dissipates after a short period of time. As we mentioned, at least for the MOG case, it is a bit difficult to discriminate between modified gravity and particle DM in this phase.

On the other hand, the second phase, which starts at intermediate times τ ≳ 200, reveals some significant differences between the modified gravity and DM models. The bar continues with an almost constant value in modified gravity, while in the cored DM models it almost constantly grows until the end of the simulation. On the other hand, in the cuspy DM model, it shows a period of constancy with a low value of bar amplitude and then, similar to the cored models, starts to grow until the end of the simulation. The other distinctive feature of bar amplitude in our models is the oscillations, which are only present in modified gravity models. This behavior might be due to the presence of different beating modes, which will be discussed in Section 4.5. However, one should be careful about the artifacts that are due to the choice of the center at which the quantities are being calculated. To be specific, let us assume that the center of the grid does not coincide with the centroid of the stellar bar. Then it is straightforward to verify that the bar amplitude given by A₂/A₀ would oscillate with roughly twice the frequency of the bar rotation (Ω_p). However, as we will show in the next subsection, this is not the case, and the oscillation period is not short enough to be one-half the rotational period of the bar. On the other hand, the GALAXY code uses McGlynn's method to find the position of the particle centroid (McGlynn 1984). With a suitable choice of the method's single parameter, this method is more sensitive to the internal condensation of the particles, where the bar is formed, and puts less weight on distant particles. Therefore, the code appropriately finds the centroid of the bar. This means that the oscillations are not artifacts and reveal a real feature. Furthermore, as we already mentioned, the grid is recentered every 16 time steps in order to minimize the numeric artifacts.

The same general behavior of the disks under modified gravity/DM halo was also presented in Roshan & Rahvar (2019), where the same gravitational models are studied for different mass profiles. Therefore, it could be concluded that, ignoring the first stage of evolution, the significant differences between bar evolution in modified gravity and DM remain unchanged when we vary the initial mass densities. Of course, the time evolution of A₂(τ) is not enough, and we still need to explore the other important properties of the stellar bar.

For a better comparison between our models, 2D projections of the particles in each model are presented in Figures 3–7. According to the face-on projection of these figures, the bar in the DM models seems to be more elongated than in modified gravity models. For the LHH model, the evolution is completed in a longer period of time. Other features apparent in these figures will be discussed in the next sections.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Another quantity that behaves differently in modified gravity and DM models is the bar pattern speed. Unlike in observations, where the required information should be derived from only one snapshot, in numerical simulations, the evolution of the disk and the rotation of its bar are measured more easily. Therefore, away from the complications of the observational methods, one can simply observe the bar rotation in the simulations to calculate the bar pattern speed ${{\rm{\Omega }}}_{p}$. The evolution of this quantity with time is demonstrated in Figure 8. From this figure, it is apparent that the pattern speed in the models of modified gravity is higher than in the DM models. This quantity decreases at a very slow rate in modified gravity, while it has a significant drop in the DM models. The decline in the bar's pattern speed in the cored DM models begins at τ ≃ 200, when the bar strength starts to increase and the second phase in the evolution of the bar amplitude has begun. This is because the halo is responsive and the bar slows down because of the dynamical friction caused by DM particles. On the other hand, for the LHH case, similar to the bar amplitude evolution, a constancy phase is visible in the pattern speed, which lasts until about τ ≃ 500, and then the second phase of the drop in pattern speed begins. In this period, the dynamical friction starts to act more effectively. However, such behavior is not observed in any mass models of the current and previous simulations of modified gravity (Tiret & Combes 2007; Roshan & Rahvar 2019). In other words, as expected, there is no dynamical friction in modified gravity models. Of course, the dynamical friction experienced by the bar via disk particles is negligible. Notice that the previously mentioned oscillations also appear in the pattern speed of the modified gravity models. Before moving on to discuss the ${ \mathcal R }$ parameter, which is of observational importance, let us discuss the LHH model with more emphasis on the amount of dynamical friction in this model.

Zoom In Zoom Out Reset image size

4.2.1. Cusp/Core Effect in the LHH Model

To better understand the behavior of the LHH model, we studied the Hernquist DM halo density during the evolution of the system. As one can see in Figure 9, at the beginning of the simulation, a cusp exists in the inner regions of the LHH model (black curves). As the simulation continues, the secular evolution in the system transforms this cusp to a core that flattens further during the simulation. It is clear from the gray curves in Figure 9 that such a transformation is absent in the initially cored Plummer halo model, which is almost unchanged during the simulation time.

Zoom In Zoom Out Reset image size

As discussed in Section 4.1, the initial rise in the bar amplitude of the LHH model is followed by a period in which the bar is weakened. According to the projected face-on snapshots of Figure 7, it is visible that during this period, there is an almost spherical configuration at the center, that is, a bulge, which might have formed due to the gravitational effect of the cuspy halo. On the other hand, it is well established that a massive bulge has stabilizing effects against global instabilities; for example, see Sellwood & Evans (2001). This symmetric configuration results in the suppression of the bar amplitude and also reduces the dynamical friction between the halo and disk particles because there is no substantial bulk motion to cause significant dynamical friction. We expect that a moving object inside a medium feels more friction than an object that rotates without any linear motion. Therefore, the pattern speed remains almost constant. On the other hand, the weak bar activity eventually destroys the cusp; see Figure 9. This mechanism has been reported in the literature to address the so-called core–cusp problem (El-Zant et al. 2004). Consequently, the bar starts to grow, as one can see in Figure 2. This directly means that the dynamical friction would rise and the pattern speed will decrease accordingly; see Figure 8.

It is worth mentioning that we also checked another cuspy model, Navarro–Frenk–White (NFW), alongside its MOG and NLG counterpart models, but with a different initial rotation curve. From these simulations, we also found similar results for the evolution of the disk under the effect of the cuspy DM halo of NFW. Because of this similarity, we omit these simulations for the sake of brevity.

Another useful quantity that is related to the bar pattern speed and has observational importance is the ${ \mathcal R }$ parameter, defined as

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

where D_L is the corotation radius and a_B is the bar semimajor axis. Observations indicate that regardless of the Hubble type, bars in disk galaxies appear to be fast rotators with $1\lt { \mathcal R }\lt 1.4$ (Aguerri et al. 2003, 2015; Cuomo et al. 2019; Garma-Oehmichen et al. 2020). However, hydrodynamical simulations in the ΛCDM paradigm report slow bars in disk galaxies with ${ \mathcal R }\gt 1.4$ (Algorry et al. 2017). This is the case also in isolated disk simulations; for example, see Debattista & Sellwood (1998). There are some proposals to address this discrepancy in the context of the DM hypothesis. For example, it is claimed that ultralight axionic DM particles suppress the dynamical friction at galactic scales and consequently leads to fast bars (Hui et al. 2017). However, the current debate on this problem in DM models has not yielded a final conclusion. On the other hand, this issue is simply resolved in modified gravity because of the absence of dynamical friction.

To calculate the ${ \mathcal R }$ parameter, at the first step one should find the corotation radius, the radius at which the rotational velocity of the pattern is equal to that of the particles (stars) in the disk. Using the pattern speed Ω_p introduced in the previous section, it is easy to find the corotation radius.

However, calculating bar length a_B is not a straightforward task, neither in real galaxies nor in simulations. Various methods have been proposed to measure this quantity. For example, visual estimation (Kormendy 1979; Martin 1995), ellipse fitting to galaxy isophotes (Wozniak et al. 1995), and Fourier decomposition on the surface brightness of the galaxy (Elmegreen & Elmegreen 1985; Ohta et al. 1990; Aguerri et al. 2000) are among the methods that have been widely used in the literature.

In order to distinguish between the behavior of our models regarding ${ \mathcal R }$, we used the Fourier decomposition method in finding a_B. In this scheme, Fourier components of the intensity (I_m) are calculated, and the bar (I_b) and interbar (I_ib) regions are introduced as

$$\begin{eqnarray}&&{I}_{b}={I}_{0}+{I}_{2}+{I}_{4}+{I}_{6}\end{eqnarray} \tag{ 8 }$$

and

$$\begin{eqnarray}&&{I}_{\mathrm{ib}}={I}_{0}-{I}_{2}+{I}_{4}-{I}_{6},\end{eqnarray} \tag{ 9 }$$

respectively, where I_m is the mth component of the Fourier decomposition. Finding the ratio of $\tfrac{{I}_{b}}{{I}_{\mathrm{ib}}}$ throughout the disk, Aguerri et al. (2000) show that the outer radius of

$$\begin{eqnarray}&&\displaystyle \frac{{I}_{b}}{{I}_{\mathrm{ib}}}\gt \displaystyle \frac{1}{2}\left[{\left(\displaystyle \frac{{I}_{b}}{{I}_{\mathrm{ib}}}\right)}_{\max }-{\left(\displaystyle \frac{{I}_{b}}{{I}_{\mathrm{ib}}}\right)}_{\min }\right]+{\left(\displaystyle \frac{{I}_{b}}{{I}_{\mathrm{ib}}}\right)}_{\min }\end{eqnarray} \tag{ 10 }$$

results in a good estimation of the bar length a_B. According to Athanassoula & Misiriotis (2002), the error in finding a_B using this method in numerical simulations is less than 4%, but it reaches $\simeq 8 \%$ for very thin bars.

By applying the above technique, the ${ \mathcal R }$ parameter could be determined for our models. However, from the projected face-on view of the models, it is clear that there are some stages, especially at the beginning of our simulations, where strong spiral waves exist in the system. Similar to the bar amplitude, these spiral arms affect the calculations of the bar length and yield a higher value for this parameter. This could result in artificially smaller values for the ${ \mathcal R }$ parameter. To avoid this effect, we consider time intervals of Δτ = 40, choose the minimum value of the bar length in each interval, and then report the related ${ \mathcal R }$ parameter (Hilmi et al. 2020). The results, namely the bar length and the ${ \mathcal R }$ parameter, are presented in Figure 10. It is clear from the left panels that in the second half of the simulation time in all of the DM models, the bar length increases with time. On the other hand, from the right panels, we see that the ${ \mathcal R }$ parameter is also an increasing quantity. This directly shows that the rate of change in the corotation radius of these models is higher than in the bar length.

Zoom In Zoom Out Reset image size

As you may see in Figure 10, the bar is in the slow regime (${ \mathcal R }\gt 1$) for all of the DM models, and also the ${ \mathcal R }$ parameter shows a growing trend. According to the previous studies, this behavior is expected because of the dynamical friction between the bar and the DM particles. However, in the modified gravity models, namely NLG and MOG, it is obvious that ${ \mathcal R }$ remains almost constant in the fast-bar regime ($1\lt { \mathcal R }\lt 1.4$), and the bar retains its rotation speed. It is apparent that the NLG model results in slightly higher values of the ${ \mathcal R }$ parameter than in MOG. It should also be mentioned that this property seems to be independent of the adopted disk/halo mass models. As seen in Figure 10, a similar trend in the LPH, LIH, and LHH models is apparent, although note that the delay in the evolution of the LHH model is also seen in the plots of bar length and ${ \mathcal R }$ parameter.

Therefore, regarding the results of pattern speed and ${ \mathcal R }$, we conclude that models under the effect of modified gravity show better compatibility with the observational results in comparison to DM halo models.

To have a more precise view of the evolution of our models, we studied the disk's behavior in the vertical direction. First, we focus on the inner radii, at which the bar dominates, and check the vertical behavior of the models with time. Then we will study the radial behavior of each model's thickness at a final time of τ ≃ 800. The two lower panels of each part in Figure 11 illustrate the time evolution of the mean $\left\langle z\right\rangle$ and rms $\sqrt{\left\langle {z}^{2}\right\rangle }$ thickness of the disk at two inner radii r = 1.1 and r = 2.1, respectively. In this figure, the top left, top right, middle left, middle right, and bottom parts present the results for NLG, MOG, LIH, LPH, and LHH, respectively. In all of the models studied in this work, except the LHH model, growth in rms thickness is visible at τ ≃ 100. Similar to the results of the previous sections, the rms thickness in the LHH model shows a delay in its evolution and starts its initial increase at about τ ≃ 200. However, in comparison to the MOG and DM models, NLG shows a sharper and more distinct increase, while for the other models, the growth happens at a slower rate. This is compatible with the higher initial bar growth rate of NLG, which was discussed in Section 4.1. Note that the evolution of the vertical behavior of the disk is linked to the bar strength evolution. Comparing the results of the rms thickness to bar amplitude in our modified gravity models, we see that immediately after the buckling instability, that is, τ ≃ 150 in NLG and τ ≃ 200 in MOG, the bars enter the oscillatory regime around a constant value. On the other hand, a second step-like behavior is visible in both the LIH and LPH models at τ ≃ 400. Such a second step-like behavior is not present in the rms curve of r = 1.1 for the LHH model, but a smooth growth in r = 2.1 is visible at around τ ≃ 800, which seems to be compatible with the slight change in the slope of the bar amplitude curve after this time.

Zoom In Zoom Out Reset image size

It should be stressed that each step-like behavior in rms thickness can be interpreted as the manifestation of the buckling instability. During this instability, the bar is weakened as the result of vertical resonances (Raha et al. 1991). It is interesting that in our DM models LPH and LIH this instability happens twice. The first buckling happens at the same time when there is a tangible reduction in the bar amplitude. However, note that for all of the DM models, the first buckling is stronger, and its effect on the bar amplitude is more explicit. For example, the rms thickness of the LIH model at r = 1.1 varies from 0.5 to about 2 during the first buckling, an enhancement by a factor of 4, but grows only from about 2.5 to 4, that is, by a factor of 1.6, in the second buckling. Consequently, the second buckling is not accompanied by a rapid reduction in the bar magnitude. Regarding the fact that A₂ is measured using the projected position of the particles in the x − y plane, in the second buckling the bar thickens in the vertical direction, and there is no substantial change in its projected width and length measured in the x − y plane.

From Figure 11, it is also apparent that in the DM models the trend of rms thickness at larger radii (dashed curves) in the second half of the simulation differs from the smaller radii (solid curves) and grows at an almost constant rate, while it stays almost constant at smaller radii. This could be because a peanut shape is formed in our DM models. Although the disk thickens during the first buckling, the peanut shape appears more explicitly after the second buckling in the LPH and LIH models. However, in the LHH model, although there is a smooth second buckling, there is a vivid peanut showing up after the first buckling. To see this fact more clearly, refer to the plots of the particles' projected positions in different planes of each model in Figures 3 and 6. In modified gravity models, there is only one buckling, and although we see a sudden increase in the rms thickness, there is no explicit twofold symmetric peanut.

The mean thickness, which is presented in the middle panel of each part in Figure 11, is a measure for the asymmetric distribution of particles around the z = 0 plane. A sudden change around the z = 0 plane in the $\left\langle z\right\rangle$ diagram happens at the same time as the growth in rms thickness.

The other quantity related to the disk's stability is the Toomre parameter Q, for which the time evolution has been illustrated in the upper panel of each part in Figure 11. We should recall that, to avoid local collapse and instability of the disk, the condition of Q > 1 should be fulfilled. From Figure 11, it can be seen that for all of our models, the Q parameter in small radii has a lower value than in larger radii. Starting from Q ≃ 1.5 at τ = 0, this parameter in MOG and NLG increases rapidly to Q ≃ 2 in R = 1.1 and remains almost constant until the end of the simulation. However, in the same radius in the LIH and LPH models, Q starts from 1.5 but reaches about 3 at the end of the simulation. On the other hand, the Q parameter in the LHH model behaves similarly in the other DM models at this radius until τ ≃ 800 and then continues to grow smoothly to about 4 at the end of its evolution.

The modified gravity and DM models show another difference at larger radii. Although they start from the same initial value for Q, NLG and MOG show a sharper growth at an earlier stage of the evolution. In the NLG model, this growth happens at about τ ≃ 50 and reaches ≃5 and remains constant until the end of the simulation. Also in the MOG model, we see that the growth starts at τ ≃ 50. However, the growth rate is not as fast as in NLG and lasts until τ ≃ 200, where the Q parameter reaches an almost constant value ≃4. On the other hand, in the DM models, this growth is smoother and starts from the second half of the simulation (τ ≃ 400 for the LPH and LIH models and τ ≃ 800 for the LHH model), where the second buckling instability happens in the disk. The final value of the Q parameter at the end of simulation reaches ≃7 in all of the DM models.

To have a better view of the overall behavior of the disk in the vertical direction, we plotted the mean and rms thickness in terms of radius for all of the models at τ = 800 in the left and right panels of Figure 12, respectively. The rms thickness in the LIH and LPH models is higher than in modified gravity models only at radii 1 ≲ R ≲ 2. This peak in rms thickness is a result of the peanut shape that occurs in the DM models. The peak in the LHH model at τ = 800 has a lower value and happens at a smaller radius. On the other hand, the rms thickness in the MOG and NLG models rises almost smoothly from the inner to outer radii and in most of the disk has higher values than in DM models. In other words, ignoring the inner parts, the disks are thicker in modified gravity models.

Zoom In Zoom Out Reset image size

From Figure 12 (left), it can also be seen that the mean thickness in NLG and MOG shows a deviation from zero. This means that, in accordance with the conclusion of Tiret & Combes (2007), our MOG and NLG models are flared and warped. This could also be seen in the edge-on projections of the disks (Figures 4 and 5).

It is necessary to mention that the vertical behavior of the modified gravity models considered in this work shows some similarities to MOND (Tiret & Combes 2007), but there are some differences. For example, the peanut shape explicitly appears in MOND simulations presented in Tiret & Combes (2007), whereas we do not confirm its existence in NLG and MOG. These differences are related to the distinct nature of these models. There is a fundamental acceleration ${a}_{0}\simeq {10}^{-10}{\rm{m}}\,{{\rm{s}}}^{-2}$ in MOND, beyond which these deviations from Newtonian dynamics appear. On the other hand, NLG's corrections originate in the nonlocal features of gravitation. There is no fundamental acceleration scale in this theory. The field equations are also different in these theories. In other words, although the gravitational potential Φ satisfies a linear differential equation in NLG and MOG, MOND's field equation is nonlinear. Therefore, in principle, we do not expect the completely same behavior for them in galactic dynamics. Further investigations on the comparison between these modified gravity models are required to reveal their viability using N-body galactic simulations.

All of the discussion presented in this subsection can be reduced to the fact that at the inner radii, in which the bar is dominating, modified gravity models lead to disks with smaller thickness; that is, the bars in DM models seem to be thicker and show stronger peanut configurations. However, the disk in modified gravity models at larger radii is more warped and flared and shows an overall higher thickness in comparison to the DM models. This is of great importance in the sense that it implies that the morphology of the disks could be different in modified gravity than in DM models. It is necessary to mention that in a similar phenomenology, it is recently claimed that the vertical structure of the galaxies could help to discriminate between modified gravity and DM (Lisanti et al. 2019a, 2019b).

As the final remark in this section, let us recall that the above-mentioned fact dealing with the thickness of the disks at small radii has also been reported in Roshan (2018) and Roshan & Rahvar (2019) for different mass profiles. Therefore, regarding the main purpose of this paper, it seems that the emergence of disks with smaller thickness at small radii compared to DM models is a model-independent feature in modified gravity models. Let us add our new finding to the previous results: disks at large radii are thicker in modified gravity models. In other words, ignoring the central region, the galactic disks seem thicker in modified gravity models compared to the standard picture.

In Figure 13 we have plotted the contours of the power spectrum with respect to the radius in the second phase of the simulations for τ > 200. This figure helps us to understand the oscillating behavior of the bar magnitude. It should be emphasized that although we use different modified gravity model parameters α and μ compared to Roshan (2018), we see the same oscillatory behavior reported in that paper. In Figure 13, the top left and top right panels belong to NLG and MOG, respectively. In both models, we see at least two distinct frequencies in the power spectrum. This means that there are two density waves with different frequencies propagating on the surface of the disk. Naturally, this yields an oscillatory behavior in the physical properties of the bar. Interestingly, the frequencies of these waves in NLG and MOG are close to each other. This is expected in that the pattern speeds of the bars in these models are also close to each other.

Zoom In Zoom Out Reset image size

In the bottom left and bottom right panels, we have shown the power spectra for the LIH and LPH models. The vertically aligned contours directly show that there is no constant frequency in the system. In other words, the frequency of the bar is varying with time. We know that it is a decreasing function with respect to time because of the dynamical friction. This feature is seen in both DM models. Note that the noisy concentration of the contours at large distances does not necessarily mean that there are several density waves because the main concentration of the baryonic matter is limited to smaller radii, as you may see in Figures 3 and 6. However, as also reported in previous studies, disks in modified gravity expand to larger radii. Therefore, the analysis of the power spectrum would be fully relevant at these radii.

As another test of the oscillatory behavior in our modified gravity models, we calculated the bar amplitude in the inner region of the disk, at r ≲ 4, where the bar mode dominates. As mentioned earlier, it is clear from Figure 13 that there are two distinct modes at radii r < 4 and r > 4 in these two models. The results are compared to the bar amplitude that was previously calculated for all of the disk particles (gray curves) in Figure 14. The value of the m = 2 mode in the outer disk is also presented in this figure. All of the values are scaled to the total number of the particles (mass) in each model. As is clear in Figure 14, the amplitude of the oscillations in the NLG and MOG models considerably changes by considering the particles inside r ≃ 4. It could be concluded that the oscillations in the bar amplitude and therefore in pattern speed are directly related to the existence of the second mode in the outer disk.

Zoom In Zoom Out Reset image size

To ensure that the results do not suffer from numeric artifacts, we have performed some tests. The conservation of energy is checked to be less than 2% until the end of the simulation. To make sure the main results are not affected by shot noise, we checked the results by changing the number of particles from 2 to 10 million. For illustration, some of the results are reported in Figure 15. As you may see in Figure 15 (left), similar to the models with 2 million particles, the initial peak in the bar amplitude happens earlier in NLG models. A similar behavior is seen in the LHH model, but the peak in LHH is slightly lower than in the NLG model, a result that is compatible with the result of the low-resolution simulation. This initial increase happens later and with a slower rate in the MOG model and then in the DM models, as was also seen previously. On the other hand, the drop of bar amplitude in the MOG model with N = 10⁷ is smoother, but at later stages, both modified gravity models reach an equal amount of bar amplitude again. Furthermore, the second phase of the increase in the bar amplitude is clearly seen in all of the DM models. Note that the constancy phase in the evolution of the LHH model is shorter and the increasing phase starts faster in comparison to the low-resolution simulation. Figure 15 (right) illustrates the evolution of pattern speed in our models with N = 10⁷ particles. As is apparent, the value of pattern speed in our modified gravity models is considerably higher than in the DM models. On the other hand, the fast drop in the LIH model and the convergence of the NLG and MOG models happen later in comparison to the N = 2 × 10⁶ case. However, it is clear that the main differences in the characteristics of DM and modified gravity models remain unaltered. Such a change in the pace of disk evolution as a result of the change in the number of particles, has also been reported in the previous studies (Ghafourian & Roshan 2017).

Zoom In Zoom Out Reset image size

Altering the time step from 0.01 to 0.005 leads to the same results as in our main simulations. This means that the time evolution is already converged and we have used a suitable choice for the time step. Furthermore, considering a higher resolution by changing the mesh dimensions to be 240 × 256 × 135, and changing lscale from 12.5 to 14.5, keeps the main results unaltered.