© ESO 2018

1. Introduction

The main goal of the Gaia mission, a European Space Agency (ESA) satellite launched in December 2013, is the production of a five-parameter astrometric catalog (i.e., including positions, parallaxes and the two components of the proper motions) at the 10–1000 μarcsecond-level (μas) for about one billion stars of our Galaxy in the magnitude range from 3 to 20.7 (Gaia Collaboration 2016).

To this end, the satellite has been designed as a scanning telescope that sweeps continuously and repeatedly the entire celestial sphere during the five years of its foreseen mission duration. The target accuracy can be reached by averaging on the ∼10³ observations per object, each at the ∼0.1 − 1 mas level, and it relies on the self-calibration capability of the instrument, as well as on complementary photometric and spectroscopic data collected on board. The latter will also make it possible to include in the Gaia catalog a multiband spectro-photometric classification and the radial velocities of the objects brighter than G ≈ 16.2 (Gaia Collaboration 2016).

The sky is scanned according to a predetermined nominal scanning law (NSL) whose parameters are fixed at the beginning of the operational phase and actively controlled by the micro-thrusters of the on-board Attitude and Orbit Control Systems (AOCS). The NSL is characterized by constant spin and precession rates, which implies that objects cross the focal plane at quasi-constant speed. As a consequence, in principle the duration of a single observation is the same for any star, and therefore the single-measurement accuracy is strongly magnitude-dependent (Fig. 1). Moreover, the scanning law parameters induce a complete scan of the celestial sphere every six months, but with a non-uniform coverage in terms of the number of single transits over a specific region, which implies that the number of observations of each object depends on its coordinates (Fig. 2). Therefore, the final astrometric accuracy of a specific object mainly depends on its magnitude and location on the sky.

Fig. 1. Approximate estimation of the Gaia single-measurement along- and across-scan accuracy (AL and AC) as function of the G magnitude. The oscillating shape up to mag 12 is due to the use of CCD gates in order to avoid saturation for bright objects. The estimation is the result of a best-fit model that averages over the whole Astrometric Field (AF).

Fig. 2. Frequency of observations as function of equatorial coordinates due to Gaia scanning law (blue ≈50–yellow ≈500).

Contrary to what usually happens for a space mission, for this mission ESA was in charge of the full development of the satellite, including the payload. The scientific community instead, organized in the Gaia Data Processing and Analysis Consortium (DPAC; Mignard & Drimmel 2007) funded by the national space agencies, has been in charge of the establishment of the data reduction pipelines of the mission. The software is developed by nine Coordination Units (CUs) while six Data Processing Centres (DPCs) spread all over Europe are commissioned of managing and running the pipelines.

The realization of the complete catalog is a complicated process in which the definition of the global reference frame is realized by a procedure called global astrometric sphere reconstruction. The latter operates with global astrometric data reduction techniques on a subset of ∼10⁷ − 10⁸ “primary stars”, similar in number to that of other modern global astrometric catalogs (Zacharias et al. 2013; Qi et al. 2015) and mainly selected at the bright end of the Gaia objects. This allows to set the astrometric parameters of these objects with respect to a standard reference system called Barycentric Celestial Reference System (BCRS; Kaplan 2005) thus providing the first materialization of such a frame of reference. The refinement of the instrument calibration needed to achieve the required catalog accuracy is carried out as part of the primary stars’ processing. Once the reference system and the calibration parameters have been established, the measurements of the remaining stellar objects (secondary stars) can be reduced by considering only their astrometric parameters as unknowns. This allows to attach them to the primary stars and therefore to densify the reference frame. Basically, a star can be included in the primaries subset when its astrometric model can be described by the classical 5 parameters. Multiple stars with too short a period, or stars with a variability too large are examples of objects that cannot belong to the primaries. The two-steps process described above, namely the one including the global sphere reconstruction and the reduction of the secondary stars, is realized within the CU3, “Core Processing”, by a pipeline called Astrometric Global Iterative Solution (AGIS; O’Mullane et al. 2011; Lindegren et al. 2012).

The global sphere reconstruction, and the materialization of a celestial reference system, has an absolute character (in the sense that it defines the reference system instead of giving the coordinates of the stars with respect to an already existing one) which implies an intrinsic difficulty at insuring the correctness of the astrometric parameters and at the same time the risk of propagating these errors everywhere in Astrophysics. This issue is well known to the scientific community, and significant effort is usually paid to both the tasks of internal verifications and of cross-checking comparison with different datasets (see, e.g. Lindegren et al. 2016; Casertano et al. 2017; Makarov et al. 2017; Frouard et al. 2018). In the case of the Gaia forerunner, HIPPARCOS, the astrometric community provided the final catalog only after having compared two sphere reconstructions realized by two independent consortia. In the case of Gaia the task is so big that it is not feasible to repeat the same approach, but the DPAC decided to constitute an Astrometric Verification Unit (AVU) within the CU3 with the goal of replicating in an independent way three specific tasks of particular importance for the sphere reconstruction, namely the Astrometric Instrument Model (AIM; Busonero 2012; Busonero et al. 2014) for the instrument and focal plane calibration, the Basic Angle Monitoring (BAM; Riva et al. 2014) which has to determine the variations of the lines of sight of the double telescope, and the Global Sphere Reconstruction (GSR) whose aim is to provide an independent sphere reconstruction and to compare it with the AGIS one.

Since its goal is to allow a cross-checking of the sphere solution, in the sense of the reference system determination, GSR does not replicate the entire AGIS pipeline, but is limited to the processing of primary sources. Moreover, GSR depends on AGIS both for the determination of such sources – even though it has the capability of providing an independent selection – and for the rejection of time intervals with noisy attitude which must not enter in the sphere reconstruction. On the other hand, GSR is in charge of the comparison task between its own solution and that of the AGIS pipeline.

In this paper we adopt the following conventions and notations:

1. in general, bold upright letters refers to 3D Euclidean vectors (x) while for basis unit vectors the notation e_â, a = 1, 2, 3 is used;

2. four-vectors are indicated by bold italic letters (x) or in index notation with Greek indexes, that is x^α, α = 0, 1, 2, 3 with 0 referring to the time coordinate;

3. the signature of the metric g_αβ is +2;

4. the symbol η_αβ denotes the Minkowskian metric diag{−1, 1, 1, 1};

5. tetrad unit vectors are identified by e_α̂, α = 0, 1, 2, 3, while e_â are the tetrad spatial axes;

6. the proper time of an observer is indicated with the Greek letter τ, while t is the coordinate time;

7. we adopt the symbol ∂_α as a shorthand notation for the derivative ∂/∂x^α.

2. Modeling the observations for the global astrometric sphere reconstruction

2.1. Geometric characterization of the observable

In a purely Euclidean geometric view, the NSL gives at each instant the orientation of the Satellite Reference System (SRS) that is of the reference triad {e_â}, a = 1, 2, 3 attached to the satellite, whose origin coincides with the barycentre of Gaia, with respect to the inertial reference system of the catalog. The position and orientation of the satellite results from the combination of three independent movements: the orbital motion of the satellite, which follows a Lissajous trajectory around the L2 point of the Sun–Earth system, the constant-rate satellite spin around the axis ${\mathbf{e}}_{\widehat{3}}$, and the precession of this axis around the Sun–satellite direction, with constant rate and keeping a constant solar aspect angle ξ. Thus in an equatorial coordinate system the SRS is completely determined by five parameters: the heliotropic angles (ξ,ν,Ω), the obliquity of the Ecliptic ε and the longitude of the nominal Sun λ_s (Fig. 3). Obviously the instantaneous orientation of the satellite depends also on the initial conditions of the NSL. Since the solar aspect angle ξ is set to 45°, the actual scanning law is fully determined by setting the remaining three degrees of freedom, namely ν, Ω and λ_s, at a given reference time.

Fig. 3. Parametrization of the Gaia Nominal Scanning Law with respect to an equatorial reference system $\{{\mathbf{E}}_{\widehat{1}},{\mathbf{E}}_{\widehat{2}},{\mathbf{E}}_{\widehat{3}}\}$. NS is the Nominal Sun.

In this geometrical model the two Fields-of-View (FoVs) of Gaia lie on the satellite’ scanning plane ${\mathbf{e}}_{\widehat{1}}-{\mathbf{e}}_{\widehat{2}}$ and are pointing symmetrically with respect to ${\mathbf{e}}_{\widehat{1}}$, separated by an (approximately) constant Basic Angle (BA) Γ. If r is the position vector of a point-like source, in the SRS the basic observables of Gaia can be represented by its along-scan measurement (AL) that is the abscissa ϕ, which is the angle between the projection of r on the instantaneous scanning plane and ${\mathbf{e}}_{\widehat{1}}$, and the across-scan measurement (AC) ζ (Fig. 4).

Fig. 4. Representation of the AL (ϕ) and AC (ζ) Gaia measurements with respect to the satellite attitude $\{{\mathbf{e}}_{\widehat{1}},{\mathbf{e}}_{\widehat{2}},{\mathbf{e}}_{\widehat{3}}\}$.

These quantities can be represented in terms of the direction cosines of r with respect to the SRS

$$\begin{array}{cc}\hfill cos\varphi & =\frac{cos{\psi }_{(\widehat{1},\mathbf{r})}}{\sqrt{1-{cos}^2{\psi }_{(\widehat{3},\mathbf{r})}}}\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill sin\zeta & =cos{\psi }_{(\widehat{3},\mathbf{r})},\hfill \end{array}$$(2)

while the direction cosine of r with respect to the axis e_â is

$$\begin{array}{c}\hfill cos{\psi }_{(\widehat{a},\mathbf{r})}=\frac{{\mathbf{e}}_{\widehat{a}}·\mathbf{r}}{{(\mathbf{r}·\mathbf{r})}^{1/2}}.\end{array}$$(3)

It is worth stressing that these modeling equations leave a potential ambiguity in the sign of the abscissa ϕ, which should be negative when the observation lies in FoV1 and positive in FoV2. However, the FoV of each observation is specified in the Gaia data, which removes this ambiguity.

2.2. Unknown parameters in the Gaia sphere reconstruction

In general e_â and r will be functions of the time of observation, of the attitude parameters and of the source parameters respectively. The latter are the classical parallax, positions and proper motions (ϖ,α,δ,μ_α,μ_δ) because in this case all the sources are primaries, while the attitude parameters have to be counted among the unknowns too, because they cannot be determined independently at the level needed for the Gaia accuracy. The attitude parameters and the explicit form of the functions, however, will depend on the specific astrometric model adopted to represent the observation.

Gaia is also a self-calibrating mission, in the sense that there are no metrologic instruments, with the only exception of the BAM, dedicated to the calibration. All the calibrations except for the monitoring of the Basic Angle variations are done using the Gaia observations only. This task is accomplished in two different ways. One is a daily calibration done by the two independent pipelines, First Look (FL) and, in the AVU subsystem, AIM; the other one is embedded in the sphere solution. The daily calibration is needed to monitor the instrument, using science data to trace directly the instrument response thanks to the repeated measurements of stars over the field, but as long as this procedure relies on a daily approach, it cannot reach the accuracy needed by the sphere reconstruction. Therefore, as it will be better specified in Sect. 2.6, AGIS and GSR have to introduce appropriate calibration parameters as additional unknowns of their models.

Finally, there exist also another class of unknowns called global parameters. Such name derives from their presence in all the rows of the system of equations that is built from the Gaia observations, as it will be explained in Sect. 3.1, and typically their number is the smallest with respect to all the other classes.

