Andrea Milani
email: milani@dm.unipi.it

Giovanni F Gronchi

Davide Farnocchia

Zoran Knežević

Robert Jedicke

Larry Denneau

Francesco Pierfederici

Topocentric Orbit Determination: Algorithms for the Next Generation Surveys

Keywords: Celestial Mechanics, Asteroids, Dynamics, Orbits

The process of calculating a good orbit from astrometric observations of the same object involves three main steps: preliminary orbit determination, least squares orbit fitting, and quality control assessing the orbit's uncertainty and reliability. For the next generation sky surveys, with much larger number density of observations, new algorithms, or at least substantial revisions of the classical ones, are needed. The classical theory of preliminary orbit algorithms was incomplete in that the consequences of the topocentric correction had not been fully studied. We show that it is possible to rigorously account for topocentric observations and that this correction may increase the number of alternate preliminary orbits without impairing the overall performance. We have developed modified least squares algorithms including the capability of fitting the orbit to a reduced number of parameters. The restricted fitting techniques can be used to improve the reliability of the orbit computing procedure when the observed arcs have small curvature.

False identification (where observations of different objects are incorrectly linked together) can be discarded with a quality control on the residuals and a 'normalization' procedure removing duplications and contradictions. We have tested our algorithms on two simulations based on the expected performance of Pan-STARRS -one of the next generation all-sky surveys. The results confirm that large sets of discoveries can be handled very efficiently resulting in good quality orbits. In these test we lost only 0.6 to 1.3% of the possible objects with a false identification rate in the range 0.02 to 0.06%.

A C C E P T E D M A N U S C R I P T 1 Introduction

The problem of preliminary orbit determination is old, with very effective solutions developed by [Laplace 1780] and [START_REF] Gauss | Theory of the Motion of the Heavenly Bodies Moving about the Sun in Conic Sections[END_REF]]. However, the methods of observing Solar System bodies have changed radically since classical times and have been changing even faster recently due to advances in digital astrometry. The question is, what needs to be improved in the classical orbit determination algorithms to handle the expected rate of data from the next generation all-sky surveys? Alternatively, what can we now use in place of the classical algorithms?

The issue is not one of computational resources because these grow at the same rate as the capability of generating astrometric data. Reliability is the main problem when handling tens of millions of detections of Solar System objects. An algorithm failing once in 1, 000 usages may have been considered reliable a few years ago, but now we must demand better performance. This is particularly important because of the strong correlation between difficulties in the orbit computation and the scientific value of the discovered object. Main Belt Asteroids (MBAs) are commonplace and their orbits are easily computed. In modern asteroid surveys only a few in a 1, 000 of the objects (to a given limiting magnitude) are the more interesting Near Earth Objects (NEOs) while a few in 100 are the equally interesting Trans Neptunian Objects (TNOs); in both cases the computation may be much more difficult for reasons explained later. Thus an algorithm that successfully computes orbits for 99% of the discoveries may still fail on a large fraction of the To determine if there are problems with the classical orbit determination techniques we must re-examine whether the assumptions invoked in the derivations of the classical techniques are still appropriate. E.g., neglecting multiple preliminary orbits and dismissing the topocentric corrections is wrong when searching for NEOs; neglecting the case in which the observed track on the sky has insignificant curvature is wrong for TNOs. Another consideration is the immensely superior computing power available to us: in the trade off between simpler and more reliable computations we almost always select the latter while the classical authors were forced to do the opposite.

Problems and solutions: preliminary orbits

It has long been known that Gauss' [START_REF] Gauss | Theory of the Motion of the Heavenly Bodies Moving about the Sun in Conic Sections[END_REF]] and Laplace's [Laplace 1780] preliminary orbit determination methods are equivalent to some level of approximation but it is important here to understand the approximations under which this is true and to check whether they are still applicable for contemporary observation. If the observations are geocentric the two methods are equivalent up to the algebraic equation of degree 8 corresponding to a quadratic approximation in time. We will show below that with topocentric observations they are not equivalent. [Crawford et al. 1930, page 99] and [Marsden 1985] have shown that Gauss' method, at least in the version of [START_REF] Merton | A modification of Gauss's method for the determination of orbits[END_REF], can be used with topocentric observations employing the same formulae but Laplace's method can not.

The difference between topocentric and geocentric observations is not

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

negligible apart from some special cases. Laplace's method can account for topocentric observations but only within an iterative procedure whose convergence can not be guaranteed. Following a suggestion by [Poincaré 1906] Laplace's method can be modified to account for topocentric observations in the degree 8 equation.

The very useful qualitative theory of [Charlier 1910] that allows us to compute the number of preliminary orbit solutions for Laplace's method does not apply to the modified method. It also does not apply to Gauss' and can not describe iterative methods. This is a problem if we are concerned with the reliability of our orbit determination algorithm. Although it is always possible to improve the orbit with iterative methods, if the first approximation provides a wrong number of solutions it is possible to completely miss the correct one!

In Section 4 we develop a new qualitative theory, fully accounting for topocentric observations, for the solutions of the equation of degree 8 for both Gauss' and the modified Laplace-Poincaré methods. We show that the number of solutions can be larger than in Charlier's theory, e.g., there may be double solutions near opposition and triple solutions at low solar elongations. Some examples are given to show that the existence of additional solutions can affect the reliability of the orbit determination for real asteroids. This progress tips the balance in favor of applying Gauss' method to preliminary orbit determination because it can handle topocentric observations in a natural way.

While these mathematical arguments are interesting an abstract qualita-

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

tive theory is useful only if it is exploited by software which can compute and keep track of all the solutions of the algebraic equations. This has not been the case to date. A tricky problem arises when, after the solution of the degree 8 equation, an iterative method is used to improve the preliminary orbit: the number of limit points of the iteration may be different from the number of starting solutions. As a compromise we use a two step procedure in which the two versions, with and without iterative improvement, are used in sequence.

Problems and solutions: weakly determined orbits

Preliminary orbit methods fail when the components of curvature in the observations are so poorly determined that even their sign is uncertain. This happens when either the observed arc is too short or the object is very distant.

In Section 5 we present how to detect when these conditions occur and how to estimate the corresponding orbit uncertainty.

To deal with these low curvature cases we use a Virtual Asteroids (VA) method. A VA is a complete orbit compatible with the available observations but by no means determined by them. A number of VA are selected, either at random [START_REF] Virtanen | Orbit computation for transneptunian objects[END_REF] or by some geometric construction [START_REF] Milani | Orbit Determination with Very Short Arcs. I Admissible Regions[END_REF], among the orbits compatible with the observations (many VA methods are described in the literature; for a recent review see [Milani 2005]). A very effective method for the problem at hand uses just one VA that is derived from the theory of the Admissible Region [START_REF] Milani | Orbit Determination with Very Short Arcs. I Admissible Regions[END_REF]].

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

The set of preliminary orbits consistent with the observations is used as a first guess for the nonlinear optimization procedure (differential corrections) that computes the nominal orbit fulfilling the principle of least squares of the observation residuals. If the preliminary orbits are well defined then at least one of them is likely to belong to the convergence domain of some least squares orbit in the differential corrections. In this case it is enough to use all the preliminary orbits as input to the differential correction routine.

If the preliminary orbits are poorly defined they may be so far from the nominal solution that the differential corrections may not converge at all or may converge to the wrong solution.

Thus, it is essential to increase the size of the convergence domain by using modified differential corrections methods. A number of these methods exist and most of them have one feature in common; the number of orbital parameters determined is less than 6. There are 4-fit methods in which 2 variables are kept fixed (e.g., [START_REF] Milani | Unbiased orbit determination for the next generation asteroid/comet surveys[END_REF]) and 5-fit or constrained least squares solutions in which one parameter is fixed (the fixed parameter need not be one of the orbital elements but may be defined in some intrinsic way [Milani et al. 2005a]).

Numerical simulations and algorithms performances

In Section 6 we describe two tests; one based on small, focused simulations with only the most difficult orbits (NEOs and TNOs), and another including all classes of solar system objects in their correct ratio (dominated by MBAs).

