© The Authors 2023

1 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, Gaia Collaboration 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 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¹. 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.

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).

2 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 6h, 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.

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

3 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 along-scan 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 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.

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.

4 Orbit fitting

4.1 Presentation

We only consider a subset of known asteroids here that have known (elliptical) orbits, as in Gaia DR3. We determined a least-squares 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.

4.2 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⁻⁸ 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 871473.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 $$\begin{array}{l}\rho ={\text{au}}_{\text{FPR}}\text{/au}\hfill \\ =149597871473.216/149597870700\hfill \\ =1.0000000051686297,\hfill \end{array}$$

which with k² ≈ (GM_⊙)_TDB, is essentially $${\left(\frac{{\left(G{M}_{\odot }\right)}_{\text{TCB}}}{{\left(G{M}_{\odot }\right)}_{\text{TDB}}}\right)}^{1/3}.$$

Now the TCB compatible values in au and au day⁻¹ can be computed from the published components X_FPR, V_FPR of the state vector with $$\begin{array}{c}{X}_{\text{TCB}}=\rho \times X_{\text{FPR}}\\ V_{\text{TCB}}=\rho \times V_{\text{FPR}}{}_{.}\end{array}$$

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, $$a_{\text{TCB}}=\rho \times a_{\text{FPR}}{}_{.}$$

We show in Sect. 5.3.1 how the Gaia TCB compatible orbits and the equivalent from JPL TDB compatible orbits are handled.

4.3 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 $${\ddot{r}}_i=-\left(G{M}_{\odot }+G{m}_i\right)\frac{r_i}{r_i^3}+{\displaystyle \sum _p{f}_p,}$$(1)

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 $${\displaystyle \sum _p{f}_{\text{p}}}=f_{\text{p}\mid \text{Grav}}+f_{\text{p|}relat}+\cdots$$(2)

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 $$f_{\text{p}\mid \text{Grav}}=-{\displaystyle \sum _{j\ne i,j\ge 1}^{n_{\text{p}}}G}{m}_j\left\{\frac{r_i-r_j}{r_{ij}^3}+\frac{r_j}{r_j^3}\right\}.$$(3)

The relativistic contribution to the equations of motion can be approximated by $$f_{p\mid relat}\simeq \frac{G{M}_{\odot }}{c^2{r}_i^3}\left\{\left[\frac{G{M}_{\odot }}{r_i}{\dot{r}}_i^2\right]r_i+4\left(r_i\cdot {\dot{r}}_i\right){\dot{r}}_i\right\}.$$(4)

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 of how the variational equations can be obtained from a general point of view.

The heliocentric positions of the asteroids are computed from numerical integration² of Eq. (1). In Gaia DR3, the Solar System model is the planetary solution INPOP10e³. 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⁴ 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: $$\sqrt{W}\text{\hspace{0.17em}d}\lambda =\sqrt{W}\left(\begin{array}{l}\text{Δ}{\alpha }_0\text{\hspace{0.17em}cos}\left({\delta }_0\right)\hfill \\ \text{Δ}{\delta }_0\hfill \\ \text{Δ}{\alpha }_1\text{\hspace{0.17em}cos}\left({\delta }_1\right)\hfill \\ \text{Δ}{\delta }_1\hfill \\ ⋮\hfill \\ \text{Δ}{\alpha }_{N-1}\text{\hspace{0.17em}cos}\left({\delta }_{N-1}\right)\hfill \\ \text{Δ}{\delta }_{N-1}\hfill \end{array}\right)=\sqrt{W}A.\text{d}q,$$(5)

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 $\sqrt{W}$ obtained through the corresponding variance or covariance matrix for the right ascension and declination of the observations. $\sqrt{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, $\sqrt{W}$ verifies ${\sqrt{W}}^T\sqrt{W}=W.$. 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 $$\gamma =\left(\begin{array}{cc}\mathrm{var}\left({\alpha }^{*}\right)& \mathrm{covar}\left({\alpha }^{*},\delta \right)\\ \mathrm{covar}\left({\alpha }^{*},\delta \right)& \mathrm{var}(\delta )\end{array}\right)=\left(\begin{array}{cc}{\sigma }_{{\alpha }^{*}{\alpha }^{*}}& {\sigma }_{{\alpha }^{*}\delta }\\ {\sigma }_{{\alpha }^{*}\delta }& {\sigma }_{\delta \delta }\end{array}\right),$$(6)

with α^* = α cos(δ) in differential quantities. Then the weights are obtained as $W={\gamma }_k^{-1}$ 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 sub-milliarcsecond 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 $$Q=\frac{1}{N_i}{\displaystyle \sum _{j=0}^{N_i-1}{(O-C)}_j^{\text{T}}}{W}_j{(O-C)}_j,$$(7)

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 $\text{d}\overline{q}={\left(A^{\text{T}}WA\right)}^{-1}{A}^{\text{T}}W\text{d}\lambda$ with corresponding variance/covariance matrix (A^TWA)⁻¹ where A^T indicates the transpose of the matrix A.

4.4 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⁻⁸ however, it can be increased in some cases. A preliminary 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 χ² computed as ${\chi }^2={\xi }_j{\gamma }_{{\xi }_j}^{-1}{\xi }_j^{\text{T}}$; with ξ the residuals of all observations evaluated with respect to the corresponding orbit and ${\gamma }_{{\xi }_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: $$c_1=\frac{\left|Q_{k+1}-Q_k\right|}{Q_{k+1}}<p,$$(8)

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, $$c_2=\Vert \Delta X\Vert =\sqrt{\frac{\Delta X^{\text{T}}M\Delta X}{N}}<p,$$(9)

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⁻⁸ , and the maximum number of iterations is set to n_max = 15. The value of p = 10⁻⁸ is somewhat arbitrary, and different values could be imposed (or even a different tolerance for each condition c₁ and c₂). 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₁ nor c₂ 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₁ or c₂ 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.

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

4.5 Post-fit residuals

The full processing resulted in 156762 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 over-whelming majority of the sources (99%) and even 76% having the ratio ℛ 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 ℛ where a sharp peak is seen around ℛ = 1 followed by a steep decline. A substantial fraction of the sources have a ratio ℛ > 1, 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. 157000) 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 or (∆ AL/σ_AL, ∆ AC/σ_AC) ≈ (0.05, −0.028), this is also the true in Gaia DR3. In the ∆ AC/σ_AC vs ∆ AL/σ_AL plot the distribution is clearly asymmetric about ∆ AC/σ_AC = 0 but appears to be much more symmetric about ∆ AL/σ_AL = 0. Again, the same is seen in the Gaia DR3 sample and is probably due to the greater uncertainties in the across scan direction. For 90% of the sources ∆ AL is in the interval [−8, +8] mas, ∆ AL/σ_AL ϵ [−1.8, +1.8].

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 in Tanga et al. (2023) 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).

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 Gmag < 11.5, the distribution is not Gaussian. This is clearer in Fig. 6. For Gmag ≥ 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.

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

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 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 quan-tiles at the 1σ level.

