© ESO, 2014

1. Introduction

Stellar proper motions have traditionally been determined by analysing the differences in position at different epochs, often separated by many decades and obtained using vastly different instruments and methods. In this process, parallaxes (and radial motions, albeit relevant to a much lesser extent) were mostly ignored.

With the advent of space astrometry, most notably the European satellite Hipparcos (1989–1993, see ESA 1997), it became necessary to treat data in a unified manner, i.e., by applying a single least-squares solution for the position, parallax, and annual proper motion. Hipparcos determined these parameters for nearly 120 000 stars¹ mostly brighter than magnitude 12, with a median uncertainty of about 1 milli-arcsecond (mas). The Tycho-2 Catalogue (Høg et al. 2000) gave additional data for 2.5 million stars observed with the Hipparcos starmappers. The re-reduction of the Hipparcos raw data (van Leeuwen 2007a,b) significantly improved the main-mission results. Today, 25 years after the launch of the satellite, these catalogues remain the main source for the astrometric parameters of these stars.

The European space astrometry mission Gaia will soon change this picture. Gaia, launched at the end of 2013, will determine the astrometric parameters of up to a billion stars between magnitude 6 and 20 with unprecedented accuracies reaching a few tens of micro-arcseconds (μas) for Gaia magnitude G ≲ 15. The vast amounts of data will be processed in a single coherent least-squares solution, which solves not only for the astrometric parameters but also for a large number of parameters describing the time-varying spacecraft attitude and the geometry of the optical instrument. Due to the very large number of parameters to be determined from the observational data the system cannot be solved directly(Bombrun et al. 2010)but has to be tackled in a block-iterative manner with the so called “Astrometric Global Iterative Solution” (AGIS). The AGIS software has been designed and implemented by groups at ESA/ESAC, Lund Observatory, and others, and is described in detail together with the fundamental algorithms and mathematical framework byLindegren et al. (2012) .

Astrometric measurements obtained in the past, even of moderate accuracy by modern standards, have lasting value as they represent a state of the Universe that is never repeated. A good example is the construction of proper motions in the Tycho-2 Catalogue using Hipparcos and century-old photographic positions. When the astrometric parameters are propagated over a long time interval, uncertainties in the tangential and radial motions accumulate to a significant positional uncertainty. Yet, long-term deviations from linear space motion (e.g., in long-period binaries) increase even more drastically with time. Such deviations might not be detectable within the time spans of the Hipparcos or Gaia missions individually, but could be detectable by combining the results of the two. Thus, although Hipparcos will soon be superseded by Gaia in terms of the expected accuracies at current epochs, its data form a unique comparison point in the past, very valuable in combination with later results. For this reason the first Gaia data release scheduled for 2016 will not only publish stellar positions and magnitudes based on the first Gaia observations, but also a combination of these observations with the Hipparcos Catalogue for all stars common between the two missions. This part of the release is called the Hundred Thousand Proper Motions project (HTPM), originally proposed by F. Mignard in a Gaia-internal technical document (Mignard 2009).

The present paper gives a recipe for the practical realisation of the HTPM project in the context of the already existing AGIS scheme for the astrometric solution of Gaia data. The proper motions in HTPM might be trivially computed from the positional differences between an early Gaia solution and the Hipparcos Catalogue – the “conventional catalogue combination” approach of Sect. 2.3. However, we argue that the more elaborate “joint solution” method described in Sect. 2.4 will have important advantages for the HTPM project, and in Sect. 3 we show how to implement it as part of AGIS. The validity and accuracy of the method is demonstrated by means of a joint solution of simulated Gaia observations of the Hipparcos stars (Sect. 4). In the final sections we discuss the limitations of the results and their validity in the light of Gaia’s full nominal mission performance, as well as possible applications of the joint solution method to other astrometric data.

The HTPM project should use the re-reduction of the raw Hipparcos data (van Leeuwen 2007b), as it represents a significant improvement over the original Hipparcos Catalogue (ESA 1997). Therefore it is also used in all our simulations. For the purpose of demonstrating the HTPM solution we regard all valid entries of the Hipparcos Catalogue as astrometrically well-behaved (effectively single) stars. Their space motions are therefore regarded as uniform (rectilinear, with constant speed) over the time interval covered by Hipparcos and Gaia. This is obviously a very simplified picture of the true content of the Hipparcos Catalogue. However, getting the solution right in this simple case is a first necessary step for any more sophisticated treatment of detected binaries and multiple stars in the Hipparcos Catalogue.

2. Theory

Combining astrometric catalogues requires that data are expressed in the same reference system and described in terms of a common kinematic model. In this section we describe the adopted model and how it is connected to the definition of the astrometric parameters. We outline the conventional approach to catalogue combination and develop the “joint solution” as an optimal generalisation of the method. We show how to detect deviations from the kinematic model or misfits between the datasets. We also outline how to reconstruct the required information from Hipparcos and how to integrate the proposed scheme in the astrometric solution algorithm of Gaia.

2.1. Kinematic model of stellar motion

The choice of astrometric parameters is a direct result of choosing a model of stellar motion. The most basic assumption is for stars to move uniformly, i.e., linearly and with constant speed, relative to the solar system barycentre (SSB). Note that this also means that the stars are assumed to be single. This is obviously not true for all of them, but a good basic assumption for most stars. During the data reduction stars that are not “well behaved” in an astrometric sense can be filtered out and treated further, e.g., by adding additional parameters for components of stellar systems or for acceleration through external influences.

A uniform space motion can be fully described by six parameters: three for the position in space at a chosen reference epoch, and three for the velocity. Traditionally, the three positional parameters are right ascension α, declination δ, and parallax ϖ relative to the SSB at the reference epoch of the catalogue. The motion is then described by three parameters, where and are the proper motions in right ascension and declination, respectively, and the third parameter μ_r is the radial motion component. The radial component is more commonly given as the radial velocity v_r in km s^-1, but in an astrometric context it is conveniently expressed as the radial proper motion (equivalent to the relative change in distance over time, or ) ${\mathit{\mu }}_{\mathit{r}}\mathrm{=}{\mathit{v}}_{\mathit{r}}\mathit{\varpi }\mathit{/}\mathit{A,}$(1)where A is the astronomical unit expressed in km yr s^-1. Only the first five parameters are classically considered astrometric parameters. Based on only a few years of observations it is usually not possible to determine the radial component from astrometry with sufficient accuracy(Dravins et al. 1999) . Hence the radial component is better determined by other techniques, i.e., from spectroscopy. For Gaia the radial component will be significant for many more stars, although the affected fraction remains very small(de Bruijne & Eilers 2012) . Even though μ_r is not determined in the astrometric solution for the vast majority of sources, it is convenient and sometimes necessary to formulate astrometric problems with the full set of six astrometric parameters, as we do in this paper. We will also show how to treat the sixth component when the radial velocity is unknown or added from spectroscopy.

2.2. Dealing with non-linearities: SMOK

When comparing and subsequently combining astrometric catalogues one needs to deal with the fact that the mapping from rectilinear to spherical coordinates is strongly non-linear. This becomes significant at the μas level when the differences in α and δ exceed some (1 μas)^1/2 ≃ 0.5 arcsec. For example, the barycentric direction traced out in α(t), δ(t) due to the proper motion will not be linear even though the star is assumed to move uniformly through space. The traditional way to deal with this is to introduce higher-order correction terms computed by Taylor expansion of the rigorous equations (e.g., Taff 1981). In this paper we take a different approach, based on the scaled modelling of kinematics (SMOK) concept described in Appendix A. For the present purpose it is sufficient to know that (α,δ) may be replaced by linear coordinates (a,d) relative to a designated, fixed comparison point, with time derivatives ȧ, ḋ representing the components of proper motion in α and δ. The six parameters a, d, ϖ, ȧ, ḋ, ṙ (where ṙ is the SMOK equivalent of the radial proper motion) provide an alternative and equivalent parametrisation of the kinematics, more convenient for the catalogue combination than the usual set α, δ, ϖ, μ_{α ∗}, μ_δ, μ_r.

2.3. Conventional catalogue combination

In the conventional catalogue combination the astrometric parameters in each catalogue are independently estimated from separate sets of observations, and the combination is done a posteriori from the individual catalogues. Let (a₁,d₁,ϖ₁) at time t₁ be the position and parallax of a star in the first catalogue, and (a₂,d₂,ϖ₂) at time t₂ the corresponding information in the second catalogue. The proper motion parameters ȧ, ḋ are then derived as the positional difference over time Δt = t₂ − t₁$\begin{array}{ccc}\mathit{ȧ}\mathrm{=}\mathrm{(}{\mathit{a}}_{\mathrm{2}}\mathrm{-}{\mathit{a}}_{\mathrm{1}}\mathrm{)}\mathit{/}\mathrm{\Delta }\mathit{t,} \mathit{ḋ}\mathrm{=}\mathrm{(}{\mathit{d}}_{\mathrm{2}}\mathrm{-}{\mathit{d}}_{\mathrm{1}}\mathrm{)}\mathit{/}\mathrm{\Delta }\mathit{t,}& & \end{array}$(2)which is possible thanks to the reformulation of the astrometric parameters in SMOK. The proper motion uncertainties are $\begin{array}{ccc}{\mathit{\sigma }}_{\mathit{ȧ}}\mathrm{=}\frac{\sqrt{{\mathit{\sigma }}_{\mathit{a}\mathrm{1}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathit{a}\mathrm{2}}^{\mathrm{2}}}}{\mathrm{\Delta }\mathit{t}}\mathit{,} {\mathit{\sigma }}_{\mathit{ḋ}}\mathrm{=}\frac{\sqrt{{\mathit{\sigma }}_{\mathit{d}\mathrm{1}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathit{d}\mathrm{2}}^{\mathrm{2}}}}{\mathrm{\Delta }\mathit{t}}\mathit{,}& & \end{array}$(3)where σ_a1 is the uncertainty of a₁, etc. The third kinematic parameter ṙ for the radial motion could in theory be derived from the (negative, relative) difference in parallax, but in practice it is derived from the spectroscopic radial velocity as discussed in Sect. 2.1.

