Evolution of spiral galaxies in nonlocal gravity

TL;DR

This study uses high-resolution simulations to compare the evolution of disk galaxies under nonlocal gravity versus standard dark matter models, revealing differences in bar formation and dynamics.

Contribution

It introduces a novel comparison between nonlocal gravity and dark matter models in galaxy evolution through detailed N-body simulations.

Findings

Higher bar instability rate in nonlocal gravity model

Bars are weaker and faster in nonlocal gravity

Differences in long-term galaxy dynamics between models

Abstract

We study the evolution of simulated disk galaxies in the context of a nonlocal theory of gravity. In this theory the appearance of the dark matter problem in cosmology and astrophysics is a manifestation of the nonlocality of the gravitational interaction. Using high-resolution N-body simulations, we investigate the dynamical evolution of disk galaxies and compare the result with the standard dark matter viewpoint. Specifically, we construct two exponential galaxy models, one in nonlocal gravity (NLG) and the other surrounded by a Plummer dark matter halo. Both systems start from the same baryonic matter distribution, particles velocities, random velocities and the initial Toomre's parameter. However, although the initial conditions are the same in both models, their long-term dynamics reveals some notable differences. For example, it turns out that the bar instability happens with higher rate in NLG model compared with the standard case. On the other hand, at the end of the simulation, we find that bars are weaker and faster in NLG compared with the standard case.

Equations

$$F_{\mu \nu }(R)=\frac{1}{R}{\int }_R{F}_{\mu \nu }(x^{\prime })d{x}^{\prime 4},$$

$$D(\mathbf{x},\omega )=ϵ(\omega )\mathbf{E}(\mathbf{x},\omega ).$$

$$D(\mathbf{x},t)=E(\mathbf{x},t)+{\int }_{-\infty }^{+\infty }G(\tau )E(\mathbf{x},t-\tau )d\tau ,$$

$$D^i(\mathbf{k},\omega )=\sum _j{ϵ}^{ij}(\mathbf{k},\omega )E^i(\mathbf{k},\omega )$$

$$D^i(\mathbf{x},t)=\sum _j\int d{x}^{\prime 3}\int d{t}^{\prime }{ϵ}^{ij}({\mathbf{x}}^{\prime },t^{\prime })E^j(\mathbf{x}-{\mathbf{x}}^{\prime },t-t^{\prime }).$$

$${\nabla }^2\varphi (\mathbf{x})=4\pi G\rho (\mathbf{x})+4\pi G\int q(\mathbf{x}-\mathbf{y})\rho (\mathbf{y}){\mathbf{d}}^3\mathbf{y},$$

$$q(r)=\frac{1}{4\pi {\lambda }_0}\frac{1+\mu r}{r^2}{e}^{-\mu r},$$

$$\phi(r)=-\frac{Gm}{r}\Big{\{}1+\alpha(1-e^{-\mu r})-\frac{\mu\alpha r}{2}\int_{\infty}^{r}\frac{e^{-\mu r^{\prime}}}{r^{\prime}}dr^{\prime}\Big{\}},$$

$$a(r)=\frac{Gm}{r^{2}}\Big{\{}1+\alpha\Big{(}1-\Big{(}1+\frac{\mu r}{2}\Big{)}e^{-\mu r}\Big{)}.\Big{\}}$$

$$\Sigma (R)=\frac{M_d}{2\pi R_d^2}{e}^{-R/R_d}$$

$$\displaystyle\rho(r)=\frac{3M_{h}}{4\pi b^{3}}\Big{[}1+\Big{(}\frac{r}{b}\Big{)}^{2}\Big{]}^{-5/2}$$

$${\sigma }_z^2(R,z)=\frac{1}{\rho (R,z)}{\int }_z^{\infty }\rho (R,z^{\prime })\frac{\partial \Phi }{\partial z^{\prime }}d{z}^{\prime }.$$

$$A_m(\tau )=\frac{1}{N}\sum _{j=1}^N{e}^{im({\varphi }_j(\tau )+\cot \gamma \ln r_j)}$$

$$A_m(\tau )={\int }_0^{\infty }{\int }_0^{2\pi }f(r,\varphi )e^{im(\varphi (\tau )+\cot \gamma \ln r)}rdrd\varphi$$

$$\mathcal{R}=\frac{D_L}{a_B}$$

Full text

galaxies: kinematics and dynamics– galaxies: spiral– instabilities– galaxies: bar growth

1 Introduction

Early numerical simulations revealed that disk galaxies are violently unstable against stellar bar formation, for example, see Miller et al (1970); Hohl (1971). This rapid bar formation was inconsistent with the observations in the 1960s. However, it turned out that a rigid spherical halo surrounding the galaxy, can substantially suppress the bar instability (Ostriker & Peebles, 1973). Currently, it is well-known that disks in live dark matter halos are more prone to bar instability compared to those in rigid halos (Athanassoula, 2002; Sellwood, 2016).

The dark matter particles distributed in a live halo can effectively slow down the bar pattern speed. More specifically the dynamical friction originated from gravitational interaction between the dark matter particles and the baryonic matter in the disk, is not negligible. This mechanism leads to slow bars at the end of dark matter halo simulations, for example in the cosmological context see Algorry et al (2017) and for isolated galaxies see Debattista & Sellwood (2000). However, this fact is inconsistent with the relevant observations on the pattern speeds. Theses observations show that spiral galaxies host fast bars (Aguerri et al, 2015). It should be noted that, for a fast bar the ratio of the corotation radius to the bar semi-major axis is smaller than $1.4$. This important issue will be discussed at greater length in the subsequent sections.

These problems along with other similar unresolved issues in the standard dark matter paradigm give rise to motivations for modified gravity approach. In this paper we study the galaxy evolution in a nonlocal gravity theory (Mashhoon, 2013, 2017) which tries to address the dark matter problem without invoking exotic particles (Rahvar & Mashhoon, 2014). Therefore, let us briefly review the similar studies in the context of other well-known modified gravity theories.

