Issue 
A&A
Volume 620, December 2018



Article Number  A40  
Number of page(s)  19  
Section  Celestial mechanics and astrometry  
DOI  https://doi.org/10.1051/00046361/201833254  
Published online  27 November 2018 
The global sphere reconstruction (GSR)
Demonstrating an independent implementation of the astrometric core solution for Gaia
^{1} INAF – Astrophysical Observatory of Torino, Pino Torinese, Italy
email: alberto.vecchiato@inaf.it
^{2} INAF – Astrophysical Observatory of Catania, Italy
^{3} EURIX, Torino, Italy
^{4} ALTEC, Torino, Italy
Received:
18
April
2018
Accepted:
4
September
2018
Context. The Gaia ESA mission will estimate the astrometric and physical data of more than one billion objects, providing the largest and most precise catalog of absolute astrometry in the history of astronomy. The core of this process, the socalled global sphere reconstruction, is represented by the reduction of a subset of these objects which will be used to define the celestial reference frame. As the HIPPARCOS mission showed, and as is inherent to all kinds of absolute measurements, possible errors in the data reduction can hardly be identified from the catalog, thus potentially introducing systematic errors in all derived work.
Aims. Following up on the lessons learned from HIPPARCOS, our aim is thus to develop an independent sphere reconstruction method that contributes to guarantee the quality of the astrometric results without fully reproducing the main processing chain.
Methods. Indeed, given the unfeasibility of a complete replica of the data reduction pipeline, an astrometric verification unit (AVU) was instituted by the Gaia Data Processing and Analysis Consortium (DPAC). One of its jobs is to implement and operate an independent global sphere reconstruction (GSR), parallel to the baseline one (AGIS, namely Astrometric Global Iterative Solution) but limited to the primary stars and for validation purposes, to compare the two results, and to report on any significant differences.
Results. Tests performed on simulated data show that GSR is able to reproduce at the subμas level the results of the AGIS demonstration run.
Conclusions. Further development is ongoing to improve on the treatment of real data and on the software modules that compare the AGIS and GSR solutions to identify possible discrepancies above the tolerance level set by the accuracy of the Gaia catalog.
Key words: astrometry / reference systems / catalogs / methods: numerical / space vehicles
© 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 fiveparameter astrometric catalog (i.e., including positions, parallaxes and the two components of the proper motions) at the 10–1000 μarcsecondlevel (μ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^{3} observations per object, each at the ∼0.1 − 1 mas level, and it relies on the selfcalibration 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 spectrophotometric 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 microthrusters of the onboard 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 quasiconstant speed. As a consequence, in principle the duration of a single observation is the same for any star, and therefore the singlemeasurement accuracy is strongly magnitudedependent (Fig. 1). Moreover, the scanning law parameters induce a complete scan of the celestial sphere every six months, but with a nonuniform 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 singlemeasurement along and acrossscan 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 bestfit 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^{7} − 10^{8} “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 twosteps 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 crosschecking 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 crosschecking 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:

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;

fourvectors 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;

the signature of the metric g_{αβ} is +2;

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

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

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

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 constantrate satellite spin around the axis , 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 . NS is the Nominal Sun. 
In this geometrical model the two FieldsofView (FoVs) of Gaia lie on the satellite’ scanning plane and are pointing symmetrically with respect to , separated by an (approximately) constant Basic Angle (BA) Γ. If r is the position vector of a pointlike source, in the SRS the basic observables of Gaia can be represented by its alongscan measurement (AL) that is the abscissa ϕ, which is the angle between the projection of r on the instantaneous scanning plane and , and the acrossscan measurement (AC) ζ (Fig. 4).
Fig. 4.
Representation of the AL (ϕ) and AC (ζ) Gaia measurements with respect to the satellite attitude . 
These quantities can be represented in terms of the direction cosines of r with respect to the SRS
while the direction cosine of r with respect to the axis e_{â} is
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 selfcalibrating 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
2.3. Gaia accuracy and the relativistic observable
The target accuracy of Gaia is at the submas 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)^{1}. 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 fourdimensional 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
where
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 fourvector 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 socalled light deflection.
Any measurement is always defined with respect to a specific observer, identified with its fourvelocity 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 fourdimensional 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_{1} and r_{2} satisfies the relation cosα = r_{1} ⋅ r_{2}/(r_{1}r_{2}), with r = (r⋅r)^{1/2}, the corresponding general relativistic expression will be
It is worth noticing that in this way the aberration effects are automatically accounted for by the inclusion of the observer’s fourvelocity in h_{αβ}.
Finally, one standard way to handle the SRS is based on the socalled tetrad formalism (de Felice & Bini 2010), which is based on the possibility of identifying a locally Lorentzian reference system associated to a specific observer , that is a set of fourdimensional axes e_{α̂}, α = 0, …, 3 called tetrad defined by the conditions
In this way, we are inducing the same 3+1 splitting of above with the addition of a set of three fourdimensional axes {e_{â}}, a = 1, 2, 3 belonging to the 3D space of u^{α}. Since by definition and , these axes are spatial (in the sense that they are orthogonal to 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
On the other hand, the AL and AC measurements of Eqs. (1) and (2) are defined in the 3D subspace of 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:

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;

find the appropriate expression of the barycentric motion of Gaia in the BCRS;

use the above fourvelocity and the BCRS metric to compute the spatial projector h_{αβ};

find the Gaia relativistic tetrad by a twosteps 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 ;

(a)

use k^{α}, h_{αβ} and in Eq. (10) to compute the needed direction cosines;
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 postNewtonian “bookkeeping” (Will 1993). In this way, considering that the typical velocities in the Solar System are ≃30 km s^{−1}, the socalled 1PN order of (ν/c)^{2} would correspond to a ∼ 10^{−8} rad in angular accuracy, that is to the mas level, allowing one to set the befitting model accuracy at the 1.5PN order (ν/c)^{3}, corresponding to 10^{−12} 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 Gaialevel 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 singlebody PPNSchwarzschild 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^{2}, of the satellite position (r_{s},θ_{s},ϕ_{s}), and of the PPNγ parameter:
First of all, it is worth stressing that the accuracy of the Gaia measurement allows for the estimation of the γ parameter as a byproduct 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 λ_{*}
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
where t_{0} 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)
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 fourvelocity of the satellite can be written as where ν_{s} is the coordinate velocity of the satellite and, using the IAU resolutions (Kaplan 2005),
Here and U is the (BCRS) Nbody gravitational potential of the Solar System at the Gaia location, which implies that the fourvelocity, 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 fourvector by setting . Then one has to fix the spatial triad by considering the following transformations (Crosta & Vecchiato 2010; Bianchi et al. 2011):

