© ESO 2021

1 Introduction

The early third data release (EDR3) of the Gaia mission (Gaia Collaboration 2021) provides the positions, proper motions, and parallaxes of more than 1 billion sources. This astrometry is relative to the International Celestial Reference Frame (ICRF; Charlot et al. 2020), which is defined by distant radio quasars. Gaia EDR3 spans more than 15 magnitudes in source brightness (more than a factor of 10⁶ in flux). The Gaia satellite uses different readout modes, called window classes, and even different exposure times, which are implemented through gatings, to handle this huge dynamic range. Different classes of observations must be calibrated to one another in order to align them with the ICRF. This is particularly true for bright stars (G ≲ 13), which have a different observation mode than faint stars and for which there are no quasars of similar brightness.

In the second Gaia data release (DR2), the bright star reference frame is rotated at ~0.15 mas yr⁻¹ relative to the faint stars and quasars (Brandt 2018; Lindegren et al. 2018)¹. In EDR3, this rotation was tentatively removed by anchoring the window class 0 (WC0) reference frame (corresponding to G ≲ 13) to the positional difference of stars between the HIPPARCOS and Gaia missions (see Sect. 4.5 of Lindegren et al. 2021). The Gaia bright star reference frame thus carries with it the 0.6 mas uncertainty in the HIPPARCOS realization of the ICRF (Kovalevsky et al. 1997) divided by the ≈25-year interval between HIPPARCOS and Gaia. Throughout the magnitude range G < 13, Gaia uses a wide range of gatings; each transition carries with it the potential for systematic offsets in the reference frame. HIPPARCOS provides almost no sources at magnitudes 11–13 to serve as references.

A proper motion bias affecting stars differently depending on their apparent magnitude can bias the eccentricity distributions inferred from a sample of binaries or can incorrectly suggest that massive stars in young associations exhibit a different kinematic behavior from that of low-mass stars. Despite representing only ~4% of the EDR3 catalog, the sources brighter than G = 13 include most O- and B-type stars within 4 kpc and represent over 70% of the subset with currently available radial velocities. A good understanding of the systematics is therefore necessary in order to fully exploit the exquisite proper motion precision provided by the Gaia astrometry.

Here, we use two independent sets of sources (resolved binaries and open clusters) to show that significant systematics are still present. We propose an additional correction to align the reference frame of the proper motions of the bright Gaia EDR3 sources with that of the fainter sources and, by extension, the ICRF. Section 2 presents the data used in this study. Section 3 investigates the dependence of the spin parameters on the magnitude and color of the sources. We conclude with Sect. 4.

2 Data used

We calibrated the proper motions of the bright stars in Gaia EDR3 to those of the faint stars without referring to external data. Our calibration sample consisted of binaries with a bright (G <14) and faint (G > 14) component and of bound star clusters whose mean proper motions we estimated using their faint members. Gaia EDR3 calibrates the faint stars to quasars of similar magnitude. This section describes our samples of binaries and open clusters, whose magnitude and sky distributions are shown in Fig. 1.

Fig. 1 Sky and magnitude distribution of the sources in our two samples. Top: equatorial coordinate distribution of the binaries and cluster members with G < 14 in healpix with Nside = 8. The center is (α,δ) = (0,0), with α increasing toward the left. Bottom: G-magnitude distribution in equal bins and with the binning used in Sect. 3.3.

2.1 Binary stars

The two components of a widely separated binary have similar motions through the Galaxy and nearly identical parallaxes. El-Badry et al. (2021) used these two facts to select pairs of stars in EDR3 whose parallaxes and proper motions are consistent with binarity. Their catalog forms the basis of our binary star sample. We sought to measure a frame rotation at the level of tens of μas yr⁻¹. For a pair of Solar mass stars at 500 pc, the orbital motion is ~ 100 μas yr⁻¹, even at a separation of 10⁴ AU. Our measurement of the frame rotation required hundreds of widely separated binaries to average over this motion.

