Issue 
A&A
Volume 649, May 2021
Gaia Early Data Release 3



Article Number  A4  
Number of page(s)  31  
Section  Astronomical instrumentation  
DOI  https://doi.org/10.1051/00046361/202039653  
Published online  28 April 2021 
Gaia Early Data Release 3
Parallax bias versus magnitude, colour, and position
^{1}
Lund Observatory, Department of Astronomy and Theoretical Physics, Lund University,
Box 43,
22100
Lund, Sweden
email: lennart@astro.lu.se
^{2}
Astronomisches RechenInstitut, Zentrum für Astronomie der Universität Heidelberg,
Mönchhofstr. 1214,
69120
Heidelberg, Germany
^{3}
HE Space Operations BV for European Space Agency (ESA), Camino bajo del Castillo, s/n, Urbanizacion Villafranca del Castillo,
Villanueva de la Cañada,
28692
Madrid, Spain
^{4}
Lohrmann Observatory, Technische Universität Dresden,
Mommsenstraße 13,
01062
Dresden, Germany
^{5}
European Space Agency (ESA), European Space Astronomy Centre (ESAC), Camino bajo del Castillo, s/n, Urbanizacion Villafranca del Castillo,
Villanueva de la Cañada,
28692
Madrid, Spain
^{6}
Vitrociset Belgium for European Space Agency (ESA), Camino bajo del Castillo, s/n, Urbanizacion Villafranca del Castillo,
Villanueva de la Cañada,
28692
Madrid, Spain
^{7}
Telespazio Vega UK Ltd for European Space Agency (ESA), Camino bajo del Castillo, s/n, Urbanizacion Villafranca del Castillo,
Villanueva de la Cañada,
28692
Madrid,
Spain
^{8}
Institute of Astronomy, University of Cambridge,
Madingley Road,
Cambridge
CB3 0HA, UK
Received:
12
October
2020
Accepted:
9
December
2020
Context. Gaia Early Data Release 3 (Gaia EDR3) gives trigonometric parallaxes for nearly 1.5 billion sources. Inspection of the EDR3 data for sources identified as quasars reveals that their parallaxes are biased, that is, they are systematically offset from the expected distribution around zero, by a few tens of microarcseconds.
Aims. We attempt to map the main dependences of the parallax bias in EDR3. In principle, this could provide a recipe for correcting the EDR3 parallaxes.
Methods. Quasars provide the most direct way for estimating the parallax bias for faint sources. In order to extend this to brighter sources and a broader range of colours, we used differential methods based on physical pairs (binaries) and sources in the Large Magellanic Cloud. The functional forms of the dependences were explored by mapping the systematic differences between EDR3 and DR2 parallaxes.
Results. The parallax bias is found to depend in a nontrivial way on (at least) the magnitude, colour, and ecliptic latitude of the source. Different dependences apply to the five and sixparameter solutions in EDR3. While it is not possible to derive a definitive recipe for the parallax correction, we give tentative expressions to be used at the researcher’s discretion and point out some possible paths towards future improvements.
Key words: astrometry / parallaxes / methods: data analysis / space vehicles: instruments / stars: distances
© ESO 2021
1 Introduction
The (early) Third Gaia Data Release (hereafter EDR3; Gaia Collaboration 2021b) provides trigonometrically determined parallaxes for approximately 1.478 billion sources in the magnitude range G ≃ 6–21. The sources include stars, unresolved binaries, compact extragalactic objects such as active galactic nuclei (AGNs), and other objects that appear roughly pointlike at the angular resolution of Gaia (~0.1 arcsec). Although Gaia in principle determines absolute parallaxes (e.g. Lindegren & Bastian 2011), imperfections in the instrument and data processing methods inevitably result in systematic errors in the published parallax values. For example, it is well known that the parallax solution is degenerate with respect to certain variations of the ‘basic angle’ between the viewing directions of the two Gaia telescopes (Butkevich et al. 2017). This means that based on the Gaia astrometric observations alone, we cannot simultaneously determine absolute parallaxes and calibrate this particular perturbation of the instrument. Conversely, when the actual instrument perturbations contain a component of this form, it will produce biased parallax values in the astrometric solution.
Since the second release of Gaia data in April 2018 (DR2; Gaia Collaboration 2018b), a number of investigations have been published that have resulted in estimates of the parallax systematics in DR2 from a variety of astrophysical objects (see e.g. Chan & Bovy 2020, and references therein). Reported values of the Gaia DR2 ‘parallax zeropoint’ (i.e. the quantity to be subtracted from the DR2 parallaxes) range from about − 30 μas to − 80 μas. While quasars in DR2 yield a median parallax of − 29 μas (Lindegren et al. 2018), the cited studies, which typically use much brighter objects, tend to give more negative values. There is consequently a strong suspicion that the parallax offset in DR2 depends on the magnitude, and possibly also on the colour of the sources (e.g. Zinn et al. 2019). Furthermore, it is known that Gaia DR2 has positiondependent offsets of the parallaxes on angular scales down to ~1 deg (Arenou et al. 2018, Lindegren et al. 2018). The systematic offset of parallaxes in DR2 could therefore be a complicated function of several variables x, including atleast the magnitude, colour, and position of the object; in the following we write this Z_{DR2} (x).
In Gaia EDR3 quasars have a median parallax of about − 17 μas, and already a simple plot of the parallaxes versus the G magnitude or colour index reveals systematic variations at a level of ~10 μas (see Sect. 4.1). Positiondependent variations are also seen on all angular scales, although with smaller amplitudes than in DR2 (Lindegren et al. 2021). Thus it is possible to define an offset function Z_{EDR3}(x) that is in general different from Z_{DR2}(x), although it may depend on the same variables x.
The letter Z adopted for these functions is a mnemonic for ‘zeropoint’. What is meant is an estimate of the bias (orsystematics) of the parallax estimate as a function of certain known variables x. The practical determination of this function is fraught with difficulties and uncertainties, and even if we knew the true parallax of every object in the catalogue, it would not be possible to define a unique function Z(x) short of tabulating the bias for every source. The function will necessarily depend on, for example, the choice of arguments in x and on their resolution and numerical representation. At best, we can achieve a prescription such that the use of ϖ_{i} − Z(x_{i}) for the parallax of source i, instead of the catalogue value ϖ_{i}, will in most cases give more accurate and consistent results in astrophysical applications.
The aim of this paper is to map some of the main dependences of the Gaia EDR3 parallax bias (zeropoint), as found in thecourse of the internal validations carried out by the Gaia astrometry team prior to the publication of the data. Because the parallax determinations in Gaia EDR3 are of two kinds, known as five and sixparameter solutions, with distinctly different properties, the results are given in the form of two functions Z_{5} (x) and Z_{6} (x) describing the bias as a function of magnitude, colour, and position for each kind of solution. In the following, the functions Z_{5} and Z_{6} always refer to EDR3, while for DR2 there is only one kind of solution, with the bias function Z_{DR2}.
Section 2 is a brief overview of some aspects of the Gaia instrument and data processing that are particularly relevant for the bias functions, such as the different observing modes depending on the magnitude of the source, and the use of colour information in the five and sixparameter solutions. In Sect. 3 we discuss the systematic differences between DR2 and EDR3 parallaxes. While the DR2 parallaxes are completely superseded by the later release, the differences could give important clues to the systematic dependences. In Sect. 4 we estimate Z_{5}, the bias function for the fiveparameter solutions in EDR3. In Sect. 5 we use a differential procedure to estimate Z_{6}, the bias function for the sixparameter solutions. A limited test of the derived bias functions is made in Sect. 6. Some possible future improvements are discussed in Sect. 7 before the findings are summarised in Sect. 8. Certain technical details are provided in the appendices.
2 Observing and processing modes for the EDR3 astrometry
The determination of biases in Gaia astrometry must rely on a comparison with external data. In principle, this process requires no prior knowledge on how the astrometry was generated, but the mapping of possible dependences is surely facilitated by some acquaintance with certain aspects of the instrument and data processing. This section provides a brief overview of relevant topics.
Fig. 1 Fraction of observations in different modes. The shaded area represents gated observations; the curves show the fraction of AF observations made in window classes WC0a, WC0b, WC1, and WC2. 
2.1 Window classes and gates
The Gaia instrument and its routine operations are described in Sects. 3 and 5 of Gaia Collaboration (2016). Further details of the initial treatment of the data and the subsequent astrometric processing can be found in Rowell et al. (2021) and Lindegren et al. (2021). The astrometric results are derived almost exclusively from observations made with the 62 CCDs occupying the central part of the focal plane assembly, the socalled astrometric field (AF; see Fig. 4 in Gaia Collaboration 2016). Of relevance here is that only small patches (‘windows’) of the CCD images around detected point sources are transmitted to the ground, and that different sampling schemes (window size, pixel binning, etc.) are used depending on the onboard estimate of the brightness of the source. Furthermore, to avoid pixel saturation for sources brighter than G ≃ 12 (where G is the magnitude in the Gaia passband; Evans et al. 2018), the CCD integration time may be reduced by means of gates (Crowley et al. 2016).
In principle, the different sampling schemes (known as window classes), and their different uses together with the gates, may require separate calibrations both in the initial treatment (image parameter determination, IPD) and in the astrometric solution. This need arises both because of subtle differences in the way the CCDs and associated electronics function depend on their mode of operation (Hambly et al. 2018), and because different data processing models are used, for example depending on whether the window contains a one or twodimensional image. Further complications are caused by unavoidable conflicts in the onboard resource management, for example overlapping and truncated windows.
In the EDR3 astrometric solution (Lindegren et al. 2021), distinct calibration models were used for observations in the four window classes designated WC0a, WC0b, WC1, and WC2. The first two provide twodimensional images with different gate usages, while WC1 and WC2 provide onedimensional images (i.e. the pixels have been binned in the other dimension) of different sizes. For the IPD a similar division was made, but without the distinction between WC1 and WC2. Figure 1 shows the fraction of AF observations in each window class as a function of G. Clearly, the window classes map out distinct intervals in G, but the transitions around G = 11, 13, and 16 are slightly fuzzy because the decision on the sampling scheme is based on the onboard realtime estimate of G, which differs from the onground estimated mean magnitude used in the plot. The transition between WC0a and WC0b occurs gradually over a whole magnitude because it involves gating thresholds that are set individually for the CCDs.
In EDR3, the astrometric parameters for a given source are typically obtained by combining some 200–500 individual AF observations. Depending on the mean magnitude of the source, the observations will be distributed among the window classes and gates in proportion to the fractions shown in Fig. 1. This explains some of the magnitudedependent effects described in later sections. As it turns out, it is also relevant whether the observations are gated or not. This creates another transition around G = 12, indicated by the shaded area in the diagram.
2.2 Five and sixparameter solutions, and the use of colour information
For a source in EDR3 that has a parallax, the astrometric parameters were determined either in a fiveparameter solution or in a sixparameter solution. The choice of solution depends on the availability of colour information, in the form of an effective wavenumber ν_{eff}, at the time when the IPD is performed. The effective wavenumber comes from the processing of BP and RP spectra in the photometric pipeline (Riello et al. 2021; De Angeli et al., in prep.) In IPD, a calibrated pointspread function (PSF) or linespread function (LSF) is fitted to the CCD samples in a window, yielding the accurate location and flux (in instrumental units) of theimage in the pixel stream. In a previous step, the PSF and LSF were calibrated as functions of several parameters, including window class and effective wavenumber.
A fiveparameter solution is computed if the source has a reliable value of ν_{eff} that can be used to select the appropriate PSF or LSF for the IPD. This means that most of the colourdependent image shifts, caused by diffraction phenomena in the telescopes and electronic effects in the CCDs, have been eliminated already in the data input to the astrometric solution (AGIS). In AGIS, the image locations for a given source are fitted by the standard astrometric model with five unknowns, α, δ, ϖ, μ_{α*}, and μ_{δ} (Lindegren et al. 2012). These parameters (for some nearby stars, also the spectroscopic radial velocity) are sufficient to describe the astrometric observations of wellbehaved sources such as single stars and quasars.
When a source does not have a reliable ν_{eff} at the time of the IPD, it will instead obtain a sixparameter solution in AGIS, where the additional (sixth) parameter is an astrometric estimate of ν_{eff} known as the pseudocolour, denoted . At the IPD stage, image locations and fluxes for the source are determined by fitting the PSF or LSF at the default effective wavenumber m^{−1}. This value isclose to the mean ν_{eff} for the faint sources that make up the majority of sources without a photometric ν_{eff}. The use of the default colour in IPD results in image locations that are biased in proportion to , where ν_{eff} is the actual (but unknown) colour. The coefficient of proportionality is a property of the instrument, known as the chromaticity; it varies with time and position in the field, but can be calibrated in the astrometric solution by means of stars of known colours (for details, see Lindegren et al. 2021). Using the calibrated chromaticity, corrections to the default colour can be estimated for the individual sources, which gives their pseudocolours. The determination of pseudocolour thus takes advantage of the (extremely small) alongscan shifts of the image centre versus wavelength, caused by optical wavefront errors and other imperfections in the astrometric instrument.
The instrument chromaticity is calibrated in AGIS by means of a special solution, using a subset of about eight million primary sources for which IPD was executed twice: once using the actual effective wavenumber ν_{eff}, known from the spectrophotometric processing, and once using the default value . For these sources it is possible to compute both fiveparameter solutions (which are the published ones) and sixparameter solutions; the latter are not published but used in Sect. 5 to derive the systematic differences between the five and sixparameter solutions.
Because the parallax bias is colour dependent, the functions Z_{5} (x) and Z_{6} (x) need to include a colour parameter among their arguments in x. Rather than using, for example, the colour index G_{BP} − G_{RP}, which is not available for all sources, the colour parameter used here is the (photometric) effective wavenumber ν_{eff} (nu_eff_used_in_astrometry) in Z_{5}, and the (astrometric) pseudocolour (pseudocolour) in Z_{6}. By definition, this colour information is available for all sources with an astrometrically determined parallax.
The use of colour information in IPD and AGIS has one additional feature that significantly affects the fiveparameter solutions with ν_{eff} > 1.72 μm^{−1} (G_{BP}−G_{RP} ≲ 0.14) or ν_{eff} < 1.24 μm^{−1} (G_{BP}−G_{RP} ≳ 3.0). The calibration of the PSF and LSF versus colour is only made for the wellpopulated interval 1.24 < ν_{eff} < 1.72 μm^{−1}, where a quadratic variation of the displacement with ν_{eff} is assumed; if the IPD requests a PSF or LSF for a wavenumber outside of this range, then the calibration at the nearest boundary is used. This ‘clamping’ of the wavenumbers guarantees that the LSF or PSF model used by the IPD is always sensible, but on the downside, it introduces some biases for sources of extreme colours. In the astrometric calibration model, the chromaticity is assumed to be linear over the entire range of wavenumbers. The combination of the two models may produce astrometric biases that are strongly nonlinear in wavenumber, and there could in particular be discontinuities in ∂Z_{5} ∕∂ν_{eff} at ν_{eff} = 1.24 and 1.72 μm^{−1}. As no clamping is used for sources with sixparameter solutions, such abrupt changes in slope are not expected for Z_{6} versus pseudocolour, but the relation can still be nonlinear from other effects.
Because ν_{eff} is computed from the detailed BP and RP spectra, while G_{BP} − G_{RP} depends on the ratio of the integrated fluxes in the two bands, there is no strict onetoone relation between the two quantities. Nevertheless, for − 0.5 ≤ G_{BP} − G_{RP} ≤ 7 the following simple formulae
represent the mean relation for stellar objects to within ± 0.007 μm^{−1} in the effective wavenumber (Lindegren et al. 2021). We note that the values of ν_{eff} used in EDR3 actually come from the spectrophotometric processing of the previous cycle, corresponding to the DR2 photometry. This is a consequence of the overall cyclic processing scheme of DPAC (Fabricius et al. 2016), in which the IPD is one of the first processes in a cycle and the photometric processing is further downstream.
Fig. 2 Celestial maps in ecliptic coordinates of some selected statistics in Gaia EDR3. Panel a: mean parallax uncertainty for sources with fiveparameter solutions and G < 14 mag. Panel b: mean correlation coefficient between parallax and pseudocolour for sources with sixparameter solutions. Panel c: smoothed median parallax of quasars, corrected for the global median offset of − 17 μas. Owing to the small number of quasars identified at small Galactic latitudes, no data are shown for  sin b  < 0.1. The maps use a Hammer–Aitoff projection in ecliptic coordinates with λ = β = 0 at the centre, ecliptic north up, and ecliptic longitude λ increasing from right to left. All three statistics exhibit some degree of systematic dependence on ecliptic latitude, which may be symmetric (as in a) or antisymmetric (as in b and c), but also even larger variations as functions of both coordinates. 
2.3 Scanning law
The two telescopes in Gaia are continuously scanning the celestial sphere according to a predefined schedule known as the scanning law. Details are given in Sect. 5.2. of Gaia Collaboration (2016). For thermal stability, it is necessary that the spin axis is kept at a fixed angle (45°) to the Sun at all times, and as a consequence, the scanning pattern is roughly symmetric about the ecliptic, at least in a statistical sense and after several years of observations. Thus, although equatorial (ICRS) coordinates are used throughout the processing and for the astrometric end products, many characteristics of the data are (approximately) aligned with ecliptic coordinates rather than equatorial.
Three examples of the ecliptic alignment are shown in Fig. 2. In the top panel a, theprecision of the parallaxes is shown to depend systematically on the ecliptic latitude (β), with 50–60% higher uncertainties along the ecliptic than near the ecliptic poles. The middle panel b shows the mean correlation coefficientbetween parallax and pseudocolour in the sixparameter solutions. This correlation is systematically positive for β ≳ 45° and negative for β ≲−45°, which is relevant for the parallax bias because an offset in the assumed colours of the sources translates into a parallax bias that is proportionalto the correlation coefficient. Although this correlation coefficient is only available for the sixparameter solutions, the correlation between the errors in colour and parallax exists also for sources with fiveparameter solutions. For example, if ν_{eff} is systematically too high in the fiveparameter solutions, the correlations will produce a more positive parallax bias in the (ecliptic) northern sky than in the south. At intermediate latitudes the correlation coefficient exhibits more complex (and larger) variations related to the scanning law. The bottom panel c is a smoothed map of quasar parallaxes, increased by a constant 17 μas to compensate for the global offset. Large regional variations are seen at middle latitudes, but a systematic difference between north and south is clearly visible for  β ≳ 45°. Although such a systematic could be produced by the correlation mechanism just described, several other explanations are possible.
The maps in panels a and b of Fig. 2 show substantial regional and local variations, especially for  β  ≲45°. These features are related to the scanning law and the (still) relatively poor coverage of the ecliptic zone obtained in the 33 months of data used for the EDR3 astrometry (global coverage is optimised for a mission length of 60 months). It can be expected that the parallax bias has regional and local variations of a similar character, and the map of quasar parallaxes (panel c) seems to confirm this, although the finer details are made invisible by the smoothing, and smallnumber statistics contribute to the variations most clearly along the Galactic zone of avoidance. In practice, it is very difficult to determine local or even regional variations in Z_{5} and Z_{6} with any degree of confidence. Even if we trust (some of) the variations seen in the quasar map, they are probably only representative for the faintsources with colours similar to those of the quasars, that is, bluer than the typical faint Galactic stars.
Thus, while a detailed mapping of Z_{5} and Z_{6} versus position is not possible, we expect on theoretical grounds a global variation with ecliptic latitude, which is also seen empirically. The only positional argument in Z_{5} and Z_{6} is therefore taken to be the ecliptic latitude, β. In the Gaia archive this coordinate is given (in degrees) as ecl_lat. Alternatively, it can be computed to within ± 43 mas using the formula (3)
In this paper, all expressions involving the ecliptic latitude are written in terms of sin β.
3 Systematic differences between EDR3 and DR2 parallaxes
Although it is not our goal to study the parallax bias for Gaia DR2, that is, the function Z_{DR2}(x), it is instructive to start this investigation by examining the systematic parallax difference between EDR3 and DR2, that is, the differential bias function (4)
In contrast to either Z_{5} or Z_{DR2}, this difference can easily be mapped in considerable detail simply by computing mean differences in parallax for sources having the same identifier (source_id) in the two releases. (We ignore here the evolution of the source list between releases, which implies that a small fraction of the sources in EDR3 may be found in DR2 under a different source_id; see Gaia Collaboration 2021b for details.) Obviously, ΔZ contains no direct quantitative information on Z_{5}, but it is reasonable to expect that it gives a good qualitative indication of the relevant dependences on G, ν_{eff}, and β. The definition of the general parametrised function Z(G, ν_{eff}, β) in Appendix A is largely guided by the shape of ΔZ.
In Eq. (4) we use Z_{5} rather than Z_{6} because there are too few sources with sixparameter solutions for G < 13 that also appear in DR2. The colour argument is ν_{eff} given by nu_eff_used_in_astrometry in EDR3, which is available for all fiveparameter solutions.
Figure 3 shows the mean parallax difference ϖ_{EDR3} − ϖ_{DR2} for sources with fiveparameter solutions in EDR3, subdivided by colour. The top panel shows the dependence on magnitude, while the middle and bottom panels show examples of the dependence on β in two different magnitude intervals. In Fig. 4 we similarly show the mean difference as a function of colour, but subdivided by magnitude and ecliptic latitude. (In this and all other diagrams with effective wavenumber or pseudocolour on the horizontal axis, the direction of the axis is reversed to follow the usual convention for colourmagnitude diagrams, that is, blue objects to the left and red to the right.) These figures and all subsequent results on ΔZ were computed using all common sources for G < 14 and a geometrically decreasing random fraction of the fainter sources. In total about 26.7 million sources were used. Mean values of the parallax differences were computed using weights proportional to .
The plots in Figs. 3 and 4 show complex dependences on G, ν_{eff}, and β, with interactions among all three arguments (i.e. the colour dependence is different depending on both G and β, etc.). We use these results to design a continuous, multidimensional function, relevant parts of which are used in Sects. 4 and 5 to represent the functions Z_{5} and Z_{6} fitted to the EDR3 data. Owing to the scarcity of data available for these fits, the general function has to be quite schematic and cannot take many of the details seen in ΔZ into account, especially for the brightest sources. With this purpose in mind, the following features are noteworthy.
When ΔZ is considered as a function of G, there are jumps at G ≃ 11.0, 12.0, 13.0, and 16.0, corresponding to the boundaries shown in Fig. 1. The transitions are not abrupt, but occur over an interval of 0.2–0.4 mag. Features at G ≲ 9 are ignored in the following because the number of sources is too small for a reliable mapping of the parallax bias. Considering ΔZ as a function of ν_{eff}, the effect of the clamping at 1.24 and 1.72 μm^{−1} (Sect. 2.2) is visible in some plots as an abrupt change in slope; within these limits, the relation is approximately linear for 1.48< ν_{eff} < 1.72 μm^{−1} and curved(cubic) for 1.24< ν_{eff} < 1.48 μm^{−1}, with no visible break at 1.48 μm^{−1}. The hook for the reddest stars (ν_{eff} < 1.18 μm^{−1}) seen for G < 16 mag is not related to any known feature of the instrument or data processing, and because it concerns very few objects, it is ignored in the following. Finally, as a function of β, ΔZ typically shows an approximately linear or quadratic dependence on sinβ.
The parametrised function described in Appendix A is designed to take these features into account. Using this form, the differential bias is approximated by (5)
where c_{j} and b_{k} are basis functions in ν_{eff} and β, specified by Eqs. (A.3) and (A.4), respectively, and q_{jk} (G) are piecewise linearfunctions of G given by Eqs. (A.2) and (A.6) in terms of the fitted coefficients z_{ijk} in Eq. (A.1).
The approximation in Eq. (5) has 13 × 5 × 3 = 195 free parameters, namely all possible combinations of i = 0…12, j = 0…4, and k = 0…2 in the fitted coefficients z_{ijk}. A simultaneous leastsquares estimation of all 195 parameters shows that many of them are poorly determined and contribute little to the overall fit. To avoid overfitting, the following procedure is used. First, all 195 parameters are estimated (ẑ_{ijk}), along with their formal uncertainties (σ_{ijk}). The parameter with the smallest standard score S_{ijk} = ẑ_{ijk}∕σ_{ijk} is then removed (set to zero), and a new fit calculated with updated uncertainties. This procedure is repeated until S_{ijk} > 2 for all the retained parameters. However, the parameters z_{ijk} with j = k = 0 are always retained at their estimated values, regardless of their scores, to avoid that the sum in Eq. (5) defaults to zero at some G independently of ν_{eff} and β.
This procedure applied to the EDR3–DR2 parallax differences yields the 137 nonzero coefficients shown in Table 1. For compactness, no uncertainties are given; all values are significant at least on the 2σ level, although some, indicated by a colon inthe table, are below 3σ. The precise values of the coefficients, as well as the subset of significant coefficients, will depend on the selection of sources used in the fit, which was increasingly nonexhaustive for G > 14. If we had used more of the faint sources in EDR3 and DR2, we would undoubtedly have obtained a better determination of the coefficients towards the faint end, and would have increased the number of significant coefficients.
As a verification of the fit, the function ΔZ defined by the coefficients in Table 1 was evaluated for each of the 26.7 million sources used in the fit, and mean values in bins of magnitude, colour, and ecliptic latitude were computed in exactly the same way as was done for the parallax differences. The results, shown by the dashed lines in Figs. 3 and 4, thus represent our multidimensional fit to the mean parallax differences shown as circles in the same diagrams. The overall fit is reasonable, but we comment on a few systematic deviations below.
First, we note that the mean parallax differences around G ≃ 7.0 and ≃ 8.0 show large (10–15 μas) positive and negative excursions that are not modelled by the fit (Fig. 3, top). This discrepancy is a direct consequence of our choice to use a linear dependence on G in the interval 6.0–10.8, which in turn is motivated by the scarcity of data for estimating Z_{5} in this interval (see Sect. 4.4.2). Furthermore, as a function of β (Fig. 3, middle and bottom) there are local excursions from the general quadratic trend on a level of ± 5 μas. These are related to localised features on the sky (cf. Fig. 2), whereas the fitted model is only meant to represent the position dependences on the largest scales. Finally, as a function of ν_{eff} (Fig. 4), the variation between 1.48 and 1.72 μm^{−1} is not linear, as assumed in Eq. (A.3), which sometimes gives deviations of 5–10 μas for moderately blue sources (ν_{eff} ≃ 1.65 μm^{−1}). In the red end the curvature of the cubic segment is not always sufficient to give a good fit around ν_{eff} ≃ 1.25 μm^{−1}.
These discrepancies could be caused by specific features in DR2, EDR3, or both. It is of course also possible that the two data sets have similar systematics that cancel in the parallax difference, and which therefore have not been identified and included in the adopted model (Appendix A). Nevertheless, as shown in the next sections, the model seems to be general and flexible enough to describe Z_{5} and Z_{6} to the level of detail permitted by the data.
Fig. 3 Mean difference in parallax between EDR3 and DR2. Top: parallax difference vs. magnitude. Middle and bottom: parallax difference vs. ecliptic latitude in two magnitude intervals. Circles connected by solid lines are weighted mean values computed in bins of variable size, with at least 1000 sources per bin; dashed lines are mean values of the fitted parametrised function ΔZ (Eq. (5)), binned as for the sources. Filled red circles are for ν_{eff} < 1.48 μm^{−1}, and open blue circles for ν_{eff} > 1.48 μm^{−1}. 
Fig. 4 Mean difference in parallax between EDR3 and DR2 as a function of effective wavenumber for several ranges of the G magnitude.Circles connected by solid lines are weighted mean values computed in bins of variable size, with at least 1000 sources per bin, and dashed lines are mean values of the fitted parametrised function ΔZ (Eq. (5)), binned as for the sources. Filled red circles are for sources with β > 0, and open blue circles for β < 0. The vertical dashed lines mark the breakpoints for the basis functions c_{j}(ν_{eff}) in Eq. (A.3), i.e. the clamping limits at 1.24 and 1.72 μm^{−1} and the midpoint at 1.48 μm^{−1}. 
Coefficients for the function ΔZ(G, ν_{eff}, β) fitted to the difference in parallax between EDR3 and DR2.
4 Fiveparameter solutions
The main methods for estimating Z_{5}(G, ν_{eff}, β) in this paper are (i) quasars (Sect. 4.1), which provide a direct estimate of the parallax bias for G ≳ 14, although only in a rather narrow interval of ν_{eff}; (ii) stars in the Large Magellanic Cloud (LMC; Sect. 4.2), which map differential variations in the bias over a range of magnitudes and colours, although only in a specific location near the south ecliptic pole; and (iii) physical pairs or binaries (Sect. 4.4), which can be used to map the differential variations for the bright stars. Additionally, we use red clump (RC) stars for a differential study of the bias in the red, where the number of quasars is too small (Appendix C). In subsequent sections, these methods are developed in full detail.
In this process no assumption is made concerning the distance to the LMC, or the absolute magnitude of the RC stars. These objects are used purely differentially, and our estimates of the parallax biases are completely anchored in the quasars, that is, placed on an absolute scale by means of their parallaxes.
4.1 Quasars
As part of EDR3, the Gaia Archive^{1} contains the table agn_cross_id, listing a total of 1 614 173 sources constituting a very clean sample of quasarlike objects, whose positions and proper motions in EDR3 formally define the reference frame of the catalogue, GaiaCRF3. The list was constructed by crossmatching the full EDR3 catalogue with 17 external catalogues of AGNs followed by a filtering based on the quality of the solutions and the astrometric parameters: the proper motions are consistent with zero, and the parallaxes with a constant offset of − 17 μas, to within five times their respective formal uncertainties. Full details of the selection procedure, and further characterisation of the sample, are given in Gaia Collaboration (2021c). In spite of the strict selection criteria, it is likely that the list contains a small number of stellar contaminants. As described below, we take this into account in a statistical way.
The quasar sample used here is a subset of agn_cross_id, consisting of 1 107 770 sources with fiveparameter solutions and effective wavenumbers (ν_{eff}) in the range 1.24 to 1.72 μm^{−1}. (For this investigation we used a preliminary version of agn_cross_id, resulting in a slightly smaller sample than the corresponding selection from the published table.) The median ν_{eff} is 1.59 μm^{−1}, with 1st and 99th percentiles at 1.44 and 1.69 μm^{−1}, which makes this sample significantly bluer than typical stars of similar magnitudes, and covering a smaller range of colours. The magnitudes range from G ≃ 13.4 to 21.0; only 97 are brighter than G = 15 and 541 are brighter than G = 16.
Figure 5 shows the weighted mean parallax plotted against G, ν_{eff}, and β. On the whole, there is a negative parallax bias: the weighted mean parallax of the sample used here is approximately − 21 μas and the median is about − 17 μas. For reference, these values are indicated by the red lines in Fig. 5. Apart from this offset, the main trends are as follows. (i) As a function of G, there is a clearly nonlinear variation with an approximately linear increase from G ≃ 17 to 20, with plateaus on either side of this interval or perhaps a decreasing trend for G > 20. (ii) As a function of ν_{eff}, the variation is approximately linear in the wellpopulated range of colours. If there is a curvature similar to what is seen in the EDR3–DR2 differences (Fig. 4) or the LMC data (Fig. 10), the interval in ν_{eff} covered by the quasars is too narrow to reveal it. (iii) As a function of β, the variation may be described by a quadratic polynomial in sinβ.
Similarly to what was observed in the EDR3–DR2 differences in Sect. 3, interactions among the three main variables are also visible in the quasar parallaxes. For example, if the quasars are binned by effective wavenumber and a quadratic polynomial a_{0} + a_{1} sinβ + a_{2} sin^{2}β is fitted to the parallaxes in each bin, it is found that a_{1} has a strong dependence on ν_{eff}, whereas no clear trend is seen for a_{2} (Fig. 6).
The trends described above, including interactions, are well approximated by the parametrised function Z(G, ν_{eff}, β) defined in Appendix A. However, the limited supports in G and ν_{eff} make it necessary toconstrain the general expression in Eq. (A.1) in several ways. Specifically, for the basis functions in magnitude, g_{i} (G) in Eq. (A.2), we only use i = 6…12 (i.e. z_{ijk} is not fitted for i < 6), and for the basis function in colour, c_{j}(ν_{eff}) in Eq. (A.3), we only use j = 0 and 1. For G < 17.5 there are not enough quasars to reliably determine a linear variation with G, as permitted by the basis functions; instead we assume that the bias is independent of G in the intervals 13.1< G < 15.9 and 16.1< G < 17.5, but allow a step around G = 16 representingthe transition from window class WC1 to WC2 (Sect. 2.1)^{2}. Additionally, the coefficients q_{12} were all found to be insignificant and were therefore constrained to zero. The resulting fit is given in Table 2.
At this point it is necessary to consider to what extent the fit is affected by a possible contamination of the quasar sample by Galactic stars. The contaminating stars will on average have higher measured parallaxes than the quasars of similar magnitude, thus biasing the fitted function towards more positive values. The effect is only expected to be important at the faint end and where the star density is high. In order to explore this, we introduce the confusion factor (6)
where D_{21} is the mean density of sources brighter than G = 21 in the vicinity of the quasar, expressed in deg^{−2}. Densities are calculated by counting the total number of sources in EDR3 in pixels of solid angle 0.8393 deg^{2} (healpix level 6), divided by the solid angle. Because the density of faint sources in EDR3 roughly doubles with each magnitude (d log_{10}D_{G}∕dG ≃ 0.3), the second term in Eq. (6) is the approximate change in density with G. Thus X is simply a convenient proxy for log_{10}D_{G}, the density of sources brighter than the quasar. For our quasar sample, X ranges fromabout 1.25 to 5.5, with a median at 3.4. Less than 1% of the quasars have X > 4.5.
When the mean quasar parallax is plotted versus X, there is a strong increase over practically the whole range of X values. Thisis not primarily caused by contamination, but rather by the positive trend of parallaxes versus G shown in Fig. 5a, transferred to X by the second term in Eq. (6). In Fig. 7 we plot instead the residual in parallax from the regression in Table 2 versus the confusion factor. Here, the mean residual is close to zero for X ≲ 3.7, after which it increases roughly linearly with X as suggestedby the dashed line. Based on this plot, we assume that the contamination bias at a particular position is proportional to (7)
with X given by Eq. (6). A dependence on position is expected not only from the varying star density, as encoded in X, but also from variations in the quality of Gaia astrometry created by the scanning law (see e.g. several plots in Lindegren et al. 2021 and Fabricius et al. 2021). Globally, the precision and number of observations improve towards the ecliptic poles, which to first order can be described as a linear dependence on sin^{2}β. We consequently model the contamination bias in the mean quasar parallax by adding the nuisance terms (8)
to the fitted model. Here, r_{0}(G) and r_{1}(G) represent the contamination bias at β = 0 and β = ±90°, respectively; both functions are piecewise linear functions of G using the basis functions g_{i}(G) in Eq. (A.2). Only the last two basis functions (i = 11 and 12) are used, as the rcoefficients are quite insignificant for G ≲ 20.
The resulting fit, including the contamination terms r_{0} and r_{1}, is given in Table 3. Compared with Table 2, where the fit did not include these terms, the global contamination bias (as shown by the difference in q_{00}) is about + 4 μas at G = 21.0 and below 1 μas at G = 20.0. Figure 8 is a map of the bias, calculated from Eq. (8) and averaged over the quasars at each location. The expected increase in bias towards the Galactic plane is very evident, but several features related to the different surveys contributing to the sample appear as well. These features probably reflect variations in the magnitude completeness, made visible through the steep increase in estimated bias towards G = 21.0.
Our best estimate of Z_{5}(G, ν_{eff}, β) from the quasar sample is therefore given by the coefficients q_{jk} in Table 3. The coefficients r_{0} and r_{1} must be ignored, as they represent the contamination bias, which should not be included in Z_{5}.
Several interesting observations can be made concerning the results in Table 3. The coefficient q_{00} represents the mean quasar parallax at the given magnitudes, averaged over the celestial sphere and reduced to the reference wavenumber, ν_{eff} = 1.48 μm^{−1}. As described in Appendix A, linear interpolation between the tabulated values should be used for other magnitudes. Its values agree rather well with the mean relation displayed in Fig. 5a. Of the remaining coefficients, the most important ones (in terms of how much they decrease the chisquare of the fit) are q_{02} and q_{11}. q_{02} describes a quadratic dependence on sinβ; the consistently negative sign means that the parallax bias is more negative towards the ecliptic poles, and this effect is strongest at the faint end. The interaction coefficient q_{11} is consistently positive and increases with magnitude; in terms of the chisquare, it is much more important than the corresponding simple coefficients q_{01} and q_{10}. At the median quasar colour, ν_{eff} = 1.59 μm^{−1}, q_{11} describes a clear north–south asymmetry of the parallaxes, which is what is seen in Fig. 5c. At the reference wavenumber 1.48 μm^{−1} the asymmetry, given by q_{01}, is smaller and less consistent. The coefficients q_{10} represent the colour gradient averaged over the celestial sphere. They are all consistent with a mean value of − 15 μas μm, which is much smaller than the slope ≃ −55 μas μm indicated by Fig. 5b. This apparent contradiction is explained by a correlation between magnitude and colour in the quasar sample: the faint quasars are on average redder than the brighter ones; together with the overall trend in Fig. 5a, this creates a stronger variation with colour in the sample as a whole than is present at a fixed magnitude. Together with the strong variation of q_{02} and q_{11} with magnitude, this illustrates the many complex dependences in the data and the difficulty of determining a unique function Z_{5} based on a limited sample with intrinsic correlations. It also shows the danger of using simple plots of the quasar parallaxes versus a single quantity such as colour for inferences on the parallax bias.
While Table 3 (ignoring r_{0} and r_{1}) in principle defines Z_{5} for any combination of arguments G, ν_{eff}, and β, it is in practice only valid in the subspace of the arguments that is well populated by the quasars. Most importantly, this does not include sources that are brighter than G ≃ 14, redder than ν_{eff} ≃ 1.48, or bluer than ν_{eff} ≃ 1.72. In order to extend Z_{5} in these directions, we resort to differential methods using physical pairs and sources in the LMC.
Fig. 5 Mean parallax of quasars binned by magnitude (panel a), effective wavenumber (panel b), and sine of ecliptic latitude (panel c). In each bin the dot is the mean parallax in EDR3 weighted by , with error bars indicating the estimated standard deviation of the weighted mean. The two red lines indicate the weighted mean value (− 21 μas) and median (− 17 μas) of the full sample. 
Coefficients for the function Z_{5}(G, ν_{eff}, β) fitted to the quasar parallaxes.
Fig. 6 Example of interactions among the dependences of quasar parallaxes on colour and ecliptic latitude. In a quadratic regression of parallax vs. sinβ, the linear coefficient (filled blue squares) exhibits a strong variation with effective wavenumber, while no such trend is shown by the quadratic coefficient (open red circles). The points have been slightly displaced sideways to avoid overlapping error bars. 
Fig. 7 Mean residual in quasar parallax after a regression on G, ν_{eff}, and β (see text), plotted against the confusion factor from Eq. (6). In each bin the dot is the mean residual weighted by , with error bars indicating the estimated standard deviation of the weighted mean. The broken dashed line is the dependence modelled by Eq. (7). 
Fig. 8 Celestial map of the mean contamination bias of the quasars, as given by Eq. (8) and Table 3. The map uses the same projection in ecliptic coordinates as the maps in Fig. 2. 
Coefficients for the extended function Z_{5}(G, ν_{eff}, β, X) fitted to the quasar parallaxes.
4.2 Large Magellanic Cloud
The distance modulus of the LMC, mag (de Grijs et al. 2014), corresponds to a parallax in the range 19.2 to 20.9 μas. The depth and inclination of the system mean that the parallaxes of individual stars have a true dispersion (and gradient) of about one μas. For our analysis we did not assume any specific mean distance to the LMC, only that the selected member sources have the same (but unknown) parallax, regardless of their colours and magnitudes. The LMC data can therefore be used to map the bias function Z_{5} (G, ν_{eff}, β) at the position of the LMC,β ≃−85°, up to an unknown additive constant.
The selection of sources in the LMC area for the analysis in this section is described in Appendix B. The sample consists of more than two million sources from EDR3, brighter than G = 19 and located within a 5° radius of the LMC centre. As discussed in the appendix, it is thought to be reasonably clean at least down to G ≃ 18. A colourmagnitude diagram (CMD) of the sample is shown in the left panel of Fig. 9.
In the right panel of Fig. 9 the same CMD is colourcoded by the median parallax at each point of the diagram. The predominantly greenish colour shows that the overall median parallax is close to zero, roughly consistent with the expected true parallax of 20 μas and a global parallax bias around − 20 μas. The dark red patch at G > 18, ν_{eff} < 1.4 is the only part of the diagram dominated by foreground stars. Several systematic deviations from the overall median are clearly visible, and cannot be attributed to foreground stars. These include the rather sharp divisions at G ≃ 13.0 and G ≃ 16.0, coinciding with the magnitude boundaries of window classes WC0b, WC1, and WC2 (Sect. 2.1 and Fig. 1); a strong difference (or gradient) in the bias for G < 13.0 between the mainsequence (at ν_{eff} ≳ 1.7) and the red supergiants (at ν_{eff} ≲ 1.4); an upturn of the bias for the bluest stars, which is stronger for G > 16.0 than for the brighter stars; and a downturn of the bias for the reddest stars, at least for G ≃ 14 to 16. For the faintest stars there appears to be a depression of the bias at intermediate colours (ν_{eff} ≃ 1.55). The cause of this depression is not known.
The systematic variation of parallax with colour and magnitude is further explored in Fig. 10, where each panel displays a different magnitude interval. The black dots show the median parallax binned by ν_{eff}. The three panels for G ≥ 13.1 show that the upturn of parallax values for the bluest stars sets in abruptly around ν_{eff} = 1.72 μm^{−1}, which is clearly related to the restriction in the range of wavenumbers to 1.24≤ ν_{eff} ≤ 1.72 μm^{−1} for the LSFand PSF calibrations, as discussed in Sect. 2.2. At the other (red) end of the interval no clear break is seen, although only a few stars in the LMC are redder than 1.24 μm^{−1}. Between 1.24 and 1.72 μm^{−1} the relation is approximately linear in the left (bluer) half of the interval, but it is clearly nonlinear in the right (redder) half, at least for G > 13.
These variations can be described by the general model in Appendix A, although several simplifications are needed for a welldetermined fit. Most importantly, because the LMC sample only covers a small area near the south ecliptic pole, the dependence of parallax on β cannot be determined; the resulting fit is valid at the mean position of the LMC, β ≃−85° or sin β ≃−0.996, but not necessarily at other locations. Modelled on Eqs. (A.5) and (A.6), but using only the basis functions in G (Eq. (A.2)) and ν_{eff} (Eq. (A.3)), the function fitted to the EDR3 parallaxes of the LMC sample is (9)
are piecewise linear functions in G, and y_{ij} are free parameters of the model, estimated by a weighted leastsquares procedure. Here and elsewhere in this paper, the leastsquares estimation is made robust against outliers by iteratively removing points that deviate by more than four times a robust estimate of the (weighted) RMS residual among all the data points. Owing to the small number of sources at the bright end and their limited coverage in ν_{eff}, the functions p_{j} (G) are forced to be constant in each of the magnitude intervals 6.0–10.8 and 11.2–11.8 (using the device described in footnote 2), and the parameters with j ≥ 2 are set to zero for i ≤ 6. At the faint end, the fit is limited to sources brighter than G = 18 in order to minimise contamination effects, and the basis functions g_{i} (G) for i = 11 and 12, which lack support for G < 18 (see Fig. A.1), are omitted in Eq. (10). The resulting fit is given in Table 4, where coefficients assumed to be zero are indicated by a dash (–). The fitted function, evaluated at representative magnitudes, is shown by the blue curves in Fig. 10.
In Table 4 the coefficients p_{0} at the different magnitudes give the mean parallax of the LMC sample reduced to the reference wavenumber ν_{eff} = 1.48 μm^{−1}. As they refer to the location of the LMC, it is not useful to compare them with the coefficients q_{00} from the quasars (Table 3), which are averages over the celestial sphere. Because we wish to avoid that our analysis depends in any way on an assumed distance modulus to the LMC, the fitted coefficients p_{0} were not further used in the determination of Z_{5}. The remaining coefficients p_{1} through p_{4} map the variation in parallax bias at the LMC location as a function of ν_{eff}. For WC0 (G ≤ 12.9) only the mean gradient in wavenumber (p_{2}) is determined and exhibits very significant variations with magnitude; the major breaks at G ≃ 11 and 13 are clearly visible in Fig. 9 (right). For WC1 and WC2 (G ≥ 13.1), the most striking feature is the relative constancy of p_{1}, p_{2}, and p_{3} with magnitude: the tabulated values all agree to within ±2σ with their weighted mean values, which are p_{1} = − 20.0 ± 1.2 μas μm, p_{2} = −1257 ± 58 μas μm^{3}, and p_{3} = +200 ± 39 μas μm. On the other hand, the coefficients p_{4}, describing the colour gradient at the blue end (ν_{eff} > 1.72 μm^{−1}), show a very clear progression with magnitude.
Fig. 9 Colourmagnitude diagrams of the LMC sample. Left: diagram colourcoded by the number of sources at a given position in the diagram. Right: diagram colourcoded by the median parallax in EDR3. 
Function ϖ_{EDR3}(G, ν_{eff}) fitted to the LMC sample.
Fig. 10 Median parallaxes for the LMC sample in Fig. 9 plotted against ν_{eff} in six magnitude intervals, as indicated in the diagrams. The orange points show the parallaxes of individual sources and give an impression of the coverage in ν_{eff} and scatter in parallax. The black dots, with error bars, show the median parallax and its uncertainty in bins of ν_{eff} with at least 20 sources per bin. The solid blue curves show the fitted combination of basis functions in Eq. (9) evaluated for a representative magnitude (G = 10.0, 11.5, 12.5, 15.0, 17.0, and 18.0) in each interval. 
4.3 Combined fit using quasars and LMC sources (G > 13)
The quasar sample contains only a few sources redder than ν_{eff} ≃ 1.44 μm^{−1} or bluer than 1.69 μm^{−1}, and therefore cannot be used to estimate q_{jk} for j = 2, 3, and 4. The LMC sample, on the other hand, gives a good determination of p_{j} for j = 1…4 in the magnitude range 13–18, but only for the specific location of the LMC. In this section we attempt to combine the two datasets in order to extend the model to the full range of colours for G > 13. This is not possible without certain additional assumptions, detailed below; however, the distance to the LMC is still left as a free parameter.
From Eq. (A.4), using sinβ_{LMC} ≃−0.996, it can be seen that the coefficients in a combined model must satisfy (11)
Only for j = 1 is a direct comparison possible; this is shown in Table 5 for G ≥ 13.1. Although the agreement between and is moderate, the differences are roughly compatible with the uncertainties. The main conclusion from the comparison is that the LMC data are vastly superior to the quasars for estimating the gradient of the bias with colour at the location of the LMC.
If a combined solution of the quasar and LMC data were attempted, in which all the parameters q_{jk} for G > 13 are freely adjusted, the result would be a model that fits both datasets very well, and is very well determined in the LMC area, but has very large uncertainties in most other locations. This would obviously not be very useful. To proceed, we need to constrain the model. We boldly make the following assumptions for WC1 and WC2 (G ≥ 13.1): (i) that thegradient with colour in the moderately blue region (ν_{eff} = 1.48–1.72 μm^{−1}) is fully described by the interaction term q_{11}, that is, q_{10} = q_{12} = 0; (ii) that the curvature with colour in the moderately red region (ν_{eff} = 1.24–1.48 μm^{−1}) is fully described by the noninteraction term q_{20}, that is, q_{21} = q_{22} = 0; (iii) that the colour gradient at the red end (ν_{eff} < 1.24 μm^{−1}) is fully described by q_{30}, that is, q_{31} = q_{32} = 0; and (iv) that the colour gradient at the blue end (ν_{eff} > 1.72 μm^{−1}) is fully described by q_{40}, that is, q_{41} = q_{42} = 0.
Assumption (i) is consistent with the observation in Sect. 4.1 that q_{12} is generally insignificant and that q_{11} is far more important than q_{10} for the overall chisquare. Assumption (ii) cannot be tested by means of the quasars, but is partially supported by the test in Appendix C, in which the red clump stars show similar curvatures at sin β = ± 0.86. Indirect support for (iii) is provided by the EDR3–DR2 differences ΔZ analysed in Sect. 3, for example the panels in Fig. 4 for G = 13–16 and 16–18, where the changes in gradient at ν_{eff} = 1.24 μm^{−1} are not distinctly different in the two hemispheres. As DR2 did not use the clamping of ν_{eff} at 1.24 and 1.72 μm^{−1} described in Sect. 2.2, it is likely that any abrupt changes in the gradient with colour seen in ΔZ at these wavenumbers are caused by the EDR3 data. The same argument may be advanced in support of (iv), although the evidence in Fig. 4 for breaks at ν_{eff} > 1.72 μm^{−1} in the northern hemisphere is very weak.
With these assumptions the coefficients q_{jk} are effectively determined by the LMC data through the conditions in Eq. (11), and it is possible to fit both datasets with a single model. To account for the offset between p_{0} and the parallax bias at the LMC location, one additional unknown must be introduced, representing the true mean parallax of the LMC; this is constrained to be independent of G. Results are given in Table 6, except for the fitted LMC parallax, which is + 22.11 ± 1.10 μas.
Comparison of the parallax gradient in colour, p_{1}, at the position of the LMC, as estimated from LMC data, quasars (QSO), and physical pairs (PP).
4.4 Physical pairs (G < 13)
After mapping Z_{5}(G, ν_{eff}, β) for G > 13 by means of the quasars and the LMC, we now turn to the brighter sources. In the LMC area the gross variation in parallax bias with colour was mapped in Sect. 4.2 (Table 4), but this local result must now be extended to the whole sphere. To this end, we used binaries (here called ‘physical pairs’), in which it can be assumed that the components have similar true parallaxes, although their magnitudes and colours maybe very different. Using the results from previous sections to anchor the parallax bias of the fainter component among the quasars, it is then possible to estimate the bias of the brighter component. Details of the method are given in Appendix D, which also describes the selection of data used in the analysis below.
4.4.1 Results for G > 10
The methodoutlined in Appendix D was applied to pairs where the magnitude of the bright component (G_{1}) is in the range 10–14 mag, and that of the faint component (G_{2}) is in the range max(G_{1}, 13.1)–G_{1} + 4 mag. The biascorrected parallax of the faint component was computed as ϖ_{2}− Z_{5}(G_{2}, ν_{eff2}, β), where ϖ_{2} and ν_{eff2} are the published EDR3 parallax and effective wavenumber of the faint component, and Z_{5} is defined by the coefficients in Table 6. The investigated magnitude interval overlaps with Table 6 for G = 13.1–14, which provides a consistency check and possibly improved estimates of the parallax bias in an interval that is poorly covered by the quasars.
The LMC data show that significant variations in bias as a function of (at least) G and ν_{eff} exist for sources brighter than G ≃ 13. Owing to the relatively small number of available physical pairs, the analysis cannot be very complex but should still include the main known or expected dependences. A severe limitation is that the scarcity of bright sources with ν_{eff} < 1.3 or > 1.6 μm^{−1} makes it practically impossible to determine q_{20}, q_{30}, and q_{40}. From the LMC data we have q_{20} ≃−1257 μas μm^{3} for G > 13.1, so that even though it cannot be usefully estimated from the physical pairs in the magnitude interval 13.1–14, consistency requires that the corresponding term is subtracted from the parallax of the brighter component before the remaining parameters are fitted. This procedure indeed reduces the sum minimised in Eq. (D.4), but not by a significant amount. A similar improvement is obtained by applying the same a priori correction for G < 13.1. We therefore assumed that q_{20} ≃−1257 μas μm^{3} throughout the range G < 14. Concerning q_{30} and q_{40} we assumed that they are zero for G < 13.1. Noting that a fit including the six coefficients q_{jk} (j ≤ 1, k ≤ 2) gives insignificant results for q_{12} at all magnitudes (lower right panel of Fig. 11), we also assumed q_{12} = 0.
Figure 11 shows results for the remaining five coefficients q_{00}, q_{01}, q_{02}, q_{10}, and q_{11}. The interaction of these terms with magnitude is fully mapped by independent fits in 35 bins of G_{1} covering the interval 10 < G_{1} < 14. All five coefficients show significant variations with G, which in most cases can be related to the boundaries discussed in Sect. 2.1.
The different symbols in Fig. 11 show the results from using different intervals in ρΔμ. The filled black circles are for ρΔμ < 2 arcsec mas yr^{−1}, which givesthe most precise estimates; the error bars are ± 1σ uncertainties obtained by bootstrap resampling. The coloured symbols (open circles, triangles, and squares) are for nonoverlapping intervals in ρΔμ. The results for the different nonoverlapping intervals in ρΔμ are to a high degree statistically independent, which gives an additional indication of the uncertainty. (They are not completely independent, as systems with more than two components may appear in more than one interval of ρΔμ.) The absence of any obvious trend in q_{00} with ρΔμ suggests that contamination bias is negligible.
More reliable estimates of the coefficients were obtained in a simultaneous fit of all the parameters, using Eq. (D.4). Here, q_{jk} were constrained to be piecewise linear functions of G with breakpoints at G = 10.8, 11.2, 11.8, 12.2, 12.9, and 13.1 (see Appendix A). For numerical stability, we required that q_{jk} are constant for G < 10.8 and > 13.1. The resulting fit is given in Table 7 and is shown by the dashed blue lines in Fig. 11.
The coefficients in the last row of Table 7 are in excellent agreement with the joint quasar and LMC results (Table 6) at G = 15.9, while the agreement is poorer at G = 13.1, where the coefficients in Table 6 are considerably more uncertain. In Table 8 we joined the two datasets by adopting the quasar/LMC results for G ≥ 15.9 and the results from the physical pairs for G ≤ 13.1, but taking q_{30} and q_{40} at G = 13.1 from Table 6 as they could not be determined from the pairs.
Coefficients for the extended function Z_{5}(G, ν_{eff}, β, X) as obtained in a combined fit to the quasar and LMC samples.
Coefficients of Z_{5}(G, ν_{eff}, β) as estimated from physical pairs for G > 10.
4.4.2 Extending the analysis to G = 6
The previous analysis of physical pairs did not reach brighter than G ≃ 10 owing to the restrictions G_{2} − G_{1} < 4 mag and G_{2} > 13.1, where the latter condition came from the necessity of using the bias function from Sect. 4.3 (valid for G > 13.1) to correct the parallax of the faint component. As a consequence, the number of available pairs rapidly decreased towards the bright end, not only because of the general scarcity of bright stars, but also because a decreasing fraction of them have faint components in the required range.
Using the coefficients in Table 7 to define a provisional parallax bias for G > 10, we can now extend the analysis to pairs as bright as G_{1} ≃ 6 by removing the second constraint, that is, by including all pairs with G_{1}< G_{2} < G_{1} + 4 mag. Not only does this extend the analysis to G > 6, but the results for G > 10 are also improved by the many more pairs included. For example, between G_{1} = 10.0 and 10.8, twice as many pairs become available for the analysis. This also allows us to remove the constraint that the coefficients are constant for G < 10. The resulting estimates are different from those in Table 7, with the most important differences seen towards the bright end. A repetition of the analysis using the improved coefficients for the biases of the faint components results in yet another slightly different set of coefficients. However, after a few more iterations, the coefficients were found to be completely stable, and hence internally consistent with the data for all useful pairs. The end result is shown in Table 9, which is our final estimate of the bias function Z_{5} (G, ν_{eff}, β) for sources withfiveparameter solutions in EDR3. Figure 21 is a visualisation of this function. The top panel of Fig. 20 shows values of the bias function for a representative selection of sources, taking into account their actual distribution in effective wavenumber at a given magnitude.
5 Sixparameter solutions
Because of their different treatments in the image parameter determination and astrometric solution, the five and sixparameter solutions have different systematics and it is necessary to consider their parallax biases separately. In principle, the same method as was used for the fiveparameter solutions could be applied to the sixparameter solutions, but in practice, this in not possible owing to the much smaller number of sixparameter solutions at all magnitudes, except for G ≳ 19 and at the very bright end (Fig. 12). Recalling (Sect. 2.2) that sixparameter solutions are used for sources that did not have reliable colour information from the BP and RP photometers in Gaia DR2, we may also note that their observationsare more often disturbed by other sources than the fiveparameter solutions, and therefore they are generally more problematic.
To circumvent the scarcity of suitable sixparameter solutions, we bootstrapped the estimation of Z_{6} in the following way on the already determined Z_{5}. For 8.16 million of the primary sources we have both five and sixparameter solutions (see Sect 2.2), and this sample was used here tomap the systematic differences in parallax between the two kinds of solution. Figure 13 shows the median of ϖ_{6} − ϖ_{5} versus G and ν_{eff}, where ϖ_{5} and ϖ_{6} are the parallaxes of a given source as obtained in the five and sixparameter solutions. The most prominent features in the plots are the positive gradient in ν_{eff} for 16 ≲ G ≲ 19.5 and 11 ≲ G ≲ 12, and an opposite gradient for the faintest stars. There are clear differences between the southern and northern ecliptic hemispheres. The systematics of ϖ_{6}− ϖ_{5} thus depend in a complex way on at least G, ν_{eff}, and β. Owing to the restrictions in the selection of primary sources for the astrometric solution, the colour region outside of the interval 1.24–1.72 μm^{−1} in ν_{eff} cannot be mapped with this sample.
In the Gaia archive, no effective wavenumber derived from photometry is provided for sources with sixparameter solutions, and many of them also lack a colour index such as G_{BP} − G_{RP}. The parallax bias function Z_{6} must therefore be expressed in terms of the pseudocolour instead of ν_{eff}. Figure 14 shows the same differences ϖ_{6} − ϖ_{5} as in Fig. 13, but plotted versus . The main trends are the same, only amplified for the faintest sources. Uncertainties in scatter some points outside of the interval 1.24–1.72 μm^{−1} especially at the faint end, where the pseudocolour becomes rather uncertain. In the southern hemisphere, the strong gradient versus at the faint end is produced by the predominantly negative correlation between parallax and pseudocolour seen in Fig. 2 (panel b).
To estimate the parallax bias function for sixparameter solutions, , the following procedure was used. First, a corrected fiveparameter parallax was calculated as (12)
for each of the ~8 million primary sources that have both kinds of solution, and for which both ν_{eff} and are available. This calculation used the coefficients in Table 9. Next, an estimate of the parallax bias in the sixparameter solutions was obtained as (13)
Finally, a robust weighted leastsquares fit of the general model in Eqs. (A.5) and (A.6) to ẑ_{6} was made, with replacing ν_{eff} in Eq. (A.3). In the fit, the data were weighted by the inverse variance, calculated as the sum of the formal variances of ϖ_{5} and ϖ_{6}. This overestimates the random errors (because the five and sixparameter solutions are positively correlated), but neglects the uncertainty of the correction Z_{5}.
Figure 15 shows the median values of ẑ_{6} in Eq. (13), plotted versus G and using the same divisions by ecliptic latitude as in Figs. 13 and 14. From Fig. 15 it is clear that Z_{6} cannot be usefully determined for m^{−1} from the available sample of common five and sixparameter solutions, at least not for G ≲ 18 (and even for fainter sources it would be very dubious, given the correlation mentioned above). For a similar reason, the results at the red end m^{−1} are also highly uncertain. Because the boundaries at 1.24 and 1.72 μm^{−1} have no special significance in the sixparameter solutions (except that the chromaticity calibration does not extend beyond these limits), it could well be that an extrapolation of the fit gives reasonable results at the more extreme colours. However, rather than relying on this, we chose to assume q_{3k} = q_{4k} = 0 in the fitted model; this effectively means that Z_{6} is clamped to its value at 1.24 or 1.72 μm^{−1} for more extreme pseudocolours. This has the added benefit of restricting Z_{6} to a finite interval even when the pseudocolour is grossly incorrect. In line with the analysis of Z_{5}, we also assumed q_{21} = q_{22} = 0.
The resulting coefficients for Z_{6} are given in Table 10. The function is visualised in Figs. 16 and 22. To facilitate comparison with the sample data in Fig. 15, the data shown in Fig. 16 were averaged over all latitudes in the middle panel, and over each hemisphere in the top and bottom panels. (To calculate the average over the whole celestial sphere, only the coefficients q_{j0} should be included; similarly, the averages in the two hemispheres are obtained from q_{j0} ± q_{j1}∕2.)
Fig. 11 Coefficients q_{jk} estimated from physical pairs as functions of G. Results for q_{12} are considered insignificant and set to zero when fitting the other five coefficients. The different symbols represent different selections onρΔμ: 0–0.5 arcsec mas yr^{−1} (red circles), 0.5–1 arcsec mas yr^{−1} (green triangles), 1–2 arcsec mas yr^{−1} (blue squares), and 0–2 arcsec mas yr^{−1} (filled black circles with lines and error bars). The dashed blue line is the global fit in Table 7. The vertical grey lines show the breakpoints for the basis functions in G defined by Eq. (A.2). 
Fig. 12 Fraction of sources in Gaia EDR3 with different kinds of solutions: fiveparameter solutions (P5 = solid red curve), sixparameter solutions (P6 = dashed blue curve), and twoparameter solutions (P2 = thin grey curve). The twoparameter solutions are ignored in this paper as they do not have parallaxes. 
Fig. 13 Median parallax difference between six and fiveparameter solutions as a function of ν_{eff} and G for the 8 million sources with both kinds of solutions. The panels show selections depending on ecliptic latitude β: southern hemisphere (panel a), all latitudes (panel b), and northern hemisphere (panel c). In this sample no sources lie outside of the interval 1.24 < ν_{eff} < 1.72 μm^{−1}. 
Fig. 14 Same data as in Fig. 13, but plotted against the astrometrically determined pseudocolour (). Uncertainties in the pseudocolour scatter some points outside of the interval 1.24–1.72 μm^{−1}. 
Fig. 15 Median ẑ_{6} as a functionof ν_{eff} and G for the sample in Figs. 13 and 14. The panels show selections depending on ecliptic latitude: southern hemisphere (panel a), all latitudes (panel b), and northern hemisphere (panel c). 
Fig. 16 Parallax bias Z_{6} according toTable 10. The panels show mean values for the southern hemisphere (panel a), all latitudes (panel b), and the northern hemisphere (panel c). 
6 Validation of the bias corrections
In this section the bias functions Z_{5} and Z_{6} defined in Tables 9 and 10 are applied to EDR3 parallaxes in order to confirm the validity of the corrections. To some extent, the same data are used as for deriving the corrections in Sects. 4.1 and 4.2, in which case the tests are not a true validation of the functions but rather a consistency check of the procedures used to derive them. Exceptions are the tests of the quasar sample with sixparameter solutions and the bright (G< 13) LMC sample, neither of which were used in the derivation, and the mean parallax of the LMC, which was a free parameter in the fit. Additional tests are described in Fabricius et al. (2021) and Gaia Collaboration (2021a).
Coefficients for the function .
6.1 Using quasars
We applied the bias corrections to the quasar sample in EDR3, expecting the mean corrected parallax to be close to zero regardless of magnitude, colour, etc. The EDR3 table agn_cross_id contains 1 215 942 sources with fiveparameter solutions and 398 231 with sixparameter solutions. In the fiveparameter case, this sample is nearly the same as was used in Sect. 4 to derive the faint part of Z_{5}. By contrast, the sixparameter sample was not previously used: as detailed in Sect. 5, Z_{6} was estimated differentially with respect to Z_{5}, using a sample of sources for which both kinds of solution were available.
Figure 17 shows the results for the fiveparameters solutions, divided according to magnitude, effective wavenumber, and ecliptic latitude. Mean values of the uncorrected parallaxes (ϖ) are shown as open black circles, and those of the corrected values (ϖ − Z_{5}) as filled blue circles. The yellow dots, showing individual uncorrected values, are mainly intended to give an impression of the distributions in G, , and sin β of the sources. Although the corrected parallaxes are not perfectly centred on zero, especially for the brightest and reddest sources, the overall improvement is fairly satisfactory.
Figure 18 shows the corresponding results for the sixparameters solutions. Although the corrected values (ϖ − Z_{6}) are clearly better than the uncorrected values, it appears that a slightly larger correction than Z_{6} might often be required. A peculiar effect noted in this sample is that the mean corrected parallax is a strong function of various goodnessoffit measures such as the renormalised unit weight error (RUWE) and excess source noise. For example, the mean corrected parallax is generally closer to zero for the subset of quasars with insignificant excess source noise (astrometric_excess_noise_sig < 2), as shown by the red squares in Fig. 18. This trend is not present in the fiveparameter sample, and we currently have no explanation for it.
6.2 Using the LMC sample
The LMC data in Table 4 for G < 13 were not used anywhere in the analysis that led to the bias estimates Z_{5} in Table 9. For example, the faint components of the physical pairs were anchored in the combined quasar–LMC solution, which was derived based only on sources with G > 13. We can therefore use the bright LMC data for a partial validation of Z_{5}. Of particular interest is the conspicuous difference in parallax bias between the blue upper main sequence and the red giants, shown in the right panel of Fig. 9 and in Fig. 10 for G = 11.2–11.8 and 12.2–12.9. Figure 19 is a CMD of the LMC sample, similar to Fig. 9, but colourcoded by the median value of ϖ − Z_{5} at each point, whereZ_{5} is defined by the coefficients in Table 9. The diagram clearly shows that Z_{5} provides a reasonable correction in most parts of the CMD, including the bright part, although there is a suggestion that the bias is more negative (by ≃ 10 μas) than as given by Z_{5} for the blue branch at G < 13. This could indicate that the sizes of q_{10} and/or q_{11}, as obtained from the physical pairs, are underestimated in Table 9. As remarked in Sect. 4.2, the dark red patch in the lower right part of the diagram (roughly G > 18, ν_{eff} < 1.4) is caused by foreground stars that dominate this part of the CMD. The unmodelled depression of the bias for the faintest stars at intermediate colours, mentioned in Sect. 4.2, is visible as a greenish patch in Fig. 19.
In Sect. 4.3 the distance the LMC was a free parameter in a fit of Z_{5} to the combined quasar and LMC data, yielding a mean parallax of + 22.11 ± 1.10 μas. This is about two standard deviations higher than commonly accepted values, for example, + 20.17 ± 0.25 μas (Pietrzyński et al. 2019; including systematic uncertainty). In absolute measure, however, the difference of about 2 μas is small compared with the regional variations of the quasar parallaxes shown for example in Fig. 2c. The angular power spectrum of quasar parallaxes in Gaia EDR3 is discussed in Lindegren et al. (2021), who estimate that the RMS variation of the parallax systematics (excluding the global offset) is about 10 μas on angular scales ≳10°. The uncertainty of ±1.1 μas from the fit in Sect. 4.3 does not take these variations into account because the LMC only probes a smaller area. We can therefore conclude that the mean corrected parallax of the LMC agrees remarkably well with the accepted value.
7 Way forward
The bias functions Z_{5} and Z_{6} provide a recipe for the systematic correction of the EDR3 parallaxes based on the particular choices of data, method, and bias models described in the preceding sections. These choices contain considerable elements of uncertainty and arbitrariness, and are no doubt also coloured by our preconceived notions. The results should therefore not be regarded as definitive. On the contrary, it is vitally important that alternative routes are explored towards better compensating for the systematics in Gaia data. Gaia Data Release 3, which is expected to be released in 2022, will be a superset of EDR3, so that investigations made based on EDR3 parallaxes will not be superseded until Data Release 4 (DR4) much later. Although DR4 is expected to be significantly better in terms of systematics, it will not be unbiased. Methodological developments made using the current data will therefore remain applicable for a long time. We point out a few possible directions for this work.
Global analysis: the approach taken here is highly heuristic, where models are built up gradually in interaction with the data analysis, using a range of different analysis tools. This is often a good way to arrive at a reasonable approximation quickly, and it naturally incorporates the already existing knowledge of the structures and difficulties inherent in the data. It probably does not give the optimal result, however. When a general model has been established, it would be advantageous to make a singe global fit to all the data. Using standard statistical techniques (e.g. Burnham & Anderson 2002), we can select the appropriate submodel by objective criteria and can evaluate the confidence intervals. This should give more precise estimates by combining datasets optimally and avoiding the accumulation of errors in the current multistep fits. It might also reveal new dependences and interactions.
More and different data: increasing the amount of data included in a global analysis could improve the precision of the bias estimation and perhaps extend its validity in magnitude–colour space. Examples of datasets that should be explored are the stars in open and globular clusters, which have the potential to cover a wide range in magnitudes and different environments, for instance crowded areas. They could also provide a more direct link between Z_{5} (for the brighter stars) and Z_{6} (for the faint members) in each cluster.
Datadriven modelling: the access to highquality astrometric, photometric, and radialvelocity data for very large stellar samples makes it possible to construct and fit statistical models that do not depend on the physical models of specific types of stars (e.g. Schönrich et al. 2019; Hogg et al. 2019). These methods depend on having a very large number of objects to beat down statistical uncertainties, and will therefore not provide a high resolution in several variables, but rather an independent overall validity check.
Bias modelling for specific applications: general bias models of the kind developed here are not necessarily the best way to handle parallax bias in specific applications. It might for example be better to include it as a free parameter directly in the physical model, for instance for luminosity calibrations. This approach was taken by many researchers using Gaia DR2 data (see Sect. 1 for some examples), and it remains a valid alternative to correcting the data.
Systematics in proper motions: A priori there is no reason to expect that the proper motions in EDR3 are less affected by systematics than the parallaxes. After all, they are jointly determined from the same observations. Fabricius et al. (2021) showed an example (their Fig. 25) in which a discontinuity of about 25 μas yr^{−1} is detected at G ≃ 13 in the EDR3 proper motions μ_{α*} of cluster stars. It is therefore clearly interesting to extend our analysis to the components of proper motion. Although less crucial than the parallax bias in most astrophysical and galactic astronomy applications of Gaia data, proper motion biases are relevant for a number of the most exacting applications, such as the search for exoplanets. Mapping these biases in clusters, however, may be nontrivial in view of internal motions, which may include systematic patterns from expansion, rotation, etc.
Feedback to Gaia calibration models: Several distinct features in Z_{5} and Z_{6}, such as the abrupt changes and reversal of gradients in colour at G ≃ 11 and 13 (Fig. 20), can be traced back to specific elements of the instrument calibration in the dataprocessing chain for EDR3. This can be used in a kind of reversed engineering to understand how the calibration models need to be improved in order to decrease systematics in future releases. This is part of the normal cyclic development work in the Gaia dataprocessing consortium, and is much helped by having tools at our disposal to evaluate systematics in the astrometric solution.
Fig. 17 Parallaxes for 1.2 million quasars with fiveparameter solutions in Gaia EDR3. Yellow dots show theindividual values plotted vs. magnitude (panel a), effective wavenumber (panel b), and sine of ecliptic latitude (panel c). Open black circles show mean values of the uncorrected parallaxes (ϖ) in bins of magnitude etc., and filled blue circles show mean values of the corrected parallaxes (ϖ − Z_{5}) in the same bins. Mean values are calculated using weights . Error bars indicate the estimated standard deviation of the weighted mean in each bin. 
Fig. 18 Parallaxes for 0.4 million quasars with sixparameter solutions in Gaia EDR3. Yellow dots show the individual values plotted vs. magnitude (panel a), pseudocolour (panel b), and sine of ecliptic latitude (panel c). Open black circles show mean values of the uncorrected parallaxes (ϖ) in bins of magnitude etc., and filled blue circles show mean values of the corrected parallaxes (ϖ − Z_{6}) in the same bins. Red open squares show mean values of the corrected parallaxes for the 78% of the quasars that have insignificant excess source noise, using a slightly different binning. Mean values are calculated using weights . Error bars indicate the estimated standard deviation of the weighted mean in each bin. 
Fig. 19 Colourmagnitude diagram of the LMC sample colourcoded by the median of the EDR3 parallax after subtracting Z_{5} as given by the cofficients in Table 9. 
Fig. 20 Parallax bias Z_{5} (top) and Z_{6} (bottom) computed according to Tables 9 and 10 for a representative sample of sources with five and sixparameter solutions in EDR3. The dots show values for the complete sample of sources (for G < 11.5) or for a random selection (G > 11.5). The colour scale indicates the mean ν_{eff} or at a given point, thus giving an impression of the mean colour dependence of the bias. The black curves show the 1st, 10th, 50th, 90th, and 99th percentiles. The thick curve is the 50th percentile or median. 
8 Summary and conclusions
We have investigated the parallax bias in Gaia EDR3 and the variation of this bias with magnitude, colour, and ecliptic latitude. The direct estimation of the bias using quasars was complemented by indirect methods using physical binaries and stars in the LMC. The indirect methods were strictly differential with respect to the quasars; in particular, no specific value was assumed for the distanceto the LMC. The expected functional form of the dependences on magnitude, colour, and ecliptic latitude was derived partly from a consideration of the known properties of the instrument and data processing, and partly from a mapping of the systematic difference between parallaxes in Gaia EDR3 and DR2 (Table 1).
Complex dependences on all three variables are evident in all the data, and the dependences are not the same for sources that have five and sixparameter solutions. For sources with fiveparameter solutions in Gaia EDR3 (astrometric_params_solved = 31) the fitted bias function, Z_{5}(G, ν_{eff}, β), is given by Eqs. (A.3)–(A.5) using the functions q_{jk} (G) obtained by linear interpolation in Table 9. Here, G is the mean magnitude of the source in the broadband Gaia photometric system (phot_g_mean_mag), ν_{eff} is the effective wavenumber (nu_eff_used_in_astrometry), and β is the ecliptic latitude (ecl_lat). A dash (−) in the table should be interpreted as zero. The function Z_{5} is illustrated in Fig. 21 and in the top panel of Fig. 20. Nominally, the bias ranges from − 94 to + 36 μas, with an RMS scatter of 18 μas when the full range of colours and magnitudes is considered; when it is weighted by the actual distribution of colours per magnitude in EDR3, as in Fig. 20 (top), the RMS scatter is about 13 μas.
For sources with sixparameter solutions in Gaia EDR3 (astrometric_params_solved = 95), the parallax bias function, , is similarly given by the functions q_{jk} (G) obtained by interpolation in Table 10. Here, is the astrometrically estimated effective wavenumber, known as pseudocolour (pseudocolour). The function Z_{6} is illustrated in Fig. 22 and in the bottom panel of Fig. 20. It is very uncertain when is < 1.24 μm^{−1} or > 1.72 μm^{−1}. The bias ranges from − 151 μas to + 130 μas, with an RMS scatter of 21 μas for the full range of colours, or 15 μas if weighted by the actual distribution of colours.
The derived relations are only applicable to the parallaxes in Gaia EDR3. Regarded as a systematic correction to the parallax, the bias function Z_{5} or Z_{6} should be subtracted from the value (parallax) given in the archive. Python implementations of both functions are available in the Gaia web pages^{3}.
As recipes for the systematic correction of the EDR3 parallaxes, these functions should be regarded as provisional and indicative. While we have reason to believe that their application will in general reduce systematics in the parallaxes, this may not always be the case. Users are urged to make their own judgement concerning the relevance of the indicated bias correction for their specific applications. Whenever possible, depending on the type and number of sources under consideration, users of EDR3 should try to derive more targeted bias estimates for their specific use cases.
It is difficult to quantify uncertainties in Z_{5} and Z_{6}. In the region of the parameter space that is well populated by the quasars (essentially G ≳ 16 and 1.4 ≲ ν_{eff} ≲ 1.7 μm^{−1}), they may be as small as a few μas, but beyond that region, uncertainties are bound to be greater because of the indirect methods we used. For redder sources (ν_{eff} ≲ 1.4 μm^{−1}, corresponding to G_{BP}−G_{RP} ≳ 1.6), both Z_{5} and Z_{6} depend critically on assumption (ii) in Sect. 4.3 that the curvature in colour, seen in the LMC, is the same over the whole sky. If this turns out to be incorrect, it could mean that the bias function for very red sources in the northern hemisphere is off by more than 40 μas. As discussed in Appendix C, we do not think this is likely, but the possibility should be kept in mind. Beyond the nonclamped interval 1.72 > ν_{eff} > 1.24 μm^{−1} (corresponding to 0.15 ≲ G_{BP} − G_{RP} ≲ 3.0) all results are in any case quite uncertain owing to the way colour information is handled in the LSF and PSF calibrations and astrometric solutions for EDR3 (Sect. 2.2).
Two unmodelled features mentioned in the paper could merit further investigation. One is the depression of the bias for G ≳ 18, ν_{eff} ≃ 1.55 μm^{−1} that is visible as a greenish patch in Fig. 19. As this region of colour–magnitude space is well covered by the quasars, the feature is probably not generally present but might be particularly strong in the LMC area. The other unmodelled feature is the dependence of the parallax bias for sixparameter solutions on excess source noise illustrated in Fig. 18. It is possible that both features are caused by to a hitherto unexplored tendency at faint magnitudes for the bias to become more negative in crowded areas, where the source excess noise or RUWE also tends to be higher.
As discussed in connection with Fig. 2 (panel c) and more extensively elsewhere (Lindegren et al. 2021), additional systematics in the EDR3 parallaxes have been identified, for example depending on position on small and intermediate angular scales. These cannot easily be mapped to any useful precision, and should rather be modelled as correlated random errors. The angularpower spectrum of quasar parallaxes presented in Sect. 5.6 of Lindegren et al. (2021) may be used to that end. Depending on the angular scale considered, the estimated RMS variation with position ranges from 5 to 26 μas, that is, of the same order of magnitude as the RMS variations in Z_{5} and Z_{6} as functions of colour or magnitude.
While itis easy enough to demonstrate that the EDR3 parallaxes contain significant systematics, it is extremely difficult to obtain accurate estimates of the bias beyond the limited region of parameter space that is well populated by the quasars. This paper does not claim to give a definitive answer but merely a rough characterisation of what we have found to be the main dependences. It is likely that better and possibly quite different estimates can be obtained in the future by means of more refined and comprehensive analysis methods. Continued exploration of the systematics is important not least in order to gain a better understanding of their causes. In the end, this will hopefully lead to much lower levels of systematics in future Gaia data releases.
Fig. 21 Parallax bias Z_{5}(G, ν_{eff}, β) according to Table 9. The panels show cuts at (a) ecliptic latitude β = − 90°, (b) β = 0°, and (c) β = +90°. 
Fig. 22 Parallax bias according to Table 10. The panels show cuts at (a) ecliptic latitude β = − 90°, (b) β = 0°, and (c) β = +90°. 
Acknowledgements
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. This work was financially supported by the Swedish National Space Agency (SNSA/Rymdstyrelsen), the European Space Agency (ESA) in the framework of the Gaia project, and the German Aerospace Agency (Deutsches Zentrum für Luft und Raumfahrt e.V., DLR) through grants 50QG0501, 50QG0601, 50QG0901, 50QG1401 and 50QG1402. We thank the Centre for Information Services and High Performance Computing (ZIH) at the Technische Universität (TU) Dresden for generous allocations of computer time; C. Babusiaux, A.G.A. Brown, J. Maíz Apellániz, A. Vallenari, and the anonymous referee for many helpful suggestions and corrections to the manuscript; and P. Ramos for implementing the Python version of the bias functions. Diagrams were produced using the astronomyoriented data handling and visualisation software TOPCAT (Taylor 2005).
Appendix A Parametrised functions
The various absolute or differential bias functions discussed in the paper are written as linear combinations of a finite set of threedimensional basis functions, with G, ν_{eff} (or ), and β as independentvariables. To allow interactions among the variables, the set of threedimensional basis functions should in the most general case be the outer product of the onedimensional basis functions on each axis. The generic function is therefore (A.1)
where g_{i}(G) are the basis functions in magnitude, c_{j}(ν_{eff}) the basis functions in colour, and b_{k}(β) the basis functions in ecliptic latitude. The coefficients z_{ijk} are the free parameters used to fit Z to the given data.
The magnitude dependence is modelled as a continuous piecewise linear function with breakpoints (knots) at γ_{0…12} = 6.0, 10.8, 11.2, 11.8, 12.2, 12.9, 13.1, 15.9, 16.1, 17.5, 19.0, 20.0, and 21.0 mag. The corresponding basis functions are (A.2)
(Fig. A.1). A linear combination of these functions may provide a reasonable approximation of variations such as those shown in the top panel of Fig. 3. In particular, the knots at G = 11.0 ± 0.2, 12.0 ± 0.2, 13.0 ± 0.1, and 16.0 ± 0.1 correspond to the transitions suggested in Fig. 1. An important property of this basis set is that at every breakpoint, exactly one basis function is = 1 and the rest are = 0. This means that for arbitrary coefficients q_{i}, the function q(G) = ∑_{j}q_{i}g_{i}(G) can be evaluated by linear interpolation among the coefficients, since q(γ_{i}) = q_{i}. This property is useful in connection with the alternative form in Eq. (A.5).
The dependence on colour is also modelled as a continuous piecewise polynomial, but of a more specific form inspired by plots like Figs. 4 and 10. Interior breakpoints are placed at ν_{eff} = 1.24, 1.48, and 1.72 μm^{−1}. The segments below 1.24, between 1.48 and 1.72, and above 1.72 are linear, while between 1.24 and 1.48 it is cubic with continuous first and second derivatives at 1.48 μm^{−1}. The basis functions are (A.3)
(Fig. A.2). Their coefficients thus represent the value at ν_{eff} = 1.48 μm^{−1}, the linear gradient between 1.24 and 1.72 μm^{−1}, the added cubic trend in the moderately red interval (from 1.24 to 1.48 μm^{−1}), and the linear gradients at the extreme colours (below 1.24 and above 1.72 μm^{−1}).
Finally, the basis functions in ecliptic latitude, (A.4)
(Fig. A.3) describe an arbitrary quadratic dependence on sin β. The term in b_{2} makes the three functions orthogonal for a uniform distribution of sources on the celestial sphere (which is also uniform in sin β).
An equivalent form of Eq. (A.1) is (A.5)
are piecewise linear in G. Equation (A.5) is useful because the functions q_{jk}(G) can be evaluated by linear interpolation among the coefficients z_{ijk}, which allows Z(G, ν_{eff}, β) to be given in the compact tabular form that we extensively use in this paper.
The coefficients z_{ijk} may be determined by standard curvefitting techniques. The problem is simple in the sense that it is linear in all the coefficients, but in practice, it is complicated by the presence of outliers and the often very incomplete coverage of the threedimensional space (G, ν_{eff}, β). Robust techniques such as L_{1}norm minimisation can be used to cope with outliers. Data coverage is more problematic and may require some judicious modification of the set of basis functions. A simple remedy could be to remove basis functions without support, for example c_{3} and c_{4} if there are too few sources of extreme colours; this is equivalent to setting the corresponding coefficients (z_{i3k} and z_{i4k}) equal to 0. The dependence on G can be simplified by adding constraints to the fit; for example, the constraint z_{i,jk} = z_{i+1,jk} will force the function q_{jk}(G) to be constant for γ_{i} ≤ G ≤ γ_{i+1}.
Fig. A.2 Basis functions c_{j}(ν_{eff}), j = 1… 4 according to Eq. (A.3). The (constant) function c_{0} = 1 is not shown. For better visibility the functions are vertically displaced by j∕2, and the amplitude of c_{2}(ν_{eff}) is increased by a factor 10. 
Appendix B Construction of the LMC sample
This appendix describes the selection of sources in the LMC area that we used for the analysis in Sect. 4.2. We also present some tests of the purity of the sample.
Adopting the LMC centre (α_{C}, δ_{C}) = (78.77°, −69.01°) as in Gaia Collaboration (2018c), we extracted all sources in Gaia EDR3 within a radius of 5° with fiveparameter solutions, G < 19, and ruwe < 1.4. A colourmagnitude (Hess) diagram of the resulting sample is shown in Fig. B.1a. The effective wavenumber ν_{eff} was used instead of a colour index on the horizontal axis, and the direction was reversed so that bluer stars are to the left and redder stars to the right. Several prominent features in this diagram are produced by Galactic foreground stars, and some additional filtering is clearly required. Because the purpose here is to study biases in parallax, it is essential that no filtering uses the actual parallax values, while the filtering already done based on position and ruwe cannot introduce a selection bias on the parallaxes. For the further selection we used the residuals in proper motion relative to a fitted, very simple kinematic model.
The positions and proper motions (including uncertainties and correlations of the latter) were converted into Galactic coordinates and then into rectangular orthographic components (x, y, μ_{x}, μ_{y}) using Eq. (2) in Gaia Collaboration (2018c), but replacing α, δ, μ_{α*}, and μ_{δ} by l, b, μ_{l*}, μ_{b}, and (α_{C}, δ_{C}) everywhere by (l_{C}, b_{C}) ≃ (279.77°, −33.77°). Using a robust (L_{1}norm minimisation) algorithm, we obtained the following linear relation in a fit including only stars brighter than G = 18: (B.1)
For each star in the sample, deviations in (μ_{x}, μ_{y}) from this model were transformed back into propermotion residuals (Δμ_{l*}, Δμ_{b}), from which the normalised deviations (B.2)
were computed. Figure B.1b shows the joint distribution of parallax and Δμ_{N} for the sample in Fig. B.1a. The selection Δμ_{N} < 10 is likely to produce a relatively clean sample of LMC stars. The colourmagnitude diagram of the filtered sample is shown in Fig. B.1c. This sample, which we used for the analysis in Sect. 4.2, contains 1457 sources with G < 13.0, 88 285 with G < 16.0, 519 203 with G < 17.5, and 2 371 761 with G < 19.0.
To validate and further quantify the cleanliness of the resulting sample, we performed the analogous selection and filtering of sources in two adjacent areas of the same size, but centred on (l_{C} ± 10°, b_{C}). Their positions and proper motions were transformed into x, y, μ_{x}, μ_{y} relative to the centre of the respective offset area, while residuals and Δμ_{N} were computed relative to the fixed values in Eq. (B.1). This should give approximately the same number andkinematic selection of Galactic stars in the three areas because their latitudes are equal and their spread in longitude is small. (Here it is essential that Galactic coordinates are used, so that the orientations of the x and yaxes relative to the Galactic plane are the same in the different areas.) The right panels d–f of Fig. B.1 show the results for one of the offset areas; plots for the other offset area are not shown but are qualitatively similar. The ϖ–Δμ_{N} plots in the middle panels confirm that the number and distribution of Galactic stars (the structure extending towards positive ϖ and high Δμ_{N}) are quite similar in the LMC and offset areas. Further statistical analysis shows that the sample in Fig. B.1c has negligible contamination by foreground stars down to G ≃ 18 at all colours, and down to G = 19 for the bluer sources (ν_{eff} ≳ 1.6).
Fig. B.1 Illustrating the selection and validation of the LMC sample. (a) Colourmagnitude diagram for the original sample centred on (l_{C}, b_{C}). (b) Joint distribution for the original sample of parallax and Δμ_{N}, the normalised propermotion difference to the fitted model. The vertical line marks the cutoff used for the selection based on proper motions. (c) Colourmagnitude diagram for the subsample with Δμ_{N} < 10. (d–f) Same as (a–c), but for a sample centred on (l_{C} − 10°, b_{C}), containing far fewer LMC stars, but roughly the same number of Galactic foreground stars as in (a–c). 
Appendix C Test using red clump stars
A crucial assumption for the joint quasar and LMC model in Sect. 4.3 is that the curvature in colour, represented by the basis function c_{2}(ν_{eff}) in Eq. (A.3), is the same over the whole celestial sphere; in other words, that interaction terms between ν_{eff} and β are negligible. This assumption is needed because too few quasars are red enough (ν_{eff} ≲ 1.35) to reliably determine the curvature at any point; instead it is derived entirely from the LMC sample, as illustrated by the curved segments in the three panels of Fig. 10 for G > 13. Such interactions, if they existed, would be represented by nonzero coefficients q_{21} and q_{22} in Table 6. The assumed invariance of the colour curvature with position rests mainly on the following analysis of the parallaxes of red clump (RC) stars in EDR3.
In the observational Hertzsprung–Russell diagram (HRD), the RC is a prominent feature made up of lowmass stars in the core heliumburning stage of their evolution. Red clump stars are widely used as standard candles because their scatterin absolute magnitudes is relatively small, especially at nearinfrared and infrared wavelengths (for a general overview, see Girardi 2016). Differences in the absolute magnitudes of RC stars are nevertheless expected, depending on factors such as their ages, effective temperatures, metallicities ([Fe/H]), and alphaelement abundances ([α/Fe]).
Fig. C.1 Hertzsprung–Russell diagrams for the union of samples N and S. Top: absolute magnitudes in the G band. Bottom: absolute magnitudes from W in Eq. (C.1). 
In the HRD of nearby giants (e.g. Fig. 10 in Gaia Collaboration 2018a), for which extinction is negligible, the RC stars occupy a narrow range of colours, approximately G_{RP} − G_{BP} = 1.2 ± 0.1, corresponding to ν_{eff} = (1.48 ± 0.02) μm^{−1}. Intrinsically,therefore, the RC stars are not nearly red enough to determine the curvature of the parallax bias versus colour, which only becomes significant for ν_{eff} ≲ 1.3 μm^{−1}. In the Galactic plane and bulge, however, it is possible to find many RC stars that are sufficiently reddened by interstellar extinction.
Fig. C.2 Differential parallax bias estimated by means of RC stars. The plots show differences between the measured parallaxes ϖ_{EDR3} and the photometric parallaxes ϖ_{RC} from Eq. (C.2) in area N (a) and S (b). Contours of constant density are shown as thin grey lines. The thick black curve traces the ridge of the feature created by the red clump stars. 
Fig. C.3 Differential parallax bias estimated by means of RC stars. The points show the locations of the ridges in Fig. C.2 for area N (filled red circles) and S (open blue squares). The curves show the expected locations at G = 17 according toTable 6 for area N (solid red) and S (longdashed blue), assuming that the curvature with colour (j = 2) is the same in both areas. If the curvature had flipped sign with β, the ridge location for area N should instead follow the upbending, shortdashed black curve. 
Because the absolute magnitudes of the RC stars are uncertain and may depend on many unknown factors, as mentioned above, our analysis of the RC parallaxes must be strictly differential and only compare samples with similar population characteristics. At the same time, we need locations at widely different ecliptic latitudes in order to see a possible variation in curvature with β. These conditions rule out the use of the inner part of the Galaxy, which only covers a limited range of ecliptic latitudes. A better strategy is to compare two areas in the Galactic plane that are symmetrically placed on either side of the Galactic centre and therefore probe similar ranges of galactocentric distances. The difference in sin β between the areas is maximised for Galactic longitudes l = ±90°, and we therefore chose areas around (l, b) = (90°, 0) or (α, δ) = (318.00°, +48.33°) (hereafter called ‘N’), and (270°, 0) or (α, δ) = (138.00°, −48.33°) (‘S’). Within a 5° radius of these points, we extracted all sources with fiveparameter solutions satisfying 13 < G < 17.5 and RUWE < 1.4 from EDR3. This gave 1.022 million sources in N and 0.686 million in S. The reason for this strong asymmetry seems to be the presence of a nearby (< 1 kpc) complex of dust clouds in the middle of the S area, possibly part of the Vela Molecular Ridge (Murphy & May 1991). While the high extinction produced by the dust clouds is desirable for our purpose, their proximity reduces the number of stars in the sample and increases their mean parallax, which is unfavourable for a precise determination of the bias. We therefore added two areas of 5° radius, centred on (l, b) = (85°, 0) and (275°, 0) to the previous selection. The resulting union sets (1.361 million sources in N; 1.561 million in S) are more similar in terms of the overall distribution of colours and distances, and still approximately symmetric in Galactic longitudes. The mean value of sin β is + 0.852 for the sources in N, and − 0.867 in S.
Figure C.1 is an HRD for the union of the two sets N and S. The absolute magnitude is (simplistically) computed using 1∕ϖ for the distance and ignoring extinction. In the top panel a the RC stars are seen as the concentration of points along a diagonal line about five magnitudes above the main sequence. In effective wavenumber, the feature starts at ν_{eff} ≃ 1.43 μm^{−1}, and it extends at least to ≃ 1.24 μm^{−1} because of the extinction reddening. In the bottom panel b the reddeningfree Wesenheit magnitude (C.1)
is used instead of G to compute the values on the vertical axis; here λ ≃ 1.9 is the slope of the reddening line for the photometric bands of Gaia (Ripepi et al. 2019). With this transformation the RC stars have an absolute magnitude M_{W} ≃−1.7 that is nearly independent of the colour. By lucky coincidence, the RC for unreddened nearby giants has a similar slope in the Gaia HRD, so the use ofW instead of G not only eliminates the effect of the reddening, but reduces the intrinsic scatter of the absolute magnitudes. Assuming that the RC stars have a fixed and known M_{W}, their photometric parallaxes can be computed as (C.2)
In Fig. C.2 we have plotted the differences ϖ_{EDR3} − ϖ_{RC}(G, G_{BP} − G_{RP}) for the sources in areas N and S versus ν_{eff}, using λ = 1.9 and M_{W} = −1.68 mag. If all the sources had absolute magnitude M_{W}, the points in these diagrams would outline the parallax bias as a function of ν_{eff}. Most of the sources are nearby mainsequence stars with more positive parallaxes; they are shown in the upper left corners of the diagrams. The RC stars form a concentration of points between 0 and − 100 μas on the vertical axis. The ridge locations, shown as the black curves in Fig. C.2, were estimated by binning the points in colour, using a bin width of 0.01 μm^{−1}, and finding the maximum of the distribution for a Gaussian kernel of 10 μas standard deviation. In Fig. C.3 the ridge locations for the two areas are plotted in the same diagram for easier comparison.
The ridge locations in Figs. C.2 and C.3 are very sensitive to the assumed values of λ and M_{W}. The value λ = 1.90 used here was adopted from Ripepi et al. (2019), while M_{W} = −1.68 was selected to give approximately the expected bias at ν_{eff} ≃ 1.4 μm^{−1} according to the analysis in Sect. 4.3. The red and blue curves in Fig. C.3 show the expected variation of the bias according to Table 6, evaluated at G = 17, where mostof the sources in the RC area are found. By adjusting M_{W}, it is possible to obtain agreement at a specific ν_{eff} for any reasonable choice of λ, but the slope and curvature of the relations defined by the points will be different. If λ is a function of the total extinction, this could also change the curvature of the relation. Little weight should therefore begiven to the circumstance that the points in Fig. C.3 roughly follow the curves computed using the values in Table 6. However, it is significant and supports the assumption made in Sect. 4.3 that the empirical relation is similar in the two areas for any reasonable M_{W} and reddening law. If, for example, the curvature instead of being independent of β were proportionalto sinβ, the expected relation in the N area would rather follow the upbending black curve in the diagram, which is clearly contradicted by the RC data.
For ν_{eff} ≲ 1.26 μm^{−1} the ridges and points in Fig. C.2 reach much more negative values on the vertical axis than expected from Table 6. Although deviations from the simplistic reddening law (constant λ) could contribute to this trend, itis probably mainly a selection effect, similar to the Malmquist bias. The strong downturn sets in for sources redder than approximately 1.26 μm^{−1}, at which point the total extinction in V is at least 4 mag and most of the sources defining the ridge are close to the faint magnitude limit of our sample at G = 17.5. A preferential selection of RC stars that are intrinsically 0.2–0.3 mag brighter than the mean population is enough to explain the discrepancy at ν_{eff} = 1.24 μm^{−1}. For ν_{eff} > 1.30 μm^{−1} this selection bias should be much smaller because the sources are on average at least one magnitude brighter. The conclusion from the analysis is that we see no evidence in the RC data for a difference in the curvature of parallax bias versus ν_{eff} between the northern and southern hemispheres.
Appendix D Data and method for physical pairs
This appendix describes the selection of data used for the analysis of physical pairs in Sect. 4.4. We also present the method applied to these data to estimate the parallax bias.
In this method only pairs with apparent separation ρ < 10 arcsec (5 × 10^{−5} rad) are considered. When their parallax is ≲20 mas, as is the case for 99.9% of the bright stars, it is not likely that the true parallax difference in the pair exceeds ± 1 μas. In the following we use subscripts 1 and 2 for the bright and faint components in a pair. When G_{2}> 13.1, the recipe derived from the quasars can be used to correct the EDR3 parallax of the faint component, ϖ_{2}, thus providing an (approximately) unbiased estimate of the true parallax of the pair. Considering a sample of physical pairs with similar G_{1}, ν_{eff1}, and β, we can then estimate the parallax bias as (D.1)
where med[x] is the sample median, which is used for robustness. A limitation of the method is that the components in a pair always have practically the same β (which is why β is not subscripted in Eq. (D.1)), so that the mapping of Z_{5} with respect to this parameter for the bright components must rely on the presumed known dependence on β for the faint components.
A main concern with the method is contamination by optical pairs, that is, chance alignments of physically unconnected sources with different true parallaxes. The use of the median in Eq. (D.1) provides good protection against outliers but is still biased if the outliers tend to fall mainly on one side of the median. This is the case for optical pairs, where the fainter component is likely to be more distant than the brighter, that is, ϖ_{1} − ϖ_{2} > 0, leading to a positive contamination bias in Eq. (D.1). Rather than eliminating the contamination bias by using a very clean sample, which may then be too small for our purpose, we adopt a heuristic approach, where a small amount of contamination is accepted and the calculated median is corrected by a statistical procedure. This requires accurate knowledge of the selection function of the sample, which in practice precludes using predefined catalogues such as the Washington Double Star Catalog (Mason et al. 2001).
Fig. D.1 Illustration of the selection of physical pairs and the contamination by optical pairs. Top panels: proper motion difference (Δμ) vs. separation (ρ) for source pairs in sample A (a) and B1 (b) before making the cut ρ > 2.2 arcsec indicated by the dashed vertical line. The three solid lines correspond to ρΔμ = 0.5, 1, and 2 arcsec mas yr^{−1}. The points are colourcoded by the median parallax difference between the components. In panel b the separation of the faint source is measured from the position of the bright source displaced by + 0.1° in declination. Bottom panels: normalised parallax difference vs. ρΔμ for the samples in the top panels, but now after the selection ρ > 2.2 arcsec. The three vertical lines correspond to the values of ρΔμ indicated in the top panels. The thin green curves show the 16th, 50th, and 84th percentiles of the distribution in normalised parallax difference. 
Instead, pairs need to be selected entirely based on EDR3 data. It is imperative that the parallax values themselves are not used anywhere in the selection process. Without risk of introducing selection biases, it is possible to define samples in terms of angular separation ρ, proper motion difference , parallax uncertainty and other quality indicators, and various photometric parameters. The main sample (denoted A) consists of pairs with G_{1} < 14, G_{1} < G_{2} < 19, and ρ < 10 arcsec. To correct for the contamination bias, we need in addition two comparison samples (B1 and B2), where faint components in the range G_{1} < G_{2} < 19 are selected within 10 arcsec of positions that are offset in declination by ± 0.1° from the bright components in sample A. We assume that B1 and B2 contain similar numbers and distributions of potential contaminants as A. The use of two comparison fields, rather than one, reduces the statistical noise introduced by the correction procedure, and their symmetrical placement around the bright component reduces the effect of local variations in star density, etc.
For the actual samples used, we started by extracting sources in Gaia EDR3 with fiveparameter solutions that satisfied (D.2)
which resulted in 14 561 255 sources. These are the bright components in sample A. Next, within a radius of 10 arcsec around each such source, we extracted faint components (physical and optical) with fiveparameter solutions and magnitudes G_{2} in the range [G_{1}, 19]. This selection is very incomplete for ρ ≲ 2 arcsec, mainly because such sources usually do not have a reliable ν_{eff} from the BP and RP photometry, and therefore lack a fiveparameter solution. For G_{2}− G_{1} ≳ 5 mag the selection is incomplete at much larger separations as well. For a welldefined selection we therefore in addition required ρ > 2.2 arcsec and G_{2} − G_{1} < 4 mag, which gave 3 336 571 faint components in sample A. No condition was imposed on σ_{ϖ} or RUWE for the faint components.
To obtain the faint components in B1 and B2, the same criteria were used as in A (that is G_{1} < G_{2} < 19, G_{2} − G_{1} < 4, and 2.2 < ρ < 10 arcsec), only with ρ measured from the offset positions. This gave 3 309 468 and 3 311 657 sources, respectively. Because B1 and B2 probably contain some faraway (ρ ≃ 0.1°) physical companions or members of clusters that include a brighter component in A, it can be inferred that sample A contains at least 26 000 physical pairs.
We note that these selections were made on a prerelease version of the EDR3 catalogue, which did not yet include the EDR3 photometric data but inherited the G magnitudes from DR2. Making the corresponding selections on the published EDR3 may therefore result in slightly different numbers.
The further selection of physical pairs was based on separation (ρ) and propermotion difference (Δμ). Low values of ρ and/or Δμ are clearly much more likely in physical pairs than in optical. However, orbital motion may give a significant proper motion difference in a physical pair, especially when the separation is small. On the other hand, the risk of selecting an optical pair is proportional to ρ^{2}, so that we can afford to increase the tolerance on Δμ for small separations. The product ρΔμ was found to be convenient for separating optical from physical pairs. This is illustrated in the top panels of Fig. D.1. Panel a shows Δμ versus ρ for sample A, but for illustration purposes, the diagram includes also pairs with ρ < 2.2 arcsec. Panel b is the corresponding plot for sample B1; the plot for B2 (not shown) is extremely similar. Both panels are colourcoded by the parallax difference in the pair, and it is obvious that A contains mainly physical pairs for ρΔμ ≲ 2 arcsec mas yr^{−1}, that is, below the topmost diagonal line. From a comparison of panels a and b it should be clear why a cut like ρ > 2.2 arcsec is needed: without it, the number of optical pairs in the top left corner of the diagram is grossly overestimated. Similar plots of G_{2}− G_{1} versus ρ motivate the selection G_{2} − G_{1} < 4 mag (cf. Fig. 6 of Lindegren et al. 2021).
The bottom panels of Fig. D.1 show the normalised parallax difference (ϖ_{1} − ϖ_{2} divided by the combined uncertainty)versus ρΔμ for sample A (panel c) and B1 (panel d), after applying the selection ρ > 2.2 arcsec. For ρΔμ ≲ 1 arcsec mas yr^{−1}, the normalised parallax differences roughly follow the expected normal distribution in sample A, whereas the distribution is much wider and displaced towards positive values in B1, and for ρΔμ ≫ 1 in both samples. Subsamples of A with a varying degree of contamination can be obtained by selecting on ρΔμ; specifically, we used the limits indicated by the three solid lines in Fig. D.1, corresponding to ρΔμ = 0.5, 1, and 2. For example, ρΔμ < 0.5 is the cleanest (and smallest) subsample, while ρΔμ < 2 gives a larger but more contaminated subsample. Increasing the upper limit on ρΔμ reduces the statistical uncertainty of the parallax bias, but could instead increase contamination bias. However, with a good procedure to correct for the contamination, the result should not depend systematically on the upper limit used.
The procedure to correct for contamination bias is illustrated in Fig. D.2. With the selection ρΔμ < 1 arcsec mas yr^{−1}, sample A contains 69 993 pairs and the median of ϖ_{1} − ϖ_{2} is + 4.1 ± 0.2 μas (uncertaintyestimated by bootstrapping). The distribution of the parallax differences in sample A is shown by the red histogram (h_{A}) in the top panel of the figure. With the same selection on ρΔμ, samples B1 and B2 contain 7226 and 7097 pairs, respectively. The mean of the distributions in B1 and B2, (h_{B1} + h_{B2})∕2, shown by the blue histogram, is clearly skewed towards positive values, as expected from the contaminants; its median parallax difference is + 29.5 ± 1.4 μas. On the assumption that the contaminants in sample A are similar, in number and distributions, as in B1 or B2, we obtain an estimate of the distribution for the true pairs in A by taking the difference between the observed distributions, that is, Δh ≡ h_{A} − (h_{B1} + h_{B2})∕2. This is shown by the shaded histogram in Fig. D.2, although only the part with positive count differences is visible owing to the logarithmic scale. Calculating the median of the difference histogram Δh, including negative counts, gives the corrected estimate + 3.4 ± 0.2 μas. The contamination correction is smaller than a μas in this case, but it can be much larger for other subsamples.
The bottom panel of Fig. D.2 shows the uncorrected and corrected medians as functions of the upper limit on ρΔμ. With increasing limit the sample sizes increase and the statistical uncertainties decrease (as shown by the error bars), but the increasing contamination bias is also apparent in the uncorrected medians. The corrected medians, on the other hand, remain virtually constant up to a limit of about 2 arcsec mas yr^{−1} on ρΔμ. This result is only meant to illustrate the principle of the contamination correction and ignores the complex dependences of both the parallax bias and the contamination bias on G, ν_{eff}, and β. The important point is that the correction enables us to benefit from the larger samples obtained with a less strict limit on ρΔμ, or alternatively, to use a finer division in magnitude or colour without sacrificing the accuracy of the estimated parallax bias.
Because we are only interested in the median of the distribution difference, it is not necessary to compute and subtract histograms; a simpler and more exact way is to use a weighted median. In an array of values {x_{i} } with (relative) weights {w_{i}}, the weighted median is the value , such that the sum of weights is the same on either side of . In our case x = ϖ_{1} − ϖ_{2}, and the weighted median is computed for the union of A, B1, and B2, setting w = 1 for pairs in sample A, and w = −0.5 for pairs in B1 and B2.
It is well known that the median minimises the sum of absolute deviations, that is, the L_{1} norm of the residuals. An alternative definition of the weighted median therefore is (D.3)
This expression can immediately be generalised to the multidimensional case, where an arbitrary function of the parameter vector z (such as thegeneral function described in Appendix A) is fitted to the data by weighted L_{1} norm minimisation, (D.4)
In all our analyses of the pairs (and for the results shown as black dots in the bottom panel of Fig. D.2), we used Eq. (D.4) with two additional refinements. First, the uncertainty of the resulting median or fit was estimated by bootstrap resampling of the union set. Secondly, the weights in Eq. (D.4) were adjusted to take the uncertainties of the parallax differences, , into account. This was done by setting the weights to 1∕σ_{Δϖ} in sample A and to − 0.5∕σ_{Δϖ} in B1 and B2. Weights proportional to are consistent with the L_{1} formalism in Eq. (D.4).
Fig. D.2 Illustrating the procedure for contamination bias correction. Top: distribution of parallax differences in sample A (thick red histogram) and for the mean of comparison samples B1 and B2 (thin blue histogram). The shaded grey histogram is the difference between the red and blue histograms. The dashed line is the corrected estimate of the parallax difference, equal to the median of the difference histogram. Bottom: uncorrected(open red circles) and corrected (filled black) estimates of the mean parallax difference vs. the cut in ρΔμ. For better visibility, the points have been slightly displaced sideways. 
References
 Arenou, F., Luri, X., Babusiaux, C., et al. 2018, A&A, 616, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Burnham, K., & Anderson, D. 2002, Model Selection and Multimodel Inference: A Practical InformationTheoretic Approach (New York: Springer) [Google Scholar]
 Butkevich, A. G., Klioner, S. A., Lindegren, L., Hobbs, D., & van Leeuwen, F. 2017, A&A, 603, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chan, V. C., & Bovy, J. 2020, MNRAS, 493, 4367 [CrossRef] [Google Scholar]
 Crowley, C., Kohley, R., Hambly, N. C., et al. 2016, A&A, 595, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 de Grijs, R., Wicker, J. E., & Bono, G. 2014, AJ, 147, 122 [Google Scholar]
 Evans, D. W., Riello, M., De Angeli, F., et al. 2018, A&A, 616, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fabricius, C., Bastian, U., Portell, J., et al. 2016, A&A, 595, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fabricius, C., Luri, X., Arenou, F., et al. 2021, A&A, 649, A5 (Gaia EDR3 SI) [Google Scholar]
 Gaia Collaboration (Prusti, T., et al.) 2016, A&A, 595, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaia Collaboration (Babusiaux, C., et al.) 2018a, A&A, 616, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaia Collaboration (Brown, A. G. A., et al.) 2018b, A&A, 616, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaia Collaboration (Helmi, A., et al.) 2018c, A&A, 616, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaia Collaboration (Antoja, T., et al.) 2021a, A&A, 649, A8 (Gaia EDR3 SI) [EDP Sciences] [Google Scholar]
 Gaia Collaboration (Brown, A. G. A., et al.) 2021b, A&A, 649, A1 (Gaia EDR3 SI) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaia Collaboration (Klioner, S. A., et al.) 2021c, A&A, 649, A9 (Gaia EDR3 SI) [EDP Sciences] [Google Scholar]
 Girardi, L. 2016, ARA&A, 54, 95 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hambly, N. C., Cropper, M., Boudreault, S., et al. 2018, A&A, 616, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hogg, D. W., Eilers, A.C., & Rix, H.W. 2019, AJ, 158, 147 [Google Scholar]
 Lindegren, L., & Bastian, U. 2011, EAS Pub. Ser., 45, 109 [CrossRef] [EDP Sciences] [Google Scholar]
 Lindegren, L., Lammers, U., Hobbs, D., et al. 2012, A&A, 538, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lindegren, L., Hernández, J., Bombrun, A., et al. 2018, A&A, 616, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lindegren, L., Klioner, S. A., Hernández, J., et al. 2021, A&A, 649, A2 (Gaia EDR3 SI) [EDP Sciences] [Google Scholar]
 Mason, B. D., Wycoff, G. L., Hartkopf, W. I., Douglass, G. G., & Worley, C. E. 2001, AJ, 122, 3466 [Google Scholar]
 Murphy, D. C., & May, J. 1991, A&A, 247, 202 [NASA ADS] [Google Scholar]
 Pietrzyński, G., Graczyk, D., Gallenne, A., et al. 2019, Nature, 567, 200 [Google Scholar]
 Riello, M., De Angeli, F., Evans, D. W., et al. 2021, A&A, 649, A3 (Gaia EDR3 SI) [Google Scholar]
 Ripepi, V., Molinaro, R., Musella, I., et al. 2019, A&A, 625, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rowell, N., Davidson, M., Lindegren, L., et al. 2021, A&A, 649, A11 (Gaia EDR3 SI) [EDP Sciences] [Google Scholar]
 Schönrich, R., McMillan, P., & Eyer, L. 2019, MNRAS, 487, 3568 [Google Scholar]
 Taylor, M. B. 2005, ASP Conf. Ser., 347, 29 [Google Scholar]
 Zinn, J. C., Pinsonneault, M. H., Huber, D., & Stello, D. 2019, ApJ, 878, 136 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Coefficients for the function ΔZ(G, ν_{eff}, β) fitted to the difference in parallax between EDR3 and DR2.
