1 Introduction

The evolution of the Universe, the kinematics, and the stability of galaxies is usually explained by the presence of Dark Matter (DM) [1]. The standard model used is the \(\Lambda \) Cold Dark Matter (\(\Lambda \)CDM) model, in which the galaxies’ rotation velocities as well as their stability is ensured by DM [2, 3]. The Newtonian inverse square gravitational attraction from baryonic matter is, however, not included in the standard model, and this has given rise to a series of proposed modifications of Newton’s inverse square law (ISL) for gravitational attraction with the intention of unifying gravity with particle physics [4,5,6,7].

The rotation velocities of galaxies deviate from the Classical Mechanics rotation velocities [8,9,10,11], this behavior is taken into account in the theories that have led to the assumption of the existence of DM in the Universe [12]. Another attempt to explain the stability and the dynamics of galaxies is the MOND theory, where the stability of the galaxies and their rotation velocities is caused by modified Newtonian dynamics at small accelerations [13]. The Newtonian dynamics can either be changed by modifying the accelerations of baryonic objects at small accelerations (MOND) or by modifying the gravitational attractions (MOGA) from objects located far away [14]. The MOND theory is reviewed in [15].

Modifications of the gravitational attraction (MOGA) have been proposed for a long time, and modified gravity are reviewed in [16,17,18]. The modification could be caused by the general relativistic gravitational lensing (GL) which deals with the effect of wave deflection caused by baryonic objects. GL is reviewed in [19, 20]. Gravitational waves were first detected in 2016 [21]. Simultaneously detection of gravitational and electromagnetic waves from the coalescence of binary neutron stars rules, however, out a class of MOND theories [22].

The kinematics of the galaxies is determined by radiation from the galaxies or by simulations. Simulations of galaxies have been performed for many decades. The dynamics of baryonic objects can be solved by the Particle-Particle/Particle-Mess (PPPM) method [23,24,25], which is a mean-field approximation where each mass unit is moving in the collective field of all the others, and the Poison equation for the PPPM grid is solved numerically. The dynamics are determined for Zeldovich’s adiabatic expansion of the Einstein-de Sitter universe [26]. Later, simulations with large-scaled computer packages with PPPM (GADGET, PHANTOM, RAMSES, AREPO) [27,28,29] are with many billions of mass units. The evolution of galaxies has been obtained from hydrodynamical large-scaled cosmological simulations [30,31,32], and with co-evolving dark matter gas and stellar objects [33, 34]. Simulations of galaxies with MOND dynamics have also been performed by (approximated) Poisson solvers and [35,36,37]. Cosmological simulations of galaxy formation are reviewed in [38].

Here, we simulate models of galaxies with the modified Newtonian dynamics, MOND, and with modified gravitational attraction, MOGA. The MOND acceleration is obtained from the interpolation function(s) [10, 13], and the modified gravitational attraction, MOGA, is obtained by replacing Newton’s gravitational inverse square attraction (ISL) with an inverse attraction (IA) for large distances. The simulations of galaxies in an expanding Universe are performed by use of a recent extension of Newton’s discrete algorithm [39,40,41,42]. The algorithm is absolutely stable and without any approximations, and the dynamics with the discrete algorithm have the same invariances as Newton’s analytic dynamics, and thus the discrete dynamics is exact in the same sense as the exact solution of Newton’s analytic second-order differential equations for classical celestial mechanics. The dynamics with Newton’s discrete algorithm is reviewed in [43], and the algorithm and the proof that Newton’s discrete dynamics has the same invariances as his analytic dynamics is given in the “Appendix.” The advantage of using Newton’s discrete algorithm is that the classical dynamics are without any approximations, and the disadvantage is that it is only possible to perform the exact long-time simulations for small ensembles of objects in bound rotations around their common mass center.

The Newtonian dynamics with MOND and with MOGA are formulated in the next section. The acceleration in MOND is modified when the acceleration of an object is below a certain threshold, \(a_0\), and independent of the individual attractions from the objects which together cause the small acceleration. This modification breaks Newton’s third law and unlike MOGA, MOND no longer conserves the momentum and angular momentum for a conservative system(see the “Appendix”, Eq. (A7), Fig. 6, and [44]), and the present simulations reveal that MOND is unstable. Section 3 presents the results of the MOND and MOGA simulations of galaxies for times corresponding to the age of the Universe. The simulations show that the MOND dynamics is unstable and releases the bound objects in the galaxies with time, whereas MOGA not only stabilizes the bound rotations of the objects in the galaxies but it also changes the rotation velocities and is qualitatively in agreement with experimentally determined rotation curves for galaxies in the Universe.

