Gaia Focused Product Release: Asteroid orbital solution. Properties and assessment

Context. We report the exploitation of a sample of epoch astrometry for 157 000 asteroids, the same object in the Gaia Data Release 3, extended over the time coverage planned for the Gaia DR4, which is not expected before the end of 2025. This data set covers more than one full orbital period for the vast majority of these asteroids. The orbital solutions are derived from the Gaia data alone over a relatively short arc compared to the observation history of many of these asteroids. Aims. The work aims to produce orbital elements for a large set of asteroids based on 66 months of accurate astrometry provided by Gaia and to assess the accuracy of these orbital solutions with a comparison to the best available orbits derived from independent observations. A second validation is performed with accurate occultation timings. Methods. We processed the raw astrometric measurements of Gaia to obtain astrometric positions of moving objects with 1D sub-mas accuracy at the bright end. For each asteroid that we matched to the data, an orbit fitting was attempted in the form of the best fit of the initial conditions at the median epoch. Results. Orbits are provided in the form of state vectors in the International Celestial Reference Frame for 156 764 asteroids, including near-Earth objects, main-belt asteroids, and Trojans. For the asteroids with the best observations, the (formal) relative uncertainty is better than 1E10. Results are compared to orbits available from the Jet Propulsion Laboratory and MPC. Their orbits are based on much longer data arcs, but from positions of lower quality. The relative differences in semi-major axes have a mean of 5E10 and a scatter of 5E9.


Introduction
In the course of its systematic scan of the sky, the ESA Gaia spacecraft has detected a wide variety of celestial sources.It will continue operating until the end of its operational life sometime in early 2025.The celestial sources are primarily stars from the Milky Way, which form the bulk of the massive Gaia catalogues that were released between 2016 and 2022 in data releases Gaia DR1, Gaia DR2, Gaia E-DR3, and Gaia DR3 proper (Gaia Collaboration 2016b,a, 2023).However, sources closer to Earth are also caught in the Gaia net, with representatives of every category of small Solar System bodies (SSSBs), such as eear-Earth objects (NEOs), main-belt asteroids (MBA), Trojans, and a few more distant objects from the Kuiper belt and from the trans-Neptunian region.In addition, Gaia data include the largest planetary moons and a sample of periodic comets that were observed not too far from their perihelion passage.
In this Gaia-coordinated Focus Product Release (Gaia FPR), the astrometry of 156 792 minor planets and 31 natural satellites is provided.However, this paper focuses only on the orbit computation of minor planets, and only this category of Solar System bodies is investigated.Likewise, photometric and spectral data are also beyond the scope of this work and have been largely presented with the Gaia DR3 results in Tanga et al. (2023).
Gaia DR2 has shown the quality of Gaia observations and their huge value for asteroids.Based on the high-precision astrometry and because the time span covered by the observations was long enough (at least ≃ 1000 days), it was possible to derive accurate orbits for the observed asteroids, as was shown in Gaia DR3 (Tanga et al. 2023), who reported a summary of the observational technique and the epoch astrometric accuracy.More precisely, Gaia DR3 was built upon 34 months of Gaia operations extending from July 2014 to May 2017.A A37, page 2 of 24 well-observed Solar System body came into the Gaia fields of view about 25-30 times on average, resulting in up to nine accurate astrometric positions at every passage 1 .Tanga et al. (2023) ranked the orbits from excellent to very approximate, with a clear difference in reliability as a function of the fraction of the orbital period covered over 34 months.
For the majority of the MBAs, this time was too short to remove strong constraints between orbital elements.In the processing reported in this paper, which yielded the Focused Product Release, the data arc can be as long as 66 months.This typically is a little longer than the orbital period of an average asteroid that orbits in the main belt between Mars and Jupiter.Therefore, the situation considerably improved as soon as we had a few pairs of data points that lie one orbital period apart from each other.This is sufficient to remove the partial degeneracy that is observed when the arc was too short.Although we have almost twice as many observations as in Gaia DR3, the main source of improvement remains the complete coverage of the orbit in the data set, which is in sharp contrast with the previous solution in Gaia DR3, which was based on only 34 months of data.The length of arc that is used is the primary factor in the improved uncertainty of the orbit.
In this paper, we aim to present and discuss the data that are released in the Gaia FPR, together with the results of the orbital fitting obtained in an extended time span, using a refined version of the fitting algorithm, with better weighting and an improved outlier rejection.Several shortcomings in the code were found and corrected in this new exploitation, and the orbits we provide now are the best we can achieve at the moment with the extremely high astrometric precision of Gaia.It is important to mention at this stage that the orbital solutions are derived from the Gaia data alone, meaning that the arc is still relatively short in comparison with the fits performed at MPC or JPL using all the available observations that so conveniently archived by the MPC, which were collected over decades or even centuries.As shown by the results, the relatively short time-span is somewhat balanced by the homogeneity of the data and the unparalleled single-observation accuracy and precision.
As mentioned above, no new photometric data in G band or reflectance spectra are released in the FPR.To associate photometry with FPR data, the G fluxes and magnitudes published in DR3 can be adopted for the corresponding transits in common.New spectrophotometric processing will be available in DR4, which is currently scheduled for the last term of 2025.Likewise, new astrometry is published in the FPR for planetary satellites, but their orbits have not yet been adjusted on these new positions.
The paper is broken down into a few sections.In Sect.2, we discuss the observation material used in the paper, and then we very briefly recall the relevant properties of the astrometric solution in Sect.3. In Sect.4, we present the principle of the orbit fitting that produced the set of state vectors at a median epoch for each asteroid.The validation and discussion of the results are taken up in detail in Sect.5, along with some recommendation for users. 1 We adopted the following terminology for the data sequencing.A passage or a transit is the crossing of a Gaia field of view.During a passage, up to nine position measurements at the CCD level are possible (primarily in the along-scan direction), which are referred to as observations, or elementary observations to remove ambiguity.An epoch or a visibility period is a set of a few successive transits within one or two days, followed by a gap of several weeks before the start of a new period.Blue shows the first 100 000 numbered asteroids with a median of 330 individual observations (91 785 orbits), and red shows the numbered asteroids with numbers ≥100 000 and a median of 190 observations (64 979 orbits).

Properties of the data set
Gaia scans the sky throughout the year in a continuous way, without virtually any significant interruption beyond unavoidable but short dead times.The principles have been described in several papers from the Gaia Collaboration, and we refer to the reliable information therein (Gaia Collaboration 2016a).The relevant facts for the orbital solutions we use here are the repeated observations of the same source, with an interval of 6-10 weeks between two sequences.A typical sequence consists of a small set of successive passages in the Preceding Field of View (pfov) and Following Field of View (ffov).The most common sequence is just a simple sequence of two passages pfov to ffov (PF) or ffov to pfov (FP) within 6 h, and then a few weeks before a new similar sequence.However, every asteroid during the lifetime of Gaia enjoyed one or two much longer sequences, such as PFPF... or FPFP... and so on.In some cases, there were more than 15 consecutive transits within 4 or 5 days in a row.For the well-observed asteroids, typically, 40 to 90 passages (field crossings) were made in 66 months, and the number of visibility periods (closely packed repeating sequences at distinct epochs) is 30 ± 10 on average.This has the strongest effect for the orbit determination because it drives the sampling in orbital phases.
During each passage, a maximum of nine astrometric positions is obtained.These make up the observations at the CCD level.The number of observations that were used (after rejections in the iterative model fitting) is shown in Fig. 1 for the asteroids that were observed best with numbers <100 000, and for the fainter and less frequently observed asteroids with numbers ≥100 000, typically with half as many transits ending with a successful Gaia detection.The transits were missed because the magnitude was too faint in some part of the orbits or because problems were found in the astrometric solution, which led to a rejection.The transition at 100 000 between the two groups is just a convenient number to obtain balanced distributions.It has no deeper significance.
Beyond the number of transits, a key parameter for the orbital solution is the time coverage, or in this context, more precisely, the fraction (<1 or >1) of the orbit that is covered between the first and last observation.That is to say, the arc length measured in orbital periods.This is shown in blue in Fig. 2 for the asteroids that were observed best with numbers <100 000, and the fainter and smaller asteroids with numbers ≥100 000 are superimposed in red.The typical coverage is about 1.2 orbital periods for the first group, and it is about 0.5 orbital periods for the Trojans.The typical coverage is shorter at about one orbital period for the fainter asteroids, with a marked tail of short arcs (less than one orbital period) that is caused not by the periods themselves, but by a systematic of the Gaia sensitivity threshold: the faint asteroids are not observed at all elongations that are theoretically allowed by the Gaia scan, but only when their distance to the Earth is small enough to have an apparent magnitude <20.7 in the Gaia white band.This on-board brightness cut is imposed by the instrument sensitivity and the ground-link data rate.An orbit is better determined when the observation arc covers a complete orbit or more, and the solution is less constrained for short arcs.This qualitative insight is confirmed in Sect.5, in which we analyse the solutions.

