© ESO 2019

1. Introduction

The Gaia Celestial Reference Frame (Gaia-CRF; Gaia Collaboration 2018a) is formally defined by the positions, as measured by Gaia, of a large number of sources that are identified as quasars. Through their cosmological distances, these objects define a kinematically non-rotating reference frame, that is, their proper motions are assumed to be zero on average. A subset of them, identified as the optical counterparts of radio sources with accurate positions in the International Celestial Reference Frame (ICRF; Ma et al. 1998) from very long baseline interferometry (VLBI) observations, are used to align the axes of the non-rotating quasar frame with the ICRF. The second release of Gaia data (DR2; Gaia Collaboration 2018b) lists 556 869 quasars whose positions at epoch J2015.5 define the optical reference frame known as Gaia-CRF2. This includes 2820 sources matched to a prototype version of ICRF3 (Jacobs et al. 2018).

The vast majority of sources in Gaia DR2 are Galactic stars with sizeable proper motions. The implicit assumption is that the positions and proper motions of the stars, and indeed the barycentric coordinates of all Gaia sources, are expressed in the same reference frame as the quasars. This is fundamental for the dynamical interpretation of the observations, which assumes the absence of the inertial (Coriolis and centrifugal) forces that would appear in a rotating frame.

Although the measurement and reduction principles of the Gaia mission have been designed to provide a globally consistent reference frame for all kinds of objects, subtle differences are inevitable as a consequence of the varying conditions under which the objects are observed. For example, the quasars defining Gaia-CRF are all faint (fewer than 1% have G <  17 mag), on average bluer than stars of comparable magnitude, and they have a very different distribution on the sky than the stars. Differences in magnitude, colour, and numerous other factors are likely to produce small shifts of the image centroids, which, if left uncalibrated, may propagate into systematic errors of the positions and proper motions.

Intricate instrument models have been set up and calibrated as part of the global astrometric reductions of Gaia data in order to eliminate such systematics. In Gaia DR2 there is nevertheless an indication that the reference frame defined by the bright stars (up to G  ≃  11 − 13) rotates with respect to the quasars at a rate of a few tenths of a milliarcsecond (mas) per year. This is illustrated in Fig. 4 of Lindegren et al. (2018) and further quantified by Brandt (2018). In Brandt’s paper the rotation shows up as a systematic offset of the proper motions of bright stars in Gaia DR2 with respect to the “HIPPARCOS–Gaia proper motions” calculated from the position differences between Gaia DR2 (at epoch 2015.5) and the HIPPARCOS catalogue (at epoch 1991.25), scaled by the epoch difference of ∼24 yr. Through the long time-baseline, the HIPPARCOS–Gaia proper motions constitute a precise set of reference values, which are moreover inertial because the positional systems of HIPPARCOS and Gaia were both aligned with the International Celestial Reference System (ICRS) at their respective epoch. The observed offset therefore points to a systematic error in the Gaia DR2 proper motions of the bright stars, equivalent to an inertial rotation of its reference frame.

The cause of this rotation is discussed in Appendix B. Briefly, it is related to the different modes in which Gaia’s CCDs are operated, depending on the magnitude of the source, and corresponding differences in the calibration models. In particular, around G  =  13, there is a transition in the on-board CCD sampling scheme from two- to one-dimensional pixel windows, causing abrupt changes in the quality of both the astrometric data and the G-band photometry (which uses the same CCDs). Examples of this are shown in Figs. 9 and B.2 of Lindegren et al. (2018), and Fig. 9 of Evans et al. (2018).

It is therefore justified to study the reference frames separately that are defined by the bright and the faint Gaia sources, and to draw the division line at G  ≃  13. Gaia-CRF2, being defined by the quasars, clearly belongs to the faint part and is by construction to very high accuracy non-rotating with respect to the ICRS; its properties are discussed elsewhere (Gaia Collaboration 2018a). The subject of this paper is the bright reference frame of Gaia, defined by the system of proper motions for sources with G ≲ 13. Available data, including the HIPPARCOS–Gaia proper motions mentioned above, are not sufficient to decide if the transition from the faint to the bright reference frame occurs abruptly at this magnitude, or more gradually over a few magnitudes; for the purpose of this paper, I will in general assume that the transition is abrupt so that all sources brighter than G = 13.0 are in the same reference frame.

Future Gaia data releases may provide proper motions that are an order of magnitude more precise than they were in DR2, and similar or even greater improvements could be obtained in the systematics. The quasars, of which many more will be found, will continue to be the main tool for examining the quality of the Gaia reference frame at faint magnitudes. Assessing its consistency at brighter magnitudes will be much harder. The HIPPARCOS–Gaia proper motions will then be of little use because their random and systematic errors are already dominated by the position errors in the HIPPARCOS celestial reference frame (HCRF), which remain unchanged. The quality of the bright reference frame can only be verified by means of the positions and proper motions of bright Gaia sources, measured to sufficient accuracy in the ICRF frame by some independent method. Most obviously, this can be done by means of differential VLBI, where the positions of radio stars are measured relative to quasars by phase-referencing techniques (Lestrade et al. 1990; Beasley & Conway 1995; Rioja & Dodson 2011; Fomalont 2012). These relative measurements are already reaching microarcsecond (μas) precision (Reid & Honma 2014).

The purpose of this paper is to explore the use of differential VLBI observations of radio stars for verification of Gaia’s bright reference frame. In Sect. 2 the required formalism is developed, whereby the positions and proper motions measured by VLBI are connected to the Gaia data. The method is tested on Gaia DR2 data in Sect. 3, using a selection of published VLBI observations, and the results and possible future improvements are discussed in Sect. 4.

2. Theory

2.1. ICRS and celestial reference frames

The ICRS is an idealised system of astronomical coordinates α, δ, whose axes are defined by convention and remain fixed with respect to distant matter in the Universe (Arias et al. 1995). The origin is at the Solar System barycentre. Any astrometric catalogue where the positions and proper motions nominally refer to the ICRS can be regarded as a practical realisation of the idealised system, and is then called a celestial reference frame (CRF). The ICRF, HCRF, and Gaia-CRF2 are examples of such reference frames, among which the ICRF has the privileged status of actually defining the conventional axes of the ICRS. In the present context it is necessary to consider that Gaia DR2 may represent (at least) two distinct reference frames, one defined by the faint quasars, and a second defined by the positions and proper motions of the bright stars in Gaia DR2.

Conceptually, a CRF can be visualised as a set of orthogonal unit vectors X, Y, Z, with origin at the Solar System barycentre, and with X pointing towards α = δ = 0, Z towards δ = +90°, and Y = Z × X. Let C = [X Y Z] be the vector triad representing the ICRF. By definition, C coincides with the axes of the ICRS and is fixed with respect to objects at cosmological distances such as the quasars. This means that the proper motions of quasars, when expressed in C, have no global component that can be interpreted as a solid-body rotation (spin) of C. Any other reference frame may have some small time-dependent offset from C described by the vector ε(t),

(1)

Thus ε(t) is the rotation of needed to align its axes with C. The sign of ε(t) is chosen for consistency with earlier publications (Lindegren & Kovalevsky 1995; Lindegren et al. 2012, 2016), where the frame offset was defined in the sense of a correction to the frame under investigation. The components of ε(t) in C or are denoted ε_X(t), ε_Y(t), ε_Z(t). The relation between the two frames is illustrated in Fig. 1.

Fig. 1. Convention for the definition of the orientation of the arbitrary frame with respect to the ICRF (C = [X Y Z]). The drawing illustrates the configuration when the orientation difference is a pure rotation about the Z axes by the positive angle εZ, i.e. ε = [0,  0,  εZ]′. The right ascension of the source at S is α in frame C and in frame .

Equation (1) is valid in the small-angle approximation, that is, ignoring terms of order ε², where ε = |ε| is the total angular offset between the two frames. This is a valid approximation in all practical cases, where ε ≲ 1 mas, or ε² <  5 × 10⁻⁹ mas (5 pas). Rigorous expressions are given in Sect. 6.1.2 of Lindegren et al. (2012).

The relation between the two frames can alternatively be expressed by means of the rotation matrix