Both tests were performed with synthetic realizations of the Pan-STARRS et al. 2007] next generation survey. Moreover, we will show that all the innovative solutions we have proposed in this paper are essential for the accurate and reliable determination of the most difficult orbits.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT [Jedicke

Equations from the Classical Theory

In this section we provide a compact summary of the basic formulae to be used in Sections 3-5 that lead to the dynamical equation and the associated equation of degree 8 for both the Gauss and Laplace methods. We also summarize Charlier's theory on the number of solutions for Laplace's method which will be compared to the new qualitative theory of Section 4.

Laplace's Method

An observation defines a unit vector ρ = (cos δ cos α, cos δ sin α, sin δ) where (α, δ) are right ascension and declination respectively. The heliocentric position of the observed body can then be written as

r = ρ + q = ρ ρ + qq
where q is the observer's position. Let s be the arc length parameter for the path described by the relative position ρ(t) and η the proper motion, then

ds dt = η = α2 cos 2 δ + δ2 ; d ds = 1 η d dt .
We use the moving orthonormal frame [Danby 1962, Sec. 7.1]

ρ , v = d ρ ds , n = ρ × v (1)
and define the geodesic curvature κ by the equation

dv ds = -ρ + κn . (2) A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
Then the relative acceleration is

d 2 ρ dt 2 = (ρ -ρη 2 ) ρ + (ρ η + 2 ρη)v + (ρη 2 κ)n (3)
and the differential equations of relative motion are

d 2 ρ dt 2 = r -q = μq q 3 - μr r 3 (4)
with the following approximations: q = q ⊕ coincides with the center of mass of the Earth and the only force operating on both the Earth and the object at r is the gravitational attraction by the Sun, (no lunar and planetary perturbations, not even indirect perturbation by the Earth itself). From

(3)•n =(4)•n d 2 ρ dt 2 • n = ρη 2 κ = μ q q • n 1 q 3 - 1 r 3 ,
which can be presented in the form

C ρ q = 1 - q 3 r 3 where C = η 2 κq 3 μq • n (5)
and is referred to as the dynamical equation in the literature 1 . C is a nondimensional quantity that can be 0 when r = q or undetermined (of the form 0/0) in the case that the O(Δt 2 ) approximation fails, i.e. the object is at an inflection point with tangent pointing at the Sun.

Given ρ, to complete the initial conditions ρ is found from (3

)•v =(4)•v -μ q • v r 3 + μ q • v q 3 = ρ η + 2 ρη . ( 6
)
The aberration can be accurately accounted for by using as epoch of the initial conditions the true time at the asteroid t = t obsρ/c, where t obs is 1 It is actually the component of the equations of motion along the normal to the path.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

the central time at which the observation ρ has been taken and at which the derivatives have been interpolated [Crawford et al. 1930, page 99].

Eq. ( 5) is the basic formula for Laplace's method using the solution in terms of either ρ or r which are not independent quantities. From the triangle formed by the vectors q, ρ, r we have the geometric equation

r 2 = ρ 2 + 2ρq cos ε + q 2 (7)
where cos ε = q • ρ is fixed by the observation direction (ε = 180 •solar elongation). By substituting ρ from eq. ( 7) in ( 5), squaring and multiplying by C 2 r 6 (C = 0, otherwise r = q, the zero circle) we obtain

P (r) = C 2 r 8 -q 2 r 6 (1 + 2C cos ε + C 2 ) + 2q 5 r 3 (1 + C cos ε) -q 8 = 0 . ( 8
)
The trivial root r = q is due to a coordinate singularity. There can be other spurious solutions of the equation ( 8) corresponding to ρ < 0 in eq. ( 5).

Charlier's Theory

The qualitative theory of [Charlier 1910] on the number of solutions is obtained by analyzing equations ( 5)- [START_REF]Region Tracklets False Objects 3-night Overhead 2-night 1-night Oppos[END_REF]. The sign of the coefficients of eq. ( 8)

are known: -(1+2C cos ε+C 2 ) < 0 and (1+C cos ε) > 0 (see [Plummer 1918]).

Thus, there are 3 changes of sign in the sequence of coefficients and ≤ 3 positive real roots. By extracting the factor (rq) we find that P (r) = (rq) P 1 (r) ; P 1 (0) = q 7 ; P 1 (q) = q 7 C (C -3 cos ε) .

The number of solutions of the polynomial equation changes where P 1 (q)

changes sign: at C = 0 ⇔ r = q and at C -3 cos ε = 0. The latter condition

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
defines the limiting curve in heliocentric polar coordinates (r, φ) by using ρ2 = r 2 + q 2 -2 r q cos φ:

4 -3 r q cos φ = q 3 r 3 . ( 9
)
INSERT FIGURE 1 Following [Charlier 1910, Charlier 1911], the number of preliminary orbit solutions can be understood using a plot of the level curves of the function C(r, ρ) as defined by the dynamical equation, in a plane with the Sun at (0, 0), the Earth at (q, 0) and the position in each half-plane defined by the bipolar coordinates (r, ρ). The limiting curve and the zero circle can be used to determine the number of solutions for an object discovered at any point of the plane 2 . There is only one solution on the right of the unlimited branches of the limiting curve near opposition. There are two solutions for every point between the unlimited branches and the zero circle. Inside the zero circle and outside the loop of the limiting curve there is only one solution. Inside that loop there are always two solutions.

Charlier's theory assumes there is always at least one preliminary orbit solution. This results from two implicit assumptions: the value of C is measured exactly from the observations (or at least to good accuracy), and the observed object exists (not being the result of a false identification). Both assumptions may fail as discussed in Sections 5 and 6.2 respectively.

A C C E P T E D M

A N U S C R I P T ACCEPTED MANUSCRIPT

Gauss' Method

Gauss' method uses 3 observations corresponding to heliocentric positions

r i = ρ i + q i i = 1, 2, 3 (10)
at times t 1 < t 2 < t 3 with t it j = O(Δt) period and assumes coplanarity

λ 1 r 1 -r 2 + λ 3 r 3 = 0 . ( 11
)
From ( 11)×r i • ĉ, where c = r i × ṙi , the coefficients λ 1 , λ 3 are obtained as triangle ratios

λ 1 = r 2 × r 3 • ĉ r 1 × r 3 • ĉ ; λ 3 = r 1 × r 2 • ĉ r 1 × r 3 • ĉ .
From (10) and ρ1 × ρ3 •(11):

ρ 2 [ ρ1 × ρ3 • ρ2 ] = ρ1 × ρ3 • [λ 1 q 1 -q 2 + λ 3 q 3 ]. ( 12
)
Next, the differences r ir 2 are expanded in powers of t ij = t it j = O(Δt). e.g. by using the f, g series formalism r i = f i r 2 + g i ṙ2 where

f i = 1 - μ 2 t 2 i2 r 3 2 + O(Δt 3 ) , g i = t i2 1 - μ 6 t 2 i2 r 3 2 + O(Δt 4 ) . ( 13
)
Then

r i × r 2 = -g i c, r 1 × r 3 = (f 1 g 3 -f 3 g 1 )c
and

λ 1 = g 3 f 1 g 3 -f 3 g 1 > 0 ; λ 3 = -g 1 f 1 g 3 -f 3 g 1 > 0 (14) f 1 g 3 -f 3 g 1 = t 31 1 - μ 6 t 2 31 r 3 2 + O(Δt 4 ) . ( 15
)
Using ( 13) and ( 15) in ( 14) we find

λ 1 = t 32 t 31 1 + μ 6r 3 2 (t 2 31 -t 2 32 ) + O(Δt 3 ) . ( 16
) A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT λ 3 = t 21 t 31 1 + μ 6r 3 2 (t 2 31 -t 2 21 ) + O(Δt 3 ) . ( 17
)
Let V = ρ1 × ρ2 • ρ3 be 3×volume of the pyramid with vertices q, r 1 , r 2 , r 3 and substituting it and ( 16), ( 17 If the terms O(Δt 4 ) are neglected, the coefficient of the 1/r 3 2 term in (18) is

B(q 1 , q 3 ) = μ 6 t 32 t 21 ρ1 × ρ3 • [(t 31 + t 32 )q 1 + (t 31 + t 21 )q 3 ]. ( 19
)
Then multiply (18) by q 3 2 /B(q 1 , q 3 ) to obtain

- V ρ 2 t 31 B(q 1 , q 3 ) q 3 2 = q 3 2 r 3 2 + A(q 1 , q 2 , q 3 ) B(q 1 , q 3 ) (20)
where

A(q 1 , q 2 , q 3 ) = q 3 2 ρ1 × ρ3 • [t 32 q 1 -t 31 q 2 + t 21 q 3 ]. ( 21
)
Let

C 0 = V t 31 q 4 2 B(q 1 , q 3 ) , h 0 = - A(q 1 , q 2 , q 3 ) B(q 1 , q 3 )
and then, by substituting into (20),

C 0 ρ 2 q 2 = h 0 - q 3 2 r 3 2 (22)
is the dynamical equation of Gauss' method, similar (but not identical) to eq. ( 5) of Laplace's method. Using (7) at time t 2 (with q 2 , ρ 2 , r 2 and ε 2 ):

P 0 (r) = C 2 0 r 8 2 -q 2 2 r 6 2 (h 2 0 + 2C 0 h 0 cos ε 2 + C 2 0 ) + 2q 5 2 r 3 2 (h 0 + C 0 cos ε 2 ) -q 8 2 = 0 (23) A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
where the sign of the coefficients is as for [START_REF]Region Tracklets False Objects 3-night Overhead 2-night 1-night Oppos[END_REF], apart from h 0 + C 0 cos ε 2 whose sign depends upon h 0 . Note that P 0 (q) = 0 such that no root can be found analytically. The number of positive roots is still ≤ 3 but a qualitative theory such as the one of Section 2.2 is not available in the literature.

After the possible values for r 2 have been found the corresponding ρ 2 values are obtained from eq. ( 22) and the velocity ṙ2 can be computed, e.g.

from the classical formulae by Gibbs [Herrick 1971, Chap. 8].

Topocentric Gauss-Laplace Methods

There is a critical difference between the methods of Gauss and Laplace.

Gauss uses a truncation (to order O(Δt 2 )) in the motion r(t) of the asteroid but the position of the observer (be it coincident with the center of the Earth or not) are used with their exact values. On the other hand, Laplace uses a truncation to the same order of the relative motion ρ(t) (see eq.(32) in Section 5.1), thus implicitly approximating the motion of the observer. This section discusses the consequences of the difference between the techniques.

Gauss-Laplace equivalence

To directly compare the two methods let us introduce in Gauss' method the same approximation to order O(Δt 2 ) in the motion of the Earth which is still assumed to coincide with the observer. The f , g series for the Earth are

q i = 1 - μ 2 t 2 i2 q 3 2 q 2 + t i2 q2 + μ 6 t 3 i2 q 3 2 3(q 2 • q2 )q 2 q 2 2 -q2 + O(Δt 4 ). ( 24
)
By using ( 24) in ( 19) we find that 24) in ( 21)

B(q 1 , q 3 ) = μ 6 t 32 t 21 ρ1 × ρ3 • [3t 31 q 2 + t 31 (t 32 -t 21 ) q2 + O(Δt 3 )]. A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT If t 32 -t 21 = t 3 + t 1 -2t 2 = 0 (i.e., the interpolation for d 2 /dt 2 is done at the central value t 2 ) then B(q 1 , q 3 ) = μ 2 t 21 t 32 t 31 ρ1 × ρ3 • q 2 (1 + O(Δt 2 )) ; otherwise, if t 2 = (t 1 + t 3 )/2 the last factor is (1 + O(Δt)). Using (
we find

A(q 1 , q 2 , q 3 ) = - μ 2 t 21 t 32 t 31 ρ1 × ρ3 • q 2 + 1 3 (t 21 -t 32 ) 3(q 2 • q2 )q 2 q 2 2 -q2 +O(Δt 5 ). If, as above, t 32 -t 21 = t 3 + t 1 -2t 2 = 0 then A(q 1 , q 2 , q 3 ) = - μ 2 t 21 t 32 t 31 ρ1 × ρ3 • q 2 (1 + O(Δt 2 ))
and we can conclude that

h 0 = - A B = 1 + O(Δt 2 ) ; otherwise, if t 2 = (t 1 + t 3 )/2 the last factor is (1 + O(Δt)). For V we need d 2 ρ dt 2 = d ˙ρ dt = d dt (ηv) = -η 2 ρ + ηv + κη 2 n (25)
to make a Taylor expansion of ρi in

t 2 ρi = ρ2 + t i2 ηv 2 + t 2 i2 2 (-η 2 ρ2 + ηv 2 + κη 2 n2 ) + O(Δt 3 ). This implies that ρ1 × ρ3 • ρ2 = 1 2 t 12 ηv 2 × t 2 32 κ η 2 n2 -t 32 η v2 × t 2 12 κ η 2 n2 • ρ2 + O(Δt 5 )
and the O(Δt 4 ) term vanishes. Thus

V = - κη 3 2 (t 12 t 2 32 -t 32 t 2 12 ) (1 + O(Δt 2 )) = κη 3 2 t 21 t 32 t 31 (1 + O(Δt 2 )) A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT C 0 = V t 31 q 4 2 B = κη 3 t 31 q 4 2 + O(Δt 3 ) μ ρ1 × ρ3 • q 2 (1 + O(Δt)) . ( 26
)
In the denominator, ρ1 × ρ3 computed to order Δt 2 is

ρ1 × ρ3 = t 31 η n2 + t 2 32 -t 2 12 2 ( η n2 -κ η 2 v2 ) + O(Δt 3 ). ( 27
) If t 32 -t 21 = t 3 + t 1 -2t 2 = 0 then C 0 = κ η 3 t 31 q 4 2 + O(Δt 3 ) μ t 31 η q 2 q2 • n2 + O(Δt 3 ) = κ η 2 q 3 2 μ q2 • n2 (1 + (OΔt 2 )) , otherwise the last factor is (1 + O(Δt)).
Thus, neglecting the difference between topocentric and geocentric observations the coefficients of the two dynamical equations ( 5) and ( 22) are the same to zero order in Δt, to order 1 if the time t 2 is the average time.

Topocentric Laplace's Method

Now let us remove the approximation that the observer sits at the center of the Earth and introduce topocentric observations into Laplace's method. The center of mass of the Earth is at q ⊕ but the observer is at q = q ⊕ +P. Let us derive the dynamical equation by also taking into account the acceleration contained in the geocentric position of the observer P(t) such that

d 2 ρ dt 2 = - μr r 3 + μq ⊕ q 3 ⊕ -P.
Multiplying by •n and using eq. ( 3)

d 2 ρ dt 2 • n = ρη 2 κ = μ q ⊕ q⊕ • n q 3 ⊕ -q ⊕ q⊕ • n r 3 -P P • n r 3 -P • n
The term P P • n/r 3 can be neglected because P/q ⊕ ≤ 4.3 × 10 -5 which is smaller than the planetary perturbations. Thus

C ρ q ⊕ = (1 -Λ n ) - q 3 ⊕ r 3 (28) A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT where C = η 2 κq 3 ⊕ μq ⊕ • n , Λ n = q 2 ⊕ P • n μq ⊕ • n = P • n (μ/q 2 ⊕ ) q⊕ • n. (29)
Note that Λ n is singular only where C is also singular. The analog of eq. ( 6), again neglecting terms of O(p/q ⊕ ), is

ρ η + 2 ρη = μ q⊕ • v q 2 ⊕ 1 -Λ v - q 3 ⊕ r 3 , Λ v = q 2 ⊕ P • v μ q⊕ • v . ( 30
)
The important fact is that Λ n and Λ v are not small. The centripetal acceleration of the observer (towards the rotation axis of the Earth) has size

| P| = Ω 2 ⊕ R ⊕ cos θ
where Ω ⊕ is the angular velocity of the Earth's rotation, R ⊕ is the radius of the Earth and θ is the observer's latitude. The maximum of | P| 3.4 cm s -2 occurs for observers located at the equator. The quantity μ/q 2 ⊕ in the denominator of Λ n is the size of the heliocentric acceleration of the Earth, 0.6 cm s -2 . Thus |Λ n | can be > 1 and the coefficient 1 -Λ n can be very different from 1; it may even be negative. Thus, without taking into account the geocentric acceleration of the observer Laplace's geocentric classical method is not a good approximation to the topocentric general case.

However, when observations from different nights are obtained from the same station at the same sidereal time the observer's acceleration cancels out and the geocentric classical Laplace's method is a good approximation.

A common technique for the implementation of Laplace's method is to apply a negative topocentric correction to obtain the geocentric observation case. When applying this correction an initial value of ρ must be assumed as a first approximation, e.g., ρ = 1 AU [Leuschner, 1913, Page 15]. If the starting value is approximately correct then an iteration cycle will eventually

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
achieve convergence. However, if the starting value is really wrong (e.g., if

the object is undergoing a close approach to the Earth) the procedure may diverge. These reliability problems discourage the use of the classical form of Laplace's method when processing large datasets, containing discoveries of different orbital classes that span a wide range of distances.

Gauss-Laplace equivalence, topocentric

When taking into account the displacement P the Taylor expansion of q i (t) of eq. ( 24) is not applicable. We need to use

q i = q 2 + t i2 q2 + t 2 i2 2 q2 + O(Δt 3 )
where q 2 (t) and its derivatives also contain P(t). By using eq. ( 27) and assuming t 21 = t 32 , eq. ( 19) and ( 21) become

B(q 1 , q 3 ) = μ η 2 t 21 t 32 t 2 31 n2 • q 2 + O(Δt 6 ) A(q 1 , q 2 , q 3 ) = q 3 2 η 2 t 21 t 32 t 2 31 n2 • q2 + O(Δt 6 ) .
Note that q2 does not appear in A at this approximation level. Thus

h 0 = - A B = - q 3 2 n2 • q2 + O(Δt 2 ) μ n2 • q 2 + O(Δt 2 )
and once again neglecting P/q ⊕ terms we find that

h 0 = - q 3 2 n2 • q⊕2 μ n2 • q 2 - q 3 2 n2 • P2 μ n2 • q 2 + O(Δt 2 ) = = q 3 2 q 3 ⊕2 - q 3 2 n2 • P2 μ n2 • q 2 + O(Δt 2 ) n2 • q 2 = q 2 n2 • q ⊕2 q 2 + P 2 q 2 = q 2 n2 • q⊕2 + O P 2 q 2 A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT h 0 = 1 - q 3 ⊕2 n2 • P2 μ n2 • q 2 + O(Δt 2 ) + O P 2 q 2 = 1 -Λ n2 + O(Δt 2 ) + O P 2 q 2
where Λ n2 is the same quantity as Λ n of eq. ( 29) computed at t = t 2 .

The conclusion is that Gauss' method used with the heliocentric positions of the observer, q i = q ⊕i +P i , is equivalent to the topocentric implementation of Laplace's method of Section 3.2 to lowest order in Δt when neglecting very small terms O(P 2 /q 2 ).

Problems in Topocentric Laplace's Method

Contrary to common belief, Laplace's and Gauss' methods are not equivalent. Gauss' method is superior because it naturally accounts for topocentric observations by using the observer's position in eq. ( 19) and ( 21). The question then arises whether we could account for topocentric observations in Laplace's method (without iterations) by introducing the term Λ n from eq. ( 29). Surprisingly, the answer is already contained in the literature in a 100 year old paper by [Poincaré 1906, pag. 177-178].

INSERT FIGURE 2

Figure 2 shows the simulated path of an approaching NEO. The apparent motion of the asteroid from night to night cannot be approximated using parabolic segments fit to a single night3 . For the geocentric path the parabolic approximation to ρ(t) as used by Laplace would be applicable.

INSERT FIGURE 3

Figure 3 shows that topocentric observations contain information beyond what is contained in the average angles and proper motion (the attributable,

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
see Section 5). Thus, to reduce the observations to the geocenter by removing the topocentric correction is not a good strategy.

Poincaré suggested computing what we call Λ n by using a value of P obtained by interpolating the values P(t i ) at the times t i of the observations (not limited to 3, one of the advantages of Laplace's method). This method could be used but its practical advantages have not yet been established.

When the observations are performed from an artificial satellite (such as the Hubble Space Telescope or, in the future, from Gaia) the acceleration P 900 cm s -2 and the Λ n and Λ v coefficients can be up to 1, 500. A few hours of observations over several orbits can produce multiple kinks (as in [Marchi et al. 2004, Figure 1]) that contain important orbital information.

Topocentric Qualitative Theory

In Gauss' method, the dynamical equation ( 22) describes the level lines C 0 = const in a bipolar coordinate system (r 2 , ρ 2 ). In rectangular heliocentric coordinates (x, y) where the x axis is along q2 (from the Sun to the observer)

we have ρ 2 = q 2 2 + x 2 + y 2 -2xq 2 and r 2 = √ x 2 + y 2 , thus we can consider the function

C 0 (x, y) = q 2 q 2 2 + x 2 + y 2 -2xq 2 h 0 - q 3 2 (x 2 + y 2 ) 3/2 . ( 31
)
For the topocentric Laplace's method, eq. ( 31) can be used to describe C as a function of (x, y), with 1 -Λ n and q ⊕ replacing h 0 and q 2 respectively.

C 0 = 0 is the zero circle r = r 0 = q/ 3 √ h 0 for h 0 > 0 and is empty otherwise. This function tends to -∞ as (x, y) → (0, 0). There is another

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
singularity in (x, y) = (q 2 , 0); as (x, y) → (q 2 , 0) we have the following behavior: for h 0 > 1, C 0 → +∞; for h 0 < 1, C 0 → -∞; for h 0 = 1 the limit of C 0 does not exist, as shown by Figure 1. The stationary points of C 0 are the pairs (x, y) with y = 0 and x such that (h 0 |x| 3q 3 2 )x = 3q 3 2 (xq 2 ). For h 0 ≤ 0 there is only a saddle (x 1 , 0), with 0 < x 1 < q 2 (Figure 4). For h 0 > 0 there is always a saddle (x 1 , 0) with x 1 < -r 0 < 0. If 0 < h 0 < 1 there are two additional solutions, x 2 and x 3 such that 0 < x 2 < q 2 < r 0 < x 3 ;

(x 2 , 0) is a saddle, (x 3 , 0) is a maximum (Figure 5). For h 0 > 1 there is no additional stationary point (Figure 6).

INSERT FIGURE 4

The number of solutions of the dynamical equation along a fixed topocentric direction can be computed as follows. We evaluate the degree 8 polynomial (23) on the zero circle

P (r 0 ) = C 2 0 q 8 h 8/3 0 1 -h 2/3 0
and take into account that P (0) < 0. By spurious we mean a root of the polynomial (23) corresponding to ρ 2 ≤ 0 in eq. ( 22).

INSERT FIGURE 5 AND FIGURE 6

For h 0 ≤ 0 there can be either 1 or 3 positive roots4 of the polynomial equation ( 23). By comparing with eq. ( 22) these roots are all spurious if C 0 ≥ 0, all lead to a preliminary orbit solution if C 0 < 0.

For 0 < h 0 < 1 there can be either 1 or 3 positive roots of the polynomial equation, but one of them must be < r 0 and is necessarily spurious for C 0 ≥ 0.

Thus there are either 1 or 3 spurious roots for C 0 ≥ 0, either 0 or 2 for C 0 < 0.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

For h 0 > 1 there can be either 1 or 3 positive roots, one of which is > r 0 and is necessarily spurious for C 0 ≤ 0. Thus there are either 1 or 3 spurious roots for C 0 ≤ 0, either 0 or 2 for C 0 > 0.

INSERT TABLE 1

In Table 1 we summarize the possible numbers of preliminary orbit solutions, for a given direction of observation ε, in the different cases, depending upon the value of h 0 and the sign of C 0 . In this Section we are not making the assumption of Charlier, that some solutions must exist, for the reasons given in Section 2.2.

This qualitative theory generalizes Charlier's, showing that the number of solutions can be quite different, e.g., we can have 2 solutions near opposition and up to 3 at low elongation. For a qualitative theory including the generalization of the limiting curve, eq. ( 9), see [Gronchi, 2007, in preparation].

Examples INSERT FIGURE 7

We would like to find examples in which the additional solutions with respect to the classical theory by Charlier are essential, i.e. cases in which the additional preliminary orbits are closer to the true orbit leading to convergence of the least squares method, the other preliminary solutions fail.

An example in which there are two solutions when observing in a direction close to the opposition is shown in An interesting feature is that the preliminary orbit using the nearer solution has residuals of the 6 observations with RMS = 66 arcsec, while the one using the farther solution has residuals with RMS = 2.5 arcsec. In this case, if only the lowest residual preliminary solution were passed to the differential corrections step the proper solution would be discarded.

INSERT FIGURE 8

It is not easy to find a good example with 3 solutions: in many cases the solution nearest to the observer has ρ 2 too small for the heliocentric 2-body approximation to be applicable. A value for ρ 2 ≤ 0.01 AU corresponds to the sphere of influence of the Earth, i.e., the region where the "perturbation" from the Earth is actually more important than the attraction from the Sun.

Thus, a solution with such a small ρ 2 must be considered spurious because the approximation used in Gauss' and Laplace's method is not valid.

To apply our arguments on the number of solutions to a real case we used detections from the first three nights

Reliability and Precision

We need to implement the algorithms discussed in this paper for the computation of preliminary orbits in a way which is reliable when applied to the large observation data sets expected in the next generation of asteroid surveys. In doing so we need to satisfy three requirements.

The first requirement is to obtain the solutions to the polynomial equations (e.g. Equation 23) in a way which is fast and reliable in providing the number of distinct real solutions. In this way we can fully exploit our understanding of the number of solutions (with topocentric observations) as described above. This is made possible using algorithms that compute the set of roots of a polynomial equation (as a complex vector), with rigorous upper bounds for the errors including roundoff. We use the algorithm by [Bini 1996] and the corresponding public domain software 5 .

The second requirement is to improve the preliminary orbit as obtained from the solutions of the degree 8 polynomial equations in such a way that it is as close as possible to the least squares solution to be later obtained by differential corrections. There is such an immense literature on this topic that we will not even attempt to provide a set of references.

Conceptually, as shown by [START_REF] Celletti | Four Classical Methods for Determining Planetary Elliptic Elements: A Comparison[END_REF], each step in the iterative procedures used to improve the preliminary orbits (which they call Gauss map 6 ) can be shown to increase the order in Δt of the approxima-

5 For the Fortran 77 version visit http://www.netlib.org/numeralgo/na10 while the Fortran 90 version is available at http://users.bigpond.net.au/amiller/pzeros.f90.

6 The classical treatises, such as [START_REF] Crawford | Determination of orbits of comets and asteroids[END_REF], use the term differential corrections for algorithms of the same class as Gauss map in [START_REF] Celletti | Four Classical Methods for Determining Planetary Elliptic Elements: A Comparison[END_REF]. We follow the terminology of the recent papers because, in modern usage, differential correc-

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
tion to the exact solutions of the 2-body equations of motion. However, [START_REF] Celletti | Dependence on the observational time intervals and domain of convergence of orbital determination methods[END_REF] have also shown that iteration of a Gauss map can diverge when the solution of the degree 8 equation is far from the fixed point of the iterative procedure, outside of its convergence domain.

The same results apply to algorithms, such as those of [Leuschner, 1913] and [START_REF] Crawford | Determination of orbits of comets and asteroids[END_REF], that improve Laplace's preliminary orbit. The difference is that in Laplace's method the first approximation is with the observations treated as geocentric (or possibly corrected with an assumed distance, [Leuschner, 1913, page 15]), while in Gauss' method (Merton's version) the first approximation properly handles topocentric observations. As mentioned earlier, this leads us to prefer Gauss' method.

We have implemented an iterative improvement algorithm for Gauss' method and found that in most cases it provides a preliminary orbit much closer to the least squares solution which is therefore a more reliable first guess for the least squares algorithms. We have found that the Gauss map diverges in a small fraction of test cases, but this behavior occurs often enough to significantly decrease the efficiency of the algorithm (see Section 6.1). In some cases the number of orbits to which the Gauss map converges is less than the number of solutions of the degree 8 equations. It can happen that one of the lost degree 8 solutions was the only one leading to a least squares solution. One method to obtain the highest efficiency without an inordinate increase in the computational cost is to run two iterations, with and without the Gauss map. In the second iteration we also accept preliminary orbits tions refers to the iterative method to solve the least squares problem.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

with comparatively large residuals (up to a RMS of 100 arcsec) to allow for significant perturbations by a third body.

The third requirement is to use modified differential corrections algorithms, with larger convergence domains, in such a way that even when the geodetic curvature and the coefficients C and C 0 of the two methods are poorly constrained by the available observations (because the arc length on the celestial sphere is too short) the rough preliminary orbit can lead to a least squares solution. This possibility is discussed in the next section.

Weak preliminary orbits

An essential difference between the classical works on preliminary orbits and the modern approach to the same problem is that the effects of the astrometric errors cannot be neglected. Since the next generation of all-sky asteroid surveys will acquire fewer observations of the objects the deviations of the observed path from a great circle may not be significant.

Uncertainty of Curvature

The explicit computation of the two components of curvature of interest for orbit determination, geodesic curvature κ and along track acceleration η, can be performed by using the properties of the orthonormal frame (1) in straightforward computation using the Riemannian structure of the unit sphere [Milani et al. 2007a, Section 6.4]. The results are

κ = 1 η 3 ( δ α -α δ) cos δ + α η 2 + ( δ) 2 sin δ = κ(α, δ, α, δ, α, δ)(32) η = 1 η α α cos 2 δ + δ δ -( α) 2 δ cos δ sin δ = η(α, δ, α, δ, α, δ) . (33) A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
Given these explicit formulae it is possible to compute the covariance matrix of the quantities (κ, η) by propagation of the covariance matrix of the angles and their derivatives with the matrix of partial derivatives for κ and η

Γ κ, η = ∂(κ, η) ∂(α, δ, α, δ, α, δ) Γ α,δ ∂(κ, η) ∂(α, δ, α, δ, α, δ) T . ( 34
)
The covariance matrix Γ α,δ for the angles and their first and second derivatives is obtained by the procedure of least squares in a fit to the individual observations as a quadratic function of time. The partials of κ and η are given below (note that the partials with respect to α are zero).

∂κ ∂δ = - 1 η 5 -2 α3 cos 2 δ sin δ δ + sin δ δ α δ2 + 2 α2 cos 2 δ sin δ α δ - -sin δ α δ3 -α5 cos 3 δ -4 α3 cos δ δ2 + α3 cos 3 δ δ2 -2 α cos δ δ4 ∂ η ∂δ = - α 2 η 3 sin(2δ) α2 α cos 2 δ + 2 δ2 α -α δ δ +2 α δ3 cos(2δ)+2 α3 δ cos 4 δ ∂κ ∂ α = 1 η 5 -α cos 3 δ 2 α δ -3 δ α + δ2 δ cos δ -α2 cos 2 δ sin δ + 2 δ2 sin δ ∂ η ∂ α = - cos δ δ η 3 -cos δ α δ + α3 sin δ cos 2 δ + 2 α sin δ δ2 + cos δ δ α ∂κ ∂ δ = - 1 η 5 cos δ α2 α cos 2 δ -2 δ2 α + 3 α δ δ -α δ sin δ α2 cos 2 δ -2 δ2 ∂ η ∂ δ = - α cos 2 δ η 3 -δ α + α3 cos δ sin δ + α δ ∂κ ∂ α = - δ cos δ η 3 ; ∂κ ∂ δ = α cos δ η 3 ; ∂ η ∂ α = α cos 2 δ η ; ∂ η ∂ δ = δ η . ( 35
)
The last four of these partials, the 2 × 2 matrix ∂(κ, η)/∂(α, δ), contribute to the principal part of the covariance of (κ, η) for short arcs (see below).

We use a full computation of the covariance matrix without approximations to assess the significance of curvature by using the formula from

A C C E P T E D M A N U S C R I P T
ACCEPTED MANUSCRIPT [Milani et al. 2007a] providing

χ 2 = κ η T Γ -1 κ, η κ η ( 36
)
and we assume that the curvature is significant if χ 2 > χ 2 min = 9.

The Infinite Distance Limit

The problem of low values of C can occur in two ways: near the zero circle and for large values of both ρ and r. On the other hand, the uncertainty in the estimates of the deviations from a great circle will depend upon the length of the observed arc (both in time Δt and in arc length ∼ η Δt). For short observed arcs it may be the case that the curvature is not significant.

Then the preliminary orbit algorithms will yield orbits which may fail as starting guesses for differential corrections.

We will now focus on the case of distant objects. We would like to estimate the magnitude of the uncertainty in the computed orbit with respect to the small parameters ν, τ, b where ν is the astrometric accuracy of the individual observations (in radians) and τ = n ⊕ Δt, b = q ⊕ /ρ are small for short observed arcs and for distant objects respectively. Note that the proper motion η for b → 0 has principal part n ⊕ b -the effect of the motion of the Earth. The uncertainty in the angles (α, δ) and their derivatives can be estimated as follows (see [Crawford et al. 1930, page 68])

Γ α,δ = O(ν) , Γ α, δ = O(ντ -1 ) , Γ α, δ = O(ντ -2 ) .
The uncertainty of the curvature components (κ, η) should be estimated by the propagation formula (34) but it can be shown that the uncertainty of

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

(δ, α, δ) contributes with lower order terms. Thus we use the estimates based upon the partials in the last line of ( 35)

∂(κ, η) ∂(α, δ) = O(b -2 ) n -2 ⊕ O(b -2 ) n -2 ⊕ O(1) O(1)
and obtain

Γ κ, η = ν O(b -4 τ -2 ) O(b -2 τ -2 )n 2 ⊕ O(b -2 τ -2 )n 2 ⊕ O(τ -2 )n 4 ⊕ .
To propagate the covariance to the variables (ρ, ρ) we use the implicit equation connecting C and ρ obtained by eliminating r from ( 7) and ( 22):

F (C, ρ) = C ρ q ⊕ + q 3 ⊕ (q 2 ⊕ + ρ 2 + 2q ⊕ ρ cos ε) 3/2 -1 + Λ n = 0 . ( 37
)
For b → 0 we have C b -1 → 1; thus C → 0 and is of the same order as b.

Although C depends upon all the variables (α, δ, α, δ, α, δ), its uncertainty mostly depends upon the uncertainty of κ and thus, ultimately, upon the difficulty in estimating the second derivatives of the angles.

Next, we compute the dependence of Γ ρ, ρ upon Γ κ, η. From the derivatives of the implicit function ρ(κ), assuming cos ε, η, n to be constant and keeping only the term of lowest order in q/ρ, we find

∂ρ ∂κ = - η 2 q 4 μ q⊕ • n ρ q ⊕ C + O q 3 ρ 3 = q ⊕ O(1) .
In the same way from (30) we deduce η = n 2 ⊕ O(b) and obtain the estimates

∂ ρ ∂κ = n ⊕ q ⊕ O(1) , ∂ ρ ∂ η = q ⊕ n ⊕ O(b -2 ) .
For the covariance matrix,

Γ ρ, ρ = ∂(ρ, ρ) ∂(κ, η) Γ κ, η ∂(ρ, ρ) ∂(κ, η) T , A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
we compute the main terms of highest order in b -1 , τ -1 as

Γ ρ, ρ = ν b -3 τ -2 q 2 ⊕ O(1) q 2 ⊕ n ⊕ O(1) q 2 ⊕ n ⊕ O(1) q 2 ⊕ n 2 ⊕ O(1) . ( 38
)
In conclusion, if (ρ, ρ) are measured in the appropriate units (AU for ρ and n ⊕ AU for ρ) their uncertainties are of the same order.

This conclusion appears different from [START_REF] Bernstein | Orbit Fitting and uncertainties for Kuiper belt objects[END_REF] who claim that for a TNO arc with low curvature the inverse distance 1/ρ can be determined in a robust way while the other variable ρ/ρ remains essentially undetermined. In fact, by propagating the covariance from eq. ( 38) to the variables (1/ρ, ρ/ρ) and expressing them in the natural units 1/q ⊕ , n ⊕ we find that the RMS of ρ/ρ is larger by a factor 1/b than the one of 1/ρ; thus there is no disagreement.

The coordinates (ρ, ρ) together with (α, δ, α, δ) form a set of Attributable Orbital Elements with the special property that the confidence region of solutions with low residuals is a very thin neighborhood of a portion of the (ρ, ρ) plane [Milani et al. 2005b, Section 3]. A similar property, but with a different plane, is shared by the Cartesian Elements. Thus these coordinates are very suitable for differential corrections when performed under conditions of quasi-linearity even for large corrections. A set of coordinates containing

(1/ρ, ρ/ρ) results in a much larger nonlinearity with corresponding increased risk of divergence.

We do agree with [START_REF] Bernstein | Orbit Fitting and uncertainties for Kuiper belt objects[END_REF]] that for a TNO observed over an arc shorter than one month there is very often an approximate degeneracy forcing the use of a constrained orbit (with only 5 free parameters). The weak direction, along which an arbitrary choice needs

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
to be made, is in the (ρ, ρ) plane, may vary and is generally not close to the ρ axis [Milani et al. 2005b, Figures 3-6]. In the context of the tests with TNO orbits described in Section 6.1, for simulated discoveries around opposition we found that the weak direction forms an angle with the ρ axis (computed with the scaling indicated by eq. ( 38)) between -31 • and +17

• while near quadrature it forms an angle with the ρ axis between -54 • and +36 • . Thus, the weak direction depends strongly upon the elongation ( [START_REF] Bernstein | Orbit Fitting and uncertainties for Kuiper belt objects[END_REF] warn that their arguments are not applicable exactly at opposition).

From Preliminary to Least Square Orbits

The procedure to compute an orbit given an observed arc with ≥ 3 nights of data (believed to belong to the same object) begins with the solution of the degree 8 equation ( 23) and ends with the differential corrections to achieve a least squares orbit with 6 solved parameters. For algorithms more efficient than the classical ones we consider up to four intermediate steps:

1. an iterative Gauss map to improve the solution of the degree 8 equation as discussed in Section 4.2;

2. adding to the preliminary orbit(s) another one, obtained from the Attributable and a value for (ρ, ρ) selected inside the Admissible Region; Intermediate step 1 has been discussed in Section 4.2.

By Attributable we mean the set of 4 variables (α, δ, α, δ) estimated at some reference time by a fit to the observations [START_REF] Milani | The Asteroid Identification Problem IV: Attributions[END_REF]. It is possible to complete an attributable to a set of orbital elements by adding the values of range and range rate (ρ, ρ) at the same time. For each attributable we can determine an Admissible Region which is a compact set in the (ρ, ρ) plane compatible with Solar System orbits [START_REF] Milani | Orbit Determination with Very Short Arcs. I Admissible Regions[END_REF].

For intermediate step 2 we distinguish two cases depending upon the topology of the Admissible Region. If it has two connected components (this occurs for distant objects observed near opposition) we select the center of symmetry of the component far from the observer. This corresponds to an orbit with 0 ≤ e < 1; note that, sometimes, a circular orbit may be incompatible with the Attributable.

If the Admissible Region is connected then we select the point along the symmetry line ρ = const at 0.8 times the maximum distance ρ compatible with e ≤ 1. This case always occurs near quadrature; if the object is distant, thus has a low proper motion η, the selected point is also far.

The selected point (ρ, ρ) in the Admissible Region completed with the Attributable provides a compatible orbit belonging to the Solar System; this is called a Virtual Asteroid (VA) [Milani 2005]. This VA method provides an additional preliminary orbit. We shall see in Section 6.1 that for TNOs this additional preliminary orbit is often required, in most cases near quadrature,

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
because the curvature is hardly significant.

Intermediate step 3 is essentially the method proposed by D. Tholen, available in his public domain software KNOBS. It has already been tested in the context of a simulation of a next generation survey in [START_REF] Milani | Unbiased orbit determination for the next generation asteroid/comet surveys[END_REF]].

