Fast particle-mesh code for Milgromian dynamics

Modified Newtonian dynamics (MOND) is a promising alternative to dark matter. To further test the theory, there is a need for fluid- and particle-dynamics simulations. The force in MOND is not a direct particle-particle interaction, but derives from a potential for which a nonlinear partial differential equation (PDE) needs to be solved. Normally, this makes the problem of simulating dynamical evolution computationally expensive. We intend to develop a fast particle-mesh (PM) code for MOND (the AQUAL formalism). We transformed the nonlinear equation for MOND into a system of linear PDEs plus one algebraic equation. An iterative scheme with the fast Fourier transform (FFT) produces successively better numerical approximations. The algorithm was tested for dynamical systems in MOND where analytical solutions are known: the two-body problem, a body with a circular ring, and a spherical distribution of particles in thermal equilibrium in the self-consistent potential. The PM code can accurately calculate the forces at subpixel scale and reproduces the analytical solutions. Four iterations are required for the potential, but when the spatial steps are small compared to the kernel width, one iteration is suffices. The use of a smoothing kernel for the accelerations is inevitable in order to eliminate the self-gravity of the point particles. Our PDE solver is $15$ to $42$ times as slow as a standard Poisson solver. However, the smoothing and particle propagation takes up most of the time above one particle per $10^3$ pixels. The FFTs, the smoothing, and the propagation part in the code can all be parallelized.

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 0 ≈ .127• c 2 Λ 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 0 .Because MOND is nonrelativistic, it applies in the regime in which |φ| c 2 and |∇φ| c 2 Λ 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).
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(Iorio , 2017)).

Theory and main result: The idea of the method 114
As the particles move in space under the influence of the grav-115 itational potential φ(r, t), they experience an acceleration given 116 by Newton's second law, (1) The baryonic mass density for the set of N point masses is Here, δ is the three-dimensional Dirac-delta function.Whereas 119 the potential is a solution of the Poisson equation in the Newton 120 theory, in MOND (the AQUAL formalism), the potential is a 121 solution to the following nonlinear partial differential equation 122 (PDE): with the boundary condition ∇φ(r, t) −→ 0 as |r| −→ ∞.
2.1.An equivalent system for AQUAL 135 We now demonstrate how to rewrite the Eq.
(3) for the potential 136 in MOND by transforming this equation into a system of coupled 137 linear PDEs for four vector fields g N , g M , H, and F, combined 138 with one algebraic equation.This system of equations is (10) The vector fields g N and g M are the Newton and MOND accel-141 eration fields.With the identification g N = −∇φ N , the two equa-142 tions in Eq. ( 5) are equivalent to the Poisson equation.Our new 143 system effectively expresses the MOND gravity field g M in the 144 Newtonian gravity field g N .Equations ( 8) and ( 10) are not inde-145 pendent equations, but are the inverse of each other.Here, F and 146 g M are the lengths of the vectors F and g M .The pair µ(x) and 147 ν(y) relate the strength of the Newtonian gravity to the strength 148 of the Milgromian gravity by 149 µ(x)x = y = F a 0 , and ν(y)y = x = g M a 0 .
The choice Eq. ( 4) for the interpolation function gives 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 g M .We thus need to identify g M = −∇φ.When we now combine Eq. ( 6) with (10), solve for g N , then substitute this solution in the first of Eq. ( 5), and use Eq. ( 7), we obtain Equation ( 9) states that the field is curl-free.Hence, g M is a conservative field with a potential.When we substitute g 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 g 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 g N and H, which must satisfy the second of Eq. ( 5) and Eq. ( 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 g 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).
3. Project g M −→ g M onto the curl-free part, as per (9).
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 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.
lensing) in the incorrect Newton-Poisson equation.This is equiv-207 alent to The apparent density ρ A is also known as the effective density 209 (see Banik & Zhao 2022).This leads to an apparent dark matter 210 density of ρ A − ρ B , also called phantom dark matter density.The 211 total baryonic mass and apparent mass in a volume W can be 212 found from the flux through its boundary, 2.2.An equivalent system for QuMOND

219
We do not consider QuMOND any further.However, the algo-220 rithm can easily be adapted using the same steps.We now describe the implementation of the method.The code 223 needs to numerically approximate the solution of the MOND 224 PDEs in Eqs. ( 5)-( 10) and propagate particle motion and a time-225 dependent mass density.The particle equations of motion Eq. ( 1) 226 are implemented with Leapfrog, Although the formulas do not explicitly show this, the velocities 228 and positions represent the values at a half time step difference.229 In principle, an actual half time step needs to be made initially 230 (and also finally) to obtain the true Leapfrog integration scheme.231 The code uses a discrete Fourier transform, so that space is 232 represented as a periodic cubic lattice of n 3 cells we call pixels.233 Hence, the scheme does not have nested grids; the grid cells are 234 all equal and do not change.However, due to the periodic na-235 ture of the simulated space, matter inside one cubic volume is 236 effectively influenced by the gravity of the neighboring copies.237 This effect is nonphysical, and a volume of empty space needs 238 to surround the matter when this effect is to be small.Ideally, 239 the center of mass is placed in the origin, and it is ensured that 240 particles stay away from the boundary of the cube.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.

