Prospects of dynamical determination of General Relativity parameter beta and solar quadrupole moment J2 with asteroid radar astronomy
Ashok K. Verma Jean-Luc Margot Adam H. Greenberg
astro-ph.EP gr-qc physics.space-ph



















































Ashok Kumar Verma
ashokverma@ucla.edu

0000-0002-7209-0004]Ashok K. Verma
Department of Earth, Planetary, and Space Sciences, University of California, Los Angeles, CA 90095, USA

0000-0001-9798-1797]Jean-Luc Margot
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

0000-0001-8834-9423]Adam H. Greenberg
Department of Physics and Astronomy, University of California, Los Angeles, CA 90095, USA
















We evaluated the prospects of quantifying the parameterized
post-Newtonian parameter β and solar quadrupole moment J_2⊙
with observations of near-Earth asteroids with large orbital
precession rates (9 to 27 arcsec century^-1).  We considered existing
optical and radar astrometry, as well as radar astrometry that can
realistically be obtained with the Arecibo planetary radar in the next
five years.  Our sensitivity calculations relied on a traditional
covariance analysis and Monte Carlo simulations.  We found that
independent estimates of β and J_2⊙ can be obtained with
precisions of 6×10^-4 and 3×10^-8,
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×10^-4 and 10^-8, respectively.  A purely dynamical
determination of solar oblateness with asteroid radar astronomy may
therefore rival the helioseismology determination.






















§ INTRODUCTION

 
The parameterized post-Newtonian (PPN) formalism is a useful framework
for testing metric theories of gravity <cit.>.  It consists of
10 dimensionless parameters that describe the general properties of
the metric.  In general relativity (GR), only 2 of the 10
parameters are non-zero. They are known as the
Eddington-Robertson-Schiff parameters γ and β.  γ
represents the amount of curvature produced by a unit mass, and
β represents the amount of nonlinearity in the superposition law
for gravity.

Several techniques have been used to place observational bounds on
these parameters <cit.>, including observations of the bending
and delay of light by spacecraft tracking <cit.> or
Very Long Baseline Interferometry <cit.>, and
fitting of ephemerides to observations of planetary positions
<cit.>.

In GR, γ and β are equal to one.  Doppler tracking of the
Cassini spacecraft has shown that γ does not differ from one by
more than 2 × 10^-5 <cit.>.  Ephemeris-based
studies prior to 2009 indicated that β - 1 does not differ from
zero by more than 10^-4 <cit.>.  More
recently, the availability of precise ranging data from the MESSENGER
Mercury orbiter <cit.> enabled improved
estimates <cit.>.  Here, we evaluate the
prospect of asteroid orbit precession measurements to place more
stringent bounds on β.  We consider Earth-based radar
observations of near-Earth asteroids with perihelion shifts larger
than 10 arcsec century^-1.

Orbital precession can also be caused by the nonuniformity of the
gravity field that results from the oblate shape of the Sun.  The
solar oblateness is characterized by the solar quadrupole moment,
J_2⊙ <cit.>.  Simultaneous estimation of
β and J_2⊙ requires that the precessional effects due to
GR and to the Sun's oblateness be disentangled.  Fortunately, GR is a
purely central effect, whereas the oblateness-induced precession has an
inclination dependence.  The two effects also have a different
distance dependence <cit.>.  As a result, observations of a
small sample of near-Earth asteroids with a variety of semi-major axes
and inclinations (Table <ref>) can in principle be used to
estimate β and J_2⊙ <cit.>.

Current estimates of the solar quadrupole moment are typically derived
on the basis of interior models of the Sun constrained by
helioseismology data <cit.>.  The current
best value from the helioseismology literature is J_2⊙ = (2.2
± 0.1) × 10^-7 <cit.>.  Dynamical estimates that do
not rely on fits to helioseismology data yield similar values of
J_2⊙ = 2.3± 0.25 × 10^-7 <cit.> and
J_2⊙ = 2.25± 0.09 × 10^-7 <cit.>.
High-precision dynamical estimates are important to validate our
understanding of the interior structure of the Sun.

Our simulations of the determination of β and J_2⊙ using a
variety of asteroid orbits suggest that independent values of β
and J_2⊙ can be obtained with satisfactory precision: with
the traditionally conservative assignment of radar uncertainties,
β can be constrained at the 6×10^-4 level and
J_2⊙ can be constrained at the 3×10^-8 level.  With
uncertainties that more closely reflect measurement errors, this
precision may be improved by a factor of ∼3. (Section
<ref>).

The outline of this paper is as follows.  In Section <ref>, we
describe our choice of target asteroids.  In Section <ref>, we
discuss the estimation of asteroid orbits with optical and radar
measurements.  Our dynamical model and data reduction procedures are
described in Section <ref> and <ref>, respectively.
Orbit determination results are presented in Section <ref>.
Simulations of the determination of β and J_2⊙ are described
in Section <ref>.