Coefficients for the function Z_{5}(G, ν_{eff}, β) fitted to the quasar parallaxes.
Coefficients for the extended function Z_{5}(G, ν_{eff}, β, X) fitted to the quasar parallaxes.
Comparison of the parallax gradient in colour, p_{1}, at the position of the LMC, as estimated from LMC data, quasars (QSO), and physical pairs (PP).
Coefficients for the extended function Z_{5}(G, ν_{eff}, β, X) as obtained in a combined fit to the quasar and LMC samples.
Coefficients of Z_{5}(G, ν_{eff}, β) as estimated from physical pairs for G > 10.
All Figures
Fig. 1 Fraction of observations in different modes. The shaded area represents gated observations; the curves show the fraction of AF observations made in window classes WC0a, WC0b, WC1, and WC2. 

In the text 
Fig. 2 Celestial maps in ecliptic coordinates of some selected statistics in Gaia EDR3. Panel a: mean parallax uncertainty for sources with fiveparameter solutions and G < 14 mag. Panel b: mean correlation coefficient between parallax and pseudocolour for sources with sixparameter solutions. Panel c: smoothed median parallax of quasars, corrected for the global median offset of − 17 μas. Owing to the small number of quasars identified at small Galactic latitudes, no data are shown for  sin b  < 0.1. The maps use a Hammer–Aitoff projection in ecliptic coordinates with λ = β = 0 at the centre, ecliptic north up, and ecliptic longitude λ increasing from right to left. All three statistics exhibit some degree of systematic dependence on ecliptic latitude, which may be symmetric (as in a) or antisymmetric (as in b and c), but also even larger variations as functions of both coordinates. 