Interpolation for the particle density function 242
The main drawback of the particle-mesh (PM) method is the lack 243 of accuracy at the short-distance scale of the pixels.Although  we also need to approximate the acceleration on particle i by the 250 field at r i , for which the value at the center of the pixel is not 251 accurate.

252
We used a Gaussian kernel to create a smoothed density in 253 order to suppresses anti-aliasing.The density of Eq. ( 2) was re-254 10 −2 10 −1 10 2 0 .5 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 m 1 /m 2 = 1.
placed by 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.

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 g 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: 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 .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.

Results: Validation using analytic solutions 284
The problem of self-gravity dominating the actual force on a

302
For this regime, Milgrom (1994Milgrom ( , 1997Milgrom ( , 2014) ) found, using the 303 method of the virial, that the kinetic energy for a bound system 304 with the center-of-mass motion at rest is with m i the individual particle masses, and M the total mass.If  We considered the two-body problem in MOND.This is a model 312 for a binary stellar system with two comparable masses with a 313 wide separation.The interparticle force is found from Eq. ( 16) 314 to be with M = m 1 + m 2 , and where r = |r 2 − r 1 | is the separation.We can rewrite the mass dependence using The equations of motion of the two particle system are 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: 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 1 and L 2 are also constant.However, the numerical error in L 1 and L 2 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 1 and r 2 , the errors disappear.For grid sizes 128 3 , 256 3 , and 512 3 , 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 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 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. ( 21)-( 22), where m 0 is the central mass and Nm is the ring mass.
and for the kinetic energy Using Eqs. ( 16)-( 18), we find that the velocity v of the reduced motion and the orbital frequency ω are given by the expression 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.

Example II: Ring galaxy
Again, we tested the code with a system that was solved analytically by Milgrom (1994): a central mass m 0 , surrounded by a ring with a radius r 2 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), The force on the particles in the ring is given by Here, M = m 0 + Nm.We can now write x/s Fig. 7. Simulated orbits for a ring of N = 100 particles with a mass m around a central mass m 0 for a mass ratio of m 0 /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.
The velocity of the ring particles and the orbital frequency are 380 found by equating the force Eq. ( 21) to The force times radius and squared velocity is plotted in Fig. 6. 382 The result of the simulation for N = 100 particles is plotted in 383 Fig. 7.If the mass of the ring is not very low compared to the 384 central mass, the system is unstable, and the ring breaks apart 385 after a few orbits.Our third application is the simulation of a spherically symmetric 388 system in thermodynamic equilibrium.This could be a globular 389 cluster with up to 10 6 stars or a galaxy cluster with up to 10 3 390 galaxies.These system are in the deep MOND regime or in the 391 crossover regime (accelerations around a 0 ).

392
A system of N point particles with equal mass m in hydro-393 static equilibrium with a temperature T has a velocity dispersion 394 given by The last equality is the thermodynamic limit of Eq. ( 16).For 396 a spherical system, the particle number density in terms of the 397 potential is given by Here, we used hydrostatic equilibrium and the ideal gas law.399 Substitution of the mass density ρ B (r) = mn(r) and Eq. ( 24) in 400 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 √ GMa 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 σ = ρ B dz in units of 3M/4πb j , j = 3, 2. .8 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).
the MOND Equation (3) gives the following equation for the po- The solution of this equation is found to be 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).
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 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), The solution given in Eqs. ( 26)-( 29) is applicable for g M a 0 .Because the maximum acceleration is found at r = b/2 2/3 , the requirement becomes GM/b 2 9a 0 /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.give x/s 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 256 3 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 ≤ 10 2 .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 ≥ 10 3 .This indicates that the retrieval of a clumpy mass distribution from weak lensing may be problematic for large N.
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 In order to obtain random realizations, our point masses were 447 found from seven random numbers ξ 1 ∈ (−1, 1), η 1 , η 2 , η 3 ∈ 448 (0, 1), ζ 1 , ζ 2 , ζ 3 ∈ (0, 2π), sampled homogeneously from these 449 intervals and the formulae 450 Fig. 12. Energies vs. time for N = 100 particles on a 256 3 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 E pot = ρ 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.
11 shows the realizations projected onto the celestial plane and 455 the resulting apparent dark matter surface density.For N = 10 3 456 or larger, the potential is nearly spherical, and it is well described 457 with our analytic solution.

458
We solved the evolution of N particles in the external po-459 tential Eq. ( 26) as a system of ODEs, and we simulated the N 460 particles in MONDian gravity using the PM code, both with the 461 same initial distribution.The two models are comparable be-462 cause both must behave identical in the thermodynamic limit 463 N −→ ∞.We also tested whether the distributions are stationary, 464 as they should be for large N.In the external potential, the ener-465 gies and angular momenta of the individual particles are strictly 466 conserved, and so is the total energy and angular momentum.