2 Modified acceleration, MOND, and modified gravitational attraction, MOGA

Newtonian dynamics for celestial objects is given by Newton’s classical dynamics

$$\begin{aligned} {\textbf {F}}_{i}(t)=m_i \text {a}_i \hat{{\textbf {a}}}_i(t) \end{aligned}$$
(1)

and his ISL law of gravitation

$$\begin{aligned} {\textbf {F}}_{i}(t)=\text {F}_i \hat{{\textbf {a}}}_i= -\sum _{j \ne i}^N {\frac{m_i m_j G}{r_{ij}^2(t)} \hat{{\textbf {r}}}_{ij}(t)} \end{aligned}$$
(2)

for the force \({\textbf {F}}_{i}(t)\) and the acceleration \( {\textbf {a}}_i(t)= \text {a}_i \hat{{\textbf {a}}}_i(t)\) for the object i in the ensemble of N objects, caused by baryonic objects j with mass \(m_j\) at distances \(r_{ij}(t)\) and at time t.

The first equation is Newton’s second law for the relation between a force acting on an object, and Newton postulated Eq. (1) in the very first part of Principia and continued by derived Eq. (1) at page 37 from a discrete analog formulation of the change in the position by a force impulse \({\textbf {F}}_{i}(t)\) at time t (see Apendix A). Much later Newton starts on page 401 in Principia by formulating four principles for philosofical rules in natural science. Rule I in an English translation reads:

We are to admit no more causes of natural things than such as are both true and sufficient to explain their appearances.

After the formulation of the philosophical principles, he succeeds in formulating his law of universal gravitation, given by Eq. (2). The law was obtained from experimental data for the positions of the planets in the Solar system and the positions of the Moon. MOND and MOGA acknowledge the philosophical principles in an attempt to explain the stability and behavior of the dynamics in the Universe. MOND by modifying the acceleration and MOGA by modifying the ISL.

The MOND theory was published in 1983 by M. Milgrom [13]. The acceleration \({\textbf {a}}_i(t)\) in classical Newtonian dynamics is modified for a small acceleration caused by the sum of interactions with the \(N-1\) baryonic objects located at large distances \(r_{ij}\) from i. The transition from the Newtonian acceleration to MOND occurs for a small acceleration \(a_i \approx a_0 =\mid {\textbf {a}}_0\mid \). This can only happen if all the attractions with object i from the other \(N-1\) baryonic objects are for large distances so that the sum of attractions results in a small acceleration \( a_i< a_0\)

$$\begin{aligned} a_i \rightarrow \mu (a_i/a_0)a_i. \end{aligned}$$
(3)

According to Bekenstein and Milgrom [14] the modification can either be performed by modifying the law of inertia (MOGA) or by modifying the acceleration (MOND). The modification, \(\mu \), of the acceleration \(\mid {\textbf {a}}_i\mid \) is obtained with the ”standard interpolation function”

$$\begin{aligned} \mu (\text {a}/a_0)=\sqrt{ \frac{1}{1+\left( \frac{a_0}{\text {a}}\right) ^2}}, \end{aligned}$$
(4)

The MOND acceleration is given by

$$\begin{aligned} \text {F}_i/m_i= \text {a}_i \sqrt{ \frac{\text {a}_i^2}{\text {a}_i^2+a_0^2}}, \end{aligned}$$
(5)

or the interpolation function proposed by [10]

$$\begin{aligned} \mu (\text {a}/a_0)=\frac{\mid \text {a}\mid }{ \mid \text {a}\mid +a_0}, \end{aligned}$$
(6)

with the modification given by

$$\begin{aligned} \text {F}_i/m_i = \text {a}_i \frac{\mid \text {a}_i\mid }{ \mid \text {a}_i \mid +a_0}, \end{aligned}$$
(7)

and acceleration

$$\begin{aligned} \text {a}_i(\text {MOND}) =\frac{ \text {F}_i}{2m_i}\left( 1+\sqrt{1+4 m_i a_0/\mid \text {F}_i\mid }\right) . \end{aligned}$$
(8)