Intermediate step 4 is fully described in [Milani et al. 2005a]. Our preferred options are to use either Cartesian or Attributable Elements scaled as described in [Milani et al. 2005a, Table 1], that is, consistently with eq. ( 38).

The steps listed above are all optional and indeed it is possible to compute good orbits in many cases without some of them. However, the steps must be linked in a suitable manner, to provide a reliable algorithm and a least squares orbit. As an example, step 1 may be used in a first iteration but omitted in a second one.

Step 2 is essential for distant objects while step 3 is used whenever the curvature is insignificant (i.e., when the observed arc is of type 1 [Milani et al. 2007a]), as determined with eq. ( 36).

Step 4 is important for weakly determined orbits where the differential corrections may diverge when starting from an initial guess with comparatively large residuals. Even step 4 may fail and cause the differential corrections to diverge. In this case the differential corrections are restarted using the outcome of the previous step. This connecting logic is an extension of the one presented in [Milani et al. 2005a, Figure 5].

Tests

We have performed a series of tests of our algorithms using a realistic Solar System Model (a catalog of orbits for synthetic objects [START_REF] Milani | Unbiased orbit determination for the next generation asteroid/comet surveys[END_REF]])

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
and a simulation of the performance of one of the next generation surveys: Pan-STARRS [Hodapp et al.(2004)]. We employ a realistic observation scheduler and instrument performance and identify which of the synthetic objects have detections above a threshold signal to noise ratio. We then add S/N-dependent astrometric error to the detections at the level expected for the Pan-STARRS survey (about 0.1 arcsec). We have not included false detections (corresponding to no synthetic object).

