EDP Sciences
Free access
Volume 565, May 2014
Article Number A56
Number of page(s) 8
Section Celestial mechanics and astrometry
DOI http://dx.doi.org/10.1051/0004-6361/201322766
Published online 08 May 2014

© ESO, 2014

1. Introduction

The astrometric determination of asteroid masses is the technique of calculating the mass of asteroids (hereafter called “perturbing asteroids” or “perturbers”) from the gravitational perturbations they exert on the motion of other asteroids (called “test asteroids”).

The usual procedure to calculate the mass of one asteroid is, first, to identify potential test asteroids by searching for close approaches with the perturbing asteroids and then computing some parameter (usually the deflection angle or the change in mean motion) characterising the effect of the perturbing asteroid, see e.g. Galád & Gray (2002). Promising cases are then further investigated to see if they yield meaningful results, since the outcome also depends on the number, distribution and quality of the available astrometric observations of the test asteroid. Asteroids not making close approaches but moving in near-commensurable orbits can also be useful for mass determinations (Kochetova 2004).

A better but more laborious way to find suitable close encounters is to integrate the orbits of test asteroids over the period of available observations with and without the perturbations of the perturbing asteroid (Michalak 2000). Significant differences in the geocentric coordinates between the two calculations indicate cases that should be investigated further.

Initially, asteroid masses were computed from their perturbations on one other asteroid, the first instance being the mass of 4 Vesta computed from the orbit of 197 Arete by Hertz (1968). When several mass values for the same asteroid are calculated from a number of test asteroids, the usual practice is to compute a weighted mean along with the associated uncertainty as the final value, as in Michalak (2000). Mass determinations where one or a few masses are determined simultaneously from a number of test asteroids have been attempted on a modest scale in the past – see Sitarski and Todororovic-Juchniewicz (1995), Goffin (2001) and others.

This paper gives the results of an attempt to simultaneously determine a large number of asteroid masses (230) from a substantial number of test asteroids (349 737). Section 2 explains the basic principles and problems of simultaneous asteroid mass determination and in Sect. 3 the details of the actual calculations are given. Section 4 briefly deals with the observations used and explains the statistical analysis of the residuals that was used to assign weights to, or reject observations. In Sect. 5 the results are presented and Sect. 6 discusses them.

Throughout this paper, asteroid masses are expressed in units of 10-10 solar masses (M) and generally given to 3 decimal places.

2. Simultaneous mass determination

The principle of the simultaneous determination of asteroid masses is simple: use the astrometric observations of N test asteroids to determine amongst others the masses of M perturbing asteroids. This direct approach has a number of distinct advantages. First, there is no need to search for close encounters to the M perturbers, assess their efficiency and select promising pairs for the actual calculations (the same can be said for near-commensurable orbits). All test asteroids are treated equal and close approaches will appear anyway when numerically integrating the N orbits. Next, one doesn’t have to compute weighted means of masses and their standard deviations. Also, the mutual correlations between the M masses are taken into account in the global solution. But all this comes at a cost: computer resources!

2.1. The normal equations

In the process of mass determination, the equations of condition of the linearised problem are usually transformed into the normal equations: (1)The coefficient matrix A is square, positive definite and symmetric about the main diagonal. For an “ordinary” mass determination (N = 1,M = 1) it is of size 7 × 7 (rows × columns). Note that the coefficients below the main diagonal do not need to be stored, but that part of the matrix is often used as temporary storage when solving the equations or computing the inverse A-1.

The normal equations for simultaneous mass determination (N > 1, M ≥ 1) are easily assembled from the normal equations of the individual test asteroids. Sitarski & Todororovic-Juchniewicz (1992) give a schematic representation for a simple case (N = 2, M = 1). In general, A is of size (6N + M) × (6N + M). While this is not prohibitive for modest applications, it is clear that A will grow rapidly in size with large values of M and especially N, and eventually exceed the available computer memory.