In the text 
Fig. 3 Mean difference in parallax between EDR3 and DR2. Top: parallax difference vs. magnitude. Middle and bottom: parallax difference vs. ecliptic latitude in two magnitude intervals. Circles connected by solid lines are weighted mean values computed in bins of variable size, with at least 1000 sources per bin; dashed lines are mean values of the fitted parametrised function ΔZ (Eq. (5)), binned as for the sources. Filled red circles are for ν_{eff} < 1.48 μm^{−1}, and open blue circles for ν_{eff} > 1.48 μm^{−1}. 

In the text 
Fig. 4 Mean difference in parallax between EDR3 and DR2 as a function of effective wavenumber for several ranges of the G magnitude.Circles connected by solid lines are weighted mean values computed in bins of variable size, with at least 1000 sources per bin, and dashed lines are mean values of the fitted parametrised function ΔZ (Eq. (5)), binned as for the sources. Filled red circles are for sources with β > 0, and open blue circles for β < 0. The vertical dashed lines mark the breakpoints for the basis functions c_{j}(ν_{eff}) in Eq. (A.3), i.e. the clamping limits at 1.24 and 1.72 μm^{−1} and the midpoint at 1.48 μm^{−1}. 

In the text 
Fig. 5 Mean parallax of quasars binned by magnitude (panel a), effective wavenumber (panel b), and sine of ecliptic latitude (panel c). In each bin the dot is the mean parallax in EDR3 weighted by , with error bars indicating the estimated standard deviation of the weighted mean. The two red lines indicate the weighted mean value (− 21 μas) and median (− 17 μas) of the full sample. 