While the proper motions are obtained by taking position differences over time, the combined parameters for position and parallax are formed as weighted means. For a this gives $\begin{array}{ccc}\mathit{â}\mathrm{=}\frac{{\mathit{a}}_{\mathrm{1}}{\mathit{\sigma }}_{\mathit{a}\mathrm{1}}^{-2}\mathrm{+}{\mathit{a}}_{\mathrm{2}}{\mathit{\sigma }}_{\mathit{a}\mathrm{2}}^{-2}}{{\mathit{\sigma }}_{\mathit{a}\mathrm{1}}^{-2}\mathrm{+}{\mathit{\sigma }}_{\mathit{a}\mathrm{2}}^{-2}}\mathit{,}& & \end{array}$(4)referring to the mean epoch of the combination $\begin{array}{ccc}\mathit{t̂}\mathit{a}\mathrm{=}\frac{{\mathit{t}}_{\mathrm{1}}{\mathit{\sigma }}_{\mathit{a}\mathrm{1}}^{-2}\mathrm{+}{\mathit{t}}_{\mathrm{2}}{\mathit{\sigma }}_{\mathit{a}\mathrm{2}}^{-2}}{{\mathit{\sigma }}_{\mathit{a}\mathrm{1}}^{-2}\mathrm{+}{\mathit{\sigma }}_{\mathit{a}\mathrm{2}}^{-2}}\mathrm{·}& & \end{array}$(5)The reference time is the optimal time in-between the two catalogues at which the position and proper motion are uncorrelated and the uncertainty of â is minimal, given by ${\mathit{\sigma }}_{\mathit{â}}^{-2}\mathrm{=}{\mathit{\sigma }}_{\mathit{a}\mathrm{1}}^{-2}\mathrm{+}{\mathit{\sigma }}_{\mathit{a}\mathrm{2}}^{-2}$. The expressions for and are analogous.

This combination scheme has some limitations, in that it does not take correlations between the astrometric parameters into account, nor the individual proper motions that may exist in each catalogue. In the next section we describe a more general approach.

2.4. Joint solution

The reduction of astrometric data is typically done using least-squares solutions, resulting in a linear system of normal equations Nx = b. Here, x is the vector of resulting astrometric parameters, N the normal equations matrix, and b a vector constructed from the residuals of the problem². The covariance C of the solution is formally given by C = N^-1.

In AGIS the observations of all well-behaved stars (“primary sources”) must be considered together in a single, very large least-squares solution (Sect. 2.7). For n primary sources, x would then be the full vector of 6n astrometric parameters, with N and b of corresponding dimensions. However, for the present exposition it is sufficient to consider one star at a time, so that x and b are of length 6 and N has dimensions 6 × 6. In practice only five of the six parameters are estimated, and N^-1 should hereafter be regarded as the inverse of the upper-left 5 × 5 submatrix³.

On the assumption that the adopted kinematic model is valid for a particular star, the matrix N and vector b encapsulate the essential information on the astrometric parameters, as determined by the least-squares solution. Thus, in order to make optimal use of the Hipparcos data for a given star there is no need to consider the individual observations of that star: all we need is contained in the “information array” [ Nb ]. In Sect. 2.6 we show how this array is reconstructed from the published Hipparcos Catalogue.

Let [ N₁b₁ ] and [ N₂b₂ ] be the information arrays for the same star as given by two independent astrometric catalogues. From the way the normal equations are calculated from observational data it is clear that the information arrays are additive, so that [ N₁b₁ ] + [ N₂b₂ ] is the information array that would have resulted from processing the two datasets together. InMichalik et al. (2012)we have proposed that the optimum combination of the catalogues is done a priori, that is by adding the corresponding arrays before solving. The result, $\begin{array}{ccc}{x̂}\mathrm{=}\mathrm{(}{{N}}_{\mathrm{1}}\mathrm{+}{{N}}_{\mathrm{2}}{\mathrm{)}}^{-1}\mathrm{(}{{b}}_{\mathrm{1}}\mathrm{+}{{b}}_{\mathrm{2}}\mathrm{)}\mathit{,}& & \end{array}$(6)is the joint solution of the astrometric parameters, with covariance Ĉ = (N₁ + N₂)^-1. The two catalogue entries for the star must use the same reference epoch and the same SMOK comparison point.

The joint solution has several advantages over the conventional combination method outlined in Sect. 2.3. Because it uses the full information in each catalogue it makes better use of the data and allows one to estimate the resulting uncertainties more accurately, taking the correlations into account. The individual proper motion information available in each catalogue is automatically incorporated in the joint proper motion. Moreover, a solution might be possible where the data in each set individually is insufficient to solve for all astrometric parameters, that is, N₁ + N₂ may be non-singular even if N₁, N₂, or both, are singular. In practice, if N₁ comes from the Hipparcos data, it will always be non-singular (since there is a Hipparcos solution), and the sum is then also non-singular. Hence it will always be possible to make a joint solution for all five astrometric parameters of the HTPM stars. Finally, the joint solution scheme is a clean and rigorous approach and can be integrated into the existing implementation of the astrometric solution for Gaia with moderate effort.

The joint solution can be seen as a multidimensional generalisation of the conventional scheme in Sect. 2.3, with N representing the weights (σ^-2) and b the astrometric parameters multiplied by their weights (e.g., aσ^-2). Then Eq. (6) is the matrix-equivalent of Eq. (4). The joint solution can also be understood in terms of Bayesian estimation theory (assuming multivariate Gaussian parameter errors), with N₁,b₁ representing the prior information, N₂,b₂ the new data, and their sums the posterior information.

2.5. Goodness of fit of the joint solution

The goodness of fit of a least-squares solution can be described in terms of the sum of the squares of the normalized post-fit residuals, $\mathit{Q}\mathrm{=}\sum _{\mathit{k}}{\left(\frac{{\mathit{\eta }}_{\mathit{k}}^{\mathrm{\left(}\mathrm{obs}\mathrm{\right)}}\mathrm{-}{\mathit{\eta }}_{\mathit{k}}^{\mathrm{\left(}\mathrm{calc}\mathrm{\right)}}}{{\mathit{\sigma }}_{\mathit{k}}}\right)}^{\mathrm{2}}\mathit{,}$(7)where ${\mathit{\eta }}_{\mathit{k}}^{\mathrm{\left(}\mathrm{obs}\mathrm{\right)}}$ and ${\mathit{\eta }}_{\mathit{k}}^{\mathrm{\left(}\mathrm{calc}\mathrm{\right)}}$ are the observed and calculated (fitted) angular focal-plane coordinates of the star in observation k, and σ_k is the standard error of the observation. Q is calculated for each star separately and is simply a function of x = (a,d,ϖ,ȧ,ḋ,ṙ)′. The least-squares solution minimizes Q and for any other parameter vector x we have $\mathit{Q}\mathrm{\left(}{x}\mathrm{\right)}\mathrm{=}\mathit{Q}\mathrm{\left(}{x̂}\mathrm{\right)}\mathrm{+}\mathrm{(}{x}\mathrm{-}{x̂}{\mathrm{)}}^{\mathrm{\prime }}{N}\mathrm{(}{x}\mathrm{-}{x̂}\mathrm{)}\mathit{.}$(8)If the kinematic model is correct and the standard errors of the observations are correctly estimated one expects the minimum value to follow the chi-square distribution with ν degrees of freedom, . Here ν = m − rank(N) is equal to the number of observations m (that is the number of terms in Eq. (7)) diminished by the rank of N. Note that this holds even if N is singular (i.e., rank(N) <n, where n is the number of fitted parameters). In the singular case is not unique, yet has a well-defined value (which may be 0 or positive).

Analogous to Eq. (8), in the joint solution we minimize the total goodness of fit, $\begin{array}{ccc} \mathit{Q}\mathrm{\left(}{x}\mathrm{\right)}& \mathrm{=}& {\mathit{Q}}_{\mathrm{1}}\mathrm{\left(}{x}\mathrm{\right)}\mathrm{+}{\mathit{Q}}_{\mathrm{2}}\mathrm{\left(}{x}\mathrm{\right)}\mathrm{=}{\mathit{Q}}_{\mathrm{1}}\mathrm{\left(}{{x̂}}_{\mathrm{1}}\mathrm{\right)}\mathrm{+}{\mathit{Q}}_{\mathrm{2}}\mathrm{\left(}{{x̂}}_{\mathrm{2}}\mathrm{\right)}\\ & & \mathrm{+}\mathrm{(}{x}\mathrm{-}{{x̂}}_{\mathrm{1}}{\mathrm{)}}^{\mathrm{\prime }}{{N}}_{\mathrm{1}}\mathrm{(}{x}\mathrm{-}{{x̂}}_{\mathrm{1}}\mathrm{)}\mathrm{+}\mathrm{(}{x}\mathrm{-}{{x̂}}_{\mathrm{2}}{\mathrm{)}}^{\mathrm{\prime }}{{N}}_{\mathrm{2}}\mathrm{(}{x}\mathrm{-}{{x̂}}_{\mathrm{2}}\mathrm{)}\mathit{.}\end{array}$(9)Here ${{x̂}}_{\mathit{i}}\mathrm{=}{{N}}_{\mathit{i}}^{-1}{{b}}_{\mathit{i}}$ is the solution obtained by using only catalogue i = 1, 2, i.e., minimizing Q_i(x), which results in the minimum value . It is readily seen that Eq. (9) is minimized precisely for the joint solution vector in Eq. (6).

Each of the four terms in Eq. (9) has a simple interpretation. The first term, , is the chi-square obtained when fitting the astrometric parameters only to the first set of data (in our case the Hipparcos data); similarly, is the chi-square obtained when fitting only to the second set of data (from Gaia). The sum of the last two terms is minimized for , and shows how much the chi-square is increased by forcing the same parameters to fit both sets of data in the joint solution. This quantity is useful for assessing whether the two datasets are mutually consistent and we therefore introduce a separate notation for it, $\mathrm{\Delta }\mathit{Q}\mathrm{=}\mathrm{(}{x̂}\mathrm{-}{{x̂}}_{\mathrm{1}}{\mathrm{)}}^{\mathrm{\prime }}{{N}}_{\mathrm{1}}\mathrm{(}{x̂}\mathrm{-}{{x̂}}_{\mathrm{1}}\mathrm{)}\mathrm{+}\mathrm{(}{x̂}\mathrm{-}{{x̂}}_{\mathrm{2}}{\mathrm{)}}^{\mathrm{\prime }}{{N}}_{\mathrm{2}}\mathrm{(}{x̂}\mathrm{-}{{x̂}}_{\mathrm{2}}\mathrm{)}\mathit{.}$(10)The two terms give the increase in chi-square due to the first and second dataset, respectively.

Long-period astrometric binaries may have significantly different proper motions at the Hipparcos and Gaia epochs, and these in turn may differ from the mean proper motion between the epochs. If the differences are significant, compared with the measurement precisions, they will result in an increased value of ΔQ. The null hypothesis, namely that the star is astrometrically well-behaved, should be rejected if ΔQ exceeds a certain critical value. In order to calculate the critical value it is necessary to know the expected distribution of ΔQ under the null hypothesis.

