The Radial Acceleration Relation and a Magnetostatic Analogy in Quasilinear MOND

Recently a remarkable relation has been demonstrated between the observed radial acceleration in disk galaxies and the acceleration predicted on the basis of baryonic matter alone. Here we study this relation within the framework of the modified gravity model MOND. The field equations of MOND automatically imply the radial acceleration relation for spherically symmetric galaxies, but for disk galaxies deviations from the relation are expected. Here we investigate whether these deviations are of sufficient magnitude to bring MOND into conflict with the observed relation. In the quasilinear formulation of MOND, to calculate the gravitational field of a given distribution of matter, an intermediate step is to calculate the"pristine field", which is a simple nonlinear function of the Newtonian field corresponding to the same distribution of matter. Hence, to the extent that the quasilinear gravitational field is approximately equal to the pristine field, the radial acceleration relation will be satisfied. We show that the difference between the quasilinear and pristine fields obeys the equations of magnetostatics, the curl of the pristine field serves as the source for the difference in the two fields, much as currents serve as sources for the magnetic field. Using the magnetostatic analogy we numerically study the difference between the pristine and quasilinear fields for simple model galaxies with a Gaussian profile. Our principal finding is that the difference between the fields is small compared to the observational uncertainties and that quasilinear MOND is therefore compatible with the observed radial acceleration relation.


I. INTRODUCTION
The problem of dark matter is one of the central problems in modern cosmology [1]. According to the conventional ΛCDM paradigm, matter constitutes less than a third of the total energy density of the Universe. The bulk of this matter is non-baryonic. Many lines of evidence point to the existence of dark matter, notably the cosmic microwave background anisotropies [2] and observations of lensing [3], but dark matter has eluded direct detection [4] and its nature remains a mystery.
Rotation curves of disk galaxies provided the earliest compelling evidence for dark matter [5]. The observation that the rotation velocity generically approaches a plateau instead of vanishing asymptotically is conventionally interpreted as evidence for a spherical dark matter halo that dominates the visible baryonic matter in the disk. Puzzling from this perspective were empirical observations such as the baryonic Tully-Fisher relation [6] that appear to suggest a strong connection between the rotation curves and the distribution of baryonic matter.
In an important recent study McGaugh et al. [7,8] demonstrated a relation between the observed radial acceleration in disk galaxies and the acceleration predicted on the basis of baryonic matter alone. This empirical relation, dubbed the "radial acceleration relation" (RAR) is based on an impressive body of data, the SPARC dataset that represents near infrared observations by the Spitzer Space Telescope and 21 cm observations by radio astronomers over several decades [9]. The radial acceleration relation has to be explained by any tenable model of dark matter. There have already been several efforts to show that this relation is compatible with ΛCDM cosmology [10][11][12]. There have also been efforts to understand this relationship using modified theories of gravity ranging from scalar tensor theories [13] to emergent gravity models [14,15].
The modified gravity theory MOND was originally introduced [16] to account for galaxy rotation and it has had some phenomenological success, e.g., providing an explanation for the baryonic Tully-Fisher relation. However MOND and its relativistic completions have wellknown problems on cluster and cosmological scales [17][18][19]. It is sometimes assumed that MOND predicts the observed radial acceleration relation, but in fact the MOND field equations imply such a relationship only for spherically-symmetric galaxies [20]. For disk galaxies deviations are to be expected. With the empirical discovery of the remarkable radial acceleration relation it becomes imperative to determine the magnitude of these deviations; that is the object of this paper. Our principle finding is that the deviations are small compared to the experimental uncertainties and that the MOND field equations are in fact compatible with the observed radial acceleration relation.
Our analysis is carried out in the framework of 'quasilinear' MOND [21]. Notwithstanding the name, quasilinear MOND is a highly nonlinear theory. In order to calculate the quasilinear gravitational field, g Q , corresponding to a distribution of matter ρ D , an intermediate step is to calculate a "pristine field", g P . The pristine field is related to the Newtonian field through a simple, nonlinear 'interpolating function', thus, if one was to plot observed versus Newtonian accelerations within pristine MOND, finding agreement with the radial acceleration relation would simply be a matter of choosing the right interpolating function [8]. To the extent that the quasilinear field is approximately equal to the pristine field, it will obey the radial acceleration relation. In order to demonstrate the approximate equality of the two fields we show that the difference b = g P − g Q obeys the equations of magnetostatics with b analogous to the magnetic field and ∇ × g P analogous to the current that sources b. This analogy provides useful qualitative insights into the deviation field b as well as a practical numerical scheme for its calculation.
Although the main focus of this paper is on modified gravity formulations of MOND it is worth noting that there is another class of models called modified inertia formulations of MOND. Although modified inertia models have received comparatively less attention [22], there has been an attempt to construct a relativistic formulation [23] and an interesting experimental test has been proposed [24]. It would be desirable in future work to explore possible contrasts between the predictions of modified inertia and modified gravity theories of MOND for galactic motion along the lines suggested by [22].
For simplicity we study simple model galaxies with a Gaussian profile. Our methods as well as our results will easily carry over to other more realistic galaxy models. The numerical analysis is most efficiently done by use of the fast Fourier transform although there are some subtleties that arise due to the long range nature of gravitational forces. Key results are presented in fig 2 which shows that the current distribution ∇ × g P for a disk galaxy resembles a Helmholtz gradient coil; fig 4 in which the dependence of the deviation b on distance from the galactic center and on the aspect ratio of the disk galaxy is systematically investigated; and fig 5 in which we display the computed radial acceleration relation for our model galaxies. Comparison of this plot with the observed radial acceleration relation leads to our conclusion.