the origin of the tetrad is located at the BCRS coordinates of the satellite barycentre, which induces a first transformation to the socalled 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;

the satellite barycentre moves with fourvelocity 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^{3}. If the fourvelocity 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 ;

the satellite orientation, that is the SRS e_{â}, can be finally obtained by a Euclidean transformation of its spatial axes, , 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) σ_{1}, σ_{2} and σ_{3} 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)^{2} order. Basically, the main requirement is that the metric has the form
where, for the null geodesic integration, we retained only the (ν/c)^{2}order perturbation term given by the gravitational potential of a massive body , so that , while the attitude model can include also higher order terms, namely 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 timedependent 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 micropropulsion system introduces a highfrequency 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 timedependent polynomial functions B(t) of degree M − 1 called BSplines (Ahlberg & Nilson 1967) whose characteristics are summarized in the following.
The expansion of the jth component of a generic representation, say S, reads
where 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 subintervals identified by a sequence of instants {τ_{k}}, k = 0, …, K called nodes, which constitutes the socalled support of the series. The BSplines are useful in our case because their support is minimal, that is B_{n}(t) ≠ 0 only in M subintervals, so if τ_{n} < t < τ_{n + 1} then,
in the sense that at each time the expansion depends only on M unknowns. For Gaia the expansion is in cubic Bsplines, that is M = 4, and therefore
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 BSplines lies in the possibility of exploiting the linearization of the observation equations to implement the socalled 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 lineofprinciple 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 longterm 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 acrossscan 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 spectrumdependent shifts dη are introduced.
In practice the AL measurement ϕ in Eq. (4) is defined as
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
in which η^{0} is the nominal AL field angle and Δη_{rfnj}, δη_{nm}, are the above mentioned calibration unknowns, namely the AL large and small scale geometric corrections, and the magnitude and spectrumdependent shifts. The indices indicate the dependencies of these parameters with respect to the instrument configuration.
The largescale 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 computed at the normalized AC pixel coordinate
where 14 ≤ μ ≤ 1979 is the pixel coordinate of the measurement and μ_{0} = 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 + μ_{0} to pixcols + μ_{0}, where pixcols is the total number of illuminated pixels. The same indexes are used to identify the dependencies of the magnitude and spectrumdependent shifts, which need a reference magnitude and frequency G_{ref} and ν_{ref} respectively.
Similarly, the AC measurement ζ of Eq. (5) is
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 leastsquares 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^{8} 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^{3} times during the mission lifetime, so the total number of observations is n_{obs} ∼ 10^{11}.
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 nonlinear 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
where , , , 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
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 or to , 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^{4}.
The resulting system of equations is intrinsically rankdeficient when astrometric and attitude parameters have to be solved at the same time. This wellknown issue comes from the invariance of the solution for a rigid threecomponents rotation and threecomponents 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^{5}
is introduced, where ϵ_{j} is the estimated standard deviation of the jth observation and ϵ_{0} is a reference value, currently assumed to the estimated uncertainty at magnitude G = 21.5. The system thus becomes
The solution in the leastsquares sense of such system provides δx̂, namely the bestfit estimation of δ x, that is used to update the catalog values x̄ to the improved estimation of the true values x̃, that is
The formal leastsquares solution of the system of Eq. (28) is δ x = [(WA)^{T}(WA)]^{−1}(WA)^{T}W b, where [(WA)^{T}(WA)]^{−1} 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 (and to )^{6}, 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^{3}. Since the normal matrix A^{T}A is n_{unk} × n_{unk} with n_{unk} ∼ 6 × 10^{8}, we can immediately realize that in the Gaia case this problem can be addressed only resorting to iterative algorithms. A bruteforce approach like this, in fact, would require about 10^{26} FLOPs a requirement which cannot be reduced to an acceptable level even taking into account the sparsity of A.
AGIS uses a blockiterative 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 attitudeinduced 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 blockiterative 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/CC++ 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 PCLSQR (Baur et al. 2008) a conjugate gradientbased algorithm, originally proposed by Paige & Saunders (1982). As for any iterative algorithm, this is equivalent to computing, at each iteration (i), an approximate solution
and then evaluating the vector of the residuals
which has to be minimized in the leastsquares 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 2norm 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 noncompatible systems it converges to a positive finite limit;