In the text 
Fig. 6 Example of interactions among the dependences of quasar parallaxes on colour and ecliptic latitude. In a quadratic regression of parallax vs. sinβ, the linear coefficient (filled blue squares) exhibits a strong variation with effective wavenumber, while no such trend is shown by the quadratic coefficient (open red circles). The points have been slightly displaced sideways to avoid overlapping error bars. 

In the text 
Fig. 7 Mean residual in quasar parallax after a regression on G, ν_{eff}, and β (see text), plotted against the confusion factor from Eq. (6). In each bin the dot is the mean residual weighted by , with error bars indicating the estimated standard deviation of the weighted mean. The broken dashed line is the dependence modelled by Eq. (7). 

In the text 
Fig. 8 Celestial map of the mean contamination bias of the quasars, as given by Eq. (8) and Table 3. The map uses the same projection in ecliptic coordinates as the maps in Fig. 2. 

In the text 
Fig. 9 Colourmagnitude diagrams of the LMC sample. Left: diagram colourcoded by the number of sources at a given position in the diagram. Right: diagram colourcoded by the median parallax in EDR3. 

In the text 
Fig. 10 Median parallaxes for the LMC sample in Fig. 9 plotted against ν_{eff} in six magnitude intervals, as indicated in the diagrams. The orange points show the parallaxes of individual sources and give an impression of the coverage in ν_{eff} and scatter in parallax. The black dots, with error bars, show the median parallax and its uncertainty in bins of ν_{eff} with at least 20 sources per bin. The solid blue curves show the fitted combination of basis functions in Eq. (9) evaluated for a representative magnitude (G = 10.0, 11.5, 12.5, 15.0, 17.0, and 18.0) in each interval. 