Modified Newtonian Dynamics (MOND) is one of the most successful theories to address the dark matter problem (Milgrom 1983; Famaey & McGaugh 2012). The local stability of disk galaxies in MOND has been studied in Milgrom (1989). There are few papers in which the evolution of spiral galaxies, with emphasis on the bar instability, has been investigated via N-body simulations, for example see Christodoulou (1991); Brada & Milgrom (1999) and Tiret & Combes (2007). The Latter provides more advanced simulations incorporating a large number of particles in the simulations. Therefore its results are more reliable than precedent papers. In this paper, it has been shown that the bar growth rate in MOND is higher than in Newtonian case. However, in MOND after reaching a maximum, the bar strength starts to decrease while in the dark matter model it does not stop increasing. In this case, the final magnitude of the bar at the end of the simulation in MOND is substantially smaller than the dark matter halo case.

Another theory which explains the flat rotation curves and the mass discrepancy in the clusters is a scalar-tensor-vector theory of gravity known as MOG in the literature (Moffat, 2006; Moffat & Rahvar, 2013, 2014). The local stability of disk galaxies in MOG has been investigated in Roshan & Abbassi (2015). On the other hand, the global stability of spiral galaxies has been studied in Ghafourian & Roshan (2017) and Roshan (2018). Both studies, using different codes, show that the bar growth rate is lower in MOG. Furthermore, the magnitude of the bar at the end of the simulation is smaller in MOG, and the pattern speed is faster than that of the standard halo model. The absence of the dynamical friction in MOG is the main reason for the latter case.

In this paper, we study the time evolution of exponential disk galaxies in the context of the nonlocal theory of gravity (NLG). This theory has been built on completely different bases compared with MOND and MOG but shares the same purpose to resolve the dark matter problem without introducing dark matter particles (Mashhoon, 2013, 2017). In this theory, gravity is similar to the electrodynamics in the non-vacuum medium and has a non-local behavior. We know this behavior from the concept of hysteresis in the magnetic materials (Jackson, 1975).

The non-locality in this modified gravity model heritages from the violation of the equivalence principle where it postulates a pointwise connection between the local frame of the accelerated observer with the background global inertial frame. Writing the general relativity in non-local form, in the weak field regime, we can write the modified Poisson equation where the potential from the all the space with a given kernel is engaged in this equation. The comparison of the dynamics from the NLG with the observation of galaxies and cluster of galaxies can explain the dynamics of large-scale structures without a need to the dark matter (Rahvar & Mashhoon, 2014).

The outline of this paper is as follows: In section 2 we briefly introduce NLG and discuss its weak filed limit. In section 3 we explain the numerical procedure to produce the initial conditions. Moreover we describe the essential ingredients of the code which evolves the point particles. We construct two models with identical initial conditions, one model in NLG without halo component, and the other in Newtonian gravity surrounded by a live dark matter halo. In section 4 we compare the galactic evolution in these models. Finally, results are summarized and discussed in section 5.

2 Nonlocal gravity: the weak field limit

It is natural to expect that relativistic effects are negligible in the secular dynamics of spiral galaxies. Excluding the central black hole, this fact is true almost everywhere in the galaxy. Therefore we need to use the weak-field limit of the theory. More specifically, in N-body simulations, we need the gravitational force between two point particles. In this section, we review the main features of NLG and its Newtonian limit.

The Lorentz invariance between the inertial observers is the central symmetry of the special relativity. In this theory, all the measurements are done locally both in the space and in the time. The idea of non-locality dates back to Bohr and Rosenfeld where they pointed out long ago that field determination might not be performed instantaneously (Bohr & Rosenfeld, 1950). This is a covariant theory, however the average value of the electromagnetic field is detemined in a finite volume of the space-time as

[TABLE]

where the size of space-time volume of $R$ might be defined by the natural volume of space as $(\hbar/mc)^{3}$ and time $\hbar/mc^{2}$.

There are also classical examples of non-locality in the electrodynamics (Jackson, 1975). For the case that the dielectric constant in a medium depends on the frequency of an external electric field, the displacement vector relates to the electric field as

[TABLE]

The Fourier transformation of this equation results in

[TABLE]

where $G(\tau)=1/2\pi\int^{\infty}_{-\infty}\left(\epsilon(\omega)-1\right)e^{i\omega\tau}d\omega$ and depending on the dielectric function, one can find the non-locality of the displacement vector in time (for detail see Chapter 7 of Jackson 1975). The extension of equation (2) is

[TABLE]

and in the Fourier space results in non-locality both in the space and time as

[TABLE]

In the superconductors where the mean-free path is comparable with the size of system, the Ohm’s law is no longer valid and dielectric tensor not only depends on the frequency but also it depends on the wavenumber.

The idea of non-locality is extended to the gravitation where for the non-inertial observers unlike to the general relativity, the equivalence of a local accelerating observer to an inertial observer is not valid. Here, we define a characteristic length of the system by $\lambda$ and a characteristic length for the acceleration of a reference frame by $L$ where the ratio of $\lambda/L$ would provide the degree of the non-locality.

We can define the characteristic length of $L$ by $L=c^{2}/a$ for linear accelerating frames and $L=c/\Omega$ for rotating systems. Now, let us calculate $\lambda/L$ for typical astrophysical and cosmological systems. (i) For the case of Earth using $a=g$ results in $\lambda/L\simeq 10^{-12}$. (ii) For the case of galaxy with the rotating velocity of $v_{rot}$, $L=c\lambda/v_{rot}$ where $\lambda$ is the size of galaxy, then $\lambda/L=v_{rot}/c\simeq 10^{-3}$. (iii), finally for the case of the Universe, the non-locality ratio is $\lambda/L=\frac{1}{2}(HR/c)^{2}$ where $H$ is the Hubble parameter and $R$ is scale of Universe in our concern. From these three scales we can see that for the larger scales and/or strong gravitational system, the non-locality is important as this effect depends on the combination of scale and the acceleration. The NLG can be written in terms of the teleparallel gravity theories since the final formalism is like the electromagnetism and we know how to make the electromagnetism non-local in a medium (Jackson, 1975).

