An Interesting Equation Leads to the MOND but Does not Rule out the Existence of Dark Matter

1. Introduction

Dark matter and dark energy are the two major puzzles in cosmology today[1]. Dark energy is proposed to explain why universe is accelerating, and dark matter is proposed to solve the abnormal rotation curve of disk galaxies. In the standard cosmological model, Λ Cold Dark Matter (ΛCDM), the component of dark energy is about 70%, the component of dark matter is about 25%, and the other component we can see is barely about 5%[2]. Although the standard cosmological model is relatively successful and accepted by most physicists, it still has some unresolved problems, such as the fine-tuning problem and the origin problem of dark energy[3]. And for dark matter, people are still not sure what the component it is[4], since none of the particles described or predicted by the Standard particle Model meet the properties of dark matter. And more importantly, the dark matter have not been detected by experiments so far. It is in this complex context that various alternative models on Dark Energy and Dark Matter have been proposed. And among these alternative models to Dark Matter, the well-known and famous model maybe the MOdified Newtonian Dynamics (MOND)[5], which is first proposed by Milgrom in 1983. Milgrom postulates that at small accelerations the usual Newtonian dynamics break down and that the law of gravity needs to be modified. To date, MOND have been widely studied and it has been shown to successfully explain a large number of astrophysical observations, such as kinematics of galaxies[6,7], dynamics of wide binary stars[8,9], radial acceleration relation of galaxies[10,11], orbital velocity of interacting galaxy pairs[12], early-universe galaxy and cluster formation[13], and galaxy-scale gravitational lensing[14,15]. However, for alternative models to Dark Energy, they are still difficulty to distinguish, and so far no alternative model has stood out, but we notice an interesting model proposed in Ref. [16], which argues that the cosmic fluid is not perfect but rather viscous and dissipative, thus introducing a viscous pressure proportional to the Hubble constant into the standard cosmological model. Although the origin of the viscous pressure is still debated, it shows that the bulk viscosity model is significantly better than the ΛCDM model in fitting the combined SNe Ia + CMB + BAO + H(z) data.

In this paper, inspired by the idea in Ref. [16] that introducing a dissipation term into the ΛCDM model to describe the expansion of universe, driven by curiosity, we also try to introduce such a dissipation term into the rotational motion equation describing the disk galaxies, and then study its effect on the rotational motion of galaxies. Surprisingly, we found a high degree of coincidence between the derived motion properties and those predicted by MOND. However, as we have known that, MOND is only an empirical model, not a fundamental theory[14]. This seems to give us a window to glimpse the nature of MOND, or in other words, even if MOND is not a correct model, but as MOND can explain so many astronomical phenomena, we can also glimpse the secrets hidden behind MOND.

2. Newtonian Dynamics with Dissipation Term

In Ref. [16], the author introduced a dissipation term into the ΛCDM model to describe the expansion of universe, making the expansion and evolution equation describing the universe becomes

$$\stackrel{••}{a}=-{\displaystyle \frac{4G{\rho }_m}{3}}a+{\displaystyle \frac{\Lambda }{3}}a+12\pi G\zeta \stackrel{•}{a}$$

(1)

where a denotes the scale factor of universe, ρ_m is the total matter density, G is the gravitational constant, Λ denotes the cosmological constant, ζ denotes the introduced viscosity coefficient, which is positive.

In Ref. [16], it shows that Eq. (1) is significantly better than the ΛCDM model in fitting the astrophysical observation data. Inspired by this, we also try to introduce such a dissipation term into the rotational motion equation describing the disk galaxies and study what will happen, even if the physical basis of this dissipation term is unclear at present, and even if the physical spatial scales of the two objects of study are different.

Similar to Eq. (1), we also assume that there is a dissipation term in the rotational motion equation describing the disk galaxy that is proportional to the velocity, and that the corresponding viscosity coefficient is also positive just as which in Eq. (1). Then, without considering relativistic effects, we can easily obtain the force equilibrium equations in the radial and circumferential directions of rotational motion

