A geometric method to locate Neptune
AA geometric method to locate Neptune
Siddharth Bhatnagar ∗ and Jayanth Vyasanakere P. † The School of Arts and Sciences, Azim Premji University, Bengaluru, 562 125, India
Jayant Murthy ‡ Indian Institute of Astrophysics, Bengaluru, 560 034, India (Dated: February 9, 2021)We develop a direct geometric method to determine the orbital parameters and mass of a planet,and we then apply the method to Neptune using high-precision data for the other planets in thesolar system. The method is direct in the sense that it does not involve curve fitting. This paper,thereby, offers a new pedagogical approach to orbital mechanics that could be valuable in a physicsclassroom.
I. INTRODUCTION
The discovery of Neptune in 1846 by observing devia-tions of Uranus from its predicted orbit remains one ofthe crowning achievements of Newtonian mechanics.
The discoverers faced a difficult inverse problem, eas-ily solvable today with electronic computers, but Her-culean in the mid-nineteenth century.
The discover-ers, Urbain Le Verrier and John Couch Adams, indepen-dently investigated Alexis Bouvard’s ephemeris tables forUranus to check their accuracy and found no computa-tional errors.
With this, it became clear that the prob-lem lay not with the tables, but rather with the planetitself! The mathematical problem of locating the preciseposition and determining the mass of the perturbing ob-ject was non-trivial;
Le Verrier and Adams proceededby linearizing the system of equations involved and solv-ing them by least squares. Computing the semi-majoraxis of the perturbing body was also a complex problem,for which an informed estimate had to be made. This wasdone using the Titius-Bode law, which held for the fiveplanets known since antiquity. Neptune has subsequentlybeen “rediscovered” with simpler methods several timesafter its real discovery. In today’s age of high-precision data, orbital state vec-tor data of most solar system bodies are available. Wehave retrieved these data for the major planets in
Carte-sian coordinates , with the origin at the instantaneous center of the Sun (see the Appendix for more details).These data, collected at a time step of 2 h, range fromthe date of Uranus’ discovery (March 1781) to the presentday (March 2020). It is clear that we have an advantageover the discoverers, in terms of the nature, quality andduration of data at hand. Using this advantage, we de-scribe a purely geometric way of locating a planet, focus-ing on Neptune.Since we work with state vector data and not the or-bital parameters of the known planets, our approach nat-urally lends itself to vector analysis. The geometricalideas are encoded into scalar (dot) products and vector(cross) products. The method emphasizes basic ideasof classical mechanics, such as Newton’s laws, Kepler’slaws, and motion in the presence of an inverse squarelaw force, but does not use curve fitting to locate theunknown planet. All these ideas are typically encoun-tered in an undergraduate physics curriculum, makingour method better suited to pedagogical purposes thanthe historical ones.
II. EQUATION OF MOTION
Our model comprises only nine objects: the Sun, theseven known planets at the time (Mercury to Uranus) andNeptune, which we consider unknown for the purpose ofdemonstrating our method. Newton’s law written for a r X i v : . [ phy s i c s . c l a ss - ph ] F e b the i th object reads d x i dt = − (cid:88) j (cid:54) = i GM j x i − x j | x i − x j | ( i = 1 to 9) , (1)where G is the universal gravitational constant, x i s and x j s stand for the position of the center of mass (COM) ofthe objects (including their moons, if any), and M j s arethe masses (again including the moons, if any). However,we need not worry about solving these coupled differen-tial equations, for we have the solution in terms of datafor all objects except Neptune, whose parameters needto be determined using these equations.The d x i /dt in Eq. (1) denote accelerations with re-spect to an inertial frame. The COM of the solar systemis a reasonably good choice for such a frame. However,the COM is not an observable, but rather a constructthat can only be obtained after having knowledge of all constituents of the system, including Neptune ( N ), whichis yet to be “discovered.” We circumvent this problemby subtracting Eq. (1) written for the Sun ( (cid:12) ) from thesame equation written for Uranus ( U ). This gives GM N (cid:18) r N − r U | r N − r U | − r N | r N | (cid:19) = V ( t ) , (2)where V ( t ) = d r U dt + G (cid:0) M (cid:12) + M U (cid:1) r U | r U | − (cid:88) j (cid:54) = (cid:12) ,U,N GM j (cid:18) r j − r U | r j − r U | − r j | r j | (cid:19) (3)and where r i = x i − x (cid:12) , which are relative coordinateswith respect to the Sun. It is instructive to write downall the terms explicitly to understand the derivation ofEqs. 2 and 3, and to note that r i − r j = x i − x j .In Eq. (3), the r i s are observable without a priori knowledge of Neptune. Thus, with Eq. (2) and the as-sembled data from Ref. 13 (as elaborated in Sec. I), V ( t )is known; see Sec. III for a physical understanding ofthis vector. The stage is now set to find the unknowns r N ( t ) and M N , which we do through a purely geometricapproach. H a L H b L H c L H d L FIG. 1. The behavior of ˆ V as Neptune (blue sphere—the oneon the outer orbit) moves around in its orbit, with Uranus(brown sphere—the one on the inner orbit) held fixed. TheSun (not shown) is at the center and the black dashed lineis the z -axis. The orbits of Uranus (brown inner orbit) andNeptune (blue outer orbit) are both shown to be slightly in-clined to the x – y plane. The red solid line denotes ˆ V . Thegreen dot-dashed line denotes vector p (Eq. (4)). (a) and (c)correspond to conjunction and opposition of the Sun–Uranus–Neptune system, respectively. (b) and (d) correspond to con-figurations, which are not conjunction or opposition, but stillsatisfy χ = 0 (see Eq. (5) and the discussion around it). Notethat this figure is a schematic for the sake of visualizing ˆ V and was not derived from any data. III. UNDERSTANDING THE VECTOR V Since our discussion revolves around the vector V , itis important to discuss what it means physically. FromEq. (2), V ( t ) would be zero if Neptune did not exist; itis that part of the acceleration of Uranus not explainedby the other planets.Consider ˆ V , which is the unit vector along V . InFig. 1, this is denoted by a red solid vector. If r U isimagined to be held fixed and r N is made to rotate an-ticlockwise once around the z axis, then, simultaneously,ˆ V will complete two anticlockwise circuits around the z axis. We will revisit ˆ V in Sec. IV.We now turn our attention to V ( t ) = | V ( t ) | , the mag-nitude of V . This is shown in Fig. 2. With V being V ( − AU / y r ) t (yr CE) FIG. 2. V , which is the magnitude of V (Eq. (3)), is shownas a function of time. The two peaks correspond to November1822 and June 1994. As argued in Sec. III, these peaks providea rough estimation of the time of conjunction of the Sun–Uranus–Neptune system. that part of Uranus’ acceleration arising exclusively fromNeptune, the peaks in Fig. 2 are expected to correspondto conjunctions. This can be confirmed by modeling theorbits of Uranus and Neptune as co-planar circles aroundthe Sun, and obtaining from Eq. (2) that V peaks dur-ing conjunction, irrespective of the radii of the orbits.Hence, it is reasonable to expect that the peaks in Fig. 2,which occurred in November 1822 and June 1994, signifyconjunctions. However, this prediction is not precise asthe planets’ orbits are neither exact circles nor exactlyco-planar. Therefore, the analysis, thus far, provides buta rough estimate for the date of conjunction. IV. DETERMINATION OF CONJUNCTIONAND OPPOSITION WITH URANUS
To obtain precise conjunction and opposition times, werecall that in a spherical polar coordinate system withthe Sun at the origin, the azimuthal angle ( φ ) for bothUranus and Neptune will be the same during conjunc-tion, and differ by π radians at opposition. From Eq. (2)and Fig. 1, it is clear that during both conjunctions and oppositions, the vector V will have the same φ coordi-nate as Uranus. This observation motivates us to definea vector p ( t ) = ˆ z × ˆ r U ( t ) , (4) -1.0-0.50.00.51.0 C O C χ t (yr CE) FIG. 3. Precise determination of the conjunction and opposi-tion of the Sun–Uranus–Neptune system. χ (Eq. (5)) is shownas a function of time. Among the several roots of χ in the du-ration shown, the ones occurring in August 1821 and March1993 determine conjunctions ( C and C ), while the one inMay 1908 determines opposition ( O ). The origin of the re-maining roots (also depicted in parts (b) and (d) of Fig. 1) isdiscussed below Eq. (5). which lies in the x – y plane and is normal to r U ( t ) (de-noted as a green dot-dashed line in Fig. 1). In terms of p ,we can hope to determine a conjunction or an oppositionby checking whether V is perpendicular to p at a giventime. In other words, we define a quantity χ , which is ameasure of the projection of ˆ V along p , as χ ( t ) = ˆ V ( t ) · p ( t ) , (5)and look for its roots.Figure 3 shows χ as a function of time. As expected, χ does vanish around the estimated times for conjunctionsand oppositions as obtained in the Sec. III by lookingat the peaks in Fig. 2. We have already seen that the φ coordinate of V is the same as that of Uranus duringboth conjunctions and oppositions. This observation canalso be used to further verify whether a given root of χ corresponds to either a conjunction or an opposition.The remaining roots of χ , which do not correspond toconjunctions or oppositions can be understood as follows.During conjunctions and oppositions, the projections ofboth the bracketed terms in Eq. (2) on p vanish. How-ever, χ can also vanish if the projections of these twoterms cancel each other. The relative dominance of theseterms flips when the system moves from a conjunction to -0.020.000.020.040.06 ξ t (yr CE) FIG. 4. The projection of V ( t ) along the normal to the orbitalplane of Uranus, quantified through ξ (Eq. (6)), is shown asa function of time. The roots of ξ in the duration shownoccur at June 1791, September 1850, June 1860, March 1932,April 1962, and May 2015. As argued in Sec. V, the rootsin September 1850 and May 2015 can be used to preciselydetermine the orbital period of Neptune. an opposition or vice versa, resulting in an extra root be-tween them. Basically, if the n th root corresponds to aconjunction, the ( n + 2)th root will be an opposition, the( n + 4)th root will again be a conjunction, and so on.In this manner, we conclude that the roots of χ inFig. 3, which happen to fall in August 1821 and March1993, are conjunctions and designate them as C and C , respectively. The root falling in May 1908 then cor-responds to an opposition, which we call O . V. DETERMINATION OF THE ORBITALPERIOD AND SEMI-MAJOR AXIS
The synodic period of Neptune with respect to Uranusis the time difference between the conjunctions C and C : T sy = 171 . T N ≈
165 years. However, thisis an approximation.We, therefore, venture to obtain T N precisely by defin-ing ξ ( t ) = ˆ V ( t ) · ˆ n U , (6)where ˆ n U is a unit vector along (cid:0) r U × d r U dt (cid:1) evaluated at our model’s epoch, March 1781 (referred to hereafter as T I ). In other words, ˆ n U is a unit normal to the orbitalplane of Uranus. The quantity ξ , a measure of projectionof V ( t ) along ˆ n U , is shown in Fig. 4 as a function oftime. When Neptune crosses the orbital plane of Uranus, ξ vanishes. However, similar to the roots of χ , betweenany two consecutive roots of this kind will exist anotherkind of root, where the components along ˆ n U of the twoterms bracketed in Eq. (2) cancel each other. Thus, withreference to Fig. 4, if the time difference between the n thand ( n + 4)th roots of ξ corresponds to T N , then thatbetween the ( n + 1)th and ( n + 5)th roots corresponds to T sy , and vice versa. In particular, the roots in June 1791and April 1962 differ by 170.8 years, which is very close to T sy . Hence, the roots in September 1850 and May 2015can be used to give T N = 164 . a ) is then obtained from T N using Kepler’s thirdlaw (assuming M N (cid:28) M (cid:12) ) to be 30.05 astronomicalunits (AU). VI. DETERMINATION OF OTHER ORBITALPARAMETERS AND MASS
To begin here, the polar angle θ and the azimuthal an-gle φ (the spherical polar coordinates) of Neptune duringa conjunction with Uranus can be determined as follows.Equation 2 can be re-expressed as r N = r U + | r N − r U | (cid:18) r N | r N | + (cid:12)(cid:12)(cid:12)(cid:12) r N − r U | r N − r U | − r N | r N | (cid:12)(cid:12)(cid:12)(cid:12) ˆ V (cid:19) . (7) As an initial guess, we can take θ for Neptune to be π/ φ of Neptune can be taken equal to that of Uranus. Byfreezing Neptune’s r coordinate to its semi-major axis a ,which has already been obtained, we have a guess for r N . With this, every quantity on the RHS of Eq. (7) isknown. Hence, the LHS of this equation can be takenas a better estimate for r N . This way, we can iterateEq. (7) to obtain φ and θ of Neptune at conjunction.A similar procedure can be carried out for an opposi-tion by re-expressing Eq. (2) as TABLE I. The converged directions (as elaborated in Sec. VI)of conjunctions (Uranus and Neptune aligned on the same sideof the Sun) and opposition (Uranus and Neptune aligned onthe opposite sides of the Sun). φ and θ are the spherical polarcoordinates, while ψ is the angle in the orbital plane from C . Date (year CE) φ (rad) θ (rad) ψ (rad) C August 1821 4.532 1.512 0 O May 1908 1.582 1.647 3.334 C March 1993 4.787 1.490 0.256 r N = | r N | (cid:18) r N − r U | r N − r U | − (cid:12)(cid:12)(cid:12)(cid:12) r N − r U | r N − r U | − r N | r N | (cid:12)(cid:12)(cid:12)(cid:12) ˆ V (cid:19) (8) and taking φ for Neptune to be π radians away fromthat of Uranus.Note that we have re-expressed Eq. (2) in two differentways: one for conjunction (Eq. (7)) and the other foropposition (Eq. (8)). These choices ensure convergenceof the iteration algorithm in the respective cases. Theconverged values of φ and θ of Neptune after iterationare given in Table I. Thus, we obtain the directions ofNeptune during C , O , and C .At this stage, we assume that the orbit of Neptuneis an ellipse with the Sun at one focus (Kepler’s firstlaw). Since the time difference between C and C isgreater than the orbital period, it is clear that C is aheadof C in the orbit. The orbital plane is characterizedby the inclination ( i ) of its normal from the z axis andthe longitude of the ascending node (Ω). These can beevaluated from the knowledge of the directions of C and C as i = 0 . . O should be very close to this plane. Anegligible component perpendicular to this plane, if any,can be dealt with by projecting O onto the plane. Let ψ be the angle measured in the plane of the ellipse from C . ψ for C and O can be easily found, since we alreadyknow the corresponding φ and θ . These values of ψ aretabulated in Table I.Using this information, the equation of the ellipse canbe written with unknowns being its eccentricity ( (cid:15) ) andthe distance of Neptune from the Sun ( r I ) at T I . These r I ( AU ) (cid:15) O : 1908 C : 1993 FIG. 5. Determination of eccentricity ( (cid:15) ) and the distance ofNeptune from the Sun ( r I ) at our model’s epoch ( T I ). Thesolid line represents the combination of (cid:15) and r I that is con-sistent with ψ during O . The dashed line is the combinationof (cid:15) and r I that is consistent with ψ during C . Their inter-section determines (cid:15) and r I of Neptune’s orbit. unknowns can be obtained by demanding that Neptune’sorbit assumes the appropriate ψ values given in Table Iat the respective times. The combination of (cid:15) and r I that are separately con-sistent with ψ during O and C is obtained as shownin Fig. 5. The intersection of these curves gives the re-quired combination of (cid:15) and r I . This gives (cid:15) = 0 . r I = 30 .
24 AU. The orbit of Neptune is thus completelydetermined.The time at which Neptune passes through the peri-helion ( T p ) and the argument of the perihelion (angularcoordinate of the perihelion in the orbital plane from theascending node, denoted by ω ) emerge as November 1892and 3.282 radians, respectively.Thus, we now know all the orbital parameters of Nep-tune, from which its location at any time can be deter-mined.The only remaining parameter is the mass of Neptune.From Eq. (2), this can be written as GM N = V / (cid:12)(cid:12)(cid:12)(cid:12) r N − r U | r N − r U | − r N | r N | (cid:12)(cid:12)(cid:12)(cid:12) . (9)Now that we know r N , we can substitute it into Eq. (9)to obtain M N . The mass should of course be a constant,but any deviation of the position obtained from the exactvalue will cause a variation in M N . Hence, we averagedthe M N values so obtained over the orbital period toarrive at M N = 1 . × kg. VII. COMPARISON WITH ACTUAL VALUES If r N ( t ) is the trajectory obtained here and r N ( t ) isthe actual trajectory of Neptune given by Ref. [13], oneway to characterize the error is% deviation = | r N ( t ) − r N ( t ) || r N ( t ) | × . (10)The deviation of r N ( t ) from r N ( t ) varies with time andis always within 1.7%. As viewed from Earth, the pre-dicted direction of Neptune always lies within one degreeof its actual direction.Another approach to testing the quality of our solutionis by characterizing the elliptical orbit through its sixelements. The inclination ( i ) of the orbital axis and thelongitude of the ascending node (Ω) specify the plane ofthe orbit. The semi-major axis ( a ) and eccentricity ( (cid:15) )specify the size and shape of the ellipse. The argument ofthe perihelion ( ω ) specifies the orientation of the ellipsein its plane. The time at which Neptune passes throughthe perihelion ( T p ) specifies the epoch.At any given time, the instantaneous position ( r N ( t ))of Neptune as seen from the Sun is already obtained fromRef. [13]. From this, the angular momentum ( L ) and theenergy ( E ) of Neptune can be calculated. While (cid:15) and a can be obtained from E and | L | , i and Ω can be obtainedfrom ˆ L . The calculation of T p and ω is more involved, andcan be determined from the eccentric anomaly. Hence,at each instant, the orbital elements can be calculated. IfNeptune and the Sun were isolated, the orbital elementswould be constants. However, due to perturbations byother planets, they do vary. We summarize these vari-ations in Table II, which also lists the orbital elementscalculated using r N ( t ) obtained from our method. It isapparent that they compare well with the correspond-ing actual values. Note that since r N ( t ) is presumed tobe perfectly elliptical, we obtain definite values for theorbital parameters, not a range of values. VIII. CONCLUSION
The discovery of Neptune was a magnificent exampleof mathematical and scientific analysis. While the oldermethods were certainly able to locate Neptune, they wereless transparent in terms of the usage of physics principlesand geometrical arguments than that advanced here.Aided with modern-age data, the method we have de-scribed to locate Neptune is simple, both conceptuallyand by means of calculation. Moreover, it is also a direct geometric method, without any curve fitting or solving ofdifferential equations. All steps in this method are basedon ideas encountered in an undergraduate curriculum,such as vector analysis and the laws of planetary mo-tion. Thus, our method offers a pedagogical alternativeto traditional methods.This method can also be applied to the case of a hith-erto undiscovered planet. For this, sound data for a longenough duration (at least until one time period of theplanet) of all the other planets in the system are required.This method can also serve to provide an initial guess formore sophisticated analyses, be it for application withinour solar system or to exoplanetary systems.
ACKNOWLEDGMENTS
The authors thank Rajaram Nityananda of Azim Pre-mji University for illuminating discussions and the threeanonymous reviewers for their helpful suggestions.
Appendix: DOWNLOADING THE DESIREDDATA
The following settings were made on Ref. 13 to gener-ate and download the data.(1) Choose “Vectors” from the “Ephemeris Type”menu. This will “generate a Cartesian state vec-tor table of any (solar system) object with respectto any major body.”
TABLE II. Comparison of Neptune’s orbital parameters and mass, obtained using the method illustrated in this paper, withthe actual values of the osculating orbit (as elaborated in Sec. VII) in the time interval considered. i (rad) Ω (rad) a (AU) (cid:15) ω (rad) T p (year CE) M N (kg)Our method 0.1124 3.9806 30.05 0.0143 3.282 Nov 1892 1 . × Actual 0.1122–0.1123 3.9769–3.9778 29.93–30.31 0.002–0.016 1.918–3.633 Aug 1855–Aug 1901 1 . × (2) Choose the object in question from the “TargetBody” menu.(3) Write “@sun” in the box provided in the “Coordi-nate Origin” menu. This will ensure that the coor-dinate origin is fixed to the center of the Sun at alltimes.(4) From the “Time Span” menu, choose a time span inthe appropriate format (specified in the same web-page next to these fields) and the necessary timestep.(5) Under the “Table Settings” menu, select vec-tor table output as “Type 2 (state vector { x,y,z,vx,vy,vz } ).” Then, choose the required out- put units from the “output units” submenu. Now,choose “body mean equator and node of date”from the “reference plane” submenu. This will setthe reference x – y plane to coincide with the Sun’smean equator. Then, choose “ICRF/J2000.0” fromthe “reference system” submenu. Leave “aberra-tions” as “Geometric states (no aberrations; instan-taneous ephemeris states).” Check the “labels,”“CSV format,” and “object page” options.(6) From the “Display/Output” menu, check the“download/save” option so that the ephemeris re-sults can be saved to a local file.Then select the “Generate Ephemeris” option and a .txtfile with the ephemeris data will be downloaded. ∗ [email protected] † [email protected] ‡ [email protected] All dates mentioned in this paper are in CE. U. J. Le Verrier, “Recherches sur les mouvements d’Uranuspar UJ Le Verrier,” Astron. Nachr. , 53–80 (1846). John Couch Adams, “Explanation of the observed irregu-larities in the motion of Uranus, on the hypothesis of dis-turbance by a more distant planet; with a determinationof the mass, orbit, and position of the disturbing body,”Mon. Not. R. Astron. Soc. , 149–152 (1846). David Rines, “The discovery of the planet Neptune,” Pop.Astron. , 482–498 (1912). W. M. Smart, “John Couch Adams and the discovery ofNeptune,” Nature (4019), 648–652 (1946). Siddharth Bhatnagar and Jayant Murthy, “Game oforbits: A gaming approach to Neptune’s discovery,”arXiv:1807.11280 (2018). Michael Martin Nieto,
The Titius-Bode Law of PlanetaryDistances: Its History and Theory , 1st ed. (Elsevier, 2014). Ernest W. Brown, “On a criterion for the prediction of anunknown planet,” Mon. Not. R. Astron. Soc. , 80–101(1931). Raymond Arthur Lyttleton, “The rediscovery of Neptune,”Vistas Astron. , 25–46 (1960). C. J. Brookes, “On a criterion for the prediction of Nep-tune,” Mon. Not. R. Astron. Soc. , 79–83 (1972). Hon Ming Lai, C. C. Lam, and Kenneth Young, “Perturba-tion of Uranus by Neptune: A modern perspective,” Am.J. Phys. (10), 946–953 (1990). Gustav Eriksson and Kevin Garcia Martin, “Discovery ofNeptune,” TRITA-SCI-GRU No. 2018-210 (KTH, Schoolof Engineering Sciences (SCI), 2018). J. D. Giorgini and JPL Solar System DynamicsGroup, see
Classical Mechanics , 3rd ed. (Pearson Education, 2007),p. 92–103. Actually, ζ = sign of the radial component of velocity at T I is another unknown parameter while solving for theellipse. With other parameters fixed, we tried both possi-bilities ( ζ = ±
1) and have kept ζ = +1 as there is no rootanalogous to that found in Fig. 5 for ζ = −−