The Einstein equations can be written in equivalent formalism so-called teleparallel gravity, namely, $GR_{\parallel}$. In this formalism the Riemannian structure of space-time has an additional Weitzenböck connection. The advantage of this formalism is that Einstein gravity equations are expressed in the form of covariant Maxwell equations. In other word, in this formalism general relativity is a gauge theory of the Abelian group of spacetime translations (Blagojević & Hehl, 2013). The advantage of this formalism is that in analogy to the electromagnetism, we can make these theories non-local. The gravitational potential in teleparallel gravity is given by the tetrad field which relates the metric of space-time to the local tangent Minkowski space as $g_{\mu\nu}(x)=e_{\mu}{}^{\hat{\alpha}}(x)e_{\nu}{}^{\hat{\beta}}(x)\eta_{\alpha\beta}$ where the connection is defined by $\Gamma^{\beta}{}_{\mu\nu}=e^{\beta}{}_{\alpha}\partial_{\nu}e^{\alpha}{}_{\mu}$. In this definition the tetrad field is globally teleparallel and the curvature of the Weitzenböck spacetime vanishes. The field equations in this formalism are expressed in terms of the gravitational field strength $C_{\mu\nu}{}^{\hat{\alpha}}=\partial_{\mu}e_{\nu}{}^{\hat{\alpha}}-\partial_{\nu}e_{\mu}{}^{\hat{\alpha}}$ which is similar to the electromagnetic tensor of $F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$. The result is similar to the Maxwell equations where according to the non-localized electromagnetic field, we can make gravity non-local in this formalism (for details see Mashhoon 2017). Actually we don’t have a microphysical description of the non-locality in GR unlike to what we have in the electromagnetism. However, we may assume the non-locality of GR as an emergent phenomenon from the quantum behavior of gravity, similar to what we have from the microphysics of the dielectric medium.

For the weak field approximation the tetrad field can be expanded around the flat space as $e^{\hat{\alpha}}{}_{\mu}=\delta^{\alpha}{}_{\mu}+\psi^{\alpha}{}_{\mu}$ where since $\psi^{\alpha}{}_{\mu}$ is the first order perturbation and in this approximation, there is no distinction between the anholonomic and holonomic indices. After detailed calculation, the modified Poisson equation in this theory can be written as

[TABLE]

where for $q(\bf{x-y})=0$ we have locality and recover the standard Poisson equation. The important characteristics of this equation is that for NLG the second term at the right hand side plays the role of dark matter. Moreover, the linearity of this equation in terms of the source is satisfied. This means that we can use the superposition of discrete sources for calculation of the overall potential for an N-body system. The kernel function has no fundamental definition and can be obtained from comparing the observational data and theory.

Here we use the kernel function of

[TABLE]

where $\lambda_{0}$ and $\mu$ are the constant parameters of this function. We substitute this kernel function in equation (6) and use Dirac-Delta function as a point-like source of gravity, the result from this integration is

[TABLE]

where we replace $\lambda_{0}$ with $\lambda_{0}=2/(\alpha\mu)$ and taking derivative from this equation we obtain the acceleration of a test mass particle as follows:

[TABLE]

We note that due the linearity of Poisson equation in terms of density, we can use the superposition of the individual particles on the overall gravitational acceleration.

This theory has been used for calculation of the acceleration of a test particle in the spiral galaxies (Rahvar & Mashhoon, 2014). The results from this analysis show the unique value of $\alpha=10.94\pm 2.56$ and $\mu=0.059\pm 0.028~{}{\text{k}pc}^{-1}$ can be used to describe the dynamics of these galaxies without need to the dark matter. Moreover, taking the adapted values for the parameters of this modified gravity model, we can also explain the dynamics of the cluster of galaxies without a need to the contribution of the dark matter.

3 Numerical method

One of the standard and high-resolution N-body codes which has been developed to study isolated galaxies is the GALAXY code Sellwood (2014). This code produces suitable initial conditions and gives the time evolution of the galactic models. The galactic models can include disk, bulge, and spherical dark matter halo. Almost all the well-known disk and halo density profiles are included in the code. A detailed description is available in the online manual Sellwood (2014). We recall that in NLG model, we do not need dark matter halos. Therefore to construct galactic models in this theory, we have modified only some parts of the code which deal with disks. As a necessary test, we measure the error in total energy during the simulations. This error is always smaller than $1$% in all the models. For simplicity, we have not included a bulge in our models.

For galaxy model in Newtonian gravity including a live halo, we use the hybrid N-body code described in (Sellwood, 2003, Appendix B). In this case two grids are used to find the gravitational field: a spherical grid for the live halo component, and a 3D cylindrical polar grid for the disk. Of course, our single component model in NLG uses only the latter grid.

3.1 Initial conditions

The initial dark matter model in Newtonian gravity is identical to the model ”EPL”111E: Exponential disk, P: Plummer sphere, L: Live halo in Roshan (2018). This model is consist of an exponential disk and a live Plummer spherical halo. The exponential profile for the disk is a satisfactory choice from the observational point of view. On the other hand, the Plummer sphere for the halo is chosen to obtain a better match between rotation curves of this model and the model in NLG. Moreover, this halo has been used in other studies where MOND’s consequences are compared with the standard dark matter paradigm, for example see Tiret & Combes (2007) and Combes (2014). Therefore, this profile enables us to compare our simulations in NLG with those of MOND. On the other hand, this halo provides rising rotation curves reminiscent of the late-type low luminosity galaxies. The exponential disk and the Plummer halo are given by the following densities respectively

[TABLE]

[TABLE]

