# Application of time transfer functions to Gaia's global astrometry - Validation on DPAC simulated Gaia-like observations

Stefano Bertone, Alberto Vecchiato, Beatrice Bucciarelli, Mariateresa Crosta, Mario G. Lattanzi, Luca Bianchi, Marie-Christine Angonin, Christophe Le Poncin-Lafitte

► To cite this version:

Stefano Bertone, Alberto Vecchiato, Beatrice Bucciarelli, Mariateresa Crosta, Mario G. Lattanzi, et al. Application of time transfer functions to Gaia's global astrometry - Validation on DPAC simulated Gaia-like observations. Astronomy & Astrophysics - A&A, 2017, 608, pp.A83. 10.1051/0004-6361/201731654. hal-01669662

HAL Id: hal-01669662

https://hal.science/hal-01669662v1

Submitted on 27 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Application of time transfer functions to Gaia's global astrometry
## Validation on DPAC simulated Gaia-like observations

Stefano Bertone<fnref n="1" />,<fnref n="2" />,*, Alberto Vecchiato<fnref n="1" />, Beatrice Bucciarelli<fnref n="1" />,<fnref n="3" />, Mariateresa Crosta<fnref n="1" />, Mario G. Lattanzi<fnref n="1" />, Luca Bianchi<fnref n="4" />, Marie-Christine Angonin<fnref n="2" />, and Christophe Le Poncin-Lafitte<fnref n="2" />

Received 26 July 2017 / Accepted 31 August 2017

## ABSTRACT

**Context.** A key objective of the ESA Gaia satellite is the realization of a quasi-inertial reference frame at visual wavelengths by means of global astrometric techniques. This requires accurate mathematical and numerical modeling of relativistic light propagation, as well as double-blind-like procedures for the internal validation of the results, before they are released to the scientific community at large.

**Aims.** We aim to specialize the time transfer functions (TTF) formalism to the case of the Gaia observer and prove its applicability to the task of global sphere reconstruction (GSR), in anticipation of its inclusion in the GSR system, already featuring the Relativistic Astrometric MODel (RAMOD) suite, as an additional semi-external validation of the forthcoming Gaia baseline astrometric solutions.

**Methods.** We extended the current GSR framework and software infrastructure (GSR2) to include TTF relativistic observation equa-
tions compatible with Gaia's operations. We used simulated data generated by the Gaia Data Processing and Analysis Consortium
(DPAC) to obtain different least-squares estimations of the full (five-parameter) stellar spheres and gauge results. These were com-
pared to analogous solutions obtained with the current RAMOD model in GSR2 (RAMOD@GSR2) and to the catalog generated with
the Gaia RELativistic Model (GREM), the model baselined for Gaia and used to generate the DPAC synthetic data.

**Results.** Linearized least-squares TTF solutions are based on spheres of about 132 000 primary stars uniformly distributed on the sky and simulated observations spanning the entire 5 yr range of Gaia's nominal operational lifetime. The statistical properties of the results compare well with those of GREM. Finally, comparisons to RAMOD@GSR2 solutions confirmed the known lower accuracy of that model and allowed us to establish firm limits on the quality of the linearization point outside of which an iteration for non-linearity is required for its proper convergence. This has proved invaluable as RAMOD@GSR2 is prepared to go into operations on real satellite data.

**Key words.** astrometry – gravitation – methods: data analysis – space vehicles: instruments

## 1. Introduction

The Gaia space satellite, operational since 2014, performs 4π (steradians) absolute astrometry, aiming at the definition of a global astrometric reference frame at visual wavelengths to unprecedented accuracies. Within this context, it will be very difficult to identify possible errors in the measurements or in the data reduction process by using external comparisons of similar accuracies across the sky as they are simply not available. For this reason, beside the main processing pipeline, the Gaia Data Processing and Analysis Consortium (DPAC) has established an astrometric verification unit (AVU) to verify crucial steps in the baseline data processing chain and report on any significant difference (Vecchiato et al. 2012). This step is essential in order to guarantee a high quality catalog to the larger scientific community.

Moreover, the processing of astrometric observations at the
μas (micro-arcsecond) accuracy required by the Gaia mission

demands to take into account relativistic effects on clock syn-
chronization, reference frame transformations as well as on
light propagation. Indeed, the behavior of electromagnetic waves
(EW) in the solar system is highly sensitive to space-time cur-
vature. Most of the available models are based on the solution
of null geodesic equations, including those currently used to
process Gaia observations: Gaia RELativistic Model (GREM;
Klioner 2003) for the main processing pipeline of the Astro-
metric Global Iterative Solution (AGIS; Lindegren et al. 2012)
and the Relativistic Astrometric MODel family (RAMOD;
Crosta et al. 2017, 2015, and references therein) for the AVU
as implemented in the Global Sphere Reconstruction (GSR;
Abbas et al. 2011) subsystem hosted at the Italian Data Process-
ing Center in Torino (DPCT, Messineo et al. 2012).

Nevertheless, other approaches exist, such as the one first
based on the Synge World Function (Le Poncin-Lafitte et al.
2004) and then improved with the use of the time transfer func-
tions (TTF; see Teyssandier & Le Poncin-Lafitte 2008). Simi-
larly to the other methods this approach provides “ready-to-
use” analytical formulas to describe light propagation in the

<footnote n="1">¹ INAF, Astrophysical Observatory of Torino, via Osservatorio 20, 10025 Pino Torinese (Torino), Italy
e-mail: stefano.bertone@aiub.unibe.ch</footnote>

<footnote n="2">² SYRTE, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Universités, UPMC Univ. Paris 06, LNE, 61 avenue de l'Observatoire, 75014 Paris, France</footnote>

<footnote n="3">³ Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan Road, 200030 Shanghai, PR China</footnote>

<footnote n="4">⁴ EURIX SRL, via Giulio Carcano 26, 10153 Torino, Italy</footnote>

<footnote>* Present address: Astronomical Institute, University of Bern, Sidlerstrasse 5, 3011 Bern, Switzerland.</footnote>
field of multiple axisymmetric bodies moving with arbitrary velocities in the post-Minkowskian framework (Hees et al. 2014a) and up to the second and third order for static perturbing bodies (Linet & Teyssandier 2013; Hees et al. 2014b). The equivalence of the TTF approach has been verified analytically with RAMOD in the static case and with GREM in the case of slowly moving monopoles (Bertone et al. 2014). However, equivalence on the analytical side does not guarantee equivalence of the numerical results (Vecchiato et al. 2014).

The task of providing the best possible confidence in the data reduction process is precisely the main responsibility of AVU. To this aim, having the possibility of reducing data with different approaches would surely increase the level of reliability from the point of view of the theoretical understanding of the data.

In this paper, we present results of the first application of the TTF approach to global astrometry in the framework of the Gaia mission, which is preparatory to its successive inclusion in the AVU/GSR operational infrastructure next to the models of the RAMOD family. In Sect. 2 we briefly recall the astrometric problem of processing Gaia observations and how it is addressed in GSR. In Sect. 3 we present our new implementation, GSR-TTF, by first describing how the TTF equations have been integrated in the GSR modeling of Gaia observations at the required accuracy. Then, we provide the equations of the coefficients required for the least-square solution of the sphere. We take advantage of the modularity of the GSR software infrastructure, which allows for an easy plugin of, for example, light propagation equations and reference frames transformations. Finally, we present in Sect. 4 the comparison of our modeling to the GREM implementation in GaiaTools (ter Linden & GaiaTools Committee 2009; Balm & GaiaTools Committee 2014) as well as solutions of a sphere of about 132 000 stars based on DPAC simulated data. Section 5 suggests possible developments and further applications for the results presented in this work and our concluding remarks.

## 2. The astrometric problem in GSR

### 2.1. Astrometric coordinates in Gaia

As shown in Fig. 1, each point of the celestial sphere can be fixed in the reference frame of the Gaia spacecraft by three direction cosines defined as

$$n_{(i)} = (\cos \alpha, \cos \beta, \cos \delta), \qquad (1)$$

with $i = 1, 2, 3$.