Let m_i and ν_i = m_i − rank(N_i) be the number of observations and degrees of freedom in catalogue i. The number of degrees of freedom in the joint solution is ν = (m₁ + m₂) − rank(N₁ + N₂). Under the null-hypothesis we have (i = 1, 2), , and consequently $\mathrm{\Delta }\mathit{Q}\mathrm{~}{\mathit{\chi }}^{\mathrm{2}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{,}$(11)where $\mathit{k}\mathrm{=}\mathit{\nu }\mathrm{-}{\mathit{\nu }}_{\mathrm{1}}\mathrm{-}{\mathit{\nu }}_{\mathrm{2}}\mathrm{=}\mathrm{rank}\mathrm{\left(}{{N}}_{\mathrm{1}}\mathrm{\right)}\mathrm{+}\mathrm{rank}\mathrm{\left(}{{N}}_{\mathrm{2}}\mathrm{\right)}\mathrm{-}\mathrm{rank}\mathrm{(}{{N}}_{\mathrm{1}}\mathrm{+}{{N}}_{\mathrm{2}}\mathrm{)}\mathit{.}$(12)In the special case when N₁, N₂, and N₁ + N₂ all have full rank (equal to n, the number of astrometric parameters) we have k = n. At a significance level of 1% the critical values of ΔQ, above which the null hypothesis should be rejected, are 15.086, 13.277, 11.345, 9.210, and 6.635 for k = 5, 4, 3, 2, and 1, respectively (e.g., Abramowitz & Stegun 2012). With this criterion only 1% of the well-behaved stars should be accidentally misclassified as not well-behaved. The expected distribution of ΔQ can be verified in the simulations which, by design, only includes well-behaved stars.

2.6. Reconstruction of N_H, b_H for the HIPPARCOS Catalogue

When using the joint solution for incorporating Hipparcos data in the solution of early Gaia data it is necessary to reconstruct the normal matrix N_H and the right hand side b_H from Hipparcos for each star. These are initially calculated for the reference epoch of the Hipparcos Catalogue (J1991.25) and later propagated to the adopted reference epoch of the joint solution (see Sect. 2.7).

Fig. 1 Relationships between catalogues during simulation runs.

Let a_H, d_H, ϖ_H, ȧ_H, ḋ_H be the astrometric parameters from the Hipparcos Catalogue after transformation into the SMOK notation (see Appendix A). The upper-left 5 × 5 submatrix of the covariance matrix can be taken without changes from the Hipparcos Catalogue (see Appendix B for details) since σ_{α ∗} = σ_a, σ_δ = σ_d, ... with sufficient accuracy at the reference epoch of the catalogue and provided that the SMOK comparison point is close enough to the astrometric parameters of the star. The sixth parameter ṙ_H and its corresponding entries in the covariance matrix need to be added from external sources or set to sensible values if not available (see below). Then the normal matrix is simply the inverse of the covariance matrix ${{N}}_{\mathrm{H}}\mathrm{=}{{C}}_{\mathrm{H}}^{-1}$ and $\begin{array}{ccc}{{b}}_{\mathrm{H}}\mathrm{=}{{N}}_{\mathrm{H}}{\mathrm{(}\begin{array}{c}\\ \end{array}\mathrm{)}}^{\mathrm{\prime }}\mathit{.}& & \end{array}$(13)ESA (1997), Vol. 1, Eq. (1.5.69) shows how to reconstruct the elements [ C₀ ] _i6 = [ C₀ ] _6i(i = 1...6), that is the sixth column and row of the covariance matrix corresponding to the radial motion μ_r. Let , , be the true values and δv_r, δϖ, δμ_r the errors. The expression in Eq. (1.5.69) for the diagonal element [ C₀ ] ₆₆ is only valid if the relative uncertainties in the radial velocity and parallax are small, i.e., , . If this is not the case we need to consider the complete expression for the calculated radial motion, $\begin{array}{ccc}{\mathit{\mu }}_{\mathit{r}}\mathrm{=}\mathit{\mu ̅}\mathit{r}\mathrm{+}\mathit{\delta }{\mathit{\mu }}_{\mathit{r}}\mathrm{=}\mathrm{(}\mathit{v̅}\mathit{r}\mathrm{+}\mathit{\delta }{\mathit{v}}_{\mathit{r}}\mathrm{)}\mathrm{(}\mathit{\varpi ̅}\mathrm{+}\mathit{\delta \varpi }\mathrm{)}\mathit{/}\mathit{A,}& & \end{array}$(14)where , leading to $\begin{array}{ccc}\mathit{\delta }{\mathit{\mu }}_{\mathit{r}}\mathrm{=}\mathrm{(}\mathit{v̅}\mathit{r}\mathit{\delta \varpi }\mathrm{+}\mathit{\varpi ̅}\mathit{\delta }{\mathit{v}}_{\mathit{r}}\mathrm{+}\mathit{\delta }{\mathit{v}}_{\mathit{r}}\mathit{\delta \varpi }\mathrm{)}\mathit{/}\mathit{A}\mathit{.}& & \end{array}$(15)Squaring and taking the expectation while assuming that the errors in parallax and radial velocity are uncorrelated gives $\begin{array}{ccc}\mathit{E}\mathrm{\left(}\mathit{\delta }{\mathit{\mu }}_{\mathit{r}}^{\mathrm{2}}\mathrm{\right)}\mathrm{=}\mathrm{(}\mathit{E}\mathrm{\left(}{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}\mathit{\delta }{\mathit{\varpi }}^{\mathrm{2}}\mathrm{\right)}\mathrm{+}\mathit{E}\mathrm{\left(}{\mathit{\varpi }}^{\mathrm{2}}\mathit{\delta }{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}\mathrm{\right)}\mathrm{+}\mathit{E}\mathrm{(}\mathit{\delta }{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}\mathit{\delta }{\mathit{\varpi }}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{)}}\mathit{/}{\mathit{A}}^{\mathrm{2}}\mathit{,}& & \end{array}$(16)where we replaced the true quantities by the observed ones. The third term is the required generalisation if v_r or ϖ is zero, or if the relative errors are large. For example, if parallax and radial motion are unknown they could be assumed to be zero with a large uncertainty. The generalized version of Eq. (1.5.69) in ESA (1997) reads $\begin{array}{ccc}{\left[{{C}}_{\mathrm{0}}\right]}_{\mathrm{66}}& \mathrm{=}& {\mathrm{(}{\mathit{v}}_{\mathit{r}\mathrm{0}}\mathit{/}{\mathit{A}}^{\mathrm{)}}}^{\mathrm{2}}{\left[{{C}}_{\mathrm{0}}\right]}_{\mathrm{33}}\mathrm{+}{\mathrm{(}{\mathit{\varpi }}_{\mathrm{0}}\mathit{/}{\mathit{A}}^{\mathrm{)}}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathit{v}\mathit{r}\mathrm{0}}^{\mathrm{2}}\mathrm{+}{\mathrm{(}{\mathit{\sigma }}_{\mathit{v}\mathit{r}\mathrm{0}}\mathit{/}{\mathit{A}}^{\mathrm{)}}}^{\mathrm{2}}{\left[{{C}}_{\mathrm{0}}\right]}_{\mathrm{33}}\mathit{,}\\ {\left[{{C}}_{\mathrm{0}}\right]}_{\mathit{i}\mathrm{6}}& \mathrm{=}& \end{array}$(17)The Hipparcos Catalogue contains numerous entries for non-single stars, for which additional parameters are given, describing deviations from uniform space motion. These additional parameters are ignored in our simulations, which regard every star as single. In the actual HTPM solution many of these stars may require more specialised off-line treatment. This is not further discussed in this paper.

2.7. Joint solution in AGIS

In reality the astrometric solution cannot be done separately for each star as described in Sect. 2.4 but must consider all the stars together with the spacecraft attitude and instrument calibration. Without prior information on the astrometric parameters this leaves the solution undetermined with respect to the reference frame. This is not the case for the joint solution, however, as the Hipparcos prior information contains positions and proper motions that are expressed in a specific reference frame, namely the Hipparcos realisation of the International Celestial Reference System (ICRS;Feissel & Mignard 1998 ). The incorporation of the Hipparcos prior in the joint solution automatically ensures that the resulting data are on the Hipparcos reference frame. If required, the data can later be transformed into a more accurate representation of the ICRS (see Sect. 5.3).

Due to the size of the data reduction problem, AGIS does not directly solve Nx = b but iteratively improves the astrometric parameters by computing the updates Δx, i.e., the difference to the current best estimate values. When incorporating Hipparcos data this requires us to also express the Hipparcos data (subscript H) as a difference to the current best estimate (subscript c). Therefore we construct $\begin{array}{ccc}\mathrm{\Delta }{{b}}_{\mathrm{H}}\mathrm{=}{{N}}_{\mathrm{H}}\mathrm{\Delta }{x}\mathrm{=}{{N}}_{\mathrm{H}}\left(\begin{array}{c}\\ \\ \\ \\ \\ \\ \end{array}\right)\mathit{.}& & \end{array}$(18)Before solving we add the corresponding matrices for the Gaia data. If no additional Gaia data would be added the solution would immediately recover the Hipparcos Catalogue parameters.

The reference epoch of the joint solution can be arbitrarily chosen. In practice the Gaia data are much better than the Hipparcos data, therefore the optimal reference epoch would always be very close to the epoch of the Gaia data alone. Assuming one releases Gaia-only data and HTPM results at the same time it might be convenient to publish both for the same reference epoch, i.e., the Gaia-only reference epoch of the data release.

3. Simulations

3.1. Logic of simulations

Simulations are based on AGISLab(Holl et al. 2012) , a small-scale version of the AGIS data reduction created and maintained at Lund Observatory. It is used to aid the development of algorithms for the astrometric data reduction of Gaia. Simulation runs are carried out in the following steps (cf. Fig. 1):

• 1.

  Creating catalogues of all the stars used in the simulation, namely the Hipparcos stars and the auxiliary stars (see below). Two catalogues are needed: a simulated “true” catalogue to generate Gaia observations and to evaluate the uncertainties of the astrometric performance, and an initial catalogue of starting values for the data reduction.

• 2.

  Simulating observations of the stars using the Nominal Scanning Law(de Bruijne et al. 2010) , including perturbations according to the expected precision of Gaia measurements.

• 3.

  Improving the astrometry of the initial catalogue through the astrometric solution (AGIS), resulting in the final catalogue. This can be done with or without incorporation of prior information from Hipparcos.

• 4.

  Evaluating the error of the resulting solution by comparing the final catalogue with the true catalogue.

Details of the first two steps are given below, while remaining steps are covered in Sect. 4.

3.2. Simulating the stellar catalogues

All catalogues consist of two parts, the Hipparcos stars and the additional auxiliary stars. The Hipparcos stars are necessary for the realisation of the HTPM scheme, and 113 396 stars are within the nominal magnitude range of Gaia (G ≃ 6–20). In order to obtain a reliable astrometric solution with a realistic modelling of the attitude constraints we find that a minimum of one million stars is needed, uniformly distributed on the sky. 886 604 auxiliary stars are therefore added to the Hipparcos stars in the solution. The astrometric results for the auxiliary stars are not included in the statistics for the HTPM performance, which is based only on the results for the Hipparcos stars. However, they contribute indirectly to the HTPM solution via the attitude.

3.2.1. Simulated “true” catalogue

The true catalogue defines the stars used for creating the simulated Gaia observations. For the real mission the true catalogue is of course not known.

To derive the Hipparcos portion of the true catalogue we assume that the true parameters deviate from the Hipparcos values by random amounts consistent with the Hipparcos covariances. The Hipparcos Catalogue is taken from CDS and contains the astrometric parameters for the reference epoch J1991.25, including their covariance matrices (Appendix B). For each star let C be its covariance matrix, L the lower triangular matrix resulting from the Cholesky decomposition C = LL′, and g a vector of six independent standard Gaussian random variables (zero mean, unit standard deviation). Then the true parameters (subscript T) are obtained by applying the error vector e = Lg to the astrometric parameters from the Hipparcos Catalogue (subscript H): $\begin{array}{ccc}{{x}}_{\mathrm{T}}& \mathrm{=}& \end{array}$(19)Since E(g) = 0, where E(...) denotes the expectation value, it follows that E(e) = 0. Moreover, since E(gg′) = I (the identity matrix), it is readily verified that e has the desired covariance E(ee′) = C. For a joint solution with simulated Gaia data the Hipparcos Catalogue needs to be propagated to the reference epoch used in the solution.

Rigorous propagation of the astrometric parameters must take into account the radial motions of the stars, for which radial velocities are needed. We use data from XHIP(Anderson & Francis 2012) , a compilation of radial velocities and other data for the Hipparcos stars from 47 different sources. We only use radial velocities with quality flag “A” or “B” in XHIP. This makes for a total of 40 171 radial velocities which are used as true values in our simulations. For the remaining Hipparcos stars we assign random radial velocities from a Gaussian distribution with v_r = 0, σ_vr = 30 km s^-1 using Eq. (19), based on the assumption that radial velocities are typically smaller than that. The radial velocity uncertainty (taken from XHIP or using 30 km s^-1) is also used to expand the 5 × 5 covariance matrix by a sixth column and row for the uncertainty and correlation of the radial motion, using Eq. (17).

For the auxiliary stars, the positions are chosen to give a random uniform distribution across the sky with a mean density of about 21 stars deg^-2, corresponding to one million stars needed for the solution. We assume magnitude G = 13 for all auxiliary stars. Since the number density of actual stars with G ≤ 13 is about 60 deg^-2 at the Galactic poles, the assumed distribution is a rather conservative estimate of the density of bright stars available for the astrometric solution. The parallaxes of the auxiliary stars are assumed to have a log-normal distribution with median parallax 2.5 mas and a standard deviation of 0.6 dex⁴. The true proper motions and radial velocities are calculated by assuming an isotropic velocity distribution relative to the Sun with a standard deviation of 30 km s^-1.

3.2.2. Initial catalogue and astrometric solution

The initial catalogue contains the starting values for the data processing. The Hipparcos portion of it is identical to the astrometric parameters read from the Hipparcos Catalogue. For the auxiliary stars the initial positions are obtained by perturbing the true positions with Gaussian noise of standard deviation 100 mas in each coordinate, while the initial parallax and proper motion are set to zero. This is similar to a real life scenario where one would assume initial stellar positions from ground based observations or the first published Gaia positions without additional knowledge on the parallax or proper motion. The astrometric values in the initial catalogue are subsequently updated by the AGIS processing, resulting in the final catalogue once the solution is found. We do not solve for the radial motion but set the radial velocity to either zero (assuming no knowledge about it) or the true value (assuming it is perfectly known). In the first case perspective acceleration may show up for some stars as discrepancies in the solution, which disappear when the true radial velocities are used instead (see Sect. 4.3).

3.2.3. Final catalogue

The final catalogue contains the astrometric parameters after data processing. The difference to the simulated true catalogue gives the final errors of the reduced data and is used to evaluate the quality of the astrometric results. In this evaluation we focus on the improvement in the astrometric parameters of the Hipparcos stars.

3.3. Simulating Gaia observations

The observations of the one million stars described above are simulated using the Nominal Scanning Law of Gaia. We neglect so called “dead time” (when no data can be accumulated for example due to orbit maintenance manoeuvres and micro-meteoroid hits), which may amount to up to 15% of the mission time. We do, however, account for the dead time originating from stellar transits coinciding with gaps between the CCDs in the focal plane, i.e., our simulations remove such observations before further processing of the data.

To account for observation noise, i.e., the expected centroiding performance of Gaia, we use a simplified noise model that ignores the gating scheme that Gaia exploits for bright star detection. This noise model assumes a constant centroiding performance for all Hipparcos stars, identical to the centroiding performance for the brightest ungated stars at magnitude 13. The typical along-scan standard error due to photon statistics is 94 μas. A second noise component is added to account for various effects, such as attitude modelling errors(Risquez et al. 2013)and uncertainties originating from geometrical calibration parameters of the spacecraft. Although this additional noise component may be correlated between individual CCD observations, we model it by quadratically adding a conservative RMS value of 300 μas to the photon statistical standard error per CCD.

Based on the current Gaia data release scenario⁵ we assume that the HTPM project will initially be based on one year of Gaia data. The simulation results presented in Sect. 4 use one year of Gaia observations centred around the adopted reference epoch J2015.0.

4. Results

4.1. Astrometric solution scenarios

Table 1 gives an overview of the four different solution scenarios investigated in this paper. The two cases called Gaia 12 do not use any prior data from the Hipparcos Catalogue, but only the 12 months of Gaia observations. The other two, called HTPM, use the Hipparcos covariances and astrometric parameters as priors in the processing of the same Gaia observations as in Gaia 12. A comparison between the HTPM and Gaia 12 scenarios thus allows one to assess the improvement brought by the Hipparcos prior information.

The scenarios are subdivided into cases A and B. In case A we assume that there is sufficient Gaia data to perform a full five-parameter astrometric solution for all stars even without the Hipparcos prior. This is an optimistic assumption, since in reality one year of data is only barely sufficient for a five-parameter solution under ideal conditions, i.e., without data gaps. Dead time as outlined before and the actual temporal distribution of observations over the year could mean that the solution must be constrained to estimate only the two positional parameters for most of the stars. We simulate this in case B by conservatively assuming that all stars for which we do not include a prior will have a two-parameter solution. In such a solution the parallaxes and proper motions are effectively assumed to be zero, which gives a large additional error component in the estimated positions⁶. While the Gaia,12-B solution is then restricted to two parameters for all stars, HTPM-B can still solve all five parameters of the Hipparcos stars. Case B might be closer to the foreseen first release of Gaia data and the first release of HTPM. Case A on the other hand demonstrates the capabilities of Gaia and HTPM once sufficient data for a full astrometric solution are available in subsequent releases of Gaia data.

4.2. Predicted astrometric accuracies of HTPM

Table 2 summarizes the results for the entire set of Hipparcos stars, and subdivided by magnitude. No results are given for the auxiliary stars, but they are similar to the results for the Hipparcos stars in the Gaia 12 scenarios. For comparison we also give the formal uncertainties from the Hipparcos Catalogue. For the positions they are given both at the original epoch J1991.25 and at the epoch J2015 of the Gaia data. It should be noted that the simulations include stars which in the Hipparcos Catalogue are described with more than five parameters, but are here treated as single stars. Excluding them from the statistics would systematically reduce the Hipparcos uncertainties in Table 2. The real HTPM solution will also include all Hipparcos stars independent of the type of solution in the Hipparcos Catalogue. A poor fit between the Gaia and Hipparcos data will then be used to filter out binary candidates for further treatment.

All Gaia 12 and HTPM uncertainties in Table 2 are derived from the distribution of the actual errors (calculated values minus true values) obtained in the solutions, using the robust scatter estimate (RSE)⁷. Rather than stating the uncertainty of α and δ separately we give the mean of the RSE in the two coordinates as the position uncertainty. Similarly the proper motion uncertainty is the mean RSE of the errors in μ_{α ∗} and μ_δ.

Proper motion.

The joint solution shows a big improvement in the proper motion uncertainties compared with the Hipparcos data. The improvement factor of HTPM compared with Hipparcos alone is 32 in case A and 25 in case B. The factors are similar because the Hipparcos position uncertainty dominates over the Gaia uncertainty in both cases. In the optimistic case A, the proper motions from the Gaia-only data are already better than Hipparcos alone, but not as good as the joint HTPM solution.

Using Eq. (3) to estimate the expected precision of the conventional combination we find in case A proper motions of 16 and 137 μas yr^-1 for the brightest and faintest magnitude bins, compared with 14 and 94 μas yr^-1 in the HTPM-A results. In case B we find 143 and 602 μas yr^-1, respectively, compared with 27 and 134 μas yr^-1 in HTPM-B. The joint solution thus gives consistently better results, as discussed in Sect. 2.4.

Fig. 2 Parallax errors in the HTPM solution for two cases. Bin width is 20 μas. In case A (full five-parameter astrometric solution for all stars, red/right histogram) the parallax errors are unbiased. In case B (two-parameter solution of the auxiliary stars, blue/left histogram) the median parallax error is − 591 μas.

Parallax.

The improved proper motions allow better to disentangle the five parameters in the joint astrometric solution (cf. Fig. 3), resulting in improved parallax uncertainties. In case A we find that the parallax uncertainties in the joint solution improve by a factor 23 compared with Hipparcos, and a factor 2 compared with Gaia 12. However, in the more realistic case B the improvement is much smaller (a factor 3 compared with Hipparcos) and the parallaxes are strongly biased as shown in Fig. 2. This bias originates from the assumption of zero parallax and proper motion in the two-parameter solution of the auxiliary stars. The true positive parallaxes result in a biased attitude, which propagates into the five-parameter solution of the Hipparcos stars making their parallaxes systematically too small. (As discussed in Sect. 5.2, this bias can be entirely avoided in later releases of Gaia data through a proper selection of primary sources.)

Position.

The extremely good Gaia observations lead to an improvement by up to a factor ~600 compared with Hipparcos positions propagated to J2015. In case A the slight improvement in the HTPM positions compared with Gaia 12 comes from the better determination of proper motion and parallax. In case B the Gaia-only positions show a high uncertainty due to the two-parameter solution which neglects the true parallaxes and proper motions of the stars. The increase in position uncertainties is especially pronounced for the fainter stars due to preferential selection of nearby high-proper motion stars in the non-survey part of the Hipparcos Catalogue, which means that their (neglected) parallaxes and proper motions are statistically much larger than for the brighter (survey) stars. In the HTPM solution for case B all five parameters are solved for the Hipparcos stars, so the sizes of their parallaxes and proper motions have no direct impact on the accuracy of the solution. However, the positional uncertainties are still much increased compared with case A, because the two-parameter solutions for the auxiliary stars degrade the attitude estimate.

Fig. 3 Distribution of the parallax and proper motion errors on a Hammer-Aitoff equatorial projection of the sky. All maps are for case A (full five-parameter solutions for all stars). Left figures: results from the 12 months’ Gaia-only simulation. Some regions of the sky are poorly observed resulting in zonal errors. Right figures: HTPM results for the same stars. The prior helps to disentangle proper motion and parallax, therefore we find a more homogeneous distribution of errors at an overall lower level. The cyan line follows the ecliptic for reference.

4.3. Goodness of fit statistics

Fig. 4 Left column: Goodness of fit values ΔQ for case A simulations. From top to bottom, the ΔQ values (grey bars) follow a χ2 distribution (red line) with five degrees of freedom. If the assumed radial velocities in the solution equal the true values, the actual and expected distributions agree perfectly. If the assumed radial velocity is unknown (set to zero) deviations from the expected distributions are seen. These outliers are caused by perspective acceleration. The markers in the quantile-quantile and scatter plots correspond to stars with radial velocities from XHIP (black dots) and to stars with random radial velocities (red crosses). The three rightmost red crosses in the scatter plots correspond to HIP 80190, HIP 80194, and HIP 67694 which have very large uncertainties in the Hipparcos Catalogue. Therefore they do not show a large ΔQ value even though they have large perspective acceleration. Right column: same plots for case B simulations (see Sect. 2.5).

As discussed in Sect. 2.5 the goodness of fit value ΔQ from Eq. (10) describes how well the joint astrometric solution fits the individual observations of both missions together. If all the observations are consistent with the kinematic model, then ΔQ is expected to follow a χ² distribution with five degrees of freedom. Larger values indicate deviations from the model, for example non-uniform motion caused by invisible companions or astrometric binaries. In the present simulations we do not include any such objects, so we expect ΔQ to follow the theoretical distribution.

The top two diagrams in the left column in Fig. 4 shows that this is indeed true in case A, if the radial velocities assumed in the solution are the true ones. The result would have been the same if the assumed radial velocities had only been wrong by a few km s^-1. If instead we assume zero radial velocities for all stars, as was done in the bottom two diagrams (while the observations were still generated with non-zero radial velocities), we find a small number of outliers. It turns out that all of them are nearby, high-velocity stars (Table 3) expected to show significant perspective acceleration, that is the change in proper motion due to the changing stellar distance and the changing angle between the line of sight and motion of the star ( Schlesinger 1917 ;van de Kamp 1977 ; Murray 1983). This perspective acceleration is not taken into account in the solution when the radial velocities are assumed to be zero, giving a mismatch between the Hipparcos data and the observed Gaia position. The positional offset due to the perspective acceleration after Δt years amounts to $\mathrm{\Delta }{\mathit{\theta }}_{\mathrm{persp}}\mathrm{=}\mathit{\mu }{\mathit{\mu }}_{\mathit{r}}\mathrm{\Delta }{\mathit{t}}^{\mathrm{2}}\mathit{,}$(20)where $\mathit{\mu }\mathrm{=}\mathrm{(}{\mathit{\mu }}_{\mathit{\alpha }\mathrm{\ast }}^{\mathrm{2}}\mathrm{+}{\mathit{\mu }}_{\mathit{\delta }}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}$ is the total proper motion. As shown in Table 3, the stars with a large ΔQ also have a large offset Δθ_persp at the Gaia epoch, compared with the positional uncertainty of the solution at that epoch.

This demonstrates that knowledge of radial velocities is required for a number of stars to avoid false positives in the detection of non-uniform space motion(de Bruijne & Eilers 2012) . It also shows that ΔQ is a useful statistic for detecting non-uniform space motion in general.

The right column in Fig. 4 shows the corresponding results in case B. Here ΔQ follows a scaled version of the expected distribution with a somewhat extended tail. The two bottom panels show that ΔQ is still a useful measure of deviations from the adopted kinematic model although it is much less sensitive than in case A. As a result only two outliers due to the perspective acceleration are found if the assumed radial velocities are set to zero. This demonstrates the strong dependency of ΔQ on the quality of the Gaia solution.

5. Discussion

5.1. Longevity of the HTPM solution: detection of binary and exoplanetary candidates

As Gaia collects further data the accuracy of the proper motions determined from Gaia data alone will eventually supersede that of HTPM. Assuming nominal mission performance and that the proper motion uncertainty scales with mission length as L^-1.5, this will happen already after 2–3 years of Gaia data have been accumulated. Still, HTPM will remain a valuable source of information as it is based on a much longer time baseline. This is relevant for long period companions which create astrometric signatures that cannot be seen in Gaia data alone. We therefore suggest that HTPM should be repeated with future Gaia releases. The goodness-of-fit of the combined solution is sensitive to small deviations of the stellar motions from the assumed (rectilinear) model. This sensitivity will dramatically increase with more Gaia data, namely when the Gaia-only proper motions become as good as the combined HTPM proper motions.

The potential for detecting faint (stellar or planetary) companions to nearby stars can be illustrated by a numerical example. Consider a 1 M_⊙ star at 10 pc distance (ϖ = 100 mas) from the Sun, with an invisible companion of mass m orbiting at a period of P ≃ 25 years (semi-major axis a ≃ 8.5 au). The astrometric signature of the companion (i.e., the angular size of the star’s orbit around their common centre of mass;Perryman 2014 ) is a_∗ ≃ aϖ(m/M_⊙) ≃ 850(m/M_⊙) mas if the orbit is seen face-on, and the instantaneous proper motion of the star relative to the centre of mass is 2πa_∗/P ≃ 200(m/M_⊙) mas yr^-1. If Hipparcos effectively measures this instantaneous proper motion which is extrapolated over Δt = 25 years, the extrapolated position from Hipparcos (with its uncertainty of about 22 mas, see Table 2) and the position observed by Gaia (with an uncertainty much lower than from Hipparcos) could differ by up to ≃ 5000(m/M_⊙) mas. Assuming that detection is possible if the position difference is at least twice as large as the positional uncertainty⁸, we find that the initial HTPM results could be sensitive to companion masses down to ≃ 10^-2M_⊙, that is brown dwarf or super-Jupiter companions.

If we instead let Gaia measure the instantaneous proper motion of the system and propagate backwards to the Hipparcos epoch, we can take advantage of the much better uncertainties of the Gaia astrometry. Two to three years of Gaia data already give proper motion uncertainties better than 30 μas yr^-1 for the bright stars, and hence extrapolated position uncertainties better than Hipparcos at its own epoch, or ≃ 0.75 mas (Table 2). Therefore the HTPM sensitivity increases roughly by a factor 30, allowing the detection of companion masses down to about 3 × 10^-4M_⊙, or Saturn-type objects at a Saturn-like distance to the host star.

This demonstrates that the results of HTPM can be used to find candidates for long period exoplanets around nearby stars, with a highly interesting companion mass range opening up with subsequent releases of Gaia data when combined with Hipparcos. These companions cannot be detected from Gaia data alone even at the full mission length, and are hard to detect through classical methods due to their long periods, low transit probability and small radial velocity signatures. Since ΔQ is sensitive to deviations from uniform space motion, whether they are seen in the Hipparcos or in the Gaia data, or both, this statistic can be used to find candidate systems in all these cases. The further exploration of the candidate systems will, however, require specialised analysis tools.

In a future publication we will explore in more detail how ΔQ can be used to identify binary and exoplanetary candidates with orbital periods of decades to centuries. Apart from the possibility to detect substellar companions for the nearest stars, this will contribute to the census of the binary population within a few hundred parsecs from the sun by filling a difficult-to-observe gap between the shorter period spectroscopic and astrometric binaries and the visually resolved long-period systems.

5.2. Two versus five parameters

When evaluating the results of our simulations, case B deserves additional attention since it is the more realistic case for the first Gaia data release, and the first simulation of this case published so far. The two-parameter solution (Gaia,12-B in Table 2) leads to a large position error of several mas. This is caused by assuming the parallax, proper and radial motion to be zero in the solution, whereas in reality they are not. The actual positional uncertainties in this case depend on the true distribution of parallaxes and proper motions for all the stars, including the auxiliary stars, which are not very well known. The numerical values given here are based on the very schematic distribution model for the auxiliary stars described in Sect. 3.2.1, and should therefore be interpreted with caution.

This position error is also relevant for the case B HTPM scenario, where the solution of the auxiliary stars is two parameters only, but where one solves all five parameters for the Hipparcos stars while incorporating prior information from the Hipparcos Catalogue. The position error of the auxiliary stars causes a poor attitude determination. This in turn leads to increased errors in the case B HTPM results (compare HTPM B and A in Table 2), with a bias in the parallax errors (see Sect. 4 and Fig. 2). For a parallax-unbiased solution it is necessary to estimate all five parameters for all stars included in the solution. Any mixture in the estimation of five and two parameters in the same solution will lead to a bias in the resulting parallaxes. This is not only true for the HTPM scenario described in this paper but also in all Gaia-only data releases. Referring to the terminology used in Sect. 6.2 ofLindegren et al. (2012) , any star for which not all five astrometric parameters can be solved must be treated as a “secondary source”, meaning that it does not contribute to the attitude determination and instrument calibration. This is necessary in order to avoid biases for the stars where all five parameters are estimated.

5.3. Frame rotation of the combined solution

For the final AGIS solution of Gaia the reference frame will be established by means of quasars, both by linking to the optical counterparts of radio (VLBI) sources defining the orientation of the International Celestial Reference Frame, and by using the zero proper motion of quasars to determine a non-rotating frame⁹. This can also be done for earlier Gaia data releases, at least for the orientation part, while the shorter time span will limit the determination of the spin. It is desirable to rotate the HTPM results into the same reference frame as used for the first Gaia data release. This must be done in two steps. First, a provisional HTPM must be computed in the Hipparcos frame (as it will be when the Hipparcos data are used as prior, see Sect. 2.7), without imposing any other constraints on the frame. This solution will contain (many) non-Hipparcos stars with only Gaia observations which include a multitude of quasars. Their positions and proper motions are used in a second step to correct the provisional HTPM (and other data in the same solution) for the estimated orientation and spin. Since the HTPM solution is integrated in AGIS, the estimation and correction of the frame can be accomplished using the procedures and tools developed for AGIS ( Lindegren et al. 2012 , Sect. 6.1).

5.4. Other applications of the joint solution method

The joint solution is applicable also to other combinations of astrometric data. Here we give two examples.

Nano-JASMINE ( Hatsutori et al. 2009 ;Yamada et al. 2013 ) is an ultra-small Japanese satellite, a technology demonstrator for the JASMINE series of near-infrared astrometry missions, scheduled for launch in 2015. It targets bright stars between magnitude 1 and 10, although the exact limits are not yet determined. Based on current performance estimates the uncertainties in stellar parameters will be similar to or slightly worse than the uncertainties of the Hipparcos data. However, the data will still be very valuable since astrometric catalogues are best at their respective epochs and Nano-JASMINE may be the only astrometric mission at its epoch observing the brightest stars in the sky. The Nano-JASMINE data can be analysed together with Hipparcos data analogously to the HTPM project to improve the proper motions of bright stars that may not be observed by Gaia(Michalik et al. 2013) .

The Tycho-2 Catalogue (Høg et al. 2000) gives positions for 2.5 million stars, derived from starmapper observations of Hipparcos. The positions at the reference epoch J1991.25 have a median internal standard error of 7 mas for stars brighter than V_T = 9 mag and 60 mas for the whole catalogue. Combining the Tycho-2 positions with Gaia data using the joint solution scheme would allow us to derive proper motions for these stars with median uncertainties of 0.3 and 2.5 mas yr^-1, respectively. This is true even in the conservative scenario (Gaia,12-B), since the major uncertainty comes from the Tycho-2 positions. In this combination the proper motions given in Tycho-2 should not be used, as they may contain systematic errors of a similar magnitude due to the incorporated old photographic material. However, the derived proper motions from a Tycho-Gaia Astrometric Solution (TGAS) could be used to correct the photographic positions in order to take advantage of a much longer temporal baseline.

6. Conclusions

We have developed the joint solution method for incorporating priors in the core astrometric solution of Gaia. The method can be used in the processing of early Gaia data to improve the proper motions of the Hipparcos stars, the so-called Hundred Thousand Proper Motions project.

Combining astrometric data from very different epochs requires careful treatment of the non-linear effects of the mapping from spherical to rectilinear coordinates and for high velocity stars due to perspective acceleration. Therefore we have introduced a scaled model of kinematics (SMOK), which allows one to handle these effects in a simple and rigorous manner.

Using simulations we have verified that HTPM, using the joint solution method, gives the expected large improvements in proper motion uncertainties for over 100 000 stars in the Hipparcos Catalogue. The predicted proper motion uncertainties range from 14 to 134 μas yr^-1 depending on the amount of Gaia data used and the stellar magnitude, about a factor 30 improvement compared with the Hipparcos uncertainties.

We have shown that HTPM also delivers improved parallaxes, which, however, may be strongly biased unless a full five-parameter solution can be obtained from Gaia-only data also for all non-Hipparcos stars. Whether these parallaxes should be published as part of an HTPM release should be decided based on the amount and quality of Gaia data available at the time.

The joint solution is applicable also to a combination of Tycho-2 positions with early Gaia data to derive parallaxes and improved proper motions for the 2.5 million stars. We suggest that this possibility of a Tycho-Gaia Astrometric Solution (TGAS) should be considered in the Gaia data release plan.

The proposed method to calculate HTPM provides a goodness-of-fit measurement ΔQ which is sensitive to deviations from the uniform linear space motion. However, accurate radial velocities are required for nearby fast moving stars in order

to avoid mistaking outliers in ΔQ for companion signatures. We recommend to publish ΔQ as well as the radial velocities used for the HTPM data reduction. This will allow further investigations of outliers which might indicate binary or exoplanetary candidates, and will permit a correction of the HTPM results if better radial velocities become available.

The full power of HTPM will not be reached with the first Gaia data, but only in subsequent releases benefiting from the increased sensitivity of ΔQ with improved Gaia results. Because of the long temporal baseline and the combination of current with historic astrometry, HTPM will remain relevant throughout the final Gaia release for the detection and measurement of binary and exoplanetary candidates.

Online material

Appendix A: Scaled modelling of kinematics (SMOK)

A formalism called SMOK is introduced in this paper to facilitate a rigorous manipulation of small (differential) quantities in the celestial coordinates. It is reminiscent of the “standard” or “tangential” coordinates in classical small-field astrometry (e.g., Murray 1983; van Altena 2013), using a gnomonic projection onto a tangent plane of the (unit) celestial sphere, but extends to three dimensions by adding the radial coordinate perpendicular to the tangent plane. This simplifies the modelling of perspective effects.

Figure A.1 illustrates the concept. In the vicinity of the star let c be a comparison point fixed with respect to the solar system barycentre (SSB). As shown in the diagrams:

• 1.

  The barycentric motion of the star is scaled by the inverse distance to c, effectively placing the star on or very close to the unit sphere.

• 2.

  Rectangular coordinates are expressed in the barycentric [ p_cq_cr_c ] system with r_c pointing towards c, and p_c, q_c in the directions of increasing right ascension and declination.

The first point eliminates the main uncertainty in the kinematic modelling of the star due to its poorly known distance. The second point allows us to express the scaled kinematic model in SMOK coordinates a, d, r that are locally aligned with α, δ, and the barycentric vector.

Fig. A.1 Two steps in the definition of SMOK coordinates. In the top diagram the motion of an object in the vicinity of the fixed point c is modelled by the function b(t) expressed in the barycentric [ xyz ] system. A scaled version of the model is constructed such that the scaled c is at unit distance from the solar system barycentre (SSB). In the bottom diagram new coordinate axes [ pcqcrc ] are chosen in the directions of increasing right ascension, declination, and distance, respectively, at the comparison point (αc,δc) being the projection of c on the unit sphere.

Up to the scale factor | c | ^-1 discussed below, the SMOK coordinate system is completely defined by the adopted comparison point (α_c,δ_c) using the orthogonal unit vectors ${{p}}_{\mathrm{c}}\mathrm{=}\left[\begin{array}{c}\\ \mathrm{-}\sin {\mathit{\alpha }}_{\mathrm{c}}\\ \cos {\mathit{\alpha }}_{\mathrm{c}}\\ \mathrm{0}\end{array}\right]\mathit{,} {{q}}_{\mathrm{c}}\mathrm{=}\left[\begin{array}{c}\\ \mathrm{-}\sin {\mathit{\delta }}_{\mathrm{c}}\cos {\mathit{\alpha }}_{\mathrm{c}}\\ \mathrm{-}\sin {\mathit{\delta }}_{\mathrm{c}}\sin {\mathit{\alpha }}_{\mathrm{c}}\\ \cos {\mathit{\delta }}_{\mathrm{c}}\end{array}\right]\mathit{,} {{r}}_{\mathrm{c}}\mathrm{=}\left[\begin{array}{c}\\ \cos {\mathit{\delta }}_{\mathrm{c}}\cos {\mathit{\alpha }}_{\mathrm{c}}\\ \cos {\mathit{\delta }}_{\mathrm{c}}\sin {\mathit{\alpha }}_{\mathrm{c}}\\ \sin {\mathit{\delta }}_{\mathrm{c}}\end{array}\right]\mathit{.}$(A.1)[ p_cq_cr_c ]is the “normal triad” at the comparison point with respect to the celestial coordinate system (Murray 1983)¹⁰. We are free to choose (α_c,δ_c) as it will best serve our purpose, but once chosen (for a particular application) it is fixed: it has no proper motion, no parallax, and no associated uncertainty. Typically (α_c,δ_c) is chosen very close to the mean position of the star.

The motion of the star in the Barycentric Celestial Reference System (BCRS) is represented by the function b(t), where b is the vector from SSB to the star as it would be observed from the SSB at time t. The scaled kinematic model s(t) = b(t) | c | ^-1 is given in SMOK coordinates as $\mathit{a}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{{p}}_{\mathrm{c}}^{\mathrm{\prime }}{s}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{,} \mathit{d}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{{q}}_{\mathrm{c}}^{\mathrm{\prime }}{s}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{,} \mathit{r}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{{r}}_{\mathrm{c}}^{\mathrm{\prime }}{s}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{,}$(A.2)and can in turn be reconstructed from the SMOK coordinates as ${s}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{{p}}_{\mathrm{c}}\mathit{a}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{+}{{q}}_{\mathrm{c}}\mathit{d}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{+}{{r}}_{\mathrm{c}}\mathit{r}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{.}$(A.3)a, d, r are dimensionless and the first two are typically small quantities (≲ 10^-4), while r is very close to unity.

The whole point of the scaled kinematic modelling is that s(t) can be described very accurately by astrometric observations, even though b(t) may be poorly known due to a high uncertainty in distance. This is possible simply by choosing the scaling such that | s(t) | = 1 at some suitable time. This works even if the distance is completely unknown, or if it is effectively infinite (as for a quasar).

The scale factor is | c | ^-1 = ϖ_c/A, where ϖ_c is the parallax of c and A the astronomical unit. The measured parallax can be regarded as an estimate of ϖ_c.

In the following we describe some typical applications of SMOK coordinates.

Appendix A.1: Uniform space motion

The simplest kinematic model is to assume that the star moves uniformly with respect to the SSB, that is ${b}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{{b}}_{\mathrm{ep}}\mathrm{+}\mathrm{(}\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{ep}}\mathrm{)}{v}\mathit{,}$(A.4)where b_ep is the barycentric position at the reference epoch t_ep, and v is the (constant) space velocity. The scaled kinematic model expressed in the BCRS is ${s}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{{s}}_{\mathrm{ep}}\mathrm{+}\mathrm{(}\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{ep}}\mathrm{)}{ṡ}\mathit{,}$(A.5)where ${{s}}_{\mathrm{ep}}\mathrm{=}{{p}}_{\mathrm{c}}\mathit{a}\mathrm{\left(}{\mathit{t}}_{\mathrm{ep}}\mathrm{\right)}\mathrm{+}{{q}}_{\mathrm{c}}\mathit{d}\mathrm{\left(}{\mathit{t}}_{\mathrm{ep}}\mathrm{\right)}\mathrm{+}{{r}}_{\mathrm{c}}\mathit{r}\mathrm{\left(}{\mathit{t}}_{\mathrm{ep}}\mathrm{\right)}\hspace{0.17em}$(A.6)and ${ṡ}\mathrm{=}{{p}}_{\mathrm{c}}\mathit{ȧ}\mathrm{+}{{q}}_{\mathrm{c}}\mathit{ḋ}\mathrm{+}{{r}}_{\mathrm{c}}\mathit{ṙ}\hspace{0.17em}$(A.7)are constant vectors. The uniform motion can also be written in SMOK coordinates as $\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\mathit{a}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\mathit{a}\mathrm{\left(}{\mathit{t}}_{\mathrm{ep}}\mathrm{\right)}\mathrm{+}\mathrm{(}\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{ep}}\mathrm{)}\mathit{ȧ}\mathit{,}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}\mathit{d}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\mathit{d}\mathrm{\left(}{\mathit{t}}_{\mathrm{ep}}\mathrm{\right)}\mathrm{+}\mathrm{(}\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{ep}}\mathrm{)}\mathit{ḋ}\mathit{,}\\ \begin{array}{c}\end{array}\mathit{r}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\mathit{r}\mathrm{\left(}{\mathit{t}}_{\mathrm{ep}}\mathrm{\right)}\mathrm{+}\mathrm{(}\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{ep}}\mathrm{)}\mathit{ṙ}\mathit{.}\end{array} \begin{array}{c}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}$(A.8)The six constants a(t_ep), d(t_ep), r(t_ep), ȧ, ḋ, ṙ are the kinematic parameters of the scaled model; however, to get the actual kinematics of the star we also need to know ϖ_c.