Astrometric solution
The astrometric processing of Solar System objects has been presented in Tanga et al. (2023, Sect. 3.4) and was described with more technical details in the Gaia DR3 documentation in Muinonen et al. (2022, Sect. 8.3.3).We refer to these papers for more detailed information.Here we just recall the key points needed for the readability of this paper and draw attention to the significant departures from the classical ground-based CCD or traditional photographic plate astrometry.As a mission destined to survey the sky several times, a scanning strategy was imposed early during the mission design.The Gaia time sampling is a direct consequence of the parameters selected to optimise this intricate scanning law.Because of the Gaia core science, the free parameters were not optimised with Solar System sources as the main targets, but for the best science return in stellar and fundamental physics.The Solar System objects are well sampled, however, just below the sky average, because their ecliptic latitude is predominantly low.This sky area is less frequently visited by Gaia than the mid latitudes.Each passage of a source in the field of view nominally includes nine crossings of the astrometric CCDs either in the preceding or the following field of view.This provides as many 1D astrometric locations in the instrument frame (pixel number) that are perfectly time tagged.These repeated measurements form the basic astrometric data from Gaia.Based on the instrument calibration and the orientation parameters (attitude), these instrument-tied positions can be transformed into nine astrometric positions on the sky, with an excellent accuracy (≈ sub-milliarcsecond to milliarcsecond level) in the scan direction and much poorer accuracy (up to 0.5 ′′ ) in the transverse direction.When it is transformed into right ascension and declination, the uncertainty appears degraded in both coordinates, but the true statistical information is kept in the covariance matrix, where correlations are close to 1.This encodes the extreme accuracy in the along-scan direction.Therefore, a standard usage of the RA, Dec positions that does not take the correlation coefficient into account would completely fail to bring the true Gaia accuracy.
The positions used in this paper to fit the orbits are at the CCD level, meaning that the number of observations is always the number of validated CCD crossings.Each field transit has a maximum of nine crossings, but the actual average is between seven and eight crossings, so that there are about seven to eight times as many observations as field transits.These are also called observations in other contexts of the Gaia data processing.Consistency checks between the multiple 1D positions in the same transit were used to accept or reject observations for inclusion in the data release, taking the expected linear motion of the Solar System object during the ≈ 40 s of the field crossing into account.All these published observations are considered in the subsequent fits, and normal points for position and velocity were not computed.
We implemented an updated version for the FPR reported in this paper for the astrometric calibration used in DR3.The global transformation that produces astrometric positions is also more accurate because it is now based on the full duration of the nominal Gaia mission of 66 months.Moreover, the FPR processing is mostly independent of the previous data releases.Even in the overlapping time between DR3 and the FPR, observations can be found in one data set but not in the other.This is a result of a different selection that was used by the internal validation procedures or residuals that are now below or above the acceptance threshold with the new orbit.These are rare occurrences, and they arise as a natural consequence of the reprocessing of the whole data set.
The precision of the astrometric measurement in the alongscan direction has been much discussed in Tanga et al. (2023).The precision that is relevant here for the orbit fitting is the combination of a random and systematic part.Its variation with G magnitude is shown in Fig. 3, where it is reproduced in the same way as in Sect.3.4.2 of the Gaia DR3 paper The elementary measurement at the CCD level is better than 1 mas for G < 18 mag and grows to reach about 10 mas at G ≃ 20 mag.This makes Gaia data outstanding and unusual in comparison with the core of the historical observations: they are extremely accurate, rather dense, homogeneous, but over a limited time range that is not too short to preclude orbit fitting.
The orbit can be fit either with observation equations expressed in the along-scan and across-scan directions, or in A37, page 4 of 24 the more usual way in RA, Dec coordinates.The use of either coordinate set is rigorously equivalent theoretically in any case provided that the 2 × 2 variance or covariance matrix is well implemented in the global weight matrix.Tests have shown that identical results are obtained with both approaches, as expected, but great care must be exercised in using the statistical information in the covariance matrix in the right manner.

Presentation
We only consider a subset of known asteroids here that have known (elliptical) orbits, as in Gaia DR3.We determined a leastsquares orbit solution for them that represents the best fit to Gaia data.Detection, follow-up, and initial orbit determination of newly discovered asteroids are made within another pipeline that was described in Carry et al. (2021).
The overall principle of the orbit computation from the Gaia astrometry was summarised in Tanga et al. (2023).As a rule, the general procedure used in the Gaia DR3 release is also used here, but the code has been updated to correct minor bugs and improve the overall consistency between the fundamental constants in their form in astronomical or metric units.

Units of time and length
Gaia time tagging is the Barycentric coordinate time (TCB) internally for the data collection and processing, and externally in the public releases.This means that, for example, the Solar System ephemerides have TCB as an independent variable and that all epoch data are given with a timing in TCB.For stellar astrometry, and even for epoch photometry, the difference between TCB and TDB (or TT), which is about 20 s, is not an issue for the user.This is no longer true for Solar System objects moving at speed of several 10 mas per second, however.This has consequences for the orbit fitting.Observations are timed in TCB, and the equations of motions must be consistent with this timescale.More specifically, the solar mass parameter must be TCB compatible.
Unfortunately, a small glitch has remained in the software between the TCB and our using of an outdated Solar System GM ⊙ equal to the square of the Gauss constant k = 0.017 202 098 95.This is equivalent to our using a unit of length that is not strictly equal to the IAU astronomical unit.Because of the constraints of the complex integrated processing of the Gaia data and the huge validation pipeline, it was not possible to correct for the scale factor in the FPR results.It seems to be a benign correction, and it is, but nothing is simple in a complex system.Fundamentally, this is a mistake, and this issue has been fixed in the software for the planned Gaia DR4.The positive side of this unfortunate bug is that noticing a systematic effect as small as 0.5 × 10 −8 in the results already means something good about the quality of the solution.
Mixing the Gauss constant and the TCB timescale results in an au FPR = 149 597 871 473.216 m instead of the IAU fundamental constant for the scale length in SI equal to 149 597 870 700.0 m.The consequence is that the Gaia published state vectors of the asteroids must be scaled to scale them properly in TCB and au, as it should have been in the first place.Afterwards, the transformation to osculating elements provides a TCB-compatible semi-major axis.More explicitly, we have ρ = au FPR /au = 149597871473.216/149597870700= 1.0000000051686297, which with k 2 ≈ (GM ⊙ ) TDB , is essentially . Now the TCB compatible values in au and au day −1 can be computed from the published components X FPR , V FPR of the state vector with It can also be added that when the FPR state vector in osculating elements is transformed before the scaling, only a single scale transformation has to be applied to the semi-major axis, We show in Sect.5.3.1 how the Gaia TCB compatible orbits and the equivalent from JPL TDB compatible orbits are handled.