A. Quasilinear MOND
In Newtonian gravity, the acceleration due to gravity g N obeys the field equations Here ρ(r) is the mass distribution producing the gravitational field. In the pristine formulation of MOND [16] a test particle however experiences an acceleration g P which is related to the Newtonian acceleration via the algebraic relation Here a 0 ∼ 10 −10 m/s 2 is the MOND acceleration scale and the interpolating function f has the asymptotic behavior f (x) ≈ 1 for x 1 and f (x) ≈ x −1/2 for x 1. Thus the observed acceleration of the test particle matches the Newtonian prediction when the Newtonian gravitational field is strong but is significantly enhanced when the Newtonian field is weak on the MOND scale. Pristine MOND exactly reproduces the radial acceleration relation; indeed the observed radial acceleration relation may be regarded as a measurement of the MOND acceleration a 0 and the interpolating function f . However, pristine MOND suffers from physical inconsistencies and lacks a Lagrangian formulation. Thus it cannot be considered a viable alternative theory of gravitation as noted in [21] and [25]. These authors introduced nonlinear MOND [25] and quasilinear MOND [21] as viable alternatives to Newtonian gravity; in this paper we focus on quasilinear MOND. Within quasilinear MOND the acceleration of the test particle is not equal to the pristine field g P but to the quasilinear field g Q which is determined by solving ∇ · g Q = ∇ · g P ≡ −4πGρ eq and ∇ × g Q = 0. (3) Invoking the Helmholtz decomposition of vector fields we see that g Q is the curl free part of g P . The divergence of g P (or more precisely ρ eq ) has the interpretation of being the combined density of dark and baryonic matter that would be required to produce the same observed acceleration of the test mass if Newtonian gravity prevailed.
To complete the formulation of quasilinear MOND we must specify the interpolating function f . A simple form that has the desired asymptotic behavior is An alternative form, following ref [7], is to take

B. Magnetostatic Analogy
For a spherically symmetric galaxy g P is curl free and hence g Q = g P . In this case it automatically follows that the observed acceleration g Q will be a simple function of the Newtonian gravitational field produced by the baryonic matter in the galaxy, namely, g Q = g N f (g N /a 0 ). However only in cases of exceptional symmetry will g P be curl free. Generically g Q is not equal to g P and hence g Q is not a local function of g N . For the radial acceleration relation to hold in quasilinear MOND we need g Q ≈ g P at least in the galactic plane.
Whether the approximate equality g Q ≈ g P holds for disc galaxies is the key question that is studied in this paper. In [20] the corresponding problem was pointed out for nonlinear MOND and an estimate was made of the size of the discrepancy between the observed acceleration and the pristine field strength in that theory. The new high quality data establishing the radial acceleration relation provides the impetus for the detailed study of quasilinear MOND reported in this paper.
To this end it is convenient to define and It then follows from eq (3) that These are precisely the equations of magnetostatics. We have already noted that ∇ · g P may be interpreted as the equivalent density of dark matter needed to reproduce MOND behavior within Newtonian gravity. Eq (8) reveals that ∇ × g P may be interpreted as a current that sources the difference between g P and g Q . Given g P eq (8) provides an efficient way to calculate b and thereby g Q as we will discuss below. Formally one can write down a solution to eq (8) expressing b in terms of j using the Biot-Savart law. Whereas the relation eq (2) is local in the sense that g P at a given point is determined by g N at the same point, in sharp contrast the relation between b and j is not local. This shows that in general g Q is not a local function of g N in quasilinear MOND.