From a geometrical point of view, Gaia measures the abscissa $\phi$ and the ordinate $\zeta$, also called along-scan (AL) and across-scan (AC) coordinate, of such a point. In particular, the AL coordinate is the angle $\phi$ between the x-axis of the spacecraft, denoted $E_x$ in Fig. 1, and the projection of the point along the great circle traced by the $E_x$ and $E_y$ axes, which identifies the instantaneous scanning direction of the satellite. The angle $\phi$ is related to the direction cosines $n_{(i)}$ by the following relations

$$\cos \phi = \frac{n_{(1)}}{\sqrt{1 - n_{(3)}^2}}, \qquad (2a)$$

$$\sin \phi = \frac{n_{(2)}}{\sqrt{1 - n_{(3)}^2}}. \qquad (2b)$$

As illustrated in Gaia Collaboration (2016) and Perryman et al. (2001), the Gaia spacecraft has two fields of view (FoV) called

<figure><img src="image_2.png" /><figcaption>Fig. 1. Fundamental angles in the Gaia reference frame $E_\alpha$. The two Fields of View (FoV) directions are indicated by $f_+$ and $f_-$ while the measured abscissa is given by $n_{(i)} = (\alpha, \beta, \delta)$ (from Vecchiato et al. 2012).</figcaption></figure>

$f_+$ and $f_-$ whose pointing directions are separated by a fixed base angle of $\Gamma = 106.5^\circ$ and which are symmetric with respect to the x-axis. Since the angular amplitude of each FoV is about $0.5^\circ$, the abscissae range is fixed in the intervals $53.25\pm0.25$ deg for $f_+$ and $-53.25\pm0.25$ deg for $f_-$. Therefore, one of Eqs. (2) is enough to determine a univocal correspondence between the value of $\phi$ and that of the direction cosines. Gaia's measurement is three times more sensitive in the AL direction than in the across one, so, in principle, the $\cos \phi$ observation equation is sufficient to solve the sphere problem, which is what we have done in this paper; however, since to first order the abscissa measurement is independent from the orientation of the $E_z$ axis, one should also build an observation equation using the AC measurement $\zeta$ to better constrain the attitude when the latter is being reconstructed along with the astrometric parameters. The abscissa is generally expressed as function of the astrometric parameters (namely, the right ascension $\alpha_*$, declination $\delta_*$, parallax $\bar{\omega}_*$, and proper motions $\mu_{\alpha*}$ and $\mu_{\delta*}$) and of the satellite attitude $a_i^{(j)}$, where the index $j$ refers to the time of observation, and $i$ to the number of parameters used to model the attitude. The latter has to be considered unknown since the satellite attitude cannot be determined by other independent measurements at the accuracy required by the mission. For a similar reason Eq. (2a) also depends on a set of instrument parameters $\{c_i\}$ to provide a sort of long-term calibration. Moreover, when working within the parametrized post-Newtonian (PPN) formalism, one should add the parameter $\gamma$ to the unknowns. A better determination of $\gamma$, which measures the amount of curvature induced by the mass-energy on space-time, shall be one of the important scientific contributions of Gaia (Perryman et al. 2001). As consequence, each AL observation can be modeled by a non-linear function of these four classes of unknowns, as

$$\cos \phi \equiv \mathcal{F}(\alpha_*, \delta_*, \bar{\omega}_*, \mu_{\alpha*}, \mu_{\delta*}, a_1^{(j)}, a_2^{(j)}, \dots, c_1, c_2, \dots, \gamma). \qquad (3)$$
## 2.2. The GSR approach to the processing of observations

The Gaia mission will perform several billion observations during its operational years, so that the problem of the so-called sphere reconstruction translates into the solution of a very large system of equations (up to $10^{10} \times 10^8$ in the case of Gaia). Solving such a big system of non-linear equations is not feasible, so the observation Eq. (3) are linearized about a convenient starting point, that is, the current best available estimate of the required unknowns. The problem is then converted in a corresponding system of linear equations

$$ - \sin \phi \Delta\phi = \frac{\partial \mathcal{F}}{\partial \alpha_{*}}|_{\vec{\alpha}_{*}} \delta \alpha_{*} + \frac{\partial \mathcal{F}}{\partial \delta_{*}}|_{\vec{\delta}_{*}} \delta \delta_{*} + \frac{\partial \mathcal{F}}{\partial \omega_{*}}|_{\vec{\omega}_{*}} \delta \omega_{*} \\ + \sum_{ij} \frac{\partial \mathcal{F}}{\partial a_{i}^{(j)}}|_{\vec{a}_{i}^{(j)}} \delta a_{i}^{(j)} + \sum_{i} \frac{\partial \mathcal{F}}{\partial c_{i}}|_{\vec{c}_{i}} \delta c_{i} \\ + \frac{\partial \mathcal{F}}{\partial \gamma}|_{\vec{\gamma}} \delta \gamma + \dots, \qquad (4) $$

where the unknowns are the corrections $\delta x$ to the a-priori values (e.g., from stellar catalogs) while the derivatives of $\mathcal{F}$ are the coefficients of the design matrix. The known-terms are then represented by the left-hand side of Eq. (4) as

$$ \sin \phi \Delta\phi = \sin(\phi_{\text{calc}}) \times (\phi_{\text{obs}} - \phi_{\text{calc}}), \qquad (5) $$

where $\phi_{\text{obs}}$ represents the observed abscissa and formally includes the measurement errors so that it can be written as $\phi_{\text{obs}} = \phi_{\text{true}} + \Delta\phi$, while $\phi_{\text{calc}}$ is the computed value given by $\arccos(\mathcal{F})$ at the starting point of the linearization (generally speaking, the value contained in the astrometric catalog).

The resulting system of equations is quite sparse since each observation refers to a single star among the millions considered in the reconstruction problem (and then only to its astrometric parameters). A similar reasoning is valid for the attitude and calibration parameters, while $\gamma$ is a global parameter in the sense that it appears in each equation of the system. The number of observations being far larger than the number of unknown parameters, the system is over-determined and can be solved by a least-squares procedure. The final goal of the entire procedure is then to get better estimates of all the intervening parameters but for that, we first need to provide an accurate model of the direction cosines of the observation.

The input of the code are packets, each containing observations characterized by: the coordinate time of observation (necessary to retrieve Gaia state vector and the planetary ephemerides at the appropriate epoch), the catalog coordinates of the observed source and all quantities used in the setup of the observation equations for the astrometric problem, which we are going to detail in the following. The values contained in the packet are then updated by the processing of the observation. For each of the processed observations the following steps are performed:

- Load the observations packets and launch the analysis routines;

- Call, for each observation, the routines defining all needed quantities for the computation of the astrometric observable;

- Define all needed vectors (star-observer, perturbing body-observer, ...) and tensors (tetrad components, metric, ...);

- Define $\phi_{\text{obs}}$ and $\phi_{\text{calc}}$, compute the known-terms (5);

- Compute the coefficients of the linearized observation Eq. (4);

- Coefficients and known-terms are associated to an observation and stored in a packet to be then used for the setup of the

observation equation and the astrometric solution in a least-squares sense.

The GSR software has been conceived to be as neutral as possible with respect to the astrometric model. In particular, as long as the latter agrees with some given input and output specifications, one can plug-in any relativistic model, which is seen as a black box by the pipeline. In the current version of GSR (GSR2, see Sect. 4) the direction cosines are provided using a relativistic model from the RAMOD family (more precisely, the version actually implemented in the operational infrastructure at DPCT is RAMOD@GSR2, Vecchiato et al., in prep.), but the modular structure of the software makes it easy to produce a GSR-compatible implementation of the TTF astrometric model, which we denote as GSR-TTF.

## 3. TTF implementation in GSR

The goal of this work is to setup the processing of astrometric observations (eventually made by Gaia) using the TTF formalism. To do it, we implement our model in the GSR software developed at Turin Observatory and we use it to generate a series of simulated observations. The result is a “GSR-TTF plug-in” which enables the pipeline to reduce the sphere with another model, which is actually more accurate than the RAMOD@GSR2 one. The TTF model can, in fact, be implemented at several different level of accuracy. We decided to limit the present one to a $(v/c)^3$ version at the same level of GREM, where $v$ is the velocity of the observer with respect to the Barycentric Celestial Reference System (BCRS). This allows to neglect all the complications coming from higher-orders of $v/c$ and to non-spherically symmetric gravity fields, which are deemed negligible for the specific problem of the global sphere reconstruction. The implementation of the astrometric observable in GSR concerns the development of both sides of Eq. (4), as we detail in this section.

