DYNAMICAL INFERENCE FROM A KINEMATIC SNAPSHOT: THE FORCE LAW IN THE SOLAR SYSTEM

The Gaia Satellite (Perryman et al. 2001) will measure positions and velocities for millions to billions of stars at varying precision. One of the principal goals of this mission is to provide the data necessary to infer the dynamical state of the Milky Way. However, there are issues in principle with inference of dynamics from a snapshot or instantaneous set of configuration and velocity measurements: the instantaneous positions and velocities have no necessary relationship with the gravitational potential or accelerations. Indeed, despite considerable literature (for example, Oort 1932; Schwarzschild 1979; Little & Tremaine 1987; Kaasalainen & Binney 1994; Johnston et al. 1999; Beloborodov & Levin 2004) there is no methodology for performing this inference that naturally handles all of the issues, including finite and non-trivial observational uncertainties or noise, missing data, non-steady aspects of the mass distribution, and the (incredibly likely) possibility that the potential is not (simply) integrable. Robust inference may not even be possible if the Milky Way has significant time dependence or is strongly chaotic or is far from showing any simple symmetries (such as axisymmetry).

Certainly there is no hope for dynamical inference on the massive scale required for the Gaia data set if we cannot perform it on much simpler, much more symmetrical, much older (in a dynamical sense), and much smaller (in a data sense) systems. In what follows, we take one of the simplest possible systems—the solar system—and the smallest possible data set—the positions and velocities of the major planets at a moment in time—and perform a complete dynamical inference. For a test system we could also have chosen the black hole at the Galactic center, where similar considerations apply. However, this system has additional issues with "missing data" because not all six phase-space coordinates are directly measurable for all stars. The solar system truly is the simplest problem in this class.

Our inferential starting point is orbital roulette (Beloborodov & Levin 2004). In this previous work, it was assumed that the orbital angles (in the action–angle formalism) are uniformly distributed. Thus, dynamical model parameters that correspond to orbital angles that suspiciously lack diversity are rejected. More specifically, the method considers a large range of possible dynamical model parameters, and computes orbital angles at every value of the parameters and a distribution statistic on those orbital angles. If the distribution statistic is designed to monotonically increase with the diversity of the angles or the flatness of the angle distribution, the "best-fit" dynamical model parameters are those that optimize the distribution statistic. This kind of approach is inherently frequentist: on many applications of these procedures to independent data sets a confidence interval captures the true dynamical parameters on a guaranteed fraction of trials. However, for any given data set the confidence interval produced might not represent a credible set of parameters.

In what follows we cast the dynamical parameters estimation as a probabilistic inference problem (a "Bayesian" approach). We adopt the same core assumption as the roulette problem: we assume that the orbital angles are uniformly distributed. However, in this framework, we must also construct prior probability distribution functions for the dynamical model parameters and conditional probabilities of the data given the model to encode the assumptions of a long-lived and bound solar system. These prior and conditional probability distributions and the data create posterior probability distribution functions for all the parameters. We use this new method to infer the gravitational force law (radial dependence and amplitude). What is new here, in the context of solar system dynamics, is that we perform this inference with only a snapshot of the kinematic state, that is, with only the positions and velocities of the planets at a single instant of time.

Of course the kinematic snapshot we employ is, in fact, a set of initial conditions for a solar system integration (Giorgini et al. 1996). These initial conditions were determined not by a single measurement at a single epoch, but are in fact the result of an optimization of a solar system integration to observed planetary positions over many decades. In the context of this paper—a demonstration of a method—it is best to think of these "data" as "simulated data" useful for testing the method. They just happen to be data that have been simulated by the analog computer we know as the solar system.

What is new here, in the context of dynamical inference, is that our method is fully probabilistic or Bayesian. This is important for future problems, such as the Gaia problem, or for inferring the mass of the black hole at the center of the Milky Way, because in these real data analysis problems, the data points come with non-trivial and highly correlated observational uncertainties, and because entire dimensions of phase space are missing or unobserved. At the Galactic center, we do not know the radial distance to anything accurately, and in the Gaia data set many of the radial velocities will not be measured. The Bayesian framework handles these real-world data issues very naturally, although in fact they are not important in the test problem we solve here. Aside from these issues of principle, it is also the case, as we will show, that the Bayesian method performs extremely well.

Of course, a lot is known about the gravitational force law in the solar system, so we do not expect, at the outset, to be surprised by our results. The first force-law inference in the solar system (Newton 1687, and also work by contemporaries, particularly Hooke, who may have priority) made use of full orbit shape determination (Kepler 1609). In this sense, Newton's problem—find the force law (from among a small, discrete group of possible force laws) that leads to elliptical orbits with the Sun at one focus—was much easier than the problem we have set for ourselves. Of course, along the way, Newton had to develop for the first time the general principles of kinematics and dynamics!