$$\left\{\begin{array}{l}\stackrel{••}{r}+{\displaystyle \frac{GM}{r^2}}={\stackrel{•}{\theta }}^{2}r+\lambda \stackrel{•}{r}\\ {\displaystyle \frac{d}{dt}}(\stackrel{•}{\theta }r\cdot r)=\lambda \stackrel{•}{\theta }r\cdot r\end{array}\right.$$

(2)

where M denotes the center mass of the disk galaxy, (r, θ) denotes the coordinates of motion and here for simplicity we only consider the motion in a two-dimensional plane, λ is the assumed viscosity coefficient. And the second formula of Eq. (2) represents that the derivative of angular momentum corresponding to the rotational motion with respect to time is equal to the rotational torque.

One may notice that, in Eq. (2) it also implies an assumption that the dissipation force is related to the mass m of the object in rotating motion, because in its first formula both sides of the equation cancel out the mass m, which is also consistent with the idea in Ref. [16]. In fact, Eq. (2) is equivalent to introducing a dissipation potential ψ into the Newtonian dynamics

$$\psi =-{\displaystyle \frac{1}{2}}\lambda m[{\stackrel{•}{r}}^2+{(r\stackrel{•}{\theta })}^2]$$

(3)

Now anyway, we can use Eq. (2) to study the rotational motion of disk galaxies. First from the second formula of Eq. (2) we can obtain

$$\stackrel{•}{\theta }{r}^2=\stackrel{•}{{\theta }_0}{r}_0^2{e}^{\lambda t}$$

(4)

where (θ₀, r₀) is the initial condition. Then substituting Eq. (4) and GM=v₀²r₀ into the first formula of Eq. (2)

$$\stackrel{••}{R}=\lambda \stackrel{•}{R}+{\displaystyle \frac{v_0^2}{r_0^2}}{e}^{2\lambda t}{\displaystyle \frac{1}{R^3}}-{\displaystyle \frac{v_0^2}{r_0^2}}{\displaystyle \frac{1}{R^2}}$$

(5)

where R=r/r₀ is a dimensionless quantity, and $v_0=\stackrel{•}{{\theta }_0}{r}_0$. Here it should be noted that, in the derivation of Eq. (5), we used the equation GM=v₀²r₀, which assumes that the object in rotating motion is in a strong gravitational field at initial, thus leading to Eq. (2) returns to the usual Newtonian dynamics.

The analytical solution of Eq. (5) is difficult to obtain, but we can obtain its numerical solution. Since there is currently no clear physical basis to help us to determine the value of λ (from the success of Newtonian dynamics, we know that λ must be an extreme small value), let us discuss the motion properties predicted by Eq. (2) by arbitrarily setting the value of λ and the initial condition of the equation. And similar to MOND, here, the motion properties we are interested in is also the relationship between the gravitational acceleration and the centripetal acceleration subjected by the object in rotating motion. And it is well known that in MOND, the relationship between the gravitational acceleration and the centripetal acceleration was assumed as

$$g_N=a\mu (a/a_0)$$

(6)

where g_N denotes the gravitational acceleration, a denotes the centripetal acceleration, μ is an empirical or undetermined interpolation function, and when a/a₀≫1, μ→1, when a/a₀≪1, μ→a/a₀, and a₀≈1.2×10^-10 m/s² is a constant.

As we know, in MOND, when the gravitational field is very weak, i.e., a/a₀≪1, there is v_u⁴=GMa₀, where v_u denotes the circumferential velocity of the object in rotating motion, which can explain the asymptotically flat rotational velocity curve of disk galaxies. Coincidentally, we also find that based on Eq. (5), no matter what the initial conditions are, as long as time evolves long enough, although the radial distance of the object in rotating motion is increasing over time, the circumferential velocity remains unchanged (here we assume the “unchanged” velocity is v_ultimate, and in the subsequent analysis, we will find that from Eq. (5) it can derive all the formulas in MOND, which in turn proves this conclusion). Thus, similar to MOND, here we define a physical quantity A₀, and it satisfies

$$v_{ultimate}^4=GM{A}_0$$

(7)

where A₀, which just to distinguish from a₀ in MOND, has the same meaning as a₀.

So next we will discuss Eq. (2) or Eq. (5) in terms of the formula form in MOND, that is

$$g_N=a\mu \text{'}(a/A_0)$$

(8)

Now we will conduct numerical analysis on Eq. (5), and Table 1 shows the value of A₀ and μ’ corresponding to different λ and initial conditions. Here we did not set units in the data in Table 1, because it is not necessary. Whether there are units in the numerical analysis does not affect properties or performance of Eq. (5).

In Table 1, the number group (1,6,7,8), (9,10,11,12), (13,14,15,16) and (19,20,21) corresponds to the same mass M, that is $GM=v_0^2{r}_0$. And from these number group we can see that A₀ is only dependent on M and λ, not on the initial conditions (v₀, r₀), as long as the initial centripetal acceleration corresponding to the initial conditions is strong enough, which is an important step for Eq. (2) yields the MOND. One might wonder why does it here say that as long as the initial centripetal acceleration is strong enough, it’s because if the initial centripetal acceleration a’₀ is too small and λ is relatively large, v_ultimate will be close to v₀, and therefore A₀≈a’₀ according to Eq. (7), which means that A₀ is dependent on M or the initial conditions (v₀, r₀), almost not on λ. In other words, if a large number of objects rotating around the center of the disk galaxy are initially in the strong gravitational potential region of the galaxy, and we assume that for this galaxy the value of λ is a constant, then, regardless of which orbits these objects are initially distributed in (i.e., their orbits are different), they will all eventually reach the same rotational circumferential velocity under the action of dissipation potential (i.e. Eq. (3)), and the fourth power of this ultimate rotational circumferential velocity (i.e., v_ultimate) is proportional to the central mass of the galaxy (i.e. Eq. (7), which also known as the famous Tully-Fisher relation [17,18,19]).

In Table 1, from the number group (1,2,3,4,5) we can see that if we assume the value of λ is a constant independent of the central mass of the disk galaxy, then, the larger the mass of a galaxy, the higher the value of A₀, which seems to echo to a model called EMOND (i.e., extend MOND, which assumes that a₀ in MOND is increased within the galaxy cluster, since that the MOND paradigm has achieved success on the galaxy scale, but it is unable to explain the galaxy clusters) [20], and we intend to further discuss this conclusion in the latter section.

In Table 1, from the number group (1,9,13) and (6,14,19) we can see that for the same galaxy, the higher the assumed value of λ, the higher the value of A₀.

At last, the most important conclusion from Table 1 and Figure 1 is that, as long as the initial centripetal acceleration corresponding to the initial conditions is strong enough, whatever the assumed value of λ is, whichever the initial orbit of the object rotating around the center of the disk galaxy is, or whatever the center mass of the disk galaxy is, μ’, the interpolation function, are or are almost the same. And from Figure 1 we can see that, when x≪1, μ’→x, and when x≫1, μ’→1, which is exactly what is assumed in the MOND’s empirical formula. Further, we find that the curve in Figure 1 can be fitted as the following function

$$\mu \text{'}(x)=0.1455{\displaystyle \frac{\sqrt{1+4x}-1}{\sqrt{1+4x}+1}}-0.09382{\displaystyle \frac{x}{\sqrt{1+x^2}}}-0.2126(1-e^{-x})+1.16092{\displaystyle \frac{x}{1+x}}$$

(9)

Unfortunately, we here cannot give a clear mathematical definition to the above statement “as long as the initial centripetal acceleration is strong enough”, mainly because for Eq. (5) we cannot obtain its analytical solution, but its numerical solution.

In above, we only introduced a viscosity coefficient λ into the Newtonian dynamics and then we have almost completely derived all the formulas assumed in MOND, especially the concrete interpolation function. However, we also obtained something different from MOND, that is, A₀ does not need to be a constant, and the disk galaxies are also gradually expanding over time at an extreme small expansion rate, just like universe is expanding.

Although the above conclusions are drawn based on the data in Table 1, the amount of data in Table 1 is far from enough. In fact, we privately expanded the amount of data in Table 1, and found that the obtained conclusions remain the same, so we did not list these lengthy data. In addition, it should be stressed again that the absence of units for the data in Table 1 does not affect the above obtained conclusions from Eq. (5), or in other words, in Table 1, whether the unit of λ is 1/s or 1/h or 1/year (the possible value of λ will be discussed in Sect. 5), and whether the unit of v₀ is m/s or km/s, and whether the unit of r₀ is m or parsec, Eq. (9) remains the same, and the relationship between A₀ and M remains the same.

3. Discussion on a₀ from Observations

MOND, as a competitive model against ΛCDM, was initially proposed with the aim of explaining the abnormal rotational velocity curves of galaxies. For MOND, the usual method is to fix the parameter a₀ and take M/L (i.e., Mass-to-light ratio) as a free parameter, then based on the observational data, an optimal fitting parameter M/L can be obtained. If we check the obtained parameter M/L is within a reasonable range of the galaxy, then it indicates that MOND can fit the rotation curve of the galaxy well. So far, many studies using the above method have shown that MOND can fit the rotation curves of many galaxies very well (for instance, in Ref. [21] the author states that no single model can fit all galaxies, but MOND has successfully fitted the largest number of galaxies). However, some other authors also attempted to fix M/L and take a₀ as a free parameter to fit the observational data of galaxies, regardless of a₀ being regarded as a universal constant in MOND. For example, in Ref. [22], Randriamampandry and Carignan studied the mass models of 15 nearby dwarf and spiral galaxies and found that galaxies with higher central surface brightness tend to favour higher values of the constant a₀, and in Ref. [23], Swaters et al. studied a sample of 27 dwarf and low surface brightness galaxies and found that lower surface brightness galaxies tend to have lower a₀. The two conclusions just coincide with the conclusion we obtained in Sect. 2, i.e., the larger the mass of a galaxy, the higher the value of A₀ (remember that A₀ and a₀ in the equation have the same meaning).

From the analysis of a large amount of literature, there has always been controversy over whether a₀ is a universal constant in MOND. On the one hand, some literature show that, fixing a₀ but allowing M/L to be a free parameter to fit the rotation curve of galaxies can indeed fit some observational data very well. On the other hand, when different galaxy samples are used to analyze whether the a₀ credible intervals for individual galaxies is compatible among these samples, some literature concluded that a₀ is not a universal constant (for instance, see Ref. [24]), and some other literature take a more aggressive approach, that is, taking a₀ as a free parameter to fit the observational data, and it is found that a₀ is related to the surface brightness of galaxies (for instance, the above mentioned Ref. [22] and Ref. [23]). The main reason for such controversy is that, on the one hand, MOND is not a fundamental theory but more of an empirical model, on the other hand, the exclusion of dark matter (if it exists) in the MOND model limits its ability to fit the observational data.

In this paper, we provide another new perspective, that is, the derived a₀ is not a constant but vary with the mass of galaxy, and the existence of dark matter is also permitted in the model (which can be used to adjust the value of a₀, since that the value of a₀ is related to the mass of galaxy or the gravitational potential). Thus, a model that allows dark matter to exist and a₀ is variable emerges. In the next section, we will see the advantages of this model over MOND in some aspects.

In addition, we know that the MOND model is also under development, such as the EMOND (i.e., Extend MOND), which is first proposed by Zhao and Famaey [20]. This model assumes that a₀ is related to the gravitational potential, and the stronger the gravitational potential, the larger the value of a₀, thereby compensating for the drawback that MOND cannot explain the galaxy clusters. The ideological of this model is consistent with the conclusion we have drawn above, although there are differences in mathematical details. We intend to discuss them in the next paper.

4. Models Comparison

To compare observation data with ΛCDM, MOND and Eq. (8), the Newtonian gravitational acceleration due to the ordinary matter (i.e., not including dark matter) must be separated from the total acceleration. Here we assume the acceleration due to baryons comes from the bulge and the disk (with no dark halo contribution), that is

$$a_B(r)={\displaystyle \frac{v_b^2(r)+v_d^2(r)}{r}}$$

(10)

where r is the distance from the center of the galaxy, v_b and v_d are the rotational velocity due to the bulge and disk respectively.

In ΛCDM, in addition to the ordinary matter, dark matter is considered, and the function of dark halo is usually expressed as [25]

$${\rho }_h(r)={\displaystyle \frac{{\rho }_0}{\frac{r}{h}(1+\frac{r^2}{h^2})}}$$

(11)

where ρ₀ and h are the scale density and scale radius of the dark halo respectively. Then the mass of dark halo can be expressed as

$$M_h(r)=4\pi {\rho }_0[\ln (1+{\displaystyle \frac{r}{h}})-{\displaystyle \frac{r}{r+h}}]$$

(12)

Then it leads to

$$v_h^2={\displaystyle \frac{G{M}_h(r)}{r}}$$

(13)

The total acceleration is then

$$g_t={\displaystyle \frac{v_b^2+v_d^2+v_h^2}{r}}$$

(14)

In MOND, there is no dark matter, and the relationship between gravitational acceleration and centripetal acceleration is expressed as Eq. (6). Since the interpolation function is empirical, here we take the “simple” function [26]

$$\mu (x)={\displaystyle \frac{x}{1+x}}$$

(15)

Then we can plot the centripetal acceleration vs. the gravitational acceleration due to baryons for the observation data in Figure 2 [27], and the fitting curves corresponding to various models are also plotted.

It should be note that the prediction in Figure 2 assumes that the galaxies being observed are isolated. Figure 2 shows that for smaller radii, i.e. large accelerations, MOND and Eq. (8) do not significant deviate from expected ΛCDM model galaxies, but that changes for lower accelerations. In addition, it indicates that the centripetal accelerations in M31 and Milky Way are smaller than the MOND prediction at lower a_B, but in line with the prediction of Eq. (8). This is mainly because there are no free parameters in the prediction of radial acceleration relation (RAR) of galaxies in MOND, while in Eq. (8), A₀ is a parameter that can vary with the mass of the galaxy, so we can adjust A₀ to fit the observation data in Figure 2.

One may wonder whether dark matter was taken into account in Eq.(8) when fitting the observation data in Figure 2. The answer is that whether you consider it or not makes no mathematical difference to the fitting results. When dark matter exists, we assume that dark matter exists near the center of the disk galaxy (similar to MOND, our model does not need the dark matter halo to be carefully adjusted to meet the galaxy’s asymptotically flat rotational velocity curve), which means that the coordinate a_B in Figure 2 should be replaced by g_t. Since dark matter exists near the center of the galaxy, and the radial position corresponding to the maximum value of a_B in Figure 2 is far away from the center of the galaxy, it is easy to prove that the g_t-axis can be obtained by adding a same value (which is determined by the mass of dark matter) to all values of the a_B-axis, or in other words, the whole shift of the a_B-axis to the right side can yield the g_t-axis. In Eq. (8), the left side is just related to g_t, and the right side is just related to a. In addition, from Table 1 it has been pointed out that the function form of μ’ has nothing to do with the central mass of the galaxy (that is, μ’ is independent of a_B and g_t), so it can be easily proved that whether there exists dark matter near the center of the galaxy does not affect the fitting result of Eq. (8) mathematically. But physically, dark matter may be needed as its mass can help to adjust the value of A₀.

5. Discussion on λ

Obviously, the value of λ assumed above is crucial to the evolution of old galaxies (but not so important for young ones). Then, an important question is, what is the value of λ? Or, as a second choice, what order of magnitude might the value of λ be. As there is currently no clear physical basis for λ, so at present the specific value of λ is difficult to obtain. But we can roughly obtain the order of magnitude of λ from the rotation curve of the galaxy.

From Eq. (4) it has vr=v₀r₀e^λt, then when v=v₀, there is r=r₀e^λt, which indicates that the radial coordinates of a rotating object will continuously increase over time (but its circumferential rotational velocity remains the same). And we know that the rotation curves of many galaxies are asymptotically flat, as shown in Figure 3, and from the viewpoint of Eq. (2) or Eq. (5), in Figure 3, when point 2 and point 1 are on the flat curve, point 2 is the evolution of point 1 over time. Therefore, we can calculate λ based on the radial coordinate ratio of point 1 to point 2 (i.e., r₁/r₀) as well as the time required for point 1 to evolve to point 2.

In Ref. [28], many rotational velocity curves of galaxies are listed (all of which show the characteristic of asymptotically flat), from these rotational velocity curves we can find that 1<r₁/r₀≤4, however, the time for evolution from point 1 to point 2 is hard to obtain, which is not only related to the mass of the galaxy, but also to the initial orbit of the object and the value of λ. Here, we take the timescale of the universe’s evolution to roughly estimate the order of magnitude of λ, that is, we assume t=10Gyr, then we can obtain that the order of magnitude of λ is about 10^-18/s, which coincidentally is the same as the order of magnitude of the present Hubble constant. It reminds us that, with such a small value of λ, the radial expansion of the galaxy itself caused by it is also extremely small, which brings great difficulties to the detection.

Another question, one may wonder whether λ must be a constant. As has been stated that, there is currently no clear physical basis for λ (Perhaps it originated from an unknown but deeper theory), so we don’t know the answer at present, however, if the rotational velocity curve of the galaxy does eventually tend to be “absolutely” flat, then λ must be a constant.

6. Summary

This paper, with no clear physical basis, just follows the idea in Ref. [16] with curious, that is, introducing a viscosity coefficient into the ΛCDM model to describe the expansion of universe, also introduces such a similar positive viscosity coefficient into the rotational motion describing the disk galaxies, and then studies what will happen. We were surprised to find that it yields all the formulas assumed in MOND, including a concrete interpolation function between the centripetal acceleration and the Newtonian acceleration, which however is empirical in MOND.

But we also obtain something different from MOND. First, a₀ is considered to be a constant in MOND, but we obtained that A₀ (here it just to distinguish from a₀ in MOND, and A₀ and a₀ in the equation have the same meaning) is not a constant but is related to the mass of the galaxy, and the more massive the galaxy, the greater the value of A₀ if the value of λ is a constant. This characteristic can help it to fit the RAR data better than MOND. Second, the paper predicts that the rotational galaxies will expand radially over time at an extreme small expansion rate, just like the expansion of universe, so it predicts a “dynamical” galaxy rather than a “static” galaxy predicted by ΛCDM and MOND.

In many literature, such as Refs. [29,30,31,32], it is found that MOND can fit the observation data of rotational galaxies better than ΛCDM, which is thought to be mainly attributed to the two extreme conditions (i.e., Newtonian and deep-MOND limits) in MOND, and has no or little relationship with the concrete form of the interpolation function in MOND. For these observation data, it is easy to know that the model in this paper can also fit them well as the two extreme conditions can also be naturally derived from our model rather than be assumed in MOND. In MOND, in order to fit these observational data, the mass-to-light ratio of the galaxy, which is the only free parameter in MOND [33], have to be adjusted, and therefore the existence of dark matter must be excluded, because the presence of dark matter can affect the mass-to-light ratio. However, the model in this paper cannot rule out the existence of dark matter, one reason is that the mass of dark matter can help to adjust the value of A₀ to better fit the RAR curves, and another reason is that the idea of this paper comes from Ref. [16], which supports the existence of dark matter. However, unlike ΛCDM, in our model, the asymptotically flat rotational velocity curve of the disk galaxy does not require the careful adjustment of dark matter (if it exists) to be achieved, which seems unnatural. It can be also achieved naturally under the viscous dynamics of the galaxy itself. So, in fact, we tend to consider that the idea of this paper to be a further optimization of ΛCDM rather than a replacement.