Appendix A.2: Relation to the usual astrometric parameters

Choosing (α_c,δ_c) to be the barycentric celestial coordinates of the star at t_ep, and ϖ_c equal to the parallax at the same epoch, we find $\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\mathit{a}\mathrm{\left(}{\mathit{t}}_{\mathrm{ep}}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathit{,}\begin{array}{c}\end{array}\mathit{d}\mathrm{\left(}{\mathit{t}}_{\mathrm{ep}}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathit{,}\begin{array}{c}\end{array}\mathit{r}\mathrm{\left(}{\mathit{t}}_{\mathrm{ep}}\mathrm{\right)}\mathrm{=}\mathrm{1}\mathit{,}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}\mathit{ȧ}\mathrm{=}{\mathit{\mu }}_{\mathit{\alpha }\mathrm{\ast }}\mathit{,}\begin{array}{c}\end{array}\mathit{ḋ}\mathrm{=}{\mathit{\mu }}_{\mathit{\delta }}\mathit{,}\begin{array}{c}\end{array}\mathit{ṙ}\mathrm{=}{\mathit{\mu }}_{\mathit{r}}\mathit{,}\end{array}\}$(A.9)where μ_{α ∗}, μ_δ are the tangential components of the barycentric proper motion at the reference epoch t_ep, and μ_r is the “radial proper motion” allowing one to take the perspective effects into account. μ_r is usually calculated from the measured radial velocity and parallax according to Eq. (1).