4.6 Orbital elements

The output of the orbit fitting for a particular asteroid consists of a 6D vector of the initial conditions at a particular 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⁻⁹−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 sttiff: 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 γ_xxO_ICRF (xx: ICRF, Gaia, JPL) measured in the ICRF equator and taken positive, like right ascension from the local γ_xx (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 the ensuing IAU Resolution Bl 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⁻⁸) 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).

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 (OICRF) is offset from the inertial equinox γJ2000. The relevant numbers are given in Table 1.

5 Analysis and orbit validation

5.1 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 156762 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.

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 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.

5.2 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 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 semi-major 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⁻¹¹ to 2 × 10⁻⁹ for the MBA, and the density peaks about 2× 10⁻¹⁰.

The actual distribution per bin of accuracy is shown in Fig. 13 for the same selections. The uncertainties are given in log-scale. The centre of the disttibution at ≈ 1.5 × 10⁻¹⁰ for 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⁻¹⁰ 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⁻⁸ and the second at ~6.9 × 10⁻⁶, while the distribution peaks at ~8 × 10⁻¹⁰ in this study. This means a gain between one to two orders of magnitude in the uncertainty, which is much larger than $\text{a}\simeq 1/\sqrt{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 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⁻⁹ 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.

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 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 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 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.

5.3 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.

5.3.1 Comparison to JPL orbits

The JPL service named Small-Body Data Base Query⁵ 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⁻¹⁰ 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 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⁻¹⁰ 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 $$\begin{array}{l}{a}_{\text{TDB}}=(1-\text{LB})a_{\text{TCB}}\\ X_{\text{TDB}}=(1-\text{LB})X_{\text{TCB}}\text{, }\\ V_{\text{TDB}}=V_{\text{TCB}}\end{array}$$(10)

as discussed in Klioner et al. (2010) or Klioner (2008). Here, L_B = 1.550519768 × 10⁻⁸ is a constant defined in the system of astronomical units and is equal to the rate (dTCB/dTDB − 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⁻¹⁰ remains between the two solutions for the first 20 000 numbered planets. Expressed in angle, when ∆α 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 expectations. The bias is even smaller at 1.3 × 10⁻¹⁰ for the first 5000 planets and grows to 8.5 × 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⁻⁹ agrees with the formal uncertainty of either catalogue, as shown in Figs. 13–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 $\Delta a/\sqrt{\sigma {(a)}_{\text{Gaia }}^2+\sigma {(a)}_{\text{JPL}}^2}$ 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.

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 19822 Gaia solutions, 261 of which are outside the plot boundaries.

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. 18 Relative uncertainty of the reference orbits computed by the JPL before propagation to the mean Gaia epoch. The bins add up to 100% for each group.

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 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.

5.3.2 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⁻⁹ in this case) and a scatter just twice as large as with JPL solution (≈1.1 × 10⁻⁸) 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.

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. 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.

5.4 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⁻⁷ or 1 × 10⁻⁶. 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⁻⁷ for the difference |∆a|/a, deviations that 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⁻⁶ 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.

5.5 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⁻⁸, with a robust Gaussian scatter of 40 and 43 mas, respectively). This is larger than the equivalent angular deviation seen with the semi-major 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.

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 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 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 Differences in argument of perihelion (as e∆ω in arcsec) between the Gaia and JPL solutions for the first 50 000 numbered minor planets.