We selected stars from the El-Badry et al. (2021) catalog using thresholds in angular separation and parallax and requiring 1 μm⁻¹ < ν_eff < 2 μm⁻¹. The latter criterion excludes six-parameter astrometric solutions (Gaia Collaboration 2021). We required stars to be separated by at least 5′′ to minimizethe primary star’s effect on the secondary’s Gaia astrometry (Lindegren et al. 2021). We then required that a star’s rotational Keplerian velocity at its projected separation not exceed 0.2 mas yr⁻¹, assuming a system mass of 1 M_⊙. This last criterion is equivalent to $$2\pi \sqrt{\frac{{\varpi }^3}{\theta }}<0.2 \text{mas},$$(1)

with parallax ϖ and angular separation θ both measured in mas. Most binaries will have a nonzero eccentricity and will preferentially be observed at apastron. Their physical separations will also be larger than their projected separations. Both of these effects make the typical orbital motion a factor of a few smaller than our formal limit of 0.2 mas yr⁻¹. Finally, we wished to select only true binaries, but spurious matches become increasingly common at wide angular separations. We applied an additional cut: that the chance alignment probability supplied by El-Badry et al. (2021) be no higher than 5%.

The above cuts produced a sample of about 160 000 binaries. We studied the dependence of the reference frame on magnitude, anchoring the reference frame to that of the faint stars (G ≳ 14). Within this sample of binaries, about 55 000 have a primary brighter than G = 14 and a secondary fainter than G = 14. These 55 000 stars, further divided by the magnitude of the primary, form the basis of the analysis in Sect. 3.

2.2 Open clusters

Our open cluster sample uses the 1903 clusters listed in Cantat-Gaudin et al. (2020) that have members brighter than G = 14. These clusters are distant and weakly bound, making them suitable calibrators for the EDR3 reference frame. We selected stars from Gaia EDR3 within the central radius of each cluster, applying thresholds of 0.3 mas in parallax and 0.5 mas yr⁻¹ in proper motion around each cluster’s listed centroids.

Since Cantat-Gaudin et al. (2020) determined the clusters’ mean astrometric parameters using Gaia DR2, we recomputed the median proper motion and parallax using our initial sample selected from EDR3. We then repeated the query on EDR3 based on those updated values. After this second iteration we obtained about 358 000 stars, of which 37 000 are brighter than G = 14. We considered the reference proper motion of each cluster to be the median proper motion of its members fainter than G = 14.

Fig. 2 Median proper motion bias computed in healpix with Nside = 8 for the sources in our sample in the magnitude range 11 < G < 13, mapped in equatorial coordinates. The center is (α, δ) = (0, 0), with α increasing toward the left. The lower panels show the results after applying the correction derived in Sect. 3 using the prescription supplied in Appendix A.

3 Aligning to the ICRF

We began by visualizing the proper motion bias across the sky using our combined sample of binaries and open cluster members. We divided the sky into healpix tiles with Nside = 8 and computed the median proper motion of stars in the magnitude range 11 < G < 13 relative to the faint (G > 14) stars. The top panels of Fig. 2 show clear evidence for frame rotation, closely resembling that seen in DR2(cf. Fig. 6 of Brandt 2018). This likely reflects the magnitude dependence of the anchoring of the WC0 reference frame to the HIPPARCOS positions at epoch 1991.25 (ESA 1997; van Leeuwen 2007). In this section, we derive the magnitude dependence of this frame rotation and supply a prescription for its removal. The lower panel of Fig. 2 shows the results of applying this correction to our full sample of 11 < G < 13 stars.

3.1 The spin

We represent the ICRF as a vector triad C = [X Y Z], where X, Y, and Z are orthogonal unit vectors pointing toward (α,δ) = (0,0), (90°,0), and (0,90°), respectively. An arbitrary reference frame $\text{tens}\tilde{C}=[\tilde{X} \tilde{Y} \tilde{Z}]$ can be obtained by rotating C around the three axes X, Y, and Z, as shown in Fig. 3.