Appendix A.3: Differential operations

Uniform space motion does not map into barycentric coordinates α(t), δ(t) that are linear functions of time. The non-linearity derives both from the curvilinear nature of spherical coordinates and from perspective foreshortening depending on the changing distance to the object. Both effects are well known and have been dealt with rigorously by several authors (e.g., Eichhorn & Rust 1970; Taff 1981). The resulting expressions are non-trivial and complicate the comparison of astrometric catalogues of different epochs. For example, approximations such as ${\mathit{\mu }}_{\mathit{\alpha }\mathrm{\ast }}\mathrm{=}\frac{\mathit{\alpha }\mathrm{\left(}{\mathit{t}}_{\mathrm{2}}\mathrm{\right)}\mathrm{-}\mathit{\alpha }\mathrm{\left(}{\mathit{t}}_{\mathrm{1}}\mathrm{\right)}}{{\mathit{t}}_{\mathrm{2}}\mathrm{-}{\mathit{t}}_{\mathrm{1}}}\cos \mathit{\delta ,} {\mathit{\mu }}_{\mathit{\delta }}\mathrm{=}\frac{\mathit{\delta }\mathrm{\left(}{\mathit{t}}_{\mathrm{2}}\mathrm{\right)}\mathrm{-}\mathit{\delta }\mathrm{\left(}{\mathit{t}}_{\mathrm{1}}\mathrm{\right)}}{{\mathit{t}}_{\mathrm{2}}\mathrm{-}{\mathit{t}}_{\mathrm{1}}}$(A.10)cannot be used when the highest accuracy is required. By contrast, the linearity of Eq. (A.8) makes it possible to write $\mathit{ȧ}\mathrm{=}\frac{\mathit{a}\mathrm{\left(}{\mathit{t}}_{\mathrm{2}}\mathrm{\right)}\mathrm{-}\mathit{a}\mathrm{\left(}{\mathit{t}}_{\mathrm{1}}\mathrm{\right)}}{{\mathit{t}}_{\mathrm{2}}\mathrm{-}{\mathit{t}}_{\mathrm{1}}}\mathit{,} \mathit{ḋ}\mathrm{=}\frac{\mathit{d}\mathrm{\left(}{\mathit{t}}_{\mathrm{2}}\mathrm{\right)}\mathrm{-}\mathit{d}\mathrm{\left(}{\mathit{t}}_{\mathrm{1}}\mathrm{\right)}}{{\mathit{t}}_{\mathrm{2}}\mathrm{-}{\mathit{t}}_{\mathrm{1}}}$(A.11)to full accuracy, provided that the same comparison point is used for both epochs. (Strictly speaking, the same scale factor must also be used, so that in general r(t₂) − r(t₁) = (t₂ − t₁)ṙ ≠ 0.) If the position at the reference epoch coincides with the comparison point used, the resulting ȧ, ḋ are the looked-for proper motion components according to Eq. (A.9); otherwise a change of comparison point is needed (see below).