Principles and algorithm
The Gaia-centric astrometric data in right ascension and declination at CCD-level are adjusted using a standard weighted linear least-squares (LLS) fit.The error model with non-diagonal covariance matrix as derived in Tanga et al. (2023) provides the weights to be applied for each observation equation.The same set of asteroids as in Gaia DR3 is considered for this new release, but over an observational time span that is roughly twice as long.
For the LLS method, the theoretical expectation is required.To spell out this function, the equations of motion are integrated simultaneously with the variational equations.In the ICRF reference frame and on the TCB timescale, the heliocentric equation of motion is expressed as where r i = |r i | and r i is the heliocentric position vector of asteroid i, GM ⊙ is the mass parameter of the Sun, and Gm i is the mass parameter of asteroid i, which can generally be neglected.The additional perturbing accelerations, added linearly, correspond to They account as in the Gaia DR3 orbital fitting for the planetary perturbations and the relativistic terms, respectively.The part due to the gravitational perturbations (by a total of n p planets, and potentially other perturbing asteroids) is given by The relativistic contribution to the equations of motion can be approximated by The variational equations are generated by the derivation of Eq. ( 1) with respect to the parameter to be corrected.Beutler (2005) and Pontriaguine (1969) give a very thorough exposition A37, page 5 of 24 of how the variational equations can be obtained from a general point of view.
The heliocentric positions of the asteroids are computed from numerical integration2 of Eq. ( 1).In Gaia DR3, the Solar System model is the planetary solution INPOP10e3 .In the current solution, the Solar System model has been updated to INPOP19a, which presents some important differences with INPOP10e.As with INPOP10e, the Earth-Moon barycentre with the seven other planets, dwarf planet Pluto, and a selection of 343 asteroids are included in the Solar System model (Fienga et al. 2013;Deram et al. 2022).For INPOP19a, the ten most massive trans-neptunian objects (TNO) have additionally been included as perturbing bodies.Furthermore, in INPOP19a, influences by the less massive TNOs are modelled by a ring.These added perturbations in INPOP19a cause a shift in the Solar System barycentre that affects the barycentric position of the Gaia satellite (Fienga et al. 2019).
No perturbation of asteroids is included yet in Eq. ( 3).This is still a limitation compared to other databases of asteroid orbits, which consider systematic perturbations by several massive asteroids (e.g.Ceres, Pallas, or Vesta) together with the major planets, or for consistency with INPOP19a.This approximation on the dynamical model introduces a small bias on the orbital solution (quantified in internal validations; see Sect.5.3.1),but almost no effect on the uncertainty of the orbits themselves.Thus we will compare and analyse the formal uncertainties of the orbital solution, not its accuracy.These systematic perturbations will be included in subsequent data releases (starting from DR4), together with other mutual perturbations, that should also allow the mass determination from a global inversion.
Because mutual interactions between asteroids are neglected, Eqs. ( 1) and ( 3) can be integrated simultaneously with the variational equations for each asteroid separately.The small corrections relative to the reference orbit are then computed at a given reference epoch by an LLS method that minimises the residuals on the astrometric position.The correction is determined for the six-dimension initial state-vector in Cartesian coordinates taken from an auxiliary file of orbital elements.In our case, it was derived from the astorb database 4 and was propagated to the reference epoch.The reference epoch for the orbital fit was taken to be the mid-point of the total time span of the observations available for the given asteroid as a proxy for the weighted mean time of the observational arc.This reference epoch is thus not common to all asteroids in the orbital data base.
To be more concise, we must solve the following linear system of equations for each asteroid source: where the vector dλ(t) = (O − C) represents the value of the difference between the observed and computed values of the measured quantity, in this case, the right ascension and declination for each of the N observations for a given asteroid.dq is the differential correction to the parameters, here the components of the initial state-vector for the asteroid under consideration, and A is the matrix formed with the partial derivatives with respect to these parameters, projected to right ascension and declination.Furthermore, all equations are weighted by the matrix √ W obtained through the corresponding variance or covariance matrix for the right ascension and declination of the observations.
√ W is computed as described in Levinger (1980).We accepted only real values and used only the positive square roots where these are needed in the formulae.Moreover, because the variance or covariance matrix for right ascension and declination of the observations is positive semi-definite and symmetric, √ W verifies We refer to Tanga et al. (2023) for a detailed explanation of the uncertainties involved in the data reduction.As explained therein, each observation of an asteroid at a given epoch, that is, a right ascension α and declination δ, has an associated variance or co-variance matrix of the form with α * = α cos (δ) in differential quantities.Then the weights are obtained as W = γ −1 k for the kth asteroid.We recall that the off-diagonal elements play a prominent role in the uncertainty deduced for the sky coordinates, as already stressed in Sect.3.
Planetary aberration as well as the gravitational light bending are then added to the computed coordinates as so-called corrections to the observations to construct the (O−C) vector in Eq. ( 5).All corrections are computed in the general relativistic framework of Gaia data reduction, considering the source at finite distance (Klioner 2003).While the relativistic light deflection effect, which can reach several milliarcseconds at moderate solar elongations, is not included in ground-based observations, it is mandatory in the treatment of Gaia data due to the submilliarcsecond astrometry in the position of the asteroids.We recall that the published right ascension and declination are corrected for the annual aberration, that is, the equivalent to the stellar aberration due only to the barycentric motion of the observer, without any contribution from the asteroid motion.
The LLS is based on the optimisation of a target function, which in our case is defined by that is, the sum of the squares of the weighted residuals (Milani & Gronchi 2010) of the superscript T meaning the transpose of the vector (O − C) j .As in the Gaia DR3 release, optimisation is performed for each asteroid individually, considering all astrometric observations N.There are no more general parameters that would couple observations from different planets in the fit.Thus, the LLS solution of the system in d q = (A T WA) −1 A T Wdλ with corresponding variance/covariance matrix (A T WA) −1 where A T indicates the transpose of the matrix A.

Iterations and convergence
The procedure of orbit correction is iterated until convergence is reached with a given tolerance p following Milani & Gronchi (2010, Sect. 5.2); initially this parameter has a value of 10 −8 however, it can be increased in some cases.A preliminary A37, page 6 of 24 outlier filtering is performed during the astrometric reduction (Sect.3), further outlier rejection at the observation level has been implemented and follows Carpino et al. (2003).
The procedure is built around a post-fit χ 2 computed as χ 2 = ξ j γ −1 ξ j ξ T j ; with ξ the residuals of all observations evaluated with respect to the corresponding orbit and γ ξ j their expected covariance matrix.The threshold for ξ is a parameter of the fitting process.The tolerance test for stopping the computation is evaluated between iteration k and k + 1 as follows: where p is a tolerance threshold, and Q k is given by Eq. ( 7) for the kth iteration.
If the next iteration is unlikely to improve the corrections, the processing can also be stopped; this can be assessed by considering the size of the last correction using an appropriate norm (see Milani & Gronchi 2010), for example, where M = A T W A is the normal matrix, ∆X is the correction vector, and N is the number of observations.In practice, this criterion is used far less frequently.In addition to the two criteria mentioned above, an absolute maximum number of iterations is also denoted n max .
The iteration starts with p = 10 −8 , and the maximum number of iterations is set to n max = 15.The value of p = 10 −8 is somewhat arbitrary, and different values could be imposed (or even a different tolerance for each condition c 1 and c 2 ).However, according to our tests, it is unlikely that this would significantly affect the overall statistical outcome of the orbit adjustment because typically, only a few iterations are required.When neither c 1 nor c 2 is lower than p after 15 iterations, the iteration is continued and p is replaced by 10 p.When another 15 iterations do not return a c 1 or c 2 that is lower than 10 p, this is repeated one more time for a last trial.In general, convergence is reached after a few iterations, and in practice, the value of p had to be changed very rarely.
Three situations are considered as a failure in the attempt to fit an orbit.First when the maximum number of iterations is reached, second if the orbital solution is not an ellipse and is rejected and third when all observations have finally been rejected by the rejection algorithm.
A total of 59 sources were rejected, but their astrometry still appears in the FPR data; these are identified with a NULL state vector in the tables.Interestingly, 31 of these 59 sources are moons of Mars, Jupiter, Saturn, Uranus, or Neptune.An attempt at orbit improvement was not even performed for these objects.Twenty-four of the remaining 28 asteroids are main-belt asteroids, three are near-Earth objects ((433) Eros, 2001 XV10, and 2003 YV37), and the last source is a Mars crosser (2006 PG1).For two of the main-belt objects, 2005 UB160 and 2001 QG115, the high rejection rate that is probably due to a larger fraction of outliers makes convergence unstable or the orbit computation based on overly short arcs unreliable.No convergence was reached.A dedicated investigation of this subset of asteroids, which consists of only 0.2‰ of the orbital catalogue, is needed to fully understand why convergence has not been reached.This is beyond our concerns in this FPR, however, because (433) Eros is a difficult case in Gaia DR3 and requires further investigations.