(2)

where the prime (′) denotes matrix transposition and the scalar product of vectors (Murray 1983). Let u be the unit vector at a certain time from the Solar System barycentre towards a celestial object. The rectangular coordinates of the vector in the two frames are given by

(3)

where (α,  δ) and are the astronomical coordinates of the object in the two frames. The column matrices in Eq. (3) are related by the matrix equation

(4)

where is the rotation matrix in Eq. (2). To verify Eq. (4), one can note that is the unit tensor (Murray 1983), and hence .

In frame the proper motions are modelled as essentially constant angular velocities, which does not permit to have a non-uniform rotation with respect to distant matter. The variation of the offset vector with time t can therefore be written

(5)

where T is an arbitrary reference epoch, and ε(T) and ω are constant vectors with components ε_X(T), ω_X, etc.

Equation (5) describes a uniform solid rotation of one frame relative to the other. The word “rotation” is ambiguous in this context because it may refer to either the instantaneous configuration ε(t) or the angular velocity ω. In the following we use “spin” for the angular velocity and “orientation” for the instantaneous configuration, specifically ε(T); for brevity, the combined or general effect may be called “rotation”.

2.2. Differences in position and proper motion

From Eqs. (2)–(4), the following first-order expressions are obtained for the coordinate differences:

(6)

(7)

The time derivative of the above gives the corresponding expressions for the proper motion differences,

(8)

(9)

where μ_α* = (dα/dt)cos δ, μ_δ = dδ/dt are the components of the proper motion in frame C and , the components in . In the small-angle approximation, one can use either set of coordinates, (α, δ) or , for the trigonometric factors in (6)–(9); the choice made here is arbitrary.

The use of equations such as Eqs. (6)–(9) for estimating the difference in orientation and spin between two astrometric catalogues has been well established in the literature for many years (among many others, e.g., Fricke 1977; Froeschle & Kovalevsky 1982; Arias et al. 1988; Brosche et al. 1991; Lestrade et al. 1995; Zhu 2000; Metz & Geffert 2004; Fedorov et al. 2011; Bobylev 2015). To estimate the difference in spin (ω), the typical procedure has been to set up Eqs. (8) and (9), using the proper motion differences for a number of sources with accurate proper motions in both catalogues (or for which the true proper motions can be assumed to be negligible), and solve the resulting overdetermined system of equations using the least-squares method. It is thus possible to estimate ω without knowing ε because the differences in α and δ between the two catalogues are of second order in Eqs. (8) and (9).

The orientation difference (ε) can be estimated by applying a similar procedure to the position differences, resulting in a set of equations like Eqs. (6) and (7) for the three unknowns ε_X(T), ε_Y(T), and ε_Z(T). If the sources have proper motion, a complication is that the position differences must be computed for a fixed common epoch (T), which may require one or both sets of positions to be propagated from their mean epoch of observation. This propagation must in addition take into account any difference in spin between the two catalogues. Except when both catalogues contain only sources with zero proper motion, it is therefore usually not possible to estimate ε(T) independently of ω.

The general procedure should consequently consider the joint estimation of ε(T) and ω. This may in fact lead to a much better determination of ω than if only proper motion differences are used. The details of the procedure are worked out below, but the basic idea is simple enough: if a set of independent positional differences are obtained for a range of epochs, the resulting Eqs. (6) and (7) will depend on both ε(T) and ω, allowing all six parameters to be determined. In the current context, this means that positional VLBI observations of Gaia sources, suitably spread out in time, will contribute to the determination of the spin. This is true even when there is only a single epoch of VLBI data per source, so that their proper motions (and parallaxes) cannot be determined purely from the VLBI observations.

The realisation that positional observations contribute to the determination of the spin is of course not new. It was implicit in several of the methods used to link the HIPPARCOS catalogue to the ICRS (Kovalevsky et al. 1997), and explicitly discussed by Walter & Sovers (2000, Chap. 7.4), who concluded that it might become desirable to revise the HIPPARCOS link if, in the future, many more radio stars obtained accurate interferometric positions.

2.3. Joint estimation of the rotation parameters

The joint estimation of ε(T) and ω from the Gaia and VLBI data for a certain set S of common sources is now considered. It is assumed that the Gaia observations refer to frame and the VLBI observations to C (=ICRS), with Eqs. (1) and (5) connecting the frames. For conciseness, the six unknown rotation parameters are written as the column matrix x = [ε_X(T),ε_Y(T),ε_Z(T),ω_X, ω_Y, ω_Z]′. The result is an estimate of x, denoted , together with its 6 × 6 covariance matrix. As long as the full covariance of is retained, the choice of T is in principle arbitrary, and for convenience, we adopt the same reference epoch as for the Gaia data, for example T = J2015.5 for Gaia DR2.

The six rotation parameters in x are not the only unknowns of the problem. Both the VLBI and the Gaia observations provide information on the astrometric parameters of the sources in the common set S, and a basic assumption is that each source i in S can only have one set of “true” astrometric parameters, here denoted by the column matrix y_i. Considering m = |S| sources, the end result of the estimation process consists of and the m estimates for i = 1…m. Effectively, each is the weighted mean of the astrometric parameters as determined by VLBI and by Gaia, after correcting the latter values for the estimated frame rotation .

As discussed in Sect. 2.4, the comparison of Gaia and VLBI measurements is potentially affected by a number of difficulties including radio–optical offsets and non-linear motions. In order to proceed with the theoretical development, these difficulties are ignored here. It is assumed that the optical and radio data refer to the same physical point source, and that this source moves through space at uniform velocity relative to the Solar System barycentre. Astrometrically, then, the radio star is completely described by the usual five parameters α, δ, ϖ (parallax), μ_α*, and μ_δ, referred to the adopted epoch T, and the radial velocity v_r, assumed to be known from spectroscopy. These parameters describe the “true” motion of the source in frame C, and differ in general both from the Gaia parameters and from those derived from the VLBI measurements. The subsequent treatment is vastly simplified if expressions are linearised around a fixed set of reference values, which is an acceptable approximation as the differences are typically much smaller than an arcsecond for a suitable choice of reference values (cf. the discussion in Sect. 2.1).

It is convenient to use the astrometric parameters as given by Gaia as reference values for the linearisation. The parameter array for source i,

(10)

thus consists of corrections to be added to the Gaia parameters. With m sources, the total number of parameters to estimate is 6 + 5m. The estimation is done using a weighted least-squares algorithm, using as “observations” the Gaia data, hereafter denoted g_i, and the VLBI measurements, denoted f_i. We now proceed to detail how these observations depend on the unknowns.

The general model of the Gaia data is

(11)

where G_i is a function mapping the model parameters to the expected Gaia data, according to the model, and γ_i is the noise. In the standard five-parameter model, g_i and γ_i are 5 × 1 column matrices. It is assumed that the Gaia data are unbiased, except for the frame rotation, and that the uncertainties are correctly represented in the Gaia catalogue, so that

(12)

where C_i is the 5 × 5 covariance matrix of the Gaia parameters for source i. Expressing both y_i and g_i differentially with respect to the Gaia values, we have g_i = 0 and the linearised version of Eq. (11) becomes

(13)

where

(14)

is the matrix containing the trigonometric factors from Eqs. (6)–(9). The third row of the matrix is zero because the parallax is unaffected by the frame rotation.

Although the linearised form of Eq. (13), with dim(y_i) = 5, is used for the rest of this paper, the more general expression, Eq. (11), should be retained for future reference, when the Gaia observations may provide additional parameters (Sect. 2.4).

The description of the VLBI data for source i is similarly written in the general form

(15)

where F_i is a function mapping the source parameters to the expected VLBI data. The rotation parameters x do not enter here because the VLBI measurements are assumed to be in the ICRF frame. The VLBI measurement errors are represented by the column matrix ν_i, with

(16)

