Exact Semianalytical Calculation of Rotation Curves with Bekenstein–Milgrom Nonrelativistic MOND

In Newtonian gravity, a density distribution ρ( r ) produces a field of accelerations

$$\begin{eqnarray}&&{{\boldsymbol{g}}}_{{\rm{N}}}[\rho ]({\boldsymbol{r}})=G\int d{\boldsymbol{r}}^{\prime} \displaystyle \frac{\rho ({\boldsymbol{r}}^{\prime} )}{| {\boldsymbol{r}}^{\prime} -{\boldsymbol{r}}{| }^{3}}({\boldsymbol{r}}^{\prime} -{\boldsymbol{r}}).\end{eqnarray} \tag{ 1 }$$

Using cylindrical coordinates in which R and z are radial and vertical distance, respectively, and ϕ is the azimuthal angle, in an axisymmetric matter distribution (ρ ≠ ρ(ϕ)), the rotation speed V_c for stars in equilibrium with the centrifugal force with Newtonian gravitation is (Chrobáková et al. 2020, Appendix A)

$$\begin{eqnarray}\begin{array}{rcl}{g}_{{\rm{N}},{\rm{R}}}(R,z) & = & \displaystyle \frac{-{V}_{c}^{2}}{R}=\displaystyle \frac{-2\,G}{R}{\displaystyle \int }_{-\infty }^{\infty }{dz}^{\prime} \\ & & \times {\displaystyle \int }_{0}^{\infty }{dR}^{\prime} \,R^{\prime} \rho (R^{\prime} ,z^{\prime} )\end{array}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray*}&&\begin{array}{l}\times [{C}_{E}\,E(k)+{C}_{K}\,K(k)]\\ {C}_{E}=\displaystyle \frac{(R^{\prime} +R)(R^{\prime} -R)+{\left(z-z^{\prime} \right)}^{2}}{[{\left(R^{\prime} -R\right)}^{2}+{\left(z-z^{\prime} \right)}^{2}]\sqrt{{\left(R^{\prime} +R\right)}^{2}+{\left(z-z^{\prime} \right)}^{2}}}\\ {C}_{K}=-\displaystyle \frac{1}{\sqrt{{\left(R^{\prime} +R\right)}^{2}+{\left(z-z^{\prime} \right)}^{2}}}\\ k=\sqrt{\displaystyle \frac{4\,R\,R^{\prime} }{{\left(R^{\prime} +R\right)}^{2}+{\left(z-z^{\prime} \right)}^{2}}},\end{array}\end{eqnarray*}$$

where K(k) and E(k) are the complete elliptical integrals of the first and second kind, respectively, and g_N,R is the Newtonian radial acceleration.

The vertical acceleration with Newtonian gravity is

$$\begin{eqnarray}&&{g}_{{\rm{N}},{\rm{z}}}(R,z)=\displaystyle \frac{-G\,z}{2\,{R}^{3/2}}{\int }_{-\infty }^{\infty }{dz}^{\prime} {\int }_{0}^{\infty }{dR}^{\prime} \displaystyle \frac{\rho (R^{\prime} ,z^{\prime} )}{\sqrt{R^{\prime} }}H(k)\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray*}&&H(k)={\int }_{0}^{2\pi }d\phi ^{\prime} \displaystyle \frac{1}{{\left[(2/{k}^{2})-\cos \phi ^{\prime} \right]}^{3/2}}\end{eqnarray*}$$

$$\begin{eqnarray*}&&k=\sqrt{\displaystyle \frac{4\,R\,R^{\prime} }{{\left(R^{\prime} +R\right)}^{2}+{\left(z-z^{\prime} \right)}^{2}}}.\end{eqnarray*}$$

The general algebraic rule of Milgrom's empirical law usually applied by astronomers to calculate rotation curves with MOND is as follows (Famaey & McGaugh 2012, Equation (7)):

$$\begin{eqnarray}&&{{\boldsymbol{g}}}_{{\rm{M}}-{\rm{R}}}=\displaystyle \frac{{{\boldsymbol{g}}}_{{\rm{N}}}}{\mu (| {{\boldsymbol{g}}}_{{\rm{M}}-{\rm{R}}}| /{a}_{0})},\end{eqnarray} \tag{ 4 }$$