where $M_{d}$ is the disk mass, and $M_{h}$ is the halo mass. Furthermore, $R_{d}$ is disk length scale, and $b$ is a characteristic radius. The disk is truncated at $r=4\,R_{d}$ and we use a cubic function in order to taper the exponential density smoothly to zero over the range $3.2\,R_{d}<R<4\,R_{d}$. The halo is truncated at $r=3\,b$. We set the mass of the disk and the halo as $M_{d}$, and $8\,M_{d}$ respectively. Furthermore we assume that $b=10\,R_{d}$. Notice that we use units such that $R_{d}=M_{d}=1$. The Newtonian gravitational constant is scaled as $G=1$. With these choices, the unit of velocity is $V_{0}=(GM_{d}/R_{d})^{1/2}$, and the time unit is written as $\tau_{0}=R_{d}/V_{0}=(R_{d}^{3}/GM)^{1/2}$. In this unit, a suitable choice for our models can be $R_{d}=2.6$ kpc and $\tau_{0}=10$ Myr, which yields $M_{d}\simeq 4\times 10^{10}M_{\odot}$ and $V_{0}\simeq 254\,\text{km}/\text{s}$.

In this paper, the time duration of the simulations is $800\tau_{0}$, and in few cases $500\tau_{0}$. The time step is $\Delta t=0.01\tau_{0}$. To obtain a better efficiency the time step is increased by successive factors of $2$ in three radial zones. Moreover to make sure that there is no dependence in the outcome on the time step $\Delta t$, we use both shorter and longer time steps, to demonstrate that our choice is appropriate. In all the models we use the cubic spline softening length (Monaghan, 1992) $\epsilon=0.16\,R_{d}$. To avoid numeric artifacts the grid is recentered every $16$ time step.

Since the initial matter distribution has been determined, the rotation curve, $v_{c}$, can be simply calculated. Let us briefly discuss the corresponding velocity dispersions. It is natural to expect that the particles, at least at $\tau=0$, move on nearly circular orbits. Although in real spiral galaxies, most of the stars are supposed to travel on nearly circular orbits, we are not constructing initial conditions identical to a real galaxy. Rather we construct a suitable background disk for studying its global stability in a transparent way. In other words, the simulations start with a disk of particles supported in equilibrium almost entirely by rotation. This cold disk enables us to observe the occurrence and evolution of the stellar bar as a pressure supported configuration. Consequently, there is a relation between radial and azimuthal dispersions as $\sigma_{\phi}=\kappa\sigma_{R}/2\Omega$, where $\Omega(R)=v_{c}/R$. On the other hand, working in the epicyclic approximation, we set the radial velocity dispersion $\sigma_{R}$ in such a way to establish the local stability criterion $Q>1$. Where $Q$ is the so-called Toomre’s stability parameter defined as $Q=\sigma_{R}/\sigma_{R,\text{crit}}$, where $\sigma_{R,\text{crit}}=3.36G\Sigma/\kappa$ and $\kappa$ is the epicyclic frequency (Toomre, 1964). We set $Q\simeq 1.5$ everywhere on the surface of the disk. This fixes the radial dependence of both radial and azimuthal velocity dispersions.

It is necessary to mention that the disk has a Gaussian density profile in the vertical direction $z$ with the scale height $z_{0}=0.05R_{d}$. To find the vertical velocities, the vertical one dimensional Jeans equation is integrated as

[TABLE]

Notice that the isotropic distribution function (DF), see Dejonghe (1987) for more details, has been used for the Plummer halo in our simulations. However, when the halo is combined with the exponential disk, one needs to recompute the equilibrium DF for the composite system. The GALAXY code uses a halo compression algorithm described by Young (1980) and Sellwood & McGaugh (2005). The algorithm is based on the fact that both radial action and the total angular momentum are conserved during compression. Finally, the particles are selected from the new DF function following the procedure studied in Debattista & Sellwood (2000).

Now, the initial conditions for the numeric galactic model in Newtonian gravity have been completely determined. Let us briefly discuss the model in NLG. Hereafter we call this model as ”ENLG” (Exponential disk in NLG). As already mentioned, there is no dark matter halo in this model. Therefore, since we have restricted ourselves to the bulge-less simulations, ENLG is a single component disk model. The disk is exponential and its physical parameters are the same as those of the model EPL. The only difference is that the gravitational force between particles is different from the standard case and is given by (9). More specifically, since the theory has two free parameters $\alpha$ and $\mu$, we are able to reproduce almost the same initial conditions produced with the two-component model EPL. However, as we shall see, the long-term evolution of the models faces significant differences.

As our final remark in this section, notice that for each responsive component we employ three million particles. For example, in EPL we employed three million particles for the disk and accordingly three million particles for the live halo. Similarly, the ENLG model has three million particles. In order to ensure that the results do not depend on the number of particles, we increase the number of particles up to $N=2\times 10^{7}$.

3.2 Rotation curves

Working in the units introduced in the previous subsection, we have plotted the rotation curves in Fig. 1. The blue curve belongs to the model ENLG, and the solid black curve is the corresponding curve for the model EPL. Furthermore, just for comparison, we have shown the rotation curve of model EM, introduced in Roshan (2018), with a dashed curve. This model is a single component exponential disk in MOG. The gravitational force in MOG and NLG includes similar correction terms. Therefore one may expect that these theories predict similar evolution for disk galaxies. However, as we shall show, this is not the case.

We emphasize again that for a meaningful comparison between different theories, the galactic models should start with the same initial conditions. From Fig. 1, we see that the initial rotation curves are similar. Although we have not plotted the initial velocity dispersions, and other functions like $Q(R)$, they also follow similar profiles.

For a given $\alpha$ and $\mu$ one may first find the rotation curve of the ENLG model, and then change the halo parameters in EPL model until the corresponding rotation curve gets almost identical to that of ENLG model. Of course, one can first find the rotation curve of EPL, and then vary $\alpha$ and $\mu$ in the model ENLG to fit the rotation curve. Using the second approach, we have fixed $\alpha=10.9$ and $\mu=4.5\times 10^{-3}R_{d}^{-1}$ in our simulations. Notice that although the $\alpha$ parameter is close to the observed values, we do not have to use observational values. In fact we have fixed $\alpha$ and $\mu$ by fitting the rotation curve of a numeric model EPL. Consequently, we do not expect that these values get identical to those obtained from fittings to real rotation curves.