The full normal Eqs. (1) can be written as: (2)where:

  • A1,1,...,AN,N are the 6 × 6 blocks of the partial derivatives for the 6 orbital elements of the N test asteroids.

  • A1,M,...,AN,M (each of size 6 × M) and AM,M (size M × M) correspond to the M masses.

  • X1,...,XN are the variables pertaining to the N asteroids, namely the corrections to the 6 orbital elements (or components of the state vector) at epoch.

  • XM are the variables pertaining to all perturbing asteroids, namely the corrections to the M masses.

The zeros represent 6 × 6 blocks of zeros corresponding to terms of second and higher order which are of course neglected in the linearised formulation of the problem. The matrix clearly contains a staggering amount of elements filled with zeros!

2.2. Solving the normal equations

One way to reduce computer storage is to store only the non-zero blocks of A: (3)When stored as a single array, the size of A is thus reduced to (6N + M) × (6 + M), or, if AM,M is stored separately, to 6N(M + 6) + M × M. The Cholesky decomposition algorithm used to invert matrix A was adapted to cope with this new form. In that way, long test calculations were done up to N = 300 000,M = 105 (this proved to be a practical upper limit). Then, following a suggestion by E. Myles Standish, a formulation was adopted that only requires a matrix of dimensions (6 + M) × (6 + M)!

The normal Eqs. (2) can be written as the pair: (4)The first step is to eliminate the Xi’s from the second equation. From the first, we get: (5)Substituting this in the second and rearranging gives: (6)With the XM’s one can now solve Eq. (5) for the Xi’s.

2.3. Standard deviations and correlations

The formal standard deviations and correlations come from the inverse of the A-matrix. Representing the components of A-1 by α’s, we have from the definition: (7)The multiplication gives: (8)The second equation of (8) gives: (9)Substituting this in the fourth yields: (10)which has already been calculated and inverted in (6).

With (9) one gets αi,i from the first equation of (8): (11)The sequence of the computation is thus: (10), (9), (11).

2.4. Test asteroids

The test asteroids are those among the first 350 000 numbered asteroids with perihelion distance smaller than 6.0 au (to exclude distant objects such as KBO’s...), a total of 349 737. Their orbits were improved using perturbations by the major planets, Ceres, Pallas and Vesta, before using their orbital elements and observations as input for the simultaneous mass determination.

The positions of the eight major planets and the Moon are taken from the JPL ephemeris DE422, which is an extended version of DE421 (Folkner et al. 2008). The coordinates of the three largest asteroids initially come from a separate orbit improvement of these bodies.

2.5. Perturbing asteroids

Selecting the perturbing asteroids is basically quite simple: just take those with the largest mass. Reliable masses have been calculated for a number of asteroids, but for most others one must rely on estimates based on the diameter and the taxonomic class. Mass estimates for these asteroids were computed assuming a bulk density of 1.80 g/cm3 for class C (including B, D, F, G, P, T and X), 2.20 g/cm3 for class S (including A, E, K, Q, R and V) and 4.20 g/cm3 for class M.

Diameters were selected from the sources given in Baer et al. (2011), or from occultation results (Broughton 2011), or else from the results of the WISE (Masieri et al. 2011) or IRAS (Tedesco et al. 1992) surveys. For 4 Vesta the value from the Dawn spacecraft was used (Russell et al. 2012). The diameter of 52 Europa is from Merline et al. (2013) and that of 349 Dembowska from Majaess et al. (2008). In case the dimensions of a tri-axial ellipsoid were given, the equivalent diameter (of a sphere of the same volume) was calculated. When none of these sources gave a diameter, a nominal value D was computed from the absolute magnitude H: In this way, mass estimates were calculated for the first 100 000 numbered minor planets and those with a mass estimate >0.005 × 10-10 M and perihelion distance q < 6.0 au, a total of 257, were withheld as perturbers for a first test run with the 349 737 test asteroids. After a few iterations it turned out that 279 Thule and the Trojans in the list didn’t yield acceptable masses, so a new list was created with q < 4.0 and these 230 asteroids were used in the final calculations.

Though it would be possible to include the masses of the major planets in the solution, they were held fixed at the DE422 values.

3. Computations