Post-fit residuals
The full processing resulted in 156 762 orbits for the FPR.In contrast to Gaia DR3, there is no post-processing of the FPR orbit determination.In the previous catalogue, a further filtering of solutions was applied based on the relative uncertainty in the semi-major axis σ a /a and arc length.We refer to Tanga et al. (2023, Sect. 3.5) for an explanation.
The time span covered by the FPR observations is roughly twice that of Gaia DR3.Although here the astrometric processing is more accurate (see Sect. 3), it was a reasonable guess to expect that the arc covered by any object would be at least as long as it was in the Gaia DR3 release.This is true for an overwhelming majority of the sources (99%) and even 76% having the ratio R of the total arc length of the FPR sample to that of the Gaia DR3 at least equal to 2. Figure 4 shows the normalised distribution of R where a sharp peak is seen around R = 2 followed by a steep decline.A substantial fraction of the sources have a ratio R > 2, while the global time coverage from Gaia has changed from 34 to 66 months.These happen to be sources with actual time coverage, between first and last observation, much below 34 months in Gaia DR3, thus increasing automatically the ratio in FPR.One can see also that there are very few sources with a ratio less than one, i.e. an unexpected longer arc in the Gaia DR3 data than in the FPR sample.This occurs when the orbit adjustment rejects a large part of the observations, resulting into a smaller orbital arc; only eight sources (over approx.157 000) belong to this group.Figure 5 displays the along scan residuals against the G magnitude, G mag .The 5σ cut-off value for rejection is visible; it is the envelope of the accepted observations, i.e. the limit between the grey and coloured points.Figure 6 shows the corresponding histograms for normalised residuals ∆AL/σ AL for sources with a G mag less than or equal to 11.5, or greater than 11.5, respectively.For accepted observations with G mag greater than 11.5, the normalised distribution is Gaussian with a mean of 0.018 and standard deviation of 1.08; for these sources the orbit adjustment is acceptable.Sources with G mag less than 11.5 do not show a Gaussian distribution, even if the rejected solutions are included; in this magnitude range the model is not entirely satisfactory, but can be improved by including some other effects, as the photocentre offset for example.For G mag > 11.5 the symmetry is striking.
Figure 7 displays the residuals in along and across scan.There is a small offset from zero, (∆ AL, ∆ AC) ≈ (0, −15) mas A37, page 7 of 24 Fig. 5. Along-scan residuals against G magnitudes.The grey points represent all observations.The coloured points are the observations that were accepted by the orbit adjustment procedure after outlier rejection chosen at |O − C| > 5 σ.For G mag < 11.5, the distribution is not Gaussian.This is clearer in Fig. 6.For G mag ≥ 11.5, the distribution is Gaussian.The orange lines in the plot are ±1 σ, the blue lines are ±2 σ, and the yellow lines are ±3 σ.The number density is given by the scale on the right.A representation of the astrometry quality over a single transit can finally be obtained by considering the dispersion of the measurements that were recorded during each passage of a source in the focal plane in the along-scan direction (Fig. 8).This figure can be compared to Fig. 12   for Gaia DR3.The figures appear nearly equivalent, with a slight improvement around the G-magnitude interval of the best performance (G mag ∼ 12-13).

Orbital elements
The output of the orbit fitting for a particular asteroid consists of a 6D vector of the initial conditions at a particular A37, page 8 of 24 reference epoch, the midpoint of the observation times, and is given in TCB.This vector combines a position and a velocity vector and is heliocentric.It is given on the ICRF axes.The transformation from the state vector to the osculating elements is in principle routine work in celestial mechanics, and fully tested software is widely available.However, here we work with sub-milliarcsecond astrometry, or with orbits at the 10 −9 -10 −10 accuracy level ( as σ(a)/a), and care must be taken everywhere to ensure that no loss of accuracy is associated with the transformation and that the underlying conventions are well understood and implemented.This is the object of this section, to explain all the conventions we used in our solution and describe where they may differ from other usages.
The transformation from the state vectors to the osculating elliptical elements is done in two steps: first, the transformation from heliocentric ICRF frame heliocentric ecliptic frame for the Cartesian coordinates and then the transformation in the ecliptic frame from state vector to osculating elements.In addition, the TCB is used consistently for Gaia while all other sources of osculating elements used TDB in their publications and this impacts the scale length.
The transformation of the state vector to the ecliptic frame depends on just one rotation matrix to connect the two frames.The main part is the rotation ϵ around the x-axis, where ϵ is the obliquity of the ecliptic.The devil lies in the details, as so often is the case.The intersection of the dynamical ecliptic with the ICRF fundamental plane has no reason to coincide with the ICRF x-axis, and this is not the case, either.In addition, the obliquity of the ecliptic relates the ecliptic (taken as a dynamical plane with a strict definition) to the celestial equator.The latter is again different from the ICRF reference plane, and the angle between the two planes of interest is not exactly the physical obliquity of the ecliptic.It is very similar, but not equal.As a consequence, and for lack of a well-agreed convention, several definitions are in use.From a single state vector, we can derive six sets of osculating elements that differ by several tens of milliarcseconds in the longitude of the ascending node or the orbital inclination.This motivated our choice to publish Gaia orbits in the form of a heliocentric state vector in the ICRF axes instead of the usual osculating elements.The former is unambiguous, and from this vector, the orbital elements can be computed so that they agree with the different choices for the ecliptic.They differ slightly with the choice of the ecliptic, although they represent the same orbit.
The issue extends beyond the published osculating elements to the users of this information.They may start from osculating elements of an asteroid to propagate to another time and compute local coordinates using the true equator of the date.Then they will call a precession or nutation package, relating the J2000 celestial equator to another epoch, or the ICRF plane to the celestial equator by using the frame bias above the precession.Therefore, the reference plane and the origin of longitudes for the osculating element are key ingredients for high-precision computation.In this context, high precision means better than 50 mas, which is the magnitude of the small angles in Table 1 or of their difference between institutes.
We have attempted to collate the various choices made for the reference plane and the origin of the longitude of node in this plane by the most relevant sources of osculating elements.Figure 9 shows the different equatorial planes and ecliptics, together with their associated obliquity.The angle ϕ is the value of γ xx O ICRF (xx: ICRF, Gaia, JPL) measured in the ICRF equator and taken positive, like right ascension from the local γ xx   Table 1.(the sign convention is opposite to Chapront et al. 2002 and similar to Hilton & Hohenkerk 2004).The so-called IERS values are taken from the SOFA software and come from the precession formulae derived in Fukushima (2003).The angle ψ is the angular distance between the intersection of the ecliptic with the ICRF plane and the intersection with the J2000 equator, positive like an ecliptic longitude.The numerical values of the angles needed to make the transformation from the ICRF to the selected ecliptic are given in Table 1.The Gaia offset was taken from Chapront et al. (2002) before the adoption of the Fukushima precession by IERS and SOFA, and has not yet been changed for continuity reasons.The ecliptic plane and origin in the ICRF equator are extremely similar (within 3 mas) to the VLBI-derived ecliptic adopted in SOFA.However, the origin of the longitude of node is not on the J2000 equator, but in the ICRF equator.
Nothing is fundamentally wrong with these multiple realisations of the ecliptic, this is just inconvenient.This all stems from the fact that the ecliptic concept is ambiguous when high accuracy is required because the plane orthogonal to the mean angular momentum of the Earth-Moon system depends on what is included in this mean, over which time it is averaged, whether long-period wiggles from planetary perturbations are removed or not, and so on.A less physical definition can appear as well with a rigid link to the ICRF based on conventional angles, in a similar way as the Galactic coordinates are related to the equatorial frame at J2000, with no attempt to recreate a physical definition accurately.
In principle, an unambiguous definition together with numerical constants have been provided in 2006 by the IAU Working Group on Precession and Ecliptic (Hilton et al. 2006) and A37, page 9 of 24 the ensuing IAU Resolution B1 in 2006.However, a unique realisation is currently not agreed, and even within the Gaia Consortium, the ecliptic has been inherited from HIPPARCOS, and change in standards in a large data-processing system is always weighed against the loss of continuity and is decided only when unavoidable.It is important for the users to be aware of these differences and to be on their guard when using and propagating orbits with a target accuracy lower than 100 mas.The differences between osculating elements may be deceptive and may not reveal real differences in the orbits.When properly transformed into state vectors in the ICRF, they may look much more similar than from their elliptic elements.
Table 1 clearly shows that all the major providers of orbits use the same definition, in which the ecliptic is very simply related to the ICRF equator by a single rotation about the ICRF X-axis, with an angle equal to the obliquity ϵ = 84381.4480′′ .We have adopted this convention in the following sections to compute the orbital elements and compare them to other solutions.However, to let the users select their preferred option, and above all, to avoid incorrect assumptions about the ecliptic that is used for the orbital elements, we decided that primary Gaia results will be published as state vectors in ICRF instead of osculating elements.Nevertheless, the covariance matrix between osculating elements is provided because it is independent (to ≈10 −8 ) of the small differences between the ecliptic frames of Table 1, and its computation from the state vector could be tricky without an easily available standard transformation routine.
Incidentally, the Gaia astrometric data for asteroids and Solar System objects themselves can provide the link between the kinematically non-rotating reference frame (ICRF) and a dynamically non-rotating reference frame (e.g.ECJ200, represented by this inertial ecliptic from the planetary solution, or an invariant plane).