467
The single-particle kinetic energies and momenta oscillate, and 468 so does the total kinetic energy and the total momentum.We creates a minimum size, we reach a subpixel spatial resolution 489 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 finitedifference gradient of this scalar function.
The runtime of our code has a time complexity of O(Nn 3 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 3 grid, the evaluation of the accelerations or/and densities starts to dominate the runtime at N = 10 3 .Therefore, in the general case, we expect that the crossover occurs at a critical density of 10 3 /128 3 ≈ 1/10 3 particles per pixel.

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 −3 , 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 threedimensional FFTs need to be performed sequentially.For each direction, we have n 2 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 ∆ρ B , ∆g N , ∆g M , ∆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

(
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 2 and 10 3 , 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).

221 3 .
The algorithm: Dealing with the self-gravity 222

Fig. 2 .
Fig.2.Simulated orbits of the two-body system projected onto the xy plane for the mass ratio m 1 /m 2 = 3/2.The mesh has a size of 128 3 , and there are 100 time steps, each with four iterations for the MOND force.Red (m 1 ) and cyan (m 2 ) 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.
244each particle is located in a single pixel, we desire accuracy that 245 is below the grid size.This influences both the field created by 246 and the acceleration on the particle.First, the field created by 247 particle i close to the mass should resemble a spherical field 248 about the point r i , which is not the center of a pixel.Second, 249

Fig. 5 .
Fig. 5. Simulated orbits of two bodies with a mass ratio m 1 /m 2 = 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.
285point particle is inherent to the PM method for MOND.The 286 effect of an imprecise subpixel resolution for the position of a 287 particle and the consequently imprecise evaluation of the accel-288 eration on that particle is always present, even for a single par-289 ticle.We therefore also need to test the accuracy of the method 290 for linear motion.We have found that in our simulations, a point 291 particle of mass m experiences nonphysical accelerations with 292 the average value of 1.08 • 10 −3 √ Gma 0 /s for a 256 3 grid (see 293 de Nijs 2023b).This small error is independent of the number of 294 pixels.295 In this section, we present results of simulations with the 296 code for three systems: the two-body problem, a ringed system, 297 and an isotropic three-dimensional system.Because the main in-298 terest is in MOND, we considered the deep MOND case.This 299 means that the interpolation functions we used in our simulation 300 are 301 µ(x) = x, ν(y) = 1 √ y .

306
there is no external field, this single equation allows us to cal-307 culate the MOND forces between particles for highly symmetric 308 cases, thereby bypassing the nonlinear PDE.We used these ana-309 lytic solutions to test our PM code.
310 4.1.Example I: Wide binary system 311 Eample III: The isothermal sphere in MOND 387 η 3 k cos ζ 3 .The Box-Muller transform was used to generate the Gaussian 451 values of the Maxwell-Boltzmann velocity distribution.Figures 452 9 and 10 show how close the realizations are to the cumulative 453 distributions for the radial coordinate and for the velocity.
469compare this with the exact system of Eqs.(1)-(3).Here, the total 470 energy, total momentum, and total angular momentum are con-471 served.The PM code simulates this system.The fluctuations in 472 the energies are shown in Fig. 12.Because the mesh in our PM 473 approach breaks spherical symmetry, angular momentum also 474 the algorithm arise from the fact that we used 478 a fixed grid.The discrete grid has n pixels in each dimension 479 with a pixel width s and has a finite volume n 3 s 3 .480 The use of FFTs in the calculation of the potentials and accel-481 eration fields generally gives rise to anti-aliasing: The functions 482 obtain nonphysical short-wavelength oscillations.This effect is 483 particularly strong when the point particles are represented with 484 isolated nonzero pixel values in the mass density, but also when 485 the CIC method is used.Smoothing with a minimum variance 486 of one pixel width s in a spherical domain of radius 4s needs to 487 be used to overcome this problem.Although the discretization 488 See Table A.1 in Appendix A for the sym-131 bols and notation.A popular choice for the interpolation func-132 tion that is consistent with all observations so far (Brouwer et al. 133 2021) is The 124 scalar function µ(x) is called the interpolation function because 125 it interpolates between the regime of high and low accelera-126 tions.At high acceleration, µ(x) −→ 1 so that the potential will 127 approach the Newton potential.At low acceleration, the deep 128 MOND regime begins, and µ(x) −→ x.An isolated and local-129 ized mass m has the asymptotic spherically symmetric potential 130 φ(r) = √ Gma 0 log |r|.
Table A.1.List of symbols and notation.Notes.Symbol and significance of the physical quantities used.Vectors and vector fields are in boldface.The length of a vector or the strength of a vector field is in normal font.Most MOND research papers simply use g for the Milgromian gravity field g M .