5.6 Comparison to stellar occupations

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⁻¹. 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.

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 σα < 1.6 × 10−10 au.

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

6 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⁻¹⁰ 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⁻⁸ 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⁻¹⁰ 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 precession-nutation 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.

Acknowledgements

We are grateful to an anonymous referee of their kind remarks and constructive comments on several key-points of the paper. Their review has helped significantly to improve the organisation and clarity of the paper. This work presents results from the European Space Agency (ESA) space mission Gaia. Gaia data are being processed by the Gaia Data Processing and Analysis Consortium (DPAC). Funding for the DPAC is provided by national institutions, in particular the institutions participating in the Gaia MultiLateral Agreement (MLA). The Gaia mission website is https://www.cosmos.esa.int/gaia. The Gaia archive website is https://archives.esac.esa.int/gaia. The full acknowledgements are available in Appendix C.

Appendix A Access FPR asteroid data in the Gaia archive

The data used in this article can be found at the Gaia mission archive of ESA⁶. By navigating to the “(Advanced) ADQL tab”, the tree of the data base can be expanded. Asteroid data are in the focused_product_release branch, in tables named sso_source (orbits) and sso_observation (epoch astrome-try).

The meaning of the fields for the astrometry is extensively explained in the Gaia DR3 documentation and the related article (Tanga et al. 2023). The orbital elements have their equivalent in the gaiadr3.sso_orbits table, with all the relevant details provided in this article. The residuals used to derive data statistics can be computed starting from the published orbits.

All the Solar System objects present in the FPR were also part of Gaia DR3. However, as the FPR processing is completely independent of the previous processing, the inverse is not strictly true and a small number of Gaia DR3 sources may be missing in the FPR. These differences are very marginal, however.

We provide below an example of an ADQL query that can be used to retrieve the state vector and semi-major axis uncertainty for all the valid orbital solutions, also excluding objects for which the heliocentric orbit cannot be computed (natural planetary satellites):

SELECT number_mp, denomination,
osc_epoch,
ESDC_ARRAY_ELEMENT(h_state_vector,1) AS x,
ESDC_ARRAY_ELEMENT(h_state_vector,2) AS y,
ESDC_ARRAY_ELEMENT(h_state_vector,3) AS z,
ESDC_ARRAY_ELEMENT(h_state_vector,4) AS vx,
ESDC_ARRAY_ELEMENT(h_state_vector,5) AS vy,
ESDC_ARRAY_ELEMENT(h_state_vector,6) AS vz,
SQRT(ESDC_ARRAY_ELEMENT
  (orbital_elements_var_covar_matrix, 1))
  AS sigma_a
FROM gaiafpr.sso_source WHERE ESDC_ARRAY_ELEMENT
(orbital_elements_var_covar_matrix,1)>0.
ORDER BY number_mp ASC

The output is sorted by asteroid number and contains 156 762 orbital solutions. The content of the whole covariance matrix can be added by a similar method. The filter for the first element of the covariance matrix is sometimes negative when the orbit solution is poor.

All the astrometric measurements for a given object can be retrieved based on their identifiers. For instance, for (704) Inter-amnia,

SELECT * FROM gaiafpr.sso_observation
WHERE number_mp=704.

To retrieve the astrometry for all asteroids associated with the best valid orbits with σ_a<r_a<1 × 10⁻¹⁰ au, we can write

SELECT astrom.*
FROM gaiafpr.sso_observation AS astrom
INNER JOIN
gaiafpr.sso_source AS so USING(number_mp)
WHERE ESDC_ARRAY_ELEMENT
(so.orbital_elements_var_covar_matrix,1)>8.
AND SQRT(ESDC_ARRAY_ELEMENT
    (so.orbital_elements_var_covar_matrix,1))<1.e-10
AND astrom.is_rejected = ‘false’

Finally, in the following example, we compute the time span covered by the observations of each asteroid,

SELECT number_mp, MAX(epoch)-MIN(epoch) AS t_delta_days
FROM (SELECT astrom.number_mp, astrom.epoch
FROM gaiafpr.sso_observation
AS astrom INNER JOIN gaiafpr.sso_source
AS so USING(number_mp)) AS subquery
WHERE(number_mp>0) GROUP BY number_mp.

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 $$r_{\text{ecl}}\text{=}{ℛ}_z\left(\psi \right){ℛ}_x(\in ){ℛ}_z\left(-\varphi \right)r_{\text{ICRF}}$$(B.1) $$v_{\text{ecl}}\text{=}{ℛ}_z\left(\psi \right){ℛ}_x(\in ){ℛ}_z\left(-\varphi \right)v_{\text{ICRF}}{}_{,}$$(B.2)