Then we assemble detections from the same observing night which could belong to the same object into tracklets (based on the angular separation and morphology of the detections). Tracklets from at least three distinct nights are then assembled into tracks. For these simulations we have used the algorithms of [START_REF] Ivezić | LSST: Comprehensive NEO Detection[END_REF]] to assemble both tracklets and tracks.

When the number density of detections per unit area is low both tracklets and tracks are (almost always) true, i.e., they contain only detections of one and the same synthetic object. When the number density is large, as expected for the next generation surveys, both tracklets and tracks can be false, i.e., containing detections belonging to different objects and/or false detections. This is why tracks need to be confirmed by computing an orbit: first a preliminary orbit, then by differential corrections another orbit which fits all the observations in the least squares sense. The structure containing the track and the derived orbit with the accessory data for quality control (covariance, weights and residuals, statistical tests) is called identification [Milani et al. 2007b].

It is important to note that we are not testing the performance of any

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
particular next generation survey (e.g. Pan-STARRS or LSST). Instead, the purpose of these tests is to measure the performance of the algorithms described in this paper according to the following criteria:

• Efficiency E: the fraction of true tracks for which good preliminary and least squares orbits were calculated.

• Accuracy A: the fraction of returned orbits that correspond to true tracks. I.e., the orbit computation should fail on false tracks (either no preliminary orbit or no least squares orbit or the residuals for the derived orbit are too large).