Following the formalism presented in Sect. 2 of Lindegren (2020), if the overall rotation ω = [ω_X ω_Y ω_Z] is small, the proper motion of a given source in the reference frame $\text{tens}\tilde{C}$ is related to its proper motion in C through the approximation: $$\left[\begin{array}{c}\widetilde{{\mu }_{\alpha *}}-{\mu }_{\alpha *}\\ \widetilde{{\mu }_{\delta }}-{\mu }_{\delta }\end{array}\right]=\left[\begin{array}{ccc}-\sin \delta \cos \alpha & -\sin \delta \sin \alpha & \cos \delta \\ \sin \alpha & -\cos \alpha & 0\end{array}\right]\left[\begin{array}{c}{\omega }_X\\ {\omega }_Y\\ {\omega }_Z\end{array}\right].$$(2)

The proper motion bias (Δμ_α*, Δμ_δ) introduced by the spin of the reference frame therefore depends on the (α, δ) coordinates of the source. Over the whole sky, Δμ_δ averages to zero, while Δμ_α* is on average proportional to ω_Z. This effect is illustrated in Fig. 2 for the stars with 11 < G < 13 and is consistent with the findings of Fabricius et al. (2021), who point out a low-significance offset of the μ_α* of bright stars, and no detectable offset in μ_δ.

Fig. 3 ICRF modeled as an orthogonal vector triad C = [X Y Z]. The transformation from C to $\text{tens}\tilde{C}=[\tilde{X} \tilde{Y} \tilde{Z}]$ is obtained by applying the spin ω = [ωX ωY ωZ], which is a combination of rotations around the axes X, Y, and Z.

3.2 Fitting the spin components

The proper motion difference between a bright source and its faint companion, or the bulk motion of the fainter cluster stars, is due to multiple factors. Part of the difference is due to measurement errors (assumed to be Gaussian) on the proper motions, which we denote as (σ_α*,b, σ_δ,b) and (σ_α*,f, σ_δ,f) for the bright and faint source, respectively. For a cluster member, σ_α*,f and σ_δ,f are negligible because they represent medians over many stars. An intrinsic velocity difference is also expected owing to the orbital motion of binaries and to internal velocity dispersion in clusters. We denote this intrinsic dispersion by σ₁.

We also accounted for the fact that both the binary and cluster star samples contain a fraction 1 − g of outliers that are background contaminants. Since our sources are distributed across the entire sky, the overall contributionof these contaminants does not skew the proper motion distribution, but it adds tails that cannot be reproduced with a single Gaussian component. Rather than performing, for example, a sigma clipping, we modeled the proper motion distribution as a two-component Gaussian mixture by adding a broader component of dispersion σ₂.

The three spin components ω_X, ω_Y, and ω_Z give a proper motion shift (Δμ_α*, Δμ_δ) at each star’s position according to Eq. (2). We optimized these three components by maximizing the likelihood given by the Gaussian mixture model $$\begin{array}{l}L={\displaystyle \prod _{\text{stars\hspace{0.05em}}i}(}\frac{g}{2\pi {\sigma }_{\alpha *,1}{\sigma }_{\delta ,1}}\exp \left[-\frac{{(D{\mu }_{\alpha *})}^2}{2{\sigma }_{\alpha *,1}^2}-\frac{{(D{\mu }_{\delta })}^2}{2{\sigma }_{\delta ,1}^2}\right]\\ +\frac{1-g}{2\pi {\sigma }_{\alpha *,2}{\sigma }_{\delta ,2}}\exp \left[-\frac{{(D{\mu }_{\alpha *})}^2}{2{\sigma }_{\alpha *,2}^2}-\frac{{(D{\mu }_{\delta })}^2}{2{\sigma }_{\delta ,2}^2}\right])\end{array}$$(3)