The MOND modifications, Eqs. (5) or (7), change the acceleration from the classical Newtonian acceleration \(\mid {\textbf {a}}\mid a_0 \) at short distances to a modified acceleration

$$\begin{aligned} \mid {\textbf {a}}_i(\text {MOND})\mid = \sqrt{ \mid \text {F}_i\mid a_0/m_i} \end{aligned}$$
(9)

for \( \mid {\textbf {a}} \mid<< a_0\).

The asymptotic modified acceleration for an isolated object i with only one gravitational interaction, \(\text {F}_i(r_{ij})=-m_im_jG/r_{ij}^{2}\), with another object, No. j, is obtained from \(\hbox {F}_i(r_{ij})\) and Eq. (9) as

$$\begin{aligned} \text {a}_i(\text {MOND})= - \frac{ \sqrt{m_jG a_0}}{r_{ij}}. \end{aligned}$$
(10)

MOND is a modification of Newtonian acceleration. But in this case, the modification might as well be formulated as a modification of Newton’s ISL law of universal gravitational attraction, where the inverse square attraction asymptotically is replaced with an inverse attraction (IA). If this modified gravitational attraction is a universal law, MOGA, the gravitational force is modified to

$$\begin{aligned} {\textbf {F}}_i(\text {MOGA})=-\sum ^N_{j \ne i}\frac{m_i m_j G}{r_{ij}^2}\left( 1+ \frac{r_{ij}}{r_0}\right) \hat{{\textbf {r}}}_{ij}(t) \end{aligned}$$
(11)

with \(r_0=\sqrt{m_jG/a_0}\).

The modified gravitational attraction, Eq. (11), is a specific example of a general modification of the force field [14, 16]. Newtonian dynamics and Newton’s discrete algorithm used in the next section are time reversible, symplectic and with the dynamical invariances: momentum, angular momentum, and energy for a conservative system [43]. MOGA with Newton’s discrete algorithm, maintains these qualities, whereas MOND does not conserve momentum and angular momentum (see “Appendix A”, Fig. 6 and [44]). The Hubble expansion of the space destroys the invariances [42], but MOND and MOGA dynamics are, however, still time reversible.

Modifications of the Newtonian ISL attraction have been proposed for a long time in an attempt to obtain the stability of galaxies by the standard model [4, 5, 7, 45, 46], and the dynamics of local group of galaxies (Milky Way and Andromeda) have been described by modified Newtonian attraction [47, 48]. For a review of extended theories of gravity see [16, 49].

The modified Newtonian gravitational potential u(r) is often modified by a Yukawa potential [48, 50, 51, 53]

$$\begin{aligned} v_{ij}(r)=u(r_{ij})[1+\alpha exp(-r_{ij}/\lambda )], \end{aligned}$$
(12)

and the corresponding modified gravitational forces are

$$\begin{aligned} \text {F}_i = -\sum ^N_{j \ne i}\frac{m_i m_j G}{r_{ij}^2}[1+\alpha exp(-r_{ij}/\lambda )(1+r_{ij}/\lambda )]. \end{aligned}$$
(13)

The parameter \(\lambda \) corresponds to the parameter \(r_0\) in MOGA, and a possible deviation from the ISL gravity was investigated for \(\lambda \) in the range \(\lambda \in [30,8000]\) nm [51,52,53], and the rotation velocities of stars in the Milky Way were used to determine a possible Yukawa correction, Eq. (12) to the gravitational attraction [7]. So far, however, there is no direct experimental evidence for deviation from the ISL gravitational attraction.

Fig. 1
figure 1

An illustrative figure of the modified accelerations a = F/m of an object as a function of the distance r to a heavy object with a mass \(M=1000m\). Length unit is in parsec and force (F = \(mMG/r^2)\) is in units of m and G. The MOND acceleration (Eq. 8) is in blue, the MOGA acceleration (Eq. 11) is in magenta, and the acceleration with the modified Yukawa attraction (Eq. 13) is in green. The accelerations are for \(r_0=\lambda \) = 1000 parsec. The Newtonian acceleration with the ISL (Eq. 2) is shown in red. The inset enlarges the differences at \(r \approx r_0\)

Full size image
Fig. 2
figure 2