### 3.1. Setup of the observation abscissae

Let us show the procedure followed to build the abscissa $\phi$ using our astrometric model. First, following the definition of astrometric observable in the tetrad formalism (Brumberg 1991) given by Hees et al. (2014b), we define the direction cosines appearing in Eq. (2a), taken at the observation point $x_B$ as

$$ n_B^{(i)} = - \frac{E_{(i)}^0 + E_{(i)}^j \hat{k}_j}{E_{(0)}^0 + E_{(0)}^j \hat{k}_j} \Big|_{x_B}, \qquad (6) $$

where we shall choose the direction triple $\hat{k}_i$ (defining the barycentric direction of light) and the tetrad $E_{(\beta)}^\alpha$ (i.e., the transformation matrix to the reference system comoving with the observer) according to the accuracy required by the Gaia mission.

The relations between the TTF and the wave vectors $k^\mu = dx^\mu/d\lambda$ at reception have been derived in Le Poncin-Lafitte et al. (2004) as

$$ (\hat{k}_i)_B = \left(\frac{k_i}{k_0}\right)_B = -c \frac{\partial \mathcal{T}_e}{\partial x_B^i} = -c \frac{\partial \mathcal{T}_r}{\partial x_B^i} \left[1 - \frac{\partial \mathcal{T}_r}{\partial t_B}\right]^{-1}. \qquad (7) $$

We put $x_A = (ct_A, x_A)$, the event of emission $\mathcal{A}$ and $x_B = (ct_B, x_B)$, the event of reception $\mathcal{B}$ of a light signal. Moreover, we define $\mathcal{T}_e$ and $\mathcal{T}_r$ as two distinct (coordinate) TTF defined as

$$ t_B - t_A = \mathcal{T}_e(t_A, x_A, x_B) = \mathcal{T}_r(t_B, x_A, x_B), \qquad (8) $$
where $\mathcal{T}_e$ and $\mathcal{T}_r$ are evaluated at the event of emission $\mathcal{A}$ and at the event of reception $\mathcal{B}$ respectively. For astrometric applications, we are only interested in $\mathcal{T}_r$.

For applications in the solar system, we can work in the weak field approximation so that we can write

$$g_{\mu\nu} = \eta_{\mu\nu} + h_{\mu\nu}, \qquad (9)$$

with $\eta_{\mu\nu} = \text{diag}(-1, +1, +1, +1)$ a Minkowskian background and $h_{\mu\nu}$ a small perturbation. Moreover, we shall consider only weak gravitational fields generated by self-gravitating extended bodies within the slow-motion, post-Newtonian approximation. So, we assume that the potentials $h_{\mu\nu}$ may be expanded as

$$ \begin{aligned} h_{00} &= h_{00}^{(2)} + O\left(\frac{1}{c^4}\right), \\ h_{0i} &= h_{0i}^{(3)} + O\left(\frac{1}{c^5}\right), \\ h_{ij} &= h_{ij}^{(2)} + O\left(\frac{1}{c^4}\right), \end{aligned} \qquad (10) $$

where $h_{\alpha\beta}^{(2n)} \sim (U/c^2)^n$, $h_{\alpha\beta}^{(2n+1)} \sim -(U/c^2)^n (v/c)$ and $U$ is the Newtonian potential.

Under these hypotheses, the expression of the TTF $\mathcal{T}_{e/r}$ in the PPN approximation has been given in Linet & Teyssandier (2002) and in Teyssandier et al. (2008) as

$$\mathcal{T}_r(\mathbf{x}_A, t_B, \mathbf{x}_B) = \frac{R_{AB}}{c} + \frac{1}{c} \Delta_r(\mathbf{x}_A, t_B, \mathbf{x}_B) + \mathcal{O}(c^{-5}), \qquad (11)$$

where $\Delta_r$ is defined as

$$\Delta_r = \frac{1}{2} R_{AB} \int_0^1 [h_{00}^{(2)} + 2N_{AB}^i h_{0i}^{(3)} + N_{AB}^j N_{AB}^i h_{ij}^{(2)}]_{z^a(\lambda)} d\lambda, \qquad (12)$$

and the integrals are taken along the Minkowskian paths $z^a(\lambda) = (x_B^0 - \lambda R_{AB}, x_B^i - \lambda R_{AB}^i)$ and we define $R_{AB}^i = x_B^i - x_A^i$, $R_{AB} = |R_{AB}^i|$ and $N^i = \frac{R_{AB}^i}{R_{AB}}$.

The derivatives needed in the definition of the wave vectors in Eq. (7) can then be computed as (see, e.g., Bertone et al. 2014; Hees et al. 2014b,a)

$$\frac{\partial \mathcal{T}_r}{\partial x_{A/B}^i} = \mp \frac{N_{AB}^i}{c} + \frac{1}{c} \frac{\partial \Delta_r}{\partial x_{A/B}^i} + \mathcal{O}(c^{-5}), \qquad (13a)$$

and

$$\frac{\partial \mathcal{T}_r}{\partial t_B} = \frac{1}{c} \frac{\partial \Delta_r}{\partial t_B} + \mathcal{O}(c^{-5}), \qquad (13b)$$

where

$$\frac{\partial \Delta_r}{\partial x_B^i} = -\frac{1}{2} \int_0^1 [R_{AB}^i \lambda m_{,0} - R_{AB}(1-\lambda)m_{,i} - \tilde{h}_i]_{z^a(\lambda)} d\lambda, \qquad (14a)$$

and

$$\frac{\partial \Delta_r}{\partial t_B} = c \frac{R_{AB}}{2} \int_0^1 [m_{,0}]_{z^a(\lambda)} d\lambda. \qquad (14b)$$

Moreover, we defined

$$ \begin{aligned} m_{,\alpha} &\equiv h_{00,\alpha}^{(2)} + 2N_{AB}^k h_{0k,\alpha}^{(3)} + N_{AB}^j N_{AB}^k h_{jk,\alpha}^{(2)} \\ \tilde{h}_i &\equiv N_{AB}^i h_{00}^{(2)} - N_{AB}^i N_{AB}^j N_{AB}^k h_{jk}^{(2)} + 2h_{0i}^{(3)} + 2N_{AB}^j h_{ij}^{(2)}. \end{aligned} $$