where μ(x) is an interpolating function. The standard interpolating function is

$$\begin{eqnarray}&&\mu (x)=\sqrt{\displaystyle \frac{1}{1+{x}^{-2}}}.\end{eqnarray} \tag{ 5 }$$

With the last two equations,

$$\begin{eqnarray}&&| {g}_{{\rm{M}}-{\rm{R}}}| =\sqrt{\displaystyle \frac{1}{2}{g}_{{\rm{N}}}^{2}+\sqrt{\displaystyle \frac{1}{4}{g}_{{\rm{N}}}^{4}+{g}_{{\rm{N}}}^{2}{a}_{0}^{2}}},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray*}&&{{\boldsymbol{g}}}_{{\rm{M}}-{\rm{R}}}={{\boldsymbol{g}}}_{{\rm{N}}}\displaystyle \frac{| {g}_{{\rm{M}}-{\rm{R}}}| }{| {g}_{{\rm{N}}}| },\end{eqnarray*}$$

as used, for instance, by Bottema & Pestaña (2015) and Lisanti et al. (2019) for rotation curve fits. The rotation speed in the plane would be ${V}_{c}=\sqrt{| {{\boldsymbol{g}}}_{{\rm{M}}-{\rm{R}}}| \,R}$.

Equation (4) leads to some consistency problems. For instance, in a two-body case, as the implied force is not symmetric in the two masses, Newton's third law does not hold, so the momentum is not conserved (Milgrom 1983a; Felten 1984; Famaey & McGaugh 2012). Precisely because of this, more sophisticated formulations of MOND were created as a modification of classical dynamics: the AQUAL model by Bekenstein & Milgrom (1984) that we will use here.

An unpleasant characteristic of MOND is the nonlinear nature of its equations. However, nonlinear field equations can be expressed as linear field equations containing a self-source term. This may be convenient for heuristic reasons and because in this manner the equations can be solved by well-known linear methods. Obviously, this does not mean that the nonlinear difficulties can be avoided, because, due to the presence of the self-source term, the equations must be solved self-consistently. However, as we shall show, this solution may be obtained with considerable precision in a few iterations.

The Poisson equation of the AQUAL solution for the potential ϕ is

$$\begin{eqnarray}&&{\rm{\nabla }}(\mu (x){\rm{\nabla }}\phi )=4\pi G\rho ,\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray*}&&x=\displaystyle \frac{| {{\boldsymbol{g}}}_{{\rm{M}}-{\rm{E}}}| ({\boldsymbol{r}})}{{a}_{0}},\end{eqnarray*}$$

where g _M−E( r ) is the exact MOND acceleration. The left-hand side of this equation may be written in the form μ(x)∇² ϕ + (∇μ(x))(∇ϕ). Since

$$\begin{eqnarray}&&{\rm{\nabla }}\mu (x)=\displaystyle \frac{1}{{a}_{0}}\displaystyle \frac{d\mu }{{dx}}(x).{\rm{\nabla }}| {{\boldsymbol{g}}}_{{\rm{M}}-{\rm{E}}}| \end{eqnarray} \tag{ 8 }$$

and

$$\begin{eqnarray}&&{{\boldsymbol{g}}}_{{\rm{M}}-{\rm{E}}}=-{\rm{\nabla }}\phi ,\end{eqnarray} \tag{ 9 }$$

we get

$$\begin{eqnarray}&&\mu (x){{\rm{\nabla }}}^{2}\phi =4\pi G\rho +\displaystyle \frac{1}{{a}_{0}}\displaystyle \frac{d\mu }{{dx}}(x){\rm{\nabla }}| {{\boldsymbol{g}}}_{{\rm{M}}-{\rm{E}}}| \,\cdot {{\boldsymbol{g}}}_{{\rm{M}}-{\rm{E}}}.\end{eqnarray} \tag{ 10 }$$

The second term in this equation is the self-source term, which can be interpreted as a density of phantom matter. Therefore, the solution is