and known covariance matrix V_i. The dimension of f_i is n_i × 1, where n_i depends on the number of VLBI measurements and their state of reduction. If the measurements have been reduced to a set of five astrometric parameters, similar to the Gaia data but referring to some epoch t_i chosen specifically for these observations, we have n_i = 5. The VLBI data could also consist of a single measurement of the topocentric position at epoch t_i, however, in which case n_i = 2; or of a sequence of topocentric positions at different epochs. In either case, F_i(y_i) involves a propagation of the source parameters from the reference epoch T to the specific epoch(s) of the VLBI data t_i. For the standard five-parameter astrometric model this propagation should be done as described in Appendix A.

Recalling that y_i = 0 represents the source parameters according to Gaia, we see that Δf_i = f_i − F_i(0) contains the differences between the actually observed VLBI data f_i and the values F_i(0) computed by propagating the Gaia parameters to the VLBI epoch. To first order in y_i, the linearised version of Eq. (15) is therefore

(17)

where is the Jacobian matrix evaluated at y_i = 0. If f_i consists of the standard n_i = 5 astrometric parameters, taken in the same order as in Eq. (10) but referring to epoch t_i, then the Jacobian is approximately given by

(18)

This expression is accurate to first order in the total proper motion over the time interval t_i − T, which in some cases could amount to many arcseconds. Because it is then not obvious that Eq. (18) is a sufficiently good approximation, it is advisable to evaluate M_i by numerical differentiation of the propagation formulae.

The generalised least-squares estimate is obtained by minimising the loss function

(19)

On the assumption of Gaussian errors, the likelihood function is proportional to exp(−Q/2), and minimising Q is then equivalent to a maximum-likelihood estimation.

Setting the partial derivative of Q with respect to each model parameter equal to zero gives the symmetric system of 6 + 5m linear equations, known as the normal equations,

(20)

(21)

These can be solved by standard numerical methods, and the inverse of the normal matrix provides an estimate of the covariance of the model parameters.

Computationally, the solution of Eqs. (20) and (21) is unproblematic as it involves only a moderate number of unknowns. In terms of numerical accuracy, it is advantageous to compute the least-squares solution using orthogonal transformations (e.g. Björck 1996) after transforming the observation equations, Eqs. (13) and (17), to an equivalent set of uncorrelated unit-weight equations. Details of this procedure are not given here.

While the least-squares problem is thus solved, there is some additional insight to be gained by further manipulation of Eqs. (19)–(21). Using Eq. (21), it is possible to write each y_i in terms of x; inserting these into Eq. (20) yields a reduced system of normal equations involving only the common parameters¹,

(22)

where

(23)

(24)

and

(25)

Solving Eq. (22) yields , and then from the m Eqs. (21). Clearly, this solution is mathematically the same as obtained directly from Eqs. (20) and (21). It is more remarkable that the covariance of (the upper left 6 × 6 submatrix of the inverse of the full normal matrix) is obtained from the reduced system as (∑_i ∈ SN_i)⁻¹.

By a similar process of eliminating the unknowns y_i, the loss function Eq. (19) can be written in terms of x as

(26)

with

(27)

The interpretation of the above equations is straightforward. Δf_i − M_iK_ix is the residual of the VLBI data with respect to the values predicted from the Gaia data, after correcting for the rotation parameters and propagating to the VLBI epoch. D_i in (25) is the covariance of Δf_i, including the contributions from the uncertainties of both the VLBI and propagated Gaia data. Equation (27) shows that minimises the sum of the squares of the VLBI residuals after normalisation by the combined uncertainties. For a given solution , we may take the quantity

(28)

as a measure of the discrepancy for source i, where n_i = dim(f_i) is the number of VLBI data points included for the source. The normalisation by n_i is essential in order to avoid penalising sources with many VLBI data points. Q_i/n_i can be interpreted as the reduced chi-square of the source, and should ideally be around unity if the astrometric model fits the source and the uncertainties are correctly estimated. Q_i/n_i is hereafter referred to as the “discrepancy measure” of the source relative to a given solution.

Equations (22)–(27) have some practical advantages over the use of orthogonal transformations to solve the least-squares problem. For the identification of outliers (Sect. 2.4), computing the solution and other statistics for a very large number of different subsets of S may be required. This can be done most efficiently by pre-computing D_i, N_i, b_i, and other quantities that do not depend on the solution.

The matrices N_i are, furthermore, useful for quantifying the amount of (Fisher) information on x contributed by each source. A source without VLBI data would formally have infinite V_i and hence N_i = 0. A source with only proper motion data from VLBI will not contribute to the estimation of ε and will consequently have zeros in the first three rows and columns of N_i. More generally, the amount of information contributed by source i to the estimation of ε(T) and ω is quantified by

(29)

respectively, where the trace is limited to the first three diagonal elements of N_i for E_i, and to the last three for Ω_i.

2.4. Modelling issues and robustness

In the preceding treatment it was assumed that the sources move through space with uniform velocity relative to the Solar System barycentre, allowing both the Gaia and VLBI measurements to be accurately modelled by five astrometric parameters per source, plus a spectroscopically determined radial velocity. This model is manifestly incorrect for a number of radio stars known to be members of binaries or more complex systems, for which the VLBI observations have determined non-linear motions (e.g. the T Tau system; Loinard et al. 2007) or even complete orbits (e.g. the Algol and UX Arietis systems; Peterson et al. 2011). A second assumption is that the centre of radio emission coincides with that of the optical emission, which is also not true for many objects with extended atmospheres, discs, jets, and other structures in the radio and/or optical images. The astrometric biases produced by these various effects range over many decades, from the undetectable to tens of milliarcseconds. As a result, the simple modelling described above will provide excellent fits for some sources and large residuals for others, with a continuum of intermediate cases.

The general method of estimation in Sect. 2.3 does permit the application of more sophisticated source models. Depending on the physical nature of a radio star, it may be possible to improve the modelling, and ultimately the accuracy of , by introducing a small number of additional unknowns, thus extending the array y_i in Eq. (10). For example, in an interacting binary it may be possible to model the offsets of the optical and radio emissions from the barycentre of the system in terms of a few geometrical elements if the main characteristics of the binary, including its period, are known from spectroscopy. Another example is the non-linear motion of a radio star perturbed by a distant companion. In this case, the comparison of radio and optical observations, taken at different epochs, is meaningful only if the non-linearity is taken into account through the addition of a few acceleration terms or possibly a complete set of orbital elements. In either case, the modelling requires an augmented parameter array y_i and more elaborate expressions for the functions G_i and F_i in Eqs. (11) and (15) than discussed in Sect. 2.3. However, the details of this are beyond the scope of this paper, where the standard model of Appendix A is used throughout.

To cope with unmodelled effects, whether they are radio–optical offsets, non-linear motions, or deviations from more elaborate models, it is imperative that the estimation procedure is robust, that is, that the result is not overly sensitive to the relatively few cases with large perturbations. The least-squares estimation of Sect. 2.3 is inherently non-robust and needs to be modified or complemented with other techniques to provide the required robustness. The strategy adopted in this paper is to identify the most problematic sources and exclude them from the solution. Consistent with the treatment in Sect. 2.3, each source with all its data is regarded as an independent entity, to be either included or rejected. Although it could happen that the model fits the data very well in one coordinate (say, α), but not in the other (δ), a reasonable standpoint is that such a source is better left out entirely. Discrepant sources can be identified by means of statistics such as Eq. (28), and by comparing solutions for different subsets of S. Details of the procedure are explained in Sect. 3.3 as it is applied to actual data.

3. Application to Gaia DR2 data

In this section we use the algorithms in Sect. 2.3 to estimate the rotation parameters of the bright Gaia DR2 reference frame, based on VLBI astrometry of radio stars collected from the literature. The primary goal is to verify, if possible, the spin detected from comparison with HIPPARCOS data (Sect. 1), but an important secondary goal is to illustrate the usefulness of positional VLBI data for estimating the spin, compared with a solution using only proper motions.

3.1. VLBI data

Recent technological advances have dramatically expanded the scope for stellar radio astrophysics (Matthews 2019). The amount of accurate VLBI astrometry that could potentially be used for validating the stellar reference frame of Gaia increases rapidly, not least thanks to a number of surveys aiming to study Galactic structure (e.g. BeSSeL, Brunthaler et al. 2011; VERA, Honma 2013; GOBELINS, Ortiz-León et al. 2018). A recent comparison of Gaia DR2 and VLBI stellar parallaxes (Xu et al. 2019) lists more than 100 targets.