Overall statistics
The same set of asteroids as was selected for Gaia DR3 was used in this processing, but over a time span of 66 months instead of the 34 months for Gaia DR3.The number of asteroids with validated astrometry that were processed this pipeline is 156 825, and we obtained converged solutions for 156 762 asteroids with an accepted state vector, or equivalently, a set of six orbital elements, at a particular reference epoch.The success rate may seem particularly high, and it is, but this results from the fact that the sources selected for this run were precisely those with a successfully computed orbit with the Gaia DR3 data, but with a shorter time range.We described the filters we applied in the course of the processing to accept a solution in Sect. 4.Although in this sample, we have numbered asteroids up to = 400 000, many of the asteroids with large numbers are not observed by Gaia or are too poorly sampled for us to compute an orbit from the small number of available positions.This feature is shown in Fig. 10, with the level of completeness as a function of the rank of the asteroid for bin of 1000 or 10 000 planets.The top plot includes up to number = 50 000 and is virtually 100% complete (99% of the first <25 000 planets with an orbit, and 97.12% of the first 50 000).The bottom plot shows a similar histogram for the full data set, showing the sharp decrease in completeness level for smaller and fainter asteroids.It reaches below 50% (5000 solutions in a bin of 10 000 planets) above number 140 000 and reaches 10% from number 300 000.Not much more improvement is expected in the future releases (Gaia DR4   and Gaia DR5) because the missing asteroids are just not detected and will remain so due to Gaia selection effect shown in Fig. 11.However, the selection boundary at number 400 000 will be abandoned, and asteroids with larger IDs will be included in the match algorithm, but they will surely be detected with a lower completeness level.

Semi-major axis
Of the means contrived to quickly grasp the global quality of a new orbit solution of a large number of asteroids, the semi-major axis remains the most significant and favoured orbital parameter.It has the desirable feature over the angular parameters of being a true geometric quantity that is independent of the reference frame.This is particularly valuable when two solutions are compared for which small differences in inclination or in longitude of node may be unconnected with the orbits, but show up just A37, page 10 of 24 as consequences of tiny rotations between reference frames.The eccentricity arguably shares this geometric nature, but the third Kepler law, linking the orbit size to the period, endows the semimajor axis with a much deeper meaning and has been taken for years as the first-rate parameter to compare solutions or to assess a solution.Above all, however, a difference in semi-major axis between two orbits has its counterpart in the mean motion and will soon show up as a secular drift in longitude, dominating any other difference from the other orbital elements and driving the accuracy of the ephemeris a few years away from the epoch.Therefore, we follow this practice here by emphasising the semi-major axis over the angular quantities, which are less significant and are too sensitive to small differences between reference frames.The plots always take the dimensionless quantity σ a /a to express the relative precision, where σ a is the formal standard deviation of the semi-major axis computed by transforming the covariance matrix of the state vector to the matrix for the osculating elements.Similarly, differences between two solutions are shown with the relative distance ∆a/a or its unsigned value ≈ amplitude).
The two plots in Fig. 12 show the relative formal uncertainty as a function of the semi-major axis for the best solutions (top) and for the whole catalogue (bottom).The overall appearance is as expected based on the observational span of 66 months: the Trojan orbits cannot be retrieved with a quality comparable to that of the bright NEOs even after scaling by the semi-major axis.For the first group, the relative uncertainty ranges from 5 × 10 −11 to 2 × 10 −9 for the MBA, and the density peaks about 2 × 10 −10 .
The actual distribution per bin of accuracy is shown in Fig. 13   the bright population with numbers <20 000 is better visible, as is the lower accuracy for the fainter group with a centre at ≈5 × 10 −10 in σ a /a.The formal uncertainty is directly derived from the fit, and the covariance matrix is scaled on the post-fit residuals.
Figure 14 compares the relative uncertainty between the Gaia DR3 and this FPR solution.It clearly shows the gain resulting from a longer and better-sampled arc for the orbit computation.The uncertainty in Gaia DR3 was typically ≃50 times larger, and the distribution was bimodal, with a second population.The second smaller peak in the histogram stems from sources with the shortest arcs (around typically 100 days).This second peak has now disappeared in the FPR as a direct consequence of a fuller coverage of the trajectories.More precisely, the first peak in the Gaia DR3 computations was at ∼2.6 × 10 −8 and the second at ∼6.9 × 10 −6 , while the distribution peaks at ∼8 × 10 −10 in this study.This means a gain between one to two orders of magnitude in the uncertainty, which is much larger than a ≃1/ √ 2 improvement that would be obtained just by doubling the number of data points used in the adjustment.It essentially results from the longer observational arc, in agreement with Desmars et al. (2013).A final plot for this section shown in Fig. 15 gives the formal uncertainty of the semi-major axis as a function of the arc length measured in orbital periods.The improvement as soon as the observation arc covers a full orbit is clearly visible in A37, page 11 of 24 the steep slope on the left in the diagram, compared to the near stationary or much slower decrease when the arc is ⪆1.2 orbital period.In this regime, the improvement will come mainly from the combination of the larger number of observations and the extended arc length, while in the small-arc regime, the photon noise is not the main source of (in)accuracy.At mission completion, with 10.6 yr of data, even the Trojans will have nearly a full orbit coverage, and their orbit uncertainty will decrease to the 1 × 10 −9 relative uncertainty.For the purpose of comparison, the light blue dots in Fig. 15 are plotted from the Gaia DR3 orbits, when only 34 months of data were used in the orbital solution, and almost no orbit had a complete orbit coverage.The improvement by a factor about 50 is outstanding in this combined plot.
The main-belt orbits are not even as good as the Trojan orbits in the 66-month solution.