The orbits of an object near a mass center in black which is thousand times heavier than the object. The dynamics was started at Perihelion at a distance to the center r(Perihelion)=5000 pc and with an excentricity \(\epsilon =0.80\). The green ellipse is the classical orbit. In blue are the first 38 corresponding (revolving) orbits with MOGA and with \(r_0=25000 \) pc. The corresponding 26 revolving orbits in red are for the MOGA with the Hubble expansion included in the dynamics [42]

Full size image

An illustrative example of the different accelerations is shown in Fig. 1. (Units are with \(m_i=G=1\) and length unit 1 corresponds to 1 pc for a galaxy. For units see [42]). The accelerations are for one object attracted by a heavy mass center with a mass that is a thousand times heavier than the object and for a modification distance \(r_0=\lambda =1000\), which corresponds to 1 kpc for the model for the Milky Way in the next section. The Yukawa modification is shown in green, the MOND acceleration is in blue, and the MOGA acceleration in magenta. The MOND and MOGA modifications are rather similar in this case, but MOGA and MOND for galaxies with many gravitational objects are, however, very different. The difference between MOGA and MOND for the dynamics of a galaxy originates from the summation of interactions between the objects. The acceleration in MOND is modified if the interactions of an object with all the other baryonic objects result in an acceleration below a certain threshold, whereas the contribution to the acceleration from the attractions of all faraway objects is modified in MOGA. The green curve is the force from the modification with the Yukawa potential, Eq. (13), and with \(\alpha =1\) and \(\lambda =\)1000. If there exist dark matter this modified force will be effective on a much shorter length scale with baryonic attractions from dark matter in the galaxies and their halos.

The mean rotation velocities of stars in galaxies are rather constant and independent of the distance to the centers of rotation as opposed to a system of rotating baryonic objects with classical Newtonian dynamics [10, 13]. This led Milgrom to propose the modification, Eq. (4) (or Eq. 6) of the Newtonian classical acceleration, and with the asymptotic modification Eq. (9). By adjusting the modified acceleration to one (isolated) baryonic object in rotation at a gravitational center he determined a value for the constant \(a_0\). Milgrom found \(a_0 \approx 1.2 \times 10^{-10}\) \(\hbox {ms}^{-2}\) to be optimal, and later investigation of the rotation curves for stars in 12 galaxies confirmed this value [10]. For a MOND modification caused by only one baryonic object with mass \(m_j\) equal to the mass of our Sun: \(m_j=m_{\text {Sun}}\)=1.989\(\times 10^{30}\) kg and with the gravitational constant \(G=6.647 \times 10^{-11}\) N \(\hbox {m}^{2}\) \(\hbox {kg}^{-2}\) the constant \(a_0\) corresponds to a distance

$$\begin{aligned} r_0 &= \sqrt{\frac{m_jG}{a_0}}= \sqrt{\frac{1.99 \times 10^{30} \times 6.65 \times 10^{-11}}{ 1.2 \times 10^{-10}}} \nonumber \\ & = 1.05 \times 10^{15} \text {m}=\frac{1.05 \times 10^{15} \text {m}}{3.09 \times 10^{16} \text {m} \text {parsec}^{-1}} \nonumber \\ & = 0.033 \text {parsec}, \end{aligned}$$
(14)

The example is for one interaction from an object with mass equal to our Sun. The modification distance is proportional with the square root of the masses, \(m_j\), and the distance is much larger for a modification caused by many heavy objects far away. The Milky Way is a barred spiral galaxy. The extension of the barred disk is 2.5–3 kpc [54], and the extension of the halos is \(\approx \) 100–300 kpc [55, 56], so the galaxy will be affected by a modification of the gravitational attractions or the accelerations.