The 41 objects considered in this study are listed in Table 1, with their corresponding Gaia DR2 identifiers in Table 2. All the sources are brighter than G = 13.0 and have full (five-parameter) astrometry in Gaia DR2. The list includes most radio stars from the early programmes, in particular, Lestrade et al. (1999), in view of their potentially high weight in the estimation of the spin components. These were complemented by a selection of more recent data mainly from the GOBELINS survey, which have observation epochs close to the Gaia DR2 epoch and therefore could contribute usefully to the estimation of the frame orientation at J2015.5. Most of the objects are young stellar objects, interacting binaries, or giants with extended atmospheres. Their celestial distribution is shown in Fig. 2.

Fig. 2. Sky distribution of the radio sources in Tables 1 and 2. Filled and open symbols mark the objects that were accepted and rejected in the baseline solution. The map is a Hammer–Aitoff projection in ICRS coordinates with α = δ = 0 in the centre, α increasing from right to left, and the Galactic equator is plotted in red.

The VLBI data have been collected from some 20 different publications as listed in Table 1. In some cases, the authors did not publish the barycentric position at mean epoch obtained in their analysis of the observations. In most of these cases, the authors did provide results from their individual VLBI sessions, however, including the dates and geocentric positions, from which the required information could be reconstructed (cf. Sect. 4.5). For the sources with reference number 5, 7, 11, 15, and 16 in the last column of Table 1, this is how the listed barycentric positions α(t) and δ(t) were derived, with their uncertainties, while the parallaxes and proper motions were taken from the cited references.

For Mira variables and red supergiants, the VLBI observations locate several maser spots in the very extended (∼100 mas) circumstellar envelopes. When a kinematic model is employed for the relative motions of the spots caused by expansion and rotation of the envelope, it is often possible to infer the position of the geometrical centre of the star (e.g. Zhang et al. 2012; Nakagawa et al. 2014). The VLBI positions of these objects given in Table 1 are not accurate enough to contribute significantly to the determination of the rotation parameters, but the systemic proper motions may contribute to the spin. For reference number 14, the positions in Table 1 were taken from SIMBAD and are not used in the solutions except to compute the trigonometric factors in Eq. (14).

For some older observation series, in particular, Lestrade et al. (1999, reference number 3), the original results have been updated using more accurate calibrator (quasar) positions from the ICRF3 catalogue², resulting in some considerable improvements. An example is Cyg X-1 (HD 226868), where the positional uncertainties at 1991.25 were reduced from (σ_α*,  σ_δ) = (1.24,  1.75) mas, as given in Lestrade et al. (1999), to the (0.308,  0.368) mas of Table 1.

In Eq. (17) the VLBI data are compared with the Gaia data propagated to the epoch (t) of the VLBI data. This was done using the formulae in Appendix A, which requires knowledge of the radial motion of the source in order to take perspective effects into account. When available, radial velocities were taken from SIMBAD (Wenger et al. 2000); otherwise, a value of zero was used. The VLBI data were assumed to be uncorrelated, that is, all matrices V_i were taken as diagonal. For the present application, the most critical correlation is that between position and proper motion, which should normally be small if t is close to the mean epoch of the VLBI measurements.

3.2. Gaia data

Gaia DR2 identifiers for the optical counterparts of the radio sources are given in Table 2. For most sources they were taken directly from SIMBAD, but in a few cases (HD 22468, T Tau, VY CMa, and AR Lac), they had to be derived by cross-matching the radio positions with the Gaia positions. Relevant optical data were retrieved from the Gaia Archive³ and the full covariance matrices C_i computed from the formal uncertainties and correlation coefficients. As these data are readily available on-line, they are not reproduced in Table 2, with the exception of the G magnitude and the re-normalised unit weight error (RUWE). The latter, computed from Archive data as described in Lindegren (2018), is a goodness-of-fit measure (formally the square root of the reduced chi-square of the astrometric solution) that should be around 1.0 for an astrometrically well-behaved source. RUWE ≳ 1.4 could indicate an astrometric binary, a (partially) resolved binary or multiple star, or an otherwise problematic source. The remaining columns in Table 2 are explained below.

It is known that the parallaxes in Gaia DR2 are systematically too small by a few tens of microarcseconds (Arenou et al. 2018). The zero-point is estimated to be about −0.03 mas for the faint quasars (Lindegren et al. 2018), but there is strong evidence that the bright stars of interest here have a more negative zero-point, around −0.05 mas (e.g. Riess et al. 2018; Zinn et al. 2018; Schönrich et al. 2019). This is broadly confirmed by the VLBI data: the median parallax difference for the 41 radio stars in Table 1 is ϖ_DR2 − ϖ_VLBI = −0.076 ± 0.025 mas. In the solutions reported below, the Gaia DR2 parallaxes of all the radio stars have been increased by 0.05 mas to take the zero-point error into account. Because parallax and the other astrometric parameters in the Gaia data are correlated, this changes the estimated rotation parameters, but only by small amounts: about 0.010 mas in ε(2015.5) and 0.005 mas yr⁻¹ in ω.

3.3. Results

Direct application of the solution method in Sect. 2.3 to the full sample of 41 sources gives a very poor fit as measured by the loss function: Q ≃ 920 915 for n ≡ ∑_in_i = 224 degrees of freedom, or a reduced chi-square of Q/n ≃ 4111. This solution also gives unrealistically high values for the spin parameters, with mas yr⁻¹. Inspection of the discrepancy measure Q_i/n_i of the individual sources shows that T Tau has by far the highest value Q_i/n_i ≃ 159 819, followed by W 40 IRS 5 at Q_i/n_i ≃ 8681, and so on. Removing T Tau from the sample and re-computing the solution and discrepancy measures gives Q ≃ 113 037 for n = ∑_in_i = 219 degrees of freedom (Q/n ≃ 516). In this solution the source with the largest discrepancy measure is W 40 IRS 5 at Q_i/n_i ≃ 8714. Removing this source as well and iterating the procedure until all sources but one have been removed gives a series of solutions with k = 0,  1, … rejected sources. The evolution of max(Q_i/n_i) and Q/n as functions of k is shown in the left panels of Fig. 3, with the corresponding orientation and spin parameters in the right panels.

Fig. 3. Evolution of discrepancy measures and rotation parameters in a series of solutions where the k most discrepant sources have been removed. a: discrepancy measure Qi/ni for the most discrepant source in the solution. Points are labelled with the name of the source. b: total discrepancy measure Q/n as a function of the number of rejected sources. c: estimated orientation parameters. d: spin parameters as functions of the number of rejected sources. Error bars are formal ±1σ uncertainties from the inverse normal matrix, i.e. not adjusted depending on the goodness-of-fit.

Errors in the rotation parameters caused by non-linear source motions and other model deficiencies are generally reduced as more outliers are removed, while the statistical (formal) uncertainties increase because fewer sources contribute to the solution. The optimum solution is a compromise between the opposite tendencies and may be found somewhere along the sequence of solutions described above. The rather smooth progression of the discrepancy measures in Fig. 3a and b gives no clear indication of the optimum k, except that it is probably in the range from 7 (removing the most obvious outliers) to ≃30 (after which Q/n <  1). The spin parameters in Fig. 3d show an abrupt change with the removal of HD 282630 at k = 11, after which the spin solution is relatively stable in the Y and Z components, but less so in X. The orientation parameters in Fig. 3c show a similar behaviour. Any k in the range from 13 to 35 in fact gives fairly consistent solutions at least in Y and Z. For the subsequent discussion, we take somewhat arbitrarily the solution at k = 15 as the baseline. This solution has Q/n = 3.55 with n = 139 degrees of freedom, and the most discrepant source is HD 290863 with Q_i/n_i ≃ 16.074.

The accepted and rejected sources and their individual discrepancy measures for the baseline solution are listed in Table 2. The rotation parameters are given in Table 3 under entry A, including the formal uncertainties calculated from the inverse normal matrix. The correlation matrix (for T = J2015.5) is

(30)