$$\begin{eqnarray}&&\phi ({\boldsymbol{r}})=\int d{\boldsymbol{r}}^{\prime} \displaystyle \frac{G\rho ({\boldsymbol{r}}^{\prime} )+\tfrac{1}{4\pi {a}_{0}}\tfrac{d\mu }{{dx}}({\boldsymbol{r}}^{\prime} ){\rm{\nabla }}| {{\boldsymbol{g}}}_{{\rm{M}}-{\rm{E}}}({\boldsymbol{r}}^{\prime} )| \,\cdot {{\boldsymbol{g}}}_{{\rm{M}}-{\rm{E}}}({\boldsymbol{r}}^{\prime} )}{\mu ({\boldsymbol{r}}^{\prime} )| {\boldsymbol{r}}-{\boldsymbol{r}}^{\prime} | }.\end{eqnarray} \tag{ 11 }$$

For a thin disk, the exact calculation of MOND acceleration g _M−E( r ) is equivalent to the Newtonian gravity calculation (equations given in Section 2.1), but setting as density

$$\begin{eqnarray}&&{\rho }^{* }({\boldsymbol{r}})=\displaystyle \frac{\rho ({\boldsymbol{r}})}{\mu (x)}+\displaystyle \frac{1}{4\pi \,G\,{a}_{0}\,\mu (x)}\displaystyle \frac{d\mu }{{dx}}(x)\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray*}&&\times \,\left|{g}_{{\rm{M}}-{\rm{E}},{\rm{R}}}({\boldsymbol{r}})\displaystyle \frac{\partial {g}_{{\rm{M}}-{\rm{E}}}({\boldsymbol{r}})}{\partial R}+{g}_{{\rm{M}}-{\rm{E}},{\rm{z}}}({\boldsymbol{r}})\displaystyle \frac{\partial {g}_{{\rm{M}}-{\rm{E}}}({\boldsymbol{r}})}{\partial z}\right|.\end{eqnarray*}$$

This equation was also derived by Milgrom (1986, Equation (4)). That is, we can calculate the acceleration iterating between Equation (1) g _M−E[ρ]( r ) = g _N[ρ^∗]( r ) and Equation (12). For the first iteration, we set g _M−E = g _M−R. In practice, when g_M−E is close to g_M−R, we only need two iterations (considering the first iteration ρ* = ρ), since g _M−E[ρ]( r ) ≈ g _N[ρ^∗[ g _M−R]]( r ) is a good approximation. We will show in Section 3.4 that even among cases with very large differences between g_M−E and g_M−R, the iterative process converges quickly in three to four iterations.

The rotation speed in the plane would again be ${V}_{c}=\sqrt{| {{\boldsymbol{g}}}_{{\rm{M}}-{\rm{E}}}| \,R}$.

By construction, g _M−E[ρ]( r ) = g _M−R[ρ]( r ) in cases of spherical symmetry in the density distribution. We will test it in the simplest case, a point-like mass in the center of the galaxy:

$$\begin{eqnarray}\sigma (R)=\left\{\begin{array}{ll}M\delta (R), & R=0\\ 0, & R\ \gt \ 0\end{array}\right\}.\end{eqnarray} \tag{ 15 }$$

In Figure 1, we show this perfect agreement (the very slight relative differences lower than 1% are due to numerical calculation errors) of both MOND algorithms in counterdistinction of Newtonian gravity.

Zoom In Zoom Out Reset image size

The most usual fit of the galactic thin disk is with a simple exponential law:

$$\begin{eqnarray}&&\sigma (R)=\displaystyle \frac{M}{2\pi {H}^{2}}\exp \left(-\displaystyle \frac{R}{H}\right),\end{eqnarray} \tag{ 16 }$$

where M is the total mass of the disk and H is its scale length. In Figures 2 and 3, we show the rotation curves for a total mass of 10¹¹ M_⊙ and scale lengths of 3 and 6 kpc, for Newtonian gravity ( g _N), MOND-rule ( g _M−R), and MOND-exact ( g _M−E), respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A distribution characterized by giving a flat rotation curve in Newtonian gravity without any extra component, such as a dark matter halo, is the Mestel disk, whose dependence with the galactocentric radius goes as (Schulz 2012)