MOND and MOGA dynamics are equal in the case of only one object attracted by a heavy mass center according to Eq. (11). Figure 2 shows the dynamics for one object in rotation around a heavy mass center with a mass one thousand times bigger than the object’s mass. The green curve is the elliptical orbit of the object with pure classical Newtonian dynamics. The dynamics was started in Perihelion with a distance to the mass center r(Perihelion)=5000, which corresponds to 5000 Parsec for the Milky Way, and with an excentricity 0.8 by which the object’s distance to the center is 45,000 in Aphelion. (The relations between the units for length l: 1 pc, time t: 1 Gyr and Hubble expansion coefficient \(H=72.1 \pm 2.0 \ \text {km} \text {s}^{-1} \text {Mpc}^{-1}\) [57] in the Universe and the corresponding units for the MD models of galaxies are: 1 pc \(\hat = \) 1 MD length unit, 1 Gyr \(\hat = 6\times 10^5\) MD time units and \( 72.1 \text {km} \text {s}^{-1} \text {Mpc}^{-1} \hat = 5.\times 10^{-8}\) in MD units. For determination of these relations see [42].) The blue curves are the first 38 revolving orbits with MOND = MOGA and with the modification distance \(r_{0}=25000\) in Eq. (11). The modified dynamics results in stable revolving orbits, but with a smaller mean distance to the mass center [39, 41]. The red curve is the first 26 corresponding orbits for MOGA dynamics with a Hubble expansion equal to the Hubble expansion of the Universe [42]. The dynamics still exhibit revolving orbits, but with an increasing mean distance. MOGA with Hubble expansion performed 64 regular orbits before the expansion released the object from the mass center, and the object performs many more orbits for stronger mass centers before its release. The Milky Way has performed \(\approx \) 60 rotations after its creation more than thirteen billion years ago and the Hubble expansion does not at all affect the stability of the revolving orbits for MOGA with heavier mass centers in times corresponding to the age of the Universe. The total mass of the Milky Way is extimated to \(2.06 \times 10^{11}-5.4 \times 10^{11}\) in unit of the mass of the Sun [58].

3 Simulations of galaxies with MOGA and with MOND

Simulations of galaxies have been performed for many decades [59], but the simulations here of models of galaxies deviate from the main part of the simulations in that they are pure Newtonian N-body simulations. Such simulations have, however, also been simulated for a long time. But one has performed a series of approximations, such as variable time steps and mean field approximations, because these simulations are very time-consuming. But the approximations ruin the exactness of the simulations. The simulations below are exact N-body Newtonian simulations without any approximations. The algorithm with the conserved dynamic invariances is given in the “Appendix” and in a recent review article about Newton’s exact discrete dynamics [43].

A galaxy and the Milky Way contain hundreds of billions of stars, and a substantial amount of baryonic gas [58, 60, 61], and it is not possible to obtain the exact dynamics with MOGA or with MOND of a galaxy with this number of objects. We have instead of simulated models of small ”galaxies” of hundred of objects in orbits around their center of gravity, and in an expanding space with various values of \(r_0\) or \(a_0\). A recent article describes how a system of baryonic objects with Newtonian discrete dynamics spontaneously creates a system with the objects in rotation about their center of gravity [41]. The algorithm is extended to also include the dynamics with the Hubble expansion of the space [42].

The algorithm is used to simulate models of a galaxy with the Hubble expansion and with MOGA and MOND, respectively. An ensemble of gravitational objects with Newtonian dynamics at time t=0 might spontaneously create a ”galaxy” system with many of the objects in a bound rotation around the center of gravity, and one can either start the simulations with MOGA or MOND at \(t=0\), or alternatively at a later time t where a Newtonian galaxy is created and it is in a rather stable state. The data reported below for twelve galaxies are started from a stable Newtonian galaxy. The twelve galaxies with the data reported below are for \(r_0= 1, 10, 100, 1000, 10,000, 100,000\) and \(a_0 = 1, 10^{-2}, 10^{-4}, 10^{-6}, 10^{-8}, 10^{-10},\), respectively, and each (stable) galaxy are simulated \(3.2 \times 10^9\) time steps corresponding to \(t=8 \times 10^6 \approx \) 13.4 Gyr or the age of the Universe. Each simulation with the \(3.2 \times 10^9\) time steps took \(\approx \) 1000 h on a fast CPU in the CPU-cluster, and the total amount of simulations for different values of \(r_0, a_0, H\), and start configurations is fifty-three.

Fig. 3
figure 3

The number of objects with mean distances \(r(\text {mean}) < 100000\) parsec to the center of mass as a function of time \(t \in [0,8 \times 10^6] \approx \) 13.4 Gyr where the galaxies are exposed to either MOND or MOGA. The red curve is for Newtonian dynamics without any modification of accelerations or attractions. The green curves are for a very weak/long range modification \( a_0= 10^{-10}\) (MOND) corresponding to \(r_0=100000\) (MOGA), and the blue curves are for a strong/short range modification \( a_0=r_0=1\). The magenta curve is MOND dynamics with \( a_0= 10^{-8}\)

Full size image

3.1 Stability of galaxies with MOGA or with MOND