A. Model and Methods
We now introduce a simple model which allows us to study the deviation from the radial acceleration relation in disk galaxies. We take the galaxies to have a Gaussian profile Here M is the mass of the galaxy and R d is its radius. ∆ is a dimensionless parameter that characterizes the aspect ratio of the galaxy; ∆ = 1 corresponds to a spherically symmetric galaxy and ∆ → 0 to a flat disk. The characteristic radial acceleration scale within Newtonian gravity for this circumstance is GM/R 2 d . The MOND field of the galaxy is therefore fully determined by the dimensionless parameters µ = GM/R 2 d a 0 and ∆. We will work in a system of units wherein R d = 1 and a 0 = 1. Then the parameter µ = GM corresponds to a measure of the mass of the galaxy in our system of units. Other simple models used to describe disk galaxies [26] can be treated using our methods and will lead to similar results; the advantage of using the Gaussian model is that it allows us to do higher resolution numerics and to rapidly explore a broader range of MOND parameters. Hence we focus on the Gaussian model here.
The first step in our analysis is to calculate g N , the Newtonian field of the galaxy. Note that it is sufficient to calculate g N for µ = 1 since by linearity g N ∝ µ. Eq (1) has the simple solutioñ in Fourier space. (A tilde above a symbol denotes the Fourier transform.) The Newtonian field can therefore be efficiently calculated using the fast Fourier transform. We assume that the galaxy is located at the center of a cube of side L and we sample functions at points on a regular cubic lattice with 2N + 1 points per edge of the cube. An amusing subtlety arises because of the long-range nature of the gravitational force. By using the Fourier transform we enforce a periodicity in the mass distribution and the gravitational field whereas in the real problem we wish to impose the boundary condition that the gravitational field vanishes far away from the galaxy. Naively one might suppose that in the limit that L is sufficiently large the results would be insensitive to boundary conditions but it is easy to demonstrate that in fact a sensitivity to boundary conditions persists in the limit L → ∞ because of the inverse square falloff of the gravitational field [27]. The remedy is to calculate not the gravitational field of the disk galaxy directly but rather of the mass distribution ρ D − ρ sph . Here ρ sph is a spherically symmetric distribution with a total mass, M , the same as the galaxy distribution ρ D . The advantage is that the gravitational field of the difference falls off like a quadrupole as 1/r 4 . This is a sufficiently rapid falloff that the result is insensitive to boundary conditions for sufficiently large L. We can then construct the gravitational field corresponding to the galaxy distribution ρ G by adding the field corresponding to the spherical distribution ρ sph which can be obtained analytically. In practice we take ρ sph to be an isotropic Gaussian with radius R d ∆ 1/3 .
Calculation of g P via eq (3) is now a simple matter of applying an algebraic function to g N . b can be calculated by solving eq (8) by Fourier methods. In order to calculate the curl of g P it is helpful to first subtract a suitable isotropic function to ensure that the results are insensitive to boundary conditions much as we do in our computation of g N . In practice we subtract the pristine field corresponding to the isotropic distribution ρ sph to facilitate the computations. Once g P and b are determined it is a simple matter to obtain g Q using eq (6).
In the results presented here we consider aspect ratios in the range 0.5 < ∆ < 1.0 and the mass parameter in the range 0.01 < µ < 10.0. Over this range of parameters it is sufficient to calculate the Newtonian field using L = 4 and N = 36. Smaller values of ∆ are more realistic for disk galaxies but would necessitate a finer grid. However in order to show the qualitative difference between pristine and quasilinear MOND, and to show that quasilinear MOND reproduces the observed radial acceleration relation it is sufficient to consider ∆ = 0.5 as discussed further in section III B; we have verified that the results extrapolate smoothly with ∆ down to ∆ = 0.1. The results are also insensitive to variations the size of the box and to the fineness of the discretization. We have also checked the accuracy of the calculated Newtonian field by using the multipole expansion to calculate the large distance asymptotic Newtonian field corresponding to the density profile eq (9). To check the accuracy of the calculated Newtonian field at short distances we have performed a non-trivial check described in appendix A. The subsequent calculation of g P is limited only by numerical roundoff since it only involves solving an algebraic relation. The calculation of g Q from g P is similar to the calculation of g N and can be checked in similar ways.