$$\begin{eqnarray}\sigma (R)=\left\{\begin{array}{ll}\displaystyle \frac{M}{2\pi {R}_{\max }R}\arccos \left(\displaystyle \frac{R}{{R}_{\max }}\right), & R\leqslant \ {R}_{\max }\\ 0, & R\gt {R}_{\max }\end{array}\right\}.\end{eqnarray} \tag{ 17 }$$

${R}_{\max }$ is the radius of the disk. For larger radii, there is not any mass. In Figure 4, we show the rotation curve for a total mass of 10¹¹ M_⊙ and maximum radius of 20 kpc, for Newtonian gravity ( g _N), MOND-rule ( g _M−R), and MOND-exact ( g _M−E), respectively.

Zoom In Zoom Out Reset image size

The Maclaurin disk, a limiting case of the Maclaurin spheroid, is applicable to the study of spiral galaxies for which the rotational velocity rises to the edge of the disk as is seen in irregular galaxies (Schulz 2009). It follows

$$\begin{eqnarray}\sigma (R)=\left\{\begin{array}{ll}\displaystyle \frac{3\,M}{2\pi \,{R}_{\max }^{2}}\sqrt{1-\left(\displaystyle \frac{R}{{R}_{\max }}\right)}, & R\leqslant \ {R}_{\max }\\ 0, & R\gt {R}_{\max }\end{array}\right\}.\end{eqnarray} \tag{ 18 }$$

${R}_{\max }$ is the radius of the disk. It is an almost constant density for $R\ll {R}_{\max }$, and it declines fast to zero for $R\lesssim {R}_{\max }$. In Figure 5, we show the rotation curve for a total mass of 10¹¹ M_⊙ and maximum radius of 20 kpc, for Newtonian gravity ( g _N), MOND-rule ( g _M−R), and MOND-exact ( g _M−E) respectively, only with two iterations of Equations (1) and (12).

Zoom In Zoom Out Reset image size

This is a case with strong differences between g _M−R and g _M−E, so we explore higher number iterations, to see that within iteration 3 or 4 it converges at R ≳ 5 kpc; see Figure 6. We see that the corrections of a higher iteration than 2 are of second order even in this case of differences between g _M−R and g _M−E, which we will not take into account from now on.

Zoom In Zoom Out Reset image size

We have seen in the previous subsections that the MOND-rule may be very inexact with respect to MOND-exact in cases of thin disks where most of the mass is concentrated in the outer parts. The most extreme case among those we have tested is the Maclaurin disk.

In all of the previous cases, we have set a total mass disk of M = 10¹¹ M_⊙. Let us define R_* as the radius where the MOND regime is significant at more than 10% in rotation speeds: $\left[\tfrac{| {V}_{{\rm{c}},\mathrm{MOND}-\mathrm{rule}}({R}_{* })-{V}_{{\rm{c}},\mathrm{Newton}}({R}_{* })| }{{V}_{{\rm{c}},\mathrm{Newton}}({R}_{* })}\right]=0.1$. We also define $r({R}_{* })\equiv \tfrac{M(R\gt {R}_{* })}{M}$ (the ratio of mass at distances larger than R_*, where MOND is significant) and d(R_*) as the maximum ∀ R > R_* of the ratio $\left[\tfrac{| {V}_{{\rm{c}},\mathrm{MOND}-\mathrm{exact}}(R)-{V}_{{\rm{c}},\mathrm{MOND}-\mathrm{rule}}(R)| }{{V}_{{\rm{c}},\mathrm{MOND}-\mathrm{rule}}(R)}\right]$. In Figures 1–5, the values of R_* are 9.8, 10.6, 0, 7.6, and 0 kpc, and r(R_*) for the five models used are 0, 0.139, 1, 0.477, and 1, respectively, whereas their respective values of d(R_*) are 0, 0.0390, 0.0592, 0.0474, and >0.70. Roughly, we can see that the highest values of d(R_*) are obtained for the values with highest mass in the external parts within the MOND regime. The error with respect to the algebraic solution is larger than 5% when more than half of the mass is in the MONDian regime, reaching in some cases >70% discrepance.