The dynamical effect of the modifications of the accelerations or gravitational forces is obtained by exposing the objects in a galaxy to MOGA or to MOND. The results reported below are obtained by exposing the objects in a Newtonian galaxy which is in a rather stable state and with an occasional release of an object before the modifications (Curve in red in Fig. 2). At the start of the modifications the system contains 460 objects with 165 with a mean distances \(r(\text {mean}) < 100000)\) parsec to the center of mass, and 295 objects with mean distances \(r\text {(mean}) > 100 000 \approx 100\) kpc. The number of objects with mean distances \(r(\text {mean}) < 100000)\) in the galaxy after the change in dynamics is shown in Fig. 3.

The stars in the halos of a galaxy are not in a stable and bound rotation, and this fact has led to the hypothesis about the existence of dark matter in the Universe. The simulations of Newtonian galaxies show, however, that the release of stars in the outer edge of a Newtonian galaxy is rare, and that the Newtonian galaxies are rather stable over time periods of many Gyr (red curve in Fig. 3) [42]. By changing the Newtonian dynamics to MOND the enhanced MOND acceleration destabilizes the galaxy and results in a release of the bound objects, whereas the MOGA dynamics has the opposite effect. The number of bound objects with \(r_0=a_0\)=1 is shown in blue, and the number with \(a_0=10^{-10}\) (MOND) and correspondingly \(r_0\)=100,000 (MOGA) is shown in green in Fig. 3. MOGA stabilizes the galaxy even for a very weak modification with \(r_0 \approx 10000- 100 000\) corresponding to a modification of the attractions at distances \(\approx \) 10–100 kpc for the Milky Way. The galaxies with MOGA still release objects, but this is very rare and the galaxies contain still more than 135 bound objects after a time \(\Delta t= 8 \times 10^6 \) corresponding to 13.4 Gyr, or the age of the Universe. The simulations are performed for \(r_0= 1 (136), 10 (143), 100 (142), 1000 (156), 10,000 (156)\) and 100, 000(142), respectively, and with the number of objects with mean distances \(r(\text {mean}) < 100 000\) at the end of the simulations indicated in the parentheses. The number of objects in the Newtonian galaxy with \(r_0 =\infty \) at the end of the simulation (red curve in Fig. 3) is \(N=\)107, so the MOGA galaxies contain significantly more bound objects than the Newtonian galaxy. The stability of the galaxy given by the number of bound objects with MOGA is not sensitive to the range \(r_0\) of the modification of the gravitational attraction.

The number of bound objects (\(r(\text {mean}) < 100000\)) with MOND dynamics is also shown in Fig. 3. The MOND dynamics releases the objects even for a very low threshold for the modified acceleration, given by \(a_0=1/r_0^2\). It was only possible to maintain some bound objects for a very weak modification of the Newtonian accelerations with \(a_0=10^{-10} \rightarrow r_0=100000\) (green MOND curve in Fig. 3). All the bound objects were released for a stronger MOND modification of the acceleration.

Galaxies with MOGA or with MOND dynamics were simulated with other start distributions of the baryonic objects and for \(H=0\) (i.e., without a Hubble expansion of the space). The simulations showed unanimously, that MOGA dynamics has a stabilizing effect on the objects in the galaxies even for a big value of \(r_0\) corresponding to that only the stars in the halos of a galaxy are affected by the modified gravitational attraction. All the MOND simulations were unstable and released the bound objects sooner or later. MOND does not conserve the angular momentum of the ensemble (see “Appendix A”, Eqs. (A7) and (A9) and Fig. 6), and the increased MOND acceleration has the opposite effect, it destabilizes the bound objects in the galaxies even for the very small value of \(a_0=10^{-10}\) and the MOND galaxies release the objects from their regular orbits in the galaxies. The same was thru for dynamics without a Hubble expansion (\(H=0\)), where the galaxies were stable for MOGA, but unstable for MOND.

3.2 Rotation velocity of the stars in a galaxy