with, for example, $$\frac{{(D{\mu }_{\alpha *})}^2}{2{\sigma }_{\alpha *,1}^2}=\frac{{({\mu }_{\alpha *,b}-{\mu }_{\alpha *,f}-\Delta {\mu }_{\alpha *})}^2}{2({\sigma }_1^2+{\sigma }_{\alpha *,b}^2+{\sigma }_{\alpha *,f}^2)}.$$(4)

For the binary sample, we took σ₁ = 100 μas yr⁻¹, which is comparable to the typical orbital motion that we expect for these stars. For the clusters, σ₁ is the quadrature sum of 10 μas yr⁻¹ (which represents the effect of small-scale correlations on the proper motion precision) and the intrinsic proper motion dispersion (which corresponds to 0.5 km s⁻¹ at the cluster distance). For both samples we used $g=\frac{1}{2}$ and σ₂ = 0.3 mas yr⁻¹. Varying g, σ₁, and σ₂ within reasonable value ranges does not change our results, but it affects the uncertainties by ~ 10%.

Maximizing Eq. (3) gives the best-fit frame rotation ω via the proper motion shifts (Δμ_α*, Δμ_δ). To assign uncertainties to these values, we used bootstrap resampling within each magnitude bin. We report the mean and standard deviations of 400 bootstrap re-samples for each magnitude bin.

3.3 Spin parameters as a function of magnitude

The best-fit spin components that we obtained in different magnitude bins are shown in Fig. 4 and listed in Table 1. The values obtained from binaries and from clusters are consistent, although the cluster sample provides slightly larger uncertainties due to a lower number of sources and larger intrinsic velocity dispersions. Both categories of sources show that the bias becomes much more significant for G < 13, clearly indicating an instrumental origin for this effect. The bottom panel of Fig. 4 shows the results obtained with the combined sample of binaries and cluster stars.

The ω_X and ω_Y componentsof the spin are very significant for the sources brighter than G = 13, while ω_Z only appears significant at G < 10.5. The ω_Y component varies strongly with magnitude, reaching up to 80 μas yr⁻¹ for stars in the magnitude range G = 12–13. This bias is five times larger than the typical proper motion error quoted in the Gaia EDR3 catalog in this magnitude range.

Figure 5 shows running medians of the difference between a star’s proper motion and that of its cluster. This bias depends on position; we divided the sky in half to avoid having it average to zero. The median bias reaches ~ 40 μas yr⁻¹ for G ≲ 13, roughly double the formal EDR3 uncertainties, but would reach higher values if computed on smaller portions of the sky. The bottom row of Fig. 5 shows a correction of the proper motions of the cluster stars based on the spin computed from binary stars. We aligned the proper motions of stars brighter than G = 13 to the ICRF through the relation $${\left[\begin{array}{c}{\mu }_{\alpha *}\\ {\mu }_{\delta }\end{array}\right]}_{\text{ICRF}}={\left[\begin{array}{c}{\mu }_{\alpha *}\\ {\mu }_{\delta }\end{array}\right]}_{\text{EDR3}}-\left[\begin{array}{ccc}-\sin \delta \cos \alpha & -\sin \delta \sin \alpha & \cos \delta \\ \sin \alpha & -\cos \alpha & 0\end{array}\right]\left[\begin{array}{c}{\omega }_X\\ {\omega }_Y\\ {\omega }_Z\end{array}\right],$$(5)

where the values of ω_X, ω_Y, and ω_Z are the best-fit spin parameters obtained in the corresponding magnitude bin. We include in Appendix A a Python function that performs the correction calibrated on the combined sample of cluster members and binaries.

For bright magnitudes (G ≲ 10), the bias that we find likely reflects errors in the alignment of HIPPARCOS with the ICRF. These were estimated to be 0.6 mas (ESA 1997), or about 25 μas yr⁻¹ if dividing by the difference between the HIPPARCOS and Gaia catalog epochs. At G = 13, the bias is similar to the rotation of (ω_X, ω_Y, ω_Z) = (–16.6, –95.0, +28.3) μas yr⁻¹ applied to bright sources (Lindegren et al. 2021), with the signs reversed. This suggests that the systematics in the bright source reference frame build gradually with magnitude away from the boundary between WC0 and WC1.