By indicating with x^S, x^A, x^C and x^G the list of source, attitude, calibration and global parameters respectively, Eqs. (1) and (2) can be formally written as

$$\begin{array}{cc}\hfill cos\varphi & =f_{\varphi }({\mathbf{x}}^{\mathrm{S}},{\mathbf{x}}^{\mathrm{A}},{\mathbf{x}}^{\mathrm{C}},{\mathbf{x}}^{\mathrm{G}})\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill sin\zeta & =f_{\zeta }({\mathbf{x}}^{\mathrm{S}},{\mathbf{x}}^{\mathrm{A}},{\mathbf{x}}^{\mathrm{C}},{\mathbf{x}}^{\mathrm{G}}).\hfill \end{array}$$(5)

2.3. Gaia accuracy and the relativistic observable

The target accuracy of Gaia is at the sub-mas level at least, moreover the catalog will be released in the BCRS and linked to the International Celestial Reference Frame (ICRF) defined by VLBI observations (Fey et al. 2015)¹. For these reasons it is necessary to use a mathematical model of the observations based on General Relativity. Basically, the task of this model is to provide a relativistically consistent formula to replace Eq. (3). This requires the relativistic equivalent of the photon’s incoming direction (r), of the satellite attitude (the SRS), and of the definition of scalar product.

In General Relativity the trajectory of a light ray connecting a source with the observer is defined by a set of four-dimensional events x^α(s), for which the condition g_αβdx^αdx^β = 0 holds, and is thus called null geodesic, parametrized by a generic affine parameter s. This trajectory can be found by integrating the geodesic equations

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2{x}^{\alpha }}{\mathrm{d}{s}^2}+{\mathrm{\Gamma }}_{\mu \nu }^{\alpha }\frac{\mathrm{d}{x}^{\mu }}{\mathrm{d}s}\frac{\mathrm{d}{x}^{\nu }}{\mathrm{d}s}=0,\end{array}$$(6)

where

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mu \nu }^{\alpha }=\frac{1}{2}{g}^{\alpha \rho }({\partial }_{\mu }{g}_{\nu \rho }+{\partial }_{\nu }{g}_{\rho \mu }-{\partial }_{\rho }{g}_{\mu \nu }),\end{array}$$(7)

and the role of r is usually played by the value of the tangent to the null geodesic at the point of observation, namely by the four-vector k^α = dx^α/ds evaluated at that point. The photon’s incoming direction computed in this way obviously depends on the relativistic effects induced by the metric and thus by the gravitational effects of the massive bodies, such as, for example, the so-called light deflection.

Any measurement is always defined with respect to a specific observer, identified with its four-velocity u^α = dx^α/dτ, which induces a natural 3+1 splitting of the spacetime, that is the identification of the space and of the time respectively associated to such observer. Specifically, u^α is tangent to the time axis of the observer, represented by its worldline, and the 3D space is the subspace orthogonal to u^α. Any spatial measurement, therefore, is performed by means of the operator h_αβ = g_αβ + u_αu_β which can project the appropriate four-dimensional quantities onto this 3D subspace, and is thus used as a replacement of the usual 3D dot product. For example, if in Euclidean geometry the angle α between two vectors r₁ and r₂ satisfies the relation cosα = r₁ ⋅ r₂/(r₁r₂), with r = (r⋅r)^1/2, the corresponding general relativistic expression will be

$$\begin{array}{c}\hfill cos\alpha =\frac{h_{\alpha \beta }{k}_1^{\alpha }{k}_2^{\beta }}{{\left(h_{\mu \nu }{k}_1^{\mu }{k}_1^{\nu }\right)}^{1/2}{\left(h_{\rho \sigma }{k}_2^{\rho }{k}_2^{\sigma }\right)}^{1/2}}·\end{array}$$(8)

It is worth noticing that in this way the aberration effects are automatically accounted for by the inclusion of the observer’s four-velocity in h_αβ.

Finally, one standard way to handle the SRS is based on the so-called tetrad formalism (de Felice & Bini 2010), which is based on the possibility of identifying a locally Lorentzian reference system associated to a specific observer $u_{\mathrm{s}}^{\alpha }$, that is a set of four-dimensional axes e_α̂, α = 0, …, 3 called tetrad defined by the conditions

$$\begin{array}{c}\hfill g_{\mu \nu }{e}_{\widehat{\alpha }}^{\mu }{e}_{\widehat{\beta }}^{\nu }={\eta }_{\alpha \beta },e_{\widehat{0}}^{\alpha }\equiv u_{\mathrm{s}}^{\alpha }.\end{array}$$(9)

In this way, we are inducing the same 3+1 splitting of above with the addition of a set of three four-dimensional axes {e_â}, a = 1, 2, 3 belonging to the 3D space of u^α. Since by definition $g_{\mu \nu }{u}^{\mu }{\mathit{e}}_{\widehat{a}}^{\nu }=g_{\mu \nu }{\mathit{e}}_{\widehat{0}}^{\mu }{\mathit{e}}_{\widehat{a}}^{\nu }=0$ and $g_{\mu \nu }{\mathit{e}}_{\widehat{a}}^{\mu }{\mathit{e}}_{\widehat{b}}^{\nu }={\eta }_{\mathit{ab}}={\delta }_{\mathit{ab}}$, these axes are spatial (in the sense that they are orthogonal to $u_{\mathrm{s}}^{\alpha }$ and thus lie on its 3D space) and orthonormal, so they provide the replacement for the Euclidean triad of Eq. (3) and can be used to define the attitude (SRS) of the satellite. The orientation of the spatial axes, in fact, is constrained only by the orthonormality condition, which implies that any spatial 3D rotation brings to a new tetrad with a different spatially oriented triad.

By combining all the above considerations we can write the general relativistic expression for the direction cosine with respect to a general spatial axis of the tetrad as

$$\begin{array}{c}\hfill cos{\psi }_{(\widehat{a},\mathit{k})}=\frac{h_{\alpha \beta }{e}_{\widehat{a}}^{\alpha }{k}^{\beta }}{{\left(h_{\mu \nu }{k}^{\mu }{k}^{\nu }\right)}^{1/2}}·\end{array}$$(10)

On the other hand, the AL and AC measurements of Eqs. (1) and (2) are defined in the 3D subspace of $u_{\mathrm{s}}^{\alpha }$ as long as the direction cosines are the relativistic ones, therefore they do not need any relativistic replacement.

We can thus summarize the general procedure applied to build the relativistic model of the Gaia observations used in GSR with the following list:

1. integrate the geodesic equations (6) for the given metric (the BCRS one, in the case of Gaia) and find an expression of k^α at the point of observation as function of the source’s coordinates;

2. find the appropriate expression $u_{\mathrm{s}}^{\alpha }$ of the barycentric motion of Gaia in the BCRS;

3. use the above four-velocity and the BCRS metric to compute the spatial projector h_αβ;

4. find the Gaia relativistic tetrad by a two-steps procedure:

   • (a)

     use Eq. (9) to define a local tetrad, whose origin is comoving with the barycentre of Gaia and whose spatial axes are kinematically parallel to those of the BCRS (“boosted tetrad”);

   • (b)

     3D rotate the spatial axes to make them coincide with the satellite orientation, thus realizing a tetrad associated to the Gaia barycentric and attitude motion $e_{\widehat{0}},e_{\widehat{a}}$;

5. use k^α, h_αβ and ${\mathit{e}}_{\widehat{a}}^{\alpha }$ in Eq. (10) to compute the needed direction cosines;

6. use Eqs. (1) and (2) to compute the Gaia measurements.

Within this general framework, the accuracy needed for the astrometric model is set by that of the Gaia measurements and of its final catalog. In practice it is useful to link the two by means of the post-Newtonian “bookkeeping” (Will 1993). In this way, considering that the typical velocities in the Solar System are ≃30 km s⁻¹, the so-called 1PN order of (ν/c)² would correspond to a ∼ 10⁻⁸ rad in angular accuracy, that is to the mas level, allowing one to set the befitting model accuracy at the 1.5PN order (ν/c)³, corresponding to 10⁻¹² rad, or ∼0.1 μas.

Actually, this has to be considered just a first, although convenient, guess. Things can be much different in the real case, and the final accuracy strongly depends on other factors, like the geometry of the observations or the observation strategy. For example, in a global problem like the sphere solution, for large regions of the sky one can safely rely on less stringent accuracy requirements, whereas relative astrometric observations close to Solar System planets would surely need a more accurate model. In this context, we also stress that Gaia-level final accuracies have also been reached in the context of differential astrometry (Brown et al. 2018).

GSR is designed to be flexible in regard to the relativistic modeling of the observable, and different relativistic models are under development to add value to the verification purposes of this pipeline. We are dealing specifically with this issue in a forthcoming publication (Vecchiato et al., in prep.) whereas in the next subsection we are concentrating on the implementation of GSR2, namely the present version of our pipeline.

2.4. Integration of the geodesic equation in the current pipeline implementation

The geodesic equations are integrated in the single-body PPN-Schwarzschild metric, where the tangent to the null geodesic k^α is a function of three constants of motion E_*, Λ_*, and λ_* characteristic of the geodesic connecting the observer with the observed object in the Schwarzschild metric², of the satellite position (r_s,θ_s,ϕ_s), and of the PPN-γ parameter:

$$\begin{array}{c}\hfill k^{\alpha }=k^{\alpha }(E_{\ast },{\mathrm{\Lambda }}_{\ast },{\lambda }_{\ast },r_{\mathrm{s}},{\theta }_{\mathrm{s}},{\varphi }_{\mathrm{s}},\gamma ).\end{array}$$(11)

First of all, it is worth stressing that the accuracy of the Gaia measurement allows for the estimation of the γ parameter as a by-product of the global sphere reconstruction (Vecchiato et al. 2003). This implies that this parameter can be treated as an unknown belonging to the vector of global parameters x^G.

Moreover, the constants of motion can be written as functions of the astrometric unknowns at the time of observation, by means of the same principle used in Vecchiato et al. (2003) and references therein, which in summary eliminates the dependence on E_* and provides two equations that implicitly relate the stellar position with Λ_* and λ_*

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill f_r({\varpi }_{\ast },r_{\mathrm{s}},{\mathrm{\Lambda }}_{\ast },\gamma )& =f_{\theta }({\theta }_{\ast },{\theta }_{\mathrm{s}},{\lambda }_{\ast }/{\mathrm{\Lambda }}_{\ast })\\ \hfill {\overline{f}}_{\theta }({\theta }_{\ast },{\theta }_{\mathrm{s}},{\mathrm{\Lambda }}_{\ast },{\lambda }_{\ast })& =f_{\varphi }({\varphi }_{\ast }-{\varphi }_{\mathrm{s}})\end{array}\end{array}$$(12)

where ϖ = a/r_* is the parallax, a = 1 AU is the parallaxes baseline, θ_* = π/2 − δ and ϕ_* ≡ α. Finally, the proper motion can be easily included in the model simply by considering the integration limits on θ_* and ϕ_* as function of the (coordinate) time, that is

$$\begin{array}{ccc}\hfill \delta & \equiv \hfill & \hfill \delta (t)=\delta (t_0)+{\mu }_{\delta }\times (t-t_0)\end{array}$$(13)

$$\begin{array}{ccc}\hfill \alpha & \equiv \hfill & \hfill \alpha (t)=\alpha (t_0)+{\mu }_{\alpha }\times (t-t_0)\end{array}$$(14)