• Goodness G: the fraction of least squares orbits close enough to the ground truth orbits to allow later recovery (e.g., in another lunation).

A speed criterion (based on CPU time) is less important because computing power grows as fast as the astrometric data7 . Still, we need to confirm that the very large data sets expected from the next generation surveys can be processed with existing and reasonable computational resources.

Small targeted tests

Since the orbits of MBAs and Jupiter Trojans are easier to compute than those of NEOs and more distant objects [START_REF] Milani | Unbiased orbit determination for the next generation asteroid/comet surveys[END_REF] we have prepared four targeted simulations: two containing only observations of NEOs and two with TNOs only. In both cases, one of the simulations covers the area near opposition and the other covers the so called sweet spots at solar elongations

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
between 60 • and 90 • . We will see that the most relevant metric is Efficiency.

Accuracy is not an issue because the number density per unit area on the sky is small (indeed, Accuracy is 100% in all tests within this Subsection).

INSERT TABLE 2

In the NEO simulations (Table 2) we have obtained very high Efficiency.

It could be improved with increased computational intensity but this could impair the Accuracy for larger data sets.

All of the cases for which an orbit was not returned, even though a true track was proposed, resulted from a failure of the preliminary orbit determination, in most cases because the degree 8 equation had only spurious roots (in the sense of Section 4.1). In other cases, some useful preliminary orbits were discarded because the RMS of the fit was large, often as high as 200