The basics of the computations have been given in many previous papers, amongst others in Goffin (1991) and Goffin (2001). The force model includes the Newtonian attraction due to the Sun, the 8 major planets and the Moon, and the M perturbing asteroids (one less when the test asteroid is one of the perturbers). The heliocentric coordinates of the major planets and the Moon are stored in a planetary file together with those of the M perturbing asteroids. Also taken into account are the first-order relativistic effects due to the Sun as given in Ivantsov (2008) and the effects of the solar oblateness.

The longest part of the calculations is the numerical integration of the N test asteroid orbits and of the differential equations of all the required partial derivatives. These calculations only became feasible with the advent of personal computers with multi-core processors.

By the very nature of the problem, the solution is obtained by successive approximations. Each iteration cycle consists of four consecutive steps, each actually corresponding to a computer programme:

  • P1: numerical integration of the N test asteroids.

  • P2: solution of the normal equations.

  • P3: recomputation of the coordinates of the perturbing asteroids.

  • P4: statistical analysis of the residuals.

The programme P1 takes care of the numerical integration of the N orbits. The mutual distances between perturbing and test asteroids are monitored and the step size is decreased when necessary. Orbital elements of all asteroids and masses of the perturbing asteroids are read in from the elements file on disk. The observations are weighted using the results of the statistical analysis, after correcting for any magnitude effects and systematic bias (see Sect. 4). The normal equations for each test asteroid are output to a disk file and all residuals (O-C) are stored in the residuals file together with relevant information (observatory code, time, visual magnitude of the asteroid, etc.).

Programme P2 reads in the separate normal equations, assembles the normal equations of the global solution, inverts the coefficient matrices and computes the 6N + M unknowns (initially 2 098 652 in total), their standard deviations and the correlation coefficients. This is actually done twice during each iteration. After the first solution, the validity of the resulting M masses is checked. Following Baer & Chesley (2008), a mass is considered acceptable when its significance, this being the ratio of the mass to its standard deviation, is larger than 2 (si = mi/σmi > + 2.0), and if the resulting density lies between 0.5 and 8.0 g/cm3. The rows and columns of the normal matrix corresponding to unacceptable masses are then modified (by filling in 0’s or 1), and the right-hand side is changed so that the second solution will force them back to their initial (estimated, and thus acceptable) values. Practise has shown that these two solutions suffice for that purpose. The improved masses and orbital elements from this solution are stored in the elements file.

The third programme (P3) uses the improved elements of the M perturbing asteroids to compute new heliocentric rectangular coordinates and copies them into the planetary file.

Programme P4 uses the residuals file to perform a statistical analysis of the performance of all observatories; the detailed description is the subject of Sect. 4. The results are stored on disk and used by P1 to assign weights to the observations. This analysis is done every 3 iterations. The sequence of the computations is thus: 3 × (P1 + P2 + P3) + P4.

4. Observations, weights and outlier rejection

4.1. Astrometric observations

The principal source of minor-planet observations are the files maintained at the Minor Planet Center (Cambridge, Mass.); the version of January 2013 was used. Throughout the years I have re-reduced most of the older observations of the first 1000 numbered asteroids made since their discovery and published in the open literature. They consist of series of visual micrometer and meridian observations, and a few series of photographic observations. Wherever possible, they were re-reduced with the aid of modern star catalogues. No further details will be given here. These observations were added to – and replaced some of – the data from the MPC. For the optical observations, those given to a precision worse than 0.01 s in right ascension or 0.1′′ in declination were not used. The total number of observations of the N test asteroids is over 89.4 million.

In the past, most researchers have used various methods for rejecting outliers and assigning weights to the observations. Some have used an iterative rejection criterion, or have assigned a-priori weights depending on the type of observation and/or the time period. Such simplifications may be easy to implement but do not reflect the real accuracy of the observations. Relevant weights must therefore be derived from a statistical analysis of the residuals resulting from orbit improvements. Since these two steps are interdependent, this is an iterative procedure.

Serious thoughts to the subject were first given by M. Carpino1 and a thorough analysis was presented in Carpino et al. (2003). However, at that time I decided to continue my own statistical analysis for two reasons. The Carpino et al. paper was not yet published when I embarked on the project of mass determinations, so no practical weighting scheme was available. Also, the many old observations which I had re-reduced are mostly not present in the MPC data set and thus I had to provide their statistics myself.