B. Results
We start with a discussion of the quantity ρ eq introduced in eq (3). If we cast a dark matter interpretation upon MOND we may regard ρ D as the density of baryonic matter and ρ eq as the combined density of baryonic and dark matter. Fig 1 shows the two densities for an isotropic galaxy with aspect ratio ∆ = 1 and mass parameter µ = 0.1. The plot of ρ D in the left panel corresponds to the Gaussian profile eq (9). The plot of ρ eq in the right panel shows three interesting features. First the density ρ eq is much greater than ρ D , consistent with the preponderance of dark matter compared to baryonic matter in the conventional ΛCDM paradigm. Second the density extends to a greater distance from the galactic center much like a dark matter halo in the conventional ΛCDM picture. Third the density has a weak divergence near the center of the galaxy (ρ eq ∝ 1/ √ r for r → 0); this is a generic feature of quasilinear MOND and it should be noted that it does not imply that the core of the galaxy has an infinite mass since the divergence is integrable. The plots in fig 1 were generated using analytic expressions for g N and g P which are available in the isotropic case. Our numerical methods can be employed to examine the behavior of ρ eq for disk galaxies. Those results are of interest but they will be reported elsewhere [28] since they are tangential to the present study.
According to the magnetostatic analogy, the quantity j = ∇ × g P may be regarded as a current that serves as the source for b, the difference between the pristine and quasilinear MOND fields. Using the product rule it is easy to see that ∇ × g P = f (g N /a 0 ) g N × ∇(g N /a 0 ). Due to the azimuthal symmetry of the galaxy it is easy to see that neither g N nor ∇g N can have an azimuthal component; it follows that j must be purely azimuthal. Thus symmetry requires that the current that sources b circulates about the axis of the galaxy. Similarly the assumed reflection symmetry of the galaxy about the x-y plane implies that the current j must be antisymmetric under reflection about the x-y plane. Thus the currents circulate in opposite senses above and below the galactic plane. This is as far as symmetry arguments go. Fig 2  shows an explicit numerical calculation of the y component of the current as a function of position in the x-y plane for a galaxy with aspect ratio ∆ = 0.5 and mass parameter µ = 0.1. The plot shows two sharp peaks and two sharp anti-peaks that are related by the azimuthal and reflection symmetries noted above. It follows that the current is concentrated into two counter-propagating circular coils that are parallel to the x-y plane, one located above the plane and the other below it, analogous to a Helmholtz gradient coil. We can easily draw upon electromagnetic intuition to draw the lines of b or to recognize that at large distance the resulting field will be quadrupolar.
We turn now to a comparison of the Newtonian, pristine MOND and quasilinear MOND gravitational fields. The upper panels of Fig.3 plot the radial component of the field as a function of the cylindrical coordinate r. Each field is computed at 36 points in the galactic plane at distances that vary uniformly from zero to two galac-tic radii. Both plots are for galaxies with aspect ratio ∆ = 0.5. The plot at left is for a galaxy that is deep in the MOND regime with parameter µ = 0.1; the plot on the right is deep in the Newtonian regime with parameter µ = 10. As expected the Newtonian field shown in blue vanishes both as r → 0 and r → ∞. We can check the large r behavior using a multipole expansion and we can also perform a non-trivial check at small r as described in Appendix A. The orange curves correspond to the pristine MOND acceleration. This is simply related to the Newtonian field via the algebraic function eq (2). The plots show that for small µ the pristine field is considerably enhanced relative to the Newtonian field but for large µ it is essentially equal to the Newtonian field. In fig 3 we have used the interpolating function in eq (4) but we have explicitly verified that the alternative interpolation eq (5) gives similar results. The black dashed curves show the radial component of the quasilinear acceleration. For this aspect ratio it appears that there is very little difference between the pristine and quasilinear fields at least in the galactic plane. We will examine this difference more closely below. The lower two panels of Fig.3 show the z component of the Newtonian, pristine, and quasilinear fields as a function of z. A key difference from the radial plots is that in this case the magnitude of the pristine acceleration is bigger than the magnitude of the quasilinear acceleration. This difference can be understood by picturing the difference between the two fields, b, as the magnetic field of a Helmholtz gradient coil. Another noteworthy feature is that the difference between the pristine and quasilinear fields is bigger along the axis than in the galactic plane, though in neither case is the difference especially large. The plots are for a galaxy with mass parameter µ = 0.1. Because the difference of the pristine and quasilinear fields is a small quantity the finiteness of the grid is visible in the jagged appearance of the plots; the deviations from a smooth fit are an estimate of the precision of the numerics. Another feature lost to low resolution is that close to the origin the radial component of b turns down and approaches zero. We have verified that with higher resolution grids this behavior is seen in the numerics. Furthermore except very near to the origin the high resolution numerics agree with the coarser plot shown in fig 4. As one might expect the difference between the pristine and quasilinear fields vanishes in the isotropic limit ∆ = 1. As the disk becomes flatter with decreasing ∆ the difference between the fields grows. The effect for ∆ = 0.1 (the approximate aspect ratio of the Milky Way, for example) is not very different from that of ∆ = 0.5. Since the latter aspect ratio can be simulated accurately with a coarser grid this is the canonical aspect ratio we have used in many of our plots. For the purposes of this paper a more faithful model of the galaxies is not necessary. It is also interesting to study how the magnitude of b varies with mass parameter µ for a fixed aspect ratio ∆. We find that the magnitude is bigger for a galaxy deep in the MOND regime than for a Newtonian galaxy; in other words it grows as µ decreases.
In fig 5 we plot the radial component of the pristine and quasilinear MOND acceleration as a function of the corresponding Newtonian acceleration for our model galaxies. In comparing this plot to the radial acceleration relation of ref [7] we note that the Newtonian acceleration corresponds to the computed baryonic acceleration in ref [7] while the pristine and quasilinear accelerations correspond to the observed acceleration of ref [7] within pristine and quasilinear MOND theories respectively. For pristine MOND by definition the observed acceleration is a simple algebraic function of the baryonic acceleration (grey curve in fig 5). For quasilinear MOND for our model galaxies there should be a distinct curve for each value of the galaxy parameters ∆ and µ. This is indeed borne out by fig 5 in which we have plotted the quasilinear predictions only for ∆ = 0.5 but for seven different values of µ. However because of the small difference between the pristine and quasilinear accelerations, all of these curves lie close to the pristine curve, and hence we arrive at our principal conclusion that quasilinear MOND is compatible with the observed radial acceleration relation. The deviations from the smooth pristine behavior are bigger for ∆ = 0.1 [29] but even in this case they remain smaller than the observational uncertainties described in ref [7]. Whether there are systematic trends buried in the observed deviations when the data are grouped according to parameters like µ and ∆, trends that can be predicted by quasilinear MOND, is a difficult problem we leave open for future work.
Although not directly related to the radial acceleration relation, we now consider oscillatory motion perpendicular to the galactic plane. A test particle moving in a circular orbit in the galactic plane will undergo such oscillations if perturbed suitably. If this oscillation frequency could be measured it would provide a clear distinction between MOND and ΛCDM. To calculate this frequency we note that in the galactic plane, by symmetry the z component of the acceleration due to gravity must vanish. Hence near the galactic plane, the acceleration due to gravity varies linearly with z and hence a particle perturbed away from the galactic plane executes simple harmonic motion. The frequency of oscillations is given by ω 2 = −∂g z /∂z. Here g z denotes the z component of g Q within quasilinear MOND and of g P in the pristine approximation to it. In the ΛCDM paradigm g z denotes the z component of the Newtonian field of visible matter and the dark matter halo combined. In fig 6 we show the oscillation frequency computed within quasilinear MOND, the pristine approximation and due to the Newtonian field of visible matter for a disk galaxy with aspect ratio ∆ = 0.5 and mass parameter µ = 1. The Newtonian field of visible matter dominates the ΛCDM contribution because it constitutes a disk whereas the dark matter is dispersed into a spherical halo. Taking the visible Newtonian frequency to be an estimate of the Note also that the deviation between gP and gQ is more pronounced in the axial plots than the radial plots; the sign of the deviation is also opposite in the two cases.
total ΛCDM oscillation frequency we see that MOND predicts a much higher frequency of oscillation than does ΛCDM. The idea of using vertical motion to distinguish MOND from ΛCDM has been extensively investigated in the literature. Recently Angus et al have fit both rotation curves and vertical velocity dispersion for a sample of galaxies from the DiskMass survey [30]. Margalit and Shaviv [31] have argued that in MOND gravity stars undergoing vertical oscillation drift in their orbits in proportion to the square of the amplitude of their vertical oscillation; for simplicity in their analysis they approximate the quasilinear MOND field by the pristine field, an approximation that is justified by the results depicted in Fig.6. Earlier Bienayme et al [32] and Nipoti et al [33] have proposed Milky Way tests to distinguish nonlinear MOND predictions from ΛCDM. Recently models of dark matter have appeared that posit the formation of a disk of dark matter; these models also predict enhanced vertical oscillation frequencies [34,35]. It is likely that future surveys providing data on both in plane and vertical motion will be needed to discriminate among these different models.