We are going to assume spherical symmetry of the solar system's force law and gravitational potential, although nothing in the general inference formalism that follows will require this. Consider a radial force law (really acceleration law) of the form

$$\begin{equation} {{\vec {\bm a}}}= -A\,\left[\frac{r}{r_0}\right]^{-\alpha }\,{\hat{r}}, \end{equation} \tag{ 1 }$$

where A is an amplitude, r is the distance from the Sun, r₀ is a distance scale (in this case we will use r₀ = 1 AU so that A can be thought of as the acceleration at Earth's orbit), α parameterizes the radial dependence, and ${\hat{r}}$ is the radial direction. In this model, the list of free parameters is

$$\begin{equation} \boldsymbol{{\omega }}\equiv \lbrace \ln A,\alpha \rbrace, \end{equation} \tag{ 2 }$$

where we have taken the logarithm of A because in inference problems, dimensioned parameters are usually best handled in the log (Jeffreys 1939; Sivia & Skilling 2006).

The potential u (potential energy per unit planet mass), radial effective potential u_eff, and binding energy per unit mass ≡ −E/m are

$$\begin{equation} u(r) = \frac{A\,r_0}{1-\alpha }\,\left[\frac{r}{r_0}\right]^{1-\alpha }, \end{equation} \tag{ 3 }$$

$$\begin{equation} u_{\mathrm{eff}}(r) = u(r) + \frac{j^2}{2\,r^2}, \end{equation} \tag{ 4 }$$

$$\begin{equation} \epsilon = - u_{\mathrm{eff}}- \frac{1}{2}\,v_r^2, \end{equation} \tag{ 5 }$$

where j² is the square of the magnitude of the planet's angular momentum per unit mass (or L²/m²), and v_r is the radial component of the velocity (the component of ${\vec{\bm v}}$ parallel to ${\vec{\bm x}}$). The perihelion and aphelion distances r_peri and r_ap are both found by setting = −u_eff. With these, we can define a radial asymmetry e as

$$\begin{equation} e \equiv \frac{r_{\mathrm{ap}}- r_{\mathrm{peri}}}{r_{\mathrm{ap}}+ r_{\mathrm{peri}}}, \end{equation} \tag{ 6 }$$

where we have called this "e" because in the Kepler–Newton α = 2 case it is the orbital eccentricity. One way of thinking of this radial asymmetry is that at any point in the space made up of the dynamical parameters $\boldsymbol{{\omega }}$ and the binding energy , the radial asymmetry e is a dimensionless description of the angular momentum magnitude.

Importantly for what follows, we can define a "radial angle" ϕ_r that increases linearly with time from perihelion passage through next perihelion passage. Any planet at radius r on an orbit with perihelion distance r_peri and aphelion distance r_ap can be assigned this angle ϕ_r by

$$\begin{equation} \phi _r\equiv \left\lbrace \begin{array}{@{}ll}\displaystyle \pi \,\frac{t(r)-t(r_{\mathrm{peri}})}{t(r_{\mathrm{ap}})-t(r_{\mathrm{peri}})} & \mbox{for $v_r > 0$} \\ [2.5ex] \displaystyle \pi +\pi \,\frac{t(r)-t(r_{\mathrm{ap}})}{t(r_{\mathrm{peri}})-t(r_{\mathrm{ap}})} & \mbox{for $v_r < 0$} \end{array}\right., \end{equation} \tag{ 7 }$$

where the first numerator is the time it takes to go from r_peri to r outbound, the first denominator is the time it takes to go from r_peri to r_ap outbound, the second numerator and denominator are the times inbound, and all time differences between two radii can be computed numerically for general values of α by integrating the inverse of the radial velocity between these radii. The first-order form of this integral has an integrable singularity at the perihelion and aphelion, which can be handled by an appropriate change of variables (e.g., Press et al. 2007). A planet observed at a set of random times spanning many orbits will be observed to have radial angles ϕ_r drawn from a flat distribution in the range 0 < ϕ_r < 2 π. This radial angle is one of the angles in the action–angle formulation of the system, which is integrable for the simple reason that it is spherically symmetric.

In what follows, we are going to use and compare several methods for inferring the force-law parameters $\boldsymbol{{\omega }}$ (the amplitude ln A and radial exponent α of the spherical force law) from an instantaneous snapshot of the positions and velocities of the eight major planets. This snapshot was taken from JPL's HORIZONS System which provides highly accurate ephemerides for solar system objects⁴ (Giorgini et al. 1996). It is an extrapolation (at the time of writing) to 2009 April 1.0, approximately 400 years after the important publication of Kepler (1609). This kinematic snapshot is given in Table 1.

Since this snapshot is obtained by integrating the positions and velocities of solar system bodies, the accuracy is limited by (1) the correctness of the dynamical model used, (2) the numerical integration of the equations of motion, and (3) the accuracy to which the initial conditions are known. It is generally believed that the dynamical model used is correct and complete, and that the numerical integration is sufficiently accurate. The main uncertainty in the ephemerides is then that due to the uncertainty in the initial conditions. The current set of initial conditions (DE405; Standish 1998) is a fit to a set of optical, radar, and VLBI observations as well as to a set of spacecraft range and Doppler points from various space missions. The uncertainties are the largest for the outer planets, since the data for these are almost entirely from optical observations (with the exception of Jupiter), and because Neptune has not been observed over a full orbit since the start of precise measurements. A comparison between the DE405 ephemerides and more recent observations shows that the positions of the inner planets are known to a fractional accuracy of approximately 10⁻⁸, while those of the outer planets are known to a fractional accuracy of 10⁻⁶–10⁻⁷ (Standish 2004). Uncertainties in the velocities are at the same fractional magnitude.

This kinematic snapshot is not, of course, a fair data set with which to perform the inference below, for the main reason that the "measured" kinematic state of the solar system is in fact the output of fitting observations with a dynamical model that assumes α = 2. For this reason, the data should be thought of as "idealized" or "simulated data" and the work must be considered a test of the method rather than a definitive inference.

The virial theorem relates the time averages 〈T〉 and 〈U〉 of a test particle's kinetic and potential energies through

$$\begin{equation} \left<T\right> = \frac{1-\alpha }{2}\,\left<U\right>, \end{equation} \tag{ 8 }$$

where α is the exponent in the radial force law. Given that a planet's potential energy is a function of the dynamical parameters $\boldsymbol{{\omega }}$, while its kinetic energy is not, the virial relation for each planet becomes a one-dimensional locus in the $\boldsymbol{{\omega }}$ space. Using kinetic and potential energies computed from the observations as a proxy for the time-averaged energies, these loci are shown in Figures 1 and 2. The fact that all eight lines cross near a single point in the space is encouraging that the system is virialized (as we expect) and that the inference will work. That the eight lines so nearly intersect in a single point is because, in fact, many of the planets are on circular orbits; the probabilistic method developed in this paper does not assume this, but it does retain the precision that is possible because of this situation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Also shown in Figures 1 and 2 is the region of parameter space in which one or more of the planets is unbound because T>U or, equivalently, < 0. In what follows, we will assign vanishing probability to regions of parameter space in which one or more planets is unbound. However, as we will see, this does not affect any of our conclusions. Also shown in Figure 1 is the region of parameter space in which one or more of the planets has r_peri < R_☉, where R_☉ is the radius of the Sun. This part of parameter space ought also to be excluded, although in practice this is not necessary for any of what follows.

In preparation for what follows, we pre-compute all of the planet radial angles ϕ_r,i—given their positions ${\vec{\bm x}}_i$ and velocities ${\vec{\bm v}}_i$—as a function of the dynamical parameters. These angles are shown in Figure 3.

Zoom In Zoom Out Reset image size

In orbital roulette (Beloborodov & Levin 2004), the idea is to compute, for each point $\boldsymbol{{\omega }}$ in parameter space, the N radial angles ϕ_r,i and analyze them statistically for being well mixed. In practice, this means applying a distribution test or multiple distribution tests to the angles and preferring parameters for which these tests are more consistent with a random or flat distribution of angles in the range 0 < ϕ_r < 2 π. Because each test provides one constraint, and in the case described here the parameter space is two dimensional, at least two qualitatively different tests are required to make localized constraints in $\boldsymbol{{\omega }}$. In addition to a test for a flat angle distribution we also apply a test for an angle–energy correlation.

About the simplest consistency test for the calculated angles is a test of the mean of the angles: Is the mean consistent with the expected mean for a uniformly distributed set of N planets? For this to perform well one must fold the angles of the inbound planets onto the interval [0, π], that is, disregard the information in the sign of the radial velocity; then, the expected mean of the angles is equal to π/2 (for details on how to test this assumption, see Beloborodov & Levin 2004). The fact that we have to perform this mapping indicates that the procedure is ad hoc. Indeed, for a uniform distribution on the circle there is no specific meaning to the perihelion and the aphelion, or to any two points, such that no real meaning can be attached to the mean of the angles between two arbitrary points on the circle.

Better, one can test the consistency of the full distribution function of the angles with a uniform distribution. This could again be done with the full [0, 2π] distributed angles or with the folded angles; the results will depend on this choice. Testing an observed distribution for consistency with an expected distribution often involves comparing cumulative distribution functions (Kolmogorov 1941). The Kolmogorov–Smirnov (K–S) test is the simplest in practice, since the distribution of the test statistic—the maximum difference between the cumulative distributions—can be approximated by an analytic function (Stephens 1970). The K–S test is by construction most sensitive to deviations near the median value; this rules out dynamical parameters at which the planets bunch up at perihelion or at aphelion, the situation in which about half are at perihelion and half at aphelion can easily dupe the test.

A statistic can be chosen that is sensitive to deviations at all values, such as the Anderson–Darling statistic (Anderson & Darling 1952). However, no approximate analytic description of the distribution of this statistic exists and in practice this distribution has to be obtained by Monte Carlo sampling (e.g., Beloborodov & Levin 2004). A statistic more appropriate to the problem at hand (although we are not primarily interested in a careful examination of the differences between different frequentist procedures) is Kuiper's statistic (Kuiper 1962). This statistic—the sum of the maximum distance of the observed cumulative distribution above and below the expected cumulative distribution—is invariant under periodic shifts of the angle and was specifically designed to test uniform distributions on the circle. The advantage over Anderson–Darling is that the asymptotic distribution of the Kuiper statistic is known (e.g., Press et al. 2007).

All of these tests for the uniformity of the distribution of the angles are shown in Figure 4. These tests can fail when for a certain combination of dynamical parameters some planets are near aphelion while other planets are near perihelion. This situation appears here in Figure 3, at α far from 2, where there are large regions in which the inner planets, especially Mercury and Venus, are near perihelion, while the outer planets, especially Uranus and Neptune, are near aphelion and vice versa. This prevents the mean angle and K–S tests from excluding those regions. The Kuiper test performs better.

Zoom In Zoom Out Reset image size

A second constraint (for the two-dimensional parameter space) comes from a second test. In regions of parameter space in which the inner planets are all near perihelion and the outer planets are all near aphelion, a significant correlation between the angles and the energies exists. This correlation is unphysical if the system is not being observed at any special time. A non-parametric test for the correlation is preferred here as the angle–energy correlation will not in general be linear. We perform a test of the angle–energy correlation using Kendall's τ (Kendall 1938). This is a rank test; it only considers the relative ordering of the angles and energies of different planets (for details on this test see Press et al. 2007). That this test is in a sense orthogonal to the tests of the uniformity of the angle distribution can be seen in Figure 4.

All of the frequentist tests permit acceptance or rejection of a dynamical model at a certain confidence level. Confidence intervals of 95% and 99% for all of the frequentist tests are shown in Figure 4. Also shown is the combination of the tests of the uniformity of the angle distribution with the test of the angle–energy correlation.

The frequentist procedures perform reasonably well in this simple problem because the number of dynamical parameters is small and the data have vanishing errors and no missing components. As the number of dynamical parameters increases the frequentist must find larger numbers of tests in order to break degeneracies. While data uncertainties can be included by sampling the error distribution of the data and combining the results (e.g., Beloborodov & Levin 2004), this will perform badly in the limit of low signal-to-noise or missing data. These difficulties are related to the fact that these procedures use only a very crude model of the data, that is, that the angles are distributed uniformly and that angle–energy correlations should be absent, which allows no room for discovery of structure in phase space. A fully Bayesian treatment of this problem can treat the phase-space distribution function as an unknown function to be inferred from the data. Modeling the full phase-space distribution function permits simultaneous inference of missing data and properly marginalized probability distributions for dynamical parameters.

Imagine that we have the three-vector positions ${\vec{\bm x}}_i$ and three-vector velocities ${\vec{\bm v}}_i$ at some time t for N planets i and a parameterized model for the gravitational acceleration law (force law per unit mass) ${\vec{a}}_{{{\omega }}}({\vec{\bm x}})$, a function of position ${\vec{\bm x}}$ and a list of parameters $\boldsymbol{{\omega }}$. We wish to obtain an estimate of the posterior probability distribution $p(\boldsymbol{{\omega }}|\lbrace {\vec{\bm x}}_i,{\vec{\bm v}}_i\rbrace)$ for the parameters, where $\lbrace {\vec{\bm x}}_i,{\vec{\bm v}}_i\rbrace$ is the set of all planet positions and velocities. We employ Bayes's theorem as follows:

$$\begin{equation} p(\boldsymbol{{\omega }}|\lbrace {\vec{\bm x}}_i,{\vec{\bm v}}_i\rbrace) = \frac{p(\lbrace {\vec{\bm x}}_i,{\vec{\bm v}}_i\rbrace |\boldsymbol{{\omega }})\,p(\boldsymbol{{\omega }})}{p(\lbrace {\vec{\bm x}}_i,{\vec{\bm v}}_i\rbrace)}, \end{equation} \tag{ 9 }$$

where, as usual, $p(\lbrace {\vec{\bm x}}_i,{\vec{\bm v}}_i\rbrace |\boldsymbol{{\omega }})$ is the likelihood or the probability distribution function for the data given the model parameters, evaluated at the observed values of the data, $p(\boldsymbol{{\omega }})$ is the prior probability distribution function for the parameters, and the denominator is (for our purposes here) a normalization constant.

For our chosen parameterization of the dynamical parameters, $\boldsymbol{{\omega }}$, a broad flat or uniform prior in the space represents a reasonable description of our (assumed) prior knowledge. The much more challenging problem is to specify the likelihood of the dynamical parameters, the conditional probability distribution function $p(\lbrace {\vec{\bm x}}_i,{\vec{\bm v}}_i\rbrace |\boldsymbol{{\omega }})$. Without detailed knowledge of how the solar system formed, this probability distribution is also a representation of our prior beliefs. Ideally our ability to learn from the observed data will not be too sensitive to these beliefs.

It is easier to express beliefs about the angle ϕ_r, radial asymmetry e, and binding energy of each planet than its position and velocity. This is because, given our assumption that the planets constitute an angle-mixed population, Jeans's theorem (Jeans 1915; Binney & Tremaine 2008) tells us that the distribution function, which is proportional to the probability of observing a planet at a certain locus in phase space, is only a function of the integrals of the motion. We wish to assign zero probability to dynamical parameters that lead to any of the computed binding energies $\epsilon ({\vec{\bm x}}_i,{\vec{\bm v}}_i,\boldsymbol{{\omega }})$ being negative (because we defined >0 to be bound). We would also like to express our prior information or assumption that the system is long-lived (in units of the dynamical time), that it is non-resonant, and that we are not seeing the system at any special time, or that the radial angles ϕ_r will be randomly distributed between 0 and 2 π. In the absence of any better information, we will try to be as agnostic as possible about the actions (conserved quantities) of the planets but extremely confident that all radial angles 0 < ϕ_r < 2 π are equally likely.

In the simple spherical or radial situation under consideration here, in which there are no missing data, we can rewrite the likelihood as a function of the planets' radial coordinates, as the orientation of the orbit does not depend on the dynamical parameters:

$$\begin{eqnarray} p({\vec{\bm x}}, {\vec{\bm v}}|\boldsymbol{{\omega }}) \propto p(r,v_r,j^2|\boldsymbol{{\omega }}) &=& |J(\ln {\epsilon },e,\phi _r;r,v_r,j^2)| \nonumber\\ &&\times p(\ln {\epsilon }, e, \phi _r|\boldsymbol{{\omega }}), \end{eqnarray} \tag{ 10 }$$

where, again, we have gone to ln  because dimensioned parameters are usually best handled in the log, and J(ln , e, ϕ_r; r, v_r, j²) is the Jacobian matrix of all the partial derivatives of (ln , e, ϕ_r) with respect to (r, v_r, j²). For spherical potentials, this Jacobian is given by

$$\begin{equation} |J(\ln {\epsilon },e,\phi _r;r,v_r,j^2)| = \frac{2\,\pi }{T_r}\,\frac{1}{\epsilon }\,\left|\frac{\partial e}{\partial j^2}\right|, \end{equation} \tag{ 11 }$$

where the derivative is evaluated at the current location in phase space and T_r is the radial period, which depends on the dynamical parameters and the integrals of the motion. For general α this radial period can be computed numerically (see Section 2). The derivative of the radial asymmetry with respect to the specific angular momentum squared can be written in terms of the perihelion and aphelion distance as

$$\begin{eqnarray} \frac{\partial e}{\partial j^2} &=& A r_0^{\alpha }\,r_{\mathrm{ap}}\, r_{\mathrm{peri}}\frac{r_{\mathrm{ap}}^{3-\alpha } - r_{\mathrm{peri}}^{3-\alpha }}{(r_{\mathrm{ap}}+r_{\mathrm{peri}})^2}\, \left[j^2-A r_0^3\left(\frac{r_{\mathrm{ap}}}{r_0}\right)^{3-\alpha }\right]^{-1}\nonumber\\ &&\times \left[j^2-A r_0^3\left(\frac{r_{\mathrm{peri}}}{r_0}\right)^{3-\alpha }\right]^{-1}. \end{eqnarray} \tag{ 12 }$$

In the special case α = 2, this Jacobian is

$$\begin{equation} |J(\ln {\epsilon },e,\phi _r;r,v_r,j^2)| \propto \epsilon ^{3/2}\,e^{-1}\,M_{\odot }^{-3},\, \end{equation} \tag{ 13 }$$

where we have dropped any terms that do not depend on the mass of the Sun.

As well as assuming that the angles are independently and uniformly distributed, we might further assume that the energy and radial asymmetry of each planet were drawn independently:

$$\begin{equation} p(\lbrace \ln {\epsilon }_i, e_i, \phi _{r,i}\rbrace |\boldsymbol{{\omega }}) = \frac{1}{(2\pi)^N} \prod _{i=1}^N p(e_i,\ln {\epsilon }_i|\boldsymbol{{\omega }}). \end{equation} \tag{ 14 }$$

However, we do not know a priori the distribution function from which radial asymmetries and binding energies were drawn. Fixing $p(e_i,\ln {\epsilon }_i|\boldsymbol{{\omega }})$ to a broad distribution would be making a strong assumption: regardless of what the data say we would continue to dogmatically believe that the distribution function was broad, which could result in poor inferences about the model as a whole.

We can assume that the planets' properties were drawn independently and identically distributed, without making strong assumptions about what that distribution was. This is achieved by introducing auxiliary parameters $\boldsymbol{{\theta }}= \lbrace \boldsymbol{{\theta }}_e, \boldsymbol{{\theta }}_{\epsilon }\rbrace$ that if known would specify the distribution function. Since we do not know these nuisance parameters, introducing a prior $p(\boldsymbol{{\theta }}|\boldsymbol{{\omega }})$ and marginalizing over them,

$$\begin{eqnarray} p(\lbrace \ln {\epsilon }_i, e_i, \phi _{r,i}\rbrace |\boldsymbol{{\omega }}) &=& \frac{1}{(2\pi)^N} \int \mathrm{d}\boldsymbol{{\theta }}\;\, p(\boldsymbol{{\theta }}|\boldsymbol{{\omega }})\nonumber\\ &&\times \prod _{i=1}^N p(e_i,\ln {\epsilon }_i|\boldsymbol{{\theta }},\boldsymbol{{\omega }}), \end{eqnarray} \tag{ 15 }$$

is then part of the inference task. With all of this in place, we apply Equation (9) to get the posterior distribution over the dynamical parameters

$$\begin{eqnarray} p(\boldsymbol{{\omega }}| \lbrace {\vec{\bm x}}_i,{\vec{\bm v}}_i\rbrace) &\propto & \left[\prod _{i=1}^N |J(\ln {\epsilon }_i,e_i,\phi _{r,i};r_i,v_{r,i},j_i^2)|\right] \nonumber\\ &&\times \int \mathrm{d}\boldsymbol{{\theta }}\;p(\boldsymbol{{\theta }},\boldsymbol{{\omega }}) \,\prod _{i=1}^N \,p(e_i,\ln {\epsilon }_i|\boldsymbol{{\theta }},\boldsymbol{{\omega }}),\quad \end{eqnarray} \tag{ 16 }$$

where each planet's value of (ln _i, e_i, ϕ_r,i) is a function of phase-space position $({\vec{\bm x}}_i,{\vec{\bm v}}_i)$ and dynamical parameters $\boldsymbol{{\omega }}$, and the Jacobian is evaluated at each planet's value of (ln _i, e_i, ϕ_r,i).

For situations in which we have missing data or noisy observations, further unknown quantities will be added to the model and marginalized over. As the model becomes more complicated it is more likely that Markov chain Monte Carlo (e.g., Neal 1993) or other approximate computational methods will be required.

6.1. Basic Method

To keep things as simple as possible, at first we model the distribution function $p(\ln {\epsilon }, e|\boldsymbol{{\omega }},\boldsymbol{{\theta }})$ as a product of a top-hat function in ln  from ln _a to ln _b with a top-hat function in e from e_a to e_b. In this context, the phase-space parameter list is

$$\begin{equation} \boldsymbol{{\theta }}= \lbrace \ln {\epsilon }_a, \ln {\epsilon }_b, e_a, e_b \rbrace. \end{equation} \tag{ 17 }$$

In what follows, we only consider values A < ln _a < ln _b < B, where A and B provide very distant (uninformative) limits (below, we will take the limit), and 0 ⩽ e_a ⩽ e_b ⩽ 1. These enforce our assumption or prior information that the system is bound.

For our initial model we also set the prior $p(\boldsymbol{{\theta }},\boldsymbol{{\omega }})$ flat or uniform in all the parameters. Now the marginalization required to compute the posterior (Equation (16)), an integral over all four of the phase-space distribution parameters in $\boldsymbol{{\theta }}$, can be performed analytically; it leaves

$$\begin{eqnarray} p(\boldsymbol{{\omega }}| \lbrace {\vec{\bm {x}}}_i,{\vec{\bm {v}}}_i\rbrace) &\propto & [\ln {\epsilon }_K - \ln {\epsilon }_L]^{2-N} \left[1-[1-e_L]^{2-N}\right.\nonumber\\ &&\left.\quad-[e_M]^{2-N} + [e_M-e_L]^{2-N}\right]\nonumber\\ &&\quad \times \prod _{i=1}^N \big|J\big(\ln {\epsilon },e,\phi _r;r_i,v_{r,i},j_i^2\big)\big|, \end{eqnarray} \tag{ 18 }$$

where ln _L is the lowest planetary binding energy at this point $\boldsymbol{{\omega }}$ in dynamical parameter space, ln _K is the highest, e_L is the lowest planetary radial asymmetry at this point $\boldsymbol{{\omega }}$, e_M is the highest, and we have taken the limit in which the range of the parameters (ln _a, ln _b) goes to infinity, or A → −∞, B → ∞. This posterior probability distribution for $\boldsymbol{{\omega }}$ represents our posterior beliefs marginalized over all possible values of {ln _a, ln _b, e_a, e_b} of the top-hat model.

6.2. Alternative Methods

In our basic method we evaluate the Jacobian |J(ln , e, ϕ_r; r_i, v_r,i, j²_i)| for each planet at the observed radius, radial velocity, and specific angular momentum, for each value of the dynamical parameters $\boldsymbol{{\omega }}$, and multiply the magnitudes of the determinants of these Jacobians together. Following Equation (18), we multiply this product with the value of the marginalized product of the distribution function and finally multiply this with the (uniform) prior on the dynamical parameters in order to obtain the posterior probability distribution for the dynamical parameters.

It is instructive to look at the magnitude of the determinant of the Jacobian for each planet as a function of the dynamical parameters in Figure 5. These factors seem to show a clear preference for the dynamical parameters of each planet to lie on its virial locus. Why is this? For nearly circular orbits the last two factors in Equation (12) are equal and are unbounded for perfectly circular orbits; this results in a natural preference for nearly circular orbits for all of the planets, i.e., for each planet to lie somewhere along its virial locus. At the same time the measured non-zero radial velocity limits the extent to which a planet's orbit can be circularized by a suitable choice of the dynamical parameters and, thus, how close to the circular singularity an individual planet can be brought. This naturally leads to a weighting scheme in which planets that are closer to circular orbits (as measured by the observed angle between the planet's velocity and the tangent to the circular orbit going through its present location) receive a higher weight in the analysis than planets on less circular orbits. This weighting scheme is good in that planets on more circular orbits are indeed more informative about the potential than planets on less circular orbits. Therefore, it may seem that the decision to model the radial asymmetries as uniform in e was fortuitous, but in fact, as we will discuss below in Section 7.1, the results obtained using the basic method are robust in that more flexible models learn that the top hat in e leads to a good interpretation of the data.