There is a discrepancy between the mean rotation velocities of the stars in the galaxies and the corresponding mean velocities of the objects in the simulated galaxies with pure Newtonian dynamics. The mean rotation velocities of the stars in a galaxy as a function of the distance to the center of rotation increase for short distances, but are rather constant at larger distances to the center of the galaxies. [10, 11, 62]. But the velocity v(r) of objects with pure Newtonian dynamics in mean declines as \(v(r) \approx (M/r)^{0.5}\) with the distance r to the mass center of the galaxy with mass M. The rotation velocity of stars in the Milky Way at the distances \(r \in [6.72, 8.40]\) kpc is \(v(r) \in [203, 240]\) km \(\hbox {s}^{-1}\) [63]. One of the results in [42] from the simulation of galaxies with pure Newtonian dynamics was, that the velocities of objects for distances \(r \in [15000,100 000]\) to the center were not located near a line with \(v(r) = (M/r)^{0.5}\), but were diffusely distributed. But a determination of the rotation velocities over a longer time interval with the root mean square (rms), obtained from 22 consecutive time intervals of the mean velocities, reveals that the distribution of the mean rotation velocities is rather Newtonian, and disagrees with the \(\approx \) constant rotation velocity in galaxies (Fig. 4).

Fig. 4
figure 4

Mean rotation velocities of the objects in the galaxy with Newtonian dynamics in red (i.e., without a modification of accelerations or attractions). Distances are in parsec and velocities are in MD units where \(v(r)=0.15\) corresponds to \(v=220\) km/s in the Milky Way. In green is the Newtonian rotation velocity \(v(r)=(M/r)^{0.5}\) for one bound object attracted by a heavy mass center. The blue line is \(v(r)=0.15\). The standard deviations for (representative) selected distances are obtained from 22 consecutive time subsets of the mean rotation velocities in the subsets. The inset shows the mean velocities at short mean distances

Full size image

The rotation velocities of the objects in a galaxy with pure Newtonian dynamics are shown in Fig. 4. The rotation velocities in red and standard deviation with magenta are obtained in the time interval \(t \in [0, 8\times 10^6]\). The green curve is the function \(v(r) = (680/r)^{0.5}\) for one object in circulation around a mass center with a mass M=680 times the mass of the object. The blue straight line is \(v=0.15\), and it corresponds to \(v=220\) km/s for stars in the Milky Way [63]. The rotation velocities decline even more rapidly with the distance to the center of rotation than given by \(v(r)=(680/r)^{0.5}\), and in disagreement with the observed mean rotation velocities of the stars in galaxies.

The dynamics with MOGA increase the rotation velocities at distances \(r> r_0\) to the center of the galaxy. The rotation velocities for MOGA dynamics and for different values of \(r_0\) are shown in the next figure. The rotation velocities are for the bound objects with the distances \(r < 100 000\). The mean velocities of the objects in Fig. 5 are determined in the time interval \(t \in [0,8\times 10^6]\). They are constant within the accuracy of the simulations already for a modified attraction with \(r_0=10000\) corresponding to 10 kpc for the Milky Way and a smaller value of the distance \(r_0\) for the onset of the modification only increases the constant mean velocity. So a modification of the ISL attraction to an inverse attraction not only stabilizes the galaxy, (Fig. 3) but also increases the rotation velocities at large distances and make them rather constant with respect to the distance to the center of rotation.

Fig. 5
figure 5

The MOGA (\(log-\)) rotation velocities of the objects in the galaxy models of the Milky Way as a function of the distances to their center of mass in the galaxies and for various values of the distance \(r_0\) for the onset of the inverse attraction. Distance is in parsec and the velocities are in MD units where a mean velocity \(\bar{v} \approx 0.15-0.25\) corresponds \(\bar{v}=220\) km \(\hbox {s}^{-1}\) in the Milky Way. The velocity interval in magenta is for velocities that correspond to rotation velocities \(v \in [0.15, 0.25]\)

Full size image