where ℛ_x,y,z(α) denotes the passive rotation of angle α about the axes x, y, or z. In matrix form, this is $${ℛ}_x\left(\alpha \right)=\left(\begin{array}{ccc}1& 0& 0\\ 0& \cos \alpha & \sin \alpha \\ 0& -\sin \alpha & \cos \alpha \end{array}\right)$$ $${ℛ}_x\left(\alpha \right)=\left(\begin{array}{ccc}\cos \alpha & \sin \alpha & 0\\ -\sin \alpha & \cos \alpha & 0\\ 0& 0& 1\end{array}\right).$$

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 forelliptical orhyperbolic orbits. However, the solar mass parametermust 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.

Appendix C Acknowledgements

This work presents results from the European Space Agency (ESA) space mission Gaia. Gaia data are being processed by the Gaia Data Processing and Analysis Consortium (DPAC). Funding for the DPAC is provided by national institutions, in particular the institutions participating in the Gaia MultiLateral Agreement (MLA). The Gaia mission website is https://www.cosmos.esa.int/gaia. The Gaia archive website is https://archives.esac.esa.int/gaia.

The Gaia mission and data processing have financially been supported by, in alphabetical order by country:

• the Algerian Centre de Recherche en Astronomie, Astro-physique et Géophysique of Bouzareah Observatory;

• the Austrian Fonds zur Förderung der wissenschaftlichen Forschung (FWF) Hertha Firnberg Programme through grants T359, P20046, and P23737;

• the BELgian federal Science Policy Office (BEL-SPO) through various PROgramme de Développement d’Expériences scientifiques (PRODEX) grants of the European Space Agency (ESA), the Research Foundation Flanders (Fonds Wetenschappelijk Onderzoek) through grant VS.091.16N, the Fonds de la Recherche Scientifique (FNRS), and the Research Council of Katholieke Univer-siteit (KU) Leuven through grant C16/18/005 (Pushing AsteRoseismology to the next level with TESS, GaiA, and the Sloan Digital Sky SurvEy – PARADISE);

• the Brazil-France exchange programmes Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) and Coordenação de Aperfeicoamento de Pessoal de Nível Superior (CAPES) - Comité Français d’Evaluation de la Coopération Universitaire et Scientifique avec le Brésil (COFECUB);

• the Chilean Agencia Nacional de Investigación y Desar-rollo (ANID) through Fondo Nacional de Desarrollo Cientí-fico y Tecnológico (FONDECYT) Regular Project 1210992 (L. Chemin);

• the National Natural Science Foundation of China (NSFC) through grants 11573054, 11703065, and 12173069, the China Scholarship Council through grant 201806040200, and the Natural Science Foundation of Shanghai through grant 21ZR1474100;

• the Tenure Track Pilot Programme of the Croatian Science Foundation and the École Polytechnique Fédérale de Lausanne and the project TTP-2018-07-1171 ‘Mining the Variable Sky’, with the funds of the Croatian-Swiss Research Programme;

• the Czech-Republic Ministry of Education, Youth, and Sports through grant LG 15010 and INTER-EXCELLENCE grant LTAUSA18093, and the Czech Space Office through ESA PECS contract 98058;

• the Danish Ministry of Science;

• the Estonian Ministry of Education and Research through grant IUT40-1;

• the European Commission’s Sixth Framework Programme through the European Leadership in Space Astrometry (ELSA) Marie Curie Research Training Network (MRTN-CT-2006-033481), through Marie Curie project PIOF-GA-2009-255267 (Space AsteroSeismology & RR Lyrae stars, SAS-RRL), and through a Marie Curie Transfer-of-Knowledge (ToK) fellowship (MTKD-CT-2004-014188); the European Commission’s Seventh Framework Programme through grant FP7-606740 (FP7-SPACE-2013-1) for the Gaia European Network for Improved data User Services (GENIUS) and through grant 264895 for the Gaia Research for European Astronomy Training (GREAT-ITN) network;

• the European Cooperation in Science and Technology (COST) through COST Action CA18104 ‘Revealing the Milky Way with Gaia (MW-Gaia)’;

• the European Research Council (ERC) through grants 320360, 647208, and 834148 and through the European Union’s Horizon 2020 research and innovation and excellent science programmes through Marie Sklodowska-Curie grants 687378 (Small Bodies: Near and Far), 682115 (Using the Magellanic Clouds to Understand the Interaction of Galaxies), 695099 (A sub-percent distance scale from binaries and Cepheids – CepBin), 716155 (Structured ACCREtion Disks – SACCRED), 745617 (Our Galaxy at full HD –Gal-HD), 895174 (The build-up and fate of self-gravitating systems in the Universe), 951549 (Sub-percent calibration of the extragalactic distance scale in the era of big surveys – UniverScale), 101004214 (Innovative Scientific Data Exploration and Exploitation Applications for Space Sciences

• EXPLORE), 101004719 (OPTICON-RadioNET Pilot), 101055318 (The 3D motion of the Interstellar Medium with ESO and ESA telescopes – ISM-FLOW), and 101063193 (Evolutionary Mechanisms in the Milky waY; the Gaia Data Release 3 revolution – EMMY);

• the European Science Foundation (ESF), in the framework of the Gaia Research for European Astronomy Training Research Network Programme (GREAT-ESF);

