MONDian three-body predictions for LISA Pathfinder
aa r X i v : . [ a s t r o - ph . C O ] A ug MONDian three-body predictions for LISA Pathfinder
Neil Bevis, ∗ Jo˜ao Magueijo, † Christian Trenkel, ‡ and Steve Kemble Theoretical Physics, Blackett Laboratory, Imperial College, London, SW7 2BZ, United Kingdom Astrium Ltd, Gunnels Wood Road, Stevenage SG1 2AS, United Kingdom (Dated: October 30, 2018)In previous work it was shown that MOND theories predict anomalously strong tidal stresses nearthe saddle points of the Newtonian gravitational potential. An analytical examination of the saddlebetween two bodies revealed a linear and a non-linear solution, valid for the outer and inner regions.Here we present a numerical algorithm for solving the MOND equations. We check the code againstthe two-body analytical solutions and explore the region transitioning between them. We thendevelop a a realistic model for the MONDian effects on the saddles of the Sun-Earth-Moon system(including further sources is straightforward). For the Sun-Earth saddle we find that the two-bodyresults are almost unchanged, with corrections increasing from full to new Moon. In contrast, theMoon saddle is an intrinsically three-body problem, but we numerically find a recipe for adapting thetwo-body solution to this case, by means of a suitable re-scaling and axis re-orientation. We explorepossible experimental scenarios for LISA Pathfinder, and the prospect of a visit to the saddle(s)at the end of the mission. Given the chaotic nature of the orbits, awareness of the full range ofthe possibilities is crucial for a realistic prediction. We conclude that even with very conservativeassumptions on the impact parameter, the accelerometers are abundantly sensitive to vindicate orrule out the theory.
I. INTRODUCTION
Modified Newtonian dynamics (MOND [1]) is a schemethat was first proposed for explaining observed dynami-cal properties of galaxies without invoking dark matter.The scheme was later incorporated into a proper the-ory with a Lagrangian formulation [2], but valid only inthe non-relativistic regime. Still later a fully covariantgravitational theory was found containing MOND phe-nomenology as a non-relativistic limit. This theory wasnamed T e V e S [3] and alternatives have been proposed(e.g. [4–6]). The observational features to be studied inthis paper depend only on their (shared) non-relativisticlimit. Constraints arising from lensing [7–9], or cosmol-ogy [10] have no bearing here and indicate issues arisingfrom the relativistic extension of these theories.The opposition between MOND and dark matterleaves considerable doubts as to how to interpret newastrophysical and cosmological data. A fair comparisonrequires re-evaluating, within each approach, the wholeset of assumptions underlying the new observations (seefor example the controversy surrounding the bullet clus-ter [11–14]). For this reason the debate would benefitfrom a direct probe, in the form of a laboratory or So-lar System experiment. Such a perspective motivateswidespread dark matter searches. The analogous “back-yard” tests of MOND include searching for anomaliesin planetary and spacecraft trajectories [15–17], strongertidal stresses in the vicinity of saddle points of the New-tonian potential [18] or Solar System manifestations of ∗ Electronic address: [email protected] † Electronic address: [email protected] ‡ Electronic address: [email protected] the MOND external field effect [19].In this paper we focus on the MONDian effects on thesaddle points of the gravitational field, with particularemphasis on the region where the Earth and Sun pullscancel (which, we stress, is not at the Lagrange pointL1). We provide realistic predictions for the future LISAPathfinder spacecraft [20], which could plausibly be redi-rected to the Earth-Sun saddle, once its primary goals atL1 have been completed. We also consider the possiblebenefits of re-direction to the saddle point near the Moon(there are not two separate Moon-Earth and Moon-Sunsaddle points) and discuss its related practical issues.As already mentioned, the predictions we calculatestem from the non-relativistic limit of T e V e S [3]. Forthis theory gravity is described by the total potentialΦ = Φ N + φ , where Φ N is the Newtonian potential and φ is an additional MOND component. The latter is phys-ically relevant only when |∇ Φ N | . a , where a is theMilgrom acceleration [1], with a ≈ − ms − in orderto fit galaxy observations without dark matter. The ex-tra field is governed by the non-linear Poisson equation: ∇ · [ µ ( κ |∇ φ | /a ) ∇ φ ] = κGρ (1)where ρ is the matter density, κ is a constant parameterand G is the underlying gravitational constant. On theleft-hand-side, µ ( y ) is an unknown function that musttend to 1 for y ≫ y for y ≪
1. Westress that this function is not the ratio of the Newto-nian and actual accelerations, which we will denote as˜ µ ( x ) with x = a/a , a function more often employedby the astronomy community. Although for sphericallysymmetric cases the conversion between ˜ µ ( x ) and µ ( y ) isstraightforward, in more general situations complicationsmay arise (see [10]). We will choose a particular form for µ ( y ) following Ref. [18] (see Eq. 5 in Section III of thepresent paper), and likewise take κ = 0 .
03 throughoutthis article. We’ll comment on other choices of µ andtheir implications for our results, in the concluding Sec-tion of this paper.In regions where |∇ Φ N | ≫ a , we have µ ≈ φ yields accelerations that are κ/ π times theNewtonian contribution. Hence by measuring the to-tal gravitational force we would measure the Newtoniangravitational constant to be (1 + κ/ π ) G and not G .However, given that here we take the value κ = 0 . G isthe standard value measured by experiment.Analytical results exist for the above equation in thecase of two point masses [18], But these are valid only incertain regimes. The “deep-MOND” solution of Ref. [18]is valid close to saddle, while the “quasi-Newtonian” solu-tion is valid for large distances. By tackling the problemusing numerical techniques we can present results accu-rate even in the gap between these two regimes. Thisis particularly important since this gap is in fact the re-gion of most practical importance for probes like LISAPathfinder.Furthermore, since the Earth-Sun saddle point lieswithin the orbit of the Moon, we must concern ourselveswith the effect of the Moon on the LISA Pathfinder mea-surements, and hence consider the more complex three-body problem. Our numerical approach enables us to dojust that, and also to explore the possible advantages ofprobing the hitherto unexplored lunar saddle point.This article is laid out as follows. In Section II wepresent the algorithm employed in our numerical code.Then in Section III we test the code by comparing its re-sults with the analytical solutions previously developedfor the two-body problem in [18]. Aware of the virtuesand limitations of the code in Section IV we explore thesolution for the full three-body problem, in a numberof configurations, and considering both the Earth-Sunand Moon-Sun saddles. Specific recommendations forthe LISA Pathfinder mission are made in Section V andmore general considerations are included in the last Sec-tion. In Appendix A we present details of our numericalalgorithm. II. NUMERICAL METHOD FOR SOLVING THENON-RELATIVISTIC EQUATIONS
In this section we present an overview of our numeri-cal approach for the solution of Eq. 1 (see Appendix Afor full details). Most importantly this involves solvingnumerically only for the region immediately surroundingthe saddle point, as shown in Fig. 1. This region containsno gravitational sources, so that Eq. (1) becomes homo-geneous . The gravitating sources make their presence This of course neglects the very small density present in inter-planetary space and any effects of the LISA Pathfinder measure-
FIG. 1: An illustration of the lattice location when calculatingfor the Earth-Sun saddle region. The calculation is performedfor the a cuboid surrounding the saddle region, but not en-closing any of the gravitating bodies. Also shown is the typeof non-uniform coordinates that we consider in order to boostthe resolution in the central regions of the lattice, i.e. nearthe saddle point. For example, the y -resolution is a functionof the y coordinate only. Not to scale. felt by the boundary conditions on our gridded volume.So long as the boundary is sufficiently far from the saddlepoint, so that 1 − µ ≪
1, then to a very good approx-imation, φ on the boundary is merely φ = κ π Φ N , i.e.a rescaled version of the Newtonian potential. We alsouse the rescaled Newtonian solution inside the box as theinitial configuration supplied to our numerical relaxationalgorithm.It turns out that solving Eq. (1) numerically in termsof φ is sub-optimal since it is only ∇ φ that appears inthis equation. However, solving in terms of either g = −∇ φ or u = − (4 πµ/κ ) ∇ φ results in there being a furtherequation added to the system in order to specify the curl-free nature of g or u /µ . Here we will use g , and so obtain: ∇ · µ g = 0 (2) ∇ × g = 0 , (3)(written in terms of u , the µ appears instead in the sec-ond equation). Fortunately we can bypass the extra workinvolved in dealing with the extra equation by choosing arelaxation algorithm for Eq. (2) that preserves the curl of g as it proceeds. Since we use a Newtonian initial config-uration, with zero curl, we obtain a solution to Eq. (3) atno additional computational cost. We represent the field g on a non-uniform lattice of discrete points, as shownin Fig. 1. Therefore, as we explain in Appendix A, inreality we obtain a solution to a finite difference versionof Eq. (2). The algorithm preserves the finite differenceversion of ∇ × g exactly (i.e. to floating-point accuracyin real life). ment probe. Interestingly, this approach is the reverse of that usedby Milgrom [21], who solved for u using the Newtoniansolution to obtain ∇ · u = 0 initially. The field u wasthen relaxed to a solution in which u /µ has zero curlwhilst not disturbing the divergence. However Milgromconsidered only 2D problems, for which there is just onecurl equation. In 3D there are three curl equations andit is beneficial to instead obtain the solution to the curlequations for free; hence our chosen algorithm. III. THE TWO-BODY PROBLEM
We first consider a two-body setup, such as the Earth-Sun system in the absence of the Moon or any other per-turbations. As we shall show later, the Moon only per-turbs slightly the results for the Earth-Sun saddle pointand hence this simpler problem provides great insightfor the three-body situation to be explored in the nextsection. Furthermore there are analytical solutions forthe two-body problem, so it presents a natural arena forchecking our numerical algorithm. Conversely our nu-merical results will fix a poorly constrained parameterwhich appears in one of the analytical solutions.In Ref. [18] solutions in two different regimes werefound for the MONDian two-body problem. The firstis valid for regions far from the saddle point and waslabelled the “quasi-Newtonian” solution; while the sec-ond is valid closer to the saddle and was dubbed the“deep-MOND” regime. Unfortunately from a practicalperspective, it would be difficult to fly a spacecraft veryclose to the saddle point, where the deep-MOND solu-tion is accurate, due to navigational difficulties but alsobecause, for example, the Earth-Sun saddle point movesas the Moon orbits (see next section) and is perturbedto a lesser extent by other solar system objects and bythe Galaxy. In contrast, for large distances from the sad-dle, where the quasi-Newtonian solution is accurate, theMONDian effects are simply too small. Hence it is in factintermediate distances, between the two regimes, whichare of most practical interest and for which our numericalresults have most value.In this section, unless otherwise specified, we will fol-low Ref. [18] and approximate the Newtonian accelera-tion (providing the MOND boundary conditions) by: g N = A (cid:18) x i − y j − z k (cid:19) . (4)where A is a constant (related to masses of the Sun andthe Earth and their separation; see [18]), the problem issymmetric about the x -axis (connecting Sun and Earth)and the saddle is at (0,0,0). This is a valid approximationin the region near the saddle point and is used in derivingthe analytical solutions. We shall use it for our numericalwork in this Section (but not in the one that follows) inorder to perform a fair comparison with the analyticalresults. We take the function µ ( κ |∇ φ | /a ) to have the formimplied by: µ p − µ = κ πa |∇ φ | . (5)This is merely for illustration purposes: it was the formchosen in Ref. [18] since it eases analytical progress (withour choice of κ = 0 .
03 also following that article). Forthe idealized case of spherically symmetric matter distri-butions this µ can be converted into the ratio of New-tonian to (total) MOND accelerations ˜ µ , which is oftendiscussed in galactic astronomy (see eg. Ref. [9, 10, 22–24]). However, galactic observations are of limited use inguiding our choice for µ since large uncertainties exist inthe matter distribution within galaxies, which in practicewill not be spherical, and a large extrapolation in accel-eration scale is required to obtain predictions for ˜ µ in theregime of interest here. The consideration of alternative µ functions and the conversion of a null result from apass of LISA Pathfinder through a saddle region into aconstraint on µ is left for future work [25]. A. Analytical solutions
The analytical solutions are most easily expressed interms of a dimensionless version of u = − (4 πµ/κ ) ∇ φ ,given by: U = κ π a u = − κ πa µ ∇ φ. (6)This may be also related to the acceleration via: g = −∇ φ = 4 πa κ (1 + U ) / U U / , (7)for the chosen µ form. The field U must be divergencefree in vacuum in order to satisfy Eq. (1), while in orderto keep φ curl-free for this µ function, U must satisfy:4(1 + U ) U ∇ × U + U × ∇ U = 0 . (8)The deep-MOND regime is characterized by U ≪ U ) term can be set to 1, whereas the quasi-Newtonian limit is characterized by U ≫
1. A transitionregion separates the two, located near the ellipsoid: x + y + z r , (9)where r is the key length-scale of problem, given by: r = 16 π a κ A . (10)We have that r ≈
381 km for the Earth-Sun saddle point(ignoring all other gravitating bodies).
1. Quasi-Newtonian solution
In the quasi-Newtonian regime, well outside ellipsoid(9), the solution can be conveniently decomposed as: U = rr N ( ψ ) + r r B ( ψ ) (11)where ( r , ψ , θ ) is a spherical coordinate system centred onthe saddle point with the two gravitating bodies locatedat ψ = 0 and ψ = π . The first term is just the Newto-nian acceleration multiplied by κ / π a , and the sec-ond “magnetic” term has finite curl and is sub-dominant.The angular profiles are given by: N ( ψ ) ≡ N r e r + N ψ e ψ (12) N r = [1 + 3 cos(2 ψ )] (13) N ψ = − sin(2 ψ ) . (14)for the Newtonian component, and by B ( ψ ) = B r ( ψ ) e r + B ψ ( ψ ) e ψ (15) B r = 25 + 3 cos 2 ψ + π √ , (16) B ψ = tan − ( √ − ψ ) + tan − ( √ ψ ) √ ψ − π √ ψ + 1sin ψ . (17)for the magnetic component (these formulae are all de-rived in [18]). Note that it is only in this regime that theNewtonian-like and curl components can be split in thismanner.
2. Deep-MOND solution
Inside the ellipsoid (9) the curl contribution becomesessential, and cannot be disentangled from the full field.A semi-analytical solution was found in Ref. [18] to be: U = C (cid:16) rr (cid:17) α − ( F ( ψ ) e r + G ( ψ ) e ψ ) (18)with α ≈ . F ( ψ ) = 0 . . ψ ) + 0 . ψ ) + . . . ,G ( ψ ) = − . ψ ) − . ψ ) + . . . , (19)which is almost identical to the Newtonian profile ( F ≈ N r and G ≈ N ψ ). These formulae were derived in [18].The normalization C is set by the boundary conditionsand matching the two solutions suggests C ≈
1. Herewe will determine C by comparison with our numericalresults. B. Numerical results and comparison to analyticalsolutions
We have first applied our numerical algorithm to thecase where the Newtonian acceleration on the latticeboundary obeys the linear approximation to the two-body problem given by Eq. (4), exactly as in the aboveanalytical results. This will allow a fair comparison be-tween the two approaches. Considering the Earth-Sunsaddle ( r ≈
381 km) and using a 257 lattice of phys-ical extent 10 000 km and central resolution ≈ . g as illustrated in Fig. 2. These canbe seen to yield a good match to the two analytical so-lutions within their respective domains and provides theappropriate interpolation between them.Due to discretization asymmetries, the position of g = is not precisely at the central site of the lattice. Thelowest value of | g | is instead found on a nearby site, butthis is just a few kilometres away and represents goodaccuracy considering the 10 000 km box size. Beforeperforming the comparison against the analytical resultswe first translate our numerical solution to take this offsetinto account, which can be important at very small r values. Without this, the gap between the thick solid anddashed lines in Fig. 2 would have been more noticeable.The numerical results enable us to determine the pre-viously poorly constrained C parameter appearing in thedeep-MOND analytical solution. To determine this wefirst converted g into U via Eq. (7) and then deter-mined the ratio of the numerical results to the C = 1analytical values for all lattice sites within bounds of r/r = 0 . → .
5, chosen to be comfortably in the deep-MOND regime. We find a ratio of C = 0 . ± .
016 ateach site, the central value from which we have used forall comparisons against the deep-MOND solution.However, g is not the key measurable quantity. We areinstead interested in the tidal stresses, to which LISAPathfinder is sensitive. When calculating the observ-able stresses we must subtract from g the unobservablerescaled Newtonian contribution (see discussion on therescaling of G in the Introduction). Hence we introducethe notation: S ij = − ∂ φ∂x i ∂x j + κ π ∂ Φ N ∂x i ∂x j , (20)for the observable MOND stress tensor, where x i = x , y ,or z . Under the linear approximation to the Newtonianacceleration field (Eq. 4), the Newtonian stress tensoris just S Nij = A diag(1 , − / , − /
2) and hence we mustsimply subtract a constant tensor from the raw MONDresults.We illustrate the nature of this subtraction in Fig. 3and Fig. 4, where it can be seen that the constant un-observable contribution dominates the full stress in thequasi-Newtonian regime (for diagonal elements of thetensor). As a result of this any inaccuracies in our nu-merical results become more important relative to the −2 −1 −2 −1 r / r κ | g r | / a numerical ψ =0numerical ψ = π analytical (deep−MOND)analytical (quasi−Newtonian)10 −2 −1 −2 −1 r / r κ | g ψ | / a numerical ψ = π /4numerical ψ =−3 π /4analytical (deep−MOND)analytical (quasi−Newtonian)10 −2 −1 −2 −1 r / r κ | g r | / a numerical ψ = π /2numerical ψ =− π /2analytical (deep−MOND)analytical (quasi−Newtonian) FIG. 2: A comparison between the numerical and analyticalresults for components of g = −∇ φ in the two-body Earth-Sun case. Results are plotted as function of r for three pairsof ψ values: 0 and π (top), π/ π/ ± π (bottom), with symmetry relating the g -component valueswithin in each pair up to the sign in the analytical case. Inthe numerical case this is also approximately true, except forvery low r where the discretization asymmetry prevents it(this serves as an estimate of the discretization errors). Notethat g ψ is zero analytically for ψ = 0 and π/ r , andthat g r tends to zero for ψ = π/ r limit. Wehave used C = 0 .
839 for the deep-MOND analytical resultsin this figure. −2−1.5−1−0.500.511.52 x 10 −13 |x| / km ∂ g i / ∂ x i / s − re−scaled Newtoniananalytical quasi−Newtoniananalytical deep−MONDnumerical (x>0)numerical (x<0) FIG. 3: The full contribution to the two-body tidal stressesfrom the MOND field and the re-scaled Newtonian contribu-tion that must be subtracted from them in order to yield theobservable component. The positive results are for ∂g x /∂x while the negative ones are for ∂g y /∂y . Results are shownfor the line y = z = 0, for which the analytical results shownare symmetric to the interchange x → − x . This symmetryis not quite realized in numerical case, as can be seen by thesmall differences between the dashed ( x >
0) and solid lines( x <
0) for the numerical solution. The numerical stressesbecome unreliable for | x | <
20 km. size of the signal, once the subtraction occurs. Further-more, when calculating the stress there are finite differ-encing errors, which become increasingly important atsmall | x | . This is evidenced by the differences betweenthe dashed ( x <
0) and solid ( x >
0) lines for the nu-merical results, and their difference relative to the deep-MOND analytical result. Hence neither regime is trivialcomputationally. However, as already noted, for LISAPathfinder we are fortunately interested in the interme-diate regime where the reliability of results is greatest.When a real spacecraft performs a saddle fly-by, it willnot pass precisely through the saddle point. We dis-cuss the likely impact parameter and trajectory for LISAPathfinder in Sec. V, but for illustrative purposes we con-sider here the form of S yy and the numerical uncertaintiesin it for a trajectory along the line y = 100 km. While thenormal Newtonian signal is much larger than that from φ , the former is simply constant (though it would slowlydrift without the linear approximation used here). TheMOND signal, on the other hand, provides a distinctivevariation as the probe passes by the saddle. The form isstable against numerical errors, although minor asymme-tries can be seen near the point of closest passage. Theseare slightly more significant when the lattice spacing isincreased by enlarging the box side to 20 000 km (fromthe 10 000 km used thus far). However, the larger boxprovides greater accuracy in the less important tails of −15 −14 −13 |x| / km S yy / s − analytical quasi−Newtoniananalytical deep−MONDnumerical (x>0)numerical (x<0) FIG. 4: The observable contribution to the ∂g y /∂y tidal stressfrom the MOND field for line y = z = 0. It can be seen thatthe the numerical results perform well even after the subtrac-tion of the unobservable re-scaled Newtonian contribution,which heavily dominates the stresses in the quasi-Newtonianregime. However, near the edge of the simulation volume thestress is attenuated due to the smooth transition of φ to meetthe rescaled-Newtonian boundary conditions. the signal (see Fig. 5).To close this Section we present results using the fulltwo-body Newtonian acceleration field for the Earth-Sunsystem, rather than the linear approximation (Eq. (4)).As can be seen in Fig. 6, the effect of the approximationis insignificant. Fig. 6 also shows the change in the formof the signal for three different impact parameters: y =25, 100 and 400 km. We see that the signal broadenswith increasing y , but more importantly, its amplitudedecreases. IV. THE THREE-BODY PROBLEM
In this section we add a third body to the problem andconsider the perturbing effect of the Moon on the Earth-Sun saddle. We also present results for the lunar saddlepoint (which is intrinsically a three-body problem, as weshall see).We must first comment that the “realistic Solar sys-tem”, as described in Ref. [18], is in fact not realistic atall, and numerical work is needed to provide even quali-tatively correct directions for a space mission. The argu-ment in Ref. [18] regarding the Earth-Moon saddle, forexample, is entirely flawed because of its use of the lin-ear approximation for the Newtonian field for the Earth-Moon system while ignoring the gravitational effects ofSun (cf. Eq. (72) in that paper). As it happens, theperturbation induced by the Sun is too large for this ap-proximation to be valid and in fact there is no true Earth-Moon saddle point. Instead as explained in Fig. 7, the −8000 −6000 −4000 −2000 0 2000 4000 6000 800010 −16 −15 −14 −13 −12 −11 −10 x / km S yy / s − numerical (10 000 km)numerical (20 000 km)analyticalStandard Newtonian −500 0 50011.52 x 10 −14 x / km S yy / s − FIG. 5: The two-body MOND stress signal S yy along the line y = 100 km for two simulation box sizes (10 000 km and20 000 km) with fixed lattice size 257 . The zoomed regionwith linear y -axis shows the minor differences between thesetwo sets of results near the point on this trajectory that isclosest to the saddle, in addition to minor asymmetries whichare due to discretization errors and more significant for thelarger box. The log-scale graph shows the tails of the signal,which are improved in the larger box. Sun-Earth-Moon system presents only two saddle points:the Sun-Earth saddle, which is more or less as the two-body analysis suggests, and a single saddle point close tothe Moon, which is intrinsically a three-body problem.But even for the latter, we shall empirically find a recipefor adapting the two-body results to the real situation.
A. The Sun-Earth saddle
For the Sun-Earth saddle, the effect of the third bodyis not drastic, because the lunar mass is approximately1 / r , is ap-proximately 381 km (see Eq. 9), so the saddle has to belocated to a precision better than this. The saddle pointis approximately 258 800 km from Earth (around two-thirds to the lunar orbital radius), and therefore even asmall perturbation is potentially important. As shownin Fig. 8, at full Moon the saddle shifts about 250 kmtowards the Sun. The magnitude of the shift remainsat approximately this level until the crescent phase ap-proaches, when the Moon begins to approach the saddle.The saddle then starts to experience a large perturba-tion and, in the 3 day period surrounding new Moon,it quickly moves over the left half of the quasi-ellipse −1500 −1000 −500 0 500 1000 1500 2468 −1500 −1000 −500 0 500 1000 1500 12 S yy / − s − −1500 −1000 −500 0 500 1000 1500−0.4−0.2 0 0.2 0.4 x / km FIG. 6: The MOND stress signal S yy along the lines y =25, 100 and 400 km (top to bottom) for the full two-bodyEarth-Sun case (thick black dashed) and when using the linearapproximation to the Newtonian acceleration (black solid). Inthe upper panel the rescaled Newtonian stress is shown (grey)for the y = 25 km case. indicated in the figure. Note that the effect near thenew Moon is significantly enhanced if the Moon is atits perigee during this period. While the perturbation isalways less than 10 000 km, and is less than 1 000 kmfor most of the lunar cycle, this cannot be ignored if aspacecraft is to be navigated through the saddle region.While the perturbation caused by the Moon generallybreaks the axial symmetry present in the two-body dis-cussion of the previous chapter, during a New or FullMoon the symmetry remains unbroken. Hence we mayconsider these two extreme cases via the same linear ap-proximation to the Newtonian acceleration, but with dif-fering A values (cf. Eq. (4)). We can then make thefollowing scaling argument. The conditions on the field U are that it must be divergence-free, satisfy Eq. (8) andmatch the boundary conditions. The first two of theseare unaffected by a rescaling of the spatial coordinates r → α r while the boundary conditions are such that: g → πa κ (cid:18) xr i − yr j − zr k (cid:19) (21)in the linear approximation of Eq. (4). Hence in scaledcoordinates r /r there is a single solution for U valid forall A . Therefore while the form of S ij stresses is unaf-fected by the change in r , the spatial extent of the signalis proportional to r and the stress magnitude scales in-versely with r . (This is a general argument to be usedlater in the Moon saddle analysis).At Full Moon the change in r is minor, increasingby about 1 km from 381 km. However at new Moon,the change is more significant, with r decreasing to FIG. 7: Maps of the Newtonian |∇ Φ N | for different con-figurations of 3 point masses, showing differing numbers ofzero solutions. Black regions denote low |∇ Φ N | values whilewhite denotes high values. The upper-left map is for the caseof equal masses arranged in an equilateral formation, whichyields 3 saddle points plus a central maximum in Φ N . Thesecond frame shows the same number of solutions, but is fora configuration much closer to the Sun-Earth-Moon system,albeit that the Earth’s mass is distributed evenly betweenthe lighter two bodies and the heavier body is slightly lighterthan the Sun. Mass is added to the heavier body in the thirdframe, causing the lower saddle and central maximum to ap-proach each other. Adding yet further mass causes them tomeet and disappear, leaving only two saddle points by thetime the solar mass is reached, as in the fourth frame. Giv-ing the lightest two bodies masses corresponding to the Earthand the Moon, then maintains this number of saddle pointsbut, as can be seen in the final frame, the lunar saddle pointis surrounded by only a small region of low acceleration andit is also significantly tilted towards Earth.
327 km . Hence the Newtonian stresses will be slightlyhigher at new Moon than in the two-body case, and veryslightly lower in the Full Moon case. Considering thesechanges at a fixed r/r value, the MOND g will be un-affected, therefore S ij at fixed r/r increases in inverseproportion to r . However, note that at new Moon, the These figures assume the Moon-Earth separation is equal to thesemi-major axis for the lunar orbit; for the lunar perigee thisvalue should be decreased by 22 km and for the apogee it shouldbe increased by 15 km. −8000 −6000 −4000 −2000 0−8000−6000−4000−200002000400060008000 ∆ x / km ∆ y / k m FIG. 8: The perturbation in the position of the Earth-Sunsaddle due to the presence of the Moon when the lunar orbitalradius is fixed at its semi-major axis (thick), with additionplots for apogee (outer) and perigee (inner) to highlight theeffect of orbital ellipticity. Negative ∆ x values denote that thesaddle moves closer to Earth, which would lie to the left in thisfigure, while the Sun would lie to the right. The saddle is onlyslightly perturbed at Full Moon, moving slightly towards theSun, but as new Moon approaches the perturbation becomeslarge and saddle moves quickly, covering the left half of thequasi-ellipse in approximately 3 days, with new Moon beingthe mid-point of this period. The perturbation at new Moonis enhanced if new Moon coincides with the lunar perigee butattenuated if at apogee. linear approximation is not so robust because the gravityfield of the Moon varies on a smaller length scale thanthat of the Earth.To consider an arbitrary Moon position, and withoutthe approximation to the Newtonian acceleration field,we require numerical methods. The results in Fig. 9,which shows three different lunar phases, indicate thatthe effect of the Moon on the MOND signal near theEarth-Sun saddle point is minor. At Full Moon the signalis essentially identical to the two-body case, since theMoon is more than twice as far from the saddle point asEarth. The saddle is displaced furthest away from the −1500 −1000 −500 0 500 1000 1500 2468 −1500 −1000 −500 0 500 1000 1500 12 S yy / − s − −1500 −1000 −500 0 500 1000 1500−0.4−0.2 0 0.2 0.4 x / km FIG. 9: The MOND stress signal S yy along the lines y = 25,100 and 400 km (top to bottom) for the three-body Earth-Suncase for different lunar phases: new Moon (thick, black, solid),Full Moon (thick, black, dashed) and when the Moon appears18 ◦ away from the Sun towards positive y (thin, black, solid).Also shown in the y = 25 km case are the rescaled Newtonianstresses (grey). Earth-Sun line when the angle between the Moon andthe Sun (measured from Earth) is approximately 18 ◦ ,and even then the line in the figure that corresponds tothis asymmetric case still has the same basic form as thetwo-body signal. Some x → − x asymmetry is instilled,but the overall conclusion is that the MOND signal forthe Earth-Sun saddle is largely unaffected by the phaseof the Moon. B. The Lunar saddle
The lunar saddle can potentially provide significantlylarger stress signals than the Earth-Sun saddle and istherefore of potential interest experimentally. However,the larger stress signal implies that the MOND region issmaller and so the challenge of navigating a space probethrough the saddle region becomes more difficult.
1. Newtonian gravity field
We must first consider the lunar saddle from a purelyNewtonian perspective, since the Newtonian physics setsthe saddle location and the boundary conditions for ourMOND calculation. It is useful to think of the lunarsaddle as a heavily perturbed Moon-Sun saddle, sincethe Earth is the least important of the three bodies inthe saddle vicinity.Although the gravitational pull of the Sun is effectivelyconstant on the scale of the lunar orbit, the Moon is ap-proximately 1 / / r for the lunar sad-dle should be around 9 times smaller than that for theMoon-Sun saddle, before third-body perturbations areconsidered.However the perturbing effect of the Earth on the lu-nar saddle is very large, as can be seen in Figures 10 and11. At new Moon the Earth helps the Moon overcomethe pull from the Sun and hence the saddle is furtherfrom the Moon than on average and the lunar stress con-tribution is smaller, yielding a larger saddle region. AtFull Moon, the Earth acts with the Sun and thereforethe saddle is closer to the Moon, the stresses are greaterand the saddle region becomes very small. Additionally,the Earth shifts the lunar saddle significantly away fromthe Moon-Sun line, by up to ≈ ◦ (measured at theMoon), although of course it lies on this line at full andnew Moon (ignoring the possible small displacement ofthe Moon off the Ecliptic). The lunar saddle is thereforeintrinsically a three-body problem.Importantly, since the Moon dominates the tidalstresses near the saddle, the form of the Newtonian accel-eration in the saddle region is still approximately that ofthe linear approximation form, just as in the Earth-Suncase where the Earth dominates the stress. However, nowthe effective “axis of symmetry” points from the Moonto the saddle location. At once this suggests a recipe foradapting the two-body MONDian solution to the lunarsaddle.
2. MOND stress signal
Given the observation that the Newtonian accelerationfield near the lunar saddle is roughly linear (as in Eq. 4),it might be expected that the MOND solution would bevery similar to the two-body case, except for re-scalingswith r . This is indeed seen in the numerical resultsdisplayed in Fig. 12, where the new Moon case is indis-tinguishable from the previous Earth-Sun results and theeven the non-axial quarter Moon case is only slightly dif-ferent, once appropriate scalings with r are introduced .For the purposes of this figure we have introduced a newcoordinate system ( x ′ , y ′ , z ), which is a rotated versionof ( x, y, z ) such that the x ′ -axis joins the Moon and thesaddle point, which still lies at (0,0,0). We have thenchosen y ′ ≈ . r and 1 . r lines for which to plot theMOND stress signals, where r ≈
81 km at new Moonand r ≈
38 km at quarter Moon. The two panels in thisfigure are hence directly comparable with the lower two
10 00020 00030 000 D i s t. f r o m M oon / k m
20 40 60 80 r ff / k m de f e c t i on / deg . θ moon / deg. FIG. 10: Upper panel: The distance of the lunar saddle fromthe Moon as a function of the lunar phase, where θ moon = 0denotes new Moon and θ moon = 180 ◦ denotes Full Moon.Middle panel: The effective r value for the lunar saddle as afunction of lunar phase. Lower panel: The deflection angle ofthe saddle away from the Moon-Sun line (towards Earth).FIG. 11: Four maps of the Newtonian acceleration strength | g N | (as before dark means small | g N | , and so a saddle, in-dicated by a X wherever it makes it clearer) in the vicin-ity of Earth for Sun-Earth-Moon system (in plane z = 0).Clockwise from top left the angle between the Moon and theSun, measured at Earth are 0 ◦ (new Moon), 18 ◦ , 90 ◦ (FirstQuarter) and 180 ◦ (Full). The Earth (0,0) and the Moon areshown by crosses. White corresponds to values greater than80 × a , black indicates values less than 25 × a , withshades of gray indicating intermediate values. Note that theEarth perturbs the Moon-Sun saddle region heavily, both inorientation and size. − (r / a ) S y ’ y ’ Lunar saddle (New)Lunar saddle (Quarter)Earth−Sun saddle (linear approx.)
FIG. 12: Stress signals for the lunar saddle, compared to theresults from the Earth-Sun saddle (in the linear approxima-tion). Lunar results are expressed in a new coordinate system( x ′ , y ′ , z ), a rotated version of ( x, y, z ) such that the x ′ -axisis the line joining the Moon and the saddle, with y ′ ≈ . r (top) and 1 . r (bottom). In the Earth-Sun case, the re-sults are for y = 100 and y = 400 km, which correspond tothe 0 . r and 1 . r , making the figure directly comparablewith the lower two panels of Fig. 6. At New moon r ≈
81 kmwhile at Quarter Moon r ≈
38 km. panels in Fig. 6. However, we have not included resultsequivalent to the upper panel there ( y ′ ≈ . r ) sinceat Quarter Moon this amounts to y ≈ . /r , so for this impact parameterthe likely stress level is largely independent of r . V. TESTING MOND WITH LISA PATHFINDER
We now explain more specifically how the results foundin this paper translate into a concrete signal, given thenavigational constraints and instrument noise propertiesof LISA Pathfinder. We consider one example here andin a companion paper [26] perform a more comprehensivestudy of the issues highlighted.
A. Overview of LISA Pathfinder Project
The LISA Pathfinder (LPF) project [20] is an Euro-pean Space Agency (ESA) mission designed to test thetechnology intended to be employed in LISA (Laser In-terferometer Space Antenna), a proposed gravity waveobservatory. While LISA itself is planned to consist ofthree spacecraft in a triangular configuration of side-length 5 × km, the Pathfinder mission will consistof two test masses within a single spacecraft and just 35 cm apart. These test masses will follow geodetic mo-tion to exceptional precision due to protection from solarradiation pressure and other unwanted non-gravitationalinfluences, while the spacecraft itself is designed to min-imize its own gravitational impact. An interferometerwill study the relative motion of the test masses suchthat LPF will be an excellent instrument for measuringtidal stresses, such as those predicted in the precedingtwo sections. The mission is presently in the “implemen-tation” phase with its launch expected in 2012.The nominal plan for the LPF mission is for it to en-ter a large-amplitude Lissajous orbit (of near halo-orbitdimensions) around the L1 Lagrange point of the Earth-Sun system, at a distance of ∼ . × km from Earth.Once the planned mission at L1 is complete, it can inprinciple be extended, by breaking from the Lissajousorbit and passing close to the Earth-Sun saddle point,or alternatively near the lunar saddle. In this section wecompare the sensitivity of LPF to the possible MONDsignal, given its navigation and measurement capabili-ties. B. Spacecraft navigation
The objective of the LPF mission does not require thespacecraft to have the large thrusters of an interplane-tary probe, since the journey out to L1 is made possibleby a propulsion module from which LPF then separates.Instead LPF will carry only low-thrust micro-propulsionsystems and the ability to navigate the probe near a sad-dle point is not something that can be taken for granted.Fortunately, the Lissajous orbit is unstable and henceLPF can be ejected from the L1 region using a fairlysmall change in its velocity. For example, a 30-day burncan yield a d v of ∼ − and put LPF on a trajec-tory that brings it close to Earth in a reasonable time-frame, of order one year. The precise orbit depends onthe timing of the burn relative to the phase of the Lis-sajous orbit, the burn magnitude, subsequent correctionmanoeuvres and any close lunar fly-bys. Without thelatter two, approaches near the Earth-Sun saddle in ap-proximately 1 yr tend to have quite large impact param-eters ∼
10 000 km but closer approaches < v manoeuvres. Due to Earth’s gravity the spacecraft willbe cruising through the saddles with a speed of the or-der of ∼ − . In what follows we will take it thatLPF can be directed through a saddle region with veloc-ity 1 kms − and an impact parameter ∼
50 km. A morecomplete study of these orbits is presented in [26].1
C. Sensitivity to MOND stress signal
Firstly, it should be noted that the design of LPF lim-its it to being sensitive only to one of the diagonal el-ements of the stress tensor, eg. S yy in our notation.Additionally, in order to yield a stable radiation pres-sure it is desirable to keep fixed the orientation of the(roughly cuboid) probe relative to the Sun, with the sideon which the solar panel is mounted facing directly at theSun. Given the vast distance of the spacecraft from theSun, the solar panel is therefore effectively aligned nor-mal to the x -axis and the spacecraft design then impliesthat S xx cannot be measured. This leaves the measure-able stress as cos ( α ) S yy + sin ( α ) S zz , where α specifiesthe orientation of the test masses relative to the y -axis.Estimates of the sensitivity of LPF as a gradiometerare normally expressed via the power spectral density ofthe stress signal. This is defined as: P ( f ) = 2 T (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)Z + T/ − T/ d t S ij ( t ) e − πift (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (22)where f is the frequency, t is the time and T is the in-tegration period. When integrated over f from 0 to in-finity this yields the mean “power” in the signal, i.e. thetime-average of S ij , with the factor of 2 present becausenegative frequencies are folded in with the positive ones.For the present article we assume that LPF will attain asensitivity in the square-root of this quantity, the ampli-tude spectral density, of 1 . × − s − / √ Hz between 1and 10 mHz [20]. Beyond this range we will take the sen-sitivity to degrade as 1 /f at lower frequencies and as f at higher frequencies. This yields a good approximationto published expectations [20].We simulate gradiometer noise under the assumptionof Gaussianity (independent Fourier mode phases) andthe above spectrum. We add this to the theoretical stresssignal, which is the combination of the standard Newto-nian and MOND components, and then apply a simplecleaning algorithm to extract a likely measured MONDsignal. Specifically, we approximately remove the New-tonian signal by performing a quadratic fit to the data.In the linear approximation, the Newtonian stress sig-nal is just a constant while the gravitational accelerationnear the saddle is so small that we may assume the probemoves with constant velocity and therefore this fit takesinto account two further orders, beyond linear, in thespatial Taylor expansion of the Newtonian acceleration.We then apply a top-hat band-pass filter to remove allfrequencies outside the range 0 . − through the Earth-Sun saddle region, passingon a trajectory with y = 50 km and z = 0, during a fullMoon. While the velocity of the probe is unlikely to bealigned with the x -axis in this manner, the details of thetrajectory orientation do not appear to be greatly impor-tant. A significant signal is clearly seen in both real-time −1 −0.5 0 0.5 1x 10 −3−2.5−2−1.5−1−0.500.511.5 x 10 −14 S yy / s − t / s −1 −0.5 0 0.5 1x 10 −6−5−4−3−2−1012 x 10 −14 S zz / s − t / s FIG. 13: Theoretical and recovered stress signals S yy and S zz for a y = 50 km, ˙ x = 1 kms − trajectory, including a real-ization of gradiometer noise and applying a simple cleaningalgorithm (see text). These results are for the three-bodyEarth-Sun saddle with the lunar phase being full Moon. and Fourier domains: Figs. 13 and 14, respectively. Ifthe actual trajectory saw LPF cross the saddle regionmore quickly, then the MOND signal would be shifted tohigher frequencies where the signal to noise is even bet-ter. An increase in the impact parameter r min reducesthe signal roughly as 1 /r min for r min ∼ r , but as 1 /r for r min ≫ r . VI. CONCLUSIONS
We solved the non-relativistic limit of T e V e S for twoand three-body saddle points using a particular µ func-tion used previously in the literature. We predicted ananomalous stress signal that is potentially detectable byfuture space-probes. We found that past analytical calcu-2 −4 −3 −2 −14 −13 −12 f / Hz AS D s − / H z − / MONDnoise
FIG. 14: The amplitude spectral density (ASD) for the S zz MOND signal for a y = 50 km, ˙ x = 1 kms − trajectory,compared to the simple noise model employed here. Notethat while the noise ASD is independent of the integrationtime, there being a single MOND event means that its ASDlessens with increasing integration time, where for the presentplot we conservatively integrate over 2 × s. lations for the two-body case are not only in agreementwith our numerical results, at least in their regimes ofapplicability, but provide a good insight into three-bodyproblems of the Sun-Earth-Moon type, where the New-tonian stresses at the saddle points are dominated by asingle body. Our results indicate that the most significanteffect of the Moon on the Earth-Sun saddle is to changeits location, while the MOND anomalous stress signal isnot greatly affected, and is stable to changes in the phaseof Moon. Furthermore, we have demonstrated that thesingle lunar saddle point near the Moon yields a stresssignal of similar form but differs in stress magnitude andspatial extent. However, for an impact parameter likelyfrom a spacecraft flyby, we note that these two factorseffectively cancel (at least in some regime).Finally, we have demonstrated that a possible exten-sion to the LISA Pathfinder mission, targeting the sad-dles, would be sensitive to the anomalous stress signalthat we have calculated. This is in part due to the dis-tinct form of the signal close to a saddle point, whichcannot be confused with the (much smoother) Newtoniansignal (or with any other signal derived from the param-eterized post-Newtonian framework). It is also due tothe likely speed of the probe through the saddle regionyielding signal variations at frequencies that minimizethe instrument noise. While this conclusion assumes theability to pass within ∼
100 km of a saddle point, webelieve that this is achievable even with the low thrustpropulsion system incorporated in the LISA Pathfinderdesign. In future work we hope to implement our studiesas part of a concerted effort to explore realistic trajecto-ries for LISA Pathfinder. We conclude with a general remark on MONDian the-ories which we hope will clarify a number of issues withthe proposed test. For all MONDian theories with a rela-tivistic formulation (and not only TeVeS, as described inthe introduction), MONDian non-relativistic effects aredue to an extra field, φ , as opposed to the gravitationalfield itself, i.e. the metric field, with Φ N associated with g . Indeed it’s very hard to covariantly modify a gravitytheory so that in the non-relativistic limit g is ruled bya non-linear Poisson equation. In the Newtonian limit(when µ ≈ φ mimics the New-tonian field and so renormalizes the gravitational con-stant ( G → G (1 + κ/ (4 π )) as discussed in the literature(e.g. [3]). To satisfy constraints, φ must be sub-dominantin this regime, and this is enforced by the parameter κ appearing in its Poisson equation (1).But this simple fact implies that one must triggerMONDian behavior in φ at Newtonian accelerations a N much higher than a . Only thus can φ ’s relative impor-tance start to increase with decreasing a N , so that bythe time a N ∼ a the field φ is not only MONDian butdominates Φ N , as required by astronomical applications.It can be easily computed that if µ turns from 1 to asingle power-law, then MONDian behavior in φ shouldbe triggered when the Newtonian acceleration a N dropsbelow (4 π/κ ) a ∼ . × a ∼ − ms − . And in-deed this is the rough Newtonian acceleration at r . Itis therefore important to realize that for LPF realisticimpact parameters we would probe the regime where φ has gone fully MONDian, but hasn’t yet dominated Φ N .The MOND signal can be detected above the Newtonianone because it has a distinctive spatial variation whereasΦ N is just a DC component.Bearing this in mind a number of issues may be clari-fied. One concerns the self-gravity of LPF. This is onlybalanced at the level of a N ∼ − ms − , so one mightthink that a test of MOND, in the regime a N ∼ a , wouldrun against the wall of self-gravity. In fact we are testingmuch higher Newtonian accelerations; for instance, for animpact parameter of 40 Km we have a N ∼ − ms − .We’d need to approach the saddle much closer than about400 meters before self-gravity became an issue (and thespacecraft itself had to be included in the computationof the location of the saddle).This should also clear the matter of the generality ofthe predictions made. There are several µ on offer forastrophysical purposes. They all have the same roughbehavior in the regime under study, where φ is fullyMONDian but hasn’t yet dominated Φ N . A negativeresult from LPF would rule out virtually all proposed µ .It would only fail to rule out very contrived µ ( y ) (neversuggested in the literature), for which two power laws areemployed in µ ( y ): a very steep one from y = 1 (felt at r )to the point where a N = a ; and the usual linear power-law for a N < a : extremely contrived. We are currentlywriting a follow up paper expanding on this matter,The high sensitivity of LISA Pathfinder enables it topotentially detect the anomalous MOND stress when the3deviations from Newtonian dynamics are tiny, in con-trast to the conditions present in outer reaches of galax-ies [7, 22]. On the other hand, constraints are yieldedby planetary orbits [15] only at much stronger accel-erations than LISA Pathfinder would probe. A saddlepoint fly-by offers the chance to study MOND where wehave an exquisite knowledge of the mass distribution andthe ability to make full calculations of the MOND accel-eration field, without assuming simplifying symmetries.Furthermore, saddle point exploration would enable aninvestigation of gravity at very low accelerations in theclose vicinity of Earth, which can be reached in a shorttime-frame and without the increased mass distributionuncertainties affecting probes sent to very large distancesfrom the Sun. A positive detection would, of course, vin-dicate the MOND paradigm, while a null result wouldprovide clean and reliable constraints upon it.In closing we note that several refinements to ournumerical calculations should be trivial. Including thequadrupole and higher multipole moments of the sourcesor the effect of Jupiter and other Solar system objects (oreven the galactic field) is straightforward, as they only af-fect the boundary conditions for our problem. Howeverthese effects are expected to be small. The shift of thesaddle location is of course non-negligible, but can be eas-ily computed with effects already studied in [18]. Oncethe new location is taken into account the change in thedetailed MOND tidal stresses can be computed just bychanging the boundary conditions of our code. However,if the effect of the Moon on the Earth-Sun predictionsfor realistic impact parameters is already so small, theseeffects can be expected to be completely negligible. Acknowledgments
We acknowledge financial support from STFC (N.B.)and thank J. Bekenstein, and C. Skordis for helpful dis-cussions. The numerical work was performed on theCOSMOS supercomputer (which is supported by STFC,HEFCE and SGI) and the Imperial College HPC facili-ties.
Appendix A: Details of numerical method
Here we briefly describe our numerical algorithm. Thisis based upon representing g on the sites of a non-uniformlattice of the form shown in Fig. 1, such that we obtaingreater spatial resolution near the saddle point. A relax-ation algorithm then cycles over each lattice site x andchanges the values of g at x and the neighboring sitesso that the divergence equation Eq. (2) is solved locally ,that is: D x = X j (cid:16) µ x g j x − µ x − j g j x − j (cid:17) = 0 . (A1) We then move to the next site and change the field val-ues so that the condition is valid there, using the newlyupdated values as we proceed. However, when enforcingthe above condition at these subsequent sites, the valueof D x at the first site will be slightly changed. We there-fore require many cycles over the whole lattice before theabove condition is closely matched globally.We stress that we use a different discretization proce-dure to Ref. [21], in that we define all three componentsof g as well as µ at the same location, that is: µ x = µ ( κ g x /a ) . (A2)This is in contrast to the more elaborate scheme used inRef. [21], in which u j is defined on the j -links betweenthe lattice sites and then µ is defined at the centers ofthe grid squares in the 2D case considered. Under thatscheme the calculation of u i /µ at any position in a 3Dcalculation would require knowledge of 33 values of u j ,whereas to determine µ x g j x here, we require merely theknowledge of the three components of g x .The presence of the non-linear function µ is a no-table complication. Furthermore its form is not uniquelyknown and therefore we desire to be able to solve for anarbitrary function. The algorithm of Ref. [21] involvessolving the curl equation for u locally, while keeping the µ values fixed at their old values. Once the µ valuesare updated using the new u , the curl equation is nolonger matched, and this slows the convergence of thealgorithm. In contrast, here we proceed by solving thedivergence equation to first order in δ g and δµ , where δ denotes the change from one step to the next. Then, asthe system converges to the solution, the terms of order δ become negligible very rapidly.Crucial in the above is to ensure that in updating thefield configuration the curl of g remains zero. We definethe our preserved discrete curl as:( ∇ × g ) k x = g j x + i − g j x r i x + i − r i x − g i x + j − g i x r j x + j − r j x , (A3)where r x is the position vector site x . In order to preservethis discrete curl at each step of the relaxation we changethe fields according to: δg j x = + C x r j x + j − r j x , (A4) δg j x − j = − C x r j x − r j x − j , (A5)where the value of C x is chosen, as below, to yieldEq. (A1) to first order. While the first equation actson all three components of g , the second acts on onlyone component per site.To determined C x , consider that to first order in δ g and δµ : δD x ≈ X j µ x δg j x + g j x δµ x − µ x − j δg j x − j − g j x − j δµ x − j r j x − r j x − j . (A6)4Writing δµ in terms of δg and substituting the abovevalues for δ g then yields: δD x C x ≈ X j " µ x ∆ j − ∆ j + + 2 X i g i x ∆ i + g j x ∆ j − d µ x d g x (A7)+ µ x − j (∆ j − ) + 2 g j x − j ∆ j − ! d µ x − j d g x − j (A8)where we have used compact notation such that:∆ j + = r j x + j − r j x (A9)∆ j − = r j x − r j x − j . (A10)We can then insert the µ derivatives for our chosen µ function. Finally, since we want D x to be zero after eachchange, we set: δD x = − D x , (A11)and hence know the C x required to obtain D x = 0, atleast to first order. (If µ had not changed during this step,the above procedure would have set D x to be exactlyzero, i.e. to all orders in δ g .)Note that in addition to preserving the curl, Eqs. (A4)and (A5) also exactly preserve the change in φ measured across the lattice as:∆ φ = X x j (cid:16) r j x + j − r j x (cid:17) g j x , (A12)where x i and x k are fixed during the summation.In practice, cycling over the lattice and solving thediscrete equation locally does not lead to rapid enoughconvergence to the solution. As D x is (approximately)zeroed at later sites, D x at earlier sites is moved slightlyaway from its desired value and a large number of it-erations of this procedure are required before each g x has converged to a good approximation. We there-fore preempt the changes to the field that will occur atother points in the cycle using the fact that these verychanges are (largely) responsible for D x being non-zero.This is achieved by a method known as successive over-relaxation (SOR, e.g. [27]) in which δg j → λδg j , where λ is the over-relaxation parameter and is larger than unity.We begin with λ = 1 and increase it once the system hasbegun to settle down, since high values of λ can initiallyresult in the RMS value of D x increasing, contrary to ourgoal.We coded the algorithm outlined in this Appendix us-ing the LATfield library [28]. [1] M. Milgrom, Astrophys. J. , 365 (1983).[2] J. Bekenstein and M. Milgrom, Astrophys. J. , 7(1984).[3] J. D. Bekenstein, Phys. Rev. D70 , 083509 (2004);Erratum-ibid.
D71 , 069901 (2005), astro-ph/0403694.[4] R. H. Sanders, Mon. Not. Roy. Astron. Soc. , 459(2005), astro-ph/0502222.[5] T. G. Zlosnik, P. G. Ferreira, and G. D. Starkman, Phys.Rev.
D74 , 044037 (2006), gr-qc/0606039.[6] T. G. Zlosnik, P. G. Ferreira, and G. D. Starkman, Phys.Rev.
D75 , 044017 (2007), astro-ph/0607411.[7] I. Ferreras, M. Sakellariadou, and M. F. Yusaf, Phys.Rev. Lett. , 031302 (2008), 0709.3189.[8] N. E. Mavromatos, M. Sakellariadou, and M. F. Yusaf,Phys. Rev.
D79 , 081301 (2009), 0901.3932.[9] I. Ferreras, N. E. Mavromatos, M. Sakellariadou, andM. F. Yusaf, Phys. Rev.
D80 , 103506 (2009), 0907.1463.[10] C. Skordis, Class. Quant. Grav. , 143001 (2009),0903.3602.[11] D. Clowe et al., Astrophys. J. , L109 (2006), astro-ph/0608407.[12] D.-C. Dai, R. Matsuo, and G. Starkman, Phys. Rev. D78 , 104004 (2008), 0806.4319.[13] G. W. Angus and S. S. McGaugh (2007), 0704.0381.[14] J. R. Brownstein and J. W. Moffat, Mon. Not. Roy. As- tron. Soc. , 29 (2007), astro-ph/0702146.[15] M. Sereno and P. Jetzer, Mon. Not. Roy. Astron. Soc. , 626 (2006), astro-ph/0606197.[16] S. G. Turyshev and V. T. Toth (2009), 0906.0399.[17] J. R. Brownstein and J. W. Moffat, Class. Quant. Grav. , 3427 (2006), gr-qc/0511026.[18] J. Bekenstein and J. Magueijo, Phys. Rev. D73 , 103513(2006), astro-ph/0602266.[19] M. Milgrom (2009), 0906.4817.[20] P. McNamara, S. Vitale, and K. Danzmann (LISA),Class. Quant. Grav. , 114034 (2008).[21] M. Milgrom, Astrophys. J. , 617 (1986).[22] B. Famaey and J. Binney, Mon. Not. Roy. Astron. Soc. , 603 (2005), astro-ph/0506723.[23] C. Nipoti, P. Londrillo, and L. Ciotti (2008), 0811.2878.[24] O. Tiret and F. Combes (2007), astro-ph/0701011.[25] N. Bevis, J. Magueijo, and A. Mozzafari, in preparation(2010).[26] C. Trenkel, S. Kemble, N. Bevis, and J. Magueijo, sub-mitted (2009).[27] Press, William and al. (1992), Numerical Recipes in C,Cambridge University Press, UK.[28] N. Bevis and M. Hindmarsh, URL