© The Authors 2024

1 Introduction

If stars and galaxies move under Newtonian gravity, dark matter is required to bind the stars to galaxies and galaxies in clusters. There is ample evidence. The mass inferred from the stellar velocities in clusters (with the virial equation; Zwicky 1933, 1937) and in galaxies (with Huygens’ equation; Rubin & Ford 1970; Rubin et al. 1980) is much higher than the mass of the luminous matter. ACDM cosmology requires six times as much matter as can be directly detected, hence, a large nonradiative matter component must dominate the Universe.

Despite a search of over half a century, no dark matter particle has been directly detected, while the list of astronomical observations that are difficult to explain with dark matter models has been growing: (i) The absence of cusps at galactic cores (McGaugh et al. 2003; de Blok 2010); (ii) bar structures at galactic cores (Chiba & Schönrich 2021; Roshan et al. 2021a,b); (iii) the low number of galaxy satellites (Bullock 2013); (iv) the absence of dynamical friction on galaxy satellites by its dark matter halo (Angus et al. 2011; Kroupa 2015; Oehm et al. 2017); (v) the shape and central brightness of the Fornax cluster (Asencio et al. 2022); (vi) the symmetry breaking of tidal tails of open star clusters (Kroupa et al. 2022); (vii) galaxies without dark matter (Moreno et al. 2022; Comerón et al. 2023); (viii) the intracluster gas does not cool (Fabian 1994); (ix) strong local lensing in galaxy clusters (Meneghetti et al. 2020); (x) the early formation of very large galaxies and clusters (Asencio et al. 2020, 2023; Labbé et al. 2023); (xi) the Hubble H₀ tension and the S₈ tension occur in standard ΛCDM cosmology (Planck Collaboration VI 2020; Haslbauer et al. 2020; Asgari et al. 2021; Riess et al. 2021); (xii) the violation of the strong equivalence principle (Blanchet & Novak 2011; Chae et al. 2020); (xiii) and the velocities in wide binary star systems (Banik & Zhao 2018; Pittordis & Sutherland 2019; Chae 2023; Hernandez 2023).

An alternative to dark matter is modified Newtonian dynamics (MOND), which was invented by Milgrom and Bekenstein (Milgrom 1983; Bekenstein & Milgrom 1984; Milgrom & Sanders 2008). Whereas Einstein’s general relativity modifies Newton when the gravity potential ϕ is deep, Milgrom’s MOND instead modifies Newton at low accelerations. The critical acceleration scale is close to the scale set by the cosmological constant Λ. It is a₀ ≈ .127⋅c²Λ^1/2. The strengths of the mutual gravity between stars, the centripetal acceleration of stars in the galaxy outskirts, and the accelerations in galaxy clusters all have values of about a₀. Because MOND is nonrelativistic, it applies in the regime in which |ϕ| ≪c² and |∇ϕ| ≲ c²Λ^1/2. The theory explains some of the listed problems and it gives a natural explanation for the observed flat rotation curves of galaxies as well as for the Tully-Fisher relation (McGaugh et al. 2000).

Although the MOND model makes definite predictions without the need for fit parameters, the model is complicated by the fact that it is not relativistic. It needs a preferred inertial reference frame (i.e., it is an aether theory) and the gravity from distant sources cause an external field effect (violating the strong equivalence principle). These effects need to be considered when the possibility of MOND is to be ruled out or confirmed.

General relativistic models that reduce to MOND at low energies introduce additional degrees of freedom. The new TeVeS theories of Skordis & Złośnik (2019, 2021) have resolved the long-standing problem that gravitational waves traveled faster than photons. This problem plagued older MOND-extensions (Bekenstein 2004).

Milgrom’s MOND modifies the equation for the gravitational potential using a nonlinear PDE. The numerical codes that have been developed to simulate the dynamical evolution in MOND are found in Brada & Milgrom (1999); Nipoti et al. (2007); Llinares et al. (2008), N-MODY by Londrillo & Nipoti (2009), and RayMOND by Candlish et al. (2014), based on multigrid discretization. The work of Angus & Diaferio (2011); Peng et al. (2016) solved for QuMOND, the quasilinear MOND model (Milgrom 2010). The code called Phantom of Ramses (Lüghausen et al. 2015; Nagesh et al. 2022) also solves QuMOND. The codes RayMOND and Phantom of Ramses are both based on RAMSES (by Teyssier 2002).

We describe an algorithm for N-body simulations in MOND. The particles can be Solar System planets, stars in a galaxy, or the galaxies in a galaxy cluster. The following main fields in astronomy are affected by Milgromian gravity:

• (i)

  Celestial mechanics. In the Solar System, the perturbation of the planetary orbits due to a possible residual MOND effect may have to be calculated. Because outer Solar System dynamics is in the crossover regime, this correction critically depends on the interpolation function, which describes the transition between MOND and normal gravity. In particular, the perihelion precession of the outer planets, where acceleration is lowest, could be affected (Pitjeva & Pitjev 2013; Paučo & Klačka 2016; Paučo 2017; Yoon & Darriba 2020). Here, the effect of MOND may lie just below the current detection limit. Binary star systems with a wide separation are in the transition of the interpolation (Hernandez et al. 2019). For both these systems a two-body model is appropriate when the external field from the Milky Way is taken into account (Iorio 2009, 2017).

• (ii)

  Dense stellar systems. These include open clusters and globular clusters, where close encounters between stars are relevant. The evolution requires inclusion of short-range gravitational interactions. Hence, individual particles must be tracked.

• (iii)

  Galactic dynamics. Examples are transient phenomena such as tidal effects between neighboring galaxies, the collisions between galaxies, the formation of tails, and the collapse or merging of smaller stellar systems. Galaxy-formation simulations in QuMOND were made by Wittenburg et al. (2020); Eappen et al. (2022).

• (iv)

  Galaxy clusters. Here, the particles are the individual galaxies. An N-body code for the dynamics with small particle numbers, between 10² and 10³, possibly supplemented with hydrodynamics for the gas component, would be the appropriate model.

• (v)

  Cosmology. Structure formation in the early Universe would be affected by MOND because the density fluctuations will initially be small, accelerations could begin with low values, and structure formation may start in the deep-MOND regime. Cosmological simulations in QuMOND were made by Katz et al. (2013); Haslbauer et al. (2020); Wittenburg et al. (2023).

2 Theory and main result: The idea of the method

As the particles move in space under the influence of the gravitational potential ϕ(r, t), they experience an acceleration given by Newton’s second law, $$m_i{\stackrel{\mathrm{..}}{r}}_i(t)=-m_i\text{∇}\varphi \left(r_i(t),t\right).$$(1)

The baryonic mass density for the set of N point masses is $${\rho }_{\text{B}}(r,t)={\displaystyle \sum _{i=1}^N{m}_i}\delta \left(r-r_i(t)\right).$$(2)

Here, δ is the three-dimensional Dirac-delta function. Whereas the potential is a solution of the Poisson equation in the Newton theory, in MOND (the AQUAL formalism), the potential is a solution to the following nonlinear partial differential equation (PDE): $$\text{∇}•\mu \left(\frac{\left|\text{∇}\varphi \right|}{a_0}\right)\text{∇}\varphi =4\pi G{\rho }_{\text{B}}(r,t),$$(3)

with the boundary condition ∇ϕ(r, t) → 0 as |r| → ∞. The scalar function µ(x) is called the interpolation function because it interpolates between the regime of high and low accelerations. At high acceleration, µ(x) → 1 so that the potential will approach the Newton potential. At low acceleration, the deep MOND regime begins, and µ(x) → x. An isolated and localized mass m has the asymptotic spherically symmetric potential $\varphi \left(r\right)=\sqrt{Gm{a}_0}$ log |r|. See Table A.1 for the symbols and notation. A popular choice for the interpolation function that is consistent with all observations so far (Brouwer et al. 2021) is $$\mu (x)=\frac{x}{\sqrt{1+x^2}}$$(4)

2.1 An equivalent system for AQUAL