Given the reduced chi-square of the solution, Q/n ≃ 3.55, it is likely that the formal uncertainties underestimate the actual errors. More realistic estimates may be obtained by bootstrap resampling of the 26 accepted sources, yielding the uncertainties in the “Adopted” entry of Table 3.

From the original sample of m = 41 sources, k = 15 were thus removed to obtain the baseline subset of m − k = 26 sources. It is not obvious that the process of successively removing the most discrepant source leads to the optimum subset in the sense that no other subset of the same size has a smaller Q/n. It is conceivable that a different procedure, for instance, starting from a smaller subset and adding the best-fitting sources (outward selection; Ben-Gal 2010), would lead to a different result. While it is impractical to test all possible combinations, an exhaustive search of the different subsets of size 26 drawn from the 33 sources with the smallest discrepancy measure did not uncover any more favourable subset. Figure 4 displays Q/n versus the spin components for these solutions. Many of them are quite different from the baseline solution in terms of the spin components, but invariably, their Q/n is then also significantly higher. This makes it credible that the spin parameters of the adopted solution are statistically significant and not the chance result of a particular combination of data for a few sources.

Fig. 4. Sensitivity of the solution to the selection of sources. The plots show the spin parameters and associated Q/n values for about four million solutions using different subsets of the sources as described in the text. The colour shows the smallest Q/n in each bin. The dashed lines mark the spin parameters for the baseline solution A with Q/n ≃ 3.55.

The solution gives improved estimates of the astrometric parameters of the sources, obtained by solving (21) for each i. These results are not tabulated, as they are practically the same as obtained from a weighted average of the VLBI and Gaia data after correcting the latter for the frame rotation and parallax zero-point. For example, the joint estimate of the parallax of Cyg X-1 is mas, very close to the weighted mean of the VLBI value, 0.547 ± 0.041 mas, and the Gaia DR2 value after correction for the zero-point, 0.472 ± 0.032 mas.

3.4. Alternative solutions

The classical way to determine ω is to solve the overdetermined system of Eqs. (8) and (9) using only the proper motion differences. In this process it is natural to assign a weight to each equation that is inversely proportional to the combined variances of the VLBI and Gaia proper motions. An equivalently weighted least-squares solution is obtained with the formalism of Sect. 2.3 simply by deleting in f_i all the VLBI data items that are not proper motions. The resulting normal equations, Eq. (22), are of course singular for the orientation parameters, but the lower-right 3 × 3 part of the equations gives the desired solution for ω. Applying this procedure to the baseline subset of 26 sources gives the result shown as solution B in Table 3. The spin parameters are reasonably close to those of the baseline solution (A), testifying to the good internal consistency of the (accepted) data. It is more interesting, however, that the formal uncertainties are significantly higher in solution B than in A: the uncertainty is a factor two higher for the Y component.

It thus appears that the positional VLBI data are at least as valuable as the proper motion data when estimating the spin, at least for the spread of VLBI epochs considered here. A direct test of this hypothesis is to make the complementary solution to what is described above, that is, deleting all proper motion items in f_i. The result is shown as solution C in Table 3. In this case, both the orientation and spin parameters are determined, and the spin is actually more precise (in terms of formal uncertainties) than in solution B where the proper motions (only) are used, although not as precise as in A. It can be noted that only 23 of the 26 sources contribute to solution C because no positional VLBI data are provided for S CrB, U Her, and RR Aql, although these sources belong to the baseline subset.

Figure 4 of Lindegren et al. (2018) suggests that the transition from the faint to the bright reference frame in Gaia DR2 does not occur abruptly at G = 13 but gradually from G ≃ 13 to ≃11 mag. Several of the sources in Table 1 have magnitudes in the transition interval and may therefore not contribute fully to the determination of the rotation parameters. In solution D the model of the Gaia data in (11) is modified so that the applied rotation is x multiplied by the magnitude-dependent function

(31)

Unfortunately, this does not improve the overall fit significantly (Q = 492.80 in D against 493.14 in A), and the discrepancy measure actually increases for four of the six accepted sources in the magnitude range 11–13. Thus the VLBI data do not clearly support the magnitude-dependent model, although it cannot be ruled out either.

3.5. Weights of the individual sources

Table 2 includes the statistics E_i and Ω_i from Eq. (29), indicating the potential weights of the sources in the orientation and spin solutions. The sources that contribute most weight to the spin solution (largest Ω_i) are AR Lac, LS I +61 303, Cyg X-1, HD 199178, and V410 Tau; all of them have Ω_i >  300 mas⁻² yr² and all are included in the baseline subset. The first three and HD 199178 mainly contribute by virtue of their relatively small positional uncertainties (≃0.3 mas) at a very early epoch, t ≃ 1992. IM Peg, one of the best-observed radio stars (Bartel et al. 2015), has a much smaller weight in this analysis because its uncertainties in Gaia DR2 are relatively large.

The main contributions to the determination of the orientation (largest E_i) come from young stellar objects in the Taurus and Orion regions that are observed as part of the GOBELINS survey (Galli et al. 2018; Kounkel et al. 2017), through their high precision and proximity in time with the Gaia DR2 epoch. The strong concentration of these sources in a small part of the sky, approximately in the +Y direction, is responsible for the relatively large uncertainty of the Y component of ε(2015.5) in Table 3.

4. Discussion

4.1. Estimated rotation parameters

The result in Table 3 yields a spin for the Gaia DR2 bright reference frame of ω  =  (−0.077, −0.096, −0.002) mas yr⁻¹, significant at the 2σ level in the Y component. Figure 3d indeed suggests that ω_Y is quite firmly fixed around −0.1 mas yr⁻¹, while the results in Z and (especially) X are less certain.

The orientation parameters ε(2015.5) are less well determined, and especially the Y component has a large uncertainty that is due to the unfavourable celestial distribution of the more recent VLBI observations included in the analysis (Sect. 3.5). Still, an orientation error of about 0.5 mas is indicated at the 2σ level in the X and Y components.

As mentioned in the Introduction, the non-zero spin of the bright reference frame of Gaia DR2 is seen also in a comparison with proper motions of HIPPARCOS stars calculated from the position differences between Gaia DR2 and the HIPPARCOS catalogue, divided by the ∼24 yr epoch difference. This comparison was made by Brandt (2018) in the course of constructing The HIPPARCOS–Gaia Catalog of Accelerations (HGCA). For 115 663 HIPPARCOS stars the HGCA gives three essentially independent sets of proper motions, namely (i) as measured by HIPPARCOS around epoch 1991.25, (ii) as measured by Gaia (in DR2) around epoch 2015.5, and (iii) as calculated from the HIPPARCOS–Gaia position differences divided by the epoch difference. The HGCA is intended for orbit fitting and for identifying candidate stars with substellar or dark companions, and to this end, Brandt (2018) made a careful cross-calibration of the three sets of proper motions in order to eliminate any systematic offsets, in particular the global rotations. For the rotation between sets (ii) and (iii) the result as given in his Table 1 is ω = ( − 0.081, −0.113, −0.038) mas yr⁻¹. Brandt used the same sign convention for this vector as in this paper, so that his result can be directly compared with our Table 3. The agreement is very good for all the solutions in Table 3; for the adopted solution, none of the component differences exceeds one standard deviation.

The comparison with Brandt (2018) rests on the assumption that the proper motions calculated from the HIPPARCOS–Gaia position differences are absolute, that is, that they are expressed in a reference frame that is non-rotating with respect to the ICRS⁴. The validity of this assumption depends on the quality of the positional reference frame of HIPPARCOS at epoch 1991.25 and of the positional reference frame for the same stars in Gaia DR2 at epoch 2015.5. If both sets of positions are aligned with the ICRS at their respective epochs, the proper motions calculated from the position differences must clearly be non-rotating with respect to the ICRS. In The HIPPARCOS and Tycho Catalogues (ESA 1997, Vol. 3, Sect. 18.7) the rms orientation error of the HIPPARCOS reference frame at epoch 1991.25 is given as 0.6 mas per axis. For the bright DR2 frame at epoch 2015.5 we may conclude from Cols. 2–4 in Table 3 that the rms error is probably smaller than 0.4 mas per axis. This gives a total rms uncertainty of mas yr⁻¹ in each component of ω.