Comparisons to other solutions
As mentioned earlier, several other sources of orbital elements are built on different observational material, different time spans, and different dynamical models.Although we compared the solution from Gaia to the Astorb (Moskovitz et al. 2022), MPC, and JPL (Park et al. 2021) solutions, only the comparisons to JPL orbits are discussed in detail in this paper, and MPC is mentioned only in passing.

Comparison to JPL orbits
The JPL service named Small-Body Data Base Query 5 allowed us to obtain orbital elements without a truncation in the output for a sample with a reasonable size.This was done for the first 50 000 numbered asteroids, which is enough to have a good sample covering a wide range of semi-major axes.The elements are provided at the epoch JD 2460000.5 TDB (25 February 2023) with all the digits, including a safety margin beyond the true precision.The epoch is offset by about 5 yr from the Gaia mean epoch in our solution, however.A propagation was therefore needed of one of the files, or for both files midway.We chose to propagate the JPL reference data to the Gaia epoch, that is, backward 5 yr.This propagation itself is the source of an additional difference between the orbits becausee it is hard to certify that an accuracy better than ≈ 1 × 10 −10 is maintained over 5 yr of numerical integration, and above all, that the Gaia and JPL dynamical models are close enough to reach this performance.This propagation for Gaia was made outside 5 See https://ssd.jpl.nasa.gov/tools/sbdb_query.html the main pipeline with an independent code.The dynamical model includes the complete form of the EIH equations and gravitational perturbations from 14 asteroids.
To ascertain that no significant deviation arose during our numerical integration, we tested this propagation on a much smaller set of asteroids with JPL orbits requested at the Gaia epoch.We were able to reproduce JPL integration at a level of 8 × 10 −10 in the differences ∆a/a after 5 yr.This is acceptable for the comparison, although we know that many of the Gaia orbits are better than this limit (at least in formal uncertainty).It is likely that these residual deviations mainly hold for the small differences between the two dynamical models (not for the same perturbing asteroids, the same masses, or the same Solar System ephemeris, etc.) and not from the numerical integration proper.
We now describe how Gaia used TCB for its time tagging, as explained in Sect.4.2.For a meaningful comparison with JPL or MPC, the same units must be used.Starting with the Gaia TCB-compatible state vector, a standard TCB-to-TDB scaling is performed as as discussed in Klioner et al. (2010) or Klioner (2008).Here, L B = 1.550519768 × 10 −8 is a constant defined in the system of astronomical units and is equal to the rate (dTCB/dT DB − 1).The scaling must not be applied again to the semi-major axis that was computed from the rescaled state vector.The effect is already included in the spatial components of the state vector.
When the two sets of osculating elements are brought to a common epoch, the differences in semi-major axes as shown in Fig. 16 are obtained in the form of the dimensionless quantity ∆a/a, with a few summary statistics.They are shown in Fig. 17 with a scatter plot for the first 50 000 numbered planets.The agreement is outstanding, even though a slight bias of 5.5 × 10 −10 remains between the two solutions for the first 20 000 numbered planets.Expressed in angle, when ∆a is taken as a positional uncertainty, which is close to the osculating epoch before the drift in mean anomaly takes over, this is equivalent to 0.1 mas.This is a very small deviation given the total independence in the construction of the two solutions.The differences in the dynamical model, in the masses, and above all, in the set of observations (none in common) drive the fit over two very distinct time spans.For a first full Gaia solution, this exceeds the most optimistic A37, page 12 of 24 expectations.The bias is even smaller at 1.3 × 10 −10 for the first 5000 planets and grows to 8.5 × 10 −10 when the full set of 50 000 is considered.This slight increase in asteroid number is hardly noticeable in the scatter in Fig. 17.The scatter about the mean in this plot is otherwise remarkably uniform for the range of selected planets.Its value at ∼5 × 10 −9 agrees with the formal uncertainty of either catalogue, as shown in Figs.[13][14][15]18.It shows that the true accuracy of the solutions is not much different from the formal accuracy.We are probably very close to the limit of what can be learnt from this comparison, which is as interesting for JPL as it is for Gaia.The normalised differences computed as ∆a/ σ(a) 2 Gaia + σ(a) 2 JPL are shown in Fig. 19 for the same set of asteroids.No change was applied to JPL-reported uncertainties to account for the propagation from the JPL to the Gaia epoch.The bias, albeit normalised, has the same origin as in Fig. 16 and is statistically significant.The scatter at ∼1.45 is not a robust estimate and takes too much weight from the tails.The inter-quartile range (IQR) for a normal distribution (IQR = 1.35σ) gives a more robust estimate of the central width at 1.607/1.35≈ 1.2 and decreases to 1.13 with the first 50 000 numbered planets and even 1.08 with the numbered planets in [25 000, 50 000].This is a very satisfactory result and indicates that the JPL uncertainties are realistic and fully compatible with the scatter of ∆a between the two solutions.Because the formal errors for Gaia are smaller than those in JPL, the normalised plot is much less sensitive to an inflation factor that should be applied to Gaia uncertainties, unless this factor is as large  as 10, in which case, the random uncertainties would be fully incompatible with the post-fit residuals.

Comparison to MPC orbits
The Minor Planet Center at Harvard (MPC) has recently provided access to its orbits with full accuracy, including the uncertainty on each element.A similar comparison was made for the MPC solution for 50 000 orbits, and it confirmed the absence of significant bias (<3 × 10 −9 in this case) and a scatter just twice as large as with JPL solution (≈1.1 × 10 −8 ) for ∆a/a.These differences between the comparison to JPL and MPC might be traced to the fact that JPL used more observations than MPC in their fit, and in particular, they used some of the Gaia DR2 observations.This good complement demonstrates the perfect consistency between three independent fits and establishes that the Gaia solution is currently better than MPC and at least at the level of JPL.Based strictly on formal uncertainties, the Gaia solution appears even better, but we remain cautious about this claim because systematics could be larger and our arc length remains limited.The real test would be to compare the Gaia orbits to accurate observations (e.g.occultations) well outside the time range of the fit.This is a work in progress.The first results presented in Sect.5.6 show that this is a work in itself, and much more work remains to be done for this goal to be reached.
A37, page 13 of 24 ) is just the small statistical right tail of the distribution of Fig. 20 and is determined by the selection threshold for the plot.The group on the right with σ a /a > 1 × 10 −7 comprises solutions with a poor fit performance and is not an issue: these orbits could have been filtered out in the pipeline and were excluded from the delivery.

Outliers from this comparison
The histogram in Fig. 20 has a vertical linear scale that prevents us from seeing the extended tails when the relative differences between Gaia and JPL orbits are greater than 1 × 10 −7 or 1 × 10 −6 .The orbital solution of Gaia has very few filters during the iterative process or in the final output.If the solution has converged, the orbit is output, regardless of the formal uncertainty.Given the possibility of very small effective arcs for some asteroids, which are also associated with the small number of detections, incorrect orbits, if rare, may be obtained in the Gaia solution.These are more likely than in JPL, which was built upon a much longer time range for the first 50 000 numbered asteroids considered in this comparison.Setting the level of anomalous differences at 1 × 10 −7 for the difference |∆a|/a, deviations exceed this threshold occur in only 110 orbits.This is indeed an astonishingly small number.It is close to 0.2% of the set and it is well accounted for in the residual analysis, as we show below.Many of these orbits may have been filtered out from the Gaia solution based on their anomalous formal uncertainty that is shown in Fig. 21.This plot shows two well separated groups.On the left, about 50 orbits with Gaia solutions that are acceptable and have a statistically significant difference with JPL.This can be viewed as the normal extended tail of the distribution, and no conclusion can be drawn without considering all solutions individually.The right part is more interesting and shows a strict relation between the departure between the two orbits (Gaia and JPL) and the increasingly lower quality of the Gaia solutions.These ≈60 orbits could be removed for the published set without damage, but it is very satisfactory to see this clear feature for such a small number of orbits.All the others except for one agree better with the JPL solution than 1.7 × 10 −6 in relative difference in semi-major axis, and they agree usually much better than this limit, as shown in Fig. 20.Essentially, there are no unexplained anomalous orbits in the set we used for the comparison.Gaia and JPL can both be proud of this.