4 Results

With the full set of initial conditions in hand, we are ready to discuss the time evolution of the models. Notice that the GALAXY code determines the gravitational forces on the cylindrical polar grid using sectoral harmonics $0\leq m\leq 8$. On the other hand, it uses the surface harmonics $0\leq l\leq 4$ to compute the gravitational force on the spherical grid. Albeit it should be noted that this numbers can be changed with the user. We emphasize that for our main model in this paper, i.e., the model ENLG, we have used the softening length $\epsilon=0.16\,R_{d}$ and grid of $193\times 224\times 125$ mesh points. In all the models the grid has been regulated to be large enough to prevent particles from escaping the grid. Moreover to make sure that the results are insensitive to our choice for $\Delta t$, we used shorter and longer values for the time step as $\Delta t=0.005$, $0.01$ and $0.03$. The results are qualitatively the same. Therefore we report the results obtained with $\Delta t=0.01$.

As already mentioned there are two main quantities which are of special importance for us, the stellar bar amplitude and the pattern speed of the bar. Both quantities are sensitive to the existence of the dark matter halo. In other words their time evolution is correlated with the physical properties of the halo. Therefore it is natural to expect serious differences in theories where there is no dark matter halo. In the following subsections, we study these quantities with details in the context of NLG and compare them with those in MOG and the standard dark matter model.

4.1 Bar growth

The amplitude of the non-axisymmetric features on the surface of the disk can be determined by the complex Fourier coefficient

[TABLE]

where $(r_{j},\phi_{j})$ are the position of the $j$th particle at time $t$ in the cylindrical polar coordinates. Moreover $\gamma$ is the pitch-angle of an $m$-fold logarithmic spiral. Logarithmic spirals make a complete basis set and consequently can be used to expand any two-dimensional function $f(r,\phi)$ as follows

[TABLE]

here $f$ is the surface density and is given by a set of delta functions. Therefore one may verify the Fourier coefficient (13). Notice that the summation is over the disk particles and does not include halo particles in the model EPL. The bar amplitude, for which $m=2$, is computed by the ratio $A_{2}/A_{0}$. For simplicity, we assume that $\cot\gamma=0$. In this case, $A_{2}$ contains both bar-like and spiral contributions. However, as we shall see, the spiral pattern exists in a short period of time compared with the time duration of the simulations, and we mostly deal with bar-like structures.

The evolution of the bar in our models, has been shown in Fig. 2. The blue curve belongs to our main model ENLG, the solid black curves are the corresponding curve for the model EPL. Furthermore, for comparison, the evolution of the bar for the model EM, has been illustrated with the black dashed curve. It is clear that the bar instability develops more rapidly in NLG compared with the standard dark matter model. More specifically, the exponential growth rate, i.e., $e^{\omega\tau}$, obtained from the top panel in Fig 2, in ENLG, EPL and EM models are $\omega\simeq 0.097$, $\omega\simeq 0.067$ and $\omega\simeq 0.058$ respectively. This means that growth rate in model ENLG is $67$% ($45$%) larger than that of the model EM (EPL). On the other hand we see that MOG has stabilizing effect, for more details on the stabilizing effects of MOG we refer the reader to Ghafourian & Roshan (2017) and Roshan (2018). It is interesting that although both theories, i.e., NLG and MOG, add almost similar Yukawa corrections to the gravitational force in the weak field limit, the onset of the bar instability and its evolution are quite different.

The rapid bar growth in NLG is reminiscent of the rapid growth in MOND. It is shown in Tiret & Combes (2007) that bar grows faster in MOND compared with the standard case. However, in MOND the maximum bar amplitude remains almost constant for a relatively long time interval ($4$ Gyr) and then decreases rapidly and stays almost constant at a much lower value, see Fig. 7 in Tiret & Combes (2007). On the other hand, in NLG the bar magnitude starts to increase at $\tau\simeq 10$, and reaches its maximum at $\tau\simeq 60$, see Fig. 2. This point is also clear from Fig. 3. In this figure we have projected the position of a small fraction of the particles on the Cartesian planes. This figure shows the evolution of the particles in the first part of the simulation where the thickness of the disk remains small. Furthermore, for better visualization the points have been smoothed, especially in the central part of the disk, using the Python module scipy.stats.gaussian_kde. Naturally, from dark red to blue color, the surface density decreases. This smoothed scatter plots are helpful to discriminate and track the evolution of the bar in a visual way. We see that at $\tau\simeq 57.6$ a spiral density wave is excited and rapidly fades its arms.

It is important to mention that in the EPL model, the bar magnitude experiences a relatively rapid growth at later times. In fact, the Newtonian bar can strengthen again (after the main drop), by exchanging angular momentum with the dark halo Athanassoula (2002).

As it is clear from Fig. 2, unlike in MOND, the bar magnitude starts to decrease and reaches its minimum value at $\tau\simeq 150$. We see that after this minimum, a new phase in the evolution is started and the bar grows with a very small slope. The second part of the evolution, where the bar’s magnitude is small and the thickness of the disk starts to increase, has been illustrated in Fig. 4. It should be noted that there are two main features apparent in figures 2 and 4. We see that the bar strength oscillates with an almost constant frequency in the second part of the simulation. Finding the power spectra of the model is helpful for revealing the reason for this oscillation. On the other hand it is clear from Fig. 4 that the thickness of the disk suddenly starts to increase with time, and it seems that the mean thickness, i.e., $\langle z\rangle$, oscillates around $z=0$ for a short period of time. In the following, we study these both cases with more detail.

4.2 Power spectra of the model ENLG

Let us start with the rapid oscillations in the magnitude of the bar. One should note that, as we already mentioned, the bar magnitude $A_{2}$ contains the contributions from spiral structures and the bar-like pattern as well. Consequently, this oscillation can be interpreted as the existence of beating/coupling between different modes. It is interesting that this oscillation exists only in modified gravity models, i.e., in the models ENLG and EM.