The rotation velocity of stars in the Milky Way can be related to the rotation velocities in the MOGA models of a galaxy. The relation between lengths is obtained from the sizes of the galaxies, and the relation between velocities is obtained from the rotation time and velocity in the Milky Way and the corresponding rotation time and velocities in the models of galaxies. The rotation velocity of stars in the Milky Way at the distances \(r \in [6.72, 8.40]\) kpc is \(v(r)=220 \pm 10\) km \(\hbox {s}^{-1}\) [63]. An estimate of the mean velocity of the objects in the models for a galaxy is associated with uncertainty. The mean rotation velocity \(\bar{v} \approx 0.15-0.25\) in the galaxy models is determined to correspond to the mean rotation velocity \(v=220\) km \(\hbox {s}^{-1}\) in the Milky Way [42], and the colored area in Fig. 5 is for mean rotation velocities in the interval \(v \in [0.15,0.25]\). The IA attractions increase systematically the mean velocities of the baryonic objects, and the rotation velocities in the MOGA galaxies only agree with the observed velocities in the Milky Way for a large value of the modification distance \(r_0\). The simulations indicate that the stability of a galaxy with a rather constant rotation velocity which corresponds to \(\approx \) the rotation velocity in the Milky Way can be obtained by a modification of the attractions for far away objects with distances corresponding to \(r \ge 10\) kpc (blue curve in Fig. 5). The estimation is with a big uncertainty, but one can, however, exclude a significantly shorter value of the modification range \(r_0\). A smaller value than \(r_0\)= 10 kpc will according to the results in Fig. 5 lead to an unrealistic high rotation velocity. The enhanced accelerations (MOND) or attractions (MOGA) for \(r_0 \le 1000 \rightarrow a_0 \ge 1.0\times 10^{-6}\) result in a too high mean rotation velocity for MOGA, and in the case of MOND, the galaxy is spontaneously destabilized.

4 Conclusion

4.1 Summary of the simulations

The Universe is more than thirteen billion years old and the Milky Way was created relatively shortly after the birth of the Universe. An appropriate time unit for the evolution of the Universe is, however, not a year but the time of one rotation of a galaxy. The Milky Way has rotated about sixty times with \(\approx \) 240 million years per rotation [63], and a Newtonian system of stars is rather stable even after the double number of rotations [42]. But the rotation velocity of the objects in the Newtonian system with ISL attraction is not constant but declines proportional to the distance r from the center of mass as \(r^{-0.5}\) (Fig. 4) and inconsistent with the rate of rotation of galaxies [10, 11, 63].

MOGA with the modification of the gravitational attraction from a Newtonian ISL attraction to an IA attraction for pairs of baryonic objects at distances \(r \ge r_0\) stabilizes the objects in the galaxy further (Figs. 2 and  3), and the modification changes the mean rotation velocity of objects in the galaxy from a classical Newtonian/Kepler behavior to a rather constant velocity to the distance to the center of rotation in the galaxy (Fig. 5). The rotation velocity in the galaxies with MOGA is in qualitative agreement with the observed rotation velocity in the Milky Way and galaxies in the Universe.

4.2 Discussion

Newton began Principia by postulating his three laws, and then he derived the second law, Eq. (1), from a discrete algorithm which today is the most used algorithm in computer simulations. The algorithm is derived from the action of a force quant [41], more than three hundred years before the formulation of Quantum Mechanics, Feynman’s path integral formulation, the General Relativity theory (GR), and quantum electrodynamics (QED). Newton’s formulation of the ISL of universal gravitation, Eq (2), was derived much later in Principia , on page 401 in the Third Book Mundi Systemate, and where Newton first formulated four philosophic principles for Natural Science. The first and in an English translation reads:

We are to admit no more causes of natural things than such as are both true and sufficient to explain their appearances.

MOGA acknowledges this principle. But MOND does not obey Newton’s third law, it does not conserve momentum and angular momentum (Fig. 6) and the MOND galaxies are unstable (Fig. 3).

The present results raise some questions: Can the dynamics of the Universe be explained alone by a modification of Newton’s laws, his universal gravitation, and QED, or do one also need explicitly to include an effect of dark matter in the dynamics?

And if a modification of the ISL law is sufficient, what causes this modification and how should this modification be formulated?

None of the two questions can of cause be answered definitively from the present investigation.

Modifications of the Newtonian ISL to an IA-like attraction have been proposed for a long time in an attempt to obtain the stability of galaxies by mean of the standard model [4, 5, 7, 45, 46]. The modifications are typically obtained by a Yukawa-force correction, Eq. (13). Milgrom’s interpolation formulae for a single pair of baryonic objects leads to an asymptotic modification of the forces, given by Eq. (11). An example of the two modifications is shown in Fig. 1 for specific values of the constants for the modifications. Another attempt to modify the dynamics of galaxies is the f(R) [64], f(T) [65,66,67] , and f(Q) [68] theories where the modified behavior is obtained by modifying Einstein’s GR theory. For a review of the f(R) modifications see the reviews[16, 49]. Here I would like to propose another possibility, viz. gravitational lensing [19,20,21] and focusing, caused by heavy centers of mass in the galaxy of the gravitational waves from objects in the galaxy that are located at long distances from the object.