where t₀ is the epoch of the astrometric catalog or an equivalent reference time. The above formulae neglect the effect on spherical coordinates at second order in proper motions (see e.g. Green 1985, p. 264), which are negligible in most cases. For example, in the catalog of more than 900 000 simulated stars considered for the demonstration run of this paper, this effect is larger than 5 μas after 5 years for just 22 objects in α and 5 objects in δ.

It is clear from this summary that, as long as the metric is a Schwarzschild one due to the gravitational pull of the Sun, the accuracy of the null geodesic integration is not the required one. At the same time it is obvious that the Gaia accuracy can be attained when the observing direction is sufficiently far from the gravitational perturbation of the Solar System bodies. One can approximately estimate these “avoidance zones” by comparing the estimated measurement accuracy with the Schwarzschild contribution to the light deflection of each single body as a function of its angular distance ψ from the source (Misner et al. 1973)

$$\begin{array}{c}\hfill \delta \psi =\frac{(1+\gamma )M}{c^2{r}_{\mathrm{s}}}\sqrt{\frac{1+cos\psi }{1-cos\psi }},\end{array}$$(15)

where M is the mass of the body and r_s its distance from the satellite. For example, the minimum angular distance that a source must have from Jupiter in order to keep the model accuracy below the ∼10 μas, level is about 10.5°.

On the other hand, the model accuracy can be enhanced by adding the estimated light deflection effects of Eq. (15) to the pure Schwarzschild model; such contribution to the light deflection has to be computed at an appropriate retarded time which takes into account the finite speed of the light. These corrections are done separately for ϕ and ζ by projecting the δψ on the AL and AC direction according to the geometrical configuration of the observation. Clearly this is in no way an exact solution, but, as shown in (Vecchiato et al., in prep.) where a more detailed description of the model will also be provided, numerical tests have proven that the accuracy of this model is much better than that of the pure Schwarzschild one.

2.5. Satellite barycentric motion and attitude model

As shown in Crosta & Vecchiato (2010), the four-velocity of the satellite can be written as $u_{\mathrm{s}}^{\alpha }=u_{\mathrm{s}}^0({\delta }_0^{\alpha }+{\nu }_{\mathrm{s}}^i{\delta }_i^{\alpha })\equiv u_{\mathrm{s}}^0\{1,{\nu }_{\mathrm{s}}^x,{\nu }_{\mathrm{s}}^{\mathit{y}},{\nu }_{\mathrm{s}}^z\}$ where ν_s is the coordinate velocity of the satellite and, using the IAU resolutions (Kaplan 2005),

$$\begin{array}{c}\hfill u_{\mathrm{s}}^0=1+\frac{U}{c^2}+\frac{1}{2}\frac{{\nu }_{\mathrm{s}}^2}{c^2}+\mathcal{O}\left(\frac{{\nu }^4}{c^4}\right)·\end{array}$$(16)

Here ${\nu }_{\mathrm{s}}^2={\left({\nu }_{\mathrm{s}}^x\right)}^2+{\left({\nu }_{\mathrm{s}}^y\right)}^2+{\left({\nu }_{\mathrm{s}}^z\right)}^2$ and U is the (BCRS) N-body gravitational potential of the Solar System at the Gaia location, which implies that the four-velocity, contrary to the null geodesic, is expressed to the right accuracy.

The construction of the attitude tetrad, as anticipated in Eq. (9), starts from this four-vector by setting ${\mathit{e}}_{\widehat{0}}^{\alpha }=u_{\mathrm{s}}^{\alpha }$. Then one has to fix the spatial triad by considering the following transformations (Crosta & Vecchiato 2010; Bianchi et al. 2011):

1. the origin of the tetrad is located at the BCRS coordinates of the satellite barycentre, which induces a first transformation to the so-called local BCRS tetrad (Bini et al. 2003). This is a tetrad at rest with respect to the BCRS, but whose spatial axes have a general relativistic contribution caused by the gravitational potential of the Solar System at that point;

2. the satellite barycentre moves with four-velocity $u_{\mathrm{s}}^{\alpha }$ in the BCRS, which induces a second transformation to the boosted tetrad with spatial axes e_â,b whose origin is comoving with the barycentre of the satellite³. If the four-velocity of the BCRS, and therefore also that of the local BCRS, is u^α, then at each instant the comoving tetrad can be computed with a special relativistic transformation, namely a boost, on the previous one whose Lorentz factor is ${\gamma }_{\mathrm{L}}={(1-{\nu }_{\mathrm{s}}^2/c^2)}^{-1/2}=-u_{\mathrm{s}}^{\alpha }{u}_{\alpha }$;

3. the satellite orientation, that is the SRS e_â, can be finally obtained by a Euclidean transformation of its spatial axes, ${\mathit{e}}_{\widehat{a}}^{\alpha }=R{e}_{\widehat{a},b}^{\alpha }$, where R is simply a 3D rotation matrix obtained by a standard attitude parametrization. AGIS uses quaternions, while GSR is based on the Modified Rodrigues Parameters (MRP) σ₁, σ₂ and σ₃ described, for example, in Schaub & Junkins (1996). Since MRP representation uses only three parameters instead of the four of the quaternions, adopting the MRP representation allows to reduce the number of attitude unknowns. Moreover, quaternions require the enforcing of a normalization condition, implemented through additional constraint equations, due to their redundant parametrization. On the other hand, such redundancy allows to avoid the problem of singularities that is present in other parametrizations (including the MRP). The scanning law of Gaia, however, is always sufficiently far from the singularity points of the MRP.

Despite their different accuracy, the null geodesic and the attitude parts stem from similar mathematical assumptions, which makes them mutually compatible at the (ν/c)² order. Basically, the main requirement is that the metric has the form

$$\begin{array}{c}\hfill g_{\alpha \beta }={\eta }_{\alpha \beta }+h_{\alpha \beta }+\mathcal{O}\left(h^2\right)\end{array}$$(17)

where, for the null geodesic integration, we retained only the (ν/c)²-order perturbation term given by the gravitational potential of a massive body $h_{00}^{(2)}=2U/c^2$, so that $h_{00}\sim h_{\mathit{ij}}=h_{00}^{(2)}{\delta }_{\mathit{ij}}+\mathcal{O}({\nu }^4/c^4)$, while the attitude model can include also higher order terms, namely $h_{0i}\sim {\nu }^3/c^3+\mathcal{O}({\nu }^5/c^5)$ which enter in the definition of the local BCRS tetrad.

Considering the problem of the attitude representation, it has to be recalled that, regardless of the specific parametrization chosen to represent the attitude matrix, each component of such parametrization has to be expressed with a finite and time-dependent number of attitude unknowns x^A. The required final accuracy of the measurements does not permit the use of a physical model for the attitude due to several factors; for example, the micro-propulsion system introduces a high-frequency noise. As in AGIS, then, an attitude parameter (that is a quaternion or an MRP component) is represented by a purely numerical expansion written as a linear combination of time-dependent polynomial functions B(t) of degree M − 1 called B-Splines (Ahlberg & Nilson 1967) whose characteristics are summarized in the following.

The expansion of the jth component of a generic representation, say S, reads

$$\begin{array}{c}\hfill S_j\left(t\right)=\sum _{n=0}^{N-1}{c}_n^{\left(j\right)}{B}_n\left(t\right),\end{array}$$(18)

where $c_n^{\left(j\right)}$ are unknown attitude coefficients to be determined. The function S_j(t) is defined in a time interval [t_beg,…,t_end] divided in K > 0 sub-intervals identified by a sequence of instants {τ_k}, k = 0, …, K called nodes, which constitutes the so-called support of the series. The B-Splines are useful in our case because their support is minimal, that is B_n(t) ≠ 0 only in M sub-intervals, so if τ_n < t < τ_{n + 1} then,

$$\begin{array}{c}\hfill S_j\left(t\right)=f(c_{n-M/2}^{\left(j\right)},\dots ,c_{n-2}^{\left(j\right)},c_{n-1}^{\left(j\right)},c_n^{\left(j\right)},c_{n+1}^{\left(j\right)},\dots ,c_{n+M/2-1}^{\left(j\right)}),\end{array}$$(19)

in the sense that at each time the expansion depends only on M unknowns. For Gaia the expansion is in cubic B-splines, that is M = 4, and therefore

$$\begin{array}{c}\hfill S_j\left(t\right)=f(c_{n-2}^{\left(j\right)},c_{n-1}^{\left(j\right)},c_n^{\left(j\right)},c_{n+1}^{\left(j\right)}).\end{array}$$(20)

Finally, having split a time segment in K − 1 intervals, the resulting number of degrees of freedom, and thus of attitude unknowns, for that segment is given by N = K + M − 1.

Another important point in using a series expansion like the B-Splines lies in the possibility of exploiting the linearization of the observation equations to implement the so-called differential attitude approach. As for the integration of the geodesics, this will be formulated in more detail in a forthcoming paper, but since the motivation and the principle of this approach is directly linked to the solution method, Sect. 3 will present a line-of-principle description from the mathematical point of view.

2.6. Instrument parameters

As anticipated in the general overview of Sect. 2, except for the BAM, Gaia has no onboard metrologic instrument, and a set of long-term calibration parameters x^C has to be introduced among the unknowns of the global sphere reconstruction. The instrument calibration parameters currently used by GSR are those described in Lindegren et al. (2012). In this model, like for the attitude, the calibration is done by a purely numerical description of the perturbations of the two field angles η and ζ on the focal plane, respectively in the along- and across-scan direction. The parameters of this model can be divided in two classes: geometric and spectrophotometric. In the former the deviations dη and dζ can be described by geometrical deviations at the level of a single CCD of the astrometric focal plane from the nominal configurations, namely shifts, shears plus rotations and distortions, while in the latter magnitude- and spectrum-dependent shifts dη are introduced.

In practice the AL measurement ϕ in Eq. (4) is defined as

$$\begin{array}{c}\hfill \varphi =\eta \pm \frac{\mathrm{\Gamma }}{2},\end{array}$$(21)

where Γ is the Basic Angle value (BA) and the sign is positive when the observation is on the preceeding Field of View (FoV2) and negative for the following Field of View (FoV1). The calibration model then puts

$$\begin{array}{c}\hfill \eta ={\eta }^0+\sum _{r=0}^2\mathrm{\Delta }{\eta }_{\mathit{rfnj}}{L}_r^{\ast }\left(\tilde{\mu }\right)+\delta {\eta }_{\mathit{nm}}+C_n^{\mathrm{mag}}(G-G_{\mathrm{ref}})+C_{\mathit{fn}}^{\mathrm{sp}}(\nu -{\nu }_{\mathrm{ref}})\end{array}$$(22)

in which η⁰ is the nominal AL field angle and Δη_rfnj, δη_nm, $C_n^{\text{mag}}\hspace{0.17em}\text{and}\hspace{0.17em}{C}_{fn}^{\text{sp}}$ are the above mentioned calibration unknowns, namely the AL large and small scale geometric corrections, and the magnitude- and spectrum-dependent shifts. The indices indicate the dependencies of these parameters with respect to the instrument configuration.

The large-scale AL parameters Δη_rfnj depend on the FoV (f), the CCD index (n) and have a temporal variation so that there exists a different set of parameters extending over a certain time interval, indexed by j. Moreover, the shift, shear plus rotation, and distortion contributions are modeled by different orders 0 ≤ r ≤ 2 of the Legendre polynomials $L_r^{\ast }\left(\tilde{\mu }\right)$ computed at the normalized AC pixel coordinate

$$\begin{array}{c}\hfill \tilde{\mu }=\frac{\mu -{\mu }_0+0.5}{1966},\end{array}$$(23)