For most observations, at Gaia accuracy we shall only consider the PN gravitational potentials of all solar system bodies as point masses. Quadrupole terms might be significant for a limited number of observations grazing giant planets (e.g., closer than 152" from Jupiter's limb at 1 µas accuracy Klioner 2003; Crosta & Mignard 2006), hence we are including them in our implementation.

Also, it has been shown (Vecchiato et al. 2003) that thanks to the satellite's high measurement accuracy, the Gaia sphere reconstruction is sensitive to variations of the $\gamma$ parameter of the PPN formalism up to the $10^{-6}-10^{-7}$ level. This processing can thus be used as a tool to improve our knowledge of such parameter, which is currently set at the $10^{-5}$ level (Will 2014). With this application in mind $\gamma$ can be straightforwardly introduced by using the PPN expression of the above metric coefficients $h_{\alpha\beta}$, where

$$h_{00}^{(2)} = \sum_{P} \frac{2w_{P}(t)}{c^2}, \qquad (15)$$

$$h_{0i}^{(3)} = -(1 + \gamma) \frac{2}{c^3} \sum_{P} w_{P}(t) v_{P}(t), \qquad (16)$$

$$h_{ij}^{(2)} = \delta_{ij} \gamma h_{00}^{(2)}, \qquad (17)$$

and $w_P(t)$ is the time-dependent gravitational potential of perturbing body $P$. This time-dependence would in principle require either to express $w$ as a retarded potential, or to compute a full integral of the geodesic. However, it has been shown in Klioner (2003) and confirmed by our analysis in Bertone et al. (2014), that at the Gaia accuracy the effects of a time-dependent perturbation due to a moving body can be taken into account simply by computing the positions of the gravitating bodies at the retarded moment of time

$$t_C = t_B - c^{-1} |\mathbf{x}_B - \mathbf{x}_P(t_B)| + \mathcal{O}(c^{-2}), \qquad (18)$$

where $t_B$ is the coordinate time of observation and $\mathbf{x}_P(t_B)$ is the position of body $P$ at $t_B$, and that the “gravitomagnetic” contributions proportional to $v_P(t)$ can be neglected.

Thus, using such constant value $\mathbf{x}_P = \mathbf{x}_P(t_C)$ in our computations yields to the following expression of the PN terms of the metric tensor:

$$h_{00} = \sum_{P} \frac{2w_{P}(t_C)}{c^2}, \qquad (19a)$$

$$h_{0i} = 0, \qquad (19b)$$

$$h_{ij} = \delta_{ij} \gamma h_{00}. \qquad (19c)$$

Let us assume that the smallest sphere containing the body has a radius equal to the equatorial radius $r_e$ of the body and that the segment joining $\mathbf{x}_A$ and $\mathbf{x}_B$ is outside this sphere. Also, let us restrict to the case of axisymmetric perturbing bodies with time-independent mass-multipoles. At any point $\mathbf{x}$ such that $r \ge r_e$, the gravitational potentials $w$ and $\mathbf{w}$ are then given (for each body $P$) by the multipole expansion (Linet & Teyssandier 2002)

$$w_P = \frac{GM}{r_P} \left[ 1 - \sum_{n=2}^{\infty} J_n \left( \frac{r_e}{r_P} \right)^n P_n \left( \frac{\mathbf{k} \cdot \mathbf{x}}{r_P} \right) \right] \quad \text{and} \quad \mathbf{w}_P = 0, \qquad (20)$$

where

$$\mathbf{r}_P = \mathbf{x}_B - \mathbf{x}_P(t_C) - \lambda R_{AB} \mathbf{N}_{AB},$$

$\mathbf{k}$ denotes the unit vector along the $x^3$-axis, the $P_n$ are the Legendre polynomials, $M$ is the mass of the body and the coefficients $J_n$ are its zonal spherical harmonics coefficients.
Under these assumptions, the monopolar and quadrupolar terms of the direction triple to be used in Eq. (6) are (Le Poncin-Lafitte & Teyssandier 2008; Bertone et al. 2013)

$$ (\hat{k}_i^B)(\mathbf{x}_A, \mathbf{x}_B) = N_{AB}^i - (\gamma + 1) \sum_P \frac{GM_P}{c^2 R_{PB}} \frac{1}{1 + N_{PA} \cdot N_{PB}} \times \left[ \frac{R_{AB}}{R_{PA}} N_{PB} - \left(1 + \frac{R_{PB}}{R_{PA}}\right) N_{AB} \right], \quad (21a) $$

$$ \begin{aligned} (\hat{k}_i^B)_{J_2}(\mathbf{x}_A, \mathbf{x}_B) = & (\gamma + 1) \frac{GM}{c^2} J_2 r_c^2 \left\{ -[\mathbf{k} \cdot (N_{PA} + N_{PB})]^2 \right. \\ & \times \left[ \frac{N_{PB} - N_{AB}}{(R_{PA} + R_{PB} - R_{AB})^3} - \frac{N_{PB} + N_{AB}}{(R_{PA} + R_{PB} + R_{AB})^3} \right] \\ & + \frac{1}{2} \left[ \frac{1 - (\mathbf{k} \cdot N_{PA})^2}{R_{PA}} + \frac{1 - (\mathbf{k} \cdot N_{PB})^2}{R_{PB}} \right] \\ & \times \left[ \frac{N_{PB} - N_{AB}}{(R_{PA} + R_{PB} - R_{AB})^2} - \frac{N_{PB} + N_{AB}}{(R_{PA} + R_{PB} + R_{AB})^2} \right] \\ & + \frac{1}{R_{PB}^3} \frac{(R_{PA} + R_{PB})R_{AB}}{R_{PA}^2} \frac{\mathbf{k} \cdot (N_{PA} + N_{PB})}{(1 + N_{PA} \cdot N_{PB})^2} \\ & \times \left[ \mathbf{k} - (\mathbf{k} \cdot N_{PB})N_{PB} \right] \\ & \left. + \frac{1}{2R_{PB}^3} \frac{R_{AB}}{R_{PA}} \frac{2(\mathbf{k} \cdot N_{PB})\mathbf{k} + [1 - 3(\mathbf{k} \cdot N_{PB})^2] N_{PB}}{1 + N_{PA} \cdot N_{PB}} \right\} \end{aligned} \quad (21b) $$

respectively, and where we used

$$ \mathbf{R}_{PX} = \mathbf{x}_X - \mathbf{x}_P(t_c), \quad (22a) $$

$$ R_{PX} = |\mathbf{R}_{PX}|, \quad N_{PX} = \frac{\mathbf{R}_{PX}}{R_{PX}}. \quad (22b) $$

An expression for the multipolar terms of higher order is given in Le Poncin-Lafitte & Teyssandier (2008) but it is not relevant for this study.

The direction triple $(\hat{k}_i^B)$, defined as the sum of Eqs. (21a) and (21b), gives the direction of the light ray at reception in the local barycentric frame. Following Eq. (6), we hence need to project this vector in the observer reference frame. Concerning this point, we have used the tetrad defined in Crosta & Vecchiato (2010) and already used in GSR, evaluated with the same metric (Eq. 19)).

The satellite reference frame defining the transformation $E_{\hat{a}}^B$ is obtained by successive transformations of the local BCRS tetrad (Bini et al. 2003)

$$ \begin{aligned} \lambda_{\hat{a}}^{\alpha} &= h_{0\alpha} \delta_{0}^{\alpha} + \left(1 - \frac{h_{0\alpha}}{2}\right) \delta_{a}^{\alpha} + O(h^2), \\ &= h_{0\alpha} \delta_{0}^{\alpha} + \left(1 - \gamma \frac{h_{00}}{2}\right) \delta_{a}^{\alpha} + O(h^2), \end{aligned} \quad (23) $$

namely the tetrad obtained from the BCRS by shifting the origin to the instantaneous barycenter of the satellite. This tetrad, and in particular the triad of four-vectors $\lambda_{\hat{a}}$ defined as the spatial part of $\lambda_{\hat{a}}^{\alpha}$, is “boosted” to the satellite rest-frame by means of an instantaneous Lorentz transformation identified by the four-velocity of the observer $u_s = \frac{(-c, v)}{\sqrt{1-v^2/c^2}}$ with respect to the local BCRS. The boosted tetrad $\tilde{\lambda}_{\hat{a}}$, namely (Crosta & Vecchiato 2010)

$$ \tilde{\lambda}_{\hat{a}}^{\alpha} = \lambda_{\hat{a}}^{\alpha} + \left(1 + 3\gamma \frac{h_{00}}{2} + \frac{1}{2} \frac{v^2}{c^2}\right) \delta_{0}^{\alpha} \frac{v^a}{c} + \frac{1}{2} \frac{v^i}{c} \delta_{i}^{\alpha} \frac{v^a}{c}, \quad (24) $$

obtained in this way represents a reference system whose origin is comoving with the barycenter of the satellite and whose spatial axes are kinematically non-rotating with respect to the BCRS. The Gaia attitude frame is finally obtained with a rotation of these spatial axes. Since the transformation is Euclidean, the rotation matrix A is exactly the attitude matrix of the satellite.

The triad resulting from these transformations, detailed in Crosta & Vecchiato (2010), establishes the Gaia relativistic attitude triad and is related to the euclidean attitude parameters by

$$ \begin{pmatrix} E_1^1 \\ E_2^1 \\ E_3^1 \end{pmatrix} = A \begin{pmatrix} \tilde{\lambda}_1^1 \\ \tilde{\lambda}_2^1 \\ \tilde{\lambda}_3^1 \end{pmatrix}. \quad (25) $$

The temporal axis of the tetrad has to be transformed alike to obtain a complete reference system attached to the observer. This result is naturally achieved by writing the transformation between the barycentric coordinate time and the observer's proper time.

Finally, introducing Eq. (6) into Eq. (2a) completes the implementation in the observation abscissa of GSR-TTF.

## 3.2. Setup of the astrometric coefficients