Appendix A.4: Changing the comparison point

Let (α₁,δ₁) and (α₂,δ₂) be different comparison points with associated triads [ p₁q₁r₁ ] and [ p₂q₂r₂ ]. If a₁(t), d₁(t), r₁(t) and a₂(t), d₂(t), r₂(t) describe the same scaled kinematics we have by Eq. (A.3) ${s}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{{p}}_{\mathrm{1}}{\mathit{a}}_{\mathrm{1}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{+}{{q}}_{\mathrm{1}}{\mathit{d}}_{\mathrm{1}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{+}{{r}}_{\mathrm{1}}{\mathit{r}}_{\mathrm{1}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{{p}}_{\mathrm{2}}{\mathit{a}}_{\mathrm{2}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{+}{{q}}_{\mathrm{2}}{\mathit{d}}_{\mathrm{2}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{+}{{r}}_{\mathrm{2}}{\mathit{r}}_{\mathrm{2}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{.}$(A.12)Thus, given a₁(t), d₁(t), r₁(t) one can compute s(t) from the first equality in Eq. (A.12), whereupon the modified functions are recovered as ${\mathit{a}}_{\mathrm{2}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{{p}}_{\mathrm{2}}^{\mathrm{\prime }}{s}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{,} {\mathit{d}}_{\mathrm{2}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{{q}}_{\mathrm{2}}^{\mathrm{\prime }}{s}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{,} {\mathit{r}}_{\mathrm{2}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{{r}}_{\mathrm{2}}^{\mathrm{\prime }}{s}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{.}$(A.13)This procedure can be applied to s(t) for any particular t as well as to linear operations on s such as differences and time derivatives.

Appendix A.5: Epoch propagation

An important application of the above formulae is for propagating the six astrometric parameters (α₁,δ₁,ϖ₁,μ_{α ∗ 1},μ_δ1,μ_r1), referring to epoch t₁, to a different epoch t₂. This can be done in the following steps:

• 1.

  Use (α₁,δ₁) as the comparison point and compute [ p₁q₁r₁ ] by Eq. (A.1). At time t₁ the SMOK parameters relative to the first comparison point are a₁(t₁) = d₁(t₁) = 0, r₁(t₁) = 1, ȧ₁ = μ_{α ∗ 1}, ḋ₁ = μ_δ1, ṙ₁ = μ_r1.

• 2.

  Calculate s(t₁) and ṡ using Eqs. (A.6)–(A.7).

• 3.

  Calculate s(t₂) by means of Eq. (A.5). Let s₂ = | s(t₂) | be its length (close to unity).

• 4.

  Calculate r₂ = s(t₂) /s₂ and hence the second comparison point (α₂,δ₂) and triad [ p₂q₂r₂ ].

• 5.

  Use Eq. (A.13) to calculate the SMOK parameters at t₂ referring to the second comparison point. For the position one trivially gets a₂(t₂) = d₂(t₂) = 0 and r₂(t₂) = s₂. For the proper motion parameters one finds $\mathit{ȧ}{}_{\mathrm{2}}\mathrm{=}{{p}}_{\mathrm{2}}^{\mathrm{\prime }}{ṡ}$, $\mathit{ḋ}{}_{\mathrm{2}}\mathrm{=}{{q}}_{\mathrm{2}}^{\mathrm{\prime }}{ṡ}$, and $\mathit{ṙ}{}_{\mathrm{2}}\mathrm{=}{{r}}_{\mathrm{2}}^{\mathrm{\prime }}{ṡ}$.

• 6.

  The astrometric parameters at epoch t₂ are α₂, δ₂, ϖ₂ = ϖ₁/s₂, μ_{α ∗ 2} = ȧ₂/s₂, μ_δ2 = ḋ₂/s₂, μ_r2 = ṙ₂/s₂.

This procedure is equivalent to the one described in Sect. 1.5.5, Vol. 1 of The Hipparcos and Tycho Catalogues (ESA 1997).

Appendix B: The HIPPARCOS Catalogue

This Appendix describes the calculation of relevant quantities from the new reduction of the Hipparcos Catalogue by van Leeuwen (2007b). Data files were retrieved from the Strasbourg astronomical Data Center (CDS) in November 2013 (catalogue I/311). These files differ slightly from the ones given on the DVD published along with the book (van Leeuwen 2007a), both in content and format, as some errors have been corrected. The data needed for every accepted catalogue entry are:

• the five astrometric parameters (α,δ,ϖ,μ_{α ∗},μ_δ);

• the 5 × 5 normal matrix N from the least-squares solution of the astrometric parameters (for a 5-parameter solution this equals the inverse of the covariance matrix C);

• the chi-square goodness-of-fit quantity Q for the 5-parameter solution of the Hipparcos data;

• the degrees of freedom ν associated with Q.

The astrometric parameters at the Hipparcos reference epoch J1991.25 are directly taken from the fields labelled RArad, DErad, Plx, pmRA, and pmDE in the main catalogue file hip2.dat. Units are [rad] for α and δ, [mas] for ϖ, and [mas yr^-1] for μ_{α ∗} and μ_δ. It is convenient to express also positional differences (such as SMOK coordinates a and d) and positional uncertainties in [mas]. The elements of N thus have units [mas^-2 yr^p], where p = 0, 1, or 2, depending on the position of the element in the matrix.

The calculation of N, Q, and ν is described hereafter in some detail as the specification of C deviates in some details from the published documentation. Clarification on certain issues was kindly provided by van Leeuwen (priv. comm.).

The number of degrees of freedom is $\mathit{\nu }\mathrm{=}{\mathit{N}}_{\mathrm{tr}}\mathrm{-}\mathit{n,}$(B.1)where N_tr is the number of field transits used (label Ntr in hip2.dat) and n is the number of parameters in the solution (see below; most stars have n = 5). The goodness-of-fit given in field F2 is the “gaussianized” chi-square(Wilson & Hilferty 1931)${\mathit{F}}_{\mathrm{2}}\mathrm{=}{\left(\frac{\mathrm{9}\mathit{\nu }}{\mathrm{2}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\left[{\left(\frac{\mathit{Q}}{\mathit{\nu }}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathrm{+}\frac{\mathrm{2}}{\mathrm{9}\mathit{\nu }}\mathrm{-}\mathrm{1}\right]$(B.2)computed from Q, the sum of the squared normalized residuals, and ν. For “good” solutions Q is expected to follow the chi-square distribution with ν degrees of freedom (Q ~ χ²(ν)), in which case F₂ approximately follows the standard normal distribution, F₂ ~ N(0,1). Thus, F₂> 3 means that Q is “too large” for the given ν at the same level of significance as the + 3σ criterion for a Gaussian variable (probability ≲ 0.0044)¹¹. Given F₂ from field F2, and ν from Eq. (B.1), it is therefore possible to reconstruct the chi-square statistic of the n-parameter solution as $\mathit{Q}\mathrm{=}\mathit{\nu }{\left[{\left(\frac{\mathrm{2}}{\mathrm{9}\mathit{\nu }}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{F}}_{\mathrm{2}}\mathrm{+}\mathrm{1}\mathrm{-}\frac{\mathrm{2}}{\mathrm{9}\mathit{\nu }}\right]}^{\mathrm{3}}\mathrm{·}$(B.3)We also introduce the square-root of the reduced chi-square, $\mathit{u}\mathrm{=}\sqrt{\mathit{Q}\mathit{/}\mathit{\nu }}\mathit{,}$(B.4)which is expected to be around 1.0 for a “good” solution (see further discussion below). u is sometimes referred to as the standard error of unit weight (Brinker & Minnick 1995).

The catalogue gives the covariance matrix in the form of an upper-diagonal “weight matrix” U such that, formally, C = (U′U)^-1. This inverse exists for all stars where a solution is given. (For the joint solution we actually need the normal matrix N = U′U, see below.) For solutions with n = 5 astrometric parameters there are n(n + 1)/2 = 15 non-zero elements in U. For some stars the solution has more than five parameters, and the main catalogue then only gives the first 15 non-zero elements, while remaining elements are given in separate tables.

Let U₁, U₂, ..., U₁₅ be the 15 values taken from the fields labelled UW in hip2.dat. The matrix U is computed as ${U}\mathrm{=}\left[\begin{array}{c}\\ {\mathit{f}}_{\mathrm{1}}{\mathit{U}}_{\mathrm{1}}& {\mathit{U}}_{\mathrm{2}}& {\mathit{U}}_{\mathrm{4}}& {\mathit{U}}_{\mathrm{7}}& {\mathit{U}}_{\mathrm{11}}\\ \mathrm{0}& {\mathit{f}}_{\mathrm{2}}{\mathit{U}}_{\mathrm{3}}& {\mathit{U}}_{\mathrm{5}}& {\mathit{U}}_{\mathrm{8}}& {\mathit{U}}_{\mathrm{12}}\\ \mathrm{0}& \mathrm{0}& {\mathit{f}}_{\mathrm{3}}{\mathit{U}}_{\mathrm{6}}& {\mathit{U}}_{\mathrm{9}}& {\mathit{U}}_{\mathrm{13}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& {\mathit{f}}_{\mathrm{4}}{\mathit{U}}_{\mathrm{10}}& {\mathit{U}}_{\mathrm{14}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& {\mathit{f}}_{\mathrm{5}}{\mathit{U}}_{\mathrm{15}}\end{array}\right]\mathit{.}$(B.5)Here f_i, i = 1...n, are scaling factors which for the CDS data must be calculated as $\begin{array}{ccc}{\mathit{f}}_{\mathrm{1}}\mathrm{=}\mathit{u}\mathit{/}{\mathit{\sigma }}_{\mathit{\alpha }}\mathit{,}& {\mathit{f}}_{\mathrm{2}}\mathrm{=}\mathit{u}\mathit{/}{\mathit{\sigma }}_{\mathit{\delta }}\mathit{,}& {\mathit{f}}_{\mathrm{3}}\mathrm{=}\mathit{u}\mathit{/}{\mathit{\sigma }}_{\mathit{\varpi }}\mathit{,}\\ {\mathit{f}}_{\mathrm{4}}\mathrm{=}\mathit{u}\mathit{/}{\mathit{\sigma }}_{\mathit{\mu \alpha }}\mathit{,}& & {\mathit{f}}_{\mathrm{5}}\mathrm{=}\mathit{u}\mathit{/}{\mathit{\sigma }}_{\mathit{\mu \delta }}\mathit{,}\end{array}$(B.6)where u is given by Eq. (B.4) and σ_· are the standard errors given in fields e_RArad through e_pmDE of hip2.dat. Equation (B.6) applies to data taken from the CDS version of the catalogue (I/311). For catalogue data on the DVD accompanying the book (van Leeuwen 2007a), scaling factors f_i = 1 apply, although those data are superseded by the CDS version.

The 5 × 5 matrix N = U′U computed using the first five rows and columns in U, as given in Eq. (B.5), contains the relevant elements of the normal matrix for any solution with n ≥ 5. Thus, for solutions with n> 5 there is no need, for the catalogue combination, to retrieve the additional elements of U from hip7p.dat, etc. The situation is different when the covariance matrix is needed: it is then necessary to compute the full n × n normal matrix N before C = N^-1 can be computed.

The normal matrix N computed as described above incorporates the formal uncertainties of the observations; as described in van Leeuwen (2007a) these are ultimately derived from the photon statistics of the raw data after careful analysis of the residuals as function of magnitude, etc. If the adopted models are correct we expect the F₂ statistic to be normally distributed with zero mean and unit standard deviation, and the standard error of unit weight, u, to be on the average equal to 1. In reality we find (for solutions with n = 5) that their distributions are skewed towards higher values, especially for the bright stars where photon noise is small and remaining calibration errors are therefore relatively more important. To account for such additional errors the published standard errors σ_{α ∗}, etc., in hip2.dat include, on a star-by-star basis, a correction factor equal to the unit weight error u obtained in its solution. This is equivalent to scaling the formal standard errors of the data used in the solution by the same factor. In order to make the computed normal matrix, covariance matrix, and goodness-of-fit statistics consistent with the published standard errors it is then necessary to apply the corresponding corrections, viz.: ${{N}}_{\mathrm{corr}}\mathrm{=}{N}{\mathit{u}}^{-2}\mathit{,} {{C}}_{\mathrm{corr}}\mathrm{=}{C}{\mathit{u}}^{\mathrm{2}}\mathit{,} {\mathit{Q}}_{\mathrm{corr}}\mathrm{=}\mathit{\nu ,} {\mathit{u}}_{\mathrm{corr}}\mathrm{=}\mathrm{1.}$(B.7)For the catalogue combination we use N_corr and Q_corr whenever u> 1, but N and Q if u ≤ 1.

Acknowledgments

We thank F. van Leeuwen for clarification on certain data items in the Hipparcos Catalogue and for providing valuable feedback as the referee. We also thank C. Fabricius for many useful comments. This work was partly carried out under ESA Contract No. 4000105564/12/NL/GE. Support from the Swedish National Space Board and the Royal Physiographic Society in Lund is gratefully acknowledged.

References

All Tables

All Figures

	Fig. 1 Relationships between catalogues during simulation runs.
In the text

	Fig. 2 Parallax errors in the HTPM solution for two cases. Bin width is 20 μas. In case A (full five-parameter astrometric solution for all stars, red/right histogram) the parallax errors are unbiased. In case B (two-parameter solution of the auxiliary stars, blue/left histogram) the median parallax error is − 591 μas.
In the text

Fig. 3 Distribution of the parallax and proper motion errors on a Hammer-Aitoff equatorial projection of the sky. All maps are for case A (full five-parameter solutions for all stars). Left figures: results from the 12 months’ Gaia-only simulation. Some regions of the sky are poorly observed resulting in zonal errors. Right figures: HTPM results for the same stars. The prior helps to disentangle proper motion and parallax, therefore we find a more homogeneous distribution of errors at an overall lower level. The cyan line follows the ecliptic for reference.

In the text

Fig. 4 Left column: Goodness of fit values ΔQ for case A simulations. From top to bottom, the ΔQ values (grey bars) follow a χ2 distribution (red line) with five degrees of freedom. If the assumed radial velocities in the solution equal the true values, the actual and expected distributions agree perfectly. If the assumed radial velocity is unknown (set to zero) deviations from the expected distributions are seen. These outliers are caused by perspective acceleration. The markers in the quantile-quantile and scatter plots correspond to stars with radial velocities from XHIP (black dots) and to stars with random radial velocities (red crosses). The three rightmost red crosses in the scatter plots correspond to HIP 80190, HIP 80194, and HIP 67694 which have very large uncertainties in the Hipparcos Catalogue. Therefore they do not show a large ΔQ value even though they have large perspective acceleration. Right column: same plots for case B simulations (see Sect. 2.5).

In the text

Fig. A.1 Two steps in the definition of SMOK coordinates. In the top diagram the motion of an object in the vicinity of the fixed point c is modelled by the function b(t) expressed in the barycentric [ xyz ] system. A scaled version of the model is constructed such that the scaled c is at unit distance from the solar system barycentre (SSB). In the bottom diagram new coordinate axes [ pcqcrc ] are chosen in the directions of increasing right ascension, declination, and distance, respectively, at the comparison point (αc,δc) being the projection of c on the unit sphere.

In the text