For revealing the existence of different modes, one may naturally find the power spectrum of the mass distribution. To do so, at each output time we first perform a one-dimensional Fourier transform of the surface density in $\phi$ (notice that the transform (14) is two-dimensional and taken over $r$ and $\phi$). The result for each mode $m$ will be a function of $r$ and $t$ as $A_{m}(r,t)$. In order to find the power in terms of $r$ and frequency $\omega$ for each sectoral harmonic $m$, we calculate the Fourier transformation in time of the coefficients $A_{m}(r,t)$ and show it by $\tilde{A}_{m}(r,\omega)$. In this case the power spectrum is given by $|\tilde{A}_{m}(r,\omega)|^{2}$.

The isocontours of the power spectrum for our main model ENLG has been shown in Fig. 5. We have shown power spectra for three different modes $m=2$, $3$ and $4$, where there are $2\times 10^{7}$ particles in the simulation. Furthermore, considering the overall behavior of bar magnitude shown in Fig. 2, we have taken the temporal Fourier transformation over two different time intervals. The top row in Fig. 5 is for the first part of the evolution $10<\tau<160$, and the bottom row is for the second part of the evolution, i.e., $160<\tau<800$.

In this figure, the solid (green) curves mark $m\Omega_{c}$ and the dashed (red) curves show $m\Omega_{c}\pm\kappa$. Where $\Omega_{c}$ and $\kappa$ denote the angular frequency for the circular motion and the Linblad epicycle frequency respectively. Notice that instead of theoretical curves, we have computed these frequencies using the initial velocities at $\tau=0$. Subsequently, one can see a rapid change in the red curves at $r\simeq 4$, which is related to the surface mass density truncating procedure explained in the section 3.

It is clear from the top row that, in the first part of the simulation there are two bar-like modes inside the corotation radius. Notice that each horizontal ridge determines a coherent density wave. In fact these $m=2$ modes, with close frequencies $\omega\simeq 0.6$ and $\omega\simeq 0.52$, start from the inner parts of the disk and stop at the corotation. One should note that the oscillations in the bar magnitude cannot be a consequence of beating between these two modes. The beat period, in this case, is large and does not coincide the period that can be inferred from Fig. 2 ($\tau_{b}\simeq 20$). Also there is a major contribution from the $m=4$ mode with frequency $\omega\simeq 0.94$. The $m=3$ mode has not an important impact on the system. In fact by plotting its magnitude it turns out that $A_{3}/A_{2}$ and $A_{3}/A_{4}$ remain small.

In the second part of the simulation, there are three $m=2$ modes inside the corotation, see bottom row in Fig. 5. We have shown them with blue dot dashed straight lines at $\omega_{1}\simeq 0.22$, $\omega_{2}\simeq 0.35$ and $\omega_{3}\simeq 0.52$. The last one is the long-lived dominant mode which was excited from the beginning of the simulation. On the other hand, we see that there are also two $m=4$ density waves with frequencies $\omega_{5}\simeq 0.87$ and $\omega_{6}\simeq 1.5$. The existence of these five waves is the main reason for the rapid oscillation in bar amplitude. Since we deal with several modes, analysis of the coupling and beating between the modes is somehow complicated. However, it seems that the beating between modes $\omega_{1}$ and $\omega_{3}$, which yields a beating frequency $\omega_{b}=\omega_{3}-\omega_{1}$ and period $\tau_{b}\simeq 20.9$, is responsible for the oscillations. This frequency is consistent with that of obtained from Fig. 2.

As our final remark in this subsection, it is interesting to mention that two bar-like beat waves can appear as an $m=4$ mode. The frequency of the resultant $m=4$ mode with frequency $\omega_{5}$, can be obtained by the sum of the frequencies of $m=2$ modes, see Masset & Tagger (1997) for more details. From this perspective, we see that the slower $m=4$ mode can be considered as the outcome of an interaction between bar-like modes with $\omega_{3}\simeq 0.52$ and $\omega_{2}\simeq 0.35$.

4.3 Vertical structure of the model ENLG

Now let us return to the vertical structure of the disk in the ENLG model. As we already mentioned, it is clear from Figs. 3 and 4, that there is a bulk oscillation around $z=0$ plane. Furthermore the rms thickness $\sqrt{\langle z^{2}\rangle}$ of the disk suddenly starts to increase at $\tau\simeq 150$. To see this point in a more quantitative way, we have plotted the time evolution of rms thickness at two different radii in the top panel in Fig. 6. More specifically, the solid and dashed curves belong to the model with $N=2\times 10^{7}$ and $N=3\times 10^{6}$ respectively. Moreover, in all the panels of Fig. 6 green and red colors belong to quantities at $R=1.1$ and $R=2.1$ respectively. It is clear that, as expected, there is no significant difference between simulations with different $N$.

In agreement with Fig. 4, we see that that $\sqrt{\langle z^{2}\rangle}$ has a rapid and step-like behavior around $\tau=150$. More specifically during the buckling, $\sqrt{\langle z^{2}\rangle}$ increases from $0.05\,R_{d}$ to $0.3\,R_{d}$. From the middle panel, it is seen that at the same time $\langle z\rangle$ oscillates around the plane of the initial disk. It is interesting to mention that the rapid rms thickness growth coincides with a substantial decrease in the bar magnitude. This behavior can be simply interpreted by considering the fact that particles resonate in the $z$ direction. In other words, the vertical oscillation can be amplified to form a peanut shape (Combes & Sanders, 1981). On the other hand, we know that the primary effect of the thickness is to reduce the growth rate of the Jeans instability in the system, see Romeo (1992). Consequently, it is natural to expect a reduction in the bar magnitude through the stabilizing effects rooted in the thickening of the disk.