§ TARGET ASTEROIDS


The per-orbit secular advance in the angular position of the perihelion is given by
<cit.>

    δω = 6π GM_⊙/a(1-e^2)c^2[(2-β+2γ)/3] 
                       + 6π/2R^2_⊙(1-3/2 sin^2i)/a^2(1-e^2)^2J_2⊙,

where ω is the argument of perihelion, GM_⊙ is the Sun's
gravitational parameter, R_⊙ is the radius of the Sun, c is
the speed of light, and a, e, and i are the semi-major axis,
eccentricity, and orbital inclination (with respect to the solar
equator) of a planetary body, respectively.  Because both GR and solar
oblateness affect perihelion precession, estimates of β and
J_2⊙ are highly correlated and it is desirable to track a
variety of solar system bodies with a range of a, e, i values to
disentangle the two effects.

Our selection of target asteroids follows the method of
<cit.>.  We select asteroids with both large perihelion
shift values and favorable observing conditions with radar (Table
<ref> and Figure <ref>).  This sample of asteroid orbits
includes a wide range of semi-major axes, eccentricities, and
inclinations, which are advantageous when simultaneously solving for
β and J_2⊙.  The predicted rates of perihelion advance,
δ̇ω̇, shown in Figure <ref> and Table
<ref> were computed assuming γ = β = 1 and
J_2⊙=2.2 × 10^-7.








§ METHODS



We first determined nominal trajectories for asteroids in our sample
with astrometric (i.e., positional) data, both optical and radar
(Table <ref>). The process involved three steps: (1) numerical
integration of each asteroid's orbit and calculation of partial
derivatives of the equations of motion with respect to the solve-for
parameters (i.e., the six components of the state vectors), (2) evaluation of 
simulated optical and radar observables
and computation of their partial derivatives with respect to the
solve-for parameters, and (3) least-squares adjustments to the solve-for 
parameters.

We used the Mission Operations and Navigation Toolkit Environment
(MONTE) software <cit.>
for orbit determination and parameter estimation.
MONTE is an astrodynamics computing platform developed by
NASA's Jet Propulsion Laboratory (JPL). MONTE is used for spacecraft
navigation and trajectory design.  MONTE has also been used for a
variety of scientific purposes, including gravity analysis
<cit.> and ephemeris generation <cit.>.






 §.§ Dynamical model



MONTE uses a variable-step Adams-Bashforth method to numerically
integrate the equations of motion and corresponding partial
derivatives. Our dynamical model includes gravitational forces from
the Sun, 8 planets, and 21 minor planets with well-determined
masses <cit.>, general relativistic effects,
and perturbations due to the oblateness of the Sun.

In addition to these forces, we have also modeled the
nongravitational Yarkovsky orbital drift.  Perihelion advance due to
GR and solar oblateness does not affect the value of the semi-major
axis, but Yarkovsky drift does.  This nongravitational effect has
been shown to affect the semi-major axes of small bodies due to the
anisotropic re-emission of absorbed
sunlight <cit.>. The change in semi-major axis
with time due to Yarkovsky orbital drift, ⟨ da/dt ⟩, was
estimated for all target asteroids with the method of
<cit.>.  The values ranged in amplitude between 4 and 50
au/My, which is plausible for kilometer-sized bodies.  Only one target (1566
Icarus) is common between our target list and the 42 Yarkovsky detections
of <cit.>, and only one target (1999 MN) is common between our
target list and the 21 Yarkovsky detections of <cit.>.  In
both cases, our Yarkovsky drift estimates are consistent with and
better constrained than prior work.

To initialize the integration process, we used a priori state vectors
extracted from the Minor Planet Center (MPC) database <cit.>.



 §.§ Existing optical and radar astrometry


We used both optical and radar astrometry to determine the nominal
trajectory of each asteroid. Optical measurements provide positional
information on the plane of the sky.  They are typically expressed as
right ascension (R.A.) and declination (decl.) in the equatorial frame of
epoch J2000.0.  We downloaded optical astrometry from the MPC
<cit.>.  We debiased optical astrometry and assigned data weights
according to the algorithm recommended by <cit.>.

Radar astrometry consists of round-trip light time, a
measurement that can provide the asteroid-observer distance, and
Doppler shift, a measurement that can provide the line-of-sight
velocity of the asteroid with respect to the observer.
Radar measurements have fractional uncertainties as small as
10^-8.  The addition of radar astrometry can decrease orbital
element uncertainties by orders of magnitude compared to an
optical-only orbit solution <cit.>.  However, the number of
radar measurements is typically small compared to the number of
optical observations (Table <ref>).

We processed a total of 12,102 optical measurements (R.A. and decl. pairs
obtained at 6051 epochs), as well as 56 range and 17 Doppler
measurements that have been published.



 §.§ Orbit determination for nominal trajectories


