Testing quasilinear modified Newtonian dynamics in the Solar System
TTesting quasilinear modified Newtonian dynamics in the Solar System
Pasquale Galianni, ∗ Martin Feix, Hongsheng Zhao, and Keith Horne SUPA, School of Physics & Astronomy, University of St Andrews, KY16 9SS, United Kingdom Department of Physics, Technion - Israel Institute of Technology, Technion City, 32000 Haifa, Israel
A unique signature of the modified Newtonian dynamics (MOND) paradigm is its peculiar be-havior in the vicinity of the points where the total Newtonian acceleration exactly cancels. In theSolar System, these are the saddle points of the gravitational potential near the planets. Typically,such points are embedded into low-acceleration bubbles where modified gravity theories `a la MONDpredict significant deviations from Newton’s laws. As has been pointed out recently, the Earth-Sunbubble may be visited by the LISA Pathfinder spacecraft in the near future, providing a uniqueoccasion to put these theories to a direct test. In this work, we present a high-precision model of theSolar System’s gravitational potential to determine accurate positions and motions of these saddlepoints and study the predicted dynamical anomalies within the framework of quasi-linear MOND.Considering the expected sensitivity of the LISA Pathfinder probe, we argue that interpolation func-tions which exhibit a “faster” transition between the two dynamical regimes have a good chanceof surviving a null result. An example of such a function is the QMOND analog of the so-calledsimple interpolating function which agrees well with much of the extragalactic phenomenology. Wehave also discovered that several of Saturn’s outermost satellites periodically intersect the Saturn-Sun bubble, providing the first example of Solar System objects that regularly undergo the MONDregime.
PACS numbers: 98.10.+z, 95.35.+d, 98.62.Dm, 95.30.Sf
I. INTRODUCTION
Since the beginning of the past century, the Solar Sys-tem (SS) has been a major lab for testing the validityof general relativity (GR). Ranging from the first evi-dence of Mercury’s perihelion precession [1] up to therecent tests of the Lense-Thirring effect [2, 3], the SShas played a key role in our understanding of the laws ofgravity. All SS tests of GR performed so far, however,have tested Einstein’s theory in the regime of moderatelystrong gravity, i.e. where the total Newtonian accelera-tion exceeds the typical values found at galactic and cos-mological scales by several orders of magnitude. At largerscales, however, it is well-known that GR fails to accountfor the dynamics of galaxies and galaxy clusters unless amassive component of dark matter (DM) [4, 5] is takeninto account. The current lack of knowledge about thecomposition and behavior of this dark component of theUniverse, which has recently increased with the study ofthe internal dynamics of dwarf galaxies [6], together withthe failing attempts to detect DM particles [7], motivatesresearch on different alternative gravity theories such as f ( R ) gravity [8], conformal gravity [9] and various for-mulations of the modified Newtonian dynamics (MOND)paradigm [10–13] like, for instance, tensor-vector-scalartheory (TeVeS) [14] and quasi-linear MOND [15].Alternative frameworks of the latter kind predict im-portant deviations from Newtonian mechanics in regionswhere the total gravitational acceleration is much smallerthan the value of a ≈ − m s − , i.e. g N (cid:28) a . In ∗ Electronic address: [email protected] most parts of the SS, the total gravitational accelera-tion exceeds this value by several orders of magnitude.As first noted in Ref. [16], however, there exist numer-ous points in the SS where the sum of all gravitationalfields exactly cancels out. These are the saddle points(SPs) of the Newtonian potential near the planets whichare close to, but do not coincide with the Lagrangian L1points (where also the centrifugal potential is taken intoaccount). These SPs are embedded into ellipsoidal low-acceleration regions (bubbles) where the gradient of thegravitational potential becomes very small, i.e. compara-ble to or smaller than a , providing an ideal environmentfor testing some of the predictions of this set of modifiedgravity theories.The LISA Pathfinder space mission (LPF) [17] is anESA project currently planned for launch in 2014. LPFwill test technologies that will be employed in the follow-ing LISA (Laser Interferometer Space Antenna) mission,which aims at the detection of gravitational waves. TheLPF will consist of two closely spaced test masses, whosedistance relative to each other will be measured by meansof an interferometer. The spacecraft will be shielded fromany non-gravitational influence such that the masses in-side the spacecraft will be in an almost perfect geodesicmotion. The LPF will therefore be able to measure tidalstresses to an unprecedented accuracy. The expected tra-jectory of LPF will be a Lissajous orbit around the L1point. With minimum corrections to its planned orbit,LPF may be redirected to the low-gravity region encas-ing the Earth-Sun or Moon-Sun SPs [18]. In previousworks, it has been claimed that the LPF instruments haveenough sensitivity to detect the extra tidal stress signalpredicted by MOND and TeVeS [16, 18]. As such, thisparticular mission is believed to provide a unique occa- a r X i v : . [ a s t r o - ph . E P ] J un sion to put such modified gravity frameworks to a directtest. The strength of the MONDian perturbations thatthe spacecraft will experience when crossing the Earth-Sun low-gravity bubble, however, depends strongly onthe transition speed between the Newtonian and MONDregimes, as well as on the intrinsic structure of whicheverMOND formulation is used. The full space of possible re-sults has not yet been explored in detail. In the presentcontribution, we use the recently formulated quasi-linearMOND framework (QMOND) [15] to make predictionsabout the expected deviations from Newtonian gravitythat would occur inside the low-gravity region encasingthe Earth-Sun and Moon-Sun SPs. Since the positionsand sizes of the low-acceleration bubbles depend on thespecific configuration of the SS at a certain instant, wealso address the problem of determining their positionsto very high precision. This is a necessary prerequisite toany realistic performance of a saddle-flyby experiment.The present paper is structured as follows: In Sec. II,we introduce the basic concepts regarding QMOND andlow-acceleration bubbles. In Sec. III, we present an accu-rate model of the SS Newtonian potential which is usedto compute the positions and three-dimensional motionsof the Earth-Sun, Moon-Sun and Saturn-Sun SPs. Fol-lowing the lines of previous work [16, 18], we then use ourNewtonian model to compute the QMOND phantom DMdistribution within the low-gravity bubbles surroundingthe SPs. In Sec. IV, we present and discuss numeri-cal solutions to the QMOND potential near the Earth-Sun SP. We further derive the extra tidal stress felt by aspacecraft near the saddle and readdress the question ofdetecting such signals with the LPF. Considering the SPbetween two massive bodies, we also derive semi-analyticsolutions for the potential which, together with a compar-ison to a full numerical treatment, are separately givenin an appendix. Finally, we conclude in Sec. V. II. LOW ACCELERATION BUBBLES IN THESOLAR SYSTEMA. Quasi-linear MOND
QMOND or (QUMOND) [15] is a recently proposedquasi-linear formulation of the original MOND paradigm[10–12]. Within this framework, the gravitational poten-tial Φ generated by a mass distribution ρ bary is a solutionof Poisson’s equation for the modified source density∆Φ = 4 πG ˆ ρ = − ∇ · [ ν ( g N /a ) g N ] , (1)where g N ≡ | g N | is the Newtonian acceleration andthe acceleration scale a ≈ . × − m/s − is Mil-grom’s constant. Here the free function ν ( y ) is relatedto the MOND interpolation function µ ( x ) according to ν ( y ) = 1 /µ ( x ) and y = xµ ( x ), where x = g/a is theMOND total gravity measured in units of a [19]. Inthe deep-MOND regime, i.e. g N (cid:28) a , we have that ν ( g N /a ) ≈ ( g N /a ) − / whereas ν → g N (cid:29) a , in which case Eq. (1)reduces to the usual Poisson equation. For later pur-poses, it is convenient to express the modified QMONDdensity ˆ ρ as the sum of the baryonic contribution ρ bary and a phantom dark matter (PDM) term ρ PDM :ˆ ρ = ρ bary + ρ PDM . (2)As has been discussed in Ref. [15], QMOND is derivablefrom an action and thus satisfies all the usual conserva-tion laws. Furthermore, it is important to note that, ashas been shown in Refs. [20, 21], our Eq. (1) is the gen-uine nonrelativistic limit of a particular class of bimetricgravity theories [22]. Contrary to TeVeS and similar pro-posed relativistic extensions, these constructions do not necessarily require a renormalization of the gravitationalconstant, meaning a difference between the bare value of G (present in the Lagrangian) or G c (the effective valueused in cosmology calculations) and the usual Newtoniangravitational constant measured on Earth. In Sec. II B,we will further elaborate on this and its implications forthe present analysis.The practical advantage of the QMOND formulationwith respect to classical MOND is the fact that it re-quires only solutions to linear partial differential equa-tions, with a further nonlinear algebraic step. In general,the QMOND density ˆ ρ can be interpreted as the mod-ified source density that yields a MOND-like potential.Using simple algebra, we can express Eq. (1) asˆ ρ = ν ( g N /a ) ρ bary − (4 πGa ) − ν (cid:48) ∇ g N · g N . (3)Near a SP (in the absence of baryonic matter), we have ρ bary = 0 and the above reduces toˆ ρ = ρ PDM = − (4 πGa ) − ν (cid:48) ∇ g N · g N , (4)where ν (cid:48) ≡ dν ( y ) dy (5)is the derivative of the interpolation function. The PDMconcept is useful to visualize the effects of nonlinear-ity in MOND using the standard Newtonian picture:For example, it was shown that, at galactic scales, thePDM may generate an effective DM disk [23, 24]. Otherwork demonstrated that placing a galaxy into an externalgraviational field leads to an hourglass-shaped distribu-tion of negative PDM density [25]. The symmetry axisof this hourglass-shaped distribution is parallel to thedirection of the external field (see, e.g., Fig. 3 of Ref.[25]). Presuming that DM particles always give rise toa positive dynamical density, this points toward one ofthe most important differences between (cold) DM andMOND at galactic scales. At solar system scales, the sit-uation is quite similar since, as we show in Sec. IV, ananalogous distribution of negative PDM is present nearthe SPs of the total Newtonian potential in the SS. B. QMOND bubbles
It has been claimed that a direct test of Newton’s lawsin the regime of ultra-low accelerations is feasible withour current technological capabilities. More precisely,it has been pointed out that a spacecraft, closely ap-proaching one of the low-gravity bubbles near the Earth,should be able to detect an additional tidal stress signalin MOND whose amplitude and shape has been predictedin the framework of TeVeS [16, 18]. Naively, one wouldexpect these bubbles to be relatively small since MONDeffects should become important only for accelerationssmaller or approximately equal to a . A simple Newto-nian calculation shows that the semi-major axis of theellipsoid defined by the condition g N = a is ≈ ≈ φ , i.e. Φ ≈ Φ N + φ ,where Φ N is the usual Newtonian potential (see App. C).Thus, even if the total potential Φ is still dominated byits Newtonian contribution, it is possible to have φ al-ready in the deep-MOND regime. This peculiarity hasbeen exploited in Ref. [16] to estimate a characteristicacceleration scale in TeVeS below which full MONDianbehaviour of φ should be triggered, a trig = (4 π/k ) a ,where k is the TeVeS scalar parameter. Dividing a trig by the tidal stress at a generic saddle point then yields acorresponding TeVeS bubble with major axis r trig = a trig T , (6)where T is the Newtonian tidal stress along the symmetryaxis in the two-body approximation, T = 2 GM ( d − r sp ) (cid:16) (cid:112) M/m (cid:17) , (7) r sp ≈ d (cid:112) m/M is the distance of the SP from the lighterbody (assuming M >> m ), and d denotes the separationbetween the two bodies. For example, the Newtoniantidal stress at the Earth-Sun SP is T ≈ × − s − .Fixing k = 0 .
03 throughout this work, we have a trig ≈ a and the above formulae allow one to estimate abubble size of 766 km for the Earth-Sun bubble, 280 kmfor the Moon-Sun, and 3 . × km for the Saturn-Sunbubble.As for the framework of QMOND and its relativisticbimetric formulation, we assume that the value of G c ,i.e. the effective gravitational constant that is relevantto cosmology (e.g. in nucleosynthesis calculations), isidentical to both the bare G (entering the theory’s La-grangian) and the value of G measured on Earth (i.e. G = G c = G bare ) [26, 27]. In the context of Milgrom’sgeneral bimetric MOND gravity [20, 21], such a situationexplicitly corresponds to setting α + β = 0 and β = 1.This leaves substantial freedom for choosing viable inter-polation functions ν ( y ) (cf. Ref. [28]), and also meansthat our Eq. (1) in Sec. II A is rigorously correct. Despite a first investigation of Friedmann-Robertson-Walker so-lutions [29], the full cosmological scenario in the frame-work of these theories remains to be worked out com-pletely, and there are presently no constraints on howfast ν ( y ) may approach the Newtonian limit. In the restof the paper, we therefore take G = G c = G bare , bearingin mind that new cosmological constraints on ν ( y ) mightarise in the future. An important consequence of thisassumption is that a characteristic bubble size akin toEq. (6) does not exist in QMOND. Indeed, the bubblesizes associated with a detectable anomalous stress sig-nal will sensitively depend on the assumed form of ν ( y ),a problem we will address in Sec. IV. However, we willnevertheless use r trig given by Eq. (6) to define bubblesizes when addressing the problem of their motions andevolution in Sec. III.A necessary prerequisite for any realistic performanceof the LPF flyby experiment is the precise knowledge ofthe positions and motions of the low-acceleration bub-bles near the Earth. Neglecting extragalactic gravita-tional fields (which typically produce accelerations of a ext << − m s − ), the total Newtonian gravitationalpotential at a point within the SS is due to the contri-butions of the Sun, the planets (as well as their satel-lites) and the gravitational field of the Milky Way. Theexquisite knowledge of the relative magnitudes of theirNewtonian potentials allows one a precise determinationof the SP positions within the SS and the shape of theisopotential surfaces nearby. An approximate estimateof perturbations induced by a constant external acceler-ation field of magnitude a p on the position of a SP isgiven by (see Ref. [16]) δr ≈ a p T . (8)Substituting values for the Earth-Sun SP, we find thatthe Moon induces a mean perturbation on the two-bodyposition of the SP of the order of δr ≈ δr ≈
10 km.However, the gravitational field generated by the MilkyWay also plays a role in determining the exact positionof the SP. Its gravitational acceleration across the SS canbe estimated assuming that the Sun’s orbit around thecenter of the Galaxy is in centrifugal equilibrium. Usingthis hypothesis, we have v /r ≈ . × − m s − (as-suming v = 220 km s − , r = 2 . × km), three ordersof magnitude less than the gravitational pull exerted byJupiter. This would translate into a positional shift ofthe SP of a few cm, and can therefore be neglected.From these back-of-the-envelope calculations, we con-clude that, within the three-body approximation the po-sition of the Earth-Sun SP can be tracked with an ac-curacy of ≈
10 km. A plausible impact parameter forthe LPF spacecraft has been estimated as b ≈
50 km[18]. As we will show in the next sections, however, thestrength of the MOND tidal stress signal decreases sub-stantially with the ”rapidity” of the transition betweenthe Newtonian and MOND regimes, which is determinedby the slope of ν ( y ). As a consequence of this, to testa wider range of interpolation functions, the LPF space-craft might need to approach the SP more closely thanpreviously thought. This motivates our choice to use amany-body code to calculate the true motion of the SPsto excellent accuracy given by our current knowledge ofthe SS. III. MOTION OF THE MOND BUBBLESA. The computational framework
We have used the symbolic calculus package
Mathe-matica to build a high accuracy model of the SS’s New-tonian potential and to compute the positions of theEarth-Sun, Moon-Sun and Saturn-Sun SPs. In com-puting the total Newtonian potential across the SS, ourmodel takes into account the planetary potential andthe contributions of the most massive natural satellites.For the positions of the planets, we use INPOP08 high-precision ephemeris [40] and adopt the most recent IAUestimates for the Gaussian gravitational constant, theplanet’s masses relative to the Sun and their gravitationalmoments. All adopted values and their correspondingreferences are listed in Table I. Furthermore, note thatthe oblateness of the Sun, Jupiter, Saturn and the Earthare taken into account when computing the potential.Assuming polar coordinates, the Newtonian potential ofa planet (centered on the planet with the equatorial planeat z = 0) is represented by the series φ N ( r, z ) = − Gmr (cid:34) − ∞ (cid:88) n =2 r n J n P n (cid:16) zr (cid:17)(cid:35) , (9)where J n the n -th gravitational moment of the planet, r = (cid:112) x + y + z and P n the Legendre polynomial oforder n . In our code the series is truncated at n = 2.As the distance of the other planets from the SPs is verylarge, for them also the terms of O ( r − ) can be neglected.Once all contributions to the Newtonian potential havebeen added together, our code calculates the positions ofthe SPs by numerically solving the equation ∇ φ N = 0using Newton’s tangents method. B. The Earth-Sun SP
In general, the SS SPs of the total Newtonian poten-tial follow non-inertial trajectories. An intuitive pictureof the characteristics of this motion can be obtained bystarting with the two-body approximation and consider-ing the effects of a perturber later. Consider two bodiesof masses M and m , respectively, with M (cid:29) m . Choosea corotating Cartesian coordinate frame centered on theorbiting object with the x -axis pointing toward the heav-ier body and call the xy -plane the ecliptic plane. In thisreference frame, the SP is stationary if the orbit of the (cid:64) days (cid:68) r s p (cid:64) k m (cid:68) FIG. 1. Distance of the Earth-Sun SP from the Earth duringa year: The gentle yearly variation is due to the eccentricityon the Earth’s orbit and its amplitude is (measured on thepicture) around ∆ r sp = 8800 km. The vertical peaks aredue to the moon’s perturbation. Their measured amplitudesoscillate between 4700 km < a < (cid:45) (cid:64) km (cid:68) Y (cid:64) k m (cid:68) FIG. 2. Epicycles described by the Earth-Sun SP around itsunperturbed position in a year (projected on the plane of theecliptic; here the xy -plane). lightest body is circular, and its distance from the centerof the orbit is fixed at r sp . If the orbit is instead ellip-tical, then a periodic oscillation of the SP along the linejoining the two bodies arises. This oscillation has thesame period as the orbiting body and an amplitude that TABLE I. Table of adopted constants: All values are expressed in SI units.
Name Symbol Value Uncertainty Reference
Gaussian gravitational constant k G ± GM (cid:12) × m s − ± × [31]Geocentric gravitational constant GM ⊕ × m s − ± × [32] M (cid:12) /M Mercury - 6.0236 × ± × [33] M (cid:12) /M Venus - 4.08523719 × ± × − [34] M (cid:12) /M Mars - 3.09870359 × ± × − [35] M (cid:12) /M Jupiter - 1.047348644 × ± × − [36] M (cid:12) /M Saturn - 3.4979018 × ± × − [37] M (cid:12) /M Uranus - 2.290298 × ± × − [38] M (cid:12) /M Neptune - 1.941226 × ± × − [39] is given by ∆ r sp = ∆ d (cid:114) mM , (10)where ∆ d is the difference between the apoapsis and pe-riapsis of the orbiting body. For the Earth-Sun SP, theabove equation gives ∆ r sp ≈ a ≈ r sp = 8800 km, veryclose to the value calculated with the help of Eq. (10).The vertical spikes caused by the Moon’s presence haveamplitudes varying between 4700 km < a <
50 100 150 200 250 300 350 (cid:45) (cid:45) (cid:45) (cid:64) days (cid:68) z (cid:64) k m (cid:68) FIG. 3. Distance of the Earth-Sun SP from the ecliptic plane:This component of the motion is due to the inclination of theMoon’s orbital plane with respect to the ecliptic. The peak-to-valley amplitude of this yearly oscillation is around 5000km. right-hand side. At new moon, the SP shifts toward theEarth, whereas it shifts toward the Sun at full moon.The epicycle’s amplitude is higher if new moon occursat perigee and lower if at apogee. The yearly motion ofthe Earth-Sun SP in the direction perpendicular to theecliptic plane is illustrated in Fig.3.
C. The Moon-Sun SP
Another SP close enough to the Earth for direct testingis the Moon-Sun saddle [18], sometimes also referred toas the “Earth-Moon” SP [16]. However, since the gravi-tational pull of the Sun on the Moon is higher than thatexerted by the Earth, it is better to look at it as a Moon-Sun SP heavily perturbed by the Earth [18]. The magni-
10 15 20 25 30 35 4024 00026 00028 00030 00032 00034 000 Time (cid:64) days (cid:68) r s p (cid:72) k m (cid:76) FIG. 4. Distance of the Moon-Sun SP from the Moon duringa sidereal month: The asymmetry of the bell-like curve is dueto the inclination of the Moon’s orbit. tude of the perturbation induced by the Earth on the po-sition of the bubble is ∆ r sp ≈ r sp ≈ (cid:45) (cid:45) (cid:45) (cid:45) (cid:45)
15 000 (cid:45)
10 000 (cid:45) (cid:64) Km (cid:68) Y (cid:64) K m (cid:68) FIG. 5. Epycicles described by the Moon-Sun SP around itsunperturbed position during a year: Here the Moon lies onthe left-hand side.
D. The Saturn-Sun SP
The Saturn-Sun low-acceleration bubble is one of thelargest in the SS. The semi-major axis of the MONDbubble is 3 . × km ≈ .
02 AU. The minimum dis-tance between Jupiter and Saturn is around 5 . × km.When this minimum is reached, Jupiter exerts a gravita-tional pull of a p ≈ . × − m s − on the Saturn-SunSP. The Newtonian two-body tidal stress at the Saturn-Sun SP calculated from Eq. (7) is T ≈ . × − s − .Using Eq. (8), we may therefore calculate the peak-to-valley amplitude of this perturbation which turns out tobe around 1 . × km. This is much lower than the vari-ation in distance due to the orbital eccentricity which is
50 100 150 200 250 300 35024 00026 00028 00030 00032 00034 00036 000 Time (cid:64) days (cid:68) r s p (cid:72) k m (cid:76) FIG. 6. Distance of the Moon-Sun SP from the Moon over ayear. (cid:45) (cid:45) (cid:45) (cid:64) days (cid:68) z (cid:64) K m (cid:68) FIG. 7. Distance of the Moon-Sun SP from the plane of theecliptic: This component of the motion is due to the incli-nation of the Moon’s orbital plane with respect to the eclip-tic. The peak-to-valley amplitude of this yearly oscillation isaround 3000 km. approximately 2 . × km. In Fig. 8, we show thedistance of the SP from Saturn during a Saturnian year.Considering the size of the Saturn-Sun MOND bub-ble, there is a considerable chance that minor bodies orsome of Saturn’s outer satellites periodically intersect it.Indeed, we have found that a subsample of the Norsegroup [41, 42], consisting of Surtur, Ymir, Fenrir andLoge, periodically go through the Saturn-Sun MONDbubble. These small satellites are characterized by retro-grade motion. The semi-major axes of their orbits varyfrom 2 . × km (Fenrir) to 2 . × km (Loge). Theirorbits are periodically intersected by the MOND bubblewhose distance from Saturn varies between 2 . × km to 2 . × km (see Fig. 8). Furthermore, one findsthat the satellites’ inclination with respect to the eclip-tic is low and varies from 3 ◦ (Surtur) to 16 ◦ (Fenrir).Since the MOND bubble stays always close to the eclipticplane, close encounters with these objects are generallyfavored. An example of such an encounter between Fen-rir and the Saturn-Sun bubble is shown in Fig. 9. For all (cid:64) days (cid:68) r s p (cid:64) AU (cid:68) FIG. 8. Distance of the Saturn-Sun saddle point from saturnduring a Saturn’s year: The amplitude of the mainly due tothe eccentricity of Saturn’s orbit is around 0 .
017 AU. (cid:64) days (cid:68) D (cid:64) AU (cid:68) FIG. 9. Distance of the Saturn’s moon Fenrir from the Saturn-Sun SP during a Saturnian year: The horizontal dashed lineshows the size of the MOND bubble according to Eq. (6).Note that 6212 days after the 2012 Jan 1, the Moon gets veryclose to the center of the bubble at a distance of only r min ≈ r trig ≈ . × km. of these satellites, we list their dates of encounter withthe Saturn-Sun bubble, their minimum distances fromthe SP as well as their crossing times in Table II.So far, the found subset of Saturn’s satellites comprisesthe only known objects in the nearby SS that periodicallyenter the MOND regime. Even though the MONDian in-fluence on the satellites’ trajectories is likely to be verysmall, it cannot be ruled out that the induced perturba-tions at each encounter may accumulate to a potentiallyobservable effect. Regarding a detailed analysis of theMONDian impact on the path of these satellites, we re-fer to a forthcoming paper. E. Error estimates
The most important sources of error on the computedpositions of the SPs are the uncertainties on the plan-
TABLE II. Norse’s satellites minimum distances to theSaturn-Sun MOND bubble.Satellite d min d min /r trig Date Crossing time (cid:0) km (cid:1) (dd/mm/yyyy) (days)Surtur 2.00 0.48 24/06/2015 64Ymir 0.65 0.16 09/02/2013 68Fenrir 1.43 0.43 03/01/2029 49Loge 0.60 0.17 16/03/2026 58 R e s i du a l s (cid:64) m (cid:68) FIG. 10. Monte-Carlo simulation output: This plot shows thedistance of the perturbed Earth-Sun SP from its unperturbedposition. The mean value of this data series is a good estima-tor of our final uncertainty on the position of the Earth-SunSP. ets’ positions and masses. We have used a Monte-Carlomethod to estimate uncertainties in the positions of theSPs. The masses of the planets are perturbed using ran-dom values drawn from normal distributions with meansat the quoted parameter values, and standard devia-tions matching those listed in Table I. For the error onthe planets’ position we use the uncertainties in the IN-POP08 ephemeris [40] which we assume to be roughly ≈ ≈ ≈
10 km forthe inner planets. With this method, we obtain an abso-lute error on the position of the Earth-Sun SP of ∆ E − S ≈ (2 − body)E − S ≈ r sp ≈ d (cid:112) m/M representingthe approximate distance of the Earth-Sun SP from theEarth-Moon system barycenter. Our results tell us thatthe major sources of error on the position of the Earth-Sun SP are the uncertainties on the GM ⊕ /GM (cid:12) ratioand on the AU. The contribution of the other planetsto the final error sum up to an amplitude of about 1 m.For the Saturn-Sun SP we obtain, in the same way, anuncertainty of ∆ Sat − S ≈ IV. QMOND POTENTIAL INSIDE THEBUBBLEA. Phantom dark matter density near the SP
As can be seen from Eq. (4), the PDM density insidea MOND bubble is entirely determined by the Newto-nian potential and the choice of a specific interpolationfunction. Using our model of the SS gravitational po-tential, we have calculated the PDM density distributionaround the Earth-Sun, Moon-Sun and Saturn-Sun SPsfor different choices of ν ( y ). In particular, we have usedthe QMOND analog of the so-called simple interpolationfunction, µ ( x ) = x/ ( x + 1), which takes the form ν ( y ) = 12 + (cid:114) y + 14 . (11)This functional form has been shown to provide a gooddescription of extragalactic phenomenology and allowsone, for instance, to fit the observed rotation curves of awide range of spiral galaxies very well [43]. In addition,we have considered the interpolation function ν ( y ) = 1 − k π + (cid:34) y + (cid:18) k π (cid:19) (cid:35) / , (12)which is the QMOND analog of the function used in Refs.[16, 18] and exhibits a behavior close to that of the func-tional form originally proposed by Bekenstein in the con-text of TeVeS [14]. Finally, we have also chosen an inter-polation function similar to Eq. (11), but giving rise toa ”faster” transition between Newtonian dynamics andMOND: ν ( y ) = 12 + (cid:18) y + 116 (cid:19) / . (13)This choice is motivated by recent claims that interpo-lation functions like Eq. (11), which still result in arelatively ”slow” transition between the two dynamicalregimes, could give rise to an anomalously high secularprecession of the outer planet’s perihelion, incompatiblewith published residuals [44, 45]. In Fig. 11, we plot thethree interpolation functions as a function of the gravi-tational field strength expressed in units of a . As canbe seen, the function defined by Eq. (12) approaches thevalue ν ( y ) = 1 (corresponding to the Newtonian regime)much more slowly than the other ones. Considering theirperformance on galactic scales, we briefly discuss all threeinterpolating functions in App. B where it is also demon-strated that the behavior of Eq. (13) at these scales isvery similar to that of the simple ν -function Eq. (11).As remarked earlier, the PDM density depends on ν (cid:48) ,which basically tells us that the PDM density within theMOND bubble is sensitive to the “transition” speed be-tween the dynamical regimes. In fact, choosing an inter-polation function with a fast enough transition, it is pos-sible to suppress non-Newtonian effects for accelerations Eq. (12)Eq. (13) Eq. (11) − − − − − y = g N / a ν ( y ) − FIG. 11. Interpolation functions ν as a function of y = g N / a :Solid, dotted and dashed lines correspond to Eq. (12), Eq.(11) and Eq. (13), respectively. Vertical dashed lines show thecharacteristic scales y = 1 and y = (4 π/k ) for k = 0 .
03. TheBekenstein-like interpolation function Eq. (12) gives rise to amuch slower transition between the two dynamical regimes. (cid:45) (cid:45)
200 0 200 40010 (cid:45) (cid:45) (cid:45) (cid:45) X (cid:72) km (cid:76) (cid:72) Π G (cid:76) (cid:215) (cid:200) Ρ pd m (cid:200) (cid:72) s (cid:76) (cid:45) FIG. 12. Logarithmic plot of the rescaled PDM density’sabsolute value 4 πGρ
PDM at a cut of y = 50 km near the Earth-Sun SP: The line styles are defined as in Fig. 11. Note thatthe PDM density decreases by several orders of magnitude forthe ν -functions with increased slope. above a . The PDM density near the SP is proportionalthe trace of the stress tensor, 4 πGρ PDM = tr( S ij ), where S ij = ∂ i ∂ j Φ, and provides a good order-of-magnitudeestimator for the maximum MOND tidal stress signalwithin a certain region - if one excludes values at specialpoints such as the SPs themselves.The above dependency on the slope of the interpolationfunction can be seen in Fig. 12 where we plot the PDMdensity at cuts of y = 50 km near the Earth-Sun SP. It isevident that the PDM density changes by several ordersof magnitude between the different choices of ν ( y ). We FIG. 13. PDM isodensity contours inside in a box of 5 × km centered on the Earth-Sun SP along the z = 0 plane andassuming ν given by Eq. (12): The positions of the Earth andthe Moon-Sun SP are marked with a large and a small dot, re-spectively. The arrow points towards the Sun and the dashedline indicates the ρ PDM = 0 contour. The PDM density isnegative within a conic region (deformed by the Moon’s pres-ence), with the symmetry axis perpendicular to the Earth-Sundirection, and positive elsewhere. The Earth is embedded intoa halo of positive PDM. Moreover, the PDM density gradientis bigger near the Moon-Sun SP. therefore expect that at the same distance (or passingtrajectory) from the SP, the corresponding tidal stresssignal will be much smaller for an interpolation functiongiven by Eq. (11) than for Eq. (12). We will confirmthis preliminary prediction by using numerical solutionsof the QMOND potential near the Earth-Sun SP below.In Fig. 13, we illustrate the PDM density distributioncalculated with the help of Eq. (12) and rescaled by afactor of 4 × πG within a box of 10 km edge lengthcentered on the Earth-Sun SP. As expected, the PDMdensity distribution around the SP exhibits cylindricalsymmetry, with the symmetry axis along the line joiningthe planet with the Sun. In the case of the Earth-SunSP, this symmetry is perturbed by the presence of theMoon. B. Multi-grid FFT solver
Since the QMOND PDM decreases quickly toward zerowhen moving away from the SP, the boundary valueproblem given in Eq.(1) may be solved by means of effi-cient numerical methods based on the Fast Fourier Trans-0form (FFT). To achieve this, we have implemented amodified version of the Poisson solver described in Ref.[46] which now incorporates a multi-grid scheme to ob-tain a higher resolution near the SP.The basic idea is the following: Beginning with a largephysical box and a coarse sampling of the underlying den-sity distribution on an equidistant grid, we solve for thepotential Φ assuming zero boundary conditions - sincewe are only interested in the purely non-Newtonian con-tribution Φ QM = Φ − Φ N . Keeping the same numberof grid points per dimension N , we then move to a gridof smaller physical size such that the previous solutionremains sufficiently accurate at the new boundaries. Afirst correction to the potential is then found by consider-ing the solution for the correction density field, which isobtained as the difference between the true PDM densityand the one calculated by interpolating the density fieldon the coarse grid. For this step, we again assume zeroboundary conditions and use a cubic spline for the in-terpolation. Finally, the above procedure is successivelyrepeated until the desired resolution in the central partsis reached, and adding the individual results yields thefull solution with the appropriate boundary conditions.Apart from the usual sanity checks (see, e.g., Ref. [46]),the method was applied to several input densities andthe corresponding results were directly compared to thoseobtained with the original single-grid solver. For our pur-poses, we have found that the relative differences betweenthese results are basically negligible, typically on the or-der of (cid:46) − . In App. A, we also present semi-analyticsolutions (both for the quasi-Newtonian and the deep-MOND regime) which are used for comparison to thefull numerical results.In the case of the Earth-Sun saddle, for instance, weset N = 448 and use three physical cubic boxes with sidesof 20000, 5000 and 1250 km, respectively. This yields amaximum resolution of approximately ∆ x max = 2 . , ∆ x max / ,
0) to avoid anynumerical difficulties arising from the divergence of thePDM density at the exact SP. Assuming these parame-ters, the numerical solution is found to be reliable for dis-tances (cid:38)
20 km from the SP (again see App. A). Compu-tationally, our spectral multi-grid method performs quitewell: Considering the above setup, the full solution andderived quantities such as the stress tensor are calculatedon a timescale of less than half an hour on an average ma-chine using a single processor.
C. Numerical results
Using the above multi-grid FFT code, we have numeri-cally solved Eq.(1) near the Earth-Sun saddle point. Theinput PDM densities have been calculated for the two-body (Earth-Sun) and the three-body (Earth-Sun-Moon)approximations. Hereafter, we adopt a reference framecentered on the SP, with the x -axis pointing along the (cid:45) (cid:45)
200 0 200 40001234567 x (cid:72) km (cid:76) S yy (cid:72) s (cid:76) (cid:45) x10 (cid:45) FIG. 14. Two-body QMOND tidal stress signal for z = 0along the lines y ≈
25 (solid line), y ≈
100 (dashed line) and y ≈
400 km (dotted line) inside the Earth-Sun MOND bubble:For the calculation, we have assumed the interpolation func-tion Eq. (12). At y ≈
25 ( y ≈ S yy ≈ . × − s − ( S yy ≈ − . × − s − ). Earth-Sun direction (increasing toward the Sun). Notethat within this reference frame, which is identical to theone used in Ref. [18], the LPF instrument will only besensitive to the S yy component of the stress tensor be-cause of certain design characteristics [47]. The resultsof our simulations for the Earth-Sun bubble are shownin Figs. 14, 15 and 16: Using the Bekenstein-like inter-polation function given by Eq. (12), Fig. 14 shows thetwo-body QMOND tidal stress signal for z = 0 at cutsof around y = 25, 100 and 400 ( ± . .
7) and a deeper centralwiggle. Note that this agrees well with the conclusions ofRef. [28] about “Type II” MOND-like theories and thesoftening effect of the TeVeS curl field.We have also computed the tidal stress signal using thesimple interpolation function Eq. (11). As expected fromthe earlier analysis of the PDM density, the tidal stresssignal decreases by approximately two orders of magni-tude in this case. From Fig. 15, one finds that the peakvalue of the tidal stress at z = 0 and y ≈
25 km is ap-proximately 8 × − s − . Including the Moon into thecalculation (see Fig. 16) has only a small impact on theestimated stress signal. We have also verified that thetidal stress signal becomes even lower ( S yy ≈ − s − at z = 0 and y ≈
25 km) if the “fast-transition” interpo-lation function Eq. (13) is used. In the next section, weshall discuss the predicted additional stress signal in thelight of the expected LPF sensitivity.1 (cid:45) (cid:45)
200 0 200 40002468 x (cid:72) km (cid:76) S yy (cid:72) s (cid:76) (cid:45) x10 (cid:45) FIG. 15. Same as Fig. 14, but now assuming the simpleinterpolation function Eq. (11): At y ≈
25 ( y ≈ S yy ≈ × − s − ( S yy ≈− . × − s − ). (cid:45) (cid:45)
200 0 200 400 (cid:45) (cid:72) km (cid:76) S yy (cid:72) s (cid:76) (cid:45) x10 (cid:45) FIG. 16. Three-body QMOND tidal stress signal within theEarth-Sun bubble for z = 0 along the line y ≈
25 km for dif-ferent lunar phases, computed using the simple interpolationfunction Eq. (11): The solid (dotted) line corresponds to thetidal stress signal at approximately full (new) moon. Notethat the presence of the Moon has only a slight impact on thesignal.
D. LPF sensitivity and the QMOND signal
Choosing the analog of the Bekenstein-like ν -functionEq. (12) in the framework of TeVeS, it was concludedthat the instrument onboard the LPF spacecraft shouldhave enough sensitivity to detect the additional MON-Dian tidal stress signal [16, 18, 28]. At a distancefrom the SP similar to the predicted impact parame-ter, this signal was estimated to be on the order of S yy ≈ − s − . As we have seen in the previous sec- tions, however, the expected signal strength is very sen-sitive to the specific form of ν ( y ), and one might questionwhether such a conclusion remains valid for other choices.In particular, we will be interested in the situation for in-terpolation functions like Eq. (11) which provide a gooddescription of the phenomenology on extragalactic scales.Apart from the instrumental characteristics, the sen-sitivity of the LPF probe to tidal stresses depends onthe frequency of the MOND signal [48]. When the LPFspacecraft crosses the Earth-Sun bubble, the frequencyscale of the MOND signal will be given by f = v/l b ,where v is the spacecraft’s velocity and l b denotes thecharacteristic length scale of variations in the MONDtidal stress. It can be seen from Figs. 14 and 15 thatthe QMOND tidal stress signal typically varies over alength of l b ≈
100 km for both choices of the interpola-tion function. Note that this scale length is the sameas the one adopted in Ref. [28]. Given a model forthe power spectral density (PSD) S h ( f ) of the noise,the expected signal-to-noise ratio ( S/N ) may be roughlyestimated from the signal’s characteristic amplitude h c , S/N ∼ h c / √ f S h [49]. To obtain a sensitivity thresholdfor the LPF, we set S/N ∼
1. Considering v = 1 . − [50] and a constant noise PSD with √ S h ≈ . × − s − / √ Hz [28], this gives f ≈ − Hz andyields a threshold on the order of 10 − s − . This shouldbe sufficient to detect the QMOND tidal stress signal forinterpolation functions like Eq. (12) whose amplitudelies between 2 × − and 8 × − s − at y = 50 km(see Fig. 14). For the simple function Eq. (11), however,the situation is different: Since the tidal stress signal de-creases by approximately two orders of magnitude, weexpect that the signal will not be detectable at the givenimpact parameter. Therefore, to have a chance of prob-ing these models, the LPF spacecraft might need to passthe Earth-Sun SP much closer than previously thought,but it is also possible that the modified signal remainscompletely undetectable.To get more insight on the above, we follow the linesof Ref. [28] and consider the signal’s amplitude spectraldensity (ASD) which is given by the square root of thePSD: P ( f ) = 2 T (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) + T/ − T/ h ( t ) exp( − πif t ) dt (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (14)where f is the frequency, h ( t ) = S yy ( vt, b,
0) correspondsto the signal, and b is the impact parameter, i.e. theminimum distance from the SP which is approached at t = 0. Fixing v = 1 . − [18, 28], we furthermorehave the integration period T which we choose as T ≈ . × s, in accordance with our numerical setup (seeSec. IV). To compare the signal ASD to that of the noise,we use the simplified LPF noise model discussed in Ref.[47] which assumes a constant baseline within 1 < f < . × − s − / √ Hz. Beyondthis range, we assume that the sensitivity degrades with1 /f and f for lower and higher frequencies, respectively.2 f [Hz] A S D ! s − / √ H z " FIG. 17.
Top panel:
The ASD of the MONDian tidal stresssignal for ν ( y ) given by Eq. (12) compared to the simplenoise model described in the text; assuming v = 1 . − ,we illustrate results for b ≈
45 (solid line), 25 (dashed line), 12(dotted line) and 1 km (dot-dashed line).
Bottom panel:
Sameas the top panel, but now assuming the simple interpolationfunction Eq. (11).
Note that although this simple noise model will suffice forobtaining realistic estimates, improved LPF noise modelshave recently been presented in Ref. [28] (see their Fig.6).For the interpolation function Eq. (12), the top panelof Fig. 17 illustrates the resulting ASDs of the signal,assuming different impact parameters between approxi-mately 45 and 1 km. Since the numerically calculatedstress becomes unreliable for distances (cid:46)
20 km from theSP (see App. A), our results for small impact parame-ters underestimate the signal at the high-frequency end,i.e. for f (cid:38) × − Hz. However, this correspondsto a frequency domain where the noise is dominant, andthus we expect this resolution effect to have no signif-icant impact on our conclusions. As can be seen fromthe figure,
S/N is typically on the order of 10 per √ Hz over a broad frequency range. Applying the techniquesof noise matched filtering (which are frequently used inthe search for gravitational waves) to the present prob-lem, it can further be shown [28, 49] that the maximal S/N , realized by an optimal filtering template, is givenby (
S/N ) opt = 2 (cid:90) ∞ df (cid:12)(cid:12)(cid:12) ˜ h ( f ) (cid:12)(cid:12)(cid:12) S h ( f ) / . (15)Here S h ( f ) again denotes the PSD of the noise and ˜ h ( f )is the Fourier transform of the signal. Using Eq. (15), the above scenario for the Bekenstein-like interpolationfunction Eq. (12) translates into integrated S/N ratiosof around 50 − b ≈
12 km and larger, this gives rise toan integrated
S/N of around 1 or less. At b ≈
45 km,for instance, we have
S/N ≈ .
5. Only in the case of b ≈ S/N crosses unity,taking a value of 2. To test models based on interpola-tion functions similar to Eq. (11), the LPF spacecraftwould therefore be required to pass the Earth-Sun SP ata distance on the order of a 1 km or less. Note, however,that since the expected integrated
S/N ratios are ratherclose to unity, these estimates are quite sensitive to thedetails of the true underlying noise model which is cur-rently unknown. Depending on the actual noise model,this could, in principal, give rise to exposing the modifiedsignal to a wider radial range or making it virtually unde-tectable. Also, considering potentially necessary impactparameters b (cid:46) S/N ratios are evenfurther suppressed, and the modified tidal stress signalwill basically become invisible to the LPF. Finally, notethat we have also considered the framework of TeVeS andaddressed the question whether there exist viable inter-polation functions which give rise to similar conclusions.For this purpose, we have used our previous choices of ν ( y ) to formally construct the corresponding functions ˆ µ for the TeVeS scalar field. A detailed description of thisprocedure and the implications of our results are sepa-rately discussed in App. C. V. CONCLUSIONS
Low-acceleration regions encasing the SPs of the grav-itational potential in the SS, may exhibit significant de-viations from Newton’s laws within the framework ofMOND or closely related theories. As the LPF space-craft may visit the Earth-Sun low-acceleration bubble inthe near future, this could provide a unique occasion toput MOND-like modifications of gravity to a direct test.While it is known that the quoted LPF sensitivity is goodenough to detect the non-Newtonian tidal stress signal incertain cases, it is still very important to explore the fullspace of theoretical models and possible outcomes of theplanned flyby experiment. Also, a necessary prerequi-site for any realistic performance of this and future ex-3periments, is the precise knowledge of the positions andmotions of such bubbles near the Earth.Using a high accuracy model of the SS Newtonian po-tential, we have determined the positions and tracked themotions of the Earth-Sun, the Moon-Sun and Saturn-SunSPs with an accuracy of ∆ E − S ≈ Sat − S ≈ b ≈
50 km, we obtain an approximately by a factor of1 . S/N ratios on the order of 1 or less. Althoughmoving to very small impact parameters on the order of (cid:46)
S/N ratios are even further sup-pressed, and the modified tidal stress signal will basicallybecome invisible to the LPF.Finally, we want to point out that our conclusions donot necessarily apply to other proposed MOND frame-works. Considering the original form of TeVeS, for in-stance, a transition corresponding to Eq. (13) yields anunacceptable interpolation function ˆ µ for the scalar field,which is briefly discussed in App. C. Moreover, furtherconstraints on the rapidity of ν ( y ) in QMOND mightarise once the full cosmological scenario of its associatedrelativistic bimetric theory has been worked out. ACKNOWLEDGMENTS
P.G. wants to thank Marco Friuli for useful discussions.The authors want to thank the anonymous referee, JoaoMagueijo and Mordehai Milgrom for useful discussions.M.F. is supported in part at the Technion by a fellowship from the Lady Davis Foundation.
Appendix A: Semi-analytic solutions near the SP
In the following, we shall consider the configurationof two bodies at distance d with masses M and m , re-spectively, where we assume that the line connecting thebodies coincides with the z -axis. Shifting the coordinateorigin to the SP and introducing spherical polar coordi-nates with z = r cos θ , the Newtonian potential in thevicinity of the saddle approximately takes the formΦ N = Cr (cid:0) − θ (cid:1) + const , (A1)where C = G d M (cid:32) M − √ M mM − m (cid:33) − + m (cid:32) − M − √ M mM − m (cid:33) − . (A2)If M (cid:29) m (as in case of the Earth-Sun SP, for instance),Eq. (A2) may be well approximated by C ≈ GM d (cid:32) (cid:114) Mm (cid:33) . (A3)Near the SP where ∆Φ N = 0, the QMOND contributionto the total gravitational potential Φ,Φ = Φ N + Φ QM , (A4)is associated with the phantom source distribution˜ ρ QM ≡ ∆Φ QM = ∇ · [ ν ( y ) ∇ Φ N ] . (A5)Although Eq. (A5) does generally not exhibit analyticsolutions, it is possible to find approximate expressionsfor certain cases, which we shall discuss in the followingsections.
1. Quasi-Newtonian domain
If the resulting dynamics is sufficiently close to theNewtonian limit, i.e. y = g N /a (cid:29)
1, we may expandthe interpolating function ν ( y ) in powers of y , ν ( y ) ≈ α y + α y + . . . . (A6)Substituting the above into Eq. (A5) and rememberingthat ∆Φ N = 0 near the saddle, we thus obtain∆Φ QM ≈ α ∇ · (cid:18) ∇ Φ N y (cid:19) + α ∇ · (cid:18) ∇ Φ N y (cid:19) + . . . ≡ ˜ ρ (1) QM + ˜ ρ (2) QM + . . . . (A7)4The linearity of Eq. (A7) allows one to solve the equationorder by order. Using Eq. (A1), one finds that the n th-order contribution to the source term can be written as˜ ρ ( n ) QM = − A n r n − θ (1 + 3 cos θ ) ( n +2) / , n ∈ N , (A8)where we have defined A n = nα n a n (2 C ) n − . (A9)Let us begin with the first-order contribution. As weare only interested in the saddle region, the approximateNewtonian potential in Eq. (A1) suggests an ansatz ofthe formΦ (1) QM ( r, θ ) = γ r β T (1) ( θ ) + const , (A10)where β and γ are assumed as constant. Inserting theabove expression into Eq. (A7) then fixes β = 1 and γ = A , leaving us with a single ordinary differentialequation for T (1) ( θ ),2 T (1) ( θ ) + cos θ sin θ ddθ T (1) ( θ ) + d dθ T (1) ( θ )= − − θ (1 + 3 cos θ ) / . (A11)For a physically relevant solution, we require that T (1) ( θ )is smooth and periodic on the interval [0 , π ] (see, e.g.,[16]). In this case, Eq. (A11) may numerically be solvedin terms of a Fourier series, yielding T (1) ( θ ) ≈ . − . θ + 0 . θ − . θ + . . . . (A12)The above series converges rather quickly and may al-ready be cut off after the first eight terms to achieve arelative accuracy of (cid:46) − .Next, we shall consider the second-order correction tothe potential Φ QM . It turns out that one cannot adoptthe same ansatz as in Eq. (A10) because it is impossibleto meet the smoothness condition of the angular partin this case. As we will see shortly, this can easily beremedied by considering an additional term proportionalto log( r/r s ) in the potential. Thus we start fromΦ (2) QM ( r, θ ) = γ (cid:20) r β T (2) ( θ ) − B log (cid:18) rr s (cid:19)(cid:21) + const . (A13)Again substituting this expression into Eq. (A7), we findthat β = 0 and γ = A . This requires the constant B tobe dimensionless and gives rise to the following equationfor T (2) ( θ ):cos θ sin θ ddθ T (2) ( θ ) + d dθ T (2) ( θ ) = − − θ (1 + 3 cos θ ) + B. (A14) Since only derivatives of T (2) ( θ ) appear in Eq. (A14),one may integrate once and finally obtain ddθ T (2) ( θ ) = 1sin θ (cid:20) θ θ − B cos θ + D − √
33 arctan (cid:16) √ θ (cid:17)(cid:35) , (A15)where D is an integration constant. Applying the peri-odicity and smoothness conditions to the above, a briefcalculation leads to D = 0 and B = 12 − √ π. (A16)Once again, an approximate solution of T (2) ( θ ) may befound by using a Fourier series approach, and the result-ing series, T (2) ( θ ) ≈ − . θ + 0 . θ − . θ + 0 . θ + . . . , (A17)quickly converges toward the real solution. In principle,one could apply the general procedure outlined in thissection also to higher orders, but for the purpose of thispaper, we shall not take these into account. Finally, notethat the semi-analytic solutions found in this section au-tomatically satisfy the desired Newtonian boundary con-ditions because of β , β <
2, meaning that the totalgravitational potential is entirely dominated by the ex-pression in Eq. (A1) when moving to radii much largerthan the true size of the bubble, i.e. r/r (cid:29)
1, where r = r trig for the interpolation function defined in Eq.(12).
2. Deep-MOND domain
In the limit of very small gravitational fields, we have y → ν ( y ) ≈ √ y , (A18)which is independent of the particularly chosen form of ν ( y ). Using Eq. (A18) in Eq. (A5), the correspondingphantom source term can be expressed as˜ ρ QM = − A √ r − θ (1 + 3 cos θ ) / , (A19)where A = 12 (cid:112) Ca . (A20)Similar to our approach in the quasi-Newtonian situa-tion, we make the ansatzΦ QM ( r, θ ) = γr β T ( θ ) + const (A21)5which is again substituted into Eq. (A5). This time weobtain β = 3 / γ = A , with the resulting differentialequation for T ( θ ) given by154 T ( θ ) + cos θ sin θ ddθ T ( θ ) + d dθ T ( θ )= − − θ (1 + 3 cos θ ) / . (A22)Analogously to the first-order contribution in the quasi-Newtonian domain, an approximate solution is deter-mined by expanding T ( θ ) into a Fourier series on theinterval [0 , π ]. A straightforward calculation then yields T ( θ ) ≈ − . − . θ + 0 . θ − . θ + . . . , (A23)where it is sufficient to consider the first six terms of Eq.(A23) for an accuracy of (cid:46) − . It is important to notethat Eq. (A21) merely constitutes a special solution toEq. (A5), and thus we still need to address the questionof boundary conditions in the deep-MOND case.Given a suitable choice of such boundary conditions,Eq. (A21) needs to be replaced byΦ QM ( r, θ ) = Ar / T ( θ ) + Φ ( r, θ ) + const , (A24)where Φ ( r, θ ) is an appropriate solution of Laplace’sequation. Using the previously assumed smoothness andperiodicity condition of the angular part and additionallyrequiring that the potential is regular at r = 0, the mostgeneral form of Φ ( r, θ ) can be expressed asΦ ( r, θ ) = ∞ (cid:88) n =1 q n C n − a n − r n P n (cos θ ) , n ∈ N , (A25)where C is defined according to Eq. (A2), P n denotes theLegendre polynomial of degree n and q n ∈ R . Since thevalidity of equation (A24) is restricted to a close vicinityof the SP at r = 0, however, we expect that only the firstfew terms in Eq. (A25) will significantly contribute tothe solution in this case. As we shall see later, this allowsone to approximate the infinite series in Eq. (A25) by acut-off one. Unfortunately, there exists no independentway of determining the coefficients q n because the saddleregion is devoid of any real sources. Thus their valuesare set by the boundary conditions in the intermediateMOND domain and need to be estimated from the fullnumerical result, which is further discussed in the nextsection.
3. Comparison to numerical results
Having derived semi-analytic solutions in both thequasi-Newtonian and the deep-MOND domain, we nowwant to compare these to the full numerical two-bodyresults of the Earth-Sun SP. To do so, however, we still
FIG. 18. Potential derivative S yy of the Earth-Sun SP eval-uated at y ≈ . z = 0 in units of km: Shown are thenumerical result for ν ( y ) given by Eq. (12) with k = 0 . need to determine the series coefficients in Eq. (A6) for aparticular choice of ν ( y ). Let us consider the expressionin Eq. (12) which is the QMOND analog of the TeVeSinterpolating function used in Ref. [18]. For y (cid:29)
1, theabove formula can be expanded as ν ( y ) = 1 + 14 (cid:18) πk (cid:19) y + O (cid:18) y (cid:19) , (A26)which reveals that the first-order term vanishes, i.e. onehas α = 0. Since the numerical approach uses differentboundary conditions (see Sec. IV B), we choose not tocompare the actual potentials, but rather their secondderivatives. In accordance with our numerical setup, weswitch to Cartesian coordinates, letting the symmetryaxis now point along the x -direction.As for a comparison to the deep-MOND solution, wemust also take care of the homogeneous contribution Φ which is determined by the potential values in the in-termediate MOND domain. The coefficients q n can becalculated by fitting the simulation results with the helpof Eq. (A24), which, of course, requires introducing acut-off in Eq. (A25). Since only the first few termssignificantly contribute in the close vicinity of the SP,we cut the series off at n = 3 and discard all remain-ing terms. Furthermore, to ensure that Eq. (A18) isapproximately satisfied, we consider only data pointswith 0 . < r/r < . k = 0 .
03 in Eq. (12) and using the corresponding6
FIG. 19. Same as Fig. 18, but now for the function ν ( y ) givenby Eq. (11). simulation data of the Earth-Sun SP, we find q ≈ . × − ,q ≈ − . × − ,q ≈ . × − . (A27)Fixing y ≈ . z = 0 in units of km, the numeri-cal two-body result for S yy ≡ ∂ y Φ QM of the Earth-SunSP is illustrated in Fig. 18. Since the numerical resultis essentially symmetric with respect to x = 0, we onlyconsider the case x >
0. As can be seen from the figures,the approximate analytic solutions are well recovered bythe full numerical result. Note that the increasing devia-tion from the quasi-Newtonian solution near the bound-ary at x = 10 km is entirely due to the different choiceof boundary conditions (see Sec. IV). The numericalsolution becomes unreliable for x (cid:46) ν ( y ) = 1 + 1 y − y + O (cid:18) y (cid:19) . (A28)From the above, one immediately identifies α = 1 and α = −
1. As is shown in Fig. 19 for S yy , the quasi-Newtonian solution provides an excellent approximationto the full numerical result in this case. On the otherhand, the deep-MOND limit is not reached at the givenresolution (∆ x max ≈ . Eq. (11)Eq. (12)Eq. (13) V c , obs NGC 1560
GasStarsStarsGasStarsGas
NGC 2403NGC 4157 V c [ k m / s ] r [kpc] FIG. 20. Rotation curves of NGC 1560 (upper panel),NGC2403 (middle panel) and NGC 4157 (lower panel): hileThe data points correspond to observations. The dashed,solid and dotted lines represent rough “fits” obtained usingEq. (12) , Eq. (11) and Eq. (13) and setting
M/L = 1(for NGC 1560 and NGC 2403) and
M/L = 1 . Appendix B: Impact of ν ( y ) on galactic scales Any functional form of ν ( y ) chosen for the SS is re-quired to be consistent with the vast amount of availableextragalactic data. The properties of Eqs. (11) and (12)relevant to these scales are well-known and have been ex-tensively discussed in the literature (see, e.g., Ref. [45]).While both prescriptions are able to describe the ob-served rotation curves of spiral galaxies, it is usually con-cluded that the data prefers the simple choice Eq. (11)(or further refined versions thereof) since it provides abetter fit to observations if one leaves the stellar mass-to-light ratio M/L as a free parameter. For Eq. (13),however, the situation is unclear. To get an idea about itsperformance, we consider the rotation curves of three testgalaxies covering different gravitational regimes: NGC1560 (low surface brightness), NGC 2403 and NGC 4157(high surface brightness). For the different choices of ν ( y ), Fig. 20 illustrates the predicted MOND rotationcurves together with the observed ones, assuming con-stant M/L ratios in the B-band. Note that since we areonly interested in the performance of Eq. (13) comparedto the others, we have not tried to actually fit the ob-served curves. Instead, to obtain a very rough match tothe data, we have used
M/L = 1 for both NGC 1560 andNGC 2403 and
M/L = 1 . g s / a ˆ µ FIG. 21. Formally constructed TeVeS interpolating functionsˆ µ for k = 0 .
03 corresponding to the ν -functions given by Eq.(12) (solid line), Eq. (11) (dashed line) and Eq. (13) (dottedline). from the figure, the behavior of Eq. (13) is very similarto that of the simple ν -function, despite their differenceson SS scales. This suggests that the interpolating func-tion Eq. (13) is not only compatible with extragalacticdata, but also recovers the successful descriptive powerof Eq. (11). Appendix C: Interpolating functions for TeVeS
So far we have only considered the situation of low-acceleration bubbles within the QMOND frameworkwhich emerges in the nonrelativistic limit of particularbimetric theories (see Sec. II). Another possibility ofrecovering MONDian dynamics in this limit is realizedin TeVeS. In the following, we want to address whetherthere exist viable interpolation functions for the TeVeSscalar field φ which give rise to similar results as foundin Sec. IV.Applying the approximations for weak fields and quasi-static systems, the gradient of the total gravitational po-tential in TeVeS approximately takes the form [14, 46] ∇ Φ = ∇ Φ N + ∇ φ, (C1)where Φ N is the Newtonian potential, obeying the usualPoisson equation, and the scalar field φ is a solution of ∇ · (ˆ µ ∇ φ ) = kGρ. (C2)Here ρ denotes the density of real matter and ˆ µ dependson the scalar field strength g s ≡ | ∇ φ | [51]. In situations where the curl field vanishes, it follows from Eq. (C2)that ∇ φ = k π ˆ µ ∇ Φ N (C3)which, if combined with Eq. (C1), leads to ∇ Φ = (cid:18) k π ˆ µ (cid:19) ∇ Φ N . (C4)Comparing the above with Eq. (1), we identify theQMOND interpolating function ν as ν = 1 + k π ˆ µ . (C5)Together with our previous assumption that ν ap-proaches unity for strong gravitational fields, the aboveleads to ˆ µ → ∞ in this limit. Such interpolation func-tions ˆ µ have been found to be problematic and may there-fore be regarded as inapplicable [52]. Avoiding this issueby imposing ˆ µ → g s /a (cid:29)
1, we obviously musthave that ν → k/ (4 π ) in the Newtonian limit, consis-tent with Eq. (C5) and resulting in an effective rescalingof the gravitational constant. Since we are here interestedin the purely non-Newtonian part of the gravitationalpotential, this can be easily achieved by formally addingthe constant k/ (4 π ) to a given expression for ν ( y ). Asa consequence, the full solutions of Eq. (1) for the suchconstructed ν -functions only differ by the term k Φ N / (4 π )from their counterparts, which leads to the same resultsas in Sec. IV after subtracting the total Newtonian con-tribution. Another way of saying this is the following:Adding a constant to ν ( y ) does not alter the effectivedensity field if the region of interest is devoid of any realphysical matter sources. Choosing the same boundaryconditions for the potential then yields identical results.Keeping this technical detail in mind, we are now readyto construct the corresponding functions ˆ µ with the helpof Eqs. (C3) and (C5). Starting with the interpolationfunction Eq. (12), a bit of algebra reveals thatˆ µ (cid:112) − ˆ µ = k πa g s , (C6)which is exactly the implicit definition of ˆ µ used in Refs.[16, 18]. Similarly, one may find expressions for the in-terpolating functions given by Eqs. (11) and (13). As-suming k = 0 .