to 300 arcsec. The VA method was not of any help, as expected, since it is intended for low curvature cases.

INSERT TABLE 3

Table 3 provides the same information as Table 2 but for TNOs. For the distant objects the preliminary orbit algorithms did not suffer a single case of failure. The very few true tracks without orbit are due to failure in quality control (described in Section 6.2).

To assess the proportion of this success that is due to the Virtual Asteroid method (which is expected to be especially effective for the low curvatures typical of TNOs) we have rerun the simulation without the VA method. The results are provided in the "No VA" column of Table 4 where it is clear that giving up the VA method would result in a significant loss of TNO discoveries

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
at opposition and even more so in the sweet spots. The conclusion is that the VA method is essential for TNOs while it is almost irrelevant for NEOs.

INSERT TABLE 4

Did the effort in reliably handling double (even triple) preliminary orbit solutions significantly improve the performance in the NEO case? We have rerun the simulations with only one preliminary orbit passed to differential corrections -the one with the lowest RMS residuals. The results (Table 4, column "1 Pre") clearly show that for NEOs near quadrature passing all the possible preliminary orbits to the differential corrections procedure is essential for maximum Efficiency.

Another test has been to truncate the algorithm after the first of the two iterations (see Section 5.3), the one with a tighter control in the RMS of the residuals for the 2-body preliminary orbit (set at 10 arcsec in these tests) and using the Gauss map. The results (column "1st It.") are that the second iteration has no effect on TNOs but is relevant for NEOs especially in the sweet spots.

