Testing MOND in the Solar System
aa r X i v : . [ a s t r o - ph . C O ] M a y Testing MOND in the Solar System
Luc BLANCHET G R ε C O , Institut d’Astrophysique de Paris, CNRS,Universit´e Pierre et Marie Curie, 98 bis boulevard Arago, 75014 Paris, France J´erˆome NOVAK
Laboratoire Univers et Th´eories, Observatoire de Paris, CNRS,Universit´e Denis Diderot, 5 place Jules Janssen, 92190 Meudon, France
The Modified Newtonian Dynamics (MOND) generically predicts a violation of the strongversion of the equivalence principle. As a result the gravitational dynamics of a systemdepends on the external gravitational field in which the system is embedded. This so-calledexternal field effect is shown to imply the existence of an anomalous quadrupolar correction,along the direction of the external galactic field, in the gravitational potential felt by planets inthe Solar System. We compute this effect by a numerical integration of the MOND equationin the presence of an external field, and deduce the secular precession of the perihelion ofplanets induced by this effect. We find that the precession effect is rather large for outergaseous planets, and in the case of Saturn is comparable to, and in some cases marginallyexcluded by published residuals of precession permitted by the best planetary ephemerides.
The Modified Newtonian Dynamics (MOND) has been proposed as an alternative to the darkmatter paradigm . At the non-relativistic level, the best formulation of MOND is the modifiedPoisson equation , ∇ · (cid:20) µ (cid:18) ga (cid:19) ∇ U (cid:21) = − πGρ , (1)where ρ is the density of ordinary (baryonic) matter, U is the gravitational potential, g = ∇ U is the gravitational field and g = | g | its ordinary Euclidean norm. The modification of thePoisson equation is encoded in the MOND function µ ( y ) of the single argument y ≡ g/a , where a = 1 . × − m / s denotes the MOND constant acceleration scale. The MOND functioninterpolates between the MOND regime corresponding to weak gravitational fields y = g/a ≪ µ ( y ) = y + o ( y ), and the Newtonian strong-field regime y ≫
1, where µ reduces to 1 so that we recover the usual Newtonian gravity.An important consequence of the non-linearity of Eq. (1) in the MOND regime, is thatthe gravitational dynamics of a system is influenced (besides the well-known tidal force) bythe external gravitational environment in which the system is embedded. This is known asthe external field effect (EFE), which has non-trivial implications for non-isolated gravitatingsystems. The EFE was conjectured to explain the dynamics of open star clusters in our galaxy ,since they do not show evidence of dark matter despite the involved weak internal accelerations(i.e. below a ). The EFE effect shows that the dynamics of these systems should actually beewtonian as a result of their immersion in the gravitational field of the Milky Way. The EFEis a rigorous prediction of the equation (1), and is best exemplified by the asymptotic behaviourof the solution of (1) far from a localised matter distribution (say, the Solar System), in thepresence of a constant external gravitational field g e (the field of the Milky Way). At largedistances r = | x | → ∞ we have U = g e · x + GM/µ e r p λ e sin θ + O (cid:18) r (cid:19) , (2)where M is the mass of the localised matter distribution, where θ is the polar angle from thedirection of the external field g e , and where we denote µ e ≡ µ ( y e ) and λ e ≡ y e µ ′ e /µ e , with y e = g e /a and µ ′ e = d µ ( y e ) / d y e . In the presence of the external field, the MOND internalpotential u ≡ U − g e · x shows a Newtonian-like fall-off ∼ r − at large distances but with aneffective gravitational constant G/µ e . a However, contrary to the Newtonian case, it exhibits anon-spherical deformation along the direction of the external field. The fact that the externalfield g e does not disappear from the internal dynamics can be interpreted as a violation of thestrong version of the equivalence principle. In two recent papers , it was shown that the imprint of the external galactic field g e on the SolarSystem (due to a violation of the strong equivalence principle) shows up not only asymptotically,but also in the inner regions of the system, where it may have implications for the motion ofplanets. This is somewhat unexpected because gravity is strong there (we have g ≫ a ) andthe dynamics should be Newtonian. However, because of the properties of the equation (1), thesolution will be given by some non-local Poisson integral, and the dynamics in the strong-fieldregion will be affected by the anomalous behaviour in the asymptotic weak-field region.We assume that the external Galactic field g e is constant over the entire Solar System. b The motion of planets of the Solar System relatively to the Sun obeys the internal gravitationalpotential u defined by u = U − g e · x , (3)which is such that lim r →∞ u = 0. Contrary to what happens in the Newtonian case, the externalfield g e does not disappear from the gravitational field equation (1) and we want to investigatenumerically its effect. The anomaly detected by a Newtonian physicist is the difference ofinternal potentials, δu = u − u N , (4)where u N denotes the ordinary Newtonian potential generated by the same ordinary matterdistribution ρ , and thus solution of the Poisson equation ∆ u N = − πGρ with the boundarycondition that lim r →∞ u N = 0. We neglect here the change in the matter distribution ρ whenconsidering MOND theory instead of Newton’s law. This is in general a good approximationbecause the gravitational field giving the hydrostatic equilibrium (and thus ρ ) is strong andMOND effects are very small. Hence u N is given by the standard Poisson integral.A short calculation shows that the anomaly obeys the Poisson equation ∆ δu = − πGρ pdm ,where ρ pdm is the density of “phantom dark matter” defined by ρ pdm = 14 πG ∇ · ( χ ∇ U ) , (5) a Recall that in the absence of the external field the MOND potential behaves like U ∼ −√ GMa ln r , showingthat there is no escape velocity from an isolated system . However since no object is truly isolated the asymptoticbehaviour of the potential is always given by (2), in the approximation where the external field is constant. b For the Milky Way field at the level of the Sun we have g e ≃ . × − m / s which happens to be slightlyabove the MOND scale, i.e. η ≡ g e /a ≃ . here we denote χ ≡ µ −
1. The phantom dark matter represents the mass density thata Newtonian physicist would attribute to dark matter. In the model , the phantom darkmatter is interpreted as the density of polarisation of some dipolar dark matter medium and thecoefficient χ represents the “gravitational susceptibility” of this dark matter medium.The Poisson equation ∆ δu = − πGρ pdm is to be solved with the boundary condition thatlim r →∞ δu = 0; hence the solution is given by the Poisson integral δu ( x , t ) = G Z d x ′ | x − x ′ | ρ pdm ( x ′ , t ) . (6)We emphasise that, contrary to the Newtonian (linear) case, the knowledge of the matter densitydistribution does not allow to obtain an analytic solution for the potential, and the solution hasto be investigated numerically. We can check that the phantom dark matter behaves like r − when r → ∞ , so the integral (6) is perfectly convergent.In the inner part of the Solar System the gravitational field is strong ( g ≫ a ) thus µ tendsto one there, and χ tends to zero. Here we adopt the extreme case where χ is exactly zero in aneighbourhood of the origin, say for r ε , so that there is no phantom dark matter for r ε ;for the full numerical integration later we shall still make this assumption by posing χ = 0inside the Sun (in particular we shall always neglect the small MOND effect at the centre ofthe Sun where gravity is vanishingly small). If ρ pdm = 0 when r ε we can directly obtainthe multipolar expansion of the anomalous term (6) about the origin by Taylor expanding theintegrand when r = | x | →
0. In this way we obtain c δu = + ∞ X l =0 ( − ) l l ! x L Q L , (7)where the multipole moments near the origin are given by Q L = G Z r>ε d x ρ pdm ∂ L (cid:18) r (cid:19) . (8)Because the integration in (8) is limited to the domain r > ε and ∂ L (1 /r ) is symmetric-trace-free(STF) there [indeed ∆(1 /r ) = 0], we deduce that the multipole moments Q L themselves areSTF. This can also be immediately inferred from the fact that ∆ δu = 0 when r ε , hencethe multipole expansion (7) must be a homogeneous solution of the Laplace equation which isregular at the origin, and is therefore necessarily made solely of STF tensors of type ˆ x L . Hencewe can replace x L in (7) by its STF projection ˆ x L . It is now clear from the non-local integral in(8) that the MONDian gravitational field (for r > r ) can influence the near-zone expansion ofthe field when r →
0. An alternative expression of the multipole moments can also be proved,either directly or by explicit transformation of the integral (8). We have Q L = − u N ( ) δ l, + ( − ) l ( ˆ ∂ L u )( ) , (9)where the Newtonian potential u N and the STF derivatives of the internal potential u are to beevaluated at the centre of the Sun.The multipole expansion (7) will be valid whenever r is much less than the MOND transitiondistance for the Solar System, defined by r = p GM/a with M the mass of the Sun and a the c Our notation is as follows: L = i · · · i l denotes a multi-index composed of l multipolar spatial indices i , · · · , i l (ranging from 1 to 3); ∂ L = ∂ i · · · ∂ i l is the product of l partial derivatives ∂ i ≡ ∂/∂x i ; x L = x i · · · x i l is theproduct of l spatial positions x i ; similarly n L = n i · · · n i l = x L /r l is the product of l unit vectors n i = x i /r ; thesymmetric-trace-free (STF) projection is indicated with a hat, for instance ˆ x L ≡ STF[ x L ], and similarly for ˆ n L and ˆ ∂ L . In the case of summed-up (dummy) multi-indices L , we do not write the l summations from 1 to 3 overtheir indices.igure 1: Left panel: profile of Q ( r ) in the Solar System, for a standard choice of function µ ( y ) [see Eq. (13a)], a = 1 . × − m.s − and g e = 1 . × − m.s − . The MOND transition radius is shown by a dash-dotted lineat r ≃ r
50 AU), where the quadrupole is almost constant.
MOND acceleration scale. This radius corresponds to the transition region where the Newtonianacceleration becomes of the order of the MOND acceleration a and therefore, MOND effectsbecome dominant. We have r ≃ So far we have elucidated the structure of the multipole expansion of the anomaly δu near theorigin. Next we resort to a numerical integration of the non-linear MOND equation (1) in orderto obtain quantitative values for the multipole moments. d The Sun being assumed to be spherically symmetric, since all the multipole moments areinduced by the presence of the external field g e in the preferred direction e , the situation isaxisymmetric and all the moments Q L will have their axis pointing in that direction e . Thuswe can define some multipole coefficients Q l by posing Q L = Q l ˆ e L , where ˆ e L denotes the STFpart of the product of l unit vectors e L = e i · · · e i l . The multipole expansion (7) reads then as δu ( r, θ ) = + ∞ X l =0 ( − ) l (2 l − r l Q l ( r ) P l (cos θ ) , (10)where P l ( z ) is the usual Legendre polynomial and θ is the angle away from the Galactic di-rection e . Although from the previous considerations the multipole coefficients Q l should beapproximately constant within the MOND transition radius r , here we compute them directlyfrom the numerical solution of (1) and shall obtain their dependence on r . With our definitionthe quadrupolar piece in the internal field is given by δu = 12 r Q ( r ) (cid:18) cos θ − (cid:19) . (11)The radial dependence of the anomaly (11) is ∝ r and can thus be separated from a quadrupolardeformation due to the Sun’s oblateness which decreases like ∝ r − .As a first result, we show in Fig. 1 the profile of the quadrupole induced by the MOND theorythrough the function Q ( r ) defined in Eq. (11). We find that this quadrupole is decreasing from d Our numerical scheme is based on the very efficient integrator of elliptic equations lorene , available fromthe website .able 1: Numerical values of the quadrupole Q together with the associated dimensionless quantity q definedby Eq. (12). All values are given near the Sun. We use different choices of the function µ ( y ) defined in Eqs. (13). MOND function µ ( y ) µ ( y ) µ ( y ) µ exp ( y ) µ TeVeS ( y ) Q [s − ] 3 . × − . × − . × − . × − . × − q .
33 0 .
19 1 . × − .
26 0 . Q ( r ) is almost constant in a large sphere surrounding the Solarsystem, as it has a relative variation lower than 10 − within 30 AU (see the zoomed region inFig. 1). We shall therefore refer to the quadrupole as a simple number, noted Q (0) or simply Q , when evaluating its influence on the orbits of Solar-system planets.On dimensional analysis we expect that the quadrupole coefficient Q should scale with theMOND acceleration a like Q = a r q ( η ) , (12)where r = p GM/a is the MOND transition radius and where the dimensionless coefficient q depends on the ratio η = g e /a between the external field and a , and on the choice of theinterpolating function µ . Our numerical results for the quadrupole are given in Table 1, fordifferent coupling functions µ ( y ). e Here we consider various cases widely used in the literature: µ n ( y ) = y n √ y n , (13a) µ exp ( y ) = 1 − e − y , (13b) µ TeVeS ( y ) = √ y − √ y + 1 . (13c)The function µ has been shown to yield good fits of galactic rotation curves ; However becauseof its slow transition to the Newtonian regime it is a priori incompatible with Solar Systemobservations. The function µ is generally called the “standard” choice and was used in fits .We include also the function µ exp having an exponentially fast transition to the Newtonianregime. The fourth choice µ TeVeS is motivated by the TeVeS theory . One should note thatnone of these functions derives from a fundamental physical principle.We have used several functions of type µ n , as defined in Eq. (13a). One can notice thatthe value of Q decreases with n , that is with a faster transition from the weak-field regimewhere µ ( y ) ∼ y , to the strong field regime where µ ( y ) ∼
1. We have been unable to determinenumerically a possible limit for Q as n goes to infinity. We investigate the consequence for the dynamics of inner planets of the Solar System of thepresence of an abnormal quadrupole moment Q oriented toward the direction e of the galacticcentre. Recall that the domain of validity of this anomaly is expected to enclose all the inner SolarSystem (for distances r . r ≈ u = u N + δu , where u N = GM/r and the perturbation function R ≡ δu is given for the quadrupole moment by Eq. (11). e Note that the quadrupole coefficient Q is found to be always positive which corresponds to a prolate elon-gation along the quadrupolar axis. e apply the standard linear perturbation equations of celestial mechanics . The unper-turbed Keplerian orbit of a planet around the Sun is described by six orbital elements. Forthese we adopt the semi-major axis a , the eccentricity e , the inclination I of the orbital plane,the mean anomaly ℓ defined by ℓ = n ( t − T ) where n = 2 π/P ( n is the mean motion, P is theorbital period and T is the instant of passage at the perihelion), the argument of the perihelion ω (or angular distance from ascending node to perihelion), and the longitude of the ascendingnode Ω. We also use the longitude of the perihelion defined by ˜ ω = ω + Ω.The perturbation function R = δu is a function of the orbital elements of the unperturbedKeplerian ellipse, say { c A } = { a, e, I, ℓ, ω, Ω } . The perturbation equations are generated bythe partial derivatives of the perturbation function with respect to the orbital elements, namely ∂R/∂c A . We express the planet’s absolute coordinates ( x, y, z ) (in some absolute Galilean frame)in terms of the orbital elements { a, e, I, ℓ, ω, Ω } by performing as usual three successive framerotations with angles Ω, I and ω , to arrive at the frame ( u, v, w ) associated with the motion,where ( u, v ) is in the orbital plane, with u in the direction of the perihelion and v oriented inthe sense of motion at perihelion. The unperturbed coordinates of the planet in this frame are u = a (cos U − e ) , (14a) v = a p − e sin U , (14b) w = 0 , (14c)where U denotes the eccentric anomaly, related to ℓ by the Kepler equation ℓ = U − e sin U .The perturbation equations provide the variations of the orbital elements d c A / d t as linear com-binations of the partial derivatives ∂R/∂c B of the perturbation function. We are interested onlyin secular effects, so we average in time the perturbation equations over one orbital period P .Denoting the time average by brackets, and transforming it to an average over the eccentricanomaly U , we have (cid:28) d c A d t (cid:29) = 1 P Z P d t d c A d t = 12 π Z π d U (1 − e cos U ) d c A d t . (15)In the following, to simplify the presentation, we shall choose the x -direction of the absoluteGalilean frame to be the direction of the galactic centre e = g e /g e . That is, we assume thatthe origin of the longitude of the ascending node Ω lies in the direction of the galactic centre.Furthermore, in order to make some estimate of the magnitude of the quadrupole effect, let usapproximate the direction of the galactic centre (which is only 5 . I = 0. In this case˜ ω = ω + Ω is the relevant angle for the argument of the perihelion. We then find the followingnon-zero evolution equations: (cid:28) d e d t (cid:29) = 5 Q e √ − e n sin(2˜ ω ) , (16a) (cid:28) d ℓ d t (cid:29) = n − Q n h e + 15(1 + e ) cos(2˜ ω ) i , (16b) (cid:28) d˜ ω d t (cid:29) = Q √ − e n h ω ) i . (16c)We recall that ˜ ω is the azimuthal angle between the direction of the perihelion and that of thegalactic centre (approximated to lie in the orbital plane). Of particular interest is the secularprecession of the perihelion h d˜ ω/ d t i due to the quadrupole effect henceforth denoted by∆ = Q √ − e n h ω ) i . (17) able 2: Results for the precession rates of planets ∆ due to the quadrupole coefficient Q . We use the values for Q for various MOND functions as computed in Table 1. Published postfit residuals of orbital precession (aftertaking into account the relativistic precession). All results are given in milli-arc-seconds per century. Quadrupolar precession rate ∆ in mas / cyMOND function Mercury Venus Earth Mars Jupiter Saturn Uranus Neptune µ ( y ) 0 .
04 0 .
02 0 . − . − .
12 5 . − .
14 7 . µ ( y ) 0 .
02 0 .
01 0 . − . − .
65 3 . − .
87 4 . µ ( y ) 2 × − − × − − × − − .
06 0 . − .
56 0 . µ exp ( y ) 0 .
03 0 .
02 0 . − . − .
88 4 . − .
01 6 . µ TeVeS ( y ) 0 .
05 0 .
02 0 . − . − .
21 5 . − .
94 8 . h d˜ ω/ d t i in mas / cyOrigin Mercury Venus Earth Mars Jupiter Saturn Uranus NeptunePitjeva − . ± − . ± . − . ± . . ± . − ± et al. − ± − ± ± .
016 0 ± . ± − ± ± · ± · Fienga et al. . ± . . ± . − . ± . − . ± . − ±
42 0 . ± .
65 - -
The precession is non-spherical, in the sense that it depends on the orientation of the orbitrelative to the galactic centre through its dependence upon the perihelion’s longitude ˜ ω . Theeffect scales with the inverse of the orbital frequency n = 2 π/P and therefore becomes moreimportant for outer planets like Saturn than for inner planets like Mercury. This is in agreementwith the fact that the quadrupole effect we are considering increases with the distance to theSun (but of course will fall down when r becomes appreciably comparable to r , see Fig. 1).Our numerical values for the quadrupole anomalous precession ∆ are reported in Table 2.As we see the quadrupolar precession ∆ is in the range of the milli-arc-second per centurywhich is not negligible. In particular it becomes interestingly large for the outer gaseous planetsof the Solar System, essentially Saturn, Uranus and Neptune. The dependence on the choice ofthe MOND function µ is noticeable only for functions µ n ( y ) defined by (13a) with large valuesof n , where the effect decreases by a factor ∼
10 between n = 2 and n = 20.We then compare in Table 2 our results to the best published postfit residuals for any pos-sible supplementary precession of planetary orbits (after the relativistic precession has beenduly taken into account), which have been obtained from global fits of the Solar System dynam-ics , , . In particular the postfit residuals obtained by the INPOP planetary ephemerides , use information from the combination of very accurate tracking data of spacecrafts orbiting dif-ferent planets. We find that the values for ∆ are smaller or much smaller than the publishedresiduals except for the planets Mars and Saturn. Very interestingly, our values are smaller orgrossly within the range of the postfit residuals for these planets. In the case of Saturn notably,the constraints seem already to exclude most of our obtained values for ∆ , except for MONDfunctions of the type µ n and given by (13a) with rather large values of n .However let us note that the INPOP ephemerides are used to detect the presence of aneventual abnormal precession, not to adjust precisely the value of that precession , . On theother hand the postfit residuals are obtained by adding by hands an excess of precession forthe planets and looking for the tolerance of the data on this excess , . But in order to reallytest the anomalous quadrupolar precession rate ∆ , one should consistently work in a MONDpicture, i.e. consider also the other effects predicted by this theory, like the precession of thenodes, the variation of the eccentricity and the inclination, and so on — see Eqs. (16). Thenone should perform a global fit of all these effects to the data; it is likely that in this way thequantitative conclusions would be different.Finally let us cautiously remark that MOND and more sophisticated theories such as TeVeS ,which are intended to describe the weak field regime of gravity (below a ), may not be extrap-olated without modification to the strong field of the Solar System. For instance it has beenrgued that a MOND interpolating function µ which performs well at fitting the rotation curvesof galaxies is given by µ defined by (13a). However this function has a rather slow transition tothe Newtonian regime, given by µ ∼ − y − when y = g/a → ∞ , which is already excludedby Solar System observations. Indeed such slow fall-off − y − predicts a constant supplementaryacceleration directed toward the Sun δg N = a (i.e. a “Pioneer” effect), which is ruled out be-cause not seen from the motion of planets. Thus it could be that the transition between MONDand the Newtonian regime is more complicated than what is modelled by Eq. (1). This is alsotrue for the dipolar dark matter model , which may only give an effective description valid inthe weak field limit and cannot be extrapolated as it stands to the Solar System. While lookingat MOND-like effects in the Solar System we should keep the previous proviso in mind. Thepotential conflict we find here with the Solar System dynamics (notably with the constraints onthe orbital precession of Saturn , ) may not necessarily invalidate those theories if they arenot “fundamental” theories but rather “phenomenological” models only pertinent in a certainregime.In any case, further studies are to be done if one wants to obtain more stringent conclusionsabout constraints imposed by Solar-system observations onto MOND-like theories. More preciseobservations could give valuable informations about an eventual EFE due to the MOND theoryand restrict the number of possible MOND functions that are compatible with the observations.More generally the influence of the Galactic field on the Solar-system dynamics through a possi-ble violation of the strong version of the equivalence principle (of which the EFE is a by-productin the case of MOND) is worth to be investigated.1. M. Milgrom. Astrophys. J. , 270:365, 1983; ibid. ibid.
Ann. Rev. Astron. Astrophys. , 40:263, 2002.3. J.D. Bekenstein and M. Milgrom.
Astrophys. J. , 286:7, 1984.4. B. Famaey, J.-P. Bruneton, and H. Zhao. Escaping from modified newtonian dynamics.
Mon. Not. Roy. Astron. Soc. , 377:L79, 2007.5. M. Milgrom. Mond effects in the inner solar system.
Mon. Not. Roy. Astron. Soc. ,399:474, 2009.6. L. Blanchet and J. Novak. External field effect of modified newtonian dynamics in thesolar system.
Mon. Not. Roy. Astron. Soc. , 412:2530, 2011.7. L. Blanchet and A. Le Tiec. Model of dark matter and dark energy based on gravitationalpolarization.
Phys. Rev. D , 78:024031, 2008.8. L. Blanchet and A. Le Tiec. Dipolar dark matter and dark energy.
Phys. Rev. D ,80:023524, 2009.9. B. Famaey and J. Binney.
Mon. Not. Roy. Astron. Soc. , 363:603, 2005.10. K. Begeman, A. Broeils, and R. Sanders.
Mon. Not. Roy. Astron. Soc. , 249:523, 1991.11. J.D. Bekenstein.
Phys. Rev. D , 70:083509, 2004.12. D. Brouwer and G.M. Clemence.
Methods of celestial mechanics . Academic Press, NewYork, 1961.13. E.V. Pitjeva.
Astron. Lett. , 31:340, 2005.14. A. Fienga, J. Laskar, P. Kuchynka, C. Leponcin-Laffite, H. Manche, and M. Gastineau.Gravity tests with inpop planetary ephemerides. 2009. In
Relativity in FundamentalAstronomy: Dynamics, Reference Frames, and Data Analysis , edited by S. Klioner et al ,Cambridge University Press, p. 159.15. A. Fienga, H. Manche, P. Kuchynka, J. Laskar, and M. Gastineau. Planetary and lu-nar ephemerides, inpop10a. 2010. The documentation of INPOP10A is available at