Zoom In Zoom Out Reset image size

The posterior probability distribution is shown in Figure 6. A strong peak is apparent near α = 2 and the Newtonian solar value for ln A. The width of the probability distribution is indicated by the 95% and 99% posterior contours. These contours are defined to enclose the smallest area that holds 95% and 99%, respectively, of the posterior probability distribution.

Zoom In Zoom Out Reset image size

In order to infer the exponent α of the force power law we perform a second marginalization of the posterior probability distribution, this time over the (for our purposes) uninteresting parameter ln A. This gives a posterior probability distribution for the parameter α,

$$\begin{equation} p(\alpha |\lbrace {\vec{\bm x}}_i,{\vec{\bm v}}_i\rbrace) = \int \mathrm{d}\ln A\, p(\boldsymbol{{\omega }}|\lbrace {\vec{\bm x}}_i,{\vec{\bm v}}_i\rbrace). \end{equation} \tag{ 19 }$$

This probability distribution is shown in Figure 7. The 95% and 99% posterior intervals in this figure are defined to exclude 2.5% and 0.5%, respectively, of the distribution on either side of the central region.

Zoom In Zoom Out Reset image size

The result of the inference is not a value for the parameters but a posterior probability distribution. A "best-fit" value for α could be obtained according to various criteria. For example, the posterior mean minimizes the expected square difference from the true parameter. However, since the posterior distribution is multi-modal, a point estimate does not capture our uncertainty. The posterior is better summarized by a credible interval, containing a given fraction of the probability mass. The 95% posterior interval is 1.989 < α < 2.052. This compares favorably with the results obtained by Newton (1687) who inferred α = 2 from a much richer data set (though a less rich model set). Modern tests constrain the value of the exponent α to a fractional accuracy of 10⁻⁸–10⁻⁹ on solar system scales (Adelberger et al. 2003; Fischbach & Talmadge 1999) using lunar laser ranging tests (Williams et al. 2004) and Keplerian tests comparing G(r)M_☉ and the rate of precession for different planets (Talmadge et al. 1988).