We also assess how much the improved differential corrections (discussed in Section 5.3) have contributed to the success of these simulations of orbit determination. The "No 4fit" column of Table 4 provides the results when the the 4-parameter fit step was not used. The results, essentially identical to the "No VA" case, indicate that the two algorithms must be used together for TNOs. Similarly, the "No LOV" columns indicates that the step with the 5-parameter least squares fit to obtain a LOV solution has a very large effect for TNOs. We get the worst results ("LSQ" column) if neither the

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
4-fit algorithm nor the LOV solutions are used and the preliminary orbits are passed directly to a full 6-parameters differential corrections.

The column labeled "Best" in Table 4 refers to the best combination of innovative and improved algorithms that we have identified. A comparison of the "Best" column with the other columns indicates that all the steps discussed in Sections 4 and 5 are essential to achieve the best results for both NEOs and TNOs.

Large scale tests

The main purpose of a large scale simulation is to measure the Accuracy.

However, Efficency and Accuracy are not independent. When there are Discordant Identifications (with some tracklets in common between different derived objects) and the orbits cannot be merged into one with all the tracklets of both (i.e. the same object has been identified twice), there is no way to choose which of the two is true. In this case we might discard both identifications which results in losing true identifications and decreasing the Efficiency. By keeping both we would decrease Accuracy while maintaining the Efficiency. Regardless of our choice, each false identification introduces permanent damage to the quality of the results8 .

Thus, with the aim of measuring the Accuracy of our algorithms we have prepared simulations for one lunation of a next generation survey both near opposition and the sweet spots. The assumed limiting magnitude was V=24

and the Solar System model was used at full density, including the over-

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
whelming majority of MBAs (11 million synthetic objects, 10 M MBAs, 269 K NEOs, 28 K TNOs). Table 5 gives the size of the realized synthetic dataset after the simulation. While the focus of this paper is those objects that have tracklets available on three nights, those that are observed for less nights are part of the problem because their tracklets can be incorrectly linked to other objects [Milani et al. 2006, figure 3].

INSERT TABLE 5

The first Accuracy problem occurs at the tracklet composition stage as shown in Table 5. Some tracklets are false because they are composed of detections belonging to different objects, but detections may appear in different tracklets so the true tracklet has probably been identified. The question is then whether the false tracklets are incorporated into the final accepted orbits. We expect that introducing false detections into the simulation will not cause an explosion of false tracklets because the noise is random in location and therefore unlikely to be spatially correlated on an image.

The second Accuracy problem occurs at the track composition stage. A track is just a hypothesis of identification to be checked by computing an orbit: at a high tracklet number density most of the tracks are false. The Overhead is the ratio between the total number of proposed and true tracks

In our simulations the Overhead was large as shown in Table 5. In the sweet spots the Overhead exceeded what was found in previous simulations [Kubica et al. 2007, Table 3].

The question is whether the orbit determination stage can produce true orbits with high Efficiency and still reject almost all the false tracks (high The most effective method to remove false tracks is a global consideration of all identifications (derived orbits). We have previously defined the normalization of lists of identifications in [Milani et al. 2005b, Section 7] and [Milani et al. 2006, Section 6]. The process removes duplications and inferior identifications but also rejects all the Discordant Identifications. This is not because they are all presumed false, indeed very often one true and one false identification are Discordant, but we do not know which is which unless one of the two has a significantly better fit. If the difference in the normalized RMS of the astrometric residuals is more than 0.25 we keep the best; otherwise we remove both and sacrifice Efficiency for Accuracy.

INSERT TABLE 7

The results of the normalization procedure are shown on the right hand side of Table 6. The false tracks can be reduced to a negligible number

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
but the Efficiency decreases as a result of the normalization as shown in the difference between the columns labeled "Eff.%" and "Eff.No.%" in Table 7.

Table 7 provides the efficiency of our algorithms as a function of orbital class. The Efficiency for NEOs and TNOs is not affected by normalization because of the lower sky-plane density of objects with their characteristic rates of motion even when embedded in a full-scale solar system population.

On the other hand, the confusion among objects with main belt rates of motion can be high and this causes our algorithm to lose a few percent of MBAs.

Nevertheless, even this problem can be solved together with recovering the few lost NEOs and comets as described below.

An analysis of the Efficiency for each of the three separate steps (track composition, orbit computation, normalization) reveals that the algorithm to generate tracks is 97.6% and 98.7% efficient at opposition and at the sweet spots, respectively, the orbit computation procedure on the proposed true tracks achieves 99.8% and 99.3% efficiency, and the normalization procedure is 98.6% and 99.4% efficient. Thus, the performance of each of the three steps is well balanced and there is not much room for improvement10 . The solution is to use a two iteration procedure.

The normalization procedure generates two outputs: the new list of identifications and the list of leftover tracklets which have not been used in the confirmed identifications. When two tracklets have detections in common, if one of the tracklets is included in a confirmed identification then the other can be ignored. Thus the set of tracklets is sharply reduced after normalization

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
which simplifies further processing. Table 8 shows that the normalization is effective in discarding false tracklets confirming with a full scale simulation what had been found with the small simulation of [START_REF] Milani | Unbiased orbit determination for the next generation asteroid/comet surveys[END_REF]].

INSERT TABLE 8

The remaining tracklets after normalization can be used as input to another iteration with different controls, perhaps using tighter requirements on residuals combined with looser thresholds on forming tracks in an effort to identify objects on hyperbolic trajectories. In subsequent iterations

Accuracy should be less of a problem because of the reduced sky-plane density of tracklets. Of course, new algorithms could be implemented in successive iterations over the tracklet data set. e.g., [Milani et al. 2005b] and [START_REF] Boattini | Low solar elongation searches for NEOs: a deep sky test and its implications for survey strategies[END_REF]]. There are only practical limits to how clever an algorithm might be and how many iterations could take place. However, to

show that the normalized Efficiency values from Table 7 are not a problem we have run an improved version of the recursive attribution algorithms of [Milani et al. 2005b] on the leftover tracklets. To control the false identifications we have used even tighter quality controls. The results are provided in Table 9 showing an almost complete recovery of the orbits that were not identified in the previous iteration.

INSERT TABLE 9

A byproduct of the procedure outlined above is the computation of normalized identifications for tracks composed of tracklets in only two nights (recall that the entire discussion above involved only tracks containing three nights of tracklets). The Efficiency and Accuracy for the 2-night linking and

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
orbit determination procedure is provided in Table 10. After this step the tracklets remaining are mostly from objects observed in just one night.

INSERT TABLE 10

Up to this point we have not assessed the quality or Goodness (Section 6) of the orbits obtained by our procedures. It is not a simple quantity to parameterize but one practical measure is to determine whether the results from one lunation can be used to attribute detections in the next (or previous) lunation. Towards this end we have run two simulations of surveys at opposition for consecutive lunations. Given the 3-night identifications for the first month we attempt to attribute to them the corresponding tracklets in the next month. The process was 99.6%, 99.7% and 99.9% efficient for objects with 1, 2 and 3 tracklets in the second lunation, respectively. There were no NEOs among the few cases of failed/incomplete attribution!

Conclusions and Future Work

The purpose of this paper is to identify efficient algorithms to compute preliminary and least squares orbits given a set of detections in a "track" (or proposed identification).

We have developed efficient and accurate algorithms by revising the classical preliminary orbit methods. The most important improvements are provisions to keep alternate solutions under control. The existence of double solutions has been known for a long time and we have shown that even triple solutions can occur. Still there is no reason this should impair the orbit determination performance.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

For the differential corrections stage that provides a least squares fit of the orbit to the detections using the preliminary orbit as the first step in the iteration, we adopted algorithms available from previous work (by ourselves and others). When the algorithms are combined with suitable control logic they significantly improve the efficiency of differential corrections even when the preliminary orbits are not close to the nominal solutions.

The third stage of orbit determination process is quality control based upon statistical analysis of the residuals. When there are only a small number of objects this may be unnecessary but, with the high detection density expected with the next generation surveys, quality control will be critical because tracklets belonging to different objects may be incorrectly identified.

Removing false identifications is not easy. We have found the method of normalization to be very effective for this purpose but, unavoidably, some true identifications are sacrificed to remove the discordant false ones. We need to select options and details of the algorithms such that the number of false identifications is kept low while the true identifications are not lost.

Although our mathematically rigorous theoretical results do not need confirmation it has been useful to test their practical performance on simulations of the next generation surveys. In this way we have shown that orbits can be computed even for the most difficult classes of orbits. We have also shown, with full density simulations including an overwhelming majority of MBAs, that the large number of objects does not result in a "false identification catastrophe". On the contrary, a large number density is compatible with a low number of lost objects provided the quality control on the residuals is Prelim. Roots Spurious h 0 ≤ 0 C 0 < 0 1 or 3 1 or 3 0 C 0 ≥ 0 0 1 or 3 1 or 3 0 < h 0 < 1 C 0 < 0 1 or 3 1 or 3 0 or 2 C 0 ≥ 0 0 or 2 1 or 3 1 or 3 h 0 > 1 C 0 ≤ 0 0 or 2 1 or 3 1 or 3 C 0 > 0 1 or 3 1 or 3 0 or 2

  like NEOs and TNOs.