where 14 ≤ μ ≤ 1979 is the pixel coordinate of the measurement and μ₀ = 14 is the pixel coordinate of the beginning of the light sensitive area of the CCD.

The small scale AL parameters δη_nm depend on the CCD index n and on the AC pixel index m, which goes from 1 + μ₀ to pixcols + μ₀, where pixcols is the total number of illuminated pixels. The same indexes are used to identify the dependencies of the magnitude- and spectrum-dependent shifts, which need a reference magnitude and frequency G_ref and ν_ref res-pectively.

Similarly, the AC measurement ζ of Eq. (5) is

$$\begin{array}{c}\hfill \zeta ={\zeta }^0+\sum _{r=0}^2\mathrm{\Delta }{\zeta }_{\mathit{rfnj}}{L}_r^{\ast }\left(\tilde{\mu }\right)+\delta {\zeta }_{\mathit{nm}}\end{array}$$(24)

where it is made evident that in the AC direction only the geometric parameters are taken into account. This calibration model introduces a degeneracy among the geometric parameters and the attitude which has to be removed by means of a set of appropriate constraint equations.

3. Reconstructing the global astrometric sphere

3.1. The linearized system of equations

The global astrometric sphere reconstruction, in the case of Gaia, requires the solution of a large, sparse and overdetermined system of linearized equations in the least-squares sense. Each observation in fact is represented by a formula either as Eq. (1) or (2), which ultimately are functions of four types of unknowns, namely x^S, x^A, x^C, x^G, and together these observations produce a system of equations whose dimension depend on the number of unknowns n_unk and of observations n_obs.

In the case of Gaia the total number of unknowns is n_unk = n^S + n^A + n^C + n^G, where the number of unknowns for the sources is n^S = 5n^*, n^* ∼ 10⁸ is the number of primary stars and n^A, n^C and n^G are the total number of attitude, instrument and global unknowns respectively. Since generally n^G ≪ n^C ≪ n^A ≪ n^*, it is reasonable to put n_unk ≲ 6n^*. Moreover, we know that on average each star is observed about 10³ times during the mission lifetime, so the total number of observations is n_obs ∼ 10¹¹.

Since n_obs ≫ n_unk the system is not only large but also overdetermined, so in principle it can be solved in the least squares sense, providing not only an estimation of the unknowns but also of their errors and correlations. Furthermore, each equation represents the observation of a single source in a specific time interval, therefore it depends only on (up to) 5 source unknowns and a comparably limited number of attitude and instrument parameters, which implies the sparseness of the system.

Finally, both Eqs. (1) and (2) are highly non-linear in their unknowns, which would make the solution of the system based on a maximum likelihood approach numerically intractable. However an approximate set of values x̄ can be given for all the parameters, thus the jth AL observation equation can be linearized around these values

$$\begin{array}{cc}\hfill -sin{\varphi }_j^{\mathrm{c}}\delta {\varphi }_j\simeq & \sum _{i=1}^{n^{\mathrm{S}}}{\frac{\partial f_{\varphi }\left(\mathbf{x}\right)}{\partial x_i^{\mathrm{S}}}|}_{\overline{\mathbf{x}}}\delta x_i^{\mathrm{S}}+\sum _{i=1}^{n^{\mathrm{A}}}{\frac{\partial f_{\varphi }\left(\mathbf{x}\right)}{\partial x_i^{\mathrm{A}}}|}_{\overline{\mathbf{x}}}\delta x_i^{\mathrm{A}}+\hfill \\ \hfill & \sum _{i=1}^{n^{\mathrm{C}}}{\frac{\partial f_{\varphi }\left(\mathbf{x}\right)}{\partial x_i^{\mathrm{C}}}|}_{\overline{\mathbf{x}}}\delta x_i^{\mathrm{C}}+\sum _{i=1}^{n^{\mathrm{G}}}{\frac{\partial f_{\varphi }\left(\mathbf{x}\right)}{\partial x_i^{\mathrm{G}}}|}_{\overline{\mathbf{x}}}\delta x_i^{\mathrm{G}}\hfill \end{array}$$(25)

where ${\varphi }_j^{\mathrm{c}}=arccos\left(f_{\varphi }\left(\overline{\mathbf{x}}\right)\right)$, $\delta {\varphi }_j={\varphi }_j^{\mathrm{o}}-{\varphi }_j^{\mathrm{c}}$, $\delta x_i=x_i^{\mathrm{true}}-{\overline{x}}_i$, ${\varphi }_j^{\mathrm{o}}$ is the abscissa of that observation measured by Gaia and x^true is the set of unknown true values. A similar expansion can be written for the across scan measurements, and the problem can be thus reduced to the solution of a linear system of observation equations

$$\begin{array}{c}\hfill \mathbf{b}=A\delta \mathbf{x}\end{array}$$(26)