In order to compute nominal asteroid trajectories, we computed the
expected values of the observables and their partial derivatives with
respect to initial state vectors.  We calculated weighted residuals by
subtracting computed measurements (C) from observed measurements
(O) and dividing the result by the corresponding observational
uncertainty (σ).  We adjusted initial state vectors with an
iterative least-squares techniques that minimized the sum of squares
of weighted residuals. Because there are 9 targets and 6 orbital
elements per asteroid in the nominal situation (γ = 1, β =1, 
J_2⊙ = 2.2 × 10^-7), we adjusted a total of 54
parameters.


We defined outliers as measurements with weighted
residuals in excess of three.  We identified and rejected 127 epochs with
outliers in the optical astrometry.  There were no outliers in the
radar astrometry.  We obtained a measure of the quality of the fit at
each iteration by computing the dimensionless rms of the weighted
residuals:

    RMS = √(1/N∑_i=1^N( O_i-C_i/σ_i)^2),

where N is the number of observations, O_i is the ith
observation, C_i is the ith computed measurement, and
σ_i is the observational uncertainty
associated with the ith observation.  We stopped the iterative
process when the change in the RMS of the weighted residuals between
two successive iterations was less than 0.01%.
RMS residuals smaller than one indicate solutions that provide good
fits to the observations (Table <ref>).






 §.§ Anticipated radar astrometry


The objectives of this study are to evaluate the precision with which
PPN parameter β and solar quadrupole moment J_2⊙ can be
determined from orbital fits constrained by existing and anticipated
optical and radar astrometry. To quantify the effect of anticipated
radar astrometry on the determination of these parameters, 
we simulated all existing optical and radar astrometry (Table <ref>) and
a number of anticipated Arecibo Observatory range measurements (Table <ref>)
with the nominal asteroid trajectories described above.
We did not attempt to simulate the effect of additional optical
astrometry, which is expected to improve the overall quality of the
fits, albeit not as powerfully as radar astrometry <cit.>.


To supplement the published astrometry with realistic anticipated
values, we used the epochs of closest approach to Earth when the
asteroids are detectable with the Arecibo radar (Table <ref>).
On the basis of prior experience, we assumed that two to four independent
data points would be collected at each future apparition.  For
apparitions in the past (identified in bold in Table <ref>),
we used the number of data points that were actually obtained.  In
total, we simulated 61 independent range measurements in addition to
the 56 published values.  For each realization in our simulations, we
added noise to the observations by randomly drawing from a Gaussian
distribution with zero mean and standard deviation equal to the
observational uncertainty.  Observational uncertainties for
observations in the future were assigned according to signal-to-noise ratio
(S/N) and experience, with values ranging between 30 and 900 m.
Uncertainties for observations in the past mirrored the actual
measurement uncertainties adopted by the observer for these data
points.









 §.§ Orbit determination with estimation of β and J_2⊙


We assigned solve-for parameters to one of two categories: local and
global.  Local parameters are specific to each asteroid, i.e., the 6
orbital elements or initial state vector (total of 9 × 6 = 54
parameters), whereas global parameters are common to all asteroids,
i.e., β and J_2⊙.  We jointly solved for these 56
parameters.

We used two independent approaches to evaluate the precision in the
determination of global parameters β and J_2⊙.  First, we
used a traditional covariance analysis (Section <ref>) as
described in <cit.>.
Second, we performed Monte Carlo simulations (Section <ref>) to
verify the results of the covariance analysis.





§ RESULTS






 §.§ Covariance analysis


A covariance analysis is a powerful technique that can be used to
evaluate the precision of solve-for parameters.  First, simulated,
noise-free measurements and their partial derivatives are computed on
the basis of nominal trajectories.  A least-squares estimation is then
performed, where the estimates logically converge on the nominal
values.  In the process, the associated covariance matrix is produced.
The expected precision of the estimated parameters is then inferred by
examining the covariance matrix. The square roots of the diagonal
elements provide the one-standard-deviation formal uncertainties.

After global  fits of 56 parameters, we obtained the following formal uncertainties:

    σ_β = 5.6 × 10^-4,


    σ_J_2⊙ =  2.7 × 10^-8,

with a correlation coefficient of -0.72.  The parameters remain
correlated because both GR and solar oblateness contribute to
perihelion precession.  However, the range of asteroid orbital
parameters (Table <ref>) helps reduce the correlation
coefficient. Consideration of the Lense-Thirring effect for the Sun 
increases our σ_β and σ_J_2⊙ estimates by 0.2% and 4%, respectively.


The expected formal uncertainty on J_2⊙ with direct dynamical
measurement of asteroids is 2.7 times the uncertainty based on
 fits to helioseismology data <cit.>.