In the text 
Fig. 11 Coefficients q_{jk} estimated from physical pairs as functions of G. Results for q_{12} are considered insignificant and set to zero when fitting the other five coefficients. The different symbols represent different selections onρΔμ: 0–0.5 arcsec mas yr^{−1} (red circles), 0.5–1 arcsec mas yr^{−1} (green triangles), 1–2 arcsec mas yr^{−1} (blue squares), and 0–2 arcsec mas yr^{−1} (filled black circles with lines and error bars). The dashed blue line is the global fit in Table 7. The vertical grey lines show the breakpoints for the basis functions in G defined by Eq. (A.2). 

In the text 
Fig. 12 Fraction of sources in Gaia EDR3 with different kinds of solutions: fiveparameter solutions (P5 = solid red curve), sixparameter solutions (P6 = dashed blue curve), and twoparameter solutions (P2 = thin grey curve). The twoparameter solutions are ignored in this paper as they do not have parallaxes. 

In the text 
Fig. 13 Median parallax difference between six and fiveparameter solutions as a function of ν_{eff} and G for the 8 million sources with both kinds of solutions. The panels show selections depending on ecliptic latitude β: southern hemisphere (panel a), all latitudes (panel b), and northern hemisphere (panel c). In this sample no sources lie outside of the interval 1.24 < ν_{eff} < 1.72 μm^{−1}. 