4.2. Grouping of observations

The complete set of observations must be divided into groups with (supposedly) homogeneous properties, i.e. by instrument (observatory, telescope), observation and reduction method, and further by any type of information that can be relevant (observer, type of measuring instrument, number of readings taken etc.). In practice, all this information is not readily available and has not been stored with the current observation format of the Minor Planet Center. Therefore a division was made per observatory code (which can thus refer to different instruments and observers) and further by the following 4 observational methods: transit and meridian (T), visual micrometer (M), photographic (P) and CCD (C).

Each of these observation groups was further split up into a number of time bins as follows. A first time bin was “filled” until its time length exceeded 2500 days or the total number of observations exceeded 100 000; then the next bin was begun. The first condition limits the time length of the bins for observatories with few observations, while the second allows for a more detailed analysis of the more productive surveys. Time bins with less than 50 observations were considered too small for analysis and were combined into the 4 groups of observatory code 500 Geocentric. This code groups all observations that have been corrected for parallax and the reference to the original observatory is thus lost anyway.

4.3. Statistical analysis

One analysis cycle consists of the following steps.

  • Using the residuals file created by P1, the observations aregrouped in time bins as decribed above.

  • Each time bin of each group is inspected for outliers. The basic idea is that each bin consists of a population of “good” data having Gaussian distribution and a population of outliers. To remove these outliers, the kurtosis parameter k is used (see Carpino et al. 2003, for a definition). Its value is equal to 3 for a normal distribution, and substantially larger when outliers are present. The calculation of k and the outlier removal is done iteratively. The residual with the largest absolute deviation from the mean μ is removed and μ and k are recomputed. This is repeated until k ≤ 3 or, as a safety measure, until k starts to increase rather than decrease.

    thumbnail Fig. 1

    Magnitude equation in right ascension for the visual micrometer observations made at 030 Arcetri (5147 observations, 1894–1919). Top graph: residuals versus time, uncorrected (μ = +1.94″, σ = 2.07″). Middle graph: residuals versus visual magnitude of the asteroid; the least-squares regression line is R = −6.92 + (0.77 ± 0.01)MV. Bottom graph: residuals versus time after magnitude correction (μ = 0.00″ = 1.86″).

    Open with DEXTER

    Table 1

    Minor planets for which no acceptable mass could be calculated.

  • The complete set of residuals of each group (observatory code + observation technique) is then analysed for a (possible) magnitude effect. This is an observational bias whereby a measurement error is correlated with the object’s brightness. Such an effect is to be expected for observations where the right ascension is measured by timing transits, such as meridian and most micrometer observations (Fricke 1979). It can also occur in photographic and CCD astrometry due to imperfections in the optical system, e.g. coma. To model the magnitude effect, a least-squares regression line of the form is computed for both right ascension and declination residuals (R is thus cosδΔα or Δδ). The quantity MV is the computed visual magnitude of the minor planet. In addition are computed: the total magnitude range ΔM of the set, the standard deviation σb of b, the correlation coefficient r (Pearson) and the statistic where N is the number of residuals. A magnitude effect is accepted as significant if: Figure 1 gives an example of a substantial magnitude effect in right ascension.

  • With outliers removed and after correcting for magnitude equation (if statistically significant) the next step is then to compute for each time bin of that group the values of the mean deviation μ (hereafter called the bias), the standard deviation σ and the standard deviation of the mean where N is the number of remaining residuals in the bin. A bias is accepted as statistically significant when it exceeds 2 times its standard deviation (| μ | /σμ > 2). For the observations of code 248 Hipparcos the bias was set equal to zero, since they are by definition on the ICRS system.

The results are applied in P1 as follows. Observations are considered as outliers, i.e. are given zero weight, if their residuals deviate more than 3σ from the mean. The weight of an observation is calculated as w = 0.1 /σ2 (unit weight thus corresponds to ). If there are n observations of the same asteroid by the same observatory (except for code 248 Hipparcos) on one night, then their weight is further divided by , so that the total weight of the n observations is not wn but . This procedure is a simple way to avoid overweighting them with respect to older measurements. The reason is that their residuals are strongly correlated because the observations were very probably made using the same reference stars.