For β, the expected formal uncertainty is
about twice the uncertainty obtained with pre-MESSENGER planetary
ephemerides <cit.>, ∼7 times the uncertainty
obtained with post-MESSENGER planetary
ephemerides <cit.>, and ∼14 times the uncertainty
obtained with MESSENGER range data <cit.>
The formal uncertainties scale linearly with the uncertainties
assigned to the measurements.  It is often the case that radar
observers assign conservative uncertainties, as evidenced by RMS
residuals or reduced chi-square metrics that are almost always smaller
than unity and most often <0.3 (Table <ref>).  Therefore, we
anticipate that the actual precision may be improved by a factor of
∼3, and the dynamical determination of J_2⊙ may be as
precise as the helioseismology determination.

In order to investigate the benefit of future observations, we also
performed covariance analyses under the assumption that observations
would stop at the end of 2017, 2019, or 2021, as opposed to 2022 in
our nominal scenario.  The results were
σ_β,2017 = 9.6  × 10^-4,
σ_β,2019 = 7.6  × 10^-4,
σ_β,2021 = 7.5  × 10^-4 and
σ_J_2⊙,2017 = 1.9  × 10^-7,
σ_J_2⊙,2019 = 4.2  × 10^-8,
σ_J_2⊙,2021 = 3.8  × 10^-8.




 §.§ Monte Carlo simulations



More robust results can be obtained by performing end-to-end
simulations that approximate the actual measurement and estimation
process.  In these analyses, integration of the trajectories and
estimation of the parameters are conducted as described in Section
<ref> with two variations.  First, we chose initial values of
the solve-for parameters that are not identical to their nominal
values.  For instance, the initial positions and velocities of all
asteroids were changed by 10 km and 0.1 ms^-1 in each direction,
respectively.  Likewise, initial values for β and J_2⊙
were changed  by 4 × 10^-4 and 5 × 10^-8, which is
approximately five times the uncertainty of recent estimates.
Second, we polluted the simulated measurements with independent noise
realizations as described in Section <ref>.

We performed 500 Monte Carlo simulations.  
After convergence of the
least-squares estimation, we compared the estimated values of
solve-for parameters with their nominal values, which produced error
estimates.
To arrive at an estimate of the uncertainties, we can fit Gaussian
distributions to the histograms of error estimates, or we can compute
the covariance matrix, as follows:

    cov(p_i, p_j) = 1/N-1∑_k=1^N(p_i^k - p_i^n)(p_j^k - p_j^n),

where N is the total number of simulations, p_i^n is the nominal
value of the ith parameter (β=1, J_2⊙ = 2.2
× 10^-7), and p_i^k is the estimated value of
the ith parameter from the kth simulation of observations.
We used Equation (<ref>) and estimated the formal uncertainties in
the solve-for parameters by computing the square root of diagonal
elements.  We found

    σ_β = 7.4 × 10^-4,


    σ_J_2⊙ =  3.7 × 10^-8,

with a correlation coefficient of -0.81.  These values confirm the
covariance analysis results.




§ CONCLUSIONS


A modest observing campaign requiring 50-60 hours of Arecibo telescope
time over the next five years can provide about 20 range measurements
of asteroids whose orbits exhibit large perihelion shift rates.  The
Arecibo Planetary Radar facility is required for these measurements
because its sensitivity is ∼20 times better than that of other
radar systems <cit.>, allowing detection of asteroids that are
not detectable elsewhere.

The Arecibo measurements will complement existing optical and radar
astrometry and enable joint orbital solutions with β and
J_2⊙ as adjustable parameters.  Independent, purely dynamical
determinations of both parameters are important because they place
bounds on theories of gravity and the interior structure the of Sun,
respectively.

Our simulation results likely under-estimated actual precision for
two reasons. First, we did not attempt to simulate the impact of
future optical astrometry nor improvements to the accuracy of star
catalogs.  Both of these effects will inevitably improve the quality
of the orbital determinations.  Second, we assumed, based on
historical evidence, that radar observers assign fairly conservative
uncertainties to their measurements, which often underestimate the
precision of the measurements by a factor of ∼3 (Table <ref>).  As a result, we
anticipate that the uncertainties of the final estimates will be close
to

    σ_β∼ 2 × 10^-4,


    σ_J_2⊙∼ 10^-8.








§ ACKNOWLEDGMENTS


A.K.V., J.L.M., and A.H.G. were supported in part by the NASA Planetary
Astronomy program under grant NNX12AG34G.  J.L.M. and A.H.G. were supported
in part by NSF Planetary Astronomy program AST-0929830 and
AST-1109772.  This work was enabled in part by the Mission Operations
and Navigation Toolkit Environment (MONTE).  MONTE is developed at the
Jet Propulsion Laboratory, which is operated by Caltech under contract
with NASA.

MONTE v124 <cit.>


















