Most spiral galaxies present an approximately flat rotation curve in the outer part, although there may be significant deviations from that behavior. A rising rotation curve is observed in the Andromeda galaxy, which was claimed to be challenging for a model with standard dark matter models or perturbations of the galactic disk by satellites (Ruiz-Granados et al. 2010). However, we see here that MOND with a Mestel or Maclaurin disk naturally gives this increase, and Newtonian for the Maclaurin disk alone; adding a dark matter halo with a strong distribution of mass in the outer parts would reinforce this trend.

A decrease of rotation curve instead of flat curve (e.g., Eilers et al. 2019; Zobnina & Zasov 2020) in the very outer disk might be more surprising for MOND, since one expects an asymptotic limit ${\mathrm{lim}}_{r\to \infty }{{\boldsymbol{g}}}_{M-R}=\tfrac{\sqrt{G\,M\,{a}_{0}}}{r}$ and ${V}_{c}=\sqrt{| {{\boldsymbol{g}}}_{M-R}| r}$. In our analyses of the different disk models, we see for the MOND-exact solution that the slope is always more negative than in the MOND-rule approximation. For instance, at R = 50 kpc we get that ${\left(\tfrac{{{dV}}_{c}}{{dR}}\right)}_{\mathrm{MOND}-\mathrm{Exact}}$ is −0.197, −0.397, −0.212, and −0.400 km s⁻¹ kpc⁻¹ for Figures 2–5, respectively, whereas ${\left(\tfrac{{{dV}}_{c}}{{dR}}\right)}_{\mathrm{MOND}-\mathrm{Rule}}$ values are −0.087, −0.254, −0.111, and −0.165 km s⁻¹ kpc⁻¹. That is, the negative slope at R = 50 kpc is multiplied in the exact solution by a factor of 1.6–2.4 with respect to the approximate rule calculation. Therefore, this factor is important and should be taken into account in a discussion about MOND plausibility in some decreasing rotation speeds.

Another caveat in the fit of rotation curves in MOND stems from the dependence of the amplitude of the rotation curve on the distance from the plane (z). In the Milky Way, Jalocha et al. (2010) observe that this dependence is strong and favors a Galaxy without a dark matter halo; however, the Jeans equation to convert azimuthal velocities into rotation speed was not fully considered. A more recent analysis by Chrobáková et al. (2020) with Gaia data, taking into account the dispersion of velocities with the Jeans equation, gives a mild or negligible dependence with z, implying that a spherical component dominates. Let us calculate here, in our exponential disk examples, the dependence on z of the rotation speed at R = 20 kpc; the results are in Figures 7 and 8. We see that, for the lowest value of H, the MOND-exact solution gives a flat dependence of ∣z∣, whereas the MOND-rule approximation predicts a slight fall-off. In a sense, the phantom mass of MOND behaves more like a spherical distribution similar to the halo, and this is something to be considered in the evaluation of the suitability of MOND to fit rotation curves away from the plane.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the example of a central point-like mass, as introduced in Section 3.1, Equation (12) gives a total density ρ*(r) with spherical symmetry (only dependent on the distance to the center r) equal to

$$\begin{eqnarray}\begin{array}{rcl}\rho * (r) & = & \rho (r)+\displaystyle \frac{1}{4\pi \,G\,{a}_{0}\,\mu (x)}\displaystyle \frac{d\mu }{{dx}}(x)\\ & & \times \left|{g}_{{\rm{M}}-{\rm{E}}}(r)\displaystyle \frac{{{dg}}_{{\rm{M}}-{\rm{E}}}(r)}{{dr}}\right|,\end{array}\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray*}&&x=\displaystyle \frac{| {g}_{{\rm{M}}-{\rm{E}}}| (r)}{{a}_{0}},\end{eqnarray*}$$