5. Results

The computations with the compact storage formulation (Eq. (3), with N = 300 000 and M = 105) had indicated that the results were sensitive to the weighting of the observations, and in particular to the recomputation of the statistics. Therefore, and in order to avoid a possible initial divergence of the iterations, the first three iterations were done without changing the weights or recomputing the perturber positions, and the recalculation of the observation statistics (P4) was only set in after the 6th iteration (and then every 3 iterations).

After the main calculation was finished, a number of test computations were performed to assess the sensitivity of the solution to the starting values of the masses.

5.1. Main calculation

At the start of the computations the masses of the 230 perturbing asteroids were set to their estimated values. After five iterations it turned out that 58 of them persistently yielded negative mass values in the first solution of the (unmodified) normal equations; they are listed in Table 1 (upper part). It was decided to remove them and to restart the computations with the 172 remaining perturbing asteroids, partly to reduce computing time2.

Due to the improvement of the N orbits and the reweighting of the observations at each iteration, the mean error of one observation decreased to 0.23″ and the total number of equations of condition used in the solution increased from 164.28 to 168.51 million. Convergence for most masses was reached after 14 iterations, i.e. after that their values started to oscillate by a few units of the 4th decimal. The mass of Pallas however showed a general decline and later started fluctuating by one or two units of the 3rd decimal. The calculations were finally stopped after the 29th iteration; the total number of linear equations used in the last iteration was 168 506 583.

Table 2

Twenty closest approaches between test and perturbing asteroids.

For 40 of the 172 perturbing asteroids no acceptable mass could be derived according to the significance and density criteria given above. They are listed in Table 1 (lower part) together with their taxonomic class, estimated diameter and mass used in the dynamical model. Table 2 gives the 20 closest approaches between test and perturbing asteroids.

As is to be expected, there are high correlations between the perturber masses and the semi-major axes of some of the test asteroids. The 10 largest values are given in Table 3; all of them involve high-numbered test asteroids. The mutual correlations between the M masses turn out to be generally low; the 10 largest values are listed in Table 4. All correlations between the 10 largest asteroid masses are 0.05.

Table 3

Ten largest correlations between the masses of perturbers and the semi-major axes of the test asteroids.

Table 4

Ten largest mutual correlations between the perturber masses.

5.2. Influence of initial assumed masses

After concluding the main calculation, the perturber masses were set to zero and the iterations were resumed without reweighting the observations, i.e. with the iteration sequence P1 + P2 + P3. After three iterations, the masses were back to the results of the main calculation to within four decimal places. Next, the masses were set to twice their assumed values and the calculations were continued. Again, three iterations sufficed. This indicates that the least-squares algorithm is rather insensitive to the starting mass values.

Then, the main calculations as descibed above were repeated twice, first with all perturber masses initially set to zero and then set to twice their estimated values, but this time with reweighting applied. The statistical analysis obtained at the end of the previous set of iterations was used at the beginning of the new test series. Convergence was deemed to be reached in both cases after 20 iterations; the number of linear equations in the last iteration steps was less than in the main calculation: −20 222 (−0.012%) and −30 237 (−0.018%) respectively.

Table 5, available at the CDS, gives the 132 asteroids with acceptable astrometric masses, 49 of which appear to be new. The table gives: minor planet number and name (asterisk added for new mass determinations), mass, standard deviation, significance, mass difference when starting from zero masses and from masses twice the assumed values, adopted equivalent diameter, taxonomic class, resulting bulk density. In this table the mass, standard deviation and significance are from the main iteration. The mass differences between the last two test iterations and the main one are also included (column Δm0 for the zero and column Δm2 for the double starting mass values).

6. Discussion

The masses given in Table 5 are presented “as such”, i.e. as the outcome of the least-squares adjustment described above. The bulk density is a quantity derived from the mass and the assumed diameter and its uncertainty is thus affected by the uncertainties of both the mass and the diameter. No physical interpretation of the densities will be attempted here. Some results conform to what is to be expected for their taxonomic class, while others are higher or lower. One should bear in mind that the taxonomic class is derived from the spectral properties of the surface of an asteroid and as such does not say much about its internal composition or porosity.