Other orbital elements
The comparison of the other orbital elements (eccentricity and the angular elements) is briefly mentioned now with additional plots.The main point is to draw attention to the details of the reference ecliptic used by either group.Any difference in the definition of the ecliptic (both the plane and the origin of longitude) would result in spurious systematic deviations that are unrelated to a real difference between the orbits.As said earlier, osculating elements of Gaia have been compared to the same ecliptic plane as JPL, with origin of longitude at the node with the ICRF equator as defined in Fig. 9.We show in Fig. 22 histograms of the difference in orbital eccentricity and inclination found between the Gaia and JPL solutions for the first 10 000 numbered asteroids.Angular units are used, with dimensionless ∆e taken as radians and converted into arcseconds (as if a∆e were a position shift and (a∆e)/a = ∆e were an angular shift).The two solutions are not biased, and the scatter about the mean is typically 5 to 10 mas or <5 × 10 −8 , with a robust Gaussian scatter of 40 and 43 mas, respectively).This is larger than the equivalent angular deviation seen with the semimajor axis, but it probably gives a good upper bound for the true quality of the orbits.The prominent tails in the inclination are not Gaussian, however.
The following comparisons in Figs.23-25 show scatter plots of the differences between Gaia and JPL for the three angular elements inclination, longitude of node, and argument of perihelion.The last two are scaled with the inclination and the eccentricity to show the equivalent of a displacement on the sky rather than the difference in the less meaningful coordinates.The three plots are very similar.They have a core within 5 mas that is surrounded by an extended tail that increases with the rank of the asteroid.An additional analysis shows that the variation in inclination has a systematic sine wave of ∼10 mas amplitude with the longitude of node.This is a clear sign of a reference-frame effect.However, this effect is hardly visible for asteroid numbers <10 000 and becomes more visible with larger numbers.A37, page 14 of 24

Comparison to stellar occultations
A final severe test of the quality of orbits published in the FPR is provided by stellar occultations.In principle, positions computed by using FPR orbits at the epoch of an observed occultation of a star present in the Gaia catalogue must match the astrometry derived from the event.
We set up our test for the objects with an orbit in the FPR and that were successfully observed by occultation.The selected data set resulted in 978 asteroids, associated with 5774 astrometric measurements by stellar occultations.Large and bright asteroids with better orbits in the MPC or JPL data bases dominate the sample because their occultations are easier to predict.In addition, because they are larger than the faint ones, their occultations last longer and are easier to catch.We propagated the FPR orbits to the epoch of each observed occultation and computed the corresponding position and its difference (O−C) in the directions of the equatorial frame with respect to the observed astrometry.We limit our study to occultations by main-belt asteroids, which are the largest majority of the sample for occultation astrometry and Gaia orbits.In the following, we discuss the total O−C, that is, the distance on the tangent sky plane between the observed and the computed position.
The global distribution of the O−C values (Fig. 26) reveals a large spread that reaches several hundred milliarcseconds.Most observations lie below ≈50 mas.This corresponds to the order of magnitude of the expected uncertainty for the dominating fraction of occultations with a small number of observers.In this situation, the cross-track astrometric error is dominated by the apparent size of the asteroid (Ferreira et al. 2022).
However, when only occultation events contemporary to the Gaia mission (same figure) are selected, a clear subset of much smaller O−C appears.It peaks at about 10 mas and has a tail of high values that is substantially suppressed.
We can then verify whether an equivalent improvement can be obtained by selecting only occultations that are related to objects with very high-quality orbits.We set a threshold on the maximum allowed σ a in such a way that the same number of occultation events as for the 2014-2020 period was selected (1819).The obtained curve shows that this selection is not sufficient to eliminate the high residuals, and the peak of the distribution does not move appreciably with respect to the overall distribution.
The apparent asteroid size (mentioned above as the main source of uncertainty) and other factors contributing to errors of occultation astrometry were factored into the total error budget following the guidelines illustrated by Herald et al. (2020).A result that is more independent of the details of this error model can be obtained by dividing the O−C values by the uncertainty associated with each position (Fig. 26, bottom panel).The obtained distribution is self-consistent and peaks around unity.The selection of the Gaia period clearly stands out as the best sample.
This evidence highlights a trend with time that is revealed by Fig. 27.Uncertainties in the occultation astrometry in the 1990s may be affected by worse statistics, but later observations are rather homogeneous (Fig. 2 in Ferreira et al. 2022).The minimum around the period of nominal operations of Gaia that was used to obtain the FPR orbits is thus not due to the properties of the errors in asteroid occultations, but is clearly related to the time interval covered by the Gaia data.Interestingly, this result allows us to quantify a degradation in the ephemeris of ∼10 mas yr −1 .This clearly only holds for the asteroid sample that is observed via occultations, which is biased towards bright asteroids.Gaia astrometry is excellent but not fully optimal for these sources.
At this stage, we also investigated this time trend for occultations involving objects with only the best orbits (same selection as in Fig. 26).We found no significant difference.This probably means that the degradation in our sample of orbits is more sensitive to the limited observing arc than to the internal differences in quality.

Conclusion
We have reported a new orbital solution for nearly 160 000 asteroids built on the Gaia astrometric measurements collected over 66 months of operation and showed its main properties in terms of accuracy and reliability.Comparisons with the best existing orbits published by the JPL or the MPC have shown the impressive consistency between these solutions, which arises from different observation sets and time coverage and dynamical modelling.We established the absence of systematic differences larger than 5 × 10 −10 in relative difference in the semi-major axes.
This Gaia FPR orbital solution is the first within the Gaia releases whose time coverage is at least equal to the orbital period of main-belt asteroids.This feature is needed to ensure the accuracy of the orbit solution.These expectations were mostly shown to be correct in this paper, with an impressive overall improvement compared to Gaia DR3 solution by a factor almost 50 (see Fig. 15), whereas the number of observations was just twice as large.
An important issue that arose in the course of this work is that institutes that publish orbital elements should agree on the reference ecliptic (same plane and same origin of longitude), while different realisations are available that differ from each other by several tens milliarcseconds.With the current quality reached at JPL, MPC, NEODyS, and Gaia, differences as small as this matter.Tiny deviations of few 1 × 10 −8 in dimensionless parameters or a few milliarcseconds in the angular elements on otherwise perfectly agreeing orbits are easily created with small differences in reference frames.With sub-milliarcsecond astrometry and orbits approaching 1 × 10 −10 in relative accuracy from JPL or Gaia, the community needs to coordinate, possibly with IAU or IERS, in order to agree about the definition and realisation of the ecliptic, physical or conventional, and about its relation to ICRF and J2000.This goes much beyond the orbits of small bodies because it must be consistent with the precessionnutation model the user will apply to obtain the positions in the equator of date.
The next step will be the future Gaia data releases, starting with Gaia DR4.They will complete the current set in two ways: (i) the photometric and spectrophotometric data, and data for comets and planetary satellites not analysed here; (ii) a much larger set of asteroids, about twice as large in number, because this has been limited in the FPR to the same set as selected for Gaia DR3, solely because the astrometric solution for the remainder was not yet available.
A few years from today, the Gaia mission will finally provide data over approximately 10 yr, which will double the time interval of the FPR data we used here, and which will improve the quality, precision, and accuracy of the orbit determination of small Solar System bodies further.With this trend, Trojan orbits are expected to be as good as those of the main-belt asteroids are today.As of June 2023, about 85% of the expected data are already safely archived on the ground for the upcoming processing.
Finally, in the following example, we compute the time span covered by the observations of each asteroid,