In order to setup the observation equation for the improvement of the a priori catalog values by a least-square procedure, we define the partial derivatives of Eq. (3) w.r.t. the astrometric parameters as in Eq. (4). The coefficients can be computed analytically since the function $F = \cos \phi_{calc}$, defined by Eq. (2a), is known. By defining $n_{(i)} = \cos \psi_{(i)}$, the coefficients of Eq. (2a) can be computed analytically as

$$ \frac{\partial \cos \phi}{\partial v} = \frac{\partial_v(\cos \psi_{(1)})}{\sqrt{1 - \cos^2 \psi_{(3)}}} + \frac{\cos \psi_{(1)} \cos \psi_{(3)} \partial_v(\cos \psi_{(3)})}{(1 - \cos^2 \psi_{(3)})^{3/2}}, \quad (26) $$

where $v$ are the so called astrometric parameters characterizing the coordinates of a star at a moment in time: the parallax $\sigma_*$, right ascension $\alpha_*$, declination $\delta_*$, and the proper motions $\mu_{\alpha_*}$ and $\mu_{\delta_*}$ in perpendicular directions. The partial derivatives appearing therein are given by

$$ \begin{aligned} \partial_v \cos \psi_{(i)} &\equiv \frac{\partial \cos \psi_{(i)}}{\partial v} \\ &= - \frac{E_{(i)}^j \partial_v \hat{k}_j (E_{(0)}^0 + E_{(0)}^l \hat{k}_l) - (E_{(i)}^0 + E_{(i)}^j \hat{k}_j) E_{(0)}^l \partial_v \hat{k}_l}{(E_{(0)}^0 + E_{(0)}^l \hat{k}_l)^2}, \end{aligned} \quad (27) $$

with

$$ \begin{aligned} \frac{\partial \hat{k}_i}{\partial v} = & \partial_v N_{AB}^i - (\gamma + 1) \sum_P \frac{GM_P}{c^2 R_{PB}} \frac{1}{1 + N_{PA} \cdot N_{PB}} \\ & \times \left\{ - \frac{\partial_v N_{PA} \cdot N_{PB}}{1 + N_{PA} \cdot N_{PB}} \times \left[ N_{AB}^i \frac{R_{AB}}{R_{PA}} - N_{AB}^i \left(1 + \frac{R_{PB}}{R_{PA}}\right) \right] \right. \\ & \left. + \frac{N_{PB}^i}{R_{PA}^2} \left(R_{PA} \partial_v R_{AB} - R_{AB} \partial_v R_{PA}\right) \right. \\ & \left. - \partial_v N_{AB}^i \left(1 + \frac{R_{PB}}{R_{PA}}\right) + N_{AB}^i R_{PB} \frac{\partial_v R_{PA}}{R_{PA}^2} \right\}. \end{aligned} \quad (28) $$
Moreover, we used the definitions (22) along with

$$ \partial_v N_{AB}^i = - \frac{(N_{AB} \times (\partial_v \mathbf{x}_A \times N_{AB}))^i}{R_{AB}}, \qquad (29a) $$

$$ \partial_v R_{AB} = - \frac{R_{AB} \cdot \partial_v \mathbf{x}_A}{R_{AB}}, \qquad (29b) $$

$$ \partial_v R_{PA} = \frac{R_{PA} \cdot \partial_v \mathbf{x}_A}{R_{PA}}, \qquad (29c) $$

$$ \partial_v N_{PA} = \frac{1}{R_{PA}} [\partial_v \mathbf{x}_A - N_{PA}(N_{PA} \cdot \partial_v \mathbf{x}_A)], \qquad (29d) $$

where $\times$ indicates a cross-product. The expression in terms of the astrometric parameters can be easily obtained by considering that the coordinate position of the source $\mathbf{x}_A(t)$ can be treated like a Euclidean vector, so that $\mathbf{x}_A(t) = (r_A \cos \alpha(t) \cos \delta(t), r_A \sin \alpha(t) \cos \delta(t), r_A \sin \delta(t))$, where $r_A = |\mathbf{x}_A|$ is the barycentric distance of the star, $\varpi = 1 UA/r_A$. Finally, the angular coordinates $\alpha(t)$ and $\delta(t)$ at a given time $t$ can be expressed as functions of the catalog positions and proper motion with, for example, $\alpha(t) = \alpha_{t0} + \mu_{\alpha,t}dt$, $dt$ being the interval of time between $t$ and the catalog reference time $t_0$. We can thus write

$$ \partial_{\varpi_*} x_A^i = - \frac{1 UA}{\varpi_A} \times n_A^i, \qquad (30a) $$

$$ \partial_{\alpha_*} x_A^i = (-r_A \sin \alpha \cos \delta, r_A \cos \alpha \cos \delta, 0), \qquad (30b) $$

$$ \partial_{\delta_*} x_A^i = (-r_A \cos \alpha \sin \delta, -r_A \sin \alpha \sin \delta, r_A \cos \delta), \qquad (30c) $$

$$ \partial_{\mu_{\alpha_*}} x_A^i = dt \times \partial_{\alpha_{\alpha_*}} x_A^i, \qquad (30d) $$

$$ \partial_{\mu_{\delta_*}} x_A^i = dt \times \partial_{\delta_{\alpha_*}} x_A^i. \qquad (30e) $$

The system is then solved following a classical least-squares method.

### 3.3. Setup of the global PPN parameter $\gamma$

Following a similar approach, we computed and implemented the partial derivatives of the astrometric observable with respect to the global PPN parameter $\gamma$, which can be used to test General Relativity using light propagation measurements.

Considering Eq. (26), we have that both the tetrad $E_{\beta}^{\alpha}$ and the direction triple $\hat{k}_B$ depend on $\gamma$ through the metric (19). Hence we get

$$ \partial_{\gamma} \cos \psi_{(i)} \equiv \frac{\partial \cos \psi_{(i)}}{\partial \gamma} = - \frac{(E_{(i)}^j \partial_{\gamma} \hat{k}_j + \partial_{\gamma} E_{(i)}^j \hat{k}_j)(E_{(0)}^0 + E_{(0)}^l \hat{k}_l) - (E_{(i)}^0 + E_{(i)}^j \hat{k}_j)E_{(0)}^l \partial_{\gamma} \hat{k}_l}{(E_{(0)}^0 + E_{(0)}^l \hat{k}_l)^2}, \qquad (31) $$

and, with Eqs. (21), (24) and (25),

$$ \frac{\partial \hat{k}_i}{\partial \gamma} = \frac{(\hat{k}_i^B) + (\hat{k}_i^B)_2}{(\gamma + 1)} \quad \text{and} \quad \partial_{\gamma} E_{(i)}^0 = \partial_{\gamma} E_{(0)}^{\alpha} = 0. \qquad (32) $$

## 4. GSR-TTF: results with Gaia simulated observations

Our implementation is tested by using observations of the Gaia satellite as simulated by DPAC using the GREM model. In particular, the dataset employed is RDS-7-F (Luri et al. 2010), a full 5 yr simulation of Gaia observations of two million stars,

approximately one million of which are Primary astrometric sources in the sense discussed in Lindegren et al. (2012). Besides the noise-free measured GREM abscissa $\eta$ (see Eq. (33) below), the data set provides: stellar coordinates, observation epoch, instantaneous satellite attitude, and the ephemerides of the solar system. These data allow us to compute the abscissa $\phi$ from a given model. As for the observation error, this can be customized. In our experiments all stars have Gaia magnitude $G = 14$, and at such brightness the expected single observation error is $\sigma_0 = 125$ µas.

The TTF model is therefore tested in two steps: by comparing the abscissae (2a) computed by TTF and GREM, and by trying some reconstructions of the global astrometric sphere. The latter had to be restricted to relatively small subsamples of stars and observations in RDS-7-F. This is because only limited resources are allocated to the test & development section of the GSR system at DPCT, the environment for running experiments with simulated data.

### 4.1. Numerical comparison of computed observables

As anticipated above, we first focused on the comparison of the abscissa $\phi$. Using the implementation described in Sect. 3.1, we applied Eq. (2a) over a set of observations from a chosen day of the given dataset using the TTF model, thus obtaining $\phi_{TTF}$. The same data are used with the tools provided by the DPAC code, which implements the GREM model, to obtain $\phi_{GREM}$. Finally, the residuals $\phi_{TTF} - \phi_{GREM}$ were computed.