Our method appears to work well: the true dynamical parameters A and α are plausible under a fairly tight posterior distribution, see Figure 6. The confidence intervals from the frequentist tests, Figure 4, are somewhat broader than the posterior intervals, but are not directly comparable. Frequentist confidence intervals do not represent beliefs about the parameters for this data set, simply a region that satisfies some sampling properties. As a result it is unsurprising that the frequentist confidence intervals appear less good if interpreted as beliefs given the data. It is difficult to use the frequentist confidence intervals to make statements about α alone: without the assumptions of the Bayesian method we cannot marginalize over the parameter A. One approach that allows frequentist guarantees is to set A adversarially, but this pessimistic view would give unreasonably large error bars for α.

That the posterior distribution $p(\alpha |\lbrace {\vec{\bm x}}_i,{\vec{\bm v}}_i\rbrace)$ turns out to be multi-modal is no surprise in the light of the virial considerations from Figures 1 and 2. The two peaks in $p(\alpha |\lbrace {\vec{\bm x}}_i,{\vec{\bm v}} _i\rbrace)$ correspond to the two main regions in parameter space in which the different virial loci cross. These virial considerations, and also the frequentist techniques, are very similar to the probabilistic approach in that they all prefer each planet to be in a non-special region of radial angle space, that is, between perihelion and aphelion; this can only happen simultaneously near the points in parameter space at which virial loci cross. The advantage of the probabilistic approach is that it explains and quantifies this reasoning, and uses it to set formal limits on the dynamical parameters.