Appendix B: Transformation of the state vector into ecliptic
The state vector resulting from the orbit fitting is a 6D vector giving the initial condition at the epoch (one per body) as a combination of the position vector and the velocity vector into a single state vector.The origin is heliocentric, and the axes are aligned to the ICRF axes.The units are au and au/day, and the scaling of Equation 10 must be applied to obtain the TDB compatible values.By reference to Fig. 9 and Table 1, the coordinates of the same vector in the heliocentric ecliptic frame are given by where R x,y,z (α) denotes the passive rotation of angle α about the axes x, y, or z.In matrix form, this is With the state vector given in the ecliptic frame (more precisely, in one of the possible choices of the ecliptic), the transformation into osculating elements is common to all the cases and is normally available with reliable routines in any group handling orbits to perform the transformation in both directions for elliptical or hyperbolic orbits.However, the solar mass parameter must also be included in in this transformation to ensure that it agrees with the mass parameter that is used in the force model.The user who needs to compute a good precision (e.g.lower than 10 mas) position in the mean equator of the date must be aware that some precession routines may automatically include some of the rotations above in the form of a frame bias to relate the ICRF frame to the J2000 equatorial frame.

Fig. 1 .
Fig. 1.Distribution of the number of observations used in the orbital solution.Blue shows the first 100 000 numbered asteroids with a median of 330 individual observations (91 785 orbits), and red shows the numbered asteroids with numbers ≥100 000 and a median of 190 observations (64 979 orbits).

Fig. 2 .
Fig. 2. Observation coverage expressed in orbital periods.Asteroid numbers <100 000 are plotted in blue, and the complementary set is shown in red.

Fig. 3 .
Fig.3.Error model in the along-scan direction for the SSO astrometry in Gaia DR3 as a function of the G magnitude.The total error is represented as given by the squared sum of the random and the systematic component.The colour represents the data density (yellow or light means a higher density).The thick line and the two thin lines on each side are the quantiles that correspond to the mean and the 1σ level.

Fig. 4 .
Fig.4.Ratio of the arc length of the FPR sources to that of the Gaia DR3 (R).The peak at ≃2 is marked, and the excursion about the mode is moderate.

Fig. 6 .
Fig.6.Histograms of the along-scan residuals normalised to the formal uncertainties for the accepted solutions with G mag ≤ 11.5 shown in blue or G mag > 11.5 shown in green.The histogram in magenta shows all sources (accepted or rejected) with G mag ≤ 11.5.A Gaussian fit to the data is shown as well.We note that for sources with G mag > 11.5, the distribution is very well represented by a Gaussian (µ = 0.018, σ = 1.08), but for G mag ≤ 11.5, the distribution is clearly not Gaussian.

Fig. 7 .
Fig. 7. Post-fit residuals of the across-scan direction vs. the along scan direction.Top: true angular values.Bottom: normalised quantities.The orange lines show the median, and the yellow lines show the 5th (left) and 95th percentile (right).The density is given in log-scale, so that only the core of the plots is populated.

Fig. 8 .
Fig. 8. Standard deviation of the values of the along-scan residuals for each CCD, computed over each transit.The colour corresponds to the point density.The lines represent the smoothed average, and the quantiles at the 1σ level.

Fig. 9 .
Fig.9.Relative positions of the ICRF and J2000 equators and of the ecliptic with the different conventions found for the ecliptic frame.The origin of the ICRF (O ICRF ) is offset from the inertial equinox γ J2000 .The relevant numbers are given in Table1.

Fig. 10 .
Fig. 10.Completeness level of successful orbital solutions as a function of the asteroid number per bin of 1000 (top) and 10 000 (bottom).Top: first 50 000 numbered asteroids, with 1000 successes achievable per bin at most.Bottom: whole set of ≈ 157 000 solutions per group of 10 000.

Fig. 11 .
Fig. 11.Number of elementary astrometric observations as function of the number of the asteroid.The black line shows the average values.Asteroids with large numbers discovered since the year 2000 are generally smaller and fainter.A large fraction cannot be seen by the Gaia telescopes and detectors, or is seen only during some limited and more favourable orbital phases.

Fig. 12 .
Fig.12.Formal uncertainty of the orbits measured by σ a /a, where a stands for the semi-major axis of the osculating orbit, and σ a shows its standard deviation.Top: First 50 000 numbered asteroids.Bottom: complete solution of 157 000 orbits.The vertical line at a = 3.15 au corresponds to a period of 66 months, which is the longest arc in the data.The vertical scale is adapted to each set.

Fig. 13 .
Fig.13.Relative formal uncertainty of the Gaia orbits resulting from the fit to the observations over 66 months, divided between the bright end (numbers <20 000) and the fainter asteroids.The bins add up to 100% for each group.

Fig. 14 .
Fig. 14.Histogram of the relative uncertainty of the semi-major axis, normalised to maximum, for all asteroids in the middle main belt.Gaia DR3 is shown in green, and this work for the FPR is shown in red.

Fig. 15 .
Fig. 15.Relative formal uncertainty of the Gaia orbits as a function of the arc-length coverage expressed in orbital periods.The black dots show this solution, and the light blue dots show the Gaia DR3 solution over 34 months.

Fig. 16 .
Fig. 16.Differences in semi-major axis (as ∆a/a) between the Gaia orbital solution and the JPL orbits transformed to the Gaia epoch for the first 20 000 asteroids.There are 19 822 Gaia solutions, 261 of which are outside the plot boundaries.

Fig. 17 .Fig. 18 .
Fig. 17.Differences in semi-major axis (as ∆a/a) between the Gaia orbital solution and the JPL orbits for the first 50 000 asteroids.The bias increases slowly from 1.3 × 10 −10 to 9.5 × 10 −10 from the left to the right.

Fig. 19 .
Fig. 19.Normalised differences in semi-major axis between the Gaia orbital solution and the JPL orbits transformed to the Gaia epoch for the first 20 000 asteroids.Compared to Fig. 16, 13 bodies lie outside the[−10, 10] interval.

Fig. 20 .
Fig. 20.Absolute values of the differences in semi-major axis between the Gaia orbital solution and the JPL orbits transformed to the Gaia epoch for the first 20 000 asteroids.Figure 10 shows 19 822 Gaia solutions with ID ≤ 20 000.The difference of 19 805 comes from the largest outliers of Fig. 21.

Fig. 21 .
Fig.21.Anomalous differences in semi-major axis between the Gaia orbital solution and the JPL orbits.The plot shows the relative difference as a function of the Gaia formal uncertainty σ a /a for all solutions with |∆a/a| > 1 × 10 −7 .The left group (σ a /a < 5 × 10 −8 ) is just the small statistical right tail of the distribution of Fig.20and is determined by the selection threshold for the plot.The group on the right with σ a /a > 1 × 10 −7 comprises solutions with a poor fit performance and is not an issue: these orbits could have been filtered out in the pipeline and were excluded from the delivery.

Fig. 22 .
Fig.22.Differences in eccentricity (expressed in arcsec) and inclination between the Gaia and JPL solutions for the first 10 000 numbered asteroids.

Fig. 23 .
Fig. 23.Differences in inclination between the Gaia and JPL solutions.The core is bounded with ±10 mas, but there are extended non-Gaussian tails.

Fig. 24 .
Fig. 24.Differences in longitude of node (as ∆Ω sin i in arcsec) between the Gaia and JPL solutions for the first 50 000 numbered minor planets.

Fig. 25 .
Fig. 25.Differences in argument of perihelion (as e∆ω in arcsec) between the Gaia and JPL solutions for the first 50 000 numbered minor planets.

Fig. 26 .
Fig. 26.Histogram of the absolute total O−C (distance on the sky) between the predicted and observed positions derived from stellar occultations.The bottom panel presents the O−C values normalised by the uncertainty of the occultation astrometry for each event.The curve labelled "2014-2020" includes only occultations around those years, while "best orbit" considers only orbits with σ a < 1.6 × 10 −10 au.

Fig. 27 .
Fig. 27.Absolute values of O−C for the occultation data set averaged on yearly bases.The error bars correspond to the standard deviation.

Table 1 .
Relevant angles defined in Fig.9that were used to relate the inertial ecliptic to the ICRF fundamental plane.