where b = {b_j}^T, j = 1, …, n_obs, is the vector of the known terms, δ x is the unknown vector and A is the n_obs × n_unk design matrix of the system whose coefficients are a_ji = (∂f/#x2202;_xi) (x̄). In the coefficients, obviously, f() = f_ϕ() or f() = f_ζ() if the observation is AL or AC respectively. The corresponding generic known term b_j is equal to $-sin{\varphi }_j^{\mathrm{c}}\delta {\varphi }_j$ or to $cos{\zeta }_j^{\mathrm{c}}\delta {\zeta }_j$, and n_obs = n_ALobs + n_ACobs, for the number of AL and AC observations n_ALobs and n_ACobs can in general be different⁴.

The resulting system of equations is intrinsically rank-deficient when astrometric and attitude parameters have to be solved at the same time. This well-known issue comes from the invariance of the solution for a rigid three-components rotation and three-components spin of the reference system, which in general requires the introduction of 6 constraint equations (see, e.g., discussions in de Felice et al. 1998, 2001; Berghea et al. 2016 about the actual need of such a constraint and their possible implementations.)

In order to take into account measurement errors, which are different for each observation, a n_obs × n_obs weight matrix⁵

$$\begin{array}{c}\hfill W=\mathrm{diag}\left(\frac{{ϵ}_0}{{ϵ}_j}\right)\end{array}$$(27)

is introduced, where ϵ_j is the estimated standard deviation of the jth observation and ϵ₀ is a reference value, currently assumed to the estimated uncertainty at magnitude G = 21.5. The system thus becomes

$$\begin{array}{c}\hfill W\mathbf{b}=(WA)\delta \mathbf{x}.\end{array}$$(28)

The solution in the least-squares sense of such system provides δx̂, namely the best-fit estimation of δ x, that is used to update the catalog values x̄ to the improved estimation of the true values x̃, that is

$$\begin{array}{c}\hfill {\mathbf{x}}^{\mathrm{true}}\simeq \tilde{\mathbf{x}}=\overline{\mathbf{x}}+\delta \widehat{\mathbf{x}}.\end{array}$$(29)

The formal least-squares solution of the system of Eq. (28) is δ x = [(WA)^T(WA)]⁻¹(WA)^TW b, where [(WA)^T(WA)]⁻¹ is also the covariance matrix providing the estimation of the variances and covariances of the unknowns.

As mentioned in Sect. 2.4, in this version of the pipeline the relativistic light deflection effect of Solar System objects different from the Sun is taken into account to improve on the model accuracy. However, it has to be stressed that this contribution is added only to the known term as a correction to ${\varphi }_j^{\mathrm{c}}$ (and to ${\zeta }_j^{\mathrm{c}}$)⁶, and it is not considered in the derivatives. This is equivalent to a correction to the observations, which are thus forced to fit more closely the pure Schwarzschild model. The numerical consequence of this approach will be explained in Sect. 5.

3.2. Solution of the linearized system of equations

It is known (see e.g. Bombrun et al. 2010 and references therein) that the computational complexity of the inversion of an N × N matrix with direct methods is ∝N³. Since the normal matrix A^TA is n_unk × n_unk with n_unk ∼ 6 × 10⁸, we can immediately realize that in the Gaia case this problem can be addressed only resorting to iterative algorithms. A brute-force approach like this, in fact, would require about 10²⁶ FLOPs a requirement which cannot be reduced to an acceptable level even taking into account the sparsity of A.

AGIS uses a block-iterative technique to solve the equation system (Lindegren et al. 2012) where, in practice, the solution is obtained by solving separately each block of unknowns (x^S, x^A, x^C and x^G) and iterating the process until convergence (Fig. 5, left). This allows an implementation as an embarassingly parallel algorithm, which is strictly needed in the pure Java environment chosen in this case. On the other hand, the computation of standard uncertainties, and covariances, of the astrometric parameters requires the inversion of the entire normal matrix, which is needed to take properly into account the attitude-induced correlations among different sources; since this task is not feasible, one must rely on some approximate covariance model to estimate such statistical parameters (Holl & Lindegren 2012; Holl et al. 2012).

Fig. 5. Schematic representation of the block-iterative procedure used by AGIS (left) and of the fully iterative one used by GSR (right). See the text for explanation.

Also lead by the initial consideration on the criticality of having an independent solution algorithm it was decided to use a fully iterative method (Fig. 5, right) to guarantee a fully general convergence mechanism of the complete system, and to allow a mixed Java/C-C++ coding of the pipeline in order to make possible the implementation of the needed parallel algorithm. Such implementation, in fact, is realized by means of a hybrid MPI/OpenMP parallel solver which runs at the CINECA supercomputing facilities (Bandieramonte et al. 2012; Becciani et al. 2014). GSR therefore implements a customized version of PC-LSQR (Baur et al. 2008) a conjugate gradient-based algorithm, originally proposed by Paige & Saunders (1982). As for any iterative algorithm, this is equivalent to computing, at each iteration (i), an approximate solution

$$\begin{array}{c}\hfill \delta {\mathbf{x}}^{\left(i\right)}={(A^{\mathrm{T}}A)}^{-1}{A}^{\mathrm{T}}{\mathbf{b}}^{(i-1)}\end{array}$$(30)

and then evaluating the vector of the residuals

$$\begin{array}{c}\hfill {\mathbf{r}}^{\left(i\right)}=\mathbf{b}-A\delta {\mathbf{x}}^{\left(i\right)}\end{array}$$(31)

which has to be minimized in the least-squares sense, according to suitable convergence conditions defined by the algorithm itself. Among the possible stopping conditions we have:

• the vector of the residuals has a 2-norm lower than a threshold value; the LSQR algorithm generates a series of residual vectors whose norm decreases monotonically; in the case of a compatible system the series goes to zero, while for non-compatible systems it converges to a positive finite limit;

• the 2-norm of A^T r⁽ⁱ⁾ (namely the norm of the residuals of the normal system A^T b = A^TAδ x) is lower than a threshold value; this is the condition used to guarantee the convergence of non-compatible systems to a solution in the least-squares sense; it is worth stressing here that the equation system to be solved for the global astrometric sphere reconstruction is non-compatible⁷, and indeed its solution is obtained in the least-squares sense, therefore this is the stopping condition that has to be reached in our case;

• the iterative estimation of the condition number of the matrix exceeds a given upper threshold;

• a fixed maximum number of iterations is reached.

The choice of the LSQR algorithm is also motivated by the possibility of further enhancing its standard implementation. Indeed, in its original definition it provides an estimation only of the diagonal elements of the covariance matrix (A^TA)⁻¹, namely the variances σ_{δ x}, but upgrades of this algorithm (see, e.g. Guo 2008; Kostina & Kostyukova 2012) would allow to estimate also any selected group of its off-diagonal elements, namely the covariances.

The GSR version of the algorithm uses a preconditioning technique, which basically consists in a renormalization of the columns of A, made to improve the convergence speed of the system. This version, moreover, is tailored on the sphere solution problem in particular for what concerns the parallelization algorithm and the memory occupancy. The latter has been optimized with respect to the classic Yale Sparse Matrix Format (Buluç et al. 2009) exploiting the almost constant number of non-zero coefficients to eliminate one of the n_obs pointer vector. Regarding the former, instead, the design matrix is built by sorting the observation by source number. In this way it is possible to distribute the stars in almost independent subsets on each Processing Element (PE), in fact the matrix-vector product, which is the core of the LSQR algorithm, can be computed by distributing on each PE an equal number n_obs/n_PE of the rows of A and of the vector b, and a number 5n^*/n_PE of astrometric unknowns. The other unknowns are then duplicated on all the PEs making the product almost communication-free (Bandieramonte et al. 2012).

3.3. Differential attitude and AC observations

Let us say that one has to compute the attitude coefficients of a generic linearized observation equation. From the above considerations we can formally write

$$\begin{array}{c}\hfill {\mathbf{x}}^A=\{{\sigma }_1(c_0^{\left(1\right)},\dots ,c_{N-1}^{\left(1\right)}),{\sigma }_2(c_0^{\left(2\right)},\dots ,c_{N-1}^{\left(2\right)}),{\sigma }_3(c_0^{\left(3\right)},\dots ,c_{N-1}^{\left(3\right)})\},\end{array}$$(32)

so that in general the catalog attitude is a set of initial values of the B-Spline expansion coefficients $\left\{{\overline{c}}_i^{\left(j\right)}\right\}$. Equation (25) can thus be written as

$$\begin{array}{cc}\hfill -sin{\varphi }_{\mathrm{calc}}\delta \varphi \simeq & \sum _{i=1}^{n^{\mathrm{S}}}{\frac{\partial f_{\varphi }\left(\mathbf{x}\right)}{\partial x_i^{\mathrm{S}}}|}_{{\mathbf{x}}_{\mathbf{0}}^{\mathrm{S}}}\delta x_i^{\mathrm{S}}+\sum _{j=1}^3{\frac{\partial f_{\varphi }\left(\mathbf{x}\right)}{\partial {\sigma }_j}|}_{{\overline{\sigma }}_j}\delta {\sigma }_j+\hfill \\ \hfill & \sum _{i=1}^{n^{\mathrm{C}}}{\frac{\partial f_{\varphi }\left(\mathbf{x}\right)}{\partial x_i^{\mathrm{C}}}|}_{\overline{\mathbf{x}}}\delta x_i^{\mathrm{C}}+\sum _{i=1}^{n^{\mathrm{G}}}{\frac{\partial f_{\varphi }\left(\mathbf{x}\right)}{\partial x_i^{\mathrm{G}}}|}_{\overline{\mathbf{x}}}\delta x_i^{\mathrm{G}}\hfill \\ \hfill =& \sum _{i=1}^{n^{\mathrm{S}}}{\frac{\partial f_{\varphi }\left(\mathbf{x}\right)}{\partial x_i^{\mathrm{S}}}|}_{{\mathbf{x}}_{\mathbf{0}}^{\mathrm{S}}}\delta x_i^{\mathrm{S}}+\sum _{j=1}^3\sum _{i=0}^{N-1}{\frac{\partial f_{\varphi }\left(\mathbf{x}\right)}{\partial {\sigma }_j}\frac{\partial {\sigma }_j}{\partial c_i^{\left(j\right)}}|}_{{\overline{c}}_i^{\left(j\right)}}\delta c_i^{\left(j\right)}\hfill \\ \hfill & \sum _{i=1}^{n^{\mathrm{C}}}{\frac{\partial f_{\varphi }\left(\mathbf{x}\right)}{\partial x_i^{\mathrm{C}}}|}_{\overline{\mathbf{x}}}\delta x_i^{\mathrm{C}}+\sum _{i=1}^{n^{\mathrm{G}}}{\frac{\partial f_{\varphi }\left(\mathbf{x}\right)}{\partial x_i^{\mathrm{G}}}|}_{\overline{\mathbf{x}}}\delta x_i^{\mathrm{G}}\hfill \\ \hfill =& \sum _{i=1}^{n^{\mathrm{S}}}{\frac{\partial f_{\varphi }\left(\mathbf{x}\right)}{\partial x_i^{\mathrm{S}}}|}_{{\mathbf{x}}_{\mathbf{0}}^{\mathrm{S}}}\delta x_i^{\mathrm{S}}+\sum _{j=1}^3\sum _{i=0}^{N-1}{\frac{\partial f_{\varphi }\left(\mathbf{x}\right)}{\partial {\sigma }_j}|}_{{\overline{\sigma }}_j}{B}_i\left(t\right)\delta c_i^{\left(j\right)}\hfill \\ \hfill & \sum _{i=1}^{n^{\mathrm{C}}}{\frac{\partial f_{\varphi }\left(\mathbf{x}\right)}{\partial x_i^{\mathrm{C}}}|}_{\overline{\mathbf{x}}}\delta x_i^{\mathrm{C}}+\sum _{i=1}^{n^{\mathrm{G}}}{\frac{\partial f_{\varphi }\left(\mathbf{x}\right)}{\partial x_i^{\mathrm{G}}}|}_{\overline{\mathbf{x}}}\delta x_i^{\mathrm{G}}\hfill \end{array}$$(33)

in which we have exploited the fact that

$$\begin{array}{c}\hfill \mathrm{d}{\sigma }_j\left(t\right)=\sum _{i=0}^{N-1}\frac{\partial {\sigma }_j}{\partial c_i^{\left(j\right)}}\mathrm{d}{c}_i^{\left(j\right)}\end{array}$$(34)

and that, from Eq. (18),

$$\begin{array}{c}\hfill \frac{\partial {\sigma }_j}{\partial c_i^{\left(j\right)}}=B_i\left(t\right).\end{array}$$(35)

In this way the updated attitude at a generic time t is

$$\begin{array}{c}\hfill {\sigma }_j(t)=\sum _{i=0}^{N-1}({\overline{c}}_i^{\left(j\right)}+\delta c_i^{(j)})B_i\left(t\right).\end{array}$$(36)

It should be observed, however, that in Eq. (33) ∂f_ϕ(x)/∂σ_j does not depend on the representation of dσ_j: it just needs a (catalog) value for the attitude parameter ${\overline{\sigma }}_j\equiv {\overline{\sigma }}_j\left(t\right)$ at the observation time t. We can therefore decide to compute this σ_j using an expansion which, in principle, may have nothing to do with the one used to represent dσ_j since our only necessity is to have a way to evaluate ${\overline{\sigma }}_j\left(t\right)$. In formulae, we could write at any time ${\sigma }_j(t)={\overline{\sigma }}_j\left(t\right)+\mathrm{d}{\sigma }_j\left(t\right)$ where

$$\begin{array}{cc}\hfill {\overline{\sigma }}_j\left(t\right)& =\sum _{m=0}^{N_{\mathrm{C}}-1}{c}_{m,\mathrm{C}}^{\left(j\right)}{B}_m\left(t\right)\hfill \end{array}$$(37)

$$\begin{array}{cc}\hfill \mathrm{d}{\sigma }_j\left(t\right)& =\sum _{i=0}^{N_{\mathrm{U}}-1}{c}_{i,\mathrm{U}}^{\left(j\right)}{B}_i\left(t\right)\hfill \end{array}$$(38)

are the catalog attitude and its update, exactly as in our original scenario, but expanded over two different supports {τ_{m, C}}, m = 0, …, N_C and {τ_{i, U}}, i = 0, …, N_U. Moreover, since the catalog attitude has to represent the starting point for the unknowns, we should have initially dσ_j(t) = 0, that is $c_{i,\mathrm{U}}^{\left(j\right)}=0$, so that in the end $c_{i,\mathrm{U}}^{\left(j\right)}=\delta c_{i,\mathrm{U}}^{\left(j\right)}$, and therefore the updated attitude is computed at any time t simply by summing the catalog and the differential attitude dσ at the same time:

$$\begin{array}{cc}\hfill {\sigma }_j(t)& ={\overline{\sigma }}_j\left(t\right)+\mathrm{d}{\sigma }_j\left(t\right)\hfill \end{array}$$(39)

$$\begin{array}{cc}& =\sum _{m=0}^{N_{\mathrm{C}}-1}{c}_{m,\mathrm{C}}^{\left(j\right)}{B}_m\left(t\right)+\sum _{i=0}^{N_{\mathrm{U}}-1}\delta c_{i,\mathrm{U}}^{\left(j\right)}{B}_i\left(t\right).\hfill \end{array}$$(40)

This facilitates the choice of the knot support in the computation of the coefficients and makes the calculation of the latter much faster. The placement of the knots, in fact, has to guarantee that each attitude parameter can be solved, a condition that can be met basically by allowing a sufficient number of observations between each pair of consecutive knots. There are, however, several reasons which could require to drop some observations after the computation of the matrix coefficients, thus inducing a rearrangement of the knot sequence; in the non-differential approach, such rearrangement would entail a recalculation of all the coefficients, while in the differential approach the latter must be computed only once, with the exception of the B_i(t) polynomials for the differential knot sequence.

One last point of the Gaia attitude reconstruction has to be stressed. Given the three coefficients $C_k^{\left(j\right)}$, j = 1, 2, 3 relative to the k-th unknown of each Rodrigues parameter σ_j ⁸, their ratios $R_k^{(i,j)}=C_k^{\left(i\right)}/C_k^{\left(j\right)}$ do not depend from the B-spline base function B_k(t). As a consequence, whatever the time dependence of two coefficients, one can always find a small time separation δt = t₂ − t₁ between two observations at t₁ and t₂ respectively which is short enough to keep $R_k^{(i,j)}$ almost constant between t₁ and t₂. As pointed out above, this “constancy” in principle does not depend on the explicit formulation of the coefficients, since it is sufficient to have time-dependent expressions repeating themselves along all the observations, however it is easy to understand that the slower the time variation, the longer this time interval can be kept. Therefore the independence of the ratios from B_k(t) contributes to increase the amount of time over which the $R_k^{(i,j)}$ remain almost the same among the different observations or, equivalently, over which the proportionality ${\mathbf{C}}_k^{\left(1\right)}\propto {\mathbf{C}}_k^{\left(2\right)}\propto {\mathbf{C}}_k^{\left(3\right)}$ holds approximately for each index k.

Since each of these values is represented into a different column of the system of equations, this quasi-proportionality produces an extremely ill-conditioned design matrix. Such an ill-conditioning problem is solved by introducing the AC observations. This problem will be discussed in more detail in a forthcoming publication (Vecchiato et al., in prep.).

4. The AVU/GSR pipeline

4.1. The overall pipeline schema

As mentioned in Sect. 1 the algorithms described so far have been implemented in the GSR pipeline, which operates in the context of the Astrometric Verification Unit of the DPAC. A detailed description of such pipeline is out of the scope of this paper, and more details can be found in Vecchiato et al. (2012) and in a forthcoming publication (Messineo et al., in prep.). Here we will thus give a synthetic summary of its main characteristics.

The GSR pipeline is composed of several scientific modules, called in sequence by the infrastructure software, which provides the overall workflow and the DB interaction functionalities (Fig. 6). In essence, GSR gets the AGIS solution and its corresponding input data. The latter are used to produce an independent sphere solution which is then compared to that of AGIS. The results of the comparison, as outcomes of statistical tests and graphics, are sent to the Gaia Main DataBase (MDB). Alerts are also foreseen if the comparison bears evidence of statistically significant differences between the two solutions.

Fig. 6. Schematic representation of the AVU/GSR pipeline.

All the scientific modules can be grouped in three main types:

• The “Pre-solver” modules, which implements the computation of the system coefficients and known terms according to the guidelines sketched in the Sects. 2 and 3.

• The “Solver” module, which computes the system solutions as described in Sect. 3.

• The “Post-solver” modules, performs all the successive operations from the catalog update to the comparison between the AGIS and GSR solutions summarized in Sect. 4.2.

As anticipated, all the modules are written in Java and run at the DPCT-ALTEC, with the exception of the Solver module which is written in C/C++ and runs at the DPCT-CINECA. The latter simply receives the coefficients and known terms for the system from ALTEC and sends back there the system solution. The pipeline concludes with the production of a report which contains all the relevant information needed to give a first evaluation of the results.

4.2. Post-Solver GSR pipeline modules

As mentioned above, GSR has to compare its results with the AGIS solution and to provide an evaluation of the differences between the two spheres. In particular, it is foreseen that possible differences larger than the expected accuracy of Gaia are investigated to assess the scientific reliability of the result.

As pointed out in Sect. 3 the sphere reconstruction problem is intrinsically rank-deficient, and as a minimum it can be solved except for a six-parameters transformation representing a rigid rotation and a spin difference of the reference system. In the equation system solved by GSR, these can be directly included as additional constraint equations, equivalent to the choice of a specific (and arbitrary) reference system which, in general, is different from that of the AGIS solution. It is therefore necessary to de-rotate one solution, that is to bring both catalogs into a common reference system, before attempting a comparison between them. Moreover, both the de-rotation and comparison procedures are done on catalogs referring to the purely spatial hypersurface of the same observer, therefore there is no need to resort to General Relativity here. The algorithms are least-squares reconstructions of purely Euclidean transformations and statistical analysis of coordinate differences.

GSR can use standard statistical algorithms, like χ² tests and Kolmogorov-Smirnov, to analyze the differences between the AGIS and GSR global astrometric sphere solutions. In addition to these, a powerful method in this context is provided by the use of a decomposition in Scalar and Vector Spherical Harmonics (VSH, see e.g. Hill 1954; Arfken & Weber 2012) as base functions to model the vector field of GSR/AGIS residuals as a series expansion, up to a suitable degree, whose coefficients can be either estimated by a standard least-squares fit or computed analytically, whenever appropriate. Once the significance of the various coefficients associated to each VSH degree is positively tested, the latter quantify the presence of a systematic error in the residuals at scale lengths of the order of 180/l deg, where l is the corresponding VSH degree, therefore their nature must be addressed.

This last algorithm is the most computationally-intensive, so its implementation is the only one besides the system solution that can require a non-embarrassingly parallel coding (to be executed at the DPCT-CINECA) in case a decomposition of order l_max ≳ 100 is required. The details of the algorithms can be found in Bucciarelli et al. (2011a,b) and in a forthcoming paper (Bucciarelli et al., in prep.), which also contains a characterization of the VSH accuracy for the Gaia case.

5. Demonstration runs on simulated data

5.1. General considerations

Similarly to what was done for AGIS (Lindegren et al. 2012), the accuracy of the AVU/GSR pipeline has to be verified before this software system enters operations. The general agreement was that GSR had to perform a “Demonstration Run test”, showing its ability to reproduce the AGIS results at a comparable accuracy level on simulated data similar to those used in the cited paper. For that purpose, the GSR demonstration run was split into two different tests:

• Test 1:

  solve the GSR astrometric sphere with perturbed Source (S), Attitude (A) and Calibration (C) parameters and noise-free observations. The aim is to give a precise assessment of the numerical accuracy of the AVU/GSR pipeline and of the GSR2 astrometric model.

• Test 2:

  solve the GSR astrometric sphere with perturbed S, A and C parameters and noisy observations, similarly to what was done for AGIS (see the cited paper above). The goal is reproducing the AGIS results at the appropriate accuracy.

Following the methodology adopted in the AGIS paper, the perturbations of the source and attitude parameters are Gaussian, while the system starts from true values for the instrument (calibration parameters). Moreover, a long-period modulation of the BA is injected in the known terms (observations) and reconstructed by means of the L₀ large-scale AL Instrument parameters (Δη) at the sub-μas level. The instrument model is used only in this scope.

As shown in Sect. 3, the global astrometric sphere reconstruction, in the case of Gaia, requires the solution in the least-squares sense of a large, sparse and overdetermined system of linearized equations. In GSR the resulting linear equation system is solved in full by utilizing a parallelized implementation of the LSQR algorithm. Provided that the starting values are sufficiently close to the true ones the system converges to its best possible least-squares solution, a statement which has to be intended in comparison with the block-iterative algorithm adopted by AGIS.

This simple framework is complicated by the discrepancy between the astrometric accuracy of the two relativistic models implemented in AGIS and the current GSR. AGIS implements the GREM astrometric model, which is accurate to the (ν/c)³ order needed to match the Gaia mission final accuracy, and the simulated data are generated with the same astrometric model. The current GSR, instead, uses a model of the RAMOD family whose (ν/c)² PPN-Schwarzschild metric, as explained in Sect. 2.4, is improved in accuracy by an approximate estimation of the contributions of the planets of the solar system and of the Moon, which are subtracted from the known terms.

This technique, however, cannot model at the required level of accuracy the observation equation of this intrinsically non-Schwarzschild problem. The accuracy of the GSR2 astrometric model, therefore, does not match that of the simulated data everywhere on the sky, and leaves unmodeled some effects at the order of the first derivative of the observable, that enter directly in the sphere reconstruction problem. Furthermore, and more importantly, it does not cope with the final accuracy expected for Gaia, especially at the bright end of the magnitude range. This approximate approach could be further improved, but at the price of a considerable mathematical complication, an effort that in any case would not be able to ensure a perfect match in terms of model accuracy.

On the other hand, it is still possible to reach (almost) the same results of the AGIS demonstration run by minimizing the impact of these modeling issues by, firstly, dropping observations too close to the planets and the Moon, and secondly, starting as close as possible to the true values. The second requirement can be met with an ad-hoc procedure, adopted for both tests constituting the demonstration run, that consists of three steps:

1. The first is a regular run of the GSR pipeline starting from input values comparable to those of AGIS.

2. The second starts from the same input values, but uses the BA reconstruction of the first step to remove the modulation from the known terms.

3. The third and final step is another run of the GSR pipeline in which the input attitude is the same as that in step 1, while the input values for the astrometric parameters are those of the solution of step 2 (the step 2 residual BA modulation is also subtracted from the known term as before).

We stress that, as long as we consider the simulated data a faithful reproduction of reality, the adopted procedure is not motivated just by the accuracy mismatch between the two astrometric models, but also by that between the accuracy of the RAMOD model and that of the Gaia measurements. Moreover, the third step takes into account the residual model-induced inaccuracy of the first order derivatives of the observable and for this reason we dub it “external iteration”, or EI run.

Indeed, it is worth noticing that this run might look like the so-called iteration for non-linearity, but it has actually another meaning. Indeed, the former is a standard procedure used when the input values are so approximate that the second order effects neglected in the linearized problems are still significantly large with respect to the measurement accuracy. Our EI, instead, is a numerical procedure that takes into account a first-order modelling accuracy issue. In short, it comes from the fact, already mentioned in Sect. 3.1, that it was chosen to model the effects of the planets only for the known terms, not including them in the coefficients of the linearized observation equation. This was done for practical reasons, and a detailed explanation, which would be too long and out of scope here, will be provided in a forthcoming paper (Vecchiato et al., in prep.). In this respect, when the current model will be replaced by a full-accuracy astrometric model (Crosta et al. 2017; Bertone et al. 2017) this ad-hoc procedure will become unnecessary. GSR, in fact, will have the same sensibility of AGIS to modelling errors and will match the Gaia measurement accuracy, so in principle it will be able to reach the required accuracy just after the first step.

The dataset utilized in these tests provides the true Gaia AL and AC measurements, which are computed from the coordinates of the simulated objects, the ephemerides of the planets and of the satellite, the NSL, and from the unperturbed instrument parameters using the implementation of the GREM model available when this dataset was produced. In comparing the results of our tests with those of AGIS, it should be considered that the former used a dataset containing 908 979 primary sources in the magnitude range 5.79 < G < 20.00, while AGIS, in its demonstration run, could use another version of the same dataset, not available anymore, that had more than 2 million primaries (Lindegren et al. 2012).

The perturbed values needed for the tests are generated with two different procedures for the sources and the attitude. The source parameters are perturbed with a Gaussian noise with zero average and a standard deviation of 20 mas for positions and annual proper motions. Moreover, negative parallaxes are not admitted, and when a random extraction produces a negative value, it is set to ϖ = 10⁻⁶ mas ≃ 4.8 × 10⁻¹⁵ rad⁹. The starting attitude values, instead, are obtained by perturbing the coefficients of the B-splines fit of the NSL (true attitude), with a Gaussian noise that produces a difference of about 10 mas between the true and perturbed orientations of the attitude axes. The initial separation between two successive nodes of the knot sequence is set to 240 s.

In test 2, the one with noisy observations, the perturbed observations are obtained by adding a Gaussian measurement noise to the true measured values. All the information needed to compute the latter are contained in the simulated dataset. The Gaussian perturbation, instead, is computed by generating a random extraction for each observation from a Gaussian distribution centered in zero and having a standard deviation corresponding to the observational error. This depends on the magnitude of the observed star, and it is obtained using the DPAC routines that implement the nominal AL/AC single-measurement error.

A star is declared solvable (for all of the 5 astrometric parameters) if it has at least 180 AF along-scan observations, and, at the same time, the difference between the observation times of its first and last observation is at least 1.5 years.

An important part of the coefficients module is the Attitude Definition Chain (ADC). The task of this piece of software is the definition of the knot sequence of the attitude, whose basic nominal separation can be adjusted to fit some constraints. Basically, there must be at least 20 observations between two adjacent knots, and the interval can be stretched up to four times the initial knot separation, namely 960 s. If after this time the minimum number of observations has not been reached yet, the B-spline sequence is segmented, thus a discontinuity in the attitude reconstruction is introduced.

The convergence condition of the LSQR algorithm are set to the most stringent requirement, namely to the machine precision accuracy. This condition is overridden if the estimated condition number of the system exceeds 10¹³, the number of iterations is larger than 50 000, or if the solver runs for more than 120 h. The preconditioning of the equation system is activated (to speed up convergence), and the six constraint equations needed to fix the intrinsic rank deficiency of the system are computed from the catalog values of “one-and-a-half stars”. Namely, we first select all stellar pairs according to the following geometric criteria:

$$\begin{array}{cc}\hfill \left|{\delta }_{1,2}\right|& <5\mathrm{deg}\hfill \end{array}$$(41)

$$\begin{array}{cc}\hfill 90-5\mathrm{deg}<|{\alpha }_2-{\alpha }_1|& <90+5\mathrm{deg}.\hfill \end{array}$$(42)

Final choice is done by choosing the brightest two among the above pairs, and the constraint equations are built by fixing α, δ, μ_α and μ_δ of the brightest and δ and μ_δ of the second brightest¹⁰.

5.2. Test 1: accuracy assessment with noise-free observations

The GSR pipeline filters out 460 primary sources that cannot be solved because they do not fit the minimum requirements explained in the previous section. Therefore, this leaves a system with 908 519 primaries and 660 599 degrees of freedom for the knot sequence. The attitude is represented in terms of MRP which implies that we have just 3 independent unknowns per knot, instead of the 4 constrained ones of the quaternion representation. In addition, the large scale instrument parameters that are estimated in the solution are associated to a total of 63 CCDs, each varying with a time scale of one month.

In this representation, therefore, the number of astrometric unknowns is 4 542 595, that of the attitude is 1 981 797, and finally, the number of AL and AC (large scale) instrument parameters is 45 360. The linearized system of equations is thus described by a design matrix of 6 569 752 columns (unknowns) and 1 330 628 523 rows (observations). In each of the three steps of this test, the solver converged to the machine-precision least-squares solution in something less than 20 000 iterations, with a condition number of ∼10⁶ (Fig. 7).

Fig. 7. Convergence plots of the first and third step. The LSQR algorithm implementation monitors the convergence status of the solution by computing two parameters, the 2-norm of the residuals vector (||r|| ≡ ||A x − b||), and the 2-norm ||AT r|| (dashed and solid lines respectively). As mentioned in Sect. 3.2, convergence to the unique least-squares solution is confirmed by the fact that the LSQR estimation of ||ATr|| is zero within the machine-precision accuracy.

Astrometric parameters. The first-step solution of the test reflects the effects of the mismodelling on the accuracy level attainable with the relativistic model implemented in GSR2, called RAMOD2d. The solution shows how this model, at first, reaches a systematic floor of ∼10 μas that, as expected given the noise-free observations, is independent from the magnitude range of the stars.

The second step provides a solution similar to the previous one because it starts from the same input catalog. Nonetheless, in this intermediate passage the BA modulation is removed from the measurements using its reconstruction from the first step. Although statistically similar to the previous one, in this way one obtains a better starting point for the final step. Indeed, jumping directly from the first to the third step, namely using the solution obtained without this “preemptive cleaning” of the known terms, would result in a less accurate final solution.

The third step, finally, gives the confirmation that the errors previously obtained are really the effect of modelling errors, since in this case the results are at the sub-μas level or better, as expected. This is also a confirmation that the adopted special procedure is able to recover the AGIS numerical accuracy, which implies that the same expectations have to apply to the second test.

Table 1 reports the results of the first and third steps. The second step is not reported because of the negligible statistical differences with respect to the first one. Following the methodology of Lindegren et al. (2012) we adopted the Robust Scatter Estimation (RSE¹¹) as a measure of the errors of the solution with respect to the true simulated values. It is worth noting that the medians of the parallaxes of the third step are one order of magnitude larger than those of the other unknowns. This residual noise can be neglected, as it is safely below the Gaia level of accuracy. Since preliminary results have shown that with the implementation of a full-accuracy astrometric model this feature disappears (Bertone et al. 2017), it can be attributed to a residual effect of the accuracy of the GSR astrometric model.

The all-sky plots of the astrometric solution for the third step are shown in Fig. 8. It is worth discussing in more detail the residual differences of these parameters. In general the largest differences show up close to the ecliptic, which is consistent with the fact that the lack of accuracy of RAMOD2d is basically due to imperfect modeling of the planets’ contribution to the null geodesic. Such discrepancy can be seen, in a plot of the GREM vs. RAMOD differences between the measurements, as a sign flip when the observing direction crosses the position of a planet. This sign flip translates in a corresponding sign flip of the astrometric parameters. Since the planets move along the ecliptic, this systematic effect is smoothed out along this plane, and we only see the sign flip in the orthogonal direction.

Fig. 8. All-sky map of the astrometric residual differences (in μas or μas yr−1) for the noise-free GSR solution after the third step.

Attitude parameters. Similar considerations hold for the attitude parameters, as shown in Table 2, where the numbers give the residual rotations around the three axes in μas. It is interesting to notice a residual, μas-level average for the e₂ axis which, as in the case of the parallaxes, is linked to the accuracy of the astrometric model.

Reconstruction of the basic angle modulation. Figure 9 shows the results of the BA modulation reconstruction for FoV1 as obtained from the L₀ (shift) Δη large scale instrument parameters. Those of FoV2 are exactly antisymmetric. The residuals are at the sub-sub-μas level (with average −0.0016 μas and standard deviation 0.015 μas for FoV1 and 0.0014 μas and standard deviation 0.015 μas for FoV2).

Fig. 9. Test 1 BA reconstruction for FoV1. That of FoV2 is not reported as it is the same with the opposite sign. The blue line represents the true modulation signal, while the red dots, which use the scale on the right side of the plot, are the differences between such signal and the final reconstruction after the three steps.

5.3. Test 2: realistic solution with noisy observations

As for the previous test, the GSR pipeline filters out 460 stars that cannot be solved because of the requirements explained in the previous section. This leaves a system with 908 519 stars and just a slightly different number (660 594) of degrees of freedom for the attitude. The large scale instrument parameters that are estimated in the solution are associated to the same total number of 63 CCDs, each varying with a time scale of one month. The number of astrometric unknowns is therefore 4 542 595, that of the attitude parameters is 1 981 782, and nothing changes for the instrument parameters (they are 45 360 as before). The system is thus described by a design matrix of 6 569 737 columns (unknowns) and 1 330 628 449 rows (observations). In each of the three steps of the demonstration run the solver converged to the machine-precision least-squares solution in something less than 42 000 iterations, with a condition number of ∼10⁶ (Fig. 10).

Fig. 10. Convergence plots of the first and third step of test 2.

Astrometric parameters. Table 3 reports, for each magnitude class, the median and the RSE for steps 1 and 3 of the GSR runs, as well as the corresponding results of the AGIS demonstration run.

As anticipated above, the first step produces an improved astrometric catalog, which however is still affected by systematics; μas-level medians and ratios of the ϖ errors to those of the other astrometric unknowns different from expectations are the numerical signatures of such systematics. In particular, it is known from both theoretical analysis and numerical simulations that the scanning law of Gaia should produce well-defined values for these ratios, specifically,

$$\begin{array}{c}\hfill {\overline{r}}_{{\alpha }^{\ast }}\equiv \frac{{\sigma }_{{\alpha }^{\ast }}}{{\sigma }_{\varpi }}=0.787,\end{array}$$(43)

and similarly r̄_δ = 0.699, r̄_μμ^∗ = 0.556, r̄_μδ = 0.496¹². These approximate estimations can thus be used as an additional check for the solutions, and the results from Table 3 show that the actual ratios are close to the expected values, and similar between them, for AGIS and for the step 3 of the GSR solution, whereas this is not the case for the step 1.

For test 2, the all-sky map plots do not show the systematics of the previous test. Indeed, in this test the remaining sub-μas systematic errors visible in the previous plots are completely negligible with respect to the Gaussian ones, which distribute uniformly on the sky. The plots for all the unknowns are quite similar, so in Fig. 11 we just show the case of the parallaxes, where the all-sky map is paired with the histogram of the normalized astrometric errors obtained by combining the data from each magnitude class. These data show that the solution of the third step, eventually, recovers the remaining modelling residuals, and that with this ad-hoc procedure GSR is able to recover an astrometric solution comparable to that of AGIS at the sub-μas level even with a less accurate astrometric model.

Fig. 11. All-sky map and normalized histogram of the parallax astrometric residual differences for the GSR solution after the third step. Units of the all-sky map are μas.

Attitude parameters. Similarly, we report in Table 4 the results for the attitude solutions from the first and third step of GSR, and those of AGIS. The numbers give the residual rotations around each axis in μas. The last column reports the ratio between the residual errors of e₂ and e₁ (y and x axes respectively), which from the geometry of the observations should be 1.34.

While the final accuracy of the astrometric parameters mostly depends on the number of observations per star, that of the attitude strongly depends on the total number of observations, and therefore on the number of stars involved in the sphere reconstruction and their magnitude. For this reason GSR, with 908 000 stars, cannot reach the same accuracy obtained by AGIS with about 2 300 000 stars. The degradation factor must be close to 1.6, which compares well to the 1.8 obtained in this test. Finally, Fig. 12 shows the histograms of the attitude errors for the x, y and z axes.

Fig. 12. Histograms of the attitude solution after the third step of test 2 for the x, y and z axes.

Reconstruction of the basic angle modulation. Figure 13 shows the results of the BA modulation reconstruction for FoV1 as obtained from the L₀ Δη large scale instrument parameters. Those of FoV2 are exactly antisymmetric. As in the case of AGIS, the errors are at the sub-μas level (with average −0.053 μas and standard deviation 0.357 μas for FoV1 and −0.047 μas and standard deviation 0.361 μas for FoV2).

Fig. 13. Test 2 BA reconstruction for FoV1. That of FoV2 is not reported as it is the same with the opposite sign. As for the previous test, the blue line represents the true modulation signal, while the red dots, which use the scale on the right side of the plot, are the differences between such signal and the final reconstruction after the three steps.

Observations’ residuals. The solution obtained after the third step is used to compute the AL and AC residuals of the condition equations. As long as the solution is close to the true values, the standard deviation of the distribution of the residuals for each star should approach the Gaussian noise used to simulate the measurement error. Table 5 reports such residuals for each magnitude class, along with the single-measurement error approximate values for the simulated data estimated with the ad-hoc empirical formulae used for Fig. 1. It has to be stressed that the values obtained with the empirical formulae are computed under the hypothesis of a Gaussian distribution with zero average for any object. For the actual residuals, the least-squares solution implies a close-to-zero average for the complete system only, and nothing can be said in general for subsets of observation equations. This means that the AL and AC residuals of each star will be generally distributed around a non-zero average, therefore the estimation by the empirical formula binned by magnitude class is likely to underestimate the real case, which is exactly what we observe in the table.

On the other hand, as we have anticipated above, we have to expect close-to-zero residuals for the complete system we actually solved (in the least-squares sense) that is the one weighted by the measurement errors via the weight matrix W, WAx = Wb. For this system we have to expect that ⟨W(Ax − b)⟩ ≃ 0. The actual result in this case is −0.50 μas, which is again at the sub-μas level as expected from the accuracy of the solution as compared to the true values.

6. Current status and future developments

Work is currently ongoing to further develop the AVU/GSR pipeline. In the following we detail the main issues that are going to be addressed in future developments.

6.1. Convergence speed and variance estimation.

As shown in the previous section, the number of iterations needed to reach a complete convergence can be very different for different situations. It is well known that such variation depends on the number of unknowns of the problem, however experience is showing, as it is clear already from the results of Test 1 and Test 2, that other factors can have a strong influence on the convergence speed. For example, this demonstration run used the Calibration parameters only to estimate a very specific signal injected in the BA, but it is expected that, when the full instrument model will be adopted, an even larger number of iterations will be needed to reach the convergence. It would be obviously an advantage to minimize such number.

One possibility is to use different kind of astrometric constraints. The full astrometric problem, unlike that represented by the AGIS approach, is intrinsically ill-conditioned, and its solution requires the introduction of six constraint equations that fix the orientation of the reference system for the positions and the proper motions. As explained in Sect. 5, currently GSR implements such equations by constraining the values of one-and-a-half stars, but other criteria are applicable. For example, it is possible to compute six “barycentric constraints” using the values of an arbitrary selection of stars. The determination of such constraints with a large set of stars usually reduces the iterations needed for the convergence.

Another possibility is simply to relax the convergence conditions. In order to guarantee the best possible solution from a purely numerical point of view, it was decided to set the convergence requirement to the most stringent value, that is to machine-precision level. However, it has been noticed that the solution δx stabilizes before the convergence parameters reach such level, and it was verified that setting it, for example, to the less demanding value of ∼10⁻¹³ would have produced an equivalent δx with a significantly reduced number of iterations. On the other hand, there is no obvious way to link the stability of the solution at the required numerical accuracy with the numerical value of ||A^Tr||, namely the residuals of the normal system that is used to estimate the convergence status of non-compatible systems. A possible practical solution is to check at regular intervals whether the variation of the solution falls below a given threshold. Technically, this option is easy to implement, since it can take advantage of a feature of our customized implementation of the LSQR algorithm, which produces intermediate solutions at regular intervals of iterations. Another option for the threshold value of the convergence condition is ǁA − Āǁ/ǁAǁ, as suggested in Paige & Saunders (1982), in which Ā are the true values of the coefficients’ matrix and the numerator can be estimated by knowing the uncertainties of the catalog values. Both these options, however, require further study to assess the reliability of the results obtained in this way. A drastic improvement in this respect, however, might have an undesired consequence.

Actually, it is also known that initially the LSQR algorithm underestimates the variances of the solution, and that such estimation improves with the number of iterations. A careful trade-off is therefore needed to get the fastest possible convergence without hampering the variances estimation or, in alternative, one might decide to resort to other approximate methods to estimate the latter. The GSR implementation of the LSQR algorithm always provides the variances along with the solution, but currently the former are used just for a first rough estimation of the quality of the solution. Further work is needed to assess the reliability of the variances estimation.

6.2. BA fit and calibration model.

As remembered in Sect. 2, the daily calibration does not guarantee a reconstruction of the instrument parameters at the final accuracy level required by the sphere reconstruction. Additional calibration parameters are thus added to the observation equations in order to recover any remaining un-calibrated residuals. However, the instrument model currently implemented in GSR is still inadequate to obtain a sphere solution coping with the accuracy goal of the Gaia catalog. Current work on real data is aimed at determining the improvements to the present instrument model needed to reach such an accuracy.

A similar reasoning has to be applied for the BA. The BAM pipeline(s) provides a daily reconstruction of the BA Variation (BAV), and it is likely that a cyclic reprocessing of these data can improve this reconstruction by providing a calibrated estimation of such variations. It is nonetheless appropriate to introduce further BA Correction (BAC) components in the observation equations to take into account a possible residual signal (Lindegren et al. 2018).

6.3. Comparison analysis

As described in Sect. 4.2, the GSR pipeline includes a Post-Solver module of comparison between the GSR and AGIS sphere reconstructions. The goal of this task is the internal validation of the primary star catalog astrometric parameters and associated formal uncertainties: this is carried out through a detailed analysis of the differences between the two solutions, thereby enlightening possible causes for discrepancies.

While the implementation of the VSH technique in the current GSR pipeline can be already used to identify and remove a residual rotation between the two catalogs, easily identified by the toroidal harmonics of degree one, the statistical robustness of the significance level test of other expansion coefficients, along with their interpretation in terms of plausible physical/geometrical effects, need a more extensive investigation, which is being addressed in the above cited forthcoming paper (Bucciarelli et al., in prep.).

Caution must be paid when the non-uniformity of the stellar distribution brakes the orthogonality of the base functions, with the consequence of giving rise to correlations among the VSH coefficients which must be taken into account; also, early truncation of the series expansion can generate spurious coefficients which result in false signal detection.

Another technique that we plan to adopt for the analysis of zonal errors is that of Infinitely Overlapping Circles (IOC; Taff et al. 1992; Bucciarelli et al. 1993): based on a statistical moving average naturally defined on the sphere; this technique is particularly easy to implement and can be successfully applied to extract local signals of scale length comparable or larger than the typical distance between neighboring sources.

6.4. Full-accuracy relativistic models

As shown in the present paper, GSR would benefit from the implementation of a relativistic model able to cope with the accuracy of the Gaia measurement. Work in this sense is in a very advanced stage, and will be reported in a forthcoming publication (Vecchaiato et al., in prep.).

6.5. Handling outliers, observations’ weighting and micro-events

Observations’ errors have to be used to build the weight matrix of the system. At present, this is done by using the expected uncertainty of the observed object (see Sect. 3.1) but, as explained in Lindegren et al. (2012), a more effective weighting can be obtained by considering both these uncertainties and the actual observations’ residuals, a procedure that in AGIS is also used to identify the outliers of the observations of a given source.

In addition to this, attitude undergoes uninterrupted perturbations because of the so-called micro-events, namely “micro-clanks” (small adjustments of the satellite structure) and “micro-hits” (due to the impact of micro-meteoroids with the Gaia satellite). These anomalies manifest themselves through non-nominal temporal variations of the AL and AC field angles, which can be interpreted as variations in the attitude scan rates. In the main pipeline they are estimated by a pre-processor running before the sphere solution, and successively removed from the observations. Failing to remove such perturbations would cause mas-level degradation to the accuracy of the sphere solution.

The above mentioned weighting procedure is currently under testing in GSR whereas, regarding the treatment of the micro-events, our pipeline is already able to remove them using the attitude corrections provided in the Gaia Main Database, while work is ongoing to implement GSR’s own procedure for micro-clanks and micro-hits.

7. Conclusions

The Astrometric Verification Unit is in charge of providing the DPAC with a pipeline able to realize a reconstruction of the global astrometric sphere independent from that of AGIS. This would allow the double-checking of the determination of the global reference system of Gaia. To this aim it is sufficient to reproduce the first stage of the process implemented by AGIS, namely the solution of the global sphere from the so-called primary sources. The absence of the secondary sources in such a pipeline, named AVU/GSR, is compensated for by the comparison task, whose goal is to identify and characterize any systematic discrepancy between the AGIS and GSR solutions larger than the expected accuracy of the Gaia catalog.

In order to guarantee a sufficient independence between the two solutions, AVU/GSR uses a different relativistic model for the astrometric observable, and a different parametrization for the attitude, based on MRP rather than quaternions. Moreover, the solution of the system of linearized equations is performed with a parallelized implementation of the LSQR full-iterative algorithm.

In this paper we showed that GSR managed to successfully reproduce the results of the demonstration run of the AGIS pipeline, as illustrated in Lindegren et al. (2012). This required the execution of two tests on simulated data for a 5 year mission duration, both solved for astrometric, attitude and instrument parameters.

The first one had the goal of assessing the sub-μas numerical accuracy of the GSR pipeline, which was obtained by using perturbed starting values and error-free observational data; the second one had to mirror the outputs of the AGIS demonstration run with perturbed values for source and attitude parameters, plus random measurement errors compatible with those expected by Gaia. Moreover, large-scale instrument calibrations were used to reconstruct a periodic signal injected in the BA value. This test was designed to produce an astrometric catalog with errors compatible with those of the final Gaia catalog.

Both tests performed according to expectations, even if the current implementation of the GSR relativistic model is less accurate than the AGIS one, its accuracy being limited by the influence of solar system bodies other than the Sun on the computation of null geodesics. Despite these limitations, the current pipeline reaches the level needed by the Gaia measurements at appropriate elongation from the perturbing object. Moreover, the model is improved by finding approximate corrections to the light deflection effect induced by the planets and the Moon, thereby extending the required accuracy to larger portions of the sky.

The influence of this issue on the final solution was investigated in this work. Test 1 revealed that GSR is critically sensible to the catalog errors because of the limited accuracy of the model; at the same time it allowed to prove that the required quality of the solution can be reached by performing an external iteration. Test 2, instead, showed that, as it was expected, the manifestation of modelling errors depends on the accuracy of Gaia measurements, and therefore at the faint end of the stellar sample the AGIS and the GSR solutions are compatible already after the first iteration, while this is true for the entire catalog only after the external iteration. This special procedure will likely be unnecessary when the ongoing work of implementing an astrometric model at the same intrinsic accuracy of GREM will be completed.

This work has put to evidence that GSR is sufficiently mature to start processing real data, and current developments aim at providing this pipeline with the features needed to meet the accuracy requirements of the Gaia solution. Some of them are made necessary by the actual behavior of the instrument and have already been tackled by AGIS; others involve a more sophisticated treatment of the comparison task, and are specific to GSR. Aim of the GSR group is to have all the needed features implemented and tested for by the next Gaia Data Release, in order to team up with AGIS and thus contribute to the production of the Gaia catalog, as foreseen in the DPAC plans.

Acknowledgments

This work was supported by the Agenzia Spaziale Italiana (ASI) through contract 2014-025-R.1.2015 to the Italian Istituto Nazionale di Astrofisica (INAF) and contract 2016-17-I.0 to the Aerospace Logistics Technology Engineering Company (ALTEC S.p.A.), and INAF. The authors wish to thank Michael Biermann (DPAC CU3 leader) and Gonzalo Gracia (DPAC Project Office Coordinator) for their support and advice during the evaluation of the results of the Demonstration Run, and the anonymous referees for their helpful comments.

References

All Tables

All Figures

	Fig. 1. Approximate estimation of the Gaia single-measurement along- and across-scan accuracy (AL and AC) as function of the G magnitude. The oscillating shape up to mag 12 is due to the use of CCD gates in order to avoid saturation for bright objects. The estimation is the result of a best-fit model that averages over the whole Astrometric Field (AF).
In the text

	Fig. 2. Frequency of observations as function of equatorial coordinates due to Gaia scanning law (blue ≈50–yellow ≈500).
In the text

	Fig. 3. Parametrization of the Gaia Nominal Scanning Law with respect to an equatorial reference system $\{{\mathbf{E}}_{\widehat{1}},{\mathbf{E}}_{\widehat{2}},{\mathbf{E}}_{\widehat{3}}\}$. NS is the Nominal Sun.
In the text

	Fig. 4. Representation of the AL (ϕ) and AC (ζ) Gaia measurements with respect to the satellite attitude $\{{\mathbf{e}}_{\widehat{1}},{\mathbf{e}}_{\widehat{2}},{\mathbf{e}}_{\widehat{3}}\}$.
In the text

	Fig. 5. Schematic representation of the block-iterative procedure used by AGIS (left) and of the fully iterative one used by GSR (right). See the text for explanation.
In the text

	Fig. 6. Schematic representation of the AVU/GSR pipeline.
In the text

Fig. 7. Convergence plots of the first and third step. The LSQR algorithm implementation monitors the convergence status of the solution by computing two parameters, the 2-norm of the residuals vector (||r|| ≡ ||A x − b||), and the 2-norm ||AT r|| (dashed and solid lines respectively). As mentioned in Sect. 3.2, convergence to the unique least-squares solution is confirmed by the fact that the LSQR estimation of ||ATr|| is zero within the machine-precision accuracy.

In the text

	Fig. 8. All-sky map of the astrometric residual differences (in μas or μas yr−1) for the noise-free GSR solution after the third step.
In the text

	Fig. 9. Test 1 BA reconstruction for FoV1. That of FoV2 is not reported as it is the same with the opposite sign. The blue line represents the true modulation signal, while the red dots, which use the scale on the right side of the plot, are the differences between such signal and the final reconstruction after the three steps.
In the text

	Fig. 10. Convergence plots of the first and third step of test 2.
In the text

	Fig. 11. All-sky map and normalized histogram of the parallax astrometric residual differences for the GSR solution after the third step. Units of the all-sky map are μas.
In the text

	Fig. 12. Histograms of the attitude solution after the third step of test 2 for the x, y and z axes.
In the text

	Fig. 13. Test 2 BA reconstruction for FoV1. That of FoV2 is not reported as it is the same with the opposite sign. As for the previous test, the blue line represents the true modulation signal, while the red dots, which use the scale on the right side of the plot, are the differences between such signal and the final reconstruction after the three steps.
In the text