The old HIPPARCOS positions thus provide an independent estimate of the spin, in good agreement with the results from the VLBI observations, and of comparable or even better accuracy. However, because its uncertainty is dominated by the unknown orientation errors of the HIPPARCOS reference frame, which will not improve, the method will be of limited value for the validation of future Gaia data releases. In contrast, the VLBI method, as discussed in Sect. 4.4, has great potential for future improvement.

The origin of the non-zero rotation parameters for the bright reference frame of Gaia DR2 is a deficiency in the specific instrument calibration model used for the DR2 astrometric solution. The relevance of this deficiency for the reference frame was not recognised at the time when the DR2 data were prepared and validated, and although its effect on the bright sources was noted (e.g. Fig. 4 in Lindegren et al. 2018), no explanation of its origin was offered. This is now understood, and the calibration model is being improved with a view to avoid a similar error in future data releases (see Appendix B).

4.2. Use of positional VLBI data

A comparison of the formal uncertainties of ω in solutions A, B, and C (Table 3) demonstrates the advantage of including positional VLBI data in a joint solution for the orientation and spin parameters. Not only do the positional data, as already noted, improve the spin solution (A is better than B for ω), but the inclusion of proper motion data also improves the orientation parameters (A is better than C for ε).

The value of positional VLBI data for the spin is especially noteworthy because it means that even single-epoch VLBI astrometry can be incorporated in the solution. This will contribute to the spin determination, especially if the data are taken at an epoch that is well separated from the Gaia epoch. This is important to keep in mind, both for the inclusion of old, possibly unpublished VLBI measurements (Sect. 4.5) and for the scheduling of future VLBI sessions specifically for the reference frame (Sect. 4.4). The latter need not be constrained by considerations of parallax factor and temporal spread of the measurements for the proper motions.

4.3. Binarity and source structure

About half of the radio sources in Table 1 are known to be binaries or members of double or multiple systems. Interacting pairs with periods from a day to tens of days include the RS CVn systems and high-mass X-ray binaries, which usually provide good fits to the single-star model unless they are perturbed by a more distant component, as is the case for σ² CrB. Binaries with periods of years to hundreds of years are more problematic unless a complete orbit can be determined. Orbits have been determined for some radio stars, but as the corresponding binary data from Gaia are not yet available, they are not included in this analysis. For other objects, the VLBI observations have detected the curvature of a long-period orbit by means of acceleration terms; this is the case for example for UX Ari, HD 283447, and T Tau, all of which obtain large discrepancy measures in the present analysis (Fig. 3a).

The radio stars included in the present analysis have not been a priori screened for known or anticipated problems with multiplicity and source structure. Apart from known binaries, several Mira variables and red supergiants have therefore been included, although they are far from ideal targets because of the issues described in Sect. 3.1. It is very likely that the optical Gaia observations are also severely affected by the extended and complex atmospheres (Chiavassa et al. 2011). This type of radio star should probably be avoided entirely for the reference frame.

A statistic that could be used for screening the sources is the RUWE given in Table 2. The RUWE measures how well the different Gaia observations that were made over a few years agree with the five-parameter single-star model. It is therefore mainly sensitive to the presence of companions at separations from a few milliarcsecond to about one arcsecond, and a relevant indicator of potentially problematic sources of that particular kind. It is not sensitive to non-linear motions, if the deviation from linear is small over the few years of the Gaia observations. A large discrepancy measure Q_i/n_i also indicates a problematic source, but of a rather different kind. In contrast to the RUWE, it is sensitive to radio–optical offsets, and it is also more sensitive to long-period perturbations if the VLBI and Gaia measurements are made at very different epochs. Figure 5 shows a weak positive correlation between the two statistics, which is to be expected because their different regimes of sensitivity overlap. However, it is clear that RUWE alone is not sufficient to find the best candidate targets for determination of the frame rotation parameters.

Fig. 5. Discrepancy measures of the sources in solution A against the RUWE in the Gaia DR2 astrometry.

Surprisingly, one source (DoAr 51) has a very large RUWE = 6.43, but nevertheless gives a good fit in the present solution (Q_i/n_i = 1.261). This object is a triple system consisting of an equal-mass pair with a period of about 8 years and a separation of approximately 56 mas at 2015.5, and a fainter tertiary component at a separation of about 790 mas (Schaefer et al. 2018). This triple configuration could well explain the high RUWE obtained with Gaia. The VLBI observations detect both components of the close pair, and the data in Table 1 refer to their centre of mass (Ortiz-León et al. 2017a). The near-coincidence of the VLBI and Gaia epochs, together with the near-coincidence of the optical photocentre with the centre of mass of the close pair and the rather large uncertainty of the Gaia proper motion (∼1 mas yr⁻¹), could explain why the discrepancy measure is not higher. With Ω_i = 3.2 mas⁻² yr², this object contributes less to the determination of ω than any of the other accepted sources in Table 2. DoAr 51 is a good example of an object for which an extended model along the lines in Sect. 2.4 could drastically increase the usefulness of the data.

4.4. Precision of future solutions

The bright reference frame in future releases of Gaia data should ideally be validated at a level compatible with the expected errors of the proper motions, which may be as small as a few μas yr⁻¹. This will require a much better accuracy in the spin parameters than is achieved in the present analysis. Clearly, the uncertainties can be reduced by including VLBI data for more radio stars, and/or using improved data for the sources already considered. As discussed in Sect. 4.2, the epoch of the added data is an important factor. If dedicated VLBI measurements are contemplated for this purpose, different scenarios can be envisaged concerning the number, distribution, and epochs of the planned observations, and it is of great interest to predict the accuracy that can be achieved in various cases. This can be done by applying the algorithm in Sect. 2.3 to simulated data sets, using the left-hand side of Eq. (22) to compute the formal precisions.

A major uncertainty in any such prediction is the extent to which binarity, source structure, and radio-optical offsets will limit the achievable accuracy. One extreme scenario is that these effects already dominate the error budget in the current analysis. Support for this could be drawn from the fact that Q_i/n_i >  1 for nearly all the accepted sources in Table 2. In this scenario it will not help much to add more and better VLBI data for the radio stars already considered, and the safest way to improve the solution may be to increase the number of sources, m, and rely on the statistical improvement by m^−1/2.

Realistically, however, the prospects are not as bleak as outlined above. Better screening of the sample, modelling of orbital motions and offsets, etc., will surely improve the results, and this process will be helped by the addition of new data both from VLBI and Gaia. The scenario at the other extreme is that such improvements will allow us to reach the formal uncertainties computed from the least-squares equations. This (optimistic) assumption is the basis for the predictions in Table 4.

The first entry in Table 4 represents the baseline solution of the present analysis (A in Table 3). For brevity, the formal uncertainties of the six rotation parameters are condensed here into two numbers, one for the orientation and one for the spin. For the next three entries, it is assumed that new VLBI observations of the same 26 accepted sources are obtained at the specified epoch t, with a precision of 0.1 mas in each coordinate. Performing the analysis as in Sect. 2.3 with the same Gaia DR2 data, we obtain the formal uncertainties in the last two columns. Remarkably, while the orientation (at T = 2015.5) is better determined with the added data, the improvement is small in the spin and practically independent of the epoch of the added VLBI measurements. This shows that the formal uncertainties of the spin in solution A are limited by the accuracy of the proper motions in Gaia DR2 rather than by the VLBI data. The slight increase in σ[ε(T)] with t arises because the effective mean epoch of the (old plus new) VLBI measurements moves away from the reference epoch T.

The middle four entries in Table 4 show predictions when Gaia data are used based on the nominal mission length L = 5 yr. It is assumed that the uncertainties of the Gaia astrometry improve as L^−1/2 for the positions and parallaxes, and as L^−3/2 for the proper motions (taking L = 1.8 yr for DR2). The covariance matrices C_i in (11) are simply scaled by the corresponding factors, leaving the correlations unchanged from DR2. The orientation parameters now refer to the corresponding reference epoch of the Gaia observations, T = 2017.1⁵. Consistent with the previous finding that the DR2 precision is the main limitation in solution A, the improved Gaia data drastically reduce the formal uncertainty of the spin parameters even without any additional VLBI data. Adding new data for the 26 stars improves the determination of the spin still further, especially if the new measurements are at a late epoch. The orientation parameters, on the other hand, are not much improved. Their uncertainties are basically limited by the VLBI position errors, typically of the order of 0.1 mas, and the small number of sources: mas.