Figure 7 .A

 7 The half line in the observing direction has two intersections with the C 0 = 0.4 level curve. The intersection point closest to the Earth leads to a useful preliminary orbit and has a 's theory with ρ 2 = 0. The more distant intersection point, at 2.2 AU from the Earth, leads to a preliminary orbit with e 10.

  (9, 11 and 12 January 2002) of observation of the asteroid 2002 AA 29 . The observations at an elongation of 111 • yield values of C 0 = 1.653 and h 0 = 1.025 such that there is only one solution with ρ 2 = 0.045 (see Figure8, left) which leads to a least squares solution with ρ 2 = 0.044. Although the value of h 0 is not very far from 1 the existence of the solution depends critically on h 0 -1 = 0. If the value of h 0 had been set to 1 there would be no solution as shown in Figure8(right).

3

  . a fit of the available observations to a 4-parameter attributable; the values of ρ and ρ are kept fixed at the previous values; 4. a fit of the available observations constrained to the Line Of Varia-), a smooth curve defined by minimization on hyperplanes orthogonal to the weak direction of the normal matrix.

  much smaller number of false tracks contain false tracklets("F.Tr." columns), that is, even the presence of a significant fraction of false tracklets affects neither Efficiency nor Accuracy.The false identifications result from combining tracklets from 2 (or 3) distinct simulated objects. With a fit passing all the quality controls we cannot a priori discard any of them: only by consulting the ground truth we know they are false. By further tightening the quality control parameters we may remove many false but also some true identifications. The values of the metric controls that we used are the result of adjustment suggested by experiments to find an acceptable balance between Accuracy and Efficiency for real survey operations.

Figure 1 :Figure 3 :Figure 5 :Figure 6 :

 1356 Figure1: Level curves of C in rectangular coordinates (solid lines) including the limiting curve (labeled) and the zero circle (dashed). For a given value of C and an observation direction (dotted) there can be either 1 or 2 solutions, e.g., for C = 0.3 there are 2.

  32 t 21 (t 31 + t 32 )q 1 + t 32 t 21 (t 31 + t 21 )q 3 ] + O(Δt 4 ) .

			) in (12), then using t 2 31 -t 2 32 = t 21 (t 31 + t 32 )
	and t 2 31 -t 2 21 = t 32 (t 31 + t 21 ) we find
		-V ρ 2 t 31 = ρ1 × ρ3 • (t 32 q 1 -t 31 q 2 + t 21 q 3 ) +	(18)
	+ ρ1 × ρ3 •	μ 2 6r 3	[t

Table 1 :

 1 Number of preliminary orbit solutions for different ranges of h 0 and C 0 . Prelim: the number of preliminary orbit solutions. Roots: the number of positive roots of the polynomial equation (23). Spurious: the number of spurious roots.

Table 2 :

 2 Performance characteristics for the NEO simulations. For each of the NEO simulations, and separately for objects observed on a different number of nights, the columns give the: [1] Total number of objects, [2] Number of Complete Identifications (containing all the tracklets belonging to the object), [3] Efficiency (defined as [2]/[1], in %), [4] Number of incomplete Identifications, [5] Fraction [4]/[1] in % of incomplete Identifications, [6] Number of objects lost (no confirmed Identification), [7] Fraction [6]/[1] in % lost.

		[1]	[2]	[3]	[4]	[5]	[6]	[7]
	Observed					Inc.		Lost
	Arc	Total Compl. Effic. Inc. Fraction Lost Fraction
	Opposition
	3-nighters	1123	1119 99.6%	0	0.0%	4	0.4%
	Sweet Spots
	3-nighters	397	389 98.0%	0	0.0%	8	2.0%

Table 3 :

 3 Performance characteristics for the TNO simulations. See Table2for column definitions.

		[1]	[2]	[3]	[4]	[5]	[6]	[7]
	Observed					Inc.		Lost
	Arc	Total Compl. Effic. Inc. Fraction Lost Fraction
	Opposition
	3-nighters	2005	2001 99.8%	3	0.15%	1	0.05%
	Sweet Spots
	3-nighters	2493	2491 99.9%	0	0.00%	2	0.08%

Table 4 :

 4 Fraction of lost 3-nighters using different algorithms. The methods used for each column are detailed in the text.

	Simulation Best 1 Pre 1st It. No VA No 4fit No LOV LSQ
	NEO Opp. 0.40% 0.50%	3.2%	0.4%	0.4%	0.4% 0.4%
	NEO Sw.	2.00% 13.4% 18.1%	2.3%	2.3%	2.8% 2.8%
	TNO Opp. 0.05% 0.05% 0.05% 38.3% 38.4%	45.7% 47.7%
	TNO Sw.	0.10% 0.10% 0.10% 75.3% 75.3%	82.6% 82.8%

Table 5 :

 5 Characteristics of the data sets for the full solar system model simulations. The columns provide the [1] survey region, [2] number of tracklets, [3] number of false tracklets, [4] number of simulated objects with observed tracklets, [5] number of simulated objects with observed tracklets on 3 different nights, [6] overhead (see text, the ratio of false to real tracks), [7] number of objects with tracklets in 2 different nights, [8]number of tracklets in only 1 night.

This plane does not correspond to a physical plane in that it also describes the points outside the ecliptic plane.

Our translation of Poincaré: It is necessary to avoid computing these quantities by starting from the law of rotation of the Earth.

All roots are counted with multiplicity; e.g., in this case a double root can occur.

Moore's law tells us that the number of elements on a chip grows exponentially with time: this applies equally to the number of pixels on a CCD and CPU speed.

False facts are highly injurious to the progress of science, for they often endure long..., C. Darwin, The Origin ofMan, 1871.

Loosening the controls for track composition would improve the Efficiency at the expense of increasing the false identifications and the losses at the normalization stage.

Acknowledgments

Milani & Gronchi are supported by the Italian Space Agency through contract 2007-XXXX, Knežević from Ministry of Science of Serbia through project 146004 "Dynamics of Celestial Bodies, Systems and Populations". Jedicke & Denneau are supported by the Panoramic Survey Telescope and Rapid Response System at the University of Hawaii's Institute for Astron-

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

Accuracy). To achieve this goal the residuals of the best fit orbits need to be submitted to a rigorous statistical quality control. Our residuals quality control algorithm uses the following 10 metrics (control values in square brackets):

• RMS of astrometric residuals divided by the assumed RMS of the observation errors (=0.1 arcsec in these simulations) [1.0] • RMS of photometric residuals in magnitudes [0.5] • bias of the residuals in RA and in DEC [1.5] • first derivatives of the residuals in RA and DEC [1.5] • second derivatives of the residuals in RA and DEC [1.5] • third derivatives of the residuals in RA and DEC [1.5] To compute the bias and derivatives of the residuals we fit them to a polynomial of degree 3 and divide the coefficients by their standard deviation as obtained from the covariance matrix of the fit 9 .

INSERT TABLE 6

The results are summarized in Tables 6 and7. Notwithstanding the tight quality controls on residuals, while processing tens of millions of proposed tracks a few thousand false tracks are found to fit all their observations well ("False" columns). The numbers are small with respect to the total number of tracks but they are not negligible as a fraction of the true tracks ("%"

9 When these algorithms are used on real data additional metrics should assess the outcome of outlier removal [START_REF] Carpino | Error Statistics of Asteroid Optical Astrometric Observations[END_REF]. For simulations this does not apply.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

tight enough and the sequence of algorithms is suitably chosen.

The performance for identification and orbit determination critically depends upon the individual algorithms and upon the pipeline design -the sequence of algorithms operating one upon the output of the last. We have used the algorithms from [START_REF] Ivezić | LSST: Comprehensive NEO Detection[END_REF]] as the first step or tracklet and track identification followed by the techniques introduced in this paper as the second step. We have mentioned the possibility of using the algorithms of [Milani et al. 2005b] as the third step. Even more complicated pipelines can be conceived but the discussion of pipeline design is beyond the scope of this paper and will be the subject of future work.

In a series of three papers ( [Milani et al. 2005b], [START_REF] Ivezić | LSST: Comprehensive NEO Detection[END_REF]] and the present one) we have defined a set of algorithms that may be used to process astrometric data for Solar System objects when the sky-plane density is much larger than it is for contemporary surveys. This will very soon be the case with Pan-STARRS [Hodapp et al.(2004)] and LSST [START_REF] Ivezić | LSST: Comprehensive NEO Detection[END_REF].

Our work on algorithm definition is a necessary step to exploit their superior survey performance and provide orbits for most observed objects. [ [START_REF] Ivezić | LSST: Comprehensive NEO Detection[END_REF]] Kubica, J., Denneau, L., Grav, T., Heasley, J., Jedicke, R., Masiero, J., Milani, A., Moore, A., Tholen, D., Wainscoat, R. J., 2007.

Efficient intra-and inter-night linking of asteroid detections using kd-trees.

Icarus 189, 151-168.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

Table 6: Accuracy results. Performance of our alogrithms before (columns 2-4) and after (columns 5-7) normalization. For each case we provide the total number of false identifications that passed our quality checks, the percentage of identifications that are false (with respect to the total number of identifications), and the number of identifications containing false tracklets.

Region

All