The Toomre’s parameter is another quantity directly related to the stability of the disk. As already mentioned this parameter accounts for the local stability of the system. To be specific, the disk is locally stable against small perturbations provided that $Q>1$. On the other hand, it is well-known that the swing amplification is extremely sensitive to the value of $Q$. When $Q$ is close to unity, this mechanism works effectively to form $m=2$ patterns. However, the amplification factor decreases rapidly with increasing $Q$ (Julian & Toomre, 1966). It also provides additional insight into the heating rate due to the activity of the bar. In the bottom panel of Fig. 6, we have plotted the evolution of the $Q$ parameter in two different radii. It is clear that $Q$ starts to increase when the bar is developed in the disk at $\tau\simeq 50$. This directly means that the occurrence of the bar instability heats the disk. Moreover, the bar influences the larger distances more effectively than the inner radii. After $\tau\simeq 75$ the parameter $Q$ gets almost constant in both radii as $Q(1.1)\simeq 2.25$ and $Q(2.1)\simeq 3.75$. It is interesting that the buckling instability which happens at $\tau\simeq 150$ does not have any impact on the $Q$ parameter. In other words, $Q$ increases due to the bar activity and not thickening of the disk by buckling instability. This directly means that the bar instability seriously heats the disk, while the short-term asymmetry around the disk plane caused by the buckling instability does not. The relatively large value for $Q$ along with the buckling mentioned above prevents further bar growth. The behavior of $Q$ with respect to time in the ENLG model is similar to what happens in MOND. It is shown in Fig. 17 of Tiret & Combes (2007) that $Q$ rises in MOND and then gets almost constant. However, it grows continuously in the dark matter model.

Comparing the vertical structure of the model ENLG with the models EPL and EM studied in Roshan (2018), one may conclude that the final rms thickness of the disks in NLG is similar to those of standard dark matter model (EPL). On the other hand, the model in MOG has a significant difference in the sense that its final rms thickness is substantially smaller than in the dark matter mode. In other words, our simulation shows that, unlike in MOG and MOND which lead to stronger peanuts Tiret & Combes (2007), the model in NLG produces similar peanuts compared with the model in standard dark matter paradigm.

4.4 Bar’s pattern speed $\Omega_{p}$

Beside the bar growth in NLG, the other important quantity which is also related to the dark matter problem is the bar’s pattern speed $\Omega_{p}(\tau)$. In this section, we put emphasis on the crucial significance of the pattern speed’s evolution in order to discriminate between dark matter particles and modified gravity. We have measured $\Omega_{p}$ in terms of time and the result has been illustrated in Fig. 7. The top and bottom panels show the evolution of $\Omega_{p}$ for $N=3\times 10^{6}$ and $N=2\times 10^{7}$ particles respectively. The blue curve belongs to the model ENLG. On the other hand solid and dashed black curves are corresponding curves in EPL and EM models.

As expected the pattern speed slows down with time in the EPL model. The dynamical friction between halo and disk particles is the main reason for this reduction. In our main model, i.e., ENLG, after $\tau\simeq 150$ the bar magnitude gets almost constant and we can assign a well defined pattern speed. As we see in Fig. 7, the pattern speed starts to decrease with a small slope. The small amplitude oscillation in the pattern speed is a result of the beating between different modes. One can easily verify that this oscillation’s frequency is consistent with the beating frequency discussed in the previous section. Since $\Omega_{p}$ is almost constant in ENLG, we can say that dynamical friction between disk particles is not effective in this case. This is also the case in the EM model shown with dashed curves in Fig. 7. It is interesting that the dynamical friction does not appear in MOND as well Tiret & Combes (2007). It should be noted that the dynamical friction not only is important in the dynamics of isolated galaxies but plays role in the interaction between galaxies. For example, it is a key mechanism to increase the merger frequency of spiral galaxies during the cosmic evolution Boylan-Kolchin et al (2008). Therefore, this issue deserves more precise investigations to reveal observable differences between standard dark matter paradigm and modified gravity.

The relevance and importance of the pattern speed in modified gravity theories get more obvious when almost all the current observations show that spiral galaxies host fast bars Aguerri et al (2015). However, dark matter simulations for both isolated galaxies and cosmological simulations in $\Lambda$CDM predict an opposite result. More specifically, these simulations show that, mainly due to the dynamical friction, simulated spiral galaxies host slow bars, for a pioneer paper that pointed out this challenge in isolated galaxies see Debattista & Sellwood (2000). Furthermore for a recent cosmological simulation that also confirms this discrepancy see Algorry et al (2017).

Therefore it seems necessary to compare our simulation with the pattern speed data. Notice that $\Omega_{p}$ cannot directly specify that whether the bar is fast or slow. Instead, it is convenient to use the ratio of the corotation radius to the bar semi-major axis $\mathcal{R}$ defined as

[TABLE]

where $D_{L}$ is the corotation radius and can be obtained by measuring the pattern speed $\Omega_{p}$. On the other hand, $a_{B}$ is the bar semi-major axis. In this case the bar is fast if $\mathcal{R}\lesssim 1.4$, and is slow if $\mathcal{R}\gtrsim 1.4$. At each time step $\tau$, we measure the pattern speed $\Omega_{p}$ and the angular velocity $\Omega(R)$. Therefore it is straightforward to find the coratation radius by solving $\Omega(D_{L})=\Omega_{p}$. As we see in Fig. 7, the pattern speed is almost constant in the ENLG model. Accordingly, the corotation radius $D_{L}$ does not move throughout the disk. Therefore the time variation of $a_{B}$ specifies the evolution of the $\mathcal{R}$ parameter.

Now our main task is to measure the bar semi-major axis. To measure $a_{B}$ in terms of time, we use a simple procedure already explained in Roshan (2018). At a given time $\tau$, working in the projected $x-y$ plane, we find a straight line parallel to the bar which crosses bar’s center. This line can be simply expressed as $y=\tan\phi(\tau)x$, where $\phi(\tau)$ is the angular displacement of the bar. Then we assume a rectangle with width $\Delta L$ and length $5$ around this line. Finally, the rectangle is divided into small elements with equal area, and the surface density at each element is calculated. We choose the bar length to be the radius at which the density of the particles is less than $B$ percent of the central element. We vary $\Delta L$ between $0.2$ and $1$, and choose three values for $B$, i.e., $10$, $20$ and $30$. Naturally, these variations cause an error bar in the final value of $a_{B}$ and consequently in $\mathcal{R}$ at each time.

