CClosed periodic orbits in anomalous gravitation
Bjorn A. Vermeersch ∗ Grenoble, France † (Dated: November 7, 2018)Newton famously showed that a gravitational force inversely proportional to the square of thedistance, F ∼ /r , formally explains Kepler’s three laws of planetary motion. But what happensto the familiar elliptical orbits if the force were to taper off with a different spatial exponent?Here we expand generic textbook treatments by a detailed geometric characterisation of the generalsolution to the equation of motion for a two-body ‘sun/planet’ system under anomalous gravitation F ∼ /r α (1 ≤ α < α , and discuss conceptual connections of the case α = 1 to Modified Newtonian Dynamics and galactic rotation curves. I. INTRODUCTION
Newton’s formulation of universal gravitation, andproof that it gives rise to Kepler’s laws of planetary mo-tion, stands as one of the great triumphs of analyticalthinking. While relativistic corrections become necessaryin certain applications, Newtonian gravity, first published332 years ago, remains firmly in the driver’s seat in spaceexploration, and continues to produce novel insights [1].Underpinning these discoveries is the central notionthat any object attracts any other object with a forceproportional to the product of the two masses and in-versely proportional to the square of their distance: F ∼ M M / (cid:107) (cid:126)r − (cid:126)r (cid:107) . Although the 1 /r spatialdependence had previously been proposed by Hooke [2],Newton was the first to demonstrate mathematically thatthis functional form gives rise to the elliptical planetaryorbits Kepler had painstakingly deduced empirically.The fact that the attractive force falls off quadraticallywith distance makes a lot of intuitive sense. Drawing par-allels with Coulomb’s electrostatic force or light emissionfrom a star, one may argue that ‘gravitons’ seemingly em-anating from the object and exerting the gravitationalpull at distance r get distributed across a surface areathat grows as r in 3D space. The stunning accuracyof celestial event predictions [3] and long list of highlysuccessful space missions [4] underpinned by Newtoniangravity models provide compelling empirical evidence.This, however, does not have to stop us from wonder-ing what would happen to Kepler’s neat ellipses if thegravitational force had a different spatial decay, say F ∼ (cid:107) (cid:126)r − (cid:126)r (cid:107) α ≤ α < What if ? , even if only hypothetically, carries ped-agogical value and often leads to thought-provoking and ∗ Email: [email protected] † The research presented here was carried out at the author’s owninitiative in his personal free time. No employer’s resources, mon-etary or otherwise, were used in its creation. entertaining insights [5]. Discussions of central forceproblems are available in text books such as [6]. A laterstudy [7] examined the generic stability and boundednessof orbits and found trajectories commonly displayed inpopular texts to be physically unfeasable. Both of thesecited works limited the mathematical description of thetrajectories to first-order Taylor perturbations aroundthe circular orbit.Here we expand and concretise these prior treatmentsby discussing the general (full-order) solution of the equa-tion of motion. In particular, we derive for a given spa-tial exponent α the emergence and detailed geometri-cal properties of closed periodic orbits as function of theinitial conditions. In addition, we show that a gravita-tional force proportional to reciprocal distance ( α = 1)has conceptual connections to Modified Newtonian Dy-namics (MOND), a prominent explanation of galactic ro-tation curves without invoking dark matter. II. EQUATION OF MOTIONA. Problem statement
We wish to study the dynamics of a two-body systemcomprised of a primary (mass M ) and satellite (mass M (cid:28) M ) under a gravitational force given by F α = G α M M (cid:107) (cid:126)r − (cid:126)r (cid:107) α , ≤ α < G α is a fractional gravitational constant with unitsm α +1 kg − s − , (cid:126)r , ( t ) denotes the time-dependent posi-tions of the two bodies, and α is a characteristic expo-nent. The limit case α = 2 corresponds to the familiarNewtonian one, for which G (cid:39) . × − m kg − s − .Before exploring the equation of motion in its generalform, we derive the radius R c and associated characteris-tics of the circular orbit that exists irrespective of the α value. These quantities will serve as normalisation valuesin the study of more complex orbits in dimensionless vari-ables. In circular motion, the gravitational force provides a r X i v : . [ phy s i c s . pop - ph ] A ug the centripetal acceleration of the satellite: F α ( (cid:107) (cid:126)r − (cid:126)r (cid:107) = R c ) = M v R c (3)from which we readily obtain the linear and angular ve-locities in terms of the orbital radius: v c = (cid:115) G α M R α − (4) ω c = v c R c = (cid:115) G α M R α +1c (5)Invoking the conservation of angular momentum M R c v c = M R ω c = constant = L (6)yields the orbital radius: R c = (cid:18) L G α M M (cid:19) / (3 − α ) (7)Finally, the orbital period is given by T c = 2 πω c = 2 π (cid:115) R α +1c G α M (8) B. Master equation
Central force fields confine the motion to a plane [7] forwhich we can introduce a polar coordinate frame ( r, ϕ )with the primary at the center r = 0. Applying Newton’ssecond law to the satellite in anomalous gravitation yieldsd r d t − r (cid:18) d ϕ d t (cid:19) = − G α M r α (9)We now invoke conservation of angular momentum M r ( ϕ ) d ϕ d t = M r ( ϕ ) ω ( ϕ ) = constant = L (10)and revert to dimensionless coordinates¯ r ≡ rR c (11)where R c is the radius of the circular orbit with sameangular moment L discussed previously. A subsequentchange of variable u ( ϕ ) ≡ r ( ϕ ) = R c r ( ϕ ) (12)finally produces d u d ϕ + u ( ϕ ) = 1 u ( ϕ ) − α (13) Without loss of generality, we can restrict ourselves toinitial conditions of the form u (0) = u ≥ , u (cid:48) (0) ≡ d u d ϕ (cid:12)(cid:12)(cid:12)(cid:12) ϕ =0 = 0 (14)because, as shown below, the general solution for u is pe-riodic in ϕ . One can therefore always perform a rotationof the coordinate axes that places a maximum of u at ϕ = 0, i.e. the orbit has a perihelion on the x axis. C. Angular and linear velocities
Conservation of angular momentum (10) readily yields ω ( ϕ ) = LM r ( ϕ ) = LM R r ( ϕ ) = ω c u ( ϕ ) (15)Points with extremal angular velocity obeyd ω d ϕ = 0 ↔ u ( ϕ ) d u d ϕ = 0 (16)As the first factor remains strictly positive for boundedorbits, the maximal and minimal angular velocities areachieved in the perihelia and aphelia respectively.The linear velocity of the satellite is given by v ( ϕ ) ≡ (cid:13)(cid:13)(cid:13)(cid:13) d (cid:126)r d t (cid:13)(cid:13)(cid:13)(cid:13) = ω ( ϕ ) (cid:115)(cid:18) d r d ϕ (cid:19) + r ( ϕ ) (17)Using (15) and v c = ω c R c , the linear velocity can beexpressed directly in terms of the inital value problemsolution u ( ϕ ) = R c /r ( ϕ ) as v ( ϕ ) = v c (cid:115)(cid:18) d u d ϕ (cid:19) + u ( ϕ ) (18)The velocity reaches an extremal value whend v d ϕ = 0 ↔ d u d ϕ · (cid:20) d u d ϕ + u ( ϕ ) (cid:21) = 0 (19)Per (13), the factor between brackets equals u α − ( ϕ ),which remains strictly positive for bounded orbits. Justlike the angular counterpart, the linear velocity is largestat perihelion and smallest at aphelion: v max = v c u , v min = v c u min (20) D. Orbital position as function of time
Rewriting (15) yieldsd ϕ = ω c u ( ϕ ) d t ⇒ t ( ϕ ) = T c π ϕ (cid:90) d ϕ (cid:48) u ( ϕ (cid:48) ) (21)This offers an explicit mapping between the polar angleof points along the orbit and the time elapsed to reachthem. Using the relation in the reverse direction providesthe position of the satellite at any given time. E. Orbital bounds
Previous stability analysis has shown that orbits gen-erated by the spatial exponents considered here (1 ≤ α ≤
2) ‘can be bounded or unbounded depending on the total(satellite) energy’ [7]. Here we establish explicit orbitalbounds for a given α as a function of u . As mentionedpreviously, the perihelion is prescribed by the initial con-ditions (14) themselves: r min = R c u max = R c u (22)The aphelion distance r max = R c /u min can be derivedfrom the conservation of total energy. Indeed, theamount of kinetic energy the satellite loses upon mov-ing from perihelion to aphelion must precisely equal thework performed against the gravitational field to increasethe distance to the primary from r min to r max : M (cid:0) v − v (cid:1) = r max (cid:90) r min G α M M r α d r (23)Invoking (4) and (20) leads to u min = 2 − u α = 2 u − u = α − (cid:0) u α − − u α − (cid:1) < α < u − u = 2 (ln u − ln u min ) α = 1 (24)Bounded orbits require strictly positive u min ( r max shouldremain finite), leading to the criterion u < (cid:18) α − (cid:19) / (3 − α ) < α ≤ ≤ u <
2, which produce bounded orbits for allexponents 1 ≤ α ≤ III. BOUNDED ORBITSA. Circular orbits
The master equation (13) supports a constant solution: u ( ϕ ) = 1 ↔ r ( ϕ ) = R c (26)which corresponds to the circular orbit discussed earlier.One easily verifies that inserting u ( ϕ ) = 1 into (15), (18)and (21) respectively recovers ω ( ϕ ) = ω c , v ( ϕ ) = v c and t ( ϕ = 2 π ) = T c as appropriate. B. Near-circular orbits
Let us consider initial conditions of type u = 1 + e with 0 < e (cid:28) u ( ϕ ) (cid:39) e cos( β e ϕ ) , β e = √ − α (28)as can be verified by substituting this form into the mas-ter equation, performing a series expansion 1 /u − α (cid:39) − (2 − α ) e cos( β e ϕ ), and equating the constant and firstharmonic terms on the left and right hand sides.In the Newtonian limit α = 2 (yielding β e = 1) thesolution becomes exact regardless the magnitude of e ,and corresponds to the familiar elliptical (0 < e <
1) orhyperbolic ( e >
1) orbits with the primary as focus.The trajectories generated for 0 < e (cid:28) ≤ α < − e (cid:46) ¯ r ( ϕ ) (cid:46) e . Interestingly,the orbits approximately follow hypotrochoids traced outby a circle with radius R c / ( β e −
1) inside a larger circlewith radius β e R c / ( β e −
1) (see Appendix A).
C. General solutions
The approximation (28) becomes progressively less ac-curate when the initial point moves inwards, as higher-order terms ought to be considered in the series expan-sion of 1 /u − α when e can no longer be considered as verysmall. However, since both d / d ϕ and fractional poweroperators conserve the periodicity of the input function,the master equation clearly supports periodic u ( ϕ ), andwe can postulate a general solution of the form u ( ϕ ) = ∞ (cid:88) n =0 A n cos( nβϕ ) (29)An approximate realisation truncated to the N lowestharmonics can conceptually be constructed as follows:1. Perform a series expansion 1 / (1 + χ ) b (cid:39) − bχ + b ( b + 1) χ − · · · of 1 /u − α up to N -th order.2. Rewrite resulting powers cos m ( nβϕ ) as linear com-binations of harmonics via de Moivre’s formula.3. Equate the coefficients of the static term and lowest N harmonics on either side.4. Enforce the initial condition: N (cid:80) n =0 A n = u .Steps 3 and 4 together produce a set of N + 2 polyno-mial equations in N + 2 unknowns, being N + 1 spectralcoefficients { A n } Nn =0 and the periodicity factor β . Ex-plicit formulae and numeric examples for N = 2 are listedin Appendix B. The procedure quickly becomes cumber-some for practical use as N increases. Accurate solutionscan instead be obtained by direct numerical integrationof the initial value problem, as will be done in Sec. IV.The presence of higher-order harmonics causes or-bits to quantitatively deviate more and more from truehypotrochoids as u − D. Closed periodic orbits
An interesting subset comprises solutions having ra-tional β = p/q ( p and q coprime integers). These yieldclosed orbits that are traversed cyclically every q revolu-tions, with extremal points forming 2 regular p -gons:perihelia : ϕ = mqp π ( m ∈ N ) (30)aphelia : ϕ = ( m + ) qp π ( m ∈ N ) (31)Each vertex is visited once per orbital cycle, in a sequencethat moves p − q point(s) clockwise along the polygons. IV. RESULTS
All results presented here were obtained by numericallysolving the initial value problem (13)–(14) through theDormand-Prince method with absolute and relative errortolerances set to 10 − . Software implementations of this4th/5th order explicit Runge-Kutta scheme are availableacross multiple applications and platforms, including the ode45 solver in Octave [8].Each orbit is plotted in dimensionless cartesian coor-dinates ( x/R c , y/R c ) i.e. normalised to the radius of thecircular orbit having the same angular momentum (dis-played for visual reference in grey shading). A. Orbital bounds of general solution
Figure 1 shows the dimensionless aphelion distance r max /R c as a function of the initial condition u for char-acteristic exponents α = 2 , . , . . . ,
1. For a given u (corresponding to the same initial normalised kinetic en-ergy v /v for all α ), the orbit is systematically moretightly bound to the primary as α becomes smaller, be-cause gravitational fields that decay less sharply in spacerequire more work to increase orbital distance. B. Periodicity of general solution
Figure 2 shows the periodicity factor β as a functionof the initial condition u for characteristic exponents α = 2 , . , . . . ,
1. Results were obtained to 4 significantdigits by determining the period 2 π/β of numerical so-lutions u ( ϕ ) computed with step size ∆ ϕ = 10 − rad.The β ( u ) dependence systematically intensifies towardslarger u values (more eccentric orbits) and more anoma-lous exponents (lower α values). We remind that theNewtionian case α = 2 is characterised by β ≡ u r m a x / R c = / u m i n α = 21.611.21.4 FIG. 1. Aphelion distance normalised to radius R c of thecircular orbit. Values observed in numerically computed or-bits (open circles) agree perfectly with the analytic result (24)derived from conservation of energy (lines). α = 2 β u FIG. 2. Periodicity factor β for the general solution u ( ϕ ) = (cid:80) A n cos( nβϕ ) of the initial value problem (13)–(14). Thetriangles on the left indicate β e = √ − α derived for near-circular orbits ( u (cid:39) β = p/q inducing closed periodic orbits. C. Geometric features
Figure 3 illustrates the close connection between theorbital shapes induced by anomalous gravitation lawsand hypotrochoids. The perihelia and aphelia of closedperiodic orbits (rational periodicity factor β ) moreoverform the vertices of regular (star) polygons (Fig. 4). −1.5 α = 1.5 , u = 1.2636 ( β ~ 8/7) x/R c y / R c α = 1.3 , u = 1.5930 ( β ~ x/R c y / R c α = 1.4 , u = 1.9188 ( β ~ x/R c y / R c −2 FIG. 3. Computed orbits (open circles) closely resemble hypotrochoids (solid lines). The quantitative agreement is best formoderate initial conditions ( u (cid:46) .
3) which induce orbits in the vicinity of the circular one (grey shading). More eccentricorbits (larger u ) have sharper lobes than their hypotrochoidal counterparts but a clear qualitative correspondence remains. −2.5
0p = p = p = p – q = α = 1.5 , u = 1.6200 ( β ~ x/R c y / R c α = 1.5 , u = 1.9131 ( β ~ x/R c y / R c α = 1.2 , u = 1.9693 ( β ~ x/R c y / R c −5
234 5 2 73 p – q = p – q = FIG. 4. Rational periodicity factors β = p/q with p and q coprime induce closed periodic orbits (black lines) whose perihelia(blue points) and aphelia (red points) form regular p -gons. A satellite revolving around the primary counter-clockwise visitseach of the vertices once per orbital cycle in a sequence that moves p − q point(s) clockwise along the polygons. D. Time-resolved closed periodic orbits
Figures 5 and 6 show various examples of closed pe-riodic orbits across a range of characteristic exponents α . Each orbit is plotted in 3 color-coded renderings thatrespectively show dimensionless angular velocity, linearvelocity, and time (normalised to circular orbit counter-parts ω c , v c and orbital period T c respectively). V. CONNECTION WITH MOND:THE CASE α = 1 Under Newtonian gravity the velocity of a star in circu-lar orbit around a galactic core reads v ( R ) = (cid:112) GM R /R where M R denotes the total galactic mass contained within the star’s orbital radius R . However, empiricalrotation curves v ( R ) do not decay as the expected 1 / √ R but become flat instead: v ( R ) ∼ constant [9]. Two mainexplanations for this this discrepancy between predictedand observed behaviour have emerged [9].On the one hand, the ‘dark matter’ hypothesis pos-tulates that galaxies possess haloes of unknown gravita-tionally interacting substance not detectable by currentmeans. This effectively adds mass to the M R value in-ferred from regular observable matter, which counteractsthe 1 / √ R velocity drop-off. Despite considerable effortsby cosmologists and particle phycisists, dark matter todate has never been observed directly, and remains a hy-pothetical substance unknown to physics [1].Modified Newtonian dynamics (MOND), on the otherhand, maintains nominal M R values but instead positsthat Newton’s second law breaks down over galacticscales, with forces becoming proportional to the squareof the acceleration when the latter is very small [10]. For-mally, MOND postulates that F = m µ (cid:18) aa (cid:19) a (32)Here µ ( x ) is an interpolation function with asympotes µ ( x (cid:28) ∼ x and µ ( x (cid:29) ∼
1, e.g. of the form µ ( x ) = 1[1 + (1 /x ) γ ] /γ , γ > a (cid:39) . × − m/s [10] is an empirical fitting con-stant acting as a characteristic acceleration below whichnon-Newtonian behaviour becomes important. Applying(32) with γ = 1 for simplicity to a circular orbit withradius R under Newtonian gravity shows that GM R R = a a + a = v /R v /R + a (34) ⇔ (cid:40) a (cid:29) a (Newtonian) : v = (cid:112) GM R /Ra (cid:28) a (deep MOND) : v = √ GM R a (35)where the latter is indeed independent of R .Alternatively, one could maintain the second law butinstead modify the gravitational force, say F ( r ) = GM M r (cid:18) r R (cid:19) (36)At small distances r (cid:28) R the force is Newtonian, whilelarge distances r (cid:29) R experience anomalous gravitationof type (1) with exponent α = 1 and fractional constant G α = G/ R . Applying Newton’s unmodified second lawto a circular orbit with radius R now yields (cid:40) R (cid:28) R ( α = 2) : v = (cid:112) GM R /RR (cid:29) R ( α = 1) : v = (cid:112) GM R / R (37) where once again the latter is independent of orbital ra-dius. Comparison with (35) shows that R = (cid:114) GM R a (38)For galactic bulge masses typically on the order of 10 Suns ( M R (cid:39) × kg) we find R (cid:39) . (cid:39) ,
000 light years (39)which is indeed comparable to the orbital distances atwhich rotation curves begin to flatten out [9].
VI. CONCLUSIONS
We deduced detailed characteristics of bounded orbitsunder anomalous gravitational forces F ∼ /r α with spa-tial exponents 1 ≤ α <
2. The familiar Keplerian ellipsesmake way for hypotrochoid-like self-intersecting trajecto-ries. Despite their convoluted appearance for most initialconditions, the orbits are governed by a general solutionthat is, just like the Newtonian one, periodic in polar an-gle. A subset of cases exhibit a ‘mathematical resonance’in which this periodic solution is traversed an integer p times for every q revolutions. These lead to fully closedperiodic orbits with the primary as centroid and periheliaand aphelia forming regular p -gons.Overall, the key differences between anomalous andNewtonian gravity could be summarised in terms of aheightened sensitivity to initial conditions for the for-mer. Indeed, moving the launch site across the boundedrange will drastically alter the orbit’s geometric appear-ance (whether it is closed or not, and if yes how many‘lobes’ it contains), while the Keplerian counterpart re-mains elliptical but simply with a different eccentricity.Some features do remain preserved, as well: all spatialexponents support circular orbits, albeit with distinctradii and orbital velocities. The latter becoming inde-pendent of the radius for a gravitational force inverselyproportional to distance, the anomalous case α = 1 iscompatible with galactic rotation curves and conceptu-ally linked to Modified Newtonian Dynamics. [1] I. Stewart, Calculating the cosmos: How mathematics un-veils the universe. (Profile Books, London, 2016).[2] A. Koyr´e, Isis What If? Serious scientific answers to ab-surd hypothetical questions. (Houghton Mifflin Harcourt,New York, 2014).[6] H. Goldstein, C. Poole, and J. Safko,
Classical Mechanics(Third Edition). (Addison Wesley, San Francisco, 2000).[7] S. Ray and J. Shamanna, arXiv:physics/0410149 (2004). [8] http://octave.sourceforge.io/octave/function/ode45.html(consulted April 2, 2018).[9] K. Begeman, A. Broeils, and R. Sanders, MNRAS ,523 (1991).[10] M. Milgrom, Astrophys. J. , 371 (1983). −3.5 0 ω / ω c α = 1 , u = 1.9336 ( β ~ 3/2) −3.5 0 3.5−3.503.5 −3.5 0 3.5−3.503.5041 v / v c t / T c ω / ω c v / v c t / T c ω / ω c v / v c t / T c ω / ω c v / v c t / T c α = 1.5 , u = 1.6200 ( β ~ 5/4) −2.5 0 2.5−2.502.5 −2.5 0 2.5−2.502.502.81 0.41.71 06.541−2 0 2−202 α = 1.3 , u = 1.5930 ( β ~ 4/3) −2 0 2−202 −2 0 2−2020.22.610.22.61 0.41.61 04.411−2.4 0 2.4−2.402.4 α = 1.6 , u = 1.5592 ( β ~ 6/5) −2.4 0 2.4−2.402.4 −2.4 0 2.4−2.402.40.41.61 07.641x/R c y / R c x/R c y / R c x/R c y / R c x/R c y / R c x/R c y / R c x/R c y / R c x/R c y / R c x/R c y / R c x/R c y / R c x/R c y / R c x/R c y / R c x/R c y / R c FIG. 5. Examples of closed periodic orbits. The color coding shows angular velocity (left column), linear velocity (middlecolumn), and time (right column), each normalised to counterparts of the circular orbit r = R c (grey shading). Time-steppedvideo animations of each of these orbits are available as supplementary media files. −1.5 ω / ω c x/R c y / R c x/R c y / R c x/R c y / R c x/R c y / R c x/R c y / R c x/R c y / R c x/R c y / R c x/R c y / R c x/R c y / R c x/R c y / R c x/R c y / R c x/R c y / R c α = 1.7 , u = 1.2636 ( β ~ −1.5 v / v c t / T c ω / ω c v / v c t / T c ω / ω c v / v c t / T c ω / ω c v / v c t / T c α = 1.1 , u = 1.7191 ( β ~ −2.8
31 0.41.81 α = 1.4 , u = 1.5284 ( β ~ −2.2 α = 1 , u = 1.5390 ( β ~ −2 FIG. 6. Examples of closed periodic orbits. The color coding shows angular velocity (left column), linear velocity (middlecolumn), and time (right column), each normalised to counterparts of the circular orbit r = R c (grey shading). Time-steppedvideo animations of each of these orbits are available as supplementary media files. Appendix A: Link with hypotrochoids
Let us consider a circle with radius R rolling inside alarger circle with radius R = βR . (Fig. 7). R dR = β R ( β -1) θ xy θ FIG. 7. Construction of hypotrochoidal curve (purple) tracedout by a point (marked × ) a distance d from the center ofa smaller circle with radius R (blue) rolling without slip-ping inside a larger circle with radius R = βR (red). Thedepicted example produces the curve shown in Fig. 3b. A point on the small circle at distance d from its centerand initially between the two circle centers traces a curveparametrised by (cid:40) x = ( β − R cos θ − d cos[( β − θ ] y = ( β − R sin θ + d sin[( β − θ ] (A1)where θ is the angle marking the position of the centerof the small circle. In polar coordinates we have r ( ϕ ) = ( β − R + d − β − R d cos[ βθ ( ϕ )]tan ϕ = ( β − R sin θ + d sin[( β − θ ]( β − R cos θ − d cos[( β − θ ] (A2)Setting R = 1 / ( β −
1) and assuming d (cid:28) (cid:40) r ( ϕ ) (cid:39) − d cos( βθ )tan ϕ (cid:39) sin θ cos θ = tan θ ↔ ϕ (cid:39) θ (A3)Meanwhile, solutions u ( ϕ ) (cid:39) e cos( βϕ ) with e (cid:28) r ( ϕ ) ≡ u ( ϕ ) (cid:39) − e cos( βϕ ) (A4) which describes the same curve as (A3) if setting d = e .For a given general solution u ( ϕ ) (not necessarily in-ducing a near-circular orbit), the hypotrochoid havingthe same periodicity factor β and same extremal radii¯ r min = 1 /u , ¯ r max = 1 /u min = 1 /u ( π/β ) is generated by R = R β = ¯ r max + ¯ r min β − , d = ¯ r max − ¯ r min Appendix B: Explicit formulae for general solutiontruncated to 2 harmonics
We seek an approximate solution of the initial valueproblem (13)–(14) of the form u ( ϕ ) (cid:39) A + A cos( βϕ ) + A cos(2 βϕ ) (B1)Executing the procedure outlined in Sec. III C producesthe following set of equations: A − α = 1 + (2 − α )(3 − α ) (cid:104) A A + A A (cid:105) − ( β − A A − α = − (2 − α ) A A (cid:104) − (3 − α ) A A (cid:105) − (4 β − A A − α = − (2 − α ) A A (cid:104) A A − (3 − α ) A A (cid:105) A + A + A = u (B2)These can be solved numerically for given pairs of gravi-tational exponent α and initial value u to obtain theperiodicity factor β and harmonic weights A , A , A . Ascould be expected, the truncated solution performs verywell for moderate initial conditions u − (cid:46) . α = 1 . u = 1 . β = 1 . − . A = 1 . − . A = 0 . . A = − . − . u = 1 . β = 1 . − . A = 1 . − . A = 0 . . A = − . − . u = 1 . β = 1 . − . A = 1 .
027 ( − . A = 0 .
482 (+0 . A = − . −−