The results are illustrated in Fig. 2, where the models are compared. The left axis marks the difference in $\mu$as between the two models while the right axis is the angular distance in degrees between a given planet and the observation, namely its elongation. Lines representing Jupiter, Saturn, and the Earth are all plotted. In particular, the periodic oscillation of the distance planet-observation illustrated in the plots is due to the Gaia scanning law (de Bruijne et al. 2010) setting a rotation period of approximately 6 h.

As expected from the analytical comparisons in Bertone et al. (2013) and Bertone et al. (2014), the plot shows that the differences between the abscissae computed with the two models are generally well below $\pm 1$ µas. The remaining sub-µas signal can be attributed to a different accounting of the retarded times in Eq. (18) or to differences in the representation of satellite attitude.

The final accuracy of a global sphere reconstruction depends on the statistical properties of the model accuracy, in the sense that possible inaccuracies affecting some observations because of its specific geometric arrangement are likely to get “blended” and “smoothed out” by other, more accurate observations featuring a more benign configuration. It is therefore appropriate to give also a statistical assessment of these differences. Figure 3 shows the distribution of abscissae differences for observations taken over a period of 40 days. The vast majority of the differences are in the range $\pm 0.2$ µas.

### 4.2. Application to the reconstruction of the celestial sphere

The main goal of the Gaia mission is to improve the quality of the stellar catalogs for coordinates, parallaxes and annual proper motions of the observed stars. This is done by evaluating the differences between the high-precision angular measurements made by the satellite and their analytical modeling as functions of the astrometric stellar parameters. As described in Sect. 2.2,
<figure><img src="image_3.png" /><figcaption>Fig. 2. Difference between the abscissae resulting from the GSR-TTF implementation and from GREM implementation in the “GaiaTools”. The left axis indicates the difference in µas between the two models – represented by the red plot; the right axis indicates the distance in degrees between a given planet and the observation – the blue, yellow and green plot representing Jupiter, Saturn and the Earth, respectively.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Abscissae differences between the GSR-TTF and GREM implementations. The vast majority of the observations show differences in the interval of ±0.2 µas, well below the accuracy required by Gaia.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Frequency of observations as function of celestial coordinates due to Gaia scanning law (blue ≈50 – yellow ≈500).</figcaption></figure>

this is essentially a minimization problem obtained by solving a large and sparse system of linearized observation equations in the least-squares sense. Indeed, the problem is largely over-determined (Gaia will provide around 700 observations for each star, see e.g. Mignard et al. 2008, or Fig. 4) which justifies the least-squares solution of the equation system.

Our implementation is then tested by using observations of the Gaia satellite simulated with the GREM model, as a replacement of the real Gaia observations, and the TTF model. In particular, for each observation it is possible to compute the true AL field angle η. Then, by definition, the abscissa of Eq. (2a) can be expressed as

$$ \phi = f \frac{\Gamma}{2} + \eta, \qquad (33) $$

where Γ, as shown in Fig. 1, is the basic angle, that is the angle between the axes of the preceding and following FoV, and f = ±1 with sign determined by the FoV in which the star is observed. In particular f takes the minus sign in the case of f-, and the positive one for f+.

The abscissae defined in this way are the values of $\phi_{\text{obs}}$ in Eq. (5), while $\phi_{\text{calc}}$ is given by the TTF model. This allows the computation of the known term of the linearized observation equation. Such an equation is completely set up by computing the first order derivatives with respect to the unknown parameters, as described in Sect. 3.2, for which we still need to define the linearization point. Finally, measurement errors must also be introduced in each observation equation.

Specifically, the observed abscissa $\phi_{\text{obs}}$ produced from the dataset should represent the value of a measurement. However, the simulated dataset provides us with the true values, in the sense specified above, that is, the GREM model is supposed to give the correct physical representation of light propagation. Moreover, the coefficients of the linearized equation are computed at specific values, which are our best estimation of the unknowns. Once again, the input catalog provides the true values, while in a more realistic situation one should take into account that our initial guesses will be perturbed by some catalog error. The worse these errors, the larger are the linearization errors which will eventually affect the solution. The above discussion, therefore, makes it clear why in our tests we follow a procedure in three steps:

1. we use the true catalog values to simulate both Gaia “measurements” ($\phi_{\text{obs}}$) and compute the right-hand-side of Eq. (4). This is equivalent to assume no measurement or catalog error on our a-priori information and allows us to test the reliability of our model (see Col. 0101 in Table 1);

2. we apply a random Gaussian error ($\mu = 0, \sigma_o = 125$ µas) to $\phi_{\text{obs}}$, but we use the true values for $\phi_{\text{calc}}$ and for the coefficients. This simulates the impact of the measurement error on the quality of the updated catalog (see Col. 0102 in Table 1), but minimizes the linearization errors on the final solution;

3. in addition to the above measurement error, we apply a Gaussian noise of $\sigma = 20$ mas to the catalog values. This introduces a contribution to the final solution due to the linearization errors of Eq. (4) (see Col. 0301 in Table 1).

We emphasize that these kind of runs are only testing the applicability of the TTF model in the case of the estimation of the astrometric parameters. The actual case of a Gaia-like mission, the attitude parameters have to be estimated as well within the global sphere reconstruction. This is because neither the usual instrumentation for the on-board attitude reconstitution nor a short-term astrometric determination can provide a sufficiently good reconstruction, given the measurement accuracy. In this proof-of-principle, however, the main goal is to test the applicability of TTF; therefore, as a first step, it is required that we compare its
<table><caption>**Table 1.** Results with GSR-TTF applied to 5 yr of observations of 132 K stars homogeneously distributed on the sky (in µas and µas/yr).</caption><thead><tr><th rowspan="2"></th><th colspan="2">0101</th><th colspan="2">0102</th><th colspan="2">0301</th></tr><tr><th>x&#772;</th><th>σ<sub>x</sub></th><th>x&#772;</th><th>σ<sub>x</sub></th><th>x&#772;</th><th>σ<sub>x</sub></th></tr></thead><tbody><tr><td>ω</td><td>0.000</td><td>0.040</td><td>-0.059</td><td>11.520</td><td>-0.033</td><td>11.480</td></tr><tr><td>α cos δ</td><td>-0.002</td><td>0.090</td><td>0.003</td><td>9.180</td><td>-0.64</td><td>9.420</td></tr><tr><td>δ</td><td>-0.009</td><td>0.100</td><td>-0.001</td><td>7.960</td><td>0.000</td><td>8.000</td></tr><tr><td>μ<sub>α</sub> cos δ</td><td>0.001</td><td>0.020</td><td>-0.011</td><td>6.440</td><td>0.021</td><td>6.430</td></tr><tr><td>μ<sub>δ</sub></td><td>-0.001</td><td>0.020</td><td>-0.014</td><td>5.680</td><td>0.010</td><td>5.660</td></tr></tbody></table>

accuracy with that of the Gaia baseline model GREM treating the attitude as perfectly known.

We therefore applied this procedure to a subset of RDS-7-F, namely, to 132 000 stars homogeneously distributed on the sky for which the five years of the planned mission provide ~92 400 000 individual observations.

The final goal of the procedure is to retrieve as much as possible the catalog values for the estimated parameters of the 132 000 processed stars. We evaluated the quality of the estimation simply by computing the residuals of the updated minus true astrometric parameters (separately for each of the five) and afterwards by calculating their average and σ.

If the two models gave exactly the same value of the abscissae, and allowing an infinite machine precision, in the 0101 case we would expect exactly a zero solution. Since computers have a finite computing precision, and numerical errors tend to cumulate, even if the first condition would be met one should expect deviations from this solution larger than the machine precision, which in our case is about $10^{-16}$. Such numbers have to be interpreted as radians, which implies that deviations cannot be smaller than $10^{-4}$ µas. Moreover, the different numerical predictions of the two models, if present, can show up both as random errors and systematic biases in the TTF reconstruction, namely in non-negligible σ and averages respectively. Within Gaia's astrometric tolerances, the two models are equivalent if in such test one obtains sub-µas values for these two quantities.

In the 0102 case we introduced a purely Gaussian measurement error, which should decrease with $\sqrt{n_{obs}}$, where $n_{obs}$ is the number of observations of a particular star. In fact, when the attitude is assumed perfectly known, along with the instrumental and global parameters, each star is completely independent from all the others, and the complete system is indeed equivalent to a set of independent systems, one for each star. Figure 5 confirms the expected behavior of the residuals for each parameter and each star. Considering that Gaia is going to measure each star about 700 times, on average, and that there are five astrometric parameters, one can roughly estimate that