the 2norm of A^{T} r^{(i)} (namely the norm of the residuals of the normal system A^{T} b = A^{T}Aδ x) is lower than a threshold value; this is the condition used to guarantee the convergence of noncompatible systems to a solution in the leastsquares sense; it is worth stressing here that the equation system to be solved for the global astrometric sphere reconstruction is noncompatible^{7}, and indeed its solution is obtained in the leastsquares 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^{T}A)^{−1}, 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 offdiagonal 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 nonzero 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 matrixvector 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 communicationfree (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
so that in general the catalog attitude is a set of initial values of the BSpline expansion coefficients . Equation (25) can thus be written as
in which we have exploited the fact that
and that, from Eq. (18),
In this way the updated attitude at a generic time t is
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 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 . In formulae, we could write at any time where
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 , so that in the end , 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:
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 nondifferential 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 , j = 1, 2, 3 relative to the kth unknown of each Rodrigues parameter σ_{j} ^{8}, their ratios do not depend from the Bspline 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_{2} − t_{1} between two observations at t_{1} and t_{2} respectively which is short enough to keep almost constant between t_{1} and t_{2}. As pointed out above, this “constancy” in principle does not depend on the explicit formulation of the coefficients, since it is sufficient to have timedependent 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 remain almost the same among the different observations or, equivalently, over which the proportionality holds approximately for each index k.
Since each of these values is represented into a different column of the system of equations, this quasiproportionality produces an extremely illconditioned design matrix. Such an illconditioning 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 “Presolver” 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 “Postsolver” 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 DPCTALTEC, with the exception of the Solver module which is written in C/C++ and runs at the DPCTCINECA. 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. PostSolver 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 rankdeficient, and as a minimum it can be solved except for a sixparameters 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 derotate one solution, that is to bring both catalogs into a common reference system, before attempting a comparison between them. Moreover, both the derotation 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 leastsquares reconstructions of purely Euclidean transformations and statistical analysis of coordinate differences.
GSR can use standard statistical algorithms, like χ^{2} tests and KolmogorovSmirnov, 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 leastsquares 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 computationallyintensive, so its implementation is the only one besides the system solution that can require a nonembarrassingly parallel coding (to be executed at the DPCTCINECA) 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 noisefree 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 longperiod modulation of the BA is injected in the known terms (observations) and reconstructed by means of the L_{0} largescale 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 leastsquares 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 leastsquares solution, a statement which has to be intended in comparison with the blockiterative 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)^{3} 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)^{2} PPNSchwarzschild 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 nonSchwarzschild 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 adhoc procedure, adopted for both tests constituting the demonstration run, that consists of three steps:

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

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.

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 modelinduced 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 socalled iteration for nonlinearity, 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 firstorder 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 fullaccuracy astrometric model (Crosta et al. 2017; Bertone et al. 2017) this adhoc 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^{−6} mas ≃ 4.8 × 10^{−15} rad^{9}. The starting attitude values, instead, are obtained by perturbing the coefficients of the Bsplines 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 singlemeasurement error.
A star is declared solvable (for all of the 5 astrometric parameters) if it has at least 180 AF alongscan 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 Bspline 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^{13}, 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 “oneandahalf stars”. Namely, we first select all stellar pairs according to the following geometric criteria:
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^{10}.
5.2. Test 1: accuracy assessment with noisefree 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 machineprecision leastsquares solution in something less than 20 000 iterations, with a condition number of ∼10^{6} (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 2norm of the residuals vector (r ≡ A x − b), and the 2norm A^{T} r (dashed and solid lines respectively). As mentioned in Sect. 3.2, convergence to the unique leastsquares solution is confirmed by the fact that the LSQR estimation of A^{T}r is zero within the machineprecision accuracy. 
Astrometric parameters. The firststep 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 noisefree 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^{11}) 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 fullaccuracy 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.
Astrometric results (estimated minus true) for test 1.
The allsky 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.
Allsky map of the astrometric residual differences (in μas or μas yr^{−1}) for the noisefree 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, μaslevel average for the e_{2} axis which, as in the case of the parallaxes, is linked to the accuracy of the astrometric model.
Attitude results (in μas) for test 1.
Reconstruction of the basic angle modulation. Figure 9 shows the results of the BA modulation reconstruction for FoV1 as obtained from the L_{0} (shift) Δη large scale instrument parameters. Those of FoV2 are exactly antisymmetric. The residuals are at the subsubμ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 machineprecision leastsquares solution in something less than 42 000 iterations, with a condition number of ∼10^{6} (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; μaslevel 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 welldefined values for these ratios, specifically,
and similarly r̄_{δ} = 0.699, r̄_{μμ∗} = 0.556, r̄_{μδ} = 0.496^{12}. 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.
Astrometric results (estimated minus true) for test 2.
For test 2, the allsky 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 allsky 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 adhoc 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.
Allsky map and normalized histogram of the parallax astrometric residual differences for the GSR solution after the third step. Units of the allsky 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_{2} and e_{1} (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.
Attitude results (in μas) for test 2.
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_{0} Δη 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 singlemeasurement error approximate values for the simulated data estimated with the adhoc 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 leastsquares solution implies a closetozero 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 nonzero 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.
AL and AC standard deviations of the residuals of the nonweighted equation system (“System” columns) after step 3 of test 2 per magnitude class compared with the singlemeasurement error used to generate the simulated data (“Est.” columns).
On the other hand, as we have anticipated above, we have to expect closetozero residuals for the complete system we actually solved (in the leastsquares 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 illconditioned, 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 oneandahalf 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 machineprecision 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^{−13} 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^{T}r, namely the residuals of the normal system that is used to estimate the convergence status of noncompatible 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 tradeoff 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 uncalibrated 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 PostSolver 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 nonuniformity 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. Fullaccuracy 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 microevents
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 socalled microevents, namely “microclanks” (small adjustments of the satellite structure) and “microhits” (due to the impact of micrometeoroids with the Gaia satellite). These anomalies manifest themselves through nonnominal 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 preprocessor running before the sphere solution, and successively removed from the observations. Failing to remove such perturbations would cause maslevel 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 microevents, 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 microclanks and microhits.
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 doublechecking 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 socalled 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 fulliterative 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 errorfree 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, largescale 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.
The ICRF is the socalled realization of the International Celestial Reference System (ICRS; Souchay & FeisselVernier 2006). Likewise, the Gaia catalog is the realization of the BCRS.
This tetrad is thus the CoMRS of the Gaia DPAC nomenclature (Lindegren et al. 2012).
Acknowledgments
This work was supported by the Agenzia Spaziale Italiana (ASI) through contract 2014025R.1.2015 to the Italian Istituto Nazionale di Astrofisica (INAF) and contract 201617I.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
 Ahlberg, J. H., & Nilson, E. N. 1967, The Theory of Splines and their Applications (Academic Press), Mathematics in Science and Engineering, 38 [Google Scholar]
 Arfken, G. B., & Weber, H. J. 2012, Mathematical Methods for Physicists, 7th edn. (Academic Press) [Google Scholar]
 Bandieramonte, M., Becciani, U., Vecchiato, A., et al. 2012, IEEE 23rd International WETICE Conference, 167 [Google Scholar]
 Baur, O., Austen, G., & Kusche, J. 2008, J. Geod., 82, 207 [NASA ADS] [CrossRef] [Google Scholar]
 Becciani, U., Sciacca, E., & Bandieramonte, M. 2014, International Conference on High Performance Computing Simulation (HPCS), 104 [Google Scholar]
 Berghea, C. T., Makarov, V. V., Frouard, J., et al. 2016, ApJ, 152, 53 [CrossRef] [Google Scholar]
 Bertone, S., Vecchiato, A., Bucciarelli, B., et al. 2017, A&A, 608, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bianchi, L., Vecchiato, A., & Bucciarelli, B. 2011, Global Sphere Reconstruction Attitude model, Tech. Rep. GAIAC3TNINAFLB00101, available from https://www.cosmos.esa.int/web/gaia/publicdpacdocuments, DPAC Livelink [Google Scholar]
 Bini, D., Crosta, M. T., & de Felice, F. 2003, Class. Quant. Grav., 20, 4695 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Bombrun, A., Lindegren, L., Holl, B., & Jordan, S. 2010, A&A, 516, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brown, T. M., Casertano, S., Strader, J., et al. 2018, ApJ, 856, L6 [NASA ADS] [CrossRef] [Google Scholar]
 Bucciarelli, B., Taff, L. G., & Lattanzi, M. G. 1993, J. Stat. Comput. Simul., 48, 29 [CrossRef] [Google Scholar]
 Bucciarelli, B., Abbas, U., Vecchiato, A., & Lattanzi, M. G. 2011a, Comparison of two Gaia sphere solutions using orthonormal bases on the sphere, Tech. Rep. GAIAC3TNINAFBB00201, available from https://www.cosmos.esa.int/web/gaia/publicdpacdocuments, DPAC Livelink [Google Scholar]
 Bucciarelli, B., Abbas, U., Vecchiato, A., & Lattanzi, M. G. 2011b, Link between two Gaia intermediate sphere solutions, Tech. Rep. GAIAC3TNINAFBB00101, available from https://www.cosmos.esa.int/web/gaia/publicdpacdocuments, DPAC Livelink [Google Scholar]
 Buluç, A., Fineman, J. T., Frigo, M., et al. 2009, in Proceedings of the Twentyfirst Annual Symposium on Parallelism in Algorithms and Architectures (New York, NY, USA: ACM), SPAA ’09, 233 [Google Scholar]
 Busonero, D. 2012, Modeling, Systems Engineering, and Project Management for Astronomy V, Proc. SPIE, 8449, 84490F [CrossRef] [Google Scholar]
 Busonero, D., Lattanzi, M. G., Gai, M., et al. 2014, Modeling, Systems Engineering, and Project Management for Astronomy VI, Proc. SPIE, 9150, 91500K [CrossRef] [Google Scholar]
 Casertano, S., Riess, A. G., Bucciarelli, B., & Lattanzi, M. G. 2017, A&A, 599, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Crosta, M., & Vecchiato, A. 2010, A&A, 509, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Crosta, M., Geralico, A., Lattanzi, M. G., & Vecchiato, A. 2017, Phys. Rev. D, 96, 104030 [NASA ADS] [CrossRef] [Google Scholar]
 de Felice, F., & Bini, D. 2010, Classical Measurements in Curved SpaceTimes (Cambridge: Oxford University Press) [CrossRef] [Google Scholar]
 de Felice, F., Lattanzi, M. G., Vecchiato, A., & Bernacca, P. L. 1998, A&A, 332, 1133 [NASA ADS] [Google Scholar]
 de Felice, F., Bucciarelli, B., Lattanzi, M. G.,& Vecchiato, A. 2001, A&A, 373, 336 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fey, A. L., Gordon, D., Jacobs, C. S., et al. 2015, AJ, 150, 58 [NASA ADS] [CrossRef] [Google Scholar]
 Frouard, J., Johnson, M., Fey, A., et al. 2018, Amer. Astron. Soc. Meet. Abstr., 231, 436.02 [NASA ADS] [Google Scholar]
 Gaia Collaboration (Prusti, T., et al.) 2016, A&A, 595, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Green, R. M. 1985, Spherical Astronomy (Cambridge: Cambridge University Press) [Google Scholar]
 Guo, R. 2008, Dimplomarbeit, University of Stuttgart, Institute of Geodesy [Google Scholar]
 Hill, E. L. 1954, Am. J. Phys., 22, 211 [NASA ADS] [CrossRef] [Google Scholar]
 Holl, B., & Lindegren, L. 2012, A&A, 543, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Holl, B., Lindegren, L., & Hobbs, D. 2012, A&A, 543, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kaplan, G. H. 2005, The IAU resolutions on astronomical reference systems, time scales, and earth rotation models: explanation and implementation, Tech. rep., U.S. Naval Observatory [Google Scholar]
 Kostina, E., & Kostyukova, O. 2012, International Series of Numerical Mathematics, 160, 197 [Google Scholar]
 Lindegren, L., Lammers, U., Hobbs, D., et al. 2012, A&A, 538, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lindegren, L., Lammers, U., Bastian, U., et al. 2016, A&A, 595, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lindegren, L., Hernández, J., Bombrun, A., et al. 2018, A&A, 616, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Makarov, V. V., Fabricius, C., & Frouard, J. 2017, ApJ, 840, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Mignard, F., & Drimmel, R. E. 2007, DPAC: Proposal for the Gaia Data Processing, Tech. Rep. GAIACDSPDPACFM03002, available from https://www.cosmos.esa.int/web/gaia/publicdpacdocuments, DPAC Livelink [Google Scholar]
 Misner, C. W., Thorne, K. S., & Wheeler, J. A. 1973, Gravitation (San Francisco: W.H. Freeman and Co.) [Google Scholar]
 O’Mullane, W., Lammers, U., Lindegren, L., Hernandez, J., & Hobbs, D. 2011, Exp. Astron., 31, 215 [NASA ADS] [CrossRef] [Google Scholar]
 Paige, C. C., & Saunders, M. A. 1982, ACM Trans. Math. Softw., 8, 43 [CrossRef] [MathSciNet] [Google Scholar]
 Qi, Z., Yu, Y., Bucciarelli, B., et al. 2015, AJ, 150, 137 [NASA ADS] [CrossRef] [Google Scholar]
 Riva, A., Gai, M., Lattanzi, M. G., et al. 2014, in RRev. Mex. Astron. Astrof., 45, 35 [NASA ADS] [Google Scholar]
 Schaub, H., & Junkins, J. L. 1996, J. Astronaut. Sci., 44, 1 [Google Scholar]
 Souchay, J., & FeisselVernier, M. 2006, IERS Technical Note, 34 [Google Scholar]
 Taff, L. G., Bucciarelli, B., & Lattanzi, M. G. 1992, ApJ, 392, 746 [NASA ADS] [CrossRef] [Google Scholar]
 Vecchiato, A., Lattanzi, M. G., Bucciarelli, B., et al. 2003, A&A, 399, 337 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
 Vecchiato, A., Abbas, U., & Bandieramonte, M. 2012, Software and Cyberinfrastructure for Astronomy II, Proc. SPIE, 8451, 84513C [CrossRef] [Google Scholar]
 Will, C. M. 1993, Theory and Experiment in Gravitational Physics (Cambridge University Press), 396 [Google Scholar]
 Zacharias, N., Finch, C. T., Girard, T. M., et al. 2013, AJ, 145, 44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Tables
AL and AC standard deviations of the residuals of the nonweighted equation system (“System” columns) after step 3 of test 2 per magnitude class compared with the singlemeasurement error used to generate the simulated data (“Est.” columns).
All Figures
Fig. 1.
Approximate estimation of the Gaia singlemeasurement along and acrossscan 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 bestfit 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 . NS is the Nominal Sun. 

In the text 
Fig. 4.
Representation of the AL (ϕ) and AC (ζ) Gaia measurements with respect to the satellite attitude . 

In the text 
Fig. 5.
Schematic representation of the blockiterative 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 2norm of the residuals vector (r ≡ A x − b), and the 2norm A^{T} r (dashed and solid lines respectively). As mentioned in Sect. 3.2, convergence to the unique leastsquares solution is confirmed by the fact that the LSQR estimation of A^{T}r is zero within the machineprecision accuracy. 

In the text 
Fig. 8.
Allsky map of the astrometric residual differences (in μas or μas yr^{−1}) for the noisefree 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.
Allsky map and normalized histogram of the parallax astrometric residual differences for the GSR solution after the third step. Units of the allsky 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.