IV. CONCLUSION
The radial acceleration relation is a plot of the observed radial acceleration in disk galaxies as a function of the acceleration that would be expected on the basis of the Newtonian gravitational field of the visible baryonic matter. Pristine MOND predicts that the relation is a smooth curve determined by the interpolating function f of the theory. Quasilinear MOND however predicts a scatter about the smooth curve resulting from Pris- tine MOND. The principal finding of this paper is that quasilinear MOND is nonetheless compatible with the observed radial acceleration relation because the predicted intrinsic scatter is small compared to the observational uncertainties. It appears a daunting task to try to relate the observed scatter about a smooth fit to the predicted intrinsic scatter. In principle one could imagine examining the deviations after sorting the galaxies in to groups according to parameters such as µ and ∆. However as already noted [29] the predicted deviations are even smaller than apparent from the plot in fig 5 and are moreover model dependent being sensitive to the unknown interpolating function f .

Appendix A: Newtonian field near galactic center
Near the origin the Newtonian potential corresponding to the distribution eq (9) can be expanded in a power series. Here we show the form of this series and a sequence of sum rules that the coefficients in the expansion must satisfy. By verifying that the sum rules are satisfied we are able to perform a non-trivial check on our numerical computations.
The Newtonian potential is related to the Newtonian field via g N = −∇φ N . We take φ N = 0 at the origin. Taking into account the azimuthal and reflection symmetries of the galaxy profile (9), to quartic order the potential must have the form φ N = A(x 2 +y 2 )+Bz 2 +C(x 2 +y 2 ) 2 +D(x 2 +y 2 )z 2 +Ez 4 +. . .
(A1) It is now a straightforward matter to compute the Laplacian of the Newtonian potential. Comparing to the corresponding series expansion of the density ρ D and making use of Poisson's equation for Newtonian gravity, For each galaxy, the accelerations are calculated at 36 points in the galactic plane at different distances from galactic center ranging from zero to two galactic radii. The accelerations are converted to units of m/s 2 assuming the value a0 = 1.2 × 10 −10 m/s 2 and are plotted on a logarithmic scale. All galaxies have ∆ = 0.5; mass parameters are µ = 0.01 (red), 0.05 (blue), 0.1 (cyan), 0.5 (purple), 1.0 (green), 5.0 (magenta) and 10.0 (orange). The grey curve shows the radial component of the pristine MOND acceleration plotted as a function of the Newtonian acceleration; the dashed line is the large field asymptote wherein the MOND accelerations equal the Newtonian. The deviations of the quasilinear MOND data from the pristine MOND curve confirms that the radial acceleration relation has intrinsic scatter in quasilinear MOND. However the scatter is small compared to observations uncertainties showing that the observed radial acceleration relation is compatible with quasilinear MOND. ∇ 2 φ N = 4πGρ, we deduce 4A + 2B = 4µ/∆ √ π, 16C + 2D = −4µ/∆ √ π, 4D + 12E = −4µ/∆ 3 √ π. (A2) Numerically the coefficients A, B, C, D and E are computed by differentiating the Newtonian field g N at the origin. We compute these derivatives in Fourier space. The numerically calculated coefficients satisfy the sum rules to a part in a thousand accuracy. The agreement confirms not only that our calculation of the Newtonian field is accurate near the origin but also that our scheme of computing space derivatives in Fourier space is accurate. This is reassuring since we do need to compute spatial derivatives in order to compute quantities of interest such as e.g. ∇ · g P and ∇ × g P .