$$ \sigma_f \approx \sqrt{\frac{5}{700}} \sigma_o \approx 10 \mu\text{as.} \quad (34) $$

Finally, in the 0301 case, we have to consider the effect of the linearization errors. It is well known that, as long as the initial guess provided by the catalog is sufficiently close to the true values, this test should give the same result of the 0102 case. It is implicit that the expression “the same result” must be intended in the numerical sense, that is “the same with respect to the required accuracy”, which in the Gaia case means to the sub-µas level.

Table 1 shows the results of these tests, while Fig. 6 shows the post-fit residuals distribution for all astrometric parameters.

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Astrometric post-fit residuals for a set of 132 000 stars as function of the number of observations for each star (red: ω, blue: α cos δ, green: δ, yellow: μ<sub>α</sub> cos δ, purple: μ<sub>δ</sub>).</figcaption></figure>

<table><caption>**Table 2.** Results with RAMOD@GSR2 applied to 5 yr of observations of 132 K stars homogeneously distributed on the sky (in µas and µas/yr).</caption><thead><tr><th rowspan="2"></th><th colspan="2">0101</th><th colspan="2">0102</th><th colspan="2">0301</th></tr><tr><th>x&#772;</th><th>σ<sub>x</sub></th><th>x&#772;</th><th>σ<sub>x</sub></th><th>x&#772;</th><th>σ<sub>x</sub></th></tr></thead><tbody><tr><td>ω</td><td>-0.123</td><td>0.270</td><td>-0.130</td><td>12.02</td><td>-0.032</td><td>12.06</td></tr><tr><td>α cos δ</td><td>-0.011</td><td>0.240</td><td>-0.012</td><td>9.610</td><td>0.519</td><td>16.950</td></tr><tr><td>δ</td><td>-0.050</td><td>0.310</td><td>-0.053</td><td>8.270</td><td>-2.943</td><td>12.600</td></tr><tr><td>μ<sub>α</sub> cos δ</td><td>-0.001</td><td>0.170</td><td>0.024</td><td>6.710</td><td>0.753</td><td>18.260</td></tr><tr><td>μ<sub>δ</sub></td><td>-0.006</td><td>0.150</td><td>0.006</td><td>5.860</td><td>1.120</td><td>10.690</td></tr></tbody></table>

Therefore, for what concerns the astrometric unknowns, we can conclude that the GREM and TTF models are equivalent. Moreover, TTF is able to recover these parameters down to the expected level of accuracy. Finally, starting from a reasonable catalog error the same astrometric accuracy can be recovered at the sub-µas level.

Table 2 reports the results of analogous runs on the same simulated datasets, but this time obtained using the RAMOD@GSR2 model, the one operational in the current version of the GSR system (GSR2) at DPCT. Runs 0101 and 0102 compare reasonably well with the corresponding ones of TTF, although the sub-µas level of the average residuals already reveals the expected intrinsic lower accuracy of the RAMOD@GSR2 model. The discrepancies appearing in run
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Histogram of post-fit residuals for the astrometric parameters affected by measurement and catalog errors (see 0301 Col. in Table 1).</figcaption></figure>

0301, the one reflecting the real case scenario, are significant at the $\mu$as level and we carefully addressed their origin.

The ultimate reason of this discrepancy is actually non-linearity effects which can be recovered by means of a so-called iteration for non-linearity, as described in detail in a forthcoming paper (Vecchiato et al., in prep.).

### 4.3. Sensitivity analysis for the global PPN parameter $\gamma$

The problem in a real Gaia-like sphere reconstruction is more complex. As already mentioned above, it includes the simultaneous estimation of satellite attitude. Besides attitude, a determination of some instrument parameters is also necessary to obtain an unbiased solution. These include, for example, long-term variations of some instrument parameters, or their fine calibration to a level which cannot be reached by daily processes. The latter, however, are unrelated to the relativistic astrometric model, and are therefore out of the scope of this paper.

On the other hand, the $\gamma$ parameter of the parametrized post-Newtonian formulation is indeed part of the relativistic model. It belongs to a class of unknowns of the sphere reconstruction problem which is called global because it enters in every single observation equation. In this case, the least-squares estimation cannot be reduced to the solution of $N$ independent $5 \times 5$ systems of equations solving for the corrections to the catalog values of the five astrometric parameters of each of the $N$ stars. Since the normal matrix is now block diagonal with 1-column border, it must be solved using a more complex algorithm. It has been shown (see, e.g., Bombrun et al. 2010) that the problem can be tackled by first solving the reduced normal equations for the $\delta\gamma$ unknown, and then forwarding that solution into the $N$ diagonal blocks to solve for the astrometric parameters. Since

for each simulated sphere solution we only have one estimate of gamma, we can study its statistical properties by using a Monte-Carlo (MC) simulation.

We performed a MC experiment by running 100 simulations for the simultaneous solution of both PPN-$\gamma$ and the five astrometric parameters. For this experiment the stellar sample (and the corresponding ensemble of satellite simulated observations) extracted from RDS-7-F had to be set to 12 000 primary sources because of limited resources allocated by the GSR2 system to the development (and test) infrastructure. Every MC solution provides an estimation $\gamma_e$ of the PPN parameter<fnref n="1" />, which can be considered a random extraction from the Gaussian sample of $\gamma$ estimations. The statistical accuracy in the determination of this parameter is thus given by the standard deviation of the distribution of the residuals with respect to the true value $\gamma_t$ (=1, the GR value), namely of the variate $\Delta\gamma = \gamma_t - \gamma_e$.

The catalog errors for all runs were set:

- for the astrometric quantities, to the same values as for the 0301 test in Sect. 4.2;

- for the PPN-$\gamma$, to $10^{-3}$.

Following the analysis presented in Vecchiato et al. (2003), it is possible to give an order-of-magnitude estimation of the expected accuracy on a single $\gamma$ measurement by applying error propagation to the formula of gravitational light deflection (e.g., Misner et al. 1973).

In the case of a Gaussian measurement error of 125 $\mu$as, as in our simulations, and a most favorable observational configuration at 45° from the Sun, the propagated error on $\gamma$ is $\sigma_{\gamma} \approx 2.5 \times 10^{-2}$. Assuming that each observation only contributes to the estimation of $\gamma$, and considering that the 12 000 primaries are observed more than $8 \times 10^6$ times in 5 yr, the improvement factor of $2.8 \times 10^3$ yields a best-accuracy value of $\sigma_{\gamma} \approx 1 \times 10^{-5}$. The outcome of the MC simulation is a distribution of residuals $\Delta\gamma = \gamma_t - \gamma_e$ centered at $-6.1 \times 10^{-5}$ and with a standard deviation $\sigma_{\gamma MC} = 1.05 \times 10^{-4}$.

Despite the apparent order-of-magnitude discrepancy with the estimated accuracy above, it can be shown that this result is in fair agreement with actual expectations. Indeed, it is well known that the $\gamma$ parameter is highly correlated with parallax, which is also estimated in our sphere reconstructions. For Gaia such correlation is similar to that of HIPPARCOS, namely $\rho \approx -0.92$ (Vecchiato et al. 2003). This means that the above prediction for $\sigma_{\gamma}$ should be increased by a factor $(1 - \rho^2)^{-1/2} \approx 2.6$, giving a revised value of $\sigma_{\gamma} \approx 2.6 \times 10^{-5}$, that is, a factor of three better than $\sigma_{\gamma MC}$.

We note, however, that the previous analysis does not take into account the apportionment of the measurement error between $\gamma$ and the other stellar parameters. Moreover, the straightforward error propagation on the light deflection formula is based on a best-case scenario that neglects other sources of errors like, for example, a non-linear effect on the $\gamma$ relative error due to the angular measurement. For these reasons, the factor of three discrepancy above appears quite reasonable.

Finally, following again the findings in (Vecchiato et al. 2003), we can extrapolate the previous value for $\sigma_{\gamma MC}$ to the case of a realistic number of observed stars. At Gaia magnitude $G = 14^2$ the number of primary stars is expected in the order