$$\begin{eqnarray*}&&{g}_{{\rm{M}}-{\rm{E}}}(r)={g}_{{\rm{M}}-{\rm{R}}}(r)=\displaystyle \frac{G\,M}{{r}^{2}\mu (x)}.\end{eqnarray*}$$

Developing this equation, and introducing Equation (5) into it, we get

$$\begin{eqnarray}&&{\rho }^{* }(r)=\rho (r)+\displaystyle \frac{M}{2\pi {r}^{3}}\displaystyle \frac{\sqrt{1+{x}^{-2}}}{2+{x}^{2}},\end{eqnarray} \tag{ 20 }$$

with ρ = 0 if r ≠ 0. In the limit of infinite distance,

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{r\to \infty }\rho * (r)=\mathop{\mathrm{lim}}\limits_{x\to 0}{\rho }^{* }(r)=\displaystyle \frac{1}{2\pi }\sqrt{\displaystyle \frac{M\,{a}_{0}}{G}}\displaystyle \frac{1}{{r}^{2}},\end{eqnarray} \tag{ 21 }$$

that is, it falls down inversely proportional to r². This is close to the decrease as r^{−1.76±0.12} obtained as the best halo fit within Newtonian gravity derived from rotation curves (Borgani et al. 1991), which was interpreted as a two-point correlation function for the galaxy background in the range of 3–350 kpc, highlighting the coincidence of the exponent (1.76) with that of the two-point correlation function among galaxies (1.77) (Borgani et al. 1991).

The total equivalent mass associated with this density ρ* within a sphere of radius r is

$$\begin{eqnarray*}\begin{array}{rcl}{M}^{* }(r) & = & 4\pi {\displaystyle \int }_{0}^{r}{dr}^{\prime} \,r{{\prime} }^{2}\rho * (r^{\prime} )\\ & = & M+\displaystyle \frac{M}{{a}_{0}}{\displaystyle \int }_{\infty }^{g(r)}{dg}^{\prime} \displaystyle \frac{1}{\mu {\left(x^{\prime} \right)}^{2}}\displaystyle \frac{d\mu }{{dx}}(x^{\prime} )\end{array}\end{eqnarray*}$$

$$\begin{eqnarray}&&=\,M\left(1+{\int }_{\infty }^{\mu (r)}d\mu ^{\prime} \,\displaystyle \frac{1}{\mu {{\prime} }^{2}}\right)=M\mu {\left(r\right)}^{-1}.\end{eqnarray} \tag{ 22 }$$

The total mass diverges as r → ∞ and μ → 0,

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{r\to \infty }{M}^{* }(r)=\sqrt{\displaystyle \frac{M\,{a}_{0}}{G}}r,\end{eqnarray} \tag{ 23 }$$

but only M is a real mass, while M^∗(r) − M stems from a fictitious density field that gravitationally behaves as Newtonian mass but does not correspond to any real mass. It is a phantom mass. One needs to imagine that each body generates a potential around it, which defines its gravitating mass (the phantom mass), and this potential is changed depending on which other gravitating bodies are around it and where they are, so the equivalence between gravitating mass and inertial mass is broken (Wu & Kroupa 2015). This concept of phantom mass as a tool to compute the MOND potential has indeed been raised to an exact concept in QUMOND, precisely motivated by earlier considerations similar to those given here but in the solar system context (Milgrom 2009). The computing of the phantom mass is the exact way to solve the QUMOND Poisson equation.

In practice, this phantom mass does not reach infinity values, because the galaxies (to be considered point-like objects at large distances) cancel their gravitational fields in regions where other galaxies have a predominant effect. Equation (19) is applicable in the volume where the acceleration of the galaxy is predominant, roughly on average up to a distance of half of the average separation among galaxies. This phantom mass would be the substitute of the nonbaryonic dark mass in the standard cosmological model, although it is not clear whether a MOND cosmology can be built, for which there are arguments in favor or against (Felten 1984; Sanders 1998).