7.1. Sensitivity Analysis

The strongly peaked Jacobian factors shown in Figure 5 depend strongly on the choice of parameterization. If the model was expressed using the radial asymmetries squared rather than the radial asymmetries, the Jacobians, some of which are shown in Figure 8, are much flatter. If the model is maintained by also transforming the prior over radial asymmetry distributions to the new parameterization, the posterior distribution over dynamical parameters will be unchanged. Setting all the radial asymmetries close to 1 is penalized elsewhere in the mathematical expression.

Zoom In Zoom Out Reset image size

However, setting the prior over radial asymmetry distributions to consider only uniform distributions was a strong assumption. If we had only considered distributions uniform in the radial asymmetry squared, this combined with the flat Jacobians in Figure 8 would not have given as tight a posterior distribution. Fortunately we can infer a suitable parameterization, or equivalently the distribution function from a richer family of possibilities, from the data.

Figure 9 shows the posterior over α from using the two more flexible priors outlined in Section 6.2. Assuming that the radial asymmetry is uniform in one of $\sqrt{e}$, e, or e² gave a 95% credible interval of 1.990 < α < 2.035. A variety of different flexible priors for the radial asymmetry distribution gave very similar results. Allowing the ln A constant in the force law to come from a flexible family of distributions changed the posterior very little. The 95% credible interval from the flexible non-parametric prior on both distributions gave a 95% credible interval of 1.991 < α < 2.040. The fine details of the posterior distribution, specifically the relative mass in the two modes, are sensitive to prior assumptions. This is unsurprising with only eight data points. However, the position of the bulk of the posterior mass, which gives the confidence interval, is surprisingly robust across a wide range of modeling choices.