We now demonstrate how to rewrite the Eq. (3) for the potential in MOND by transforming this equation into a system of coupled linear PDEs for four vector fields 𝑔_N, 𝑔_M, H, and F, combined with one algebraic equation. This system of equations is $$\{\begin{array}{lll}\text{∇}•g_{\text{N}}=-4\pi G{\rho }_{\text{B}},\hfill & \text{∇}\times g_{\text{N}}=0,\hfill & \left(5\right)\hfill \\ F=g_{\text{N}}+H,\hfill & \hfill & \left(6\right)\hfill \\ \hfill & \text{∇}•H=0,\hfill & \left(7\right)\hfill \\ g_{\text{M}}=v\left(\frac{F}{a_0}\right)F,\hfill & \hfill & \left(8\right)\hfill \\ \hfill & \text{∇}\times g_{\text{M}}=0,\hfill & \left(9\right)\hfill \\ F=\mu \left(\frac{g_{\text{M}}}{a_0}\right)g_{\text{M}}.\hfill & \hfill & \left(10\right)\hfill \end{array}$$(5)–(10)

The vector fields 𝑔_N and 𝑔_M are the Newton and MOND acceleration fields. With the identification 𝑔_N=−∇ϕ_N, the two equations in Eq. (5) are equivalent to the Poisson equation. Our new system effectively expresses the MOND gravity field 𝑔_M in the Newtonian gravity field 𝑔_N. Equations (8) and (10) are not independent equations, but are the inverse of each other. Here, F and 𝑔_M are the lengths of the vectors F and 𝑔_M. The pair µ(x) and v(y) relate the strength of the Newtonian gravity to the strength of the Milgromian gravity by $$\mu (x)x=y=\frac{F}{a_0},\text{\hspace{1em} and \hspace{1em}}v(y)y=x=\frac{g_{\text{M}}}{a_0}\text{. }$$

The choice Eq. (4) for the interpolation function gives $$v(y)=\sqrt{\frac{1}{2}+\frac{1}{2}\sqrt{1+\frac{4}{y^2}}}.$$

Here follows the proof that the system of Eqs. (5)–(10) is equivalent to Eq. (3). Instead of an equation for the potential ϕ, we have an equation for the acceleration field 𝑔_M. We thus need to identify 𝑔_M = −∇ϕ. When we now combine Eq. (6) with Eq. (10), solve for 𝑔_N, then substitute this solution in the first of Eq. (5), and use Eq. (7), we obtain $$-4\pi G{\rho }_{\text{B}}\stackrel{(5)}{=}\mathbf{\nabla }•g_{\text{N}}\stackrel{(6)}{=}\mathbf{\nabla }•F-\mathit{\nabla }•H\stackrel{(7)}{=}\mathbf{\nabla }•F\stackrel{(10)}{=}\mathbf{\nabla }•\mu \left(\frac{g_{\text{M}}}{a_0}\right)g_{\text{M}}.$$

Equation (9) states that the field is curl-free. Hence, 𝑔_M is a conservative field with a potential. When we substitute 𝑔_M = −∇ϕ, we find Eq. (3). This proves that the system Eqs. (5)–(10) implies Eq. (3). For the converse, we note that when we have a solution of Eq. (3), we can calculate 𝑔_M. This automatically satisfies Eq. (9). We can then define F with Eq. (10). We then let Eq. (6) be its Helmholtz decomposition. This defines the functions 𝑔_N and H, which must satisfy the second of Eqs. (5) and (7), stating that the first component is curl-free and the second component is divergence-free. This completes the proof.

In any numerical code that uses this system, the linear PDEs can be solved using standard numerical methods for linear equations. The nonlinear equation requires a special solver. We have devised the following numerical iteration scheme that approaches the solution of the system (5)–(10). We first solve Eq. (5) for the Newtonian gravity field 𝑔_N, which we store in memory. Then we take H = 0 for the initial value. The iteration scheme consists of the following steps:

• 1.

  Calculate F with Eq. (6).

• 2.

  Calculate g_M with Eq. (8).

• 2.

  Project 𝑔_M → 𝑔_M|| onto the curl-free part, as per Eq. (9).

• 4.

  Calculate F with Eq. (10).

• 5.

  Calculate H with Eq. (6).

• 6.

  Project H → H_⊥ onto the divergence-free part, as per Eq. (7).

• 7.

  Repeat: Go to step 1.

The calculation in steps 3 and 6 can quickly be made in the Fourier domain (see Appendix C for the equations). The use of the fast Fourier transform (FFT) allows us to perform these steps fast.

The iterative procedure converges to a numerically stable solution. Because we start with a zero H field, its change after the iterations allows us to estimate the error in the numerical solution. It was found that the error initially drops with a factor of 10 per iteration, but after about 15 iterations, 5 iterations are required for the error to drop with a factor 10, until machine precision is reached (see Platschorre 2019, where it was used for the first time). One clue to the explanation for the fast convergence is the following. For a spherical field, we have H = 0 (see Bekenstein & Milgrom 1984). Because the masses are spherical, close to a mass, the solution is approximately spherically symmetric. Far away from a collection of masses, the solution becomes spherical as well. Thus, near the masses and away from the masses, the initial value is correct. In the iterative steps, the field in the intermediate regions is changed. Here, H is not zero.

The apparent matter (baryonic plus dark matter) distribution can be calculated by substituting the correct potential (the potential that is consistent with the observed particle motions and lensing) in the incorrect Newton-Poisson equation. This is equivalent to $${\rho }_{\text{B}}(r)=\frac{-1}{4\pi G}\mathbf{\nabla }•g_{\text{N}},{\rho }_{\text{A}}(r)=\frac{-1}{4\pi G}\mathbf{\nabla }•g_{\text{M}}.$$(11)

The apparent density ρ_A is also known as the effective density (see Banik & Zhao 2022). This leads to an apparent dark matter density of ρ_A − ρ_B, also called phantom dark matter density. The total baryonic mass and apparent mass in a volume W can be found from the flux through its boundary, $$M_{\text{B}}(W)=-{\displaystyle ∯{}_{\partial W}{g}_{\text{N}}•\text{d}S,M_{\text{A}}(W)=-}{\displaystyle ∯{}_{\partial W}{g}_{\text{M}}•\text{d}S.}$$(12)

2.2 An equivalent system for QuMOND

An alternative to the AQUAL formalism of MOND is quasilinear MOND or QuMOND (Milgrom 2010; Candlish et al. 2014; Lüghausen et al. 2015; Nagesh et al. 2022). This model can be cast as a similar system $$\{\begin{array}{ccc}\mathbf{\nabla }•g_{\text{N}}=-4\pi G{\rho }_{\text{B}},& \mathbf{\nabla }\times g_{\text{N}}=0,& \\ g_{\text{M}}+H=v\left(\frac{g_{\text{N}}}{a_0}\right)g_{\text{N}},& \mathbf{\nabla }\times g_{\text{M}}=0,& \mathbf{\nabla }•H=0.\end{array}$$

We do not consider QuMOND any further. However, the algorithm can easily be adapted using the same steps.

3 The algorithm: Dealing with the self-gravity

We now describe the implementation of the method. The code needs to numerically approximate the solution of the MOND PDEs in Eqs. (5)–(10) and propagate particle motion and a time-dependent mass density. The particle equations of motion Eq. (1) are implemented with Leapfrog, $$v_i\to v_i+g_i(t)\Delta t,r_i\to r_i+v_i\Delta t,i=1\dots N.$$(13)

Although the formulas do not explicitly show this, the velocities and positions represent the values at a half time step difference. In principle, an actual half time step needs to be made initially (and also finally) to obtain the true Leapfrog integration scheme.

The code uses a discrete Fourier transform, so that space is represented as a periodic cubic lattice of n³ cells we call pixels. Hence, the scheme does not have nested grids; the grid cells are all equal and do not change. However, due to the periodic nature of the simulated space, matter inside one cubic volume is effectively influenced by the gravity of the neighboring copies. This effect is nonphysical, and a volume of empty space needs to surround the matter when this effect is to be small. Ideally, the center of mass is placed in the origin, and it is ensured that particles stay away from the boundary of the cube.

3.1 Interpolation for the particle density function

The main drawback of the particle-mesh (PM) method is the lack of accuracy at the short-distance scale of the pixels. Although each particle is located in a single pixel, we desire accuracy that is below the grid size. This influences both the field created by and the acceleration on the particle. First, the field created by particle i close to the mass should resemble a spherical field about the point r_i, which is not the center of a pixel. Second, we also need to approximate the acceleration on particle i by the field at r_i, for which the value at the center of the pixel is not accurate.

We used a Gaussian kernel to create a smoothed density in order to suppresses anti-aliasing. The density of Eq. (2) was replaced by $${\rho }_{\text{B}}(r)=\frac{1}{s^3{(2\pi )}^{3/2}}{\displaystyle \sum _{i=1}^N{m}_i}\exp \left(\frac{-{\left|r-r_i\right|}^2}{2s^2}\right).$$(14)

Hence, each of the N point particles is represented by the Gaussian in Eq. (14) centered about the particle position. We chose a standard deviation with a width of one pixel s. The Gaussian is then effectively a cloud-in-cell (CIC) method of order 12 (see Hockney & Eastwood 1988). This is explained in Appendix B. Because Gaussians drop off fast, we implemented the calculation of the density by only addressing the pixels that are within a distance of 4s from the central pixel containing the particle. To do this, we stored a list with the 251 coordinates of the ball-shaped domain with a radius of 4s. This is illustrated in Fig. 1. Moreover, we created a lookup table for the exponential function.

Fig. 1 Anti-aliasing needed to suppress the self-gravity. A cross section of the Gaussian density Eq. (14) with one-pixel variance about (blue) points in the middle grid cell, representing a 12th order weight function. We only account for the pixels inside a sphere with a radius of 4s (dashed circles). This results in a domain of 251 cells (inside the contours). The linear motion of a point particle across the middle cell (left to top right) shows that the simulated density changes smoothly between panels.

3.2 Interpolation for the forces on the particles

Next, we describe how we approximated the acceleration of particle i. Again, we wished to know the acceleration at a point r_i inside some pixel with subpixel accuracy. However, the acceleration 𝑔_B is only calculated at the pixel positions. We would therefore need to interpolate between the values at different pixels. For the force on particle i, located at r_i, we used the following smoothed average: $$g_i=\frac{-1}{s^3{(2\pi )}^{3/2}}{\displaystyle \iiint \mathbf{\nabla }}\varphi (r)\exp \left(\frac{-{\left|r-r_i\right|}^2}{2s^2}\right)\text{d}x \text{d}y \text{d}z.$$(15)

On the mesh, the gradient of the potential ∇ϕ in this expression is calculated with the finite-difference gradient, and the integral in this expression becomes a pixel sum. We only summed over pixels within a distance 4s from the pixel containing r_i.

The smoothing methods discussed in this section will generally improve the accuracy in any PM method. However, in the case of MOND, this smoothing is essential for dealing with the self-gravity. The software is available on GitHub¹ under the MIT license and is archived in the ASCL (de Nijs et al. 2023).

4 Results: Validation using analytic solutions

The problem of self-gravity dominating the actual force on a point particle is inherent to the PM method for MOND. The effect of an imprecise subpixel resolution for the position of a particle and the consequently imprecise evaluation of the acceleration on that particle is always present, even for a single particle. We therefore also need to test the accuracy of the method for linear motion. We have found that in our simulations, a point particle of mass m experiences nonphysical accelerations with the average value of 1.08⋅${10}^{-3}\sqrt{Gm{a}_0}/s$ for a 256³ grid (see de Nijs 2023). This small error is independent of the number of pixels.

In this section, we present results of simulations with the code for three systems: the two-body problem, a ringed system, and an isotropic three-dimensional system. Because the main interest is in MOND, we considered the deep MOND case. This means that the interpolation functions we used in our simulation are $$\mu (x)=x,v(y)=\frac{1}{\sqrt{y}}.$$

For this regime, Milgrom (1994, 1997, 2014) found, using the method of the virial, that the kinetic energy for a bound system with the center-of-mass motion at rest is $$E_{\text{kin}}=\frac{\overline{M\mathit{v}•\mathit{v}}}{2}=\frac{M\overline{v^2}}{2}=\frac{\sqrt{G{a}_0}}{3}\left(M^{3/2}-{\displaystyle \sum _i{m}_i^{3/2}}\right),$$(16)

with m_i the individual particle masses, and M the total mass. If there is no external field, this single equation allows us to calculate the MOND forces between particles for highly symmetric cases, thereby bypassing the nonlinear PDE. We used these analytic solutions to test our PM code.

4.1 Example I: Wide binary system

We considered the two-body problem in MOND. This is a model for a binary stellar system with two comparable masses with a wide separation. The interparticle force is found from Eq. (16) to be $$F_{12}=\frac{2\sqrt{G{a}_0}}{3r}\left(M^{3/2}-m_1^{3/2}-m_2^{3/2}\right),$$(17)

with M = m₁ + m₂, and where r = |r₂ − r₁| is the separation. We can rewrite the mass dependence using $$\begin{array}{l}{M}^{3/2}-m_1^{3/2}-m_2^{3/2}=m_1\left(\sqrt{M}-\sqrt{m_1}\right)+m_2\left(\sqrt{M}-\sqrt{m_2}\right)\\ =m_1{m}_2\left(\frac{1}{\sqrt{M}+\sqrt{m_1}}+\frac{1}{\sqrt{M}+\sqrt{m_2}}\right).\end{array}$$(18)

The equations of motion of the two particle system are $$\{\begin{array}{l}{\stackrel{.}{r}}_1=v_1,\hfill \\ {\stackrel{.}{r}}_2=v_2,\hfill \\ {\stackrel{.}{v}}_1=\frac{2\sqrt{G{a}_0}}{3}\frac{M^{3/2}-m_1^{3/2}-m_2^{3/2}}{m_1}\frac{r_2-r_1}{{\left|r_1-r_2\right|}^2},\hfill \\ {\stackrel{.}{v}}_2=\frac{2\sqrt{G{a}_0}}{3}\frac{M^{3/2}-m_1^{3/2}-m_2^{3/2}}{m_2}\frac{r_1-r_2}{{\left|r_1-r_2\right|}^2}.\hfill \end{array}$$(19)

A numerical solution of this system of ordinary differential equations (ODE) is shown in Fig. 2. Here, it is compared with the results of our PM code, showing that the two methods agree reasonably well, but diverge ultimately. The system has the following constants of motion: $$\begin{array}{l}P=m_1{\upsilon }_1+m_2{\upsilon }_2,\hfill \\ R(t)-\frac{P}{M}t=\frac{m_1{r}_1+m_2{r}_2}{M}-\frac{P}{M}t,\hfill \\ r(t)\times P,\hfill \\ L=m_1{r}_1\times {\upsilon }_1+m_2{r}_2\times {\upsilon }_2,\hfill \\ E=\frac{m_1{v}_1^2+m_2{v}_2^2}{2}+2\sqrt{G{a}_0}\frac{M^{3/2}-m_1^{3/2}-m_2^{3/2}}{3}\log \left|r_1-r_2\right|.\hfill \end{array}$$

Because the PM is cubic, which breaks the local spherical symmetry, we expect that the total angular momentum does not remain constant in the simulations. How well angular momentum is conserved by our PM code is shown in Fig. 3. If the center of mass is located at the origin, R = 0, the single-particle angular momenta L₁ and L₂ are also constant. However, the numerical error in L₁ and L₂ becomes large when the particle is within a distance of one pixel from the origin. The reason is that the center of mass may have drifted away from the origin by one pixel. When we subtract R from r₁ and r₂, the errors disappear. For grid sizes 128³, 256³, and 512³, the error in total angular momentum is 4%, 2%, and 1%, respectively.

The solution with the two particles in circular orbit about the center of mass is $$\begin{array}{ll}{r}_1=\frac{m_{2}r}{M}(i\cos \omega t+j\sin \omega t),\hfill & r_2=\frac{-m_{1}r}{M}(i\cos \omega t+j\sin \omega t),\hfill \\ {\upsilon }_1=\frac{-m_{2}v}{M}(i\sin \omega t-j\cos \omega t),\hfill & {\upsilon }_2=\frac{m_{1}v}{M}(i\sin \omega t-j\cos \omega t),\hfill \end{array}$$

where r = |r₂ − r₁| = r₁ + r₂ is the separation, υ = |υ₂ − υ₁| is the relative velocity, and ω is the orbital frequency. With Huygens’ law of the centripetal force applied to the motion of the reduced coordinate r, we obtain for the force $$\frac{m_1{v}_1^2}{r_1}=\frac{m_2{v}_2^2}{r_2}=\frac{m_1{m}_2{v}^2}{Mr}=F_{12},$$

and for the kinetic energy $$\frac{m_1{v}_1^2+m_2{v}_2^2}{2}=\frac{m_1{m}_2{v}^2}{2M}=\frac{M\overline{v^2}}{2}=E_{\text{kin}}.$$

Using Eqs. (16)–(18), we find that the velocity υ of the reduced motion and the orbital frequency ω are given by the expression $$v^2={(\omega r)}^2=\frac{M^2\overline{v^2}}{m_1{m}_2}=\frac{2\sqrt{GM{a}_0}}{3}\left(\frac{1}{1+\sqrt{m_1/M}}+\frac{1}{1+\sqrt{m_2/M}}\right).$$(20)

The velocities are independent of the separation. The theoretical dependence of the velocity on the mass ratio is plotted in Fig. 4. Figure 5 shows the numerical simulation of our PM code for two particles with the initial relative velocity from Eq. (20). It convincingly shows that motion is circular, which verifies Milgrom’s two-particle force law. It also illustrates that we can obtain reasonable subpixel resolution.

Fig. 2 Simulated orbits of the two-body system projected onto the xy plane for the mass ratio m1/m2 = 3/2. The mesh has a size of 1283, and there are 100 time steps, each with four iterations for the MOND force. Red (m1) and cyan (m2) are the orbits obtained from numerical integration of Eqs. (19). Purple and blue are the orbits calculated with the PM code. The difference between the PM code and the ODEs is a slow dephasing of the orbits.

Fig. 3 Angular momentum vs. time t for the simulation shown in Fig. 2. These are constants of motion of the ODE system Eq. (19) and numerically constant (with machine precision, not shown) when integrated with Leapfrog. In the PM code, the total angular momentum (blue) is well conserved, but the single-particle angular momenta (orange and red) have large errors when the distance to the origin is smaller than one pixel.

Fig. 4 Force times radius (blue) and squared velocity (orange) vs. mass ratio for the two-body system. The horizontal axis is logarithmic, making the curves symmetric about m1 /m2 = 1.

4.2 Example II: Ring galaxy

Again, we tested the code with a system that was solved analytically by Milgrom (1994): a central mass m₀, surrounded by a ring with a radius r₂ with N particles with a mass m (i.e., the system has N + 1 particles). A beautiful realization of such a ring galaxy is Hoag’s Object (Schweizer et al. 1987). The central mass is at rest at the origin. We created a rotating ring of N particles on a regular N-gon (N ≥ 2), $$\begin{array}{l}{r}_i=r_{2}i\cos (\omega t+2\pi i/N)+r_{2}j\sin (\omega t+2\pi i/N),\hfill \\ {\upsilon }_i=-vi\sin (\omega t+2\pi i/N)+vj\cos (\omega t+2\pi i/N).\hfill \end{array}$$

The force on the particles in the ring is given by $$F_2=\frac{2\sqrt{G{a}_0}}{3N{r}_2}\left(M^{3/2}-m_0^{3/2}-N{m}^{3/2}\right).$$(21)

Here, M = m₀ + N_m. We can now write $$\begin{array}{cc}{M}^{3/2}-m_0^{3/2}& =m_0\left(\sqrt{M}-\sqrt{m_0}\right)+Nm\sqrt{M}\\ & =MNm\left(\frac{1}{\sqrt{M}+\sqrt{m_0}}+\frac{\sqrt{m_0}}{M}\right).\end{array}$$

The velocity of the ring particles and the orbital frequency are found by equating the force Eq. (21) to F₂ = mυ²/r₂ = mω²r₂, $$v^2=\frac{M\overline{v^2}}{Nm}=\frac{2\sqrt{GM{a}_0}}{3}\left(\frac{1}{1+\sqrt{m_0/M}}+\sqrt{m_0/M}-\sqrt{m/M}\right).$$(22)

The force times radius and squared velocity is plotted in Fig. 6. The result of the simulation for N = 100 particles is plotted in Fig. 7. If the mass of the ring is not very low compared to the central mass, the system is unstable, and the ring breaks apart after a few orbits.

Fig. 5 Simulated orbits of two bodies with a mass ratio m1/m2 = 3/2 as in Fig. 2, but with 400 time steps. The initial velocities are given by Eq. (20) and should result in circular orbits. These are plotted in Fig 4. This simulation thus confirms the force formula Eq. (17) for the deep MOND case and validates the algorithm. It also illustrates how well the algorithm performs on the subpixel scale.

4.3 Eample III: The isothermal sphere in MOND

Our third application is the simulation of a spherically symmetric system in thermodynamic equilibrium. This could be a globular cluster with up to 10⁶ stars or a galaxy cluster with up to 10³ galaxies. These system are in the deep MOND regime or in the crossover regime (accelerations around a₀)·

A system of N point particles with equal mass m in hydrostatic equilibrium with a temperature T has a velocity dispersion given by $$\frac{3k_{\text{B}}T}{2m}=\frac{\overline{v•v}}{2}=\frac{\overline{v^2}}{2}=\frac{\sqrt{GM{a}_0}}{3}.$$(23)

The last equality is the thermodynamic limit of Eq. (16). For a spherical system, the particle number density in terms of the potential is given by $$\frac{n(r)}{n(0)}=\exp \left(\frac{-m\varphi (r)}{k_{\text{B}}T}\right)=\exp \left(\frac{-3\varphi (r)}{\overline{v^2}}\right)=\exp \left(\frac{-9\varphi (r)}{2\sqrt{GM{a}_0}}\right).$$(24)

Here, we used hydrostatic equilibrium and the ideal gas law. Substitution of the mass density ρ_B,(r) = mn(r) and Eq. (24) in the MOND Eq. (3) gives the following equation for the potential: $$\frac{1}{a_0{r}^2}\frac{\text{d}}{\text{d}r}{\left(r\frac{\text{d}\varphi }{\text{d}r}\right)}^2=4\pi G{\rho }_{\text{B}}(0)\exp \left(\frac{-9\varphi (r)}{2\sqrt{GM{a}_0}}\right).$$(25)

The solution of this equation is found to be $$\varphi (r)=\frac{2\sqrt{GM{a}_0}}{3}\log \left(1+\frac{r^{3/2}}{b^{3/2}}\right),$$(26) $$g_{\text{M}}(r)=-r\sqrt{\frac{GM{a}_0}{b^{3}r}}{\left(1+\frac{r^{3/2}}{b^{3/2}}\right)}^{-1},$$(27) $$M(r)=\frac{M{r}^3}{b^3}{\left(1+\frac{r^{3/2}}{b^{3/2}}\right)}^{-2},$$(28) $${\rho }_{\text{B}}(r)=\frac{3M}{4\pi b^3}{\left(1+\frac{r^{3/2}}{b^{3/2}}\right)}^{-3}.$$(29)

Here, b is a parameter with the dimension of a length. When the function in Eq. (26) is substituted in Eq. (25), we can solve for this parameter b in terms of the mass M and the central density ρ_B(0). Hence, we have found a two-parameter family of solutions. The function M(r) in Eq. (28) represents the mass inside a sphere of radius r, and M = lim_r→∞ M(r) is the total mass. The physical meaning of the parameter b becomes clear when we calculate the mass of a homogeneous sphere with radius b, with the constant density equal to the value ρ(0) at the center of the isothermal distribution, given in Eq. (29). This mass is precisely equal to the total mass $$M=4\pi \underset{0}{\overset{\infty }{{\displaystyle \int }}}{\rho }_{\text{B}}(r)r^2\text{d}r=4\pi \underset{0}{\overset{b}{{\displaystyle \int }}}{\rho }_{\text{B}}(0)r^2\text{d}r=\frac{4\pi b^3}{3}{\rho }_{\text{B}}(0).$$

The density is finite at the origin and also drops off sufficiently fast to zero to allow normalization: No cutoff is needed to keep the total mass finite. These physically convenient properties are absent from the isothermal solution in Newtonian gravity. The potential and density and the resulting deflection potential and surface-mass density are plotted in Fig. 8. Kinetic energy and potential energy are, with Eq. (23), $$E_{\text{kin}}=\frac{M\overline{{\upsilon }^2}}{2},E_{\text{pot}}=4\pi \underset{0}{\overset{\infty }{{\displaystyle \int }}}{\rho }_{\text{B}}(r)\varphi (r)r^2\text{d}r=\frac{3M\overline{{\upsilon }^2}}{2}.$$

The solution given in Eqs. (26)–(29) is applicable for 𝑔_M ≪ a₀· Because the maximum acceleration is found at r = b/2^2/3, the requirement becomes GM/b² ≪ 9a₀/2^4/3. The Virgo cluster, for example, satisfies this. Equation (29) could be a realistic description for spherical gas clouds, elliptical galaxies, and clusters.

We derived the formula for the apparent mass density, which would give the apparent dark matter distribution function. Equations (11)–(12) give $$\begin{array}{l}{M}_{\text{A}}(r)=\sqrt{\frac{M{a}_0}{G}}r{\left(1+\frac{b^{3/2}}{r^{3/2}}\right)}^{-1},\hfill \\ {\rho }_{\text{A}}(r)=\frac{1}{4\pi }\sqrt{\frac{M{a}_0}{G{b}^{3}r}}\left(\frac{5}{2}+\frac{r^{3/2}}{b^{3/2}}\right){\left(1+\frac{r^{3/2}}{b^{3/2}}\right)}^{-2}.\hfill \end{array}$$

The total apparent mass diverges as r →∞. Whereas the actual density ρ_B is smooth at the center, the apparent mass density ρ_A has a cusp at the core. The isothermal solution shows typical deep MOND behavior: Only the total baryonic mass M determines the velocity distribution, with a central density that can be varied independently. The lower the density at the center, the larger the system system size and the higher the apparent dark matter mass.

We now study the dynamics of the isothermal solution with our PM code. We need the cumulative distribution function for the radial probability density. According to Eq. (28), this is $$\frac{M(r)}{M}={\left(1+\frac{b^{3/2}}{r^{3/2}}\right)}^{-2}.$$

In order to obtain random realizations, our point masses were found from seven random numbers ξ_ı ∈ (−1,1), η₁,η₂,η₃ ∈ (0,1), ζ₁,ζ₂,ζ₃ ∈ (0,2π), sampled homogeneously from these intervals and the formulae $$\begin{array}{l}{r}_i=b\left(\sqrt{1-{\xi }_1^2}\left(i\cos {\zeta }_1+j\sin {\zeta }_1\right)+k{\xi }_1\right){\left(\frac{1}{\sqrt{{\eta }_1}}-1\right)}^{-2/3},\\ v_i=\sqrt{-\frac{2}{3}\overline{v^2}\log {\eta }_2}\left(i\cos {\zeta }_2+j\sin {\zeta }_2\right)+\sqrt{-\frac{2}{3}\overline{v^2}\log {\eta }_3}k\cos {\zeta }_3.\end{array}$$(30)

The Box-Muller transform was used to generate the Gaussian values of the Maxwell-Boltzmann velocity distribution. Figures 9 and 10 show how close the realizations are to the cumulative distributions for the radial coordinate and for the velocity. Figure 11 shows the realizations projected onto the celestial plane and the resulting apparent dark matter surface density. For N = 10³ or larger, the potential is nearly spherical, and it is well described with our analytic solution.

We solved the evolution of N particles in the external potential Eq. (26) as a system of ODEs, and we simulated the N particles in MONDian gravity using the PM code, both with the same initial distribution. The two models are comparable because both must behave identical in the thermodynamic limit N → ∞. We also tested whether the distributions are stationary, as they should be for large N. In the external potential, the energies and angular momenta of the individual particles are strictly conserved, and so is the total energy and angular momentum. The single-particle kinetic energies and momenta oscillate, and so does the total kinetic energy and the total momentum. We compare this with the exact system of Eqs. (1)–(3). Here, the total energy, total momentum, and total angular momentum are conserved. The PM code simulates this system. The fluctuations in the energies are shown in Fig. 12. Because the mesh in our PM approach breaks spherical symmetry, angular momentum also fluctuates slightly.

Fig. 6 Force times radius (lower curves) and squared velocity (higher curves) vs. the mass ratio for particles on a regular N-gon orbiting a central body in circular motion. From Eqs. (2l)–(22), where m0 is the central mass and Nm is the ring mass.

Fig. 7 Simulated orbits for a ring of N = 100 particles with a mass m around a central mass m0 for a mass ratio of m0/Nm = 3/2 as in Fig. 2. The orbits of the first 80 time steps of 10 particles are shown in color. In gray we plot 576 steps for all particles. The initial velocities are given by Eq. (22), plotted in Fig. 6, which in theory results in circular orbits. Although this confirms the force formula Eq. (21), numerical round-off errors on the subpixel scale propagate and cause the orbits to diverge.

Fig. 8 Potentials and mass distributions for the isothermal state. In teal and orange we show the potential ϕ(r) as given by Eq. (26) and the deflection potential ψ=∫ϕdz for the line of sight passing the center at distance r, determining weak gravitational lensing. The asymptotes are shown (in gray) and the units are $\sqrt{GM{a}_0}{b}^j$, j = 0,1. In blue and red we show the density ρB(r) as given by Eq. (29) and the surface mass density σ=∫ρBdz in units of 3M/4πbj, j = 3,2.

Fig. 9 Cumulative mass M(r) of the spherical isothermal galaxy vs. r. The (free) parameter b determines the galaxy size. The (blue) curve is the theoretical result given by Eq. (28). The realizations (orange, red, and purple) are calculated from Eq. (30).

Fig. 10 Cumulative velocity distribution for the spherical isothermal galaxy. This is the fraction #(≤v)/N of particles with velocities below v vs. the value v of this upper bound. The realizations (orange, red, and purple) are calculated from Eq. (30).

5 Discussion

5.1 Accuracy of the method

The inaccuracies in the algorithm arise from the fact that we used a fixed grid. The discrete grid has n pixels in each dimension with a pixel width s and has a finite volume n³s³.

The use of FFTs in the calculation of the potentials and acceleration fields generally gives rise to anti-aliasing: The functions obtain nonphysical short-wavelength oscillations. This effect is particularly strong when the point particles are represented with isolated nonzero pixel values in the mass density, but also when the CIC method is used. Smoothing with a minimum variance of one pixel width s in a spherical domain of radius 4 s needs to be used to overcome this problem. Although the discretization creates a minimum size, we reach a subpixel spatial resolution by using interpolation. We found that the retrieval of the forces on the particles with a Gaussian smoothing gives the most accurate results. To further limit the aliasing and increase accuracy, we did not use the numerical acceleration field, but calculated the potential, and we derived the forces instead using the finite-difference gradient of this scalar function.

The runtime of our code has a time complexity of O(Nn³ log n), hence it scales linearly in the number of particles and almost linearly in the number of pixels because we made use of the FFT. The runtime as a function of particle count is plotted in Fig. (13). For a 128³ grid, the evaluation of the accelerations or/and densities starts to dominate the runtime at N = 10³. Therefore, in the general case, we expect that the crossover occurs at a critical density of 10³/128³ ≈ 1/10³ particles per pixel.

5.2 Improvements of the code

The code could be made faster, allowing simulations with higher resolution/larger volumes, with smaller time steps/longer simulated times, or with more particles. We list below a few ideas for possible improvements that we did not implement.

• (i)

  The code is currently written in Python, using the FFTW library, which is a C library for the FFTs. The code could be entirely rewritten in a compiler programming language such as C, and compiled into fast-running machine code.

• (ii)

  A primitive data type can be chosen for the floats with half-precision, used by GPUs (with the cuFFTW library). The machine precision would be 10⁻³, which is close to the unavoidable error due to self-gravity. Its implementation should reduce memory use and speed up the code.

• (iii)

  Both the potential solver and the particle propagator in the code can be parallelized. The calculation of the potential requires 42 FFTs (3 scalar plus 13 vector FFTs). These three-dimensional FFTs need to be performed sequentially. For each direction, we have n² normal FFTs of a vector with n elements. These can be performed in parallel. The increments of positions and velocities of the N particles are independent, hence their calculations can also be made in parallel. Parallelizing the updating of the density is less straightforward. If the distance between two particles is smaller than 8 pixels, the density needs to be updated one particle at a time.

• (iv)

  If the density function changes only slightly at each time step, the potential solver needs just one iteration when it is seeded with the previous solution. This strategy can be employed when the spatial steps are smaller than the width of the weight functions (i.e. for small time steps), or when nearby particles have mostly overlapping weight functions (i.e. for high number densities). Densities are also smooth when the matter is modeled as a fluid. We found that the PDE solver indeed sped up by a factor three.

• (v)

  If changes in the fields at the time steps are small, we can calculate the increments $\text{Δ}{\rho }_{\text{B}},\text{Δ}{g}_{\text{N}},\text{Δ}{g}_{\text{M}},\text{Δ}F$, ∆F, and ∆H, instead of the full fields. This will improve the accuracy and reduce the required number of significant digits so that half-precision can be used. For the nonlinear Eqs. (8) and (10), we can use the linearizations

$$\begin{array}{l}\text{Δ}{g}_{\text{M}}(r)=v\left(\frac{F}{a_0}\right)\text{Δ}F+v^{\prime }\left(\frac{F}{a_0}\right)\frac{F•\text{Δ}F}{a_{0}F}F,\\ \text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}Δ}F(r)=\mu \left(\frac{g_{\text{M}}}{a_0}\right)\text{Δ}{g}_{\text{M}}+{\mu }^{\prime }\left(\frac{g_{\text{M}}}{a_0}\right)\frac{g_{\text{M}}•\text{Δ}{g}_{\text{M}}}{a_0{g}_{\text{M}}}{g}_{\text{M}}.\end{array}$$

5.3 Applications targeted to special systems

We validated the PM code by simulating the two-body system and ring systems. For applications of systems with a high symmetry like this, however, the code could be adapted to make it more accurate and much faster. For example, a binary system can be simulated using ODEs and does not require a PDE solver, as shown in Sect. 4.1 for the deep MOND system Eq. (19). For the general case, the interpolation function μ is different, but the system is still governed by a two-body force that is a function of distance and the masses. One could now use the PDE solver of our code to calculate the force. This numerical force function can then be used in the ODE system for the actual simulations of the binary.

In absence of an external field, the two-body system has a perfectly cylindrical symmetry: The potential is a function of the distance to and height on the axis of symmetry. By writing the PDE in cylindrical coordinates and using the Hankel transform, we limit the number of transforms in the PDE solver. The evaluation of the force function could thus be calculated fast, allowing one to test various choices of interpolation functions. For other systems with cylindrical symmetry, such as the motion of a test particle in a disk galaxy, this method can also be used to enhance the accuracy and speed (see Platschorre 2019, for a study of rotation curves using the Hankel transform in this way).

The external field can be included in our PM code by fixing a plane at maximum distance from the particles in the cube where a permanent dipolar mass distribution is placed, that is, a layer of pixels with mass M and an adjacent pixel layer of mass −M.

When the external field is included, one can again evaluate the force function for the two-body problem numerically. Now, the force also depends on the polar angle (the angle between the direction of the external field and the vector r₂ − r₁), and not just on the interparticle distance.

Fig. 11 Surface-mass densities σ(x,y) and deflection potentials ψ(x,y) for a globular cluster from the central 128 × 128 part of a simulated cube of 2563 pixels. The yellow dots are the stars that define the actual density. The resulting apparent dark matter surface density is plotted in blue according to Eq. (11). Negative value are shown in red (Milgrom 1986). They are found for N ≤ 102. The equipotential curves of the deflection potential for Newtonian and Milgromian gravity are plotted in green and cyan. Although the densities are far from circularly symmetric, the deflection potential is nearly circular for N ≥ 103. This indicates that the retrieval of a clumpy mass distribution from weak lensing may be problematic for large N.

Fig. 12 Energies vs. time for N = 100 particles on a 2563 grid. The kinetic, potential, and total energy for the self-consistent interacting system of the particles from the PM code are shown as solid curves (red, purple, and blue). The kinetic, potential, and total energy from the ODE system of the particles in the external potential are shown as dashed curves (orange, violet, and teal). The fluctuations in the total energies are purely numerical, and the other energies have thermal fluctuations. The dashed gray curve is Epot = ∫ ρBϕdxdydz with the PM matter density and Eq. (26), hence substituting a finite N solution in the potential obtained in the thermodynamic limit. We subtracted the self-energies.

Fig. 13 Runtime vs. particle count N for a 1283 grid with the Python code on an Intel Core i7-9750H. The time spend on the FFTs (blue) is roughly constant, and the time for the evaluation of the particle accelerations and/or densities grows linear with N (orange, red, and purple). The total runtime shows a crossover around 103 particles (black). The dashed lines are guides with slopes of zero and one.

5.4 Extensions and future research directions

Extensions to the code could make the code useful for more realistic modeling and other fields. Some examples are listed below.

• (i)

  Close encounters. If the ratio Gm/s² is much higher than a₀, the field on the pixel scale is Newtonian. It may be possible to obtain subpixel resolution for this case by calculating the near field with a particle-particle component.

• (ii)

  Hydrodynamics. In galaxy clusters, the intergalactic gas component dominates the total mass contribution from the galaxies. One could implement Euler’s equations for inviscid flow inside the MOND acceleration field. The evolution of the density ρ_B, pressure ρ, and flow field u can be calculated from the system $$\begin{array}{l}\{\begin{array}{l}{\dot{\rho }}_{\text{B}}=u•{ℱ}^{-1}\left(\text{i}k{\widehat{\rho }}_{\text{B}}\right)+{\rho }_{\text{B}}{ℱ}^{-1}(\text{i}k•\widehat{u}),\hfill \\ \dot{p}=u•{ℱ}^{-1}(\text{i}k\widehat{p})+\gamma p{ℱ}^{-1}(\text{i}k•\widehat{u}),\hfill \\ \dot{u}={ℱ}^{-1}(\text{i}k\widehat{ℰ})-u\times {ℱ}^{-1}(\text{i}k\times \widehat{u})+{ℱ}^{-1}(\text{i}k\widehat{p})/{\rho }_{\text{B}}+g_{\text{M}}.\hfill \end{array}\\ \\ \end{array}$$

• Here, $\widehat{ℰ}=ℱ(u•u)/2$ is the transformed specific kinetic energy, and γ is the polytropic constant. These are the Euler equations in Lamb form, where the derivatives are calculated using Fourier transforms. By running these on a GPU with cuFFTW in half-precision floats, one gains a speed-up of at least one order of magnitude.

• (iii)

  Cosmology. In simulations of structure formation in the early Universe, one can introduce comoving coordinates to model the expansion of the Universe. This amounts to the equations for particle motion (Peebles 1993) $${\dot{r}}_i=v_i/a,{\dot{v}}_i=g_{\text{M}}\left(r_i\right)-v_i\dot{a}/a,$$

• and the density in Eq. (5) can be replaced with $a{\rho }_{\text{B}}(r)-a\overline{{\rho }_{\text{B}}}$, where the scale factor of the expanding Universe is given by a(t).

6 Conclusions

We described, coded, and tested an algorithm for the simulation of the dynamical evolution of N bodies for AQUAL-MOND. MOND is an alternative model for the gravity between stars, galaxies, and galaxy clusters, where standard theory requires CDM. We now list our final conclusions:

• (i)

  The MOND Eq. (3) is shown to be equivalent to the system Eqs. (5)−(10) of linear PDEs plus an algebraic equation. Therefore, linear PDE solvers can be used;

• (ii)

  Because we have a linear system, we can use Fourier transforms to solve it. Our PM code relies on 42 FFTs, where a standard Poisson solver needs only two. The use of FFTs makes it faster than the finite-element codes N-MODY and RayMOND (Londrillo & Nipoti 2009; Candlish et al. 2014), especially for large grids;

• (iii)

  Four iterations are required to obtain an accuracy of 1%;

• (iv)

  Self-gravity is the dominant cause of the numerical error. We needed a Gaussian smoothing kernel of one pixel width over a spherical domain with a radius of four pixels. The standard CIC method used in Poisson solvers cannot be used;

• (v)

  When the number of particles per pixel exceeds 1/10³, particle propagation takes more time than the PDE solver;

• (vi)

  We derived the analytic isothermal state Eqs. (26)−(29) for the spherical deep MOND case. It is a two-parameter family, with a mass and a length scale. Because the deflection potential deviates significantly from that for a point mass, the solution could be relevant for elliptic galaxies and clusters;

• (vii)

  We tested our PM code for systems where deep MOND analytic solutions (by Milgrom 1994) are known and found good agreement. We also tested the isothermal state by creating a Maxwell-Boltzmann distribution with density Eq. (29) and verified that the system is stationary.

Promising applications of the N-body MOND code include the numerical comparison with actual astrophysical observations on the Solar System, wide-binaries, and galaxy clusters.

Appendix A Table of symbols

Appendix B Fourier transforms

In the text, we did not distinguish much between continuous functions on ℝ³ and the discrete representation on (ℤ/nℤ)³ for the n × n × n grid. We simply considered the discrete functions as an approximation of the continuous functions on a finite interval. Let the physical dimension of the cubic pixels be s, and let the cubic volume have length L = ns. We assume that the functions smoothly drop to zero (or a constant) for $\left|x\right|>\frac{L}{2},\left|y\right|>\frac{L}{2},\left|z\right|>\frac{L}{2}$. The position coordinates that correspond to the grid points take on the values $$r=\left(\begin{array}{l}x\hfill \\ y\hfill \\ z\hfill \end{array}\right)=\left(\begin{array}{l}i\hfill \\ j\hfill \\ k\hfill \end{array}\right)s,i,j,k\in \left\{-\frac{n}{2},\dots ,\frac{n}{2}-1\right\}.$$

Fourier space has the elements k as reciprocal vectors. It has the grid values $$k=\left(\begin{array}{l}{k}_1\hfill \\ k_2\hfill \\ k_3\hfill \end{array}\right)=\left(\begin{array}{l}i\hfill \\ j\hfill \\ k\hfill \end{array}\right)\frac{2\pi }{L},i,j,k\in \left\{-\frac{n}{2},\dots ,\frac{n}{2}-1\right\}.$$

The scalar function ϕ(r) is expressed in its Fourier transform $ℱ(\varphi )=\widehat{\varphi }(k)$ via the expansion in plane waves, that is, the inverse Fourier transform $$\varphi (r)=\frac{1}{{(2\pi )}^3}{\displaystyle \iiint \widehat{\varphi }}(k){\text{e}}^{\text{i}k•r}\text{d}{k}_1 \text{d}{k}_2 \text{d}{k}_3\approx \frac{1}{L^3}{\displaystyle \sum _k\widehat{\varphi }}(k){\text{e}}^{\text{i}k•r},$$

and the Fourier coefficients, that is, the direct Fourier transform is $$\widehat{\varphi }(k)={\displaystyle \iiint \varphi }(r){\text{e}}^{-\text{i}k•r}\text{d}x \text{d}y \text{d}z\approx s^3{\displaystyle \sum _r\varphi }(r){\text{e}}^{-\text{i}k•r}.$$

In our PM algorithm, we modeled point particles with a smooth density centered around the points r_i. The standard CIC method uses the n-fold convolution of the box function box(x/s)box(y/s)box(z/s) for the nth-order method (Hockney & Eastwood 1988). Here, the box function is defined by $$\mathrm{box}(y)=\{\begin{array}{ll}1\hfill & \text{ if \hspace{1em}}\left|y\right|<\frac{1}{2},\hfill \\ 0\hfill & \text{ else }\hfill \end{array}.$$

For large n, the convolution in one dimension and its Fourier transform may be approximated by $$\begin{array}{l}\frac{1}{s^n}\underset{n\text{ times }}{\underbrace{\mathrm{box}\left(\frac{x}{s}\right)*\cdots *\mathrm{box}\left(\frac{x}{s}\right)}}\approx \frac{1}{s\sqrt{\pi n/6}}\exp \left(\frac{-6x^2}{n{s}^2}\right),\text{ \hspace{0.17em}\hspace{0.17em}for \hspace{0.17em}\hspace{0.17em}}n\gg 1\text{, }\\ {\left(\mathrm{sinc}\frac{s{k}_1}{2}\right)}^n\approx \exp \left(\frac{-n{s}^2{k}_1^2}{24}\right),\text{\hspace{1em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em} \hspace{0.17em}for \hspace{0.17em}\hspace{0.17em}}n\gg 1.\end{array}$$

We therefore used the following density, with order n = 12: $$\rho (r)=\frac{m_i}{s^3{(2\pi )}^{3/2}}\exp \left(\frac{-{\left|r-r_i\right|}^2}{2s^2}\right).$$

The Fourier transform of this single-particle density is $$\widehat{\rho }(k)=m_i\exp \left(\frac{-s^2|k{|}^2}{2}\right){\text{e}}^{-\text{i}k•r_i}.$$

A standard deviation of a single pixel in position space implies a drop off of below ${\text{e}}^{-{\pi }^2/2}=.007$ at the boundary of the reciprocal space. Furthermore, the density is very nearly spherically symmetric. The difference between the 12th-order CIC function and the Gaussian with the same variance in three dimensions is 10⁻⁷ at most.

Appendix C Flow scheme of the PM code

This is the initialization procedure.

• 1.

  Create a lookup table for the exponential function.

• 2.

  Create a table of the 251 coordinates in a sphere of radius 4.

• 3.

  Create a n³ matrix with the vector function k.

• 4.

  Create a n³ matrix for the inverse Laplace operator, $$\text{A}=\{\begin{array}{lll}\frac{1}{k•k}\hfill & \text{ if }\hfill & k\ne 0,\hfill \\ 0\hfill & \text{ if }\hfill & k=0.\hfill \end{array}$$

• 5.

  Clear a n³ matrix for the density: ρ_B(r) = 0.

• 6.

  Loop through the particles i = 1 … N.

  • Set the mass m_i, and initialize the position r_i and velocity v_i.

  • Next i. End particle loop.

This is the main loop for time steps in the simulation.

• 1.

  Loop through all particles i = 1 … N (in case n³ > 251 N).

  • Remove particle i from the mass density by filling the pixels with zeros according to the list of 251 pixels.

  • Next i. End particle loop.

• 2.

  Loop through all particles i = 1 … N.

  • Implement the position step of Leapfrog: r_i → r_i + v_i∆t.

  • Add particle i to the mass density, with the list of 251 pixels and the lookup table: ρ_B(r) → ρ_B(r) + ρ_i(r), using Eq. (14).

  • Next i. End particle loop.

• 3.

  FFT the density: ${\widehat{\rho }}_{\text{B}}(k)=ℱ\left({\rho }_{\text{B}}\right)$.

• 4.

  Calculate ${\widehat{\varphi }}_{\text{N}}(k)=4\pi G\text{A}{\widehat{\rho }}_{\text{B}}$.

• 5.

  Inverse FFT the potential: ${\varphi }_{\text{N}}(r)={ℱ}^{-1}\left({\widehat{\varphi }}_{\text{N}}\right)$.

• 6.

  Find Newton gravity g_N(r) from ϕ_N with finite differences.

• 7.

  Initialize the mass flux field by F(r) = g_N.

• 8.

  Start the iteration loop.

  • Calculate the gravity field: $g_{\text{M}}(r)=v\left(\sqrt{F•F}/a_0\right)F$.

  • FFT the gravity field: ${\widehat{g}}_{\text{M}}(k)=ℱ\left(g_{\text{M}}\right)$.

  • Calculate the potential in Fourier space: $\widehat{\varphi }(k)=\text{iA}k•{\widehat{g}}_{\text{M}}$.

  • Make the field conservative: ${\widehat{g}}_{\text{M}}(k)=-\text{i}k\widehat{\varphi }$.

  • Exit the iteration loop at the fourth iteration.

  • Inverse FFT the gravity field: $g_{\text{M}}(r)={ℱ}^{-1}\left({\widehat{g}}_{\text{M}}\right)$.

  • Calculate the field: $F(r)=\mu \left(\sqrt{g_{\text{M}}•g_{\text{M}}}/a_0\right)g_{\text{M}}$.

  • Calculate the magnetic component, H(r) = F − g_N.

  • FFT the magnetic field: $\widehat{H}(k)=ℱ(H)$.

  • Make the field divergence free: $\widehat{H}(k)\to \widehat{H}-k\text{A}k•\widehat{H}$.

  • Inverse FFT the magnetic field: $H(r)={ℱ}^{-1}(\widehat{H})$.

  • Calculate the field with F(r) = g_N + H.

  • Next iteration.

• 9.

  Inverse FFT the potential: $\varphi (r)={ℱ}^{-1}(\widehat{\varphi })$.

• 10.

  Loop through all particles i = 1 … N.

  • Evaluate g_i, with the list of 251 pixels and the lookup table using Eq. (15).

  • Implement the velocity step of Leapfrog: v_i→ v_i + g_i∆t.

  • Next i. End particle loop.

• 11.

  Update time t → t + ∆t, and

• 12.

  Go to the next time step.

References

All Tables

All Figures

Fig. 1 Anti-aliasing needed to suppress the self-gravity. A cross section of the Gaussian density Eq. (14) with one-pixel variance about (blue) points in the middle grid cell, representing a 12th order weight function. We only account for the pixels inside a sphere with a radius of 4s (dashed circles). This results in a domain of 251 cells (inside the contours). The linear motion of a point particle across the middle cell (left to top right) shows that the simulated density changes smoothly between panels.

In the text

Fig. 2 Simulated orbits of the two-body system projected onto the xy plane for the mass ratio m1/m2 = 3/2. The mesh has a size of 1283, and there are 100 time steps, each with four iterations for the MOND force. Red (m1) and cyan (m2) are the orbits obtained from numerical integration of Eqs. (19). Purple and blue are the orbits calculated with the PM code. The difference between the PM code and the ODEs is a slow dephasing of the orbits.

In the text

Fig. 3 Angular momentum vs. time t for the simulation shown in Fig. 2. These are constants of motion of the ODE system Eq. (19) and numerically constant (with machine precision, not shown) when integrated with Leapfrog. In the PM code, the total angular momentum (blue) is well conserved, but the single-particle angular momenta (orange and red) have large errors when the distance to the origin is smaller than one pixel.

In the text

	Fig. 4 Force times radius (blue) and squared velocity (orange) vs. mass ratio for the two-body system. The horizontal axis is logarithmic, making the curves symmetric about m1 /m2 = 1.
In the text

Fig. 5 Simulated orbits of two bodies with a mass ratio m1/m2 = 3/2 as in Fig. 2, but with 400 time steps. The initial velocities are given by Eq. (20) and should result in circular orbits. These are plotted in Fig 4. This simulation thus confirms the force formula Eq. (17) for the deep MOND case and validates the algorithm. It also illustrates how well the algorithm performs on the subpixel scale.

In the text

	Fig. 6 Force times radius (lower curves) and squared velocity (higher curves) vs. the mass ratio for particles on a regular N-gon orbiting a central body in circular motion. From Eqs. (2l)–(22), where m0 is the central mass and Nm is the ring mass.
In the text

Fig. 7 Simulated orbits for a ring of N = 100 particles with a mass m around a central mass m0 for a mass ratio of m0/Nm = 3/2 as in Fig. 2. The orbits of the first 80 time steps of 10 particles are shown in color. In gray we plot 576 steps for all particles. The initial velocities are given by Eq. (22), plotted in Fig. 6, which in theory results in circular orbits. Although this confirms the force formula Eq. (21), numerical round-off errors on the subpixel scale propagate and cause the orbits to diverge.

In the text

Fig. 8 Potentials and mass distributions for the isothermal state. In teal and orange we show the potential ϕ(r) as given by Eq. (26) and the deflection potential ψ=∫ϕdz for the line of sight passing the center at distance r, determining weak gravitational lensing. The asymptotes are shown (in gray) and the units are $\sqrt{GM{a}_0}{b}^j$, j = 0,1. In blue and red we show the density ρB(r) as given by Eq. (29) and the surface mass density σ=∫ρBdz in units of 3M/4πbj, j = 3,2.

In the text

	Fig. 9 Cumulative mass M(r) of the spherical isothermal galaxy vs. r. The (free) parameter b determines the galaxy size. The (blue) curve is the theoretical result given by Eq. (28). The realizations (orange, red, and purple) are calculated from Eq. (30).
In the text

	Fig. 10 Cumulative velocity distribution for the spherical isothermal galaxy. This is the fraction #(≤v)/N of particles with velocities below v vs. the value v of this upper bound. The realizations (orange, red, and purple) are calculated from Eq. (30).
In the text

Fig. 11 Surface-mass densities σ(x,y) and deflection potentials ψ(x,y) for a globular cluster from the central 128 × 128 part of a simulated cube of 2563 pixels. The yellow dots are the stars that define the actual density. The resulting apparent dark matter surface density is plotted in blue according to Eq. (11). Negative value are shown in red (Milgrom 1986). They are found for N ≤ 102. The equipotential curves of the deflection potential for Newtonian and Milgromian gravity are plotted in green and cyan. Although the densities are far from circularly symmetric, the deflection potential is nearly circular for N ≥ 103. This indicates that the retrieval of a clumpy mass distribution from weak lensing may be problematic for large N.

In the text

Fig. 12 Energies vs. time for N = 100 particles on a 2563 grid. The kinetic, potential, and total energy for the self-consistent interacting system of the particles from the PM code are shown as solid curves (red, purple, and blue). The kinetic, potential, and total energy from the ODE system of the particles in the external potential are shown as dashed curves (orange, violet, and teal). The fluctuations in the total energies are purely numerical, and the other energies have thermal fluctuations. The dashed gray curve is Epot = ∫ ρBϕdxdydz with the PM matter density and Eq. (26), hence substituting a finite N solution in the potential obtained in the thermodynamic limit. We subtracted the self-energies.

In the text

Fig. 13 Runtime vs. particle count N for a 1283 grid with the Python code on an Intel Core i7-9750H. The time spend on the FFTs (blue) is roughly constant, and the time for the evaluation of the particle accelerations and/or densities grows linear with N (orange, red, and purple). The total runtime shows a crossover around 103 particles (black). The dashed lines are guides with slopes of zero and one.

In the text