• the European Space Agency (ESA) in the framework of the Gaia project, through the Plan for European Cooperating States (PECS) programme through contracts C98090 and 4000106398/12/NL/KML for Hungary, through contract 4000115263/15/NL/IB for Germany, through PROgramme de Développement d’Expériences scientifiques (PRODEX) grants 4000132054 for Hungary and through contract 4000132226/20/ES/CM;

• the Academy of Finland through grants 299543, 307157, 325805, 328654, 336546, and 345115 and the Magnus Ehrn-rooth Foundation;

• the French Centre National d’Études Spatiales (CNES), the Agence Nationale de la Recherche (ANR) through grant ANR-10-IDEX-0001-02 for the ‘Investissements d’avenir’ programme, through grant ANR-15-CE31-0007 for project ‘Modelling the Milky Way in the Gaia era’ (MOD4Gaia), through grant ANR-14-CE33-0014-01 for project ‘The Milky Way disc formation in the Gaia era’ (ARCHEOGAL), through grant ANR-15-CE31-0012-01 for project ‘Unlocking the potential of Cepheids as primary distance calibrators’ (UnlockCepheids), through grant ANR-19-CE31-0017 for project ‘Secular evolution of galaxies’ (SEGAL), and through grant ANR-18-CE31-0006 for project ‘Galactic Dark Matter’ (GaDaMa), the Centre National de la Recherche Scientifique (CNRS) and its SNO Gaia of the Institut des Sciences de l’Univers (INSU), its Programmes Nationaux: Cosmologie et Galaxies (PNCG), Gravitation Références Astronomie Métrologie (PNGRAM), Planétolo-gie (PNP), Physique et Chimie du Milieu Interstellaire (PCMI), and Physique Stellaire (PNPS), supported by INSU along with the Institut National de Physique (INP) and the Institut National de Physique nucléaire et de Physique des Particules (IN2P3), and co-funded by CNES; the ‘Action Fédératrice Gaia’ of the Observatoire de Paris, and the Région de Franche-Comté;

• the German Aerospace Agency (Deutsches Zentrum für Luft- und Raumfahrt e.V., DLR) through grants 50QG0501, 50QG0601, 50QG0602, 50QG0701, 50QG0901, 50QG1001, 50QG1101, 50QG1401, 50QG1402, 50QG1403, 50QG1404, 50QG1904, 50QG2101, 50QG2102, and 50QG2202, and the Centre for Information Services and High Performance Computing (ZIH) at the Technische Universität Dresden for generous allocations of computer time;

• the Hungarian Academy of Sciences through the János Bolyai Research Scholarship (G. Marton and Z. Nagy), the Lendület Programme grants LP2014-17 and LP2018-7 and the Hungarian National Research, Development, and Innovation Office (NKFIH) through grant KKP-137523 (‘Seis-moLab’);

• the Science Foundation Ireland (SFI) through a Royal Society - SFI University Research Fellowship (M. Fraser);

• the Israel Ministry of Science and Technology through grant 3-18143 and the Israel Science Foundation (ISF) through grant 1404/22;

• the Agenzia Spaziale Italiana (ASI) through contracts I/037/08/0, I/058/10/0, 2014-025-R.0, 2014-025-R.1.2015, and 2018-24-HH.0 and its addendum 2018-24-HH.1-2022 to the Italian Istituto Nazionale di Astrofisica (INAF), contract 2014-049-R.0/1/2, 2022-14-HH.0 to INAF for the Space Science Data Centre (SSDC, formerly known as the ASI Science Data Center, ASDC), contracts I/008/10/0, 2013/030/I.0, 2013-030-I.0.1-2015, and 2016-17-I.0 to the Aerospace Logistics Technology Engineering Company (ALTEC S.p.A.), INAF, and the Italian Ministry of Education, University, and Research (Ministero dell’Istruzione, dell’Università e della Ricerca) through the Premiale project ‘MIning The Cosmos Big Data and Innovative Italian Technology for Frontier Astrophysics and Cosmology’ (MITiC);

• the Netherlands Organisation for Scientific Research (NWO) through grant NWO-M-614.061.414, through a VICI grant (A. Helmi), and through a Spinoza prize (A. Helmi), and the Netherlands Research School for Astronomy (NOVA);

• the Polish National Science Centre through HARMONIA grant 2018/30/M/ST9/00311 and DAINA grant 2017/27/L/ ST9/03221 and the Ministry of Science and Higher Education (MNiSW) through grant DIR/WK/2018/12;

• the Portuguese Fundação para a Ciência e a Tecnologia (FCT) through national funds, grants 2022.06962.PTDC and 2022.03993.PTDC, and work contract DL 57/2016/CP1364/ CT0006, grants UIDB/04434/2020 and UIDP/04434/2020 for the Instituto de Astrofisica e Ciências do Espaço (IA), grants UIDB/00408/2020 and UIDP/00408/2020 for LASIGE, and grants UIDB/00099/2020 and UIDP/00099/ 2020 for the Centro de Astrofisica e Gravitação (CENTRA);

• the Slovenian Research Agency through grant P1-0188;