We could say that the idea of dependence on relative accelerations would be impossible to sustain based on theoretical grounds. There is a historical discussion already from the times of Helmholtz: a field cannot have a huge N degrees of freedom as would be required if the force were to depend on all of the relative accelerations with each of the N particles of the gravitational interaction. A dependence on the absolute acceleration makes more sense. Nonetheless, the concept of absolute acceleration is also problematic, since this requires an absolute frame of reference in the universe. Which is this absolute reference system? The cosmic microwave background radiation? But then the center of each galaxy has some acceleration with respect to that system. On the other hand, if we put the absolute reference frame in the center of a given galaxy, how can we understand the motions in other galaxies, which would have non-MONDian accelerations with respect to the first one? It is not clear either whether the acceleration is in the comoving or physical cosmological frame, and as mentioned, it is not even clear whether we may have a MOND cosmology (Felten 1984; Sanders 1998). For considering the comoving frame, we would need a well-understood metric, equivalent to the one derived from general relativity. If we considered a Newtonian-like approach and paid attention only to the physical accelerations with an expansion of the universe equivalent to the standard model, we would have a relative redshift drift among galaxies with $\dot{z}\approx {H}_{0}\,z$ at low z (Bolejko et al. 2019), so the relative acceleration would be $g=c\dot{z}\approx 5.7{a}_{0}z$, which is not negligible and can be considered in the Newtonian regime when the distance between galaxies is high (z ≳ 0.17). Something remains unclear with respect to the concept of absolute acceleration. Certainly, the concepts of MOND are slippery, but we may forget about the conceptual theoretical problems and see whether the phenomenological rules can be applied at least within one galaxy.

There are some attempts to clarify the question of the superposition of fields, distinguishing the external field effect (Milgrom 1983a; Famaey & McGaugh 2012, Section 6.3) and the internal field effect and claiming that the internal accelerations of the subsystem are irrelevant to how that subsystem responds to an external field. Only the field of the parent system at the position of the center of mass of the subsystem is relevant to that. Only the center of mass of such systems matters to determining their orbits in MOND, not their internal structure, nor the magnitude of their internal accelerations. MOND has to be described by a nonlinear theory. This basically means that the acceleration endowed by two bodies to a third is not the (vectorial) sum of the individual accelerations produced by each separately.

Bekenstein & Milgrom (1984) derived mathematically two consequences of the superposition of fields: (1) that the acceleration of the center of mass in a system much smaller than the source of the external field (e.g., a star subject to the external field produced by the rest of the galaxy) follows the acceleration imposed by the external object when the radius of the system trends to infinity (Bekenstein & Milgrom 1984, Section IV); and (2) if the external field is Newtonian, assuming also a large enough radius, the internal field of an arbitrary mass in the system is Newtonian even though the internal accelerations are much smaller than a₀ (Bekenstein & Milgrom 1984, Section V). Nonetheless, it is not clear what can be considered external or internal. Given an atom that feels two gravitational accelerations, how can it distinguish whether the force comes from a nearby source or a distant source? Or is MOND also dependent on distance apart from the dependence on the acceleration? In principle, it is not. Certainly, a particle only feels a total acceleration without distinguishing where it comes from. The distinction between internal and external field is not something a particle is aware of. This was indeed the apparent paradox that one of us introduced (López-Corredoira 2018). However, these suspicions of contradiction were not correct.

From our analysis, the explanations of these cases are as follows: On Earth, for instance, the phantom mass (the second term in Equation (12)) of the fields created by the atom or the small particle is canceled, due to the action of the strong field of Earth in which it is embedded. Also, for a multiple-star system in the inner disk of a galaxy within the Newtonian regime, similar argumentations can be given. Therefore, we cannot explore MOND effects in binary stars or globular clusters in the inner Galaxy, as apparently found in some observations (e.g., Scarpa et al. 2017), unless this MOND logic does not apply for some other reason.

However, on the interior of a star in the outer disk of a galaxy, the self-gravity cannot cancel the MONDian phantom mass created by the center of the galaxy because almost all of the space filled by this mass is far from the volume of the Newtonian regime owing to self-gravity of the star, which is negligible in comparison with the total volume. Therefore, MONDian acceleration is applied over all of the atoms of the star, and consequently the center of mass of the star follows a MONDian dynamics. On a binary or multiple-star system in the outer disk of a galaxy, similar argumentation can be given.