<footnote n="1">¹ Actually, each least-squares solution provides the adjustment $\delta\gamma_e$ to the $\gamma$ catalog value $\gamma_{cat}$ in the sense $\gamma_e = \gamma_{cat} + \delta\gamma_e$.</footnote>

<footnote n="2">² i.e., the stellar magnitude corresponding to the level of observing error we simulated.</footnote>
of millions; if we then apply their improvement factor for the case of $1 \times 10^6 G = 14$ stars we get $\sigma_{\gamma MC} = 1.05 \times 10^{-4}/80 \approx 1.3 \times 10^{-6}$, a value that compares well with more recent re- estimations of the Gaia potential for the measurement of PPN-$\gamma$ (Mignard & Klioner 2010; de Bruijne 2012).

## 5. Conclusions and perspectives

In this paper we present the latest developments of GSR-TTF, an implementation of the TTF relativistic model within the current version of the Global Sphere Reconstruction software infrastruc- ture for the reduction of Gaia astrometric observations (GSR2) running at the DPAC DPCT. Specifically, we provide the first re- sults of GSR-TTF as applied to five years of DPAC simulated observations.

Light propagation modeling in GSR-TTF is based on the fully relativistic (post-Minkowskian) background of the TTF. This means it can be easily expanded to take into account smaller relativistic effects on light propagation if a higher accuracy is necessary. For instance, the impact of moving axisymmetric bod- ies could be easily implemented by adding their formulation as given, for example, in Hees et al. (2014a).

We first provide analytical equations for both relativistic light propagation and aberration due to Gaia attitude and motion. Then, we show that GSR-TTF results for the direction cosines, that is, the observed direction, of a set of sources are equivalent to those of the GREM implementation in Gaia-Tools (the model and software responsible for the processing of Gaia observations, see, e.g., Balm & GaiaTools Committee 2014; ter Linden & GaiaTools Committee 2009) at the $\mu$s level or better.

The main goal of the Gaia mission is the best ever determination of positions, parallaxes and annual proper motions of all detected stars. Hence, we tested the ability of GSR-TTF to perform such a task on ~132 000 stars by using DPAC simulated observations over 5 yr. The statistical analysis of the results proves GSR-TTF is able to recover astrometric parameters with the expected accuracy when dealing with both realistic measurement and catalog errors. The estimate of instrumental parameters as well as corrections to the Gaia nominal attitude have not been included in this study but are available as part of the GSR framework and software infrastructure and will be experimented with the upcoming versions of the GSR-TTF implementation. Next, we addressed the combined estimate of astrometric parameters with the global PPN-parameter $\gamma$ via a Monte-Carlo experiment. Our results appear quite consistent with the most recent works on measuring PPN-$\gamma$ via Gaia's astrometry.

GSR-TTF is then a powerful processing tool to contribute to the Gaia final catalog, for example, as additional validation resource within the Global Sphere Reconstruction system of AVU. This has already been the case, having allowed to uncover the need of an iteration for non-linearity when solving the sphere with the lower accuracy RAMOD model (RAMOD@GSR2) currently in the operational infrastructure at DPCT.

**Acknowledgements.** The authors acknowledge the financial support of UIF/UFI (French-Italian University) program, CNRS/GRAM, CNES, and of the ASI contract to INAF 2014-025-R.1.2015 (Gaia Mission – The Italian Participation in DPAC). B. Bucciarelli wishes to thank the support of the Chinese Academy of Science President’s International Fellowship initiative (PIFI) for 2017. We wish to thank our DPCT colleagues for their constant support throughout this work. We also thank DPAC-CU2 for making available the Gaia simulated data used in this paper.

## References

Abbas, U., Vecchiato, A., Bucciarelli, B., Lattanzi, M. G., & Morbidelli, R. 2011, EAS Publ. Ser., 45, 127

Balm, P., & GaiaTools Committee. 2014. Livelink Technical Note, Gaia-C1-UG-ESAC-MTL-015-02

Bertone, S., Le Poncin-Lafitte, C., Crosta, M., et al. 2013, in SF2A-2013: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, eds. L. Cambresy, F. Martins, E. Nuss, & A. Palacios, 155

Bertone, S., Minazzoli, O., Crosta, M., et al. 2014, Class. Quant. Grav., 31, 015021

Bini, D., Crosta, M. T., & de Felice, F. 2003, Class. Quant. Grav., 20, 4695

Bombrun, A., Lindegren, L., Holl, B., & Jordan, S. 2010, A&A, 516, A77

Brumberg, V. A. 1991, Essential relativistic celestial mechanics (Bristol, England and New York, Adam Hilger), 271

Crosta, M., & Vecchiato, A. 2010, A&A, 509, A37

Crosta, M., Vecchiato, A., de Felice, F., & Lattanzi, M. G. 2015, Class. Quant. Grav., 32, 165008

Crosta, M., Geralico, A., Vecchiato, A., & Lattanzi, M. 2017, Phys. Rev. D, submitted

Crosta, M. T., & Mignard, F. 2006, Class. Quant. Grav., 23, 4853

de Bruijne, J., Siddiqui, H., Lammers, U., et al. 2010, in IAU Symp. 261, eds. S. A. Klioner, P. K. Seidelmann, & M. H. Soffel, 331

de Bruijne, J. H. J. 2012, Ap&SS, 341, 31

Gaia Collaboration (Prusti, T., et al.) 2016, A&A, 595, A1

Hees, A., Bertone, S., & Le Poncin-Lafitte, C. 2014a, Phys. Rev. D, 90, 084020

Hees, A., Bertone, S., & Le Poncin-Lafitte, C. 2014b, Phys. Rev. D, 89, 064045

Klioner, S. A. 2003, AJ, 125, 1580

Le Poncin-Lafitte, C., Linet, B., & Teyssandier, P. 2004, Class. Quant. Grav., 21, 4463

Le Poncin-Lafitte, C., & Teyssandier, P. 2008, Phys. Rev. D, 77, 044029

Lindegren, L., Lammers, U., Hobbs, D., et al. 2012, A&A, 538, A78

Linet, B., & Teyssandier, P. 2002, Phys. Rev. D, 66, 024045

Linet, B., & Teyssandier, P. 2013, Class. Quant. Grav., 30, 175008

Luri, X., Martinez, O., & Isasi, Y. 2010, Livelink Technical Note, Gaia-C2-TR-UB-XL-020-01

Messineo, R., Morbidelli, R., Martino, M., et al. 2012, in Software and Cyberinfrastructure for Astronomy II, Proc. SPIE, 8451, 84510E

Mignard, F., & Klioner, S. A. 2010, in Relativity in Fundamental Astronomy: Dynamics, Reference Frames, and Data Analysis, eds. S. A. Klioner, P. K. Seidelmann, & M. H. Soffel, IAU Symp., 261, 306

Mignard, F., Bailer-Jones, C., Bastian, U., et al. 2008, in IAU Symp. 248, eds. W. J. Jin, I. Platais, & M. A. C. Perryman, 224

Misner, C., Thorne, K., & Wheeler, J. 1973, Gravitation (San Francisco: W. H. Freeman)

Perryman, M. A. C., de Boer, K. S., Gilmore, G., et al. 2001, A&A, 369, 339

ter Linden, M., & GaiaTools Committee. 2009. Livelink Technical Note, Gaia-C1-SP-ESAC-MTL-019-01

Teyssandier, P., & Le Poncin-Lafitte, C. 2008, Class. Quant. Grav., 25, 145020

Teyssandier, P., Le Poncin-Lafitte, C., & Linet, B. 2008, in Lasers, Clocks and Drag-Free Control: Exploration of Relativistic Gravity in Space, eds. H. Dittus, C. Lammerzahl, & S. G. Turyshev, Astrophys. Space Sci. Lib., 349, 153

Vecchiato, A., Lattanzi, M. G., Bucciarelli, B., et al. 2003, A&A, 399, 337

Vecchiato, A., Abbas, U., Bandieramonte, M., et al. 2012, in SPIE Conf. Ser., 8451

Vecchiato, A., Gai, M., Lattanzi, M. G., et al. 2014, J. Phys. Conf. Ser., 490, 012241

Will, C. M. 2014, Liv. Rev. Relat., 17, 4