• the Spanish Ministry of Economy (MINECO/FEDER, UE), the Spanish Ministry of Science and Innovation (MCIN), the Spanish Ministry of Education, Culture, and Sports, and the Spanish Government through grants BES-2016-078499, BES-2017-083126, BES-C-2017-0085, ESP2016-80079-C2-1-R, FPU16/03827, RTI2018-095076-B-C22, PID2021-122842OB-C22, PDC2021-121059-C22, and TIN2015-65316-P (‘Computatión de Altas Prestaciones VII’), the Juan de la Cierva Incorporatión Programme (FJCI-2015-2671 and IJC2019-04862-I for F. Anders), the Severo Ochoa Centre of Excellence Programme (SEV2015-0493) and MCIN/AEI/10.13039/501100011033/ EU FEDER and Next Generation EU/PRTR (PRTR-C17.I1); the European Union through European Regional Development Fund ‘A way of making Europe’ through grants PID2021-122842OB-C21, PID2021-125451NA-I00, CNS2022-13523 and RTI2018-095076-B-C21, the Institute of Cosmos Sciences University of Barcelona (ICCUB, Unidad de Excelencia ‘Mar a de Maeztu’) through grant CEX2019-000918-M, the University of Barcelona’s official doctoral programme for the development of an R+D+i project through an Ajuts de Personal Investigador en For-mació (APIF) grant, the Spanish Virtual Observatory project funded by MCIN/AEI/10.13039/501100011033/ through grant PID2020-112949GB-I00; the Centro de Investigatión en Tecnologias de la Informatión y las Comunicaciones (CITIC), funded by the Xunta de Galicia through the collaboration agreement to reinforce CIGUS research centers, research consolidation grant ED431B 2021/36 and scholarships from Xunta de Galicia and the EUESF ED481A-2019/155 and ED481A 2021/296; the Red Española de Supercomputación (RES) computer resources at MareNostrum, the Barcelona Supercomputing Centre -Centro Nacional de Supercomputación (BSC-CNS) through activities AECT-2017-2-0002, AECT-2017-3-0006, AECT-2018-1-0017, AECT-2018-2-0013, AECT-2018-3-0011, AECT-2019-1-0010, AECT-2019-2-0014, AECT-2019-3-0003, AECT-2020-1-0004, and DATA-2020-1-0010, the Departament d’nnovació, Universitats i Empresa de la Generalitat de Catalunya through grant 2014-SGR-1051 for project ‘Models de Programació i Entorns d’Execució Parallels’ (MPEXPAR), and Ramon y Cajal Fellowships RYC2018-025968-I, RYC2021-031683-I and RYC2021-033762-I, funded by MICIN/AEI/10.13039/501100011033 and by the European Union NextGenerationEU/PRTR and the European Science Foundation (‘Investing in your future’); the Port d’Informació Cientifica (PIC), through a collaboration between the Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT) and the Institut de F sica d’Altes Energies (IFAE), supported by the call for grants for Scientific and Technical Equipment 2021 of the State Program for Knowledge Generation and Scientific and Technological Strengthening of the R+D+i System, financed by MCIN/AEI/ 10.13039/501100011033 and the EU NextGeneration/PRTR (Hadoop Cluster for the comprehensive management of massive scientific data, reference EQC2021-007479-P);

• the Swedish National Space Agency (SNSA/ Rymdstyrelsen);

• the Swiss State Secretariat for Education, Research, and Innovation through the Swiss Activités Nationales Complémentaires and the Swiss National Science Foundation through an Eccellenza Professorial Fellowship (award PCEFP2_194638 for R. Anderson);

• the United Kingdom Particle Physics and Astronomy Research Council (PPARC), the United Kingdom Science and Technology Facilities Council (STFC), and the United Kingdom Space Agency (UKSA) through the following grants to the University of Bristol, Brunel University London, the Open University, the University of Cambridge, the University of Edinburgh, the University of Leicester, the Mullard Space Sciences Laboratory of University College London, and the United Kingdom Rutherford Appleton Laboratory (RAL):

  PP/D006503/1, PP/D006511/1, PP/ D006546/1, PP/D006570/1, PP/D006791/1, ST/I000852/1,
  ST/J005045/1,   ST/K00056X/1,   ST/K000209/1,
  ST/K000756/1,   ST/K000578/1,   ST/L002388/1,
  ST/L006553/1,   ST/L006561/1,   ST/N000595/1,
  ST/N000641/1,   ST/N000978/1,   ST/N001117/1,
  ST/S000089/1,   ST/S000976/1,   ST/S000984/1,
  ST/S001123/1, ST/S001948/1, ST/S001980/1, ST/S002103/1,
  ST/V000969/1,   ST/W002469/1,   ST/W002493/1,
  ST/W002671/1,   ST/W002809/1,   EP/V520342/1,
  ST/X00158X/1,   ST/X001601/1,   ST/X001636/1,
  ST/X001687/1, ST/X002667/1, ST/X002683/1 and
  ST/X002969/1.

The Gaia project and data processing have made use of:

• the Set of Identifications, Measurements, and Bibliography for Astronomical Data (SIMBAD, Wenger et al. 2000), the ‘Aladin sky atlas’ (Bonnarel et al. 2000; Boch & Fer-nique 2014), and the VizieR catalogue access tool (Ochsenbein et al. 2000), all operated at the Centre de Données astronomiques de Strasbourg (CDS);

• the National Aeronautics and Space Administration (NASA) Astrophysics Data System (ADS);

• the SPace ENVironment Information System (SPENVIS), initiated by the Space Environment and Effects Section (TEC-EES) of ESA and developed by the Belgian Institute for Space Aeronomy (BIRA-IASB) under ESA contract through ESA’s General Support Technologies Programme (GSTP), administered by the BELgian federal Science Policy Office (BELSPO);

• the software products TOPCAT, STIL, and STILTS (Taylor 2005, Taylor 2006);

• Matplotlib (Hunter 2007);

• IPython (Pérez & Granger 2007);

• Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration 2018);

• R (R Core Team 2013);

• the HEALPix package (Górski et al. 2005, http://healpix.sourceforge.net/);

• Vaex (Breddels & Veljanoski 2018);

• the Hipparcos-2 catalogue (van Leeuwen 2007). The Hip-parcos and Tycho catalogues were constructed under the responsibility of large scientific teams collaborating with ESA. The Consortia Leaders were Lennart Lindegren (Lund, Sweden: NDAC) and Jean Kovalevsky (Grasse, France: FAST), together responsible for the Hipparcos Catalogue; Erik Høg (Copenhagen, Denmark: TDAC) responsible for the Tycho Catalogue; and Catherine Turon (Meudon, France: INCA) responsible for the Hipparcos Input Catalogue (HIC);

• the Tycho-2 catalogue (Høg et al. 2000), the construction of which was supported by the Velux Foundation of 1981 and the Danish Space Board;

• The Tycho double star catalogue (TDSC, Fabricius et al. 2002), based on observations made with the ESA Hippar-cos astrometry satellite, as supported by the Danish Space Board and the United States Naval Observatory through their double-star programme;

• data products from the Two Micron All Sky Survey (2MASS, Skrutskie et al. 2006), which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center (IPAC) / California Institute of Technology, funded by the National Aeronautics and Space Administration (NASA) and the National Science Foundation (NSF) of the USA;

• the ninth data release of the AAVSO Photometric All-Sky Survey (APASS, Henden et al. 2016), funded by the Robert Martin Ayers Sciences Fund;

• the first data release of the Pan-STARRS survey (Chambers et al. 2016; Magnier et al. 2020a; Waters et al.2020; Magnier et al. 2020 c,b; Flewelling et al. 2020). The Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, the Queen’s University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration (NASA) through grant NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation through grant AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation;

• the second release of the Guide Star Catalogue (GSC2.3, Lasker et al. 2008). The Guide Star Catalogue II is a joint project of the Space Telescope Science Institute (STScI) and the Osservatorio Astrofisico di Torino (OATo). STScI is operated by the Association of Universities for Research in Astronomy (AURA), for the National Aeronautics and Space Administration (NASA) under contract NAS5-26555. OATo is operated by the Italian National Institute for Astrophysics (INAF). Additional support was provided by the European Southern Observatory (ESO), the Space Telescope European Coordinating Facility (STECF), the International GEMINI project, and the European Space Agency (ESA) Astrophysics Division (nowadays SCI-S);

• the eXtended, Large (XL) version of the catalogue of Positions and Proper Motions (PPM-XL, Roeser et al. 2010);

• data products from the Wide-field Infrared Survey Explorer (WISE), which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, and NEO-WISE, which is a project of the Jet Propulsion Laboratory/California Institute of Technology. WISE and NEO-WISE are funded by the National Aeronautics and Space Administration (NASA);

• the first data release of the United States Naval Observatory (USNO) Robotic Astrometric Telescope (URAT-1, Zacharias et al. 2015);

• the fourth data release of the United States Naval Observatory (USNO) CCD Astrograph Catalogue (UCAC-4, Zacharias et al. 2013);

• the sixth and final data release of the Radial Velocity Experiment (RAVE DR6, Steinmetz et al. 2020a,b). Funding for RAVE has been provided by the Leibniz Institute for Astrophysics Potsdam (AIP), the Australian Astronomical Observatory, the Australian National University, the Australian Research Council, the French National Research Agency, the German Research Foundation (SPP 1177 and SFB 881), the European Research Council (ERC-StG 240271 Galac-tica), the Istituto Nazionale di Astrofisica at Padova, the Johns Hopkins University, the National Science Foundation of the USA (AST-0908326), the W.M. Keck foundation, the Macquarie University, the Netherlands Research School for Astronomy, the Natural Sciences and Engineering Research Council of Canada, the Slovenian Research Agency, the Swiss National Science Foundation, the Science & Technology Facilities Council of the UK, Opticon, Strasbourg Observatory, and the Universities of Basel, Groningen, Heidelberg, and Sydney. The RAVE website is at https://www.rave-survey.org/;

• the first data release of the Large sky Area Multi-Object Fibre Spectroscopic Telescope (LAMOST DR1, Luo et al. 2015);

• the K2 Ecliptic Plane Input Catalogue (EPIC, Huber et al. 2016);

• the ninth data release of the Sloan Digitial Sky Survey (SDSS DR9, Ahn et al. 2012). Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the United States Department of Energy Office of Science. The SDSS-III website is http://www.sdss3.org/. SDSS-III is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS-III Collaboration including the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Instituto de Astrofísica de Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University;

• the thirteenth release of the Sloan Digital Sky Survey (SDSS DR13, Albareti et al. 2017). Funding for SDSS-IV has been provided by the Alfred P. Sloan Foundation, the United States Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is https://www.sdss.org/. SDSS-IV is managed by the Astrophys-ical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University;

• the second release of the SkyMapper catalogue (SkyMap-per DR2, Onken et al. 2019, Digital Object Identifier 10.25914/5ce60d31ce759). The national facility capability for SkyMapper has been funded through grant LE130100104 from the Australian Research Council (ARC) Linkage Infrastructure, Equipment, and Facilities (LIEF) programme, awarded to the University of Sydney, the Australian National University, Swinburne University of Technology, the University of Queensland, the University of Western Australia, the University of Melbourne, Curtin University of Technology, Monash University, and the Australian Astronomical Observatory. SkyMapper is owned and operated by The Australian National University’s Research School of Astronomy and Astrophysics. The survey data were processed and provided by the SkyMapper Team at the Australian National University. The SkyMapper node of the All-Sky Virtual Observatory (ASVO) is hosted at the National Computational Infrastructure (NCI). Development and support the SkyMapper node of the ASVO has been funded in part by Astronomy Australia Limited (AAL) and the Australian Government through the Commonwealth’s Education Investment Fund (EIF) and National Collaborative Research Infrastructure Strategy (NCRIS), particularly the National eResearch Collaboration Tools and Resources (NeCTAR) and the Australian National Data Service Projects (ANDS);

• the Gaia-ESO Public Spectroscopic Survey (GES, Gilmore et al. 2022; Randich et al. 2022). The Gaia-ESO Survey is based on data products from observations made with ESO Telescopes at the La Silla Paranal Observatory under programme ID 188.B-3002. Public data releases are available through the ESO Science Portal. The project has received funding from the Leverhulme Trust (project RPG-2012-541), the European Research Council (project ERC-2012-AdG 320360-Gaia-ESO-MW), and the Istituto Nazionale di Astrofisica, INAF (2012: CRA 1.05.01.09.16; 2013: CRA 1.05.06.02.07).

The GBOT programme (GBOT) uses observations collected at (i) the European Organisation for Astronomical Research in the Southern Hemisphere (ESO) with the VLT Survey Telescope (VST), under ESO programmes 092.B-0165, 093.B-0236, 094.B-0181, 095.B-0046, 096.B-0162, 097.B-0304, 098.B-0030, 099.B-0034, 0100.B-0131, 0101.B-0156, 0102.B-0174, 0103.B-0165, 0104.B-0081, 0106.20ZA.001 (OmegaCam), 0106.20ZA.002 (FORS2), 0108.21YF; and under INAF programs 110.256C, 112.266Q; and (ii) the Liverpool Telescope, which is operated on the island of La Palma by Liverpool John Moores University in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias with financial support from the United Kingdom Science and Technology Facilities Council, and (iii) telescopes of the Las Cumbres Observatory Global Telescope Network.

In addition to the currently active DPAC (and ESA science) authors of the peer-reviewed papers accompanying the data release, there are large numbers of former DPAC members who made significant contributions to the (preparations of the) data processing. In addition to the DPAC consortium, past and present, there are numerous people, mostly in ESA and in industry, who have made or continue to make essential contributions to Gaia, for instance those employed in science and mission operations or in the design, manufacturing, integration, and testing of the spacecraft and its modules, subsystems, and units. Many of those will remain unnamed yet spent countless hours, occasionally during nights, weekends, and public holidays, in cold offices and dark clean rooms. They are acknowledged in the Gaia Documentation

References

All Tables

All Figures

	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).
In the text

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

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.

In the text

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

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 Gmag < 11.5, the distribution is not Gaussian. This is clearer in Fig. 6. For Gmag ≥ 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.

In the text

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

In the text

	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.
In the text

	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 quan-tiles at the 1σ level.
In the text

	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 (OICRF) is offset from the inertial equinox γJ2000. The relevant numbers are given in Table 1.
In the text

	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.
In the text

	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.
In the text

	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.
In the text

	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.
In the text

	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.
In the text

	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.
In the text

	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 19822 Gaia solutions, 261 of which are outside the plot boundaries.
In the text

	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.
In the text

	Fig. 18 Relative uncertainty of the reference orbits computed by the JPL before propagation to the mean Gaia epoch. The bins add up to 100% for each group.
In the text

	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.
In the text

	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.
In the text

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. 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.

In the text

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

	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.
In the text

	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.
In the text

	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.
In the text

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 σα < 1.6 × 10−10 au.

In the text

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