The result has been illustrated in Fig. 8 for the model ENLG. In the left panel, we have shown the bar semi-major axis $a_{B}$. We see that bar starts with a relatively large length and then decreases to reach a minimum at $\tau\simeq 200$ right after the buckling instability. Then $a_{B}$ starts to grow and this growth continues to the end of the simulation. On the other hand, we have plotted $\mathcal{R}$ in the right panel. We see that due to the existence of a minimum in $a_{B}$, there is a maximum in the time evolution of $\mathcal{R}$. In other words, the growth of the bar length together with the constancy of the corotation radius enforces $\mathcal{R}$ to decrease toward the border $\mathcal{R}=1.4$. Notice that observed bars lie in the interval $0.9\lesssim\mathcal{R}\lesssim 1.4$.

It should be stressed that in the EPL model, dynamical friction reduces the pattern speed and consequently the corotation radius increases with time. Therefore $\mathcal{R}$ is an increasing function of time and reaches $\mathcal{R}\simeq 4$ at the end of the simulation, see Roshan (2018). On the other hand, the model EM in MOG has a similar behavior as in ENLG. In other words, these simulations show that isolated galaxies in NLG and MOG lead to faster bars compared to the model in standard dark matter paradigm. Of course, the ratio $\mathcal{R}$ obtained in our simulations for ENLG (and MOG) are still higher than the relevant observations. Nonetheless, they are much smaller than similar simulations with dark matter halo. Therefore one may conclude that we see a better agreement in the context of modified gravity theories. However, it should be emphasized that the model studied in this paper is an ideal model in the sense that we have not included the bulge component. Consequently, still more investigations with more realistic galaxy models are required to shed light on this important issue.

As our final remark, we emphasize that the dynamics of the pattern speed is a rare situation in which modified gravity and particles dark matter behave completely opposite to each other. We know that modified gravity theories are presented to play the same role as the dark matter particles to explain the rotation curves of spiral galaxies or the cosmological structure formation. But this is not the case when studying the bar evolution, and it seems that some features of particles, like the dynamical friction, cannot be recovered by modifications in the gravitational theory.

5 Discussion and Conclusion

In this paper, we studied the evolution of a numerical galactic model in the context of a nonlocal gravity theory. More specifically, we have modified and tested the GALAXY code to include the linear effects of NLG. Our numeric model, called ENLG, consists of an exponential disk, and its initial properties are identical to the model EPL. The model EPL consists of a spherical Plummer halo and an exponential disk. Therefore we have constructed two identical models one in NLG without dark matter halo and the other in standard paradigm including a spherical halo. The main purpose of this study is to compare the evolution and dynamics of the stellar bar in the above mentioned theories. Since the gravitational force is somehow similar to that of MOG, for the sake of completeness, we have also compared the results with MOG.

Although models start with the same initial conditions, the dynamics of the bar is substantially different. For example, we found that, at least for the exponential disks, the bar instability is stronger in NLG. More specifically, the stellar bar in nonlocal grows with a higher rate compared with the dark matter model and MOG. We recall that a rigid halo can suppress the bar instability and a live halo, in principle, can trigger and even enhance the instability (Athanassoula, 2002; Sellwood, 2016). From this perspective, the NLG effects appear like a live halo.

When the bar instability occurs in ENLG, and accordingly the bar reaches its maximum strength, the resonance of the particles vertical to the plane of the disk, increases the thickness. In other words, buckling instability occurs and finally, a peanut shape appears in the simulation. The peanut is meaningfully stronger than in MOG and is comparable to that of in dark matter halo model.

We also studied the long-term evolution of the stellar bar in the model ENLG. Unlike in the dark matter model, the pattern speed does not decrease with time. More specifically the dynamical friction is not effective in the ENLG disk, while it is a key feature in the dark matter halo model EPL. Finally, we measured the ratio of corotation radius to the bar semi-major axis, i.e., $\mathcal{R}$, to compare the pattern speed with the observations. We found that after $\tau\simeq 200$, $\mathcal{R}$ starts to decrease and at the end of the simulation reaches the border $\mathcal{R}=1.4$. While in the dark matter model this parameter increases with time and reaches much higher values. In other words, NLG predicts faster bars in contrary to the standard case which leads to slower bars. This totally opposite prediction along with the relevant observations may help to distinguish between modified gravity and dark matter particles.

As our final remark in this paper, for future works, it should be noted that more realistic galactic models including initial spherical bulge are required to obtain a reliable comparison between N-body simulations in modified gravity and observations. Furthermore, more accurate analysis like 2D line-of-sight kinematics of the simulations would be helpful to derive the main properties of bars and peanuts in a way suitable for comparison with observations, for example see Iannuzzi & Athanassoula (2015); Molaeinezhad et al (2016). To the best of our knowledge, such analysis has not been done in modified theories of gravity. On the other hand, the cosmological simulations have not been developed in the context of modified theories of gravity without cold dark matter. Therefore a huge numerical and theoretical effort is still required to decide about the viability of these theories in the dynamics of isolated galaxies and in the cosmological galaxy formation scenario.

The main part of the calculations have been done on ARGO cluster in International Center for Theoretical Physics (ICTP). MR would like to thank ICTP for giving access to ARGO cluster. MR would like to appreciate Bahram Mashhoon for valuable discussions on the foundation of the non-local theory of gravity. Furthermore, we are grateful to Shahram Abbassi for providing us with a high performance computer. Also we would like to thank referee for his/her useful comments that improved this work. This work is supported by Ferdowsi University of Mashhad under Grant No. 2/47321 (02/05/1397).