The mass and size of 4 Vesta is accurately known from the orbit of the Dawn spacecraft. Russell et al. (2012) give a value of 2.59076 ± 0.00001 × 1020 kg corresponding to 1.30260 ± 0.00001 × 10-10 M, an equivalent diameter of 523.15 km and a density of 3.456 g/cm3. The result from this study, to 4 decimal places, is 1.2995 ± 0.0015 yielding a density of 3.448 g/cm3, in good agreement with the Dawn result.

It would be a constructive achievement to develop an objective criterion to exclude test asteroids that do not significantly contribute to the solution. The correlations between the masses and the semi-major axes come to mind as an obvious choice. Using the sum of the absolute values of the M correlation coefficients as an indicator for the usefulness of a test asteroid, test calculations were made with N = 10,000 and M = 10 (the 10 largest masses). With a limit for the sum of 0.01, about 15% were excluded and the average difference in the 10 masses was 0.0003; this became 48% and 0.0012 with 0.02 as the limit. It seems worthwhile to test this criterion on the full set of test asteroids.

In general, the two test calculations clearly show that the most influential factor is the treatment of the observations (weighting and outlier rejection). In that respect, it should be remarked that the observations reduced with the USNO A1.0, A 2.0 and B 1.0 star catalogues have not been corrected for the systematic biases as described in Chesley et al. (2010) or in Röser et al. (2010). It should be worthwile to include these in future work.

The planetary ephemeris uncertainties are negligible in this case, certainly the positional errors. The planetary masses come from spacecraft encounters, and DE421/422 is up-to-date. The best asteroid-produced planetary mass of Jupiter, for example, was inaccurate by about 6 parts in 106; the spacecraft-produced values are two orders of magnitude more accurate.

The Gaia satellite, launched at the end of 2013, will in the course of its 5 year nominal duration systematically observe ~300 000 asteroids brighter than V = 20 with an accuracy of 0.3–5 mas. This will enable the determination of the masses of about 150 individual asteroids from close encounters, ~100 with an accuracy of 50% and ~50 with an accuracy of 10% (Mouret et al. 2008). The Gaia star catalogue will enable the reduction of future, and the re-reduction of old observations with a very high level of accuracy (about 10 mas, see Desmars et al. 2013), thus considerably increasing the accuracy of simultaneous asteroid mass determinations.

Finally, Table 6, available at the CDS, gives an overview of mass results from recent determinations for those minor planets given in Table 5 and published in roughly the last ten years. It gives: minor planet number and name, mass, standard deviation, reference. For a more complete list including older determinations, the reader is referred to Zielenbach (2011). Besides dynamical determinations using asteroids as test objects, there are also determinations from natural satellites and from planetary ephemerides. Note that these use completely different types of observations.


Then on his web site; not available any more now.


For the calculations with 230 asteroid masses, the total reported CPU time of the programme P1 was about 96.5 h per iteration, slightly over 59 h with 172 perturbers.


I am most indebted to E. Myles Standish for his keen interest in the subject, his willingness to critically review this paper and for bringing in his expertise where my mathematical knowledge failed.


All Tables

Table 1

Minor planets for which no acceptable mass could be calculated.

Table 2

Twenty closest approaches between test and perturbing asteroids.

Table 3

Ten largest correlations between the masses of perturbers and the semi-major axes of the test asteroids.

Table 4

Ten largest mutual correlations between the perturber masses.

All Figures

thumbnail Fig. 1

Magnitude equation in right ascension for the visual micrometer observations made at 030 Arcetri (5147 observations, 1894–1919). Top graph: residuals versus time, uncorrected (μ = +1.94″, σ = 2.07″). Middle graph: residuals versus visual magnitude of the asteroid; the least-squares regression line is R = −6.92 + (0.77 ± 0.01)MV. Bottom graph: residuals versus time after magnitude correction (μ = 0.00″ = 1.86″).

Open with DEXTER
In the text