Prospects of dynamical determination of General Relativity parameter beta and solar quadrupole moment J2 with asteroid radar astronomy
DDraft version November 12, 2018
Typeset using L A TEX twocolumn style in AASTeX61
PROSPECTS OF DYNAMICAL DETERMINATION OF GENERAL RELATIVITY PARAMETER β ANDSOLAR QUADRUPOLE MOMENT J (cid:12) WITH ASTEROID RADAR ASTRONOMY
Ashok K. Verma, Jean-Luc Margot,
1, 2 and Adam H. Greenberg Department of Earth, Planetary, and Space Sciences, University of California, Los Angeles, CA 90095, USA Department of Physics and Astronomy, University of California, Los Angeles, CA 90095, USA
Submitted to ApJABSTRACTWe evaluated the prospects of quantifying the parameterized post-Newtonian parameter β and solar quadrupolemoment J (cid:12) with observations of near-Earth asteroids with large orbital precession rates (9 to 27 arcsec century − ).We considered existing optical and radar astrometry, as well as radar astrometry that can realistically be obtainedwith the Arecibo planetary radar in the next five years. Our sensitivity calculations relied on a traditional covarianceanalysis and Monte Carlo simulations. We found that independent estimates of β and J (cid:12) can be obtained withprecisions of 6 × − and 3 × − , respectively. Because we assumed rather conservative observational uncertainties,as is the usual practice when reporting radar astrometry, it is likely that the actual precision will be closer to 2 × − and 10 − , respectively. A purely dynamical determination of solar oblateness with asteroid radar astronomy maytherefore rival the helioseismology determination. Keywords: astrometry — gravitation — minor planets, asteroids: general — relativistic processes —Sun: fundamental parameters — techniques: radar astronomy
Corresponding author: Ashok Kumar [email protected] a r X i v : . [ a s t r o - ph . E P ] A ug Verma et al. 2017 INTRODUCTIONThe parameterized post-Newtonian (PPN) formal-ism is a useful framework for testing metric theoriesof gravity (Will 2014). It consists of 10 dimension-less parameters that describe the general properties ofthe metric. In general relativity (GR), only 2 of the10 parameters are non-zero. They are known as theEddington − Robertson − Schiff parameters γ and β . γ represents the amount of curvature produced by a unitmass, and β represents the amount of nonlinearity inthe superposition law for gravity.Several techniques have been used to place obser-vational bounds on these parameters (Will 2014), in-cluding observations of the bending and delay of lightby spacecraft tracking (Bertotti et al. 2003, e.g.,) orVery Long Baseline Interferometry (e.g., Lambert &Le Poncin-Lafitte 2009), and fitting of ephemerides toobservations of planetary positions (e.g., Folkner 2009;Fienga et al. 2011; Verma et al. 2014; Fienga et al. 2015).In GR, γ and β are equal to one. Doppler track-ing of the Cassini spacecraft has shown that γ doesnot differ from one by more than 2 × − (Bertottiet al. 2003). Ephemeris-based studies prior to 2009 indi-cated that β − − (Folkner 2009; Pitjeva & Pitjev 2014). More re-cently, the availability of precise ranging data from theMESSENGER Mercury orbiter (Solomon et al. 2001)enabled improved estimates (Verma et al. 2014; Fiengaet al. 2015; Park et al. 2017). Here, we evaluate theprospect of asteroid orbit precession measurements toplace more stringent bounds on β . We consider Earth-based radar observations of near-Earth asteroids withperihelion shifts larger than 10 arcsec century − .Orbital precession can also be caused by the nonuni-formity of the gravity field that results from the oblateshape of the Sun. The solar oblateness is character-ized by the solar quadrupole moment, J (cid:12) (e.g., Kaula2000). Simultaneous estimation of β and J (cid:12) requiresthat the precessional effects due to GR and to the Sun’soblateness be disentangled. Fortunately, GR is a purelycentral effect, whereas the oblateness-induced precessionhas an inclination dependence. The two effects also havea different distance dependence (Misner et al. 1973). Asa result, observations of a small sample of near-Earthasteroids with a variety of semi-major axes and inclina-tions (Table 1) can in principle be used to estimate β and J (cid:12) (Margot 2003; Margot & Giorgini 2009).Current estimates of the solar quadrupole momentare typically derived on the basis of interior modelsof the Sun constrained by helioseismology data (e.g.,Mecheri et al. 2004; Antia et al. 2008). The currentbest value from the helioseismology literature is J (cid:12) = (2 . ± . × − (Will 2014). Dynamical estimates thatdo not rely on fits to helioseismology data yield similarvalues of J (cid:12) = 2 . ± . × − (Fienga et al. 2015)and J (cid:12) = 2 . ± . × − (Park et al. 2017). High-precision dynamical estimates are important to validateour understanding of the interior structure of the Sun.Our simulations of the determination of β and J (cid:12) using a variety of asteroid orbits suggest that indepen-dent values of β and J (cid:12) can be obtained with satisfac-tory precision: with the traditionally conservative as-signment of radar uncertainties, β can be constrained atthe 6 × − level and J (cid:12) can be constrained at the3 × − level. With uncertainties that more closelyreflect measurement errors, this precision may be im-proved by a factor of ∼
3. (Section 4).The outline of this paper is as follows. In Section 2,we describe our choice of target asteroids. In Section 3,we discuss the estimation of asteroid orbits with opti-cal and radar measurements. Our dynamical model anddata reduction procedures are described in Section 3.1and 3.2, respectively. Orbit determination results arepresented in Section 3.3. Simulations of the determina-tion of β and J (cid:12) are described in Section 4. TARGET ASTEROIDSThe per-orbit secular advance in the angular positionof the perihelion is given by (Misner et al. 1973) δω = 6 πGM (cid:12) a (1 − e ) c (cid:20) (2 − β + 2 γ )3 (cid:21) + 6 π R (cid:12) (1 − / i ) a (1 − e ) J (cid:12) , (1)where ω is the argument of perihelion, GM (cid:12) is the Sun’sgravitational parameter, R (cid:12) is the radius of the Sun, c is the speed of light, and a , e , and i are the semi-majoraxis, eccentricity, and orbital inclination (with respectto the solar equator) of a planetary body, respectively.Because both GR and solar oblateness affect perihelionprecession, estimates of β and J (cid:12) are highly correlatedand it is desirable to track a variety of solar systembodies with a range of a , e , i values to disentangle thetwo effects.Our selection of target asteroids follows the method ofMargot (2003). We select asteroids with both large per-ihelion shift values and favorable observing conditionswith radar (Table 1 and Figure 1). This sample of as-teroid orbits includes a wide range of semi-major axes,eccentricities, and inclinations, which are advantageouswhen simultaneously solving for β and J (cid:12) . The pre-dicted rates of perihelion advance, ˙ δω , shown in Figure1 and Table 1 were computed assuming γ = β = 1 and J (cid:12) = 2 . × − . METHODS rospects of determination of GR parameter β with asteroid radar astronomy Table 1.
Selected asteroids and orbital elements: Semima-jor Axis ( a ), Eccentricity ( e ), and Inclination with Respectto the Ecliptic ( i ec ) and Sun’s equator ( i eq ). Target a (au) e i ec (deg) i eq (deg) ˙ δω ( (cid:48)(cid:48) cy − )1566 Icarus 1.078 0.827 22.9 15.8 10.11998 TU3 0.787 0.484 5.41 3.41 9.111999 KW4 0.642 0.688 38.9 46.0 22.11999 MN 0.674 0.665 2.02 5.25 18.52000 BD19 0.876 0.895 25.7 28.0 26.92000 EE14 0.662 0.533 26.5 26.1 15.02001 YE4 0.677 0.541 4.82 11.0 14.42004 KH17 0.712 0.499 22.1 14.9 12.02006 CJ 0.676 0.755 10.3 16.1 23.7 Note.
The predicted rate of perihelion advance in arcseccentury − ( (cid:48)(cid:48) cy − ), ˙ δω , was computed using Equation (1). We first determined nominal trajectories for asteroidsin our sample with astrometric (i.e., positional) data,both optical and radar (Table 2). The process involvedthree steps: (1) numerical integration of each asteroid’sorbit and calculation of partial derivatives of the equa-tions of motion with respect to the solve-for parameters(i.e., the six components of the state vectors), (2) eval-uation of simulated optical and radar observables andcomputation of their partial derivatives with respect tothe solve-for parameters, and (3) least-squares adjust-ments to the solve-for parameters.We used the Mission Operations and NavigationToolkit Environment (MONTE) software (Evans et al.2016, MONTE v124) for orbit determination and pa-rameter estimation. MONTE is an astrodynamics com-puting platform developed by NASA’s Jet PropulsionLaboratory (JPL). MONTE is used for spacecraft nav-igation and trajectory design. MONTE has also beenused for a variety of scientific purposes, including grav-ity analysis (Verma & Margot 2016) and ephemerisgeneration (Greenberg et al. 2017).3.1.
Dynamical model
MONTE uses a variable-step Adams-Bashforthmethod to numerically integrate the equations of motionand corresponding partial derivatives. Our dynamicalmodel includes gravitational forces from the Sun, 8 plan-ets, and 21 minor planets with well-determined masses(Konopliv et al. 2011), general relativistic effects, andperturbations due to the oblateness of the Sun.In addition to these forces, we have also modeled thenongravitational Yarkovsky orbital drift. Perihelion ad-vance due to GR and solar oblateness does not affect the value of the semi-major axis, but Yarkovsky driftdoes. This nongravitational effect has been shown toaffect the semi-major axes of small bodies due to theanisotropic re-emission of absorbed sunlight (e.g., Bot-tke et al. 2006). The change in semi-major axis withtime due to Yarkovsky orbital drift, (cid:104) da/dt (cid:105) , was esti-mated for all target asteroids with the method of Green-berg et al. (2017). The values ranged in amplitude be-tween 4 and 50 au/My, which is plausible for kilometer-sized bodies. Only one target (1566 Icarus) is commonbetween our target list and the 42 Yarkovsky detectionsof Nugent et al. (2012), and only one target (1999 MN)is common between our target list and the 21 Yarkovskydetections of Farnocchia et al. (2013). In both cases, ourYarkovsky drift estimates are consistent with and betterconstrained than prior work.To initialize the integration process, we used a prioristate vectors extracted from the Minor Planet Center(MPC) database (Minor Planet Center 2017).3.2.
Existing optical and radar astrometry
We used both optical and radar astrometry to de-termine the nominal trajectory of each asteroid. Op-tical measurements provide positional information onthe plane of the sky. They are typically expressed asright ascension (R.A.) and declination (decl.) in theequatorial frame of epoch J2000.0. We downloaded op-tical astrometry from the MPC (Minor Planet Center2017). We debiased optical astrometry and assigneddata weights according to the algorithm recommendedby Farnocchia et al. (2015).Radar astrometry consists of round-trip light time,a measurement that can provide the asteroid − observerdistance, and Doppler shift, a measurement that canprovide the line-of-sight velocity of the asteroid with re-spect to the observer. Radar measurements have frac-tional uncertainties as small as 10 − . The addition ofradar astrometry can decrease orbital element uncer-tainties by orders of magnitude compared to an optical-only orbit solution (Ostro & Giorgini 2004). However,the number of radar measurements is typically smallcompared to the number of optical observations (Table2).We processed a total of 12,102 optical measurements(R.A. and decl. pairs obtained at 6051 epochs), as wellas 56 range and 17 Doppler measurements that havebeen published.3.3. Orbit determination for nominal trajectories
In order to compute nominal asteroid trajectories, wecomputed the expected values of the observables andtheir partial derivatives with respect to initial state vec-tors. We calculated weighted residuals by subtracting
Verma et al. 2017
Semi-major Axis [au] E cc e n t r i c i t y " / c y " / c y " / c y " / c y Semi-major Axis [au] I n c li n a t i o n [ d e g ] Icarus1998 TU3 1999 KW41999 MN 2000 BD192000 EE14 2001 YE42004 KH17 2006 CJ
Figure 1.
Distribution of asteroid orbital elements for asteroids in our sample. The corresponding rates of perihelion shift,predicted with Equation (1), are shown as contour lines. computed measurements ( C ) from observed measure-ments ( O ) and dividing the result by the correspond-ing observational uncertainty ( σ ). We adjusted initialstate vectors with an iterative least-squares techniquesthat minimized the sum of squares of weighted residuals.Because there are 9 targets and 6 orbital elements perasteroid in the nominal situation ( γ = 1 , β = 1 , J (cid:12) =2 . × − ), we adjusted a total of 54 parameters.We defined outliers as measurements with weightedresiduals in excess of three. We identified and re-jected 127 epochs with outliers in the optical astrome-try. There were no outliers in the radar astrometry. Weobtained a measure of the quality of the fit at each itera-tion by computing the dimensionless rms of the weightedresiduals: RMS = (cid:118)(cid:117)(cid:117)(cid:116) N N (cid:88) i =1 (cid:18) O i − C i σ i (cid:19) , (2)where N is the number of observations, O i is the i th ob-servation, C i is the i th computed measurement, and σ i is the observational uncertainty associated with the i thobservation. We stopped the iterative process when thechange in the RMS of the weighted residuals betweentwo successive iterations was less than 0.01%. RMSresiduals smaller than one indicate solutions that pro-vide good fits to the observations (Table 2).3.4. Anticipated radar astrometry
The objectives of this study are to evaluate the preci-sion with which PPN parameter β and solar quadrupolemoment J (cid:12) can be determined from orbital fits con-strained by existing and anticipated optical and radarastrometry. To quantify the effect of anticipated radar astrometry on the determination of these parameters, wesimulated all existing optical and radar astrometry (Ta-ble 2) and a number of anticipated Arecibo Observatoryrange measurements (Table 3) with the nominal aster-oid trajectories described above. We did not attemptto simulate the effect of additional optical astrometry,which is expected to improve the overall quality of thefits, albeit not as powerfully as radar astrometry (Ostro& Giorgini 2004).To supplement the published astrometry with realis-tic anticipated values, we used the epochs of closest ap-proach to Earth when the asteroids are detectable withthe Arecibo radar (Table 3). On the basis of prior expe-rience, we assumed that two to four independent datapoints would be collected at each future apparition. Forapparitions in the past (identified in bold in Table 3),we used the number of data points that were actuallyobtained. In total, we simulated 61 independent rangemeasurements in addition to the 56 published values.For each realization in our simulations, we added noiseto the observations by randomly drawing from a Gaus-sian distribution with zero mean and standard deviationequal to the observational uncertainty. Observationaluncertainties for observations in the future were assignedaccording to signal-to-noise ratio (S/N) and experience,with values ranging between 30 and 900 m. Uncertain-ties for observations in the past mirrored the actual mea-surement uncertainties adopted by the observer for thesedata points.3.5. Orbit determination with estimation of β and J (cid:12) We assigned solve-for parameters to one of two cate-gories: local and global. Local parameters are specific toeach asteroid, i.e., the 6 orbital elements or initial state rospects of determination of GR parameter β with asteroid radar astronomy Table 2.
Selected asteroids and corresponding observations: Observational Interval, Number of Optical Pairs (R.A. and Decl.)of Observations, and Number of Published Range and Doppler Observations.
Target Observational interval N opt N rng N dop RMS opt
RMS rng
RMS dop − − ... ... ... ... − ... ... ... ... − ... ... − ... ... ... ... − ... ... ... ... − − ... ... − Note.
The last three columns provides the post-fit root-mean-square of weighted residuals.
Table 3.
Selected asteroids and simulated observations:Years of Close Earth Approaches (yyyy − Target Year of close approach N range Uncertainties (m)1998 TU3 , 19 5 75–9001999 KW4 , 17, 18, 19, 20 12 40–3001999 MN , , , 20 10 300–3752000 EE14 , , 21, 22 11 300–6002001 YE4 , , 21 10 30–6002004 KH17 , , 22 9 60–300 Note.
Years highlighted in bold correspond to epochs forwhich data have already been collected. The nextdetectable approach of 1566 Icarus is not until 2024. vector (total of 9 × β and J (cid:12) .We jointly solved for these 56 parameters.We used two independent approaches to evaluate theprecision in the determination of global parameters β and J (cid:12) . First, we used a traditional covariance analysis(Section 4.1) as described in Bierman (1977). Second,we performed Monte Carlo simulations (Section 4.2) toverify the results of the covariance analysis. RESULTS4.1.
Covariance analysis
A covariance analysis is a powerful technique that canbe used to evaluate the precision of solve-for parameters.First, simulated, noise-free measurements and their par-tial derivatives are computed on the basis of nominal tra- jectories. A least-squares estimation is then performed,where the estimates logically converge on the nominalvalues. In the process, the associated covariance matrixis produced. The expected precision of the estimatedparameters is then inferred by examining the covariancematrix. The square roots of the diagonal elements pro-vide the one-standard-deviation formal uncertainties.After global fits of 56 parameters, we obtained thefollowing formal uncertainties: σ β = 5 . × − , (3) σ J (cid:12) = 2 . × − , (4)with a correlation coefficient of -0.72. The parametersremain correlated because both GR and solar oblatenesscontribute to perihelion precession. However, the rangeof asteroid orbital parameters (Table 1) helps reducethe correlation coefficient. Consideration of the Lense-Thirring effect for the Sun increases our σ β and σ J (cid:12) estimates by 0.2% and 4%, respectively.The expected formal uncertainty on J (cid:12) with di-rect dynamical measurement of asteroids is 2.7 timesthe uncertainty based on fits to helioseismology data(Antia et al. 2008). For β , the expected formal uncer-tainty is about twice the uncertainty obtained with pre-MESSENGER planetary ephemerides (Konopliv et al.2011), ∼ ∼
14 times theuncertainty obtained with MESSENGER range data(Park et al. 2017) The formal uncertainties scale linearlywith the uncertainties assigned to the measurements. Itis often the case that radar observers assign conservativeuncertainties, as evidenced by RMS residuals or reducedchi-square metrics that are almost always smaller thanunity and most often < . Verma et al. 2017 anticipate that the actual precision may be improved bya factor of ∼
3, and the dynamical determination of J (cid:12) may be as precise as the helioseismology determination.In order to investigate the benefit of future obser-vations, we also performed covariance analyses underthe assumption that observations would stop at theend of 2017, 2019, or 2021, as opposed to 2022 inour nominal scenario. The results were σ β, =9 . × − , σ β, = 7 . × − , σ β, = 7 . × − and σ J (cid:12) , = 1 . × − , σ J (cid:12) , = 4 . × − , σ J (cid:12) , = 3 . × − .4.2. Monte Carlo simulations
More robust results can be obtained by perform-ing end-to-end simulations that approximate the actualmeasurement and estimation process. In these analyses,integration of the trajectories and estimation of the pa-rameters are conducted as described in Section 3 withtwo variations. First, we chose initial values of the solve-for parameters that are not identical to their nominalvalues. For instance, the initial positions and velocitiesof all asteroids were changed by 10 km and 0.1 ms − ineach direction, respectively. Likewise, initial values for β and J (cid:12) were changed by 4 × − and 5 × − , whichis approximately five times the uncertainty of recent es-timates. Second, we polluted the simulated measure-ments with independent noise realizations as describedin Section 3.We performed 500 Monte Carlo simulations. Afterconvergence of the least-squares estimation, we com-pared the estimated values of solve-for parameters withtheir nominal values, which produced error estimates.To arrive at an estimate of the uncertainties, we canfit Gaussian distributions to the histograms of error es-timates, or we can compute the covariance matrix, asfollows:cov( p i , p j ) = 1 N − N (cid:88) k =1 ( p ki − p ni )( p kj − p nj ) , (5)where N is the total number of simulations, p ni is thenominal value of the i th parameter ( β = 1, J (cid:12) =2 . × − ), and p ki is the estimated value of the i th parameter from the k th simulation of observations. Weused Equation (5) and estimated the formal uncertain-ties in the solve-for parameters by computing the square root of diagonal elements. We found σ β = 7 . × − , (6) σ J (cid:12) = 3 . × − , (7)with a correlation coefficient of -0.81. These values con-firm the covariance analysis results. CONCLUSIONSA modest observing campaign requiring 50 −
60 hoursof Arecibo telescope time over the next five years canprovide about 20 range measurements of asteroids whoseorbits exhibit large perihelion shift rates. The AreciboPlanetary Radar facility is required for these measure-ments because its sensitivity is ∼
20 times better thanthat of other radar systems (Naidu et al. 2016), allowingdetection of asteroids that are not detectable elsewhere.The Arecibo measurements will complement existingoptical and radar astrometry and enable joint orbital so-lutions with β and J (cid:12) as adjustable parameters. Inde-pendent, purely dynamical determinations of both pa-rameters are important because they place bounds ontheories of gravity and the interior structure the of Sun,respectively.Our simulation results likely under-estimated actualprecision for two reasons. First, we did not attemptto simulate the impact of future optical astrometry norimprovements to the accuracy of star catalogs. Bothof these effects will inevitably improve the quality ofthe orbital determinations. Second, we assumed, basedon historical evidence, that radar observers assign fairlyconservative uncertainties to their measurements, whichoften underestimate the precision of the measurementsby a factor of ∼ σ β ∼ × − , (8) σ J (cid:12) ∼ − . (9)ACKNOWLEDGMENTSA.K.V., J.L.M., and A.H.G. were supported in partby the NASA Planetary Astronomy program undergrant NNX12AG34G. J.L.M. and A.H.G. were sup-ported in part by NSF Planetary Astronomy programAST-0929830 and AST-1109772. This work was en-abled in part by the Mission Operations and NavigationToolkit Environment (MONTE). MONTE is developedat the Jet Propulsion Laboratory, which is operated byCaltech under contract with NASA. Software:
MONTE v124 (Evans et al. 2016) rospects of determination of GR parameter β with asteroid radar astronomy Antia, H. M., Chitre, S. M., & Gough, D. O. 2008, A&A,477, 657Bertotti, B., Iess, L., & Tortora, P. 2003, nature, 425, 374Bierman, G. J. 1977, Factorization Methods for DiscreteSequential Estimation (volume 128, 241 pages AcademicPress, New York, NY, 1977)Bottke, Jr., W. F., Vokrouhlick´y, D., Rubincam, D. P., &Nesvorn´y, D. 2006, Annual Review of Earth andPlanetary Sciences, 34, 157Evans, S., Taber, W., Drain, T., et al. 2016, in The 6thInternational Conference on Astrodynamics Tools andTechniques (ICATT), International Conference onAstrodynamics Tools and Techniques, Darmstadt,Germany. https://indico.esa.int/indico/event/111/session/30/contribution/177/material/paper/0.pdfhttps://indico.esa.int/indico/event/111/session/30/contribution/177/material/paper/0.pdf