In the text 
Fig. 14 Same data as in Fig. 13, but plotted against the astrometrically determined pseudocolour (). Uncertainties in the pseudocolour scatter some points outside of the interval 1.24–1.72 μm^{−1}. 

In the text 
Fig. 15 Median ẑ_{6} as a functionof ν_{eff} and G for the sample in Figs. 13 and 14. The panels show selections depending on ecliptic latitude: southern hemisphere (panel a), all latitudes (panel b), and northern hemisphere (panel c). 

In the text 
Fig. 16 Parallax bias Z_{6} according toTable 10. The panels show mean values for the southern hemisphere (panel a), all latitudes (panel b), and the northern hemisphere (panel c). 

In the text 
Fig. 17 Parallaxes for 1.2 million quasars with fiveparameter solutions in Gaia EDR3. Yellow dots show theindividual values plotted vs. magnitude (panel a), effective wavenumber (panel b), and sine of ecliptic latitude (panel c). Open black circles show mean values of the uncorrected parallaxes (ϖ) in bins of magnitude etc., and filled blue circles show mean values of the corrected parallaxes (ϖ − Z_{5}) in the same bins. Mean values are calculated using weights . Error bars indicate the estimated standard deviation of the weighted mean in each bin. 

In the text 
Fig. 18 Parallaxes for 0.4 million quasars with sixparameter solutions in Gaia EDR3. Yellow dots show the individual values plotted vs. magnitude (panel a), pseudocolour (panel b), and sine of ecliptic latitude (panel c). Open black circles show mean values of the uncorrected parallaxes (ϖ) in bins of magnitude etc., and filled blue circles show mean values of the corrected parallaxes (ϖ − Z_{6}) in the same bins. Red open squares show mean values of the corrected parallaxes for the 78% of the quasars that have insignificant excess source noise, using a slightly different binning. Mean values are calculated using weights . Error bars indicate the estimated standard deviation of the weighted mean in each bin. 