The last four entries in Table 4 are for an extended Gaia mission covering a full decade of observations (L = 10 yr). Here the spin parameters receive another boost in precision, while there is almost no improvement in the orientation parameters. The σ[ε(T)] has a shallow minimum at t ≃ 2024 because the effective mean epoch of the VLBI measurements is then close to the assumed reference epoch, T = 2019.6.

Taking into account that the two extreme scenarios considered above are probably to some extent true, we conclude from this crude assessment that a combination of actions should be taken to ensure that a robust and accurate estimate of the rotation parameters can be derived at the time when the final Gaia results become available. These actions should include compiling and recalibrating past VLBI measurements, as well as securing new observations of both old targets and as many new as possible.

4.5. A plea to VLBI observers

Most of the observational VLBI programmes used in this study address problems in Galactic or stellar astrophysics, and are primarily concerned with the parallaxes and proper motions of the radio stars, or of their orbits, surface structures, and similar. Consequently, many of the publications do not provide the (barycentric) positions that the authors must have derived along with the parallaxes and proper motions by fitting an astrometric model to their positional measurements. This is not a problem as long as the authors provide the individual measurements, and their times, on which the fit was based. For this reason, as mentioned in Sect. 3.1, several of the barycentric positions in Table 1 were derived by the present author by fitting the standard model (Appendix A) to the published VLBI measurements. Although this procedure is not without advantages (see below), it is of course simpler if the full solution is provided by the observers.

In this context, it is useful to stress the unique historical value of positional data. In contrast to the parallaxes, for example, which can be re-determined at a future time, positional measurements can never be repeated, and for a number of applications their value only increases with time. A plea to observers making high-precision VLBI astrometry is therefore that they publish the full result of their astrometric fits, including the barycentric position and corresponding epoch.

It is nevertheless good practice to publish the individual position measurements for future uses as well. This allows alternative models, which may include acceleration terms or orbital parameters for example, to be fitted in combination with other data. The general method described in Sect. 2.3 is readily adapted to the use of individual VLBI measurements, and this has the advantage that otherwise neglected correlations are fully taken into account.

Relative astrometry using phase-referencing techniques are converted into absolute coordinates using an assumed position in the ICRS of the reference source (calibrator), usually a quasar. To first order, small errors in the calibrator position directly transfer to the measured coordinates of the target (Reid & Honma 2014). This gives a constant offset of no consequence when the parallax and proper motion of the target are fitted to the data, but it might be important for the present application. It is customary to specify the calibrators used in phase-referencing observations, and sometimes also their assumed positions. Knowing the identities and adopted positions of the calibrators is indeed highly desirable because it allows the target positions to be corrected when improved calibrator positions become available (cf. Sect. 3.1).

5. Conclusions

This paper provides a rigorous mathematical framework for estimating the orientation and spin of the Gaia reference frame, in which the Gaia data are optimally combined with VLBI measurements of bright radio sources. The simultaneous estimation of the orientation (ε) and spin (ω) is essential for achieving the best accuracy. The method takes full advantage of past and future single-epoch VLBI measurements of Gaia sources for the determination of the spin. Independent estimates of their proper motions from VLBI can be incorporated in the solution, but are not required by the method.

Applied to published VLBI data for a sample of 41 bright (G ≤ 13 mag) radio sources, the method indicates that the bright reference frame of Gaia DR2 rotates by about 0.1 mas yr⁻¹ relative to the faint quasar frame. The solution retains 26 of the investigated radio sources, while 15 are rejected based on a statistical discrepancy measure sensitive to binarity and source structure. The adopted non-zero result in Table 3 is significant at the 2σ level. The orientation error of ≃0.5 mas of the bright reference frame with respect to ICRS also appears to be significant.

This determination of the spin is consistent with independent estimates obtained from a comparison of the extrapolated Gaia DR2 positions with the HIPPARCOS catalogue at epoch J1991.25 (Lindegren et al. 2018; Brandt 2018).

The accuracy of the present study is limited by the relatively small number of radio stars included, by the uncertainties of the Gaia DR2 proper motions, and by issues related to the astrophysical nature of the sources.

The origin of the spin of the bright reference frame of Gaia DR2 is understood, and measures have been taken to to avoid this problem in future data releases (see Appendix B). Nevertheless, it is important that the consistency of future reference frames can be validated across the full range of magnitudes, and the present method offers such a possibility for the bright stars. As many as possible of the already existing VLBI measurements of radio stars should be used for this purpose, but it is very desirable to complement this by re-observing many of them in the coming years, and if possible, add new targets to the list for improved robustness. The use and recalibration of old, possibly unpublished data should be pursued. In this context, the unique historical value of positional VLBI measurements needs to be stressed. As argued in Sect. 4.5, observers should ensure that relevant intermediate data and meta-information are preserved for optimal future uses of their data.

Acknowledgments

I wish to thank S. Klioner for many stimulating discussions on this and related topics, and for a critical reading of the first draft. I also thank the anonymous referee and A. Bombrun, U. Bastian, M. Biermann, J. Castañeda, M. Davidson, C. Fabricius, E. Gerlach, J. Hernández, D. Hobbs, R. Le Poole, P. McMillan, M. Ramos-Lerate, and N. Rowell for valuable comments and advice. Support from the Swedish National Space Board is gratefully acknowledged. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. The work has also made use of the SIMBAD database, operated at CDS, Strasbourg, France (Wenger et al. 2000). Calculations were made using MATLAB by The MathWorks, Inc. Diagrams were produced using the astronomy-oriented data handling and visualisation software TOPCAT (Taylor 2005).

References

Appendix A: Standard model of stellar motion

The propagation of stellar astrometric data in time is normally based on a set of approximations, assumed to be valid for the observed motion of a single star, or the centre of mass of a binary or multiple system, over a limited interval of decades to centuries. Referred to in The HIPPARCOS and Tycho Catalogues (ESA 1997, Vol. 1, Sect. 1.2.8) as “The Standard Model of Stellar Motion”, these approximations have been adopted as the basis for astrometric reductions at least since the days of Schlesinger (1917) and are still used for the analysis of Gaia data (Lindegren et al. 2012). Among the most important explicit or implicit assumptions of the model are (i) that the source moves with constant velocity vector relative to the Solar System barycentre, (ii) that light-time effects beyond the Solar System can be ignored, and (iii) that aberration effects caused by the curvature of the Galactic orbits can be ignored. This is not the place to discuss the validity of these assumptions⁶, but because the model, formally represented by the function F_i(y_i) in (15), is central to the paper, it may be useful to summarise the main steps of the calculation. The Barycentric Celestial Reference System (BCRS) is used throughout, with times expressed in seconds or Julian years of barycentric coordinate time (TCB), and distances in km or au.

A.1. Propagation of the astrometric parameters

We first consider the propagation of the astrometric parameters from epoch T to t. Let α, δ, ϖ, μ_α*, μ_δ, and v_r be the astrometric parameters and radial velocity at the original epoch T, and α(t), etc. their values at epoch t. The first step is to compute the barycentric unit vector towards the source at time T,

(A.1)

and the unit vectors in the directions of increasing α and δ,

(A.2)

Next we compute

(A.3)

where

(A.4)

is the astronomical unit. The vectors r and m are such that for constant space velocity, the barycentric vector to the source at time t is proportional to

(A.5)

The astrometric parameters at epoch t are thus recovered from s(t) and m after rescaling and performing the inverse operations of (A.1)–(A.3):

(A.6)

(A.7)

(A.8)

(A.9)

(A.10)

(A.11)

(A.12)

where

(A.13)

are the updated vectors along +α and +δ.

A.2. Coordinate direction to the source