Zoom In Zoom Out Reset image size

As a warning, it is possible to obtain poor results through bad modeling choices. The reasonable posterior from our basic method, Figure 7, was to some extent due to chance. The family of possible radial asymmetry distributions is inflexible as it forces this distribution to be flat over some range in radial asymmetry. Using the basic method it is not possible to learn the true distribution, which is not flat in radial asymmetry, from data, because the truth is assigned zero prior probability. We were fortunate to choose a parameterization in which the inflexible family was not too unreasonable. The non-parametric approach is much more flexible in that its weak smoothness constraints assign non-zero prior probabilities to all reasonable distribution functions such that the true one could be picked out with enough data. However, care must also be taken with non-parametric modeling. An unsmooth Dirichlet prior over bin heights in a non-parametric prior over distributions gave a broad and biased distribution over α.

7.2. Future Work

First of all we thank Scott Tremaine for many inspiring discussions that helped shape this paper. It is also a pleasure to thank Kyle Cranmer, Dustin Lang, Yuri Levin, Phil Marshall, Bill Press, and Hans-Walter Rix for very helpful discussions and the anonymous referee for valuable comments. This research was partially supported by NASA (ADP grant NNX08AJ48G), the NSF (grant AST-0908357), and a Research Fellowship of the Alexander von Humboldt Foundation. This project made use of the HORIZONS System provided by the Solar System Dynamics Group of the Jet Propulsion Laboratory, the NASA Astrophysics Data System, and the open-source Python modules scipy, numpy, and matplotlib.