In the text 
Fig. 19 Colourmagnitude diagram of the LMC sample colourcoded by the median of the EDR3 parallax after subtracting Z_{5} as given by the cofficients in Table 9. 

In the text 
Fig. 20 Parallax bias Z_{5} (top) and Z_{6} (bottom) computed according to Tables 9 and 10 for a representative sample of sources with five and sixparameter solutions in EDR3. The dots show values for the complete sample of sources (for G < 11.5) or for a random selection (G > 11.5). The colour scale indicates the mean ν_{eff} or at a given point, thus giving an impression of the mean colour dependence of the bias. The black curves show the 1st, 10th, 50th, 90th, and 99th percentiles. The thick curve is the 50th percentile or median. 

In the text 
Fig. 21 Parallax bias Z_{5}(G, ν_{eff}, β) according to Table 9. The panels show cuts at (a) ecliptic latitude β = − 90°, (b) β = 0°, and (c) β = +90°. 

In the text 
Fig. 22 Parallax bias according to Table 10. The panels show cuts at (a) ecliptic latitude β = − 90°, (b) β = 0°, and (c) β = +90°. 

In the text 
Fig. A.1 Basis functions g_{i}(G), i = 0… 12 according toEq. (A.2). 

In the text 
Fig. A.2 Basis functions c_{j}(ν_{eff}), j = 1… 4 according to Eq. (A.3). The (constant) function c_{0} = 1 is not shown. For better visibility the functions are vertically displaced by j∕2, and the amplitude of c_{2}(ν_{eff}) is increased by a factor 10. 