We now turn to the calculation of the position of the source as obtained from a single VLBI measurement at time t. As opposed to the barycentric direction r(t) discussed above, we now require the topocentric direction from the observer towards the source. In VLBI astrometry, geometric delays are calculated from the source positions and station coordinates expressed in the BCRS frame and corrected for the gravitational delay caused by bodies in the Solar System (for details, see Sovers et al. 1998). The astrometric position measured by VLBI therefore corresponds to the geometric direction from the observer towards the target at the time of observation, unaffected by stellar aberration and gravitational deflection. This direction, also known as the coordinate direction (Murray 1983), is here denoted .

The modelling of in terms of the astrometric parameters is simple because it only involves a shift of origin from the Solar System barycentre to the observer. With b(t) denoting the position of the observer at the time of observation, expressed as BCRS coordinates in au, we have

(A.14)

where angular brackets signify vector normalisation, ⟨a⟩=a/|a|. One small complication in Eq. (A.14) is that the position of the source should be evaluated for the barycentric time t_B obtained by adding the Römer delay to the time of observation,

(A.15)

where c is the speed of light. For an observer on the Earth, the Römer delay is at most about 500 s. Neglecting the delay produces an error equal to the proper motion of the source over this time interval, which could amount to 0.17 mas for Barnard’s star. While the effect is thus negligible for most stars, it is always safer to take it into account. On the other hand, it is an acceptable approximation to use r instead of r(t) in (A.15). (For the position of the observer, b(t), we can normally take the centre of the Earth, which is readily available from standard ephemerides, although the diurnal parallax of the nearest star is as much as 32 μas). The celestial coordinates of the source are finally obtained from the X, Y, Z components of in analogy with Eqs. (A.7) and (A.8).

Appendix B: Origin of the rotation of the bright reference frame of Gaia DR2

As mentioned in Sect. 4.1, the deviating reference frame for bright (G ≲ 13) sources in Gaia DR2 is caused by a deficiency in the astrometric instrument calibration model used for DR2. Although a detailed discussion of the effect is beyond the scope of this paper, it may be useful to outline the basic mechanism and how the effect can be avoided in future data releases. Further details may be found in the documentation of forthcoming releases.

The astrometric calibration model used for Gaia DR2 is described in Sect. 3.3 of Lindegren et al. (2018). It contains several terms that depend on the “window class” (WC) of the individual CCD observation. Window classes WC0/1/2 are different schemes for sampling the pixels around a detected source (for details, see Sect. 3.3.5 in Gaia Collaboration 2016). The WC is decided by an on-board algorithm, based mainly on the G magnitude of the source estimated as it enters the field of view. WC0 (G ≲ 13) provides a small two-dimensional image of the source, from which both the along- and across-scan coordinates can be derived. For WC1 (13 ≲ G ≲ 16) and WC2 (G ≳ 16), the image is marginalised in the across-scan direction, retaining a one-dimensional array, of different lengths in WC1 and WC2, from which only the along-scan coordinate can be derived.

In the cyclic data-processing scheme adopted for Gaia (see Sect. 7.2 of Gaia Collaboration 2016), the astrometric global iterative solution (AGIS) is closely connected to a preceding step known as the intermediate data update (IDU). The IDU provides image centroids to AGIS, which are one- or two-dimensional depending on the WC; in turn, AGIS furnishes the IDU with the improved astrometry needed to calibrate the line- and point-spread functions used by the IDU for centroiding. AGIS and IDU are thus inseparable parts of a larger astrometric task, and the calibrations discussed below should be understood as including the relevant parts of the IDU calibrations, also done per WC.

The use of different calibration models for the three window classes effectively means that there are three separate instruments to calibrate, one for each WC⁷. Because the WC essentially depends on the magnitude of the source, most sources are always observed in the same WC. Disregarding, for the moment, the sources that are observed in more than one WC, we have three disjoint subsets of sources, each subset being used to calibrate one WC. If the calibrations are allowed to be arbitrary functions of time, each subset could have its own reference frame, with an arbitrary offset in orientation and spin, without causing any inconsistency in the observations. All that is needed is that the time-dependent calibration of each WC at any time exactly matches the positional offsets of the sources in the corresponding subset, caused by the rotation of its frame. Although there is no particular reason why the calibration should vary in exactly this way, it is inevitable that unrelated variations caused by model imperfections contain some such component.

In reality, there is a good deal of overlap between the window classes, such that many sources around G ≃ 13 are sometimes observed in WC0, and other times in WC1, and similarly around G ≃ 16 for WC1 and 2. This should in principle suffice to guarantee a consistent reference frame across all magnitudes because it is not possible to obtain consistent solutions for the positions and proper motions of the multi-WC sources unless these parameters are expressed in the same reference frame for all the observations. Because there is no evidence of a discontinuity in spin at G ≃ 16 (Fig. 4 in Lindegren et al. 2018), the overlap indeed seems to have provided a good connection between the reference frames of WC1 and 2. As shown in the same figure, however, it did not work for the WC0/1 transition, or perhaps a more gradual change was created from G ≃ 13 to ≃11 (Sect. 3.4).

Of relevance here is that the modelling of observations in WC0 is, for a number of reasons, much more difficult than in WC1 and 2. The two-dimensional images in WC0 require the point-spread function to be modelled in both the along- and across-scan directions, in contrast to the simpler line-spread functions used in WC1 and 2. The across-scan profile has additional dependences on the precession-induced drift rate, resulting in a more complex calibration. Furthermore, to increase the dynamic range of Gaia’s CCDs, stars brighter than G ≃ 12 normally obtain reduced integration time through the activation of TDI blocking gates (Crowley et al. 2016), and the gates also require separate calibrations. In spite of the gates, a substantial fraction of the WC0 observations are affected by pixel saturation, which causes additional issues. In Gaia DR2 the calibration of WC0 was in fact far from satisfactory, as is readily seen in a plot of the rms post-fit residual versus magnitude (Fig. 9 in Lindegren et al. 2018). This prevented a consistent astrometric modelling of the multi-WC sources across the WC0/1 boundary, which may have been the direct cause of the observed rotation.

To avoid a similar problem in future Gaia data releases, it is necessary to improve the modelling of WC0 observations. This is already part of ongoing activities towards the next releases. It will also help to maximise the number of multi-WC sources in the primary astrometric solution. Additionally, the calibration models for the different window classes may be constrained not to contain any time-dependent components representing a rotation difference between the corresponding reference frames. Although these steps should ensure a consistent reference frame for all window classes, it remains important that the consistency can be confirmed, for example by means of VLBI observations.

All Tables

All Figures

	Fig. 1. Convention for the definition of the orientation of the arbitrary frame with respect to the ICRF (C = [X Y Z]). The drawing illustrates the configuration when the orientation difference is a pure rotation about the Z axes by the positive angle εZ, i.e. ε = [0,  0,  εZ]′. The right ascension of the source at S is α in frame C and in frame .
In the text

	Fig. 2. Sky distribution of the radio sources in Tables 1 and 2. Filled and open symbols mark the objects that were accepted and rejected in the baseline solution. The map is a Hammer–Aitoff projection in ICRS coordinates with α = δ = 0 in the centre, α increasing from right to left, and the Galactic equator is plotted in red.
In the text

Fig. 3. Evolution of discrepancy measures and rotation parameters in a series of solutions where the k most discrepant sources have been removed. a: discrepancy measure Qi/ni for the most discrepant source in the solution. Points are labelled with the name of the source. b: total discrepancy measure Q/n as a function of the number of rejected sources. c: estimated orientation parameters. d: spin parameters as functions of the number of rejected sources. Error bars are formal ±1σ uncertainties from the inverse normal matrix, i.e. not adjusted depending on the goodness-of-fit.

In the text

	Fig. 4. Sensitivity of the solution to the selection of sources. The plots show the spin parameters and associated Q/n values for about four million solutions using different subsets of the sources as described in the text. The colour shows the smallest Q/n in each bin. The dashed lines mark the spin parameters for the baseline solution A with Q/n ≃ 3.55.
In the text

	Fig. 5. Discrepancy measures of the sources in solution A against the RUWE in the Gaia DR2 astrometry.
In the text