The radial acceleration relation and a magnetostatic analogy in quasilinear MOND

In Newtonian gravity, the acceleration due to gravity g_N obeys the field equations

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {{\bf{g}}}_{N}=-4\pi G\rho ({\bf{r}})\hspace{3mm}\mathrm{and}\hspace{3mm}{\rm{\nabla }}\times {{\bf{g}}}_{N}=0.\end{eqnarray} \tag{ 1 }$$

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

$$\begin{eqnarray}&&{{\bf{g}}}_{P}=f\left(\displaystyle \frac{{g}_{N}}{{a}_{0}}\right){{\bf{g}}}_{N}.\end{eqnarray} \tag{ 2 }$$

Here a₀ ∼ 10⁻¹⁰ m s⁻² 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 RAR; indeed the observed RAR may be regarded as a measurement of the MOND acceleration a₀ 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

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {{\bf{g}}}_{Q}={\rm{\nabla }}\cdot {{\bf{g}}}_{P}\equiv -4\pi G{\rho }_{\mathrm{eq}}\hspace{3mm}\mathrm{and}\hspace{3mm}{\rm{\nabla }}\times {{\bf{g}}}_{Q}=0.\end{eqnarray} \tag{ 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

$$\begin{eqnarray}&&f(x)=\sqrt{1+\displaystyle \frac{1}{x}}.\end{eqnarray} \tag{ 4 }$$

An alternative form, following [7], is to take

$$\begin{eqnarray}&&f(x)=\displaystyle \frac{1}{1-\exp (-\sqrt{x})}.\end{eqnarray} \tag{ 5 }$$

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, ${{\bf{g}}}_{Q}={{\bf{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 RAR 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 RAR provides the impetus for the detailed study of quasilinear MOND reported in this paper.

To this end it is convenient to define

$$\begin{eqnarray}&&{\bf{b}}={{\bf{g}}}_{P}-{{\bf{g}}}_{Q}\end{eqnarray} \tag{ 6 }$$

and

$$\begin{eqnarray}&&{\bf{j}}={\rm{\nabla }}\times {{\bf{g}}}_{P}.\end{eqnarray} \tag{ 7 }$$

It then follows from equation (3) that

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {\bf{b}}=0\hspace{3mm}\mathrm{and}\hspace{3mm}{\rm{\nabla }}\times {\bf{b}}={\bf{j}}.\end{eqnarray} \tag{ 8 }$$

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. Equation (8) reveals that ∇ × g_P may be interpreted as a current that sources the difference between g_P and g_Q.

Given g_P equation (8) provides an efficient way to calculate b and thereby ${{\bf{g}}}_{Q}$ as we will discuss below. Formally one can write down a solution to equation (8) expressing b in terms of j using the Biot–Savart law. Whereas the relation equation (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.

We now introduce a simple model which allows us to study the deviation from the RAR in disk galaxies. We take the galaxies to have a Gaussian profile

$$\begin{eqnarray}&&{\rho }_{D}=\displaystyle \frac{M}{{\pi }^{3/2}{\rm{\Delta }}{R}_{d}^{3}}\exp \left[-\displaystyle \frac{({x}^{2}+{y}^{2})}{{R}_{d}^{2}}\right]\exp \left[-\displaystyle \frac{{z}^{2}}{{{\rm{\Delta }}}^{2}{R}_{d}^{2}}\right].\end{eqnarray} \tag{ 9 }$$

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 ${\rm{\Delta }}\to 0$ to a flat disk. The characteristic radial acceleration scale within Newtonian gravity for this circumstance is ${GM}/{R}_{d}^{2}$. The MOND field of the Galaxy is therefore fully determined by the dimensionless parameters $\mu ={GM}/{R}_{d}^{2}{a}_{0}$ and Δ. We will work in a system of units wherein R_d = 1 and a₀ = 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 ∝ μ. Equation (1) has the simple solution

$$\begin{eqnarray}&&{\tilde{{\bf{g}}}}_{N}({\bf{k}})={\rm{i}}4\pi G\ {\bf{k}}\displaystyle \frac{1}{{k}^{2}}{\tilde{\rho }}_{D}({\bf{k}})\end{eqnarray} \tag{ 10 }$$

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\to \infty$ because of the inverse square falloff of the gravitational field⁴ . 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⁴. 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 equation (3) is now a simple matter of applying an algebraic function to g_N. b can be calculated by solving equation (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 equation (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 RAR it is sufficient to consider Δ = 0.5 as discussed further in section 3.2; 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 equation (9). To check the accuracy of the calculated Newtonian field at short distances we have performed a non-trivial check described in the appendix. 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.

We start with a discussion of the quantity ρ_eq introduced in equation (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. Figure 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 equation (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 (${\rho }_{\mathrm{eq}}\propto 1/\sqrt{r}$ for $r\to 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 figure 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 [27] since they are tangential to the present study.

Zoom In Zoom Out Reset image size

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 ${\rm{\nabla }}\times {{\bf{g}}}_{P}=f^{\prime} ({g}_{N}/{a}_{0})\ {{\bf{g}}}_{N}\times {\rm{\nabla }}({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. Figure 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.

Zoom In Zoom Out Reset image size

We turn now to a comparison of the Newtonian, pristine MOND and quasilinear MOND gravitational fields. The upper panels of figure 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 galactic 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\to 0$ and $r\to \infty$. 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 the appendix. The orange curves correspond to the pristine MOND acceleration. This is simply related to the Newtonian field via the algebraic function equation (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 figure 3 we have used the interpolating function in equation (4) but we have explicitly verified that the alternative interpolation equation (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 figure 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.

Zoom In Zoom Out Reset image size

Figure 4 shows the radial component of b the difference between the pristine and quasilinear fields plotted as a function of r. 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 figure 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.

Zoom In Zoom Out Reset image size

In figure 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 RAR of [7] we note that the Newtonian acceleration corresponds to the computed baryonic acceleration in [7] while the pristine and quasilinear accelerations correspond to the observed acceleration of [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 figure 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 figure 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 RAR. The deviations from the smooth pristine behavior are bigger for Δ = 0.1⁵ but even in this case they remain smaller than the observational uncertainties described in [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.

Zoom In Zoom Out Reset image size

Although not directly related to the RAR, 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 ${\omega }^{2}=-\partial {g}_{z}/\partial 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 figure 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 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 [28]. Margalit and Shaviv [29] 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 figure 6. Earlier Bienayme et al [30] and Nipoti et al [31] 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 [32, 33]. It is likely that future surveys providing data on both in plane and vertical motion will be needed to discriminate among these different models.

Zoom In Zoom Out Reset image size