In the text 
Fig. A.3 Basis functions b_{k}(β), k = 0… 2 according toEq. (A.4). 

In the text 
Fig. B.1 Illustrating the selection and validation of the LMC sample. (a) Colourmagnitude diagram for the original sample centred on (l_{C}, b_{C}). (b) Joint distribution for the original sample of parallax and Δμ_{N}, the normalised propermotion difference to the fitted model. The vertical line marks the cutoff used for the selection based on proper motions. (c) Colourmagnitude diagram for the subsample with Δμ_{N} < 10. (d–f) Same as (a–c), but for a sample centred on (l_{C} − 10°, b_{C}), containing far fewer LMC stars, but roughly the same number of Galactic foreground stars as in (a–c). 

In the text 
Fig. C.1 Hertzsprung–Russell diagrams for the union of samples N and S. Top: absolute magnitudes in the G band. Bottom: absolute magnitudes from W in Eq. (C.1). 

In the text 
Fig. C.2 Differential parallax bias estimated by means of RC stars. The plots show differences between the measured parallaxes ϖ_{EDR3} and the photometric parallaxes ϖ_{RC} from Eq. (C.2) in area N (a) and S (b). Contours of constant density are shown as thin grey lines. The thick black curve traces the ridge of the feature created by the red clump stars. 

In the text 
Fig. C.3 Differential parallax bias estimated by means of RC stars. The points show the locations of the ridges in Fig. C.2 for area N (filled red circles) and S (open blue squares). The curves show the expected locations at G = 17 according toTable 6 for area N (solid red) and S (longdashed blue), assuming that the curvature with colour (j = 2) is the same in both areas. If the curvature had flipped sign with β, the ridge location for area N should instead follow the upbending, shortdashed black curve. 

In the text 
Fig. D.1 Illustration of the selection of physical pairs and the contamination by optical pairs. Top panels: proper motion difference (Δμ) vs. separation (ρ) for source pairs in sample A (a) and B1 (b) before making the cut ρ > 2.2 arcsec indicated by the dashed vertical line. The three solid lines correspond to ρΔμ = 0.5, 1, and 2 arcsec mas yr^{−1}. The points are colourcoded by the median parallax difference between the components. In panel b the separation of the faint source is measured from the position of the bright source displaced by + 0.1° in declination. Bottom panels: normalised parallax difference vs. ρΔμ for the samples in the top panels, but now after the selection ρ > 2.2 arcsec. The three vertical lines correspond to the values of ρΔμ indicated in the top panels. The thin green curves show the 16th, 50th, and 84th percentiles of the distribution in normalised parallax difference. 

In the text 
Fig. D.2 Illustrating the procedure for contamination bias correction. Top: distribution of parallax differences in sample A (thick red histogram) and for the mean of comparison samples B1 and B2 (thin blue histogram). The shaded grey histogram is the difference between the red and blue histograms. The dashed line is the corrected estimate of the parallax difference, equal to the median of the difference histogram. Bottom: uncorrected(open red circles) and corrected (filled black) estimates of the mean parallax difference vs. the cut in ρΔμ. For better visibility, the points have been slightly displaced sideways. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.