Fig. 4 Components of the best-fit frame rotation derived from a sample of wide-separation binaries (top), open cluster members (middle), and both samples combined (bottom). The uncertainties are derived from bootstrap resampling within each magnitude bin.

Fig. 5 Proper motion bias before and after correction. Top panels: running medians of the differences between proper motions of bright open cluster members and the proper motions of their host clusters as traced by the faint stars. Bottom panels: same running medians after correcting the proper motions of the bright stars using the spin of the reference frame derived from binary stars. The shaded areas correspond to 1σ and 2σ intervals computed using bootstrap resampling.

3.4 Evidence for a color dependence

In order to calibrate color-dependent spread functions, the Gaia astrometric processing relies on prior color information given by the effective wavenumberν_eff. This is available in the EDR3 catalog for each source as nu_eff_used_in_astrometry. To investigate a color dependence of ω, we allowed the three spin parameters to vary linearly with ν_eff by expressing them as $$\{\begin{array}{l}{\omega }_X={\omega }_{X,0}+c_X({\nu }_{\text{eff}}-1.5)\hfill \\ {\omega }_Y={\omega }_{Y,0}+c_Y({\nu }_{\text{eff}}-1.5)\hfill \\ {\omega }_Z={\omega }_{Z,0}+c_Z({\nu }_{\text{eff}}-1.5)\hfill \end{array}\text{}.$$(6)

This adds three free parameters to the fitting procedure described in Sect. 3.2.

The best-fit values we obtain for these coefficients in each magnitude bin are shown in Fig. 6, along with the uncertainty estimated by bootstrapping. We obtain negative values for c_X and c_Y, indicating that the bias is stronger for smaller values of ν_eff (corresponding to redder sources). We find mostly positive values for c_Z.

To assess whether these values could be the result of statistical fluctuations, and whether c_X, c_Y, and c_Z are constant with magnitude, we applied χ² statistics. In Fig. 6, we report the p values corresponding to the hypothesis that the true value of c_X, c_Y, or c_Z is zero (p_no), and to the hypothesis that they are nonzero and constant with G magnitude (p_dep). We obtain small values for p_no, providing conclusive evidence that ω_X varies with ν_eff, as well as marginal evidence that ω_Y and ω_Z do so as well. The p_dep values indicate that, within the statistical precision allowed by our sample, c_X, c_Y, and c_Z do not appear to vary with G magnitude.

We next tested whether this color correction, computed mainly from the binary stars, improves the proper motion residuals for the clusters members. The binary stars and cluster members have systematically different ν_eff, with the bright cluster stars being more distant and bluer than the bright components of nearby wide binaries. The cluster stars also show a larger dispersion of ν_eff values.

We find that adding a color correction to the magnitude correction given by Eq. (5) slightly increases the dispersion of the cluster stars’ proper motion residuals (i.e., it degrades the quality of the correction). We hypothesize that a linear dependence of the spin parameters on ν_eff is a poor approximation for the wider range of colors represented by bright open cluster members: Eq. (5) tends to overcorrect these stars. Unfortunately, our sample is insufficient to fully characterize the dependence of the reference frame on color. We therefore recommend only the magnitude correction given by Eq. (5) but advise that color-dependent frame rotations are present at ~ 10 μas yr⁻¹ for a star somewhat bluer or redder than the Gaia median of ν_eff ≈ 1.5.

Fig. 6 Coefficients obtained in different magnitude bins when fitting color-dependent spin parameters as defined in Eq. (6). The shaded areas correspond to 1σ and 2σ uncertainties obtained by bootstrapping, and pno and pdep are defined in the text. While cX in particular is significantly nonzero, we do not recommend the application of a color correction proportional to νeff. Our sample lacks the statistical power to fully characterize the chromaticity of the spin.

4 Summary and conclusion

The large magnitude range covered by the Gaia observations requires the use of different gatings and readout modes. The reference frame of the faint sources can be calibrated directly to the ICRF using distant radio quasars that have an optical counterpart.Due to the unprecedented astrometric capabilities of the Gaia spacecraft, no reference external sample of sufficient quality is available for brighter sources; the bright observations are calibrated to the fainter ones to align them with the ICRF.

In this paper we have shown that the reference frame of the proper motions of the bright EDR3 sources (G < 13) rotates with respect to that of the faint EDR3 sources. The resulting proper motion bias reaches 80 μas yr⁻¹ in the magnitude range G = 11–13, five times larger than the nominal uncertainty listed in the Gaia EDR3 catalog for these sources. We have also shown evidence for a second-order dependence of the spin parameters on the color of the source when modeled as a linear function of the effective wavenumber ν_eff. A color correction proportional to ν_eff, however, proves unsatisfactory over the wide range of colors present in our open cluster sample. We recommend the use of only a magnitude-dependent correction and caution that color-dependent systematics of up to ~ 10 μas yr⁻¹ will still be present.

Our main result is a magnitude-based correction to be applied to the Gaia EDR3 proper motions of bright stars in order to align them with the ICRF. The correction can be computed using Eq. (5) with the spin parameters listed in Table 1. Appendix A provides sample Python code to realize this correction for a list of stars in Gaia EDR3.

Acknowledgments

We thank L. Lindegren, C. Fabricius, and F. Arenou for their feedback and helpful comments, as well as the referee and editor who helped improve this manuscript. TCG acknowledges support by the Spanish Ministry of Science, Innovation and University (MICIU/FEDER, UE) through grants RTI2018-095076-B-C21 and the Institute of Cosmos Sciences University of Barcelona (ICCUB, Unidad de Excelencia ‘María de Maeztu’) through grant CEX2019-000918-M. 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.

Appendix A Python implementation of the frame spin correction

This function applies Eq. (5) using the values listed in Table 1.

References

All Tables

All Figures

	Fig. 1 Sky and magnitude distribution of the sources in our two samples. Top: equatorial coordinate distribution of the binaries and cluster members with G < 14 in healpix with Nside = 8. The center is (α,δ) = (0,0), with α increasing toward the left. Bottom: G-magnitude distribution in equal bins and with the binning used in Sect. 3.3.
In the text

	Fig. 2 Median proper motion bias computed in healpix with Nside = 8 for the sources in our sample in the magnitude range 11 < G < 13, mapped in equatorial coordinates. The center is (α, δ) = (0, 0), with α increasing toward the left. The lower panels show the results after applying the correction derived in Sect. 3 using the prescription supplied in Appendix A.
In the text

	Fig. 3 ICRF modeled as an orthogonal vector triad C = [X Y Z]. The transformation from C to $\text{tens}\tilde{C}=[\tilde{X} \tilde{Y} \tilde{Z}]$ is obtained by applying the spin ω = [ωX ωY ωZ], which is a combination of rotations around the axes X, Y, and Z.
In the text

	Fig. 4 Components of the best-fit frame rotation derived from a sample of wide-separation binaries (top), open cluster members (middle), and both samples combined (bottom). The uncertainties are derived from bootstrap resampling within each magnitude bin.
In the text

Fig. 5 Proper motion bias before and after correction. Top panels: running medians of the differences between proper motions of bright open cluster members and the proper motions of their host clusters as traced by the faint stars. Bottom panels: same running medians after correcting the proper motions of the bright stars using the spin of the reference frame derived from binary stars. The shaded areas correspond to 1σ and 2σ intervals computed using bootstrap resampling.

In the text

Fig. 6 Coefficients obtained in different magnitude bins when fitting color-dependent spin parameters as defined in Eq. (6). The shaded areas correspond to 1σ and 2σ uncertainties obtained by bootstrapping, and pno and pdep are defined in the text. While cX in particular is significantly nonzero, we do not recommend the application of a color correction proportional to νeff. Our sample lacks the statistical power to fully characterize the chromaticity of the spin.

In the text