

\section{Introduction} 
\label{sec:intro}

\textit{Gaia} Early Data Release 3 (EDR3; \citeads{EDR3-DPACP-130}) contains provisional 
astrometric and photometric data for more than 1.8~billion ($1.8\times 10^9$) sources
based on the first 34~months of observations made by the European Space Agency's 
\textit{Gaia} mission (\citeads{2016A&A...595A...1G}) since the start of the nominal 
operations in July 2014. The astrometric data in EDR3 include the five astrometric 
parameters (position, parallax, and proper motion) for 1.468~billion sources, and the 
approximate positions at epoch J2016.0 for an additional 344~million mostly faint 
sources. All sources have magnitudes in \textit{Gaia}'s unfiltered photometric 
passband $G$, and 1.544~billion have two-colour photometry in the passbands 
$G_\text{BP}$ and $G_\text{RP}$ defined by the blue and red photometers 
(BP and RP; \citealt{EDR3-DPACP-117}). The magnitudes of the well-observed sources range 
from $G=6$ to 21. All data are publicly available in the online \textit{Gaia} Archive at 
\href{https://archives.esac.esa.int/gaia}{\tt https://archives.esac.esa.int/gaia}.

The EDR3 is a subset of the full \textit{Gaia} Data Release 3 (DR3), planned for the 
first half of 2022. The full release will provide a much wider set of data, including
detailed spectrophotometric and variability information, additional astrometric 
data on non-single and extended objects, and the radial velocities, object classification, 
and astrophysical parameters for many sources. However, the basic astrometric 
information on the \textit{Gaia} DR3 sources, obtained by treating all of them as 
single stars, has already been provided in EDR3 and will not change for DR3.
  
This paper gives an overview of the processing leading up to the EDR3 astrometry,
as well as of the main characteristics of the astrometric results. Further details are provided 
in the online documentation of the \textit{Gaia} Archive and in specialised papers. In 
particular, the celestial reference frame of \textit{Gaia} (E)DR3 is described in 
\citet{EDR3-DPACP-133}, the parallax bias (zero point) is discussed in
\citet{EDR3-DPACP-132}, and the overall properties of the release are reviewed in
\citet{EDR3-DPACP-126}. A general description of the \textit{Gaia} mission can be found
in \citetads{2016A&A...595A...1G}.

The core astrometric solution for \textit{Gaia}, known as AGIS (astrometric global 
iterative solution), was comprehensively described in the pre-launch paper by 
\citetads{2012A&A...538A..78L}. This remains a useful general reference for AGIS 
in spite of the many modifications and improvements introduced since 2012. 
We also refer frequently to \citetads{2018A&A...616A...2L},
which describes the astrometric solution for \textit{Gaia} DR2.


\begin{figure*}
\center
  \includegraphics[width=16cm]{fig01.pdf} 
    \caption{Main steps of the EDR3 astrometry processing and their place in the
    cyclic processing scheme of DPAC. \textit{Gaia} EDR3 (and DR3) are generated in
    the third processing cycle (cycle~3). The stretches of observational data processed 
    in the different cycles are indicated by thick horizontal lines. The boxes connected by arrows 
    show the sequence of processing steps and their interdependencies, but they are not placed
    chronologically on the timeline. Only steps directly relevant for the astrometry are shown,
    leaving out most of the complexities of the full DPAC processing. No details are given for DR1.
    The first month of 
    the nominal mission, with observations made in the ecliptic pole scanning law (EPSL) mode, 
    was not used for the astrometry in DR2 and EDR3, but may be incorporated in later releases. 
    The Whitehead eclipse avoidance manoeuvre (WEAM) on 16 July 2019 marks the beginning
    of the extended mission. In the first year of the extended mission (data segments DS6 
    and DS7), scanning was made in the reversed precession mode (Sect.~\ref{sec:acRate}). 
    The processes SDM, CALIPD, AGIS, and PhotPipe are explained in Sect.~\ref{sec:tasks}.
    }
    \label{fig:timeline}
\end{figure*}




\section{Overview of the astrometric processing}
\label{sec:overview}

\subsection{Main processing tasks}
\label{sec:tasks}

In the cyclic processing scheme adopted by the \textit{Gaia} Data Processing and Analysis 
Consortium (DPAC; \citeads{2016A&A...595A...1G}), EDR3 and DR3 are products of the 
third processing cycle, using observations in the first four data segments called DS0--DS3
in Fig.~\ref{fig:timeline}. The data segments are just a convenient, but essentially
arbitrary division of the raw data by acquisition time. As suggested in the figure, the cycles 
treat successively larger chunks of the raw data by including additional segments, but in 
every cycle the old segments are always reprocessed together with the new ones. This 
iterative reprocessing of earlier data segments is necessary in order to achieve the 
uniformly best treatment of all the data, and a consistent assignment of source 
identifiers to the on-board detections.

In Fig.~\ref{fig:timeline}, the boxes labelled PhotPipe represent the complex photometric 
processing described elsewhere \citep{EDR3-DPACP-117,DPACP-118,DPACP-119,DPACP-120}.
This is not part of the astrometric processing as such but is included in the diagram because 
the photometric information, and in particular the colour information encoded in the 
effective wavenumbers ($\nu_\text{eff}$; Sect.~\ref{sec:nuEff}) calculated in PhotPipe, are 
needed for calibrating the colour-dependent line-spread and point-spread functions (LSF 
and PSF) of the astrometric instrument. Because PhotPipe runs essentially in parallel with 
the astrometric solution (AGIS), this implies that the astrometric processing in cycle~$N$ 
must use photometric information from cycle~$N{-}1$.

The boxes in Fig.~\ref{fig:timeline} labelled SDM, CALIPD, and AGIS represent the three 
main stages in the processing of the raw CCD (charge-coupled device) data that are of 
immediate relevance for the astrometry:

In the first stage, the SDM (source, detection-classifier, and cross-match) aims to identify all on-board
detections belonging to the same source and assign a unique source identifier (\gacs{source\_id}) 
to each such cluster of detections \citep{EDR3-DPACP-124}. An important part of the process is the 
identification of spurious detections, created for example by the diffraction spikes of bright stars 
\citepads{2016A&A...595A...3F}. 
Because the updated source list and table of links to the (genuine) detections created by the 
SDM is used by all subsequent processes, this is one of the first tasks to be executed in a cycle.
The source list from the previous cycle is a starting point for the task, but the new data and 
improved reconstruction of the satellite attitude (a key element in translating observed 
transit times into positions) unavoidably require some of the old sources to be split or merged,
in addition to creating entirely new ones. For example, EDR3 contains many pairs of sources 
(most of which are genuine binaries) that are separated by less than 0.4~arcsec, where DR2 had only one. 
Such cases could lead to the assignment of new source identifiers for both components.
The auxiliary table \gacs{dr2\_neighbourhood} helps to trace the evolution of source identifiers.

The second stage, CALIPD, consists of two parts, calibration (CAL) and image parameter 
determination (IPD). In CAL, the LSF (for one-dimensional observations) and PSF 
(for two-dimensional observations) are calibrated as functions of time, colour, 
and several other variables in order to take into account the optical imperfections of the 
instrument and their temporal evolution \citep{EDR3-DPACP-73}. The LSF and PSF 
describe the shape of the image profile for a point source as well as the small 
displacement caused by chromatic effects (Sect.~\ref{sec:nuEff}).
In IPD, the LSF or PSF relevant for a particular observation is fitted to the sampled CCD image, 
yielding precise estimates of its one- or two-dimensional location in the pixel stream 
and of the total flux of the image in the $G$ band \citepads{2016A&A...595A...3F}. 
The resulting image locations constitute the main input data for the astrometric solution,
while the flux estimates are used for the determination of $G$ magnitudes in PhotPipe.
Whereas CAL only uses a small fraction of the available observations for the LSF and PSF
calibrations, IPD is applied to all observations in the skymapper (SM) and
astrometric field (AF).
 
In the third stage, AGIS performs a simultaneous least-squares
estimation of the attitude, instrument calibration, and the five astrometric parameters for
a subset of well-behaved primary sources (about 14.3~million in cycle~3). The calibration
includes corrections for effects that are not accounted for in the CALIPD, or only partially
corrected at that stage. The comprehensive pre-launch description of AGIS in 
\citetads{2012A&A...538A..78L} is complemented by specifics of the
current models in Sect.~\ref{sec:models}.

SDM and CALIPD belong to the intermediate data update (IDU) system, which includes 
several additional tasks such as astrophysical background estimation 
\citepads{2016A&A...595A...3F} and electronic calibrations \citepads{2018A&A...616A..15H}.
Compared with DR2, several major improvement of the IDU have been introduced with
cycle~3. In SDM the treatment of high-proper motion stars, very bright stars, variable 
sources, and close pairs has been much improved \citep{EDR3-DPACP-124}. In CALIPD the 
image profiles (LSF and PSF) are no longer assumed to be independent of time and
colour, as was the case in cycle~2, and a much more realistic two-dimensional model (PSF) is 
used \citep{EDR3-DPACP-73}. Moreover, as described in Sect.~\ref{sec:nuEff} and
Fig.~\ref{fig:timeline}, the CALIPD and AGIS tasks are for the first time iterated in 
order that CALIPD may benefit from the improved astrometry, attitude, and instrument 
calibration obtained by including the new data segment (DS3) in AGIS.



\subsection{Observations used}
\label{sec:data}

\textit{Gaia} EDR3 is based on data collected from the start of the nominal observations 
on 25 July 2014 (10:30 UTC) until 28 May 2017 (08:45 UTC), or 1038~days
(data segments DS0--DS3 in Fig.~\ref{fig:timeline}).
Similarly to the astrometric solution for DR2 \citepads{2018A&A...616A...2L}, this
solution did not use the observations in the first month of the operational phase, when
the special ecliptic pole scanning law (EPSL) was employed. The data for the astrometry 
therefore start on 22 August 2014 (21:00 UTC) and cover 1009~days or 2.76~yr, with some 
interruptions mentioned below. 

The time coverage for this solution is therefore about one year longer than the astrometric 
solution for \textit{Gaia} DR2, which covered 640~days or 1.75~year. The expected 
improvement from the added data and longer time baseline scales as $T^{-1/2}$ for the
parallaxes and positions at the mean epoch of observation, and as $T^{-3/2}$ for the
proper motions; thus uncertainties should be smaller by a factor 0.80 for the parallaxes
and positions, and by a factor 0.51 for the proper motions. As shown in Sect.~\ref{sec:form}, the
median ratios of the formal uncertainties are slightly better than this thanks to additional
improvements in the instrument and attitude modelling. The reference epoch J2016.0 used 
for the astrometry in \textit{Gaia} EDR3 (Sect.~\ref{sec:source}) is close to the mid-point of 
the observations.

The on-board mission timeline (OBMT) is conveniently used to label on-board events; it is
expressed as the number of nominal revolutions of exactly 21\,600~s (6~h) on-board time from 
an arbitrary origin.%
\footnote{The rubidium atomic clock on board of \textit{Gaia} does not count SI seconds because 
it is a free-running oscillator with some (very small) time-dependent frequency error. This is 
calibrated in a special part of the data processing (see Sect.~\ref{sec:aux}), but ignored when giving 
intervals in OBMT.}
The approximate relation between OBMT (in revolutions) and barycentric
coordinate time (TCB, in Julian years) at \textit{Gaia} is
\begin{equation}\label{eq:obmt}
\text{TCB} \simeq \text{J}2015.0 + (\text{OBMT} - 1717.6256~\text{rev})/(1461~\text{rev~yr}^{-1}) \, ,
\end{equation} 
or as a Julian Date,
\begin{equation}\label{eq:obmtJD}
\text{JD}2457023.75 + (\text{OBMT} - 1717.6256~\text{rev})/(4~\text{rev~d}^{-1}) \, .
\end{equation} 
The nominal observations start at OBMT 1078.38~rev. The astrometric solution used
data in the interval OBMT 1192.13--5230.09~rev (J2014.64032--J2017.40415), with 
major gaps as listed in Table~\ref{tab:gaps}.

\begin{table}
\caption{Major gaps and events affecting the astrometric solution.
\label{tab:gaps}}
\small
\begin{tabular}{ccrl}
\hline\hline
\noalign{\smallskip}
$t_\text{beg}$ & $t_\text{end}$ & length & description\\
\noalign{\smallskip}
\hline
\noalign{\smallskip}
1220.400 & 1225.200 & 4.800 & VPU reset \\
1316.490 & 1389.113 & 72.623 & decontamination {\#}4 \\
1443.800 & 1444.200 & 0.400 & refocus (following FoV) \\
1653.800 & 1660.000 & 6.200 & PAA anomaly \\
1820.900 & 1830.000 & 9.100 & PAA anomaly \\
2094.000 & 2099.000 & 5.000 & PDHU anomaly \\
2179.125 & 2191.000 & 11.875 & VPU software update \\
2192.252 & 2195.218 & 2.967 & observation gap\\
2238.000 & 2242.000 & 4.000 & unknown \\
2322.300 &	2401.559 & 79.259 & decontamination {\#}5 \\
2405.967 & 2408.643 & 2.676 & observation gap\\
2408.935 & 2409.968 & 1.033 & observation gap\\
2574.640 &	2575.400 & 0.760 & refocus (preceding FoV) \\
2954.200 & 2958.000 & 3.800 & moon eclipse \\
3045.133 & 3049.000 & 3.867 & PDHU anomaly\\
3603.250 & 3605.227 & 1.976 & observation gap\\
3646.800 & 3650.000 & 3.200 & moon eclipse \\
3663.700 & 3667.000 & 3.300 & moon eclipse \\
4074.210 & 4076.063 & 1.853 & observation gap\\
4112.463 & 4180.000 & 67.537 & decontamination {\#}6 \\
4271.753 & 4275.200 & 3.445 & PAA anomaly \\
4477.441 & 4481.000 & 3.550 & PAA anomaly \\
4512.502 & 4515.000 & 2.498 & PAA anomaly \\
4545.144 & 4548.000 & 2.856 & PAA anomaly \\
5078.547 & 5080.600 & 2.053 & STR anomaly \\
\noalign{\smallskip}
\hline
\end{tabular}
\tablefoot{The table lists gaps longer than 1~revolution (0.25~day) and some other events of relevance
for the calibration model. $t_\text{beg}$ and $t_\text{end}$ are the start and end times of the
gap in OBMT revolutions (see Eq.~\ref{eq:obmt}). The third column is the length of the gap in 
revolutions. Abbreviations: 
VPU = video processing unit,
FoV = field of view,
PAA = phased array antenna,
PDHU = payload data handling unit,
STR = star tracker.
}
\end{table}

\subsection{Use of colour information in CALIPD and AGIS}
\label{sec:nuEff}

In the focal plane of an all-reflecting telescope, free of wavefront aberrations, the 
point-spread function (PSF) is completely symmetric. Although the width 
of the PSF increases with wavelength, because of diffraction, its position does 
not change and is consequently
independent of the spectral composition of the light (achromatic). This is no 
longer true for a real instrument like \textit{Gaia}. Inevitable coma-like wavefront 
errors produce asymmetric PSFs, in which both the shape and location 
depend on the spectrum. Subtle wavelength-dependent effects can also be 
introduced by the CCD detector itself. We use `chromaticity' as a generic term 
for these several effects, but especially for the variation of the PSF location 
with colour. Chromaticity creates colour-dependent biases in the astrometric
results, unless it is properly calibrated and corrected for in the processing.

Chromaticity should ideally be completely eliminated already in CALIPD, so that
the astrometric solution (AGIS) would not need to care about the sources having 
different colours. This requires (i) that in CAL both the shape and location of the 
LSF or PSF are accurately calibrated as functions of the spectral energy 
distribution (multiplied by the wavelength passband); and (ii) that in IPD the location 
and flux of the image are estimated using the correct profile, depending on 
the actual spectrum of the source in each observation. The astrometric parameters 
determined in the subsequent AGIS solution will then be free from chromatic 
biases. There is a certain circularity here: To achieve (i), CAL must be able to identify 
the point in the image profile that corresponds to the achromatic centre of the source, 
and this can only be done by means of the (achromatic) astrometric parameters 
determined by AGIS. This strong interdependency between CALIPD and AGIS is
dealt with by executing the two tasks alternately, which motivates the sequence 
CALIPD~3.1, AGIS~3.1, CALIPD~3.2, AGIS~3.2 in Fig.~\ref{fig:timeline}. The 
CALIPD/AGIS sequence should ideally be iterated until convergence, but in cycle~3
only two iterations (3.1 and 3.2) were made. This appears to be sufficient in practice,
because AGIS is able to eliminate most of the chromatic effects left uncorrected in IPD
via the colour-dependent terms in the AGIS calibration model (Sect.~\ref{sec:cal}).

In cycle~3 two simplifying assumptions are made, both of which may be relaxed at
some future time. The first is that the spectral information needed for the chromaticity
correction is fully encoded
in the effective wavenumber, defined as $\nu_\text{eff}=\langle\lambda^{-1}\rangle$. 
Here $\lambda$ is the wavelength, and angular brackets denote a mean value 
weighted by the detected photon flux per unit wavelength interval. 
This quantity was chosen, in preference to (say) the effective wavelength or 
colour index, based on pre-launch studies using the properties of the \textit{Gaia} 
instrument as expected at the time. According to these studies, the effective wavenumber 
provides a good one-dimensional parametrisation of chromaticity for ordinary 
stellar spectra, but may not be enough to describe shifts at the few $\mu$as 
level in atypical cases such as quasar spectra. Thus, more complex dependencies on 
the source spectrum may have to be considered in the future, but for the time being 
we use $\nu_\text{eff}$ as defined.


The second assumption is that the spectrum (or effective wavenumber) is the 
same in all observations of a given source. Although this is a sufficiently good 
approximation for most sources, it may prevent us from reaching the full potential 
of \textit{Gaia} for some variable objects. The remedy is simple in principle, namely
to use the actual colour of the source at each observation, but this may require 
an additional iteration over PhotPipe and the variability analysis 
\citepads{2018A&A...618A..30H}.

Even with the simplifications mentioned above, the conditions for eliminating 
chromaticity in CALIPD are not fully met in cycle~3. The main obstacle is that many 
sources do not have reliable colour information that can be used to select the 
appropriate image profiles for the IPD. The effective wavenumbers used in 
CALIPD~3.1 and 3.2 were calculated in PhotPipe~2 directly from the sampled and 
calibrated mean BP and RP spectra, and are given in EDR3 as 
\gacs{nu\_eff\_used\_in\_astrometry}. The analysis of the BP and RP spectra 
is very challenging in crowded areas and at the faintest magnitudes, owing to 
the blending of overlapping spectra and the difficulty to estimate the background 
accurately \citep{EDR3-DPACP-118}. A strict filtering on the quality of 
$\nu_\text{eff}$ was adopted in order to avoid that biases in the photometric
colour might propagate into the astrometry. Of particular concern was the BP+RP 
flux excess issue \citepads{2018A&A...616A...4E}, which in DR2 tended to make faint 
sources in crowded areas too blue. As a result of the adopted filtering, about
two thirds of the sources in EDR3 do not have a valid \gacs{nu\_eff\_used\_in\_astrometry}.
The situation is more favourable for brighter sources, where, for example, only 
12\% of the sources with $G<18$~mag lack a valid $\nu_\text{eff}$.

For the many sources without a valid $\nu_\text{eff}$, special procedures were 
used both in IPD and AGIS. In the IPD, image parameters 
were estimated by fitting the calibrated LSF or PSF for the default wavenumber
$\nu_\text{eff}^\text{\,def}=1.43~\mu$m$^{-1}$. This value was chosen to be
close to the mean $\nu_\text{eff}$ of faint sources, for which the default value 
is mostly used; thus, the averaged error introduced by the procedure is minimised. 
In AGIS, a six-parameter solution was computed for these sources, where the sixth
unknown, after the standard five astrometric parameters, is the pseudocolour.
This quantity, denoted $\hat{\nu}_\text{eff}$, is an astrometric estimate of the 
effective wavenumber $\nu_\text{eff}$. In order to estimate the pseudocolour, 
it is assumed that the chromatic shift of the image location caused by using the 
wrong colour (that is, the default colour) in IPD is linearly proportional to 
${\nu}_\text{eff}-\nu_\text{eff}^\text{\,def}$. The constant of proportionality 
is a property of the instrument that can be determined in a special AGIS calibration 
solution using sources for which $\nu_\text{eff}$ is known (step~\ref{step:Cprime}
in Sect.~\ref{sec:steps}). 

The sources with six-parameter solutions in EDR3 are identified by the flag  
$\gacs{astrometric\_params\_solved}=95$. The estimated $\hat{\nu}_\text{eff}$ 
(expressed in $\mu$m$^{-1}$) is given as \gacs{pseudocolour}; like the other 
astrometric parameters it comes with a formal uncertainty (\gacs{pseudocolour\_error}) 
and correlation coefficients (\gacs{ra\_pseudocolour\_corr}, etc.). The full $6\times 6$ 
covariance matrix can thus be reconstructed, which makes it possible to compute 
improved estimates of the astrometric parameters if a better estimate of the colour  
than the pseudocolour is available (see Appendix~\ref{sec:col6p}). We note that 
\gacs{nu\_eff\_used\_in\_astrometry} is not given for the sources with six-parameter 
solutions.

Conversely, sources with a standard five-parameter solution 
($\gacs{astrometric\_params\_solved}=31$) have the field
\gacs{nu\_eff\_used\_in\_astrometry} set, but no \gacs{pseudocolour}. 
Neither colour field is set for sources that have only a position in EDR3
($\gacs{astrometric\_params\_solved}=3$). 

\begin{figure}
\centering
  \includegraphics[width=0.9\hsize]{fig02.pdf} 
    \caption{Relation between the colour index and effective wavenumber for a
    random sample of 1.5~million sources in EDR3 brighter than $G=18$. The
    dashed curve is the approximate mean relation in Eq.~(\ref{eq:nuEff1}) or (\ref{eq:nuEff2}).}
    \label{fig:nuEffVsBpRp}
\end{figure}

The relation between the colour index  $G_\text{BP}-G_\text{RP}$ (\gacs{bp\_rp})
and effective wavenumber $\nu_\text{eff}$ (\gacs{nu\_eff\_used\_in\_astrometry}) in EDR3
is illustrated in Fig.~\ref{fig:nuEffVsBpRp}. As shown by the diagram, there is no unique 
one-to-one relation between the two colour parameters. One reason is that $\nu_\text{eff}$ 
for cycle~3 was computed in the previous cycle (by PhotPipe~2), and is therefore not 
completely consistent with other photometric data in EDR3, including the colour indices. 
But the main reason for the scatter is the very different methods of computation 
($\nu_\text{eff}$ as a weighted sum over the sampled BP and RP spectra, 
$G_\text{BP}-G_\text{RP}$ from the integrated BP and RP fluxes), which give slightly 
different results depending on the detailed spectra. When an approximate relation
is needed, the following analytical formulae may be useful:
\begin{gather}\label{eq:nuEff1}
\nu_\text{eff} \simeq 1.76 - \frac{1.61}{\pi}\,\text{atan}\,
\Bigl(0.531(G_\text{BP}-G_\text{RP})\Bigr)\quad \mu\text{m}^{-1}\, , \\
G_\text{BP}-G_\text{RP} \simeq \frac{1}{0.531}\,\text{tan}\,
\Biggl(\frac{\pi}{1.61}\,(1.76-\nu_\text{eff})\Biggr)\quad \text{mag} \, . \label{eq:nuEff2}
\end{gather}
For $-0.5\le G_\text{BP}-G_\text{RP}\le 7$ they represent the mean relation for stellar
objects to within $\pm 0.007~\mu$m$^{-1}$ in the effective wavenumber. 
The atan/tan functions conveniently describe the non-linear relation to a useful approximation,
and has the additional advantage that $\nu_\text{eff}$ is restricted to the physically plausible
interval $[0.955,~2.565]~\mu$m$^{-1}$ for arbitrarily large (positive or negative) colour 
indices. 


\subsection{Auxiliary data}
\label{sec:aux}

The processing of \textit{Gaia} data aims at producing the most accurate astrometric catalogue 
consistent with the observations, using a minimum of external auxiliary data. Some external data 
are nevertheless needed, for example to align the catalogue with the celestial reference system 
and correct for stellar aberration. The main auxiliary data used in the processing are described below. 

\paragraph{Reference frame} The orientation of the axes of the International Celestial Reference System 
(ICRS) is conventionally defined by means of the accurate positions for extragalactic radio sources observed
by very long baseline interferometry. As of 1~January 2019, the defining list is the third realisation of 
the International Celestial Reference Frame (ICRF3; \citeads{2020A&A...ICRF3}) containing 4588
radio sources. The orientation of \textit{Gaia}-CRF3, the celestial reference frame of \textit{Gaia} EDR3
\citep{EDR3-DPACP-133}, was fixed by means of 2269 ICRF3 S/X sources, for which optical counterparts have 
been identified in EDR3 and which have a valid colour information \gacs{nu\_eff\_used\_in\_astrometry}.
In order to correct a specific problem identified with the bright reference frame of EDR3 
(Sect.~\ref{sec:adHoc}) we also make use of the positional reference frame of \textsc{Hipparcos} at 
epoch J1991.25 as defined by the revised \textsc{Hipparcos} catalogue \citepads{2007ASSL..350.....V}.  

\paragraph{Ephemerides} Accurate barycentric ephemerides of \textit{Gaia} and of all the major bodies 
in the solar system, as well as for some moons and minor planets, are needed in order to interpret the
directions observed by \textit{Gaia} in terms of astrometric parameters defined in the 
barycentric system. The solar system ephemeris used for EDR3 
is the INPOP10e provided by the IMCCE \citepads{2016NSTIM.104.....F}. 
The orbit of \textit{Gaia} was determined at the Mission Operations Centre (MOC)
located at ESOC (Darmstadt, Germany), using conventional Doppler and range tracking 
as well as Delta-Differential One-way Range (Delta-DOR) measurements, the latter
using two tracking stations and calibrated by simultaneous observations of a quasar 
with known position.

The elementary along-scan (AL) astrometric observation is the precise time, 
$t_\text{obs}$, when the centre of a
stellar image crosses the calibrated fiducial line on the CCD. This time is initially given as an on-board
time (OBT), that is the number of nanoseconds counted by the on-board rubidium clock from an arbitrary 
origin, but must be transformed to the coordinate time (TCB) of the event before it can be used in the
astrometric solution. This transformation, known as the time ephemeris, is derived from an analysis of 
time couples (the OBT of a signal generated on board and the reading of the ground-station clocks
when it was received at the ground station), using a sophisticated model that takes into account \textit{Gaia}'s position relative to the Earth, Earth orientation parameters, relativistic effects in the signal propagation,
the influence of the Earth's troposphere, differences between the ground-station clocks and UTC,  
etc.\ \citepads{2017SSRv..212.1423K}. 

\paragraph{Basic-angle corrector} 
The basic angle monitor (BAM) is an interferometric device measuring short-term
($\lesssim 1$~day) variations of the basic angle at $\mu$as precision 
\citepads{2016SPIE.9904E..2DM}. BAM measurements are available since before the start of
nominal operations and throughout the entire period of observations used for EDR3. They
were processed off-line using the same methods as for DR2 
(Sect.~2.4 in \citeads{2018A&A...616A...2L}), resulting in a table of basic-angle
jumps (with the estimated time and amplitude of each jump) and, in between the jumps, 
a continuous function of time represented by a spline. 
The jumps and spline together define the function $\Delta\Gamma(t)$ in Eq.~(\ref{eq:calAL}).

\section{Models}
\label{sec:models}

\subsection{Source model}
\label{sec:source}

The astrometric processing for \textit{Gaia} EDR3 is based on a consistent 
theory of relativistic astronomical reference systems \citepads{2003AJ....126.2687S}. 
The primary coordinate system is the Barycentric Celestial Reference System 
(BCRS) with origin at the solar system barycentre and axes aligned with the 
International Celestial Reference System (ICRS). The time-like coordinate of 
the BCRS is the barycentric coordinate time (TCB). The
\textit{Gaia} relativity model (\citeads{2003AJ....125.1580K}; \citeads{2004PhRvD..69l4001K})
provides a rigorous general-relativistic modelling of astrometric observations.

For the purpose of deriving the main astrometric results in EDR3, it is assumed that
all sources outside of the solar system move with uniform velocity relative to the solar 
system barycentre. Thus, non-linear motions caused by binarity and other perturbations 
are presently ignored, but will be taken into account in future \textit{Gaia} releases.
In the present model, which we refer to as the standard model of stellar motion 
(\citeads{1997ESASP1200.....E}; \citeads{2020A&A...633A...1L}), the motion of the source is 
completely specified by six kinematic
parameters, conventionally taken to be the standard five astrometric parameters
($\alpha$, $\delta$, $\varpi$, $\mu_{\alpha*}{\,=\,}\mu_\alpha\cos\delta$, $\mu_\delta$)
and the radial velocity ($v_r$). All parameters refer to the adopted reference epoch, 
which for the EDR3 astrometry is
$\text{J2016.0}=\text{JD~2457\,389.0 (TCB)}=\text{1 January 2016, 12:00:00}$ (TCB).
This is exactly 0.5~Julian year (182.625~days) later than the reference epoch 
J2015.5 adopted for \textit{Gaia} DR2.

In spite of the well-known fact that a large fraction of the stars in the solar neighbourhood are members
of double and multiple systems, the standard model of stellar motion is very often a good 
model for the observed motions of stars in our Galaxy, and further away, at least over the
relatively short time span covered  by \textit{Gaia}'s observations. In practice, only 
$\sim$10\% of the stars may have proper motions that are noticeably non-linear over 
a few years (cf.\ \citeads{2005ESASP.576...97S}). One reason for this is the extremely
wide range of periods in physical systems, which means that most of them either have
too long periods to show significant curvature over a short time, or they are so close and 
have such short periods that their photocentric wiggles are small and average out over a 
few years. The standard model is also very often an excellent approximation for extragalactic
sources such as active galactic nuclei (AGNs) or quasars. The astrometric solution for 
\textit{Gaia} relies heavily on the lucky circumstance that the motions of most point-like 
sources in the sky can be accurately represented by this simple model. 

The standard model takes into account perspective acceleration through terms depending 
on the radial velocity $v_r$. In \textit{Gaia} DR2 this effect was only considered for some 
50 nearby \textsc{Hipparcos} sources; for EDR3 it is taken into account whenever possible, 
using radial-velocity data from \textit{Gaia}'s radial-velocity spectrometer 
(RVS; \citealt{2018A&A...616A...6S}) as provided in \textit{Gaia} DR2. For a small number of nearby 
stars (mainly white dwarfs), this was complemented with radial velocities from the literature.
Apart from the change in reference epoch and the more frequent use of radial velocity
data, the source model for EDR3 is exactly the same as was used for DR2.

In the standard model, the radial velocity is needed, in addition to the usual five astrometric parameters,
for a complete specification of the six-dimensional phase space vector of a nearby star. 
Because of this, $v_r$ (or $\mu_r= v_r\varpi/A_\text{u}$, where $A_\text{u}$ is the astronomical
unit) is sometimes called the sixth astrometric parameter. This is potentially confusing
in connection with the six-parameter solutions discussed in Sect.~\ref{sec:nuEff} and
elsewhere, where the sixth parameter is the pseudocolour $\hat{\nu}_\text{eff}$, that is 
the astrometrically estimated effective wavenumber (colour) of the source. In contrast
to the pseudocolour, the radial velocity is never estimated from \textit{Gaia} data
in any of the solutions discussed here, although it will be possible in the future for a small 
number of nearby high-velocity stars \citepads{1999A&A...348.1040D}.

\subsection{Attitude model}
\label{sec:attitude}

The attitude model for \textit{Gaia} EDR3 is the same as was used for DR2,
except that AL observations made in window class WC0b (see Sect.~\ref{sec:cal})
were not used for the attitude determination. The attitude model includes a 
pre-computed AL corrective attitude that removes much of the rapid attitude 
irregularities created by micro-clanks and high-frequency thruster noise. 
We refer to Sect.~3.2 of \citetads{2018A&A...616A...2L} for a description
of the DR2 model.


\begin{figure*}
\center
  \includegraphics[width=16cm]{fig03.pdf} 
\caption{Layout of CCDs in \textit{Gaia}'s focal plane. Star images move 
from right to left in the diagram, as indicated in the lower part of the drawing by the 
nominal paths of two images, one in the preceding FoV (PFoV) and one in the 
following FoV (FFoV). The along-scan (AL) and across-scan (AC) directions are 
indicated in the top left corner. To the right, one of the CCDs is shown magnified, with 
the fiducial observation lines indicated for selected gates ($g$). Also indicated is the
AC pixel coordinate $\mu$, running from 13.5 to 1979.5 across the image area of
each CCD. The skymappers (SM1, SM2) provide source image detection and FoV 
discrimination, but their measurements are not used in the astrometric solution. The 
astrometric field (AF1--AF9) provides accurate AL measurements and (for two-dimensional 
windows) AC positions. Other CCDs are used for the blue and red photometers 
(BP, RP), the radial-velocity spectrometer (RVS), wavefront sensing (WFS), and
basic-angle monitoring (BAM). One of the CCD strips (AF3) illustrates the system for 
labelling individual CCDs by strip and row index. The origin of the field angles 
$(\eta,\zeta)$ is at different physical locations on the CCDs in the two fields.
(Adapted from \citeads{2012A&A...538A..78L}.)\label{fig:fpa}}
\end{figure*}

 
\begin{table*}[t]
\caption{Summary of the astrometric calibration model and number of calibration
parameters in the astrometric solution for \textit{Gaia} EDR3.\label{tab:cal}}
\small
\begin{tabular}{rllcccccccccccccr}
\hline\hline
\noalign{\smallskip}
&& \multicolumn{1}{l}{Basis functions} 
& \multicolumn{13}{c}{---------------------- Multiplicity of dependencies ----------------------} & 
\multicolumn{1}{c}{Number of} \\ 
\multicolumn{2}{l}{Effect ($i$) and brief description} & 
\multicolumn{1}{l}{$K_{lm}(\tilde{\mu},\tilde{t})$} & \multicolumn{1}{c}{$K_{lm}$} 
& \multicolumn{1}{c}{$j$} & $f$ & \multicolumn{1}{c}{$n$} & 
$g$ & $b$ & $w$ & \multicolumn{1}{c}{$\nu_\text{eff}$} & $G$ & $S$ & $\phi$ 
& \multicolumn{1}{c}{$\Delta t$} & $\dot{\zeta}$ &
\multicolumn{1}{c}{parameters} \\
\noalign{\smallskip}
\hline
\noalign{\smallskip}
1& AL large-scale geometric & $lm=00,10,20,01$ & 4 & 310 & 2 & 62 & -- & -- & 4 & -- & -- & -- & -- & -- & -- &  615\,040\\ 
2& AL medium-scale gate & $lm=00,10$ & 2 & 19 & 2 & 62 & 8 & 9 & -- & -- & -- & -- & -- & -- & -- & 339\,264 \\ 
3& AL large-scale colour & $lm=00,01$ & 2 & 19 & 2 & 62 & -- & -- & 4 & 1 & -- & -- & -- & -- & -- & 18\,848 \\ 
4& AL large-scale saturation & $lm=00$ & 1 & 19 & 2 & 62 & -- & -- & 2 & -- & -- & 1 & -- & -- & -- & 4\,712 \\ 
5& AL large-scale subpixel & $lm=00$ & 1 & 19 & 2 & 62 & -- & -- & 4 & -- & -- & -- & 2 & -- & -- & 18\,848 \\ 
6& AL large-scale CTI & $lm=00$ & 1 & 19 & 2 & 62 & -- & -- & 4 & -- & -- & -- & -- & 4 & -- &  37\,696\\ 
7& AL large-scale AC rate & $lm=00$ & 1 & 19 & 2 & 62 & -- & -- & 1 & -- & -- & -- & -- & -- & 1 &  2\,356\\ 
\noalign{\smallskip}
\hline
\noalign{\smallskip}
8& AC large-scale geometric & $lm=00,10,20,01$ & 4 & 19 & 2 & 62 & -- & -- & 2 & -- & -- & -- & -- & -- & -- &  18\,848\\ 
9& AC large-scale gate & $lm=00$ & 1 & 19 & 2 & 62 & 8 & -- & -- & -- & -- & -- & -- & -- & -- & 18\,848 \\ 
10& AC large-scale colour & $lm=00$ & 1 & 19 & 2 & 62 & -- & -- & 2 & 1 & -- & -- & -- & -- & -- & 4\,712 \\ 
11& AC large-scale magnitude & $lm=00$ & 1 & 19 & 2 & 62 & -- & -- & 2 & -- & 1 & -- & -- & -- & -- & 4\,712 \\ 
12& AC large-scale saturation & $lm=00$ & 1 & 19 & 2 & 62 & -- & -- & 2 & -- & -- & 1 & -- & -- & -- & 4\,712 \\ 
\noalign{\smallskip}
\hline
\end{tabular}
\tablefoot{The column Basis functions lists the combinations of indices $l$ and $m$
used to model variations with AC coordinate on a CCD ($\tilde{\mu}$) and with time within 
a time granule ($\tilde{t}$). Multiplicity of dependencies gives the number of distinct 
functions or values for each dependency, or a dash if there is no dependency: 
basis functions ($K_{lm}$, Eq.~\ref{eq:Klm}), granule index ($j$), 
field index ($f$), CCD index ($n$), gate ($g$), stitch block ($b$), window class ($w$), 
effective wavenumber ($\nu_\text{eff}$), magnitude ($G$), saturation ($S$), 
subpixel phase ($\phi$), time since last charge
injection ($\Delta t$), and AC scan rate ($\dot{\zeta}$). The last column 
is the product of multiplicities, equal to the number of calibration parameters for the effect.
The SMs are not considered here.}
\end{table*}

\subsection{Calibration model}
\label{sec:cal}

The astrometric calibration model for \textit{Gaia} EDR3 is similar to the one
used for DR2, as described in Sect.~3.3 of \citetads{2018A&A...616A...2L}, but
with additional dependencies described below. 
The general principles of the calibration model are described in Sect.~3.4 of 
\citetads{2012A&A...538A..78L}, and only a few basic concepts are recalled here.
At any time, the attitude represents a solid-body rotation from the celestial reference
system to \textit{Gaia}'s scanning reference system (SRS), nominally fixed with respect
to the CCDs as viewed through the two FoVs (preceding and following).
Within a FoV, directions with respect to the SRS are usually expressed by means of the
field angles $(\eta,\zeta)$, with origin at the nominal centre of the FoV (Fig.~\ref{fig:fpa}). 
According to the scanning law, stellar images traverse the FoV in the direction of decreasing 
$\eta$ (at the AL rate of approximately 60~arcsec~s$^{-1}$) and at approximately constant 
$\zeta$ (the AC rate is at most $\pm 0.18$~arcsec~s$^{-1}$). The fundamental AL measurement 
used for the astrometry is the precise time when an image transits across a fiducial 
`observation line' line fixed to the CCD (Fig.~\ref{fig:fpa}, right). The astrometric calibration 
of the instrument (as opposed to the LSF and PSF calibrations by CAL) is essentially a specification 
of the location of the observation line in field angles, that is of the functions $\eta(\mu)$ and 
$\zeta(\mu)$, where $\mu$ is the AC pixel coordinate.  

More precisely, the AL and AC calibration functions are written as the 
sums of the nominal calibrations and several `effects' that describe the dependence on 
various quantities, such as time, CCD, and FoV (see Eqs.~\ref{eq:calAL} and \ref{eq:calAC}). 
The effects used in the EDR3 calibration model in the AF
are listed in Table~\ref{tab:cal}.
The skymappers (SM1 and SM2 in Fig.~\ref{fig:fpa}) obtain a similar, but simpler, 
calibration. However, the SM observations are not at all used in the astrometric solution, 
and their calibration is not discussed in this paper.

Compared with the corresponding table for the 
DR2 model (Table~2 in \citeads{2018A&A...616A...2L}), Table~\ref{tab:cal} contains effects with
several new dependencies ($S$, $\phi$, $\Delta t$, $\dot{\zeta}$). Their introduction
in the model was motivated by systematic trends seen in the residuals from preliminary 
solutions, in which the calibration model did not include the effects. 
The complete set of dependencies is as follows.
\begin{itemize}
\item
AC pixel coordinate $\mu$ on the CCD, which is a continuous value running
from 13.5 to 1979.5 across the AC extent of the CCD image area (Fig.~\ref{fig:fpa}).
The offset by 13~pixels allows for the presence of pre-scan pixel data.
\item
Time $t$, divided into granules such that $t_j\le t<t_{j+1}$ in the granule indexed 
by $j$. Two different time axes are used, with 310 and 19 granules spanning the length 
of the data; the typical duration of the granules is, respectively, about 3~d and 63~d. 
\item 
FoV index $f$, specifying preceding or following FoV. We use the
convention $f=+1$ in the PFoV and $f=-1$ in the FFoV.
\item
CCD $n$, with 62 different values in the astrometric field.
\item
Gate $g$, taking eight different values with $g=0$ for ungated observations
(Fig.~\ref{fig:obsFreqVsWcGate}). The number of active TDI lines is 4500 for $g=0$,
2900 ($g=12$), 2048 ($g=11$), 1024 ($g=10$), 512 ($g=9$), 256 ($g=8$), 
128 ($g=7$), and 16 ($g=4$). Gates 1--3, 5, and 6 are not used in normal operations.
\item
Stitch block $b$, with nine different values in the AC direction of a CCD. $b$ is uniquely
defined by the AC pixel coordinate through $b=\lfloor (\mu+128.5)/250\rfloor$,
where $\lfloor\,\rfloor$ is the floor function.
\item
Window class $w$, with four values. In the DR2 calibration model, three window classes 
WC0, WC1, and WC2 were used, approximately corresponding to magnitude ranges
$G\lesssim 13$, $13\lesssim G\lesssim 16$, and $16\lesssim G$, respectively. (The
WC represents the CCD sampling scheme chosen at detection time, depending 
mainly on the real-time estimate of the magnitude derived from the SM observation, 
but also on several other factors. There is consequently no strict relation between the 
mean calibrated $G$ magnitude, given in the catalogue, and the WC.) 
In the EDR3 model, WC0 was further 
subdivided into WC0a  (for $G\lesssim 11$) and WC0b (for $11\lesssim G\lesssim 13$),
see Fig.~\ref{fig:obsFreqVsWcGate}.
\item
Effective wavenumber $\nu_\text{eff}$ is the photon-weighted inverse wavelength, 
calculated from the BP and RP spectra in the photometric processing \citep{EDR3-DPACP-118} 
and expressed in $\mu\text{m}^{-1}$. The cyclic processing scheme 
adopted by DPAC implies that the
$\nu_\text{eff}$ used for the EDR3 astrometry was generated in the preceding cycle,
corresponding to DR2 photometry, and is sometimes missing or inconsistent with the 
EDR3 photometry. The actual values used for the astrometry (and IDU pre-processing)
is given in the \textit{Gaia} Archive as \gacs{nu\_eff\_used\_in\_astrometry}. For sources
without a reliable $\nu_\text{eff}$ a special calibration was employed (step~\ref{step:Cprime}
in Sect.~\ref{sec:steps}).   
\item
Magnitude $G$: Like the effective wavenumber, the magnitude used in the astrometric
processing was derived from DR2, but since the differences are generally small and only the
(less critical) AC calibration depends on $G$, the actual value used is not given in the 
Archive.  
\item
Saturation $S$: This is a flag produced by the IPD as part of the IDU pre-processing. 
It is set to 1 if the raw observed sample exceeds 
a pre-defined conservative threshold, as determined from early mission data, for the 
CCD column and sample binning; otherwise $S=0$. 
The astrometric effects of the saturation are only calibrated for WC0a and WC0b.
\item
Subpixel phase $\phi$: This is $2\pi$ times the fractional part of the precise observation 
time $t$, as determined by the IPD and expressed in TDI periods of on-board time. (The TDI 
period is the time it takes to shift the charges on the CCDs by one pixel AL, or approximately
0.982~ms.) Inaccuracies in the LSF and PSF calibrations used for the IPD may result in
systematic AL errors that are periodic functions of $\phi$.
\item
Time since the last charge injection $\Delta t$: To minimise the effects of charge transfer
inefficiency (CTI) in the CCDs, charge injections are made at regular time intervals
of 2000~TDI periods. CTI may cause systematic AL shifts of the image centroids, which
increase with $\Delta t$.    
\item
Across-scan (AC) rate $\dot{\zeta}$: The nominal scanning law of \textit{Gaia} 
\citepads{2016A&A...595A...1G} produces a quasi-periodic ($\simeq 6$~h period) variation 
of the AC rate, with an amplitude of approximately 173~mas~s$^{-1}$. Imperfections 
in the PSF modelling may result in systematic AL errors that depend on the AC rate.
This dependence is only calibrated for observations in WC0b using gates 11, 12, and 0
(that is, for a CCD exposure time of about 2.0, 2.8, or 4.4~s).
\end{itemize}
Within a time granule, the variation with $t$ and $\mu$ is modelled as a 
linear combination of basis functions
\begin{equation}\label{eq:Klm}
K_{lm}(\tilde{\mu},\tilde{t})=\tilde{P}_l(\tilde{\mu})\tilde{P}_m(\tilde{t}\,) \, ,
\end{equation}
where $\tilde{P}_l(x)$, $\tilde{P}_m(x)$ are the shifted Legendre polynomials%
\footnote{The shifted Legendre polynomials $\tilde{P}_n(x)$ are related to the 
(ordinary) Legendre polynomials $P_n(x)$ by $\tilde{P}_n(x)=P_n(2x-1)$.
Specifically, $\tilde{P}_0(x)=1$, $\tilde{P}_1(x)=2x-1$, and 
$\tilde{P}_2(x)=6x^2-6x+1$. The shifted Legendre polynomials are orthogonal 
on $0\le x\le 1$.\label{fn2}}
of degree $l$ and $m$,
$\tilde{\mu}=(\mu-13.5)/1966$ is the normalised AC pixel coordinate, and 
$\tilde{t}=(t-t_j)/(t_{j+1}-t_j)$ the normalised time within granule $j$.
The third and fourth columns in Table~\ref{tab:cal} list the combination of indices 
$l$ and $m$ used for a particular effect and the number of basis functions 
$K_{lm}$. Most of the effects only
use $lm=00$, meaning that the effect is modelled as constant with $t$ and $\mu$
for a given combination of the other indices and variables.

Each combination of indices $l$, $m$, $j$, $f$, $n$, $g$, $b$, and $w$ indicated in 
Table~\ref{tab:cal} is a `calibration unit' and receives an independent calibration.
Within a calibration unit, effect $i$ is a linear combination of products 
$\smash{K_{lm}(\tilde{\mu},\tilde{t})\,\Psi^{(i)}_k(x)}$, where $K_{lm}$ describes the dependence
on $t$ and $\mu$ according to Eq.~(\ref{eq:Klm}), and $\smash{\Psi^{(i)}_k(x)}$ ($k=0,~1,~\dots$) 
describe the dependence on some other variable $x$, which could be
$\nu_\text{eff}$, $G$, $S$, $\phi$, $\Delta t$, or $\dot{\zeta}$. The relevant
functions are:
\begin{gather}
\label{eq:PhiNuEff}
\Psi^{(3,\,10)}_0(\nu_\text{eff})=\nu_\text{eff}-1.43~\mu\text{m}^{-1} \, , \\
\label{eq:PhiG}
\Psi^{(11)}_0(G)=G-12.6\, , \\
\label{eq:PhiS}
\Psi^{(4,\,12)}_0(S)=S\, , \\ 
\label{eq:PhiPhi}
\Psi^{(5)}_0(\phi)=\cos\phi\, ,\quad \Psi^{(5)}_1(\phi)=\sin\phi\, , \\
\label{eq:PhiDeltaT}
\Psi^{(6)}_k(\Delta t) = a_k - \exp(-\Delta t/\tau_k) \, ,\quad k=0\dots 3\, , \\
\label{eq:PhiDotZeta}
\Psi^{(7)}_0(\dot{\zeta}) = |\,\dot{\zeta}\,|^2\, .
\end{gather}
The function $\Psi^{(11)}_0(G)$ in Eq.~(\ref{eq:PhiG}) is only used for ungated 
AC observations; otherwise, it is set to 1.
Equation~(\ref{eq:PhiPhi}) describes a periodic variation with subpixel phase $\phi$.
In Eq.~(\ref{eq:PhiDeltaT}), the variation with time since the last charge injection 
$\Delta t$ is assumed to be a linear combination of four exponentials,
with e-folding times $\tau_k=10$, 100, 500, and 2000~TDI periods. The constants 
$a_k=(\tau_k/2000)[1-\exp(-2000/\tau_k)]$ are such that the mean
value of $\smash{\Psi^{(6)}_k(\Delta t)}$ over $0\le\Delta t\le 2000$ is zero. 
This means that the mean displacement of the images caused by the CTI is 
not taken out by this calibration, only its variation with $\Delta t$. The resulting 
calibration parameters are thus mainly interesting as diagnostics of the effect 
(Fig.~\ref{fig:cti}). The quadratic dependence on the AC rate $\dot{\zeta}$ 
in Eq.~(\ref{eq:PhiDotZeta}) models a possible bias caused by the AC smearing of 
the PSF; this effect was not well modelled in the PSF calibration for EDR3 
(cf.\ Sect.~\ref{sec:acRate} and Appendix~\ref{sec:plxAc}). 
Formally, $\Psi^{(1,2,8,9)}=1$ for the effects that only depend on $t$ and $\mu$.

The complete AL calibration model is 
\begin{equation}\label{eq:calAL}
\begin{split}
\eta(\mu,t,\nu_\text{eff},\dots) &= \eta_{ng}^{(0)}(\mu) + \frac{f}{2}\Delta\Gamma(t)\\
&+ \sum_{i=1}^7 \sum_{l,\,m,\,k} \Delta\eta^{(i)}_{lmk} K_{lm}\,\Psi^{(i)}_k + \delta\eta(f,t,\eta,\zeta) \, ,
\end{split}
\end{equation}
where the first term is the nominal location of the fiducial observation line for CCD $n$,
gate $g$ 
(Eq.~14 in \citeads{2012A&A...538A..78L}); the second contains the basic angle correction
$\Delta\Gamma(t)$ derived from BAM data (Sect.~\ref{sec:aux}); the third is the sum of the 
seven effects $i$ in the upper part of Table~\ref{tab:cal}, with calibration parameters 
$\Delta\eta^{(i)}_{lmk}$; and the last term is the spin-related distortion model fitted as 
global parameters (Sect.~\ref{sec:glob}). 
For brevity, the dependences on $f$, $n$, $g$, $b$, and $w$ and the arguments of 
$K_{lm}$ and $\smash{\Psi^{(i)}_k}$ have been suppressed.
This gives a total of 1\,036\,764 AL parameters (not counting the spin-related distortion
parameters), which is more than three times as many as used for the DR2 calibration model
(see Table~1 in \citeads{2018A&A...616A...2L}). Besides the longer time interval 
covered by the data, this reflects the more complex modelling made necessary (and
possible) thanks to the generally improved quality of the input data and resulting
solution. Several new effects have been introduced (saturation, subpixel, CTI, and 
AC rate), and the AL large-scale geometric calibration now depends also on the 
window class. Some of these effects should eventually be taken out by the LSF and PSF 
calibrations, but that was not yet possible in the present cycle.



\begin{figure}
\centering
  \includegraphics[width=8cm]{fig04.pdf}
  \caption{Relative frequency of observations in the various
  combinations of window class and gate, as a function of magnitude.
  The four blocks represent the four window classes (WC); within each WC
  the eight stripes represent (from top to bottom) gate number 4, 7, 8, 9, 10, 
  11, 12, and 0. The graph was constructed from a random 1\% sample of
  the AF observations of the primary sources. The faint sources observed in 
  WC0a at gate~0 are the Calibration Faint Stars, a small fraction of faint 
  observations receiving full-pixel resolution windows for calibration purposes
  \citepads{2016A&A...595A...1G}.}
  \label{fig:obsFreqVsWcGate}
\end{figure}


The AC calibration model is 
\begin{equation}\label{eq:calAC}
\begin{split}
\zeta(\mu,t,\nu_\text{eff},\dots) = \zeta_{fng}^{(0)}(\mu) 
&+ \sum_{i=8}^{12} \sum_{l,\,m,\,k} \Delta\zeta^{(i)}_{lmk} K_{lm}\,\Psi^{(i)}_k \\
\Bigl[~&+~\delta\zeta(f,t,\eta,\zeta)~\Bigr] \, ,
\end{split}
\end{equation}
where $\zeta^{(0)}_{fng}(\mu)$ is the nominal calibration and $\Delta\zeta^{(i)}_{lmk}$ 
the calibration parameter for the five effects in the bottom part of Table~\ref{tab:cal}.
The last term is the spin-related distortion in AC; but as explained in Sect.~\ref{sec:glob}
this is not used for EDR3, and is therefore put within square brackets in Eq.~(\ref{eq:calAC}).  
This gives a total of 51\,832 AC calibration parameters, which is 10\% smaller than 
in DR2, in spite of the longer time period covered. The main reason for this is that the 
dependences on time and AC coordinate were found to be overly complicated in DR2 
and have been simplified.

Selected results of the astrometric calibration are given in Appendix~\ref{sec:calData}.

\subsection{Spin-related distortion model}
\label{sec:glob}

Thermo-mechanical perturbations of the instrument over time
scales close to and below the rotational period of 6~h
present a special problem for the AGIS calibrations. Such variations
will be called `quick' below. It is known
that it is impossible to fully calibrate the quick variations of the
instrument \citepads{2017A&A...603A..45B}. The reason for this is 
a degeneracy between the source parameters, attitude parameters, 
and calibration parameters for the quick variations of the instrument. 
This is known as the VBAC degeneracy (Velocity error and
effective Basic Angle Calibration).%
\footnote{The mention of velocity here may seem surprising. For an astrometric 
satellite such as \textit{Gaia}, the observational effects of a 
small error in the translational velocity, as used in the modelling of stellar aberration, 
is found to be indistinguishable from a certain combination of errors in the attitude and in the 
basic angle. This aspect of the data processing plays no role for EDR3, where velocity 
is taken to be known (Sect.~\ref{sec:aux}).} 
A special case of this is the well-known degeneracy between the global parallax zero 
point and a specific form of basic-angle and attitude variations \citepads{2017A&A...603A..45B},
but the VBAC degeneracy is much more general: Any time-dependent distortion of 
the celestial positions is observationally indistinguishable from some specific 
combination of attitude errors and quick instrument variations. In the design 
considerations for \textit{Gaia} it was indeed a fundamental requirement that 
the instrument must either be extremely stable on time scales shorter than a few times 
the spin period, or have the means to monitor the variations continuously to a very high 
precision.

However, the VBAC degeneracy does not imply that arbitrary quick instrument variations
are degenerate with the source and attitude parameters. On the contrary, most such variations 
are not degenerate, and therefore in principle possible to calibrate from the astrometric 
observations themselves, that is, without the need for special procedures or metrology 
devices like the BAM. Over the last decade, 
a considerable effort has been put into investigating how, and to what extent, quick 
variations of the \textit{Gaia} instrument can be calibrated from the astrometric observations.
The spin-related distortion model presented here is a limited version of more general 
models that may be used for future \textit{Gaia} releases. It nevertheless represents 
a significant advance over the DR2 model (Sect.~3.4 in \citeads{2018A&A...616A...2L}).
Although the distortion model logically belongs to the instrument calibration model,
it is fitted as part of the global block in AGIS for purely implementation-technical reasons. 

\subsubsection{General model}
\label{sec:vbacfoc}

We begin by formulating the spin-related distortion model in its most general form.
At any moment of time an arbitrary distortion of the AL and AC
field angles $(\eta,\zeta)$ for a source in a given FoV $f$
can be represented as a two-dimensional expansion over a family of
orthogonal functions $\Phi_{lm}(\eta,\zeta)$:
\begin{eqnarray}\label{expansion-general-eta}
\delta\eta(f,t,\eta,\zeta)
=\sum\limits_{l\ge0}\sum\limits_{m\ge0}\;\delta\eta_f^{lm}(t)\;\Phi_{lm}(\eta,\zeta)\,,\\
\label{expansion-general-zeta}
\delta\zeta(f,t,\eta,\zeta)
=\sum\limits_{l\ge0}\sum\limits_{m\ge0}\;\delta\zeta_f^{lm}(t)\;\Phi_{lm}(\eta,\zeta)\,.
\end{eqnarray}
\noindent
For the nominally rectangular AF of \textit{Gaia}, a convenient set of orthogonal 
functions are the products of the Legendre polynomials $P_k(x)$ for the AL and AC coordinates: 
\begin{equation}\label{Phi-lm-Gaia}
\Phi_{lm}(\eta,\zeta)=P_l(\tilde\eta)\,P_m(\tilde\zeta)\,,
\end{equation}
where $\tilde\eta$ and $\tilde\zeta$ are the field angles $\eta$
and $\zeta$ linearly scaled to the interval $[-1,1]$. 
(Although not apparent in these equations, the scaling of $\zeta$ is actually different in the
two FoVs, owing to the offset of the AC origins indicated in Fig.~\ref{fig:fpa}.)

We use Eqs.~(\ref{expansion-general-eta}) and (\ref{expansion-general-zeta}) to model the
spin-related distortion terms in Eqs.~(\ref{eq:calAL}) and (\ref{eq:calAC}). Their time-variations
are defined by the functions $\smash{\delta\eta_f^{lm}(t)}$ and $\smash{\delta\zeta_f^{lm}(t)}$
to be specified below. 

Full-scale simulations of the AGIS solution have shown that the terms in Eqs.~(\ref{expansion-general-eta}) 
and (\ref{expansion-general-zeta}) with $l+m\ge1$ are not degenerate with the source parameters and attitude, 
and can therefore safely be determined from the observations. For example, simulated variations 
in which $\smash{\delta\eta_f^{lm}(t)}$ and $\smash{\delta\zeta_f^{lm}(t)}$ for $1\le l+m\le 5$ were 
represented as B-splines, with a knot interval of 10~min and random coefficients, could be completely 
recovered with no rank-deficiency and only a moderate slowdown of the convergence of the 
iterative solution. The terms with $l=m=0$, on the other hand, involve the VBAC degeneracy,
and the time variation of these terms need to be chosen with care in order not to jeopardise 
the source parameters. It is therefore natural to split the further specification of the general 
spin-related distortion model in two parts, corresponding to VBAC (for $l=m=0$) and 
FOC (for $l+m\ge 1$). Here, FOC stands for Focal length and Optical distortion Calibration, 
since the calibration for $l+m\ge 1$ obviously covers also a variation of the focal length 
of the instrument.

The general model in Eqs.~(\ref{expansion-general-eta}) and (\ref{expansion-general-zeta})
can also describe slow variations of the instrument and is in principle degenerate with certain 
parts of the calibration model described in Sect.~\ref{sec:cal}. As this degeneracy does 
not involve the source parameters it is harmless for the astrometry, but since it could slow down 
the convergence of the AGIS iterations it should nevertheless be avoided if practically feasible. 
One such case is the near-degeneracy mentioned below between effect~7 in Table~\ref{tab:cal} 
and the Fourier terms of order $p=2$ in Eqs.~(\ref{eq:calVF3b}) and (\ref{eq:calVF3c}).


\subsubsection{FOC}
\label{sec:foc}

Because the FOC calibration ($l+m\ge 1$) has no degeneracy with the source and attitude 
parameters, we are quite free to choose the maximum degree $l+m$ of the basis functions
$\Phi_{lm}(\eta,\zeta)$ and the parametrisation of the time dependency of the functions
$\smash{\delta\eta_f^{lm}(t)}$ and $\smash{\delta\zeta_f^{lm}(t)}$. Of course, it is always
necessary to limit the degree and time resolution so that the number of parameters is reasonable
in relation to the number of observations, keeping the overall solution numerically and 
practically tractable.
For example, because the AF has nine strips of CCDs (see Fig.~\ref{fig:fpa}),
it is not numerically feasible to have the AL degree $l>9$. 
Numerous test solutions using EDR3 data were made to explore some 
of the many possible options, and the configuration finally adopted for AGIS~3.2 is in some 
sense the best one found in the limited time available. 

One conclusion from the test solutions was that 
the FOC correction in AC does not bring any improvement at this stage and it is therefore not
used for EDR3; hence the bracketed term in Eq.~(\ref{eq:calAC}). It was also found that the
polynomials in Eq.~(\ref{expansion-general-eta}) can be restricted to $1\le l+m\le 3$, which
gives 18 coefficients $\smash{\delta\eta_f^{lm}(t)}$ to be considered (nine per FoV). Similarly to the AL
calibration model in Sect.~\ref{sec:cal}, it was found necessary to determine FOC separately
for each window class. For WC0a, WC1, and WC2, the ten coefficients with $1\le l+m\le 2$ 
were determined as cubic splines with a knot interval of 20~min. The remaining eight 
coefficients with $l+m=3$ were fitted as Fourier polynomials
\begin{equation}\label{eq:calVF3b}
\delta\eta_f^{lm}(t) = d(t)^{-2} \sum_{p=1}^8 \Bigl(c_{f\!plm}\cos p\Omega(t) 
+ s_{f\!plm}\sin p\Omega(t)\Bigr)\, ,
\end{equation}
where $d(t)$ is the Sun--\textit{Gaia} distance in au, $\Omega(t)$ is
the heliotropic spin phase \citepads{2018A&A...616A...2L}, and 
$c_{f\!plm}$ and $s_{f\!plm}$ are the free parameters fitted to the data.
The scaling by $d(t)^{-2}$ here and in the following equations accounts for the 
variation in solar irradiance.
For WC0b, all 18 coefficients $\smash{\eta_f^{lm}(t)}$ were fitted as
Fourier polynomials, as in Eq.~(\ref{eq:calVF3b}), but omitting the terms with
$p=2$ to avoid the near-degeneracy with effect~7 in Table~\ref{tab:cal}.
Furthermore, all Fourier polynomials for FOC were fitted independently for
the two time intervals before and after OBMT~4513~rev.  That moment of
time (one of the gaps in Table~\ref{tab:gaps}) was found to be a boundary 
between slightly different behaviours of the residuals in test solutions; 
ultimately, this behaviour can be traced back to a particular change of 
LSF and PSF models in CALIPD 3.1 at that moment.
The resulting FOC model has a total of 2\,033\,184 parameters.


\subsubsection{VBAC}
\label{sec:vbac}

The terms in Eqs.~(\ref{expansion-general-eta}) and (\ref{expansion-general-zeta}) with 
$l=m=0$ represent the distortion averaged over each FoV. It is readily seen that this is 
equivalent to a combination of four time dependent variations, namely,
(i) of the AL attitude by 
$\smash{\frac{1}{2}\bigl(\delta\eta_{+1}^{0,0}+\delta\eta_{-1}^{0,0}\bigr)}$;
(ii) of the AC attitude in the PFoV by
$\smash{\delta\zeta_{+1}^{0,0}}$;
(iii) of the AC attitude in the FFoV by
$\smash{\delta\zeta_{-1}^{0,0}}$; and
(iv) of the basic angle by
$\smash{\delta\Gamma=\delta\eta_{+1}^{0,0}-\delta\eta_{-1}^{0,0}}$.
The flexibility of the attitude modelling means that the first three variations are completely
degenerate with the attitude determination, and should not be further considered in the 
VBAC model. Therefore, the only variation to consider for $l=m=0$ is a time-dependent 
basic angle variation, $\delta\Gamma(t)$. 

$\delta\Gamma(t)$ can be regarded as an additive correction to the basic angle variation 
$\Delta\Gamma(t)$ in Eq.~(\ref{eq:calAL}) that comes from the analysis of BAM data 
(Sect. \ref{sec:aux}). 
It should be recalled that $\Delta\Gamma(t)$ includes both basic angle jumps (due to
sudden structural changes in the optics) and a smooth representation of the basic angle 
variations between jumps, including a very good approximation of the quick variations.
However, because the CCDs for the BAM are located outside of the AF
(Fig.~\ref{fig:fpa}), we cannot assume that the variations measured with the BAM are 
fully representative for the whole FoV -- indeed, in the presence of FOC distortion this
is not to be expected. Moreover, the BAM device itself may be subject to perturbations that
are not relevant for the astrometric observations. For these reasons it is highly desirable
to calibrate as much as possible of the basic angle variations directly from the astrometric
observations, which can be done with VBAC. Owing to the VBAC degeneracy there are nevertheless 
components of the basic angle variations that cannot be determined from the observations,
and the BAM signal remains indispensable as the only handle we may have on those components.

For EDR3, the same representation of $\delta\Gamma(t)$ was 
used as for DR2 (Eq.~10 in \citeads{2018A&A...616A...2L}), but split in two parts, 
\begin{equation}\label{eq:deltaGamma-split}
\delta\Gamma(t)=\delta\Gamma_\text{A}(t)+\delta\Gamma_\text{B}(t)\,,
\end{equation}
with
\begin{multline}\label{eq:calVF3c}
\delta\Gamma_\text{A}(t) = d(t)^{-2} 
\sum_{p=2}^8\, \Bigl[\delta C_{p,0}+(t-t_\text{ep})\delta C_{p,1}\Bigr]\cos p\Omega(t) \\
+ d(t)^{-2} 
\sum_{p=1}^8\, \Bigl[\delta S_{p,0}+(t-t_\text{ep})\delta S_{p,1}\Bigr]\sin p\Omega(t) \,  
\end{multline}
and
\begin{equation}\label{eq:deltaGamma-C10}
\delta\Gamma_\text{B}(t) = d(t)^{-2} \Bigl[\delta C_{1,0}+(t-t_\text{ep})\delta C_{1,1}\Bigr]\cos\Omega(t)\,.
\end{equation}
Here $d(t)$ and $\Omega(t)$ have the same meaning as in Eq.~(\ref{eq:calVF3b}),
$t_\text{ep}=\text{J}2016.0$, and $\delta C_{p,q}$, $\delta S_{p,q}$ are the constant 
coefficients determined from the data. The split in Eq.~(\ref{eq:deltaGamma-split}) 
is motivated by the near-degeneracy of $\delta C_{1,0}$ with a global parallax shift
(\citeads{2017A&A...603A..45B}; \citeads{2018A&A...616A...2L}),
which necessitates a special treatment of this term; this is deferred till Sect.~\ref{sec:C10}.
The parameter $\delta C_{1,1}$ is included in $\delta\Gamma_\text{B}(t)$ only because it
naturally belongs together with $\delta C_{1,0}$; it is not strongly correlated with other
parameters and could instead have been fitted with the other VBAC parameters in 
Eq.~(\ref{eq:calVF3c}).

The representation of $\delta\Gamma_\text{A}(t)$ in
Eq.~(\ref{eq:calVF3c}) contains 30 parameters describing linear
variations of the $d(t)^{-2}$-scaled Fourier coefficients in
$\Omega(t)$. Analysis of the test AGIS solutions and their residuals has shown
that the effective basic angle variations obtained with this model are substantially 
different for the different window classes. A separate set of Fourier coefficients was
therefore fitted for each window class.
Moreover, similarly as for the Fourier coefficients
in the FOC model, separate fits were made for the time intervals before and after 
OBMT~4513~rev, and the coefficients for $p=2$ were omitted for WC0b.
The resulting model for $\delta\Gamma_\text{A}(t)$ has a total of 232 parameters.


\subsubsection{Treatment of the near-degeneracy with parallax ($\delta C_{1,0}$)}
\label{sec:C10}

Here we consider the VBAC correction $\delta\Gamma_\text{B}(t)$ in 
Eq.~(\ref{eq:deltaGamma-C10}), containing the two additional parameters
$\delta C_{1,0}$ and $\delta C_{1,1}$. Unlike the parameters in $\delta\Gamma_\text{A}(t)$,
which were fitted per window class and separately before and after OBMT~4513~rev, 
there is only a single set of these two parameters. They are fitted using all data except WC0b.

As already mentioned, the parameter $\delta C_{1,0}$ cannot be easily fitted in an
iterative solution like AGIS because it is highly correlated to a global shift of all parallaxes
\citepads{2017A&A...603A..45B}. However, this also means that if the correction $\Delta\Gamma(t)$ 
to the basic angle derived from BAM data has an error described by $\delta C_{1,0}$, there will be
a global shift of the parallaxes. Owing to the profound scientific importance of the parallax
zero point, every effort should be made to avoid such an error. To this end a
method has been developed to calibrate $\delta C_{1,0}$ directly from the astrometric 
observations of \textit{Gaia}. The method was thoroughly tested in a series of detailed
end-to-end simulations of the iterative solution, which demonstrated the feasibility of the 
method and probed the limits of its applicability. It was tested with cycle~2 data (but 
not used in the solution for \textit{Gaia} DR2), and finally employed in the primary astrometric 
solution for EDR3. Full details of the method will be published elsewhere; here we describe 
only its most important elements.

The possibility to fit $\delta C_{1,0}$ 
is based on the small but not completely negligible differences between
heliocentric and barycentric quantities, and the fact that the solar irradiance and parallax
factor scale differently with the varying distance from the Sun or solar system barycentre.
Based on physical considerations, the model for the basic angle variations in 
Eqs.~(\ref{eq:calVF3c}) and (\ref{eq:deltaGamma-C10}) scales as $d(t)^{-2}$ and is
periodic in $\Omega(t)$, where $d(t)$ and $\Omega(t)$ are heliocentric, that is, reckoned 
with respect to the Sun. On the other hand, the AL parallax effect depends on the 
corresponding barycentric quantities $d_\text{b}(t)$ and $\Omega_\text{b}(t)$ measured
relative to the solar system barycentre. More precisely, 
the ability to determine absolute parallaxes depends on the AL parallax 
factor being different in the two FoVs (see Fig.~2 in \citeads{2016A&A...595A...1G}). 
Relevant for the parallax zero point is therefore the 
differential AL parallax factor $d_\text{b}(t)\sin\xi_\text{b}(t)\cos\Omega_\text{b}(t)$
(see footnote~\ref{fn:Omega} in Appendix~\ref{sec:plxAc}), where $\xi_\text{b}(t)$ 
is the angle between \textit{Gaia}'s spin axis and the direction to the barycentre.
Thus, in our model for $\delta\Gamma(t)$, only the term proportional to 
$\cos\Omega(t)$, that is the one containing $\delta C_{1,0}$, has a strong 
correlation with the parallax zero point.
The differences between the heliocentric and barycentric 
quantities, of the order of 0.01~au and 0.01~rad, and the annual variations in
$d$ and $d_\text{b}$, by about $\pm 1.7$\%, all contribute 
towards a decorrelation of the parallax zero point from $\delta C_{1,0}$. 

For the actual cycle~3 data, the correlation coefficient between $\delta C_{1,0}$ 
and the parallax zero point is $\simeq\,$0.99992. Such a high degree of correlation 
(collinearity) in a least-squares estimation problem would normally be considered 
crippling, but it need not be so if the number of observations is very high, which it is 
in this case, and the modelling is sufficiently accurate, which we strive for. 
If all the $\sim$85~million unknowns in the primary astrometric solution could be 
obtained by direct solution of the normal equations, a valid solution for $\delta C_{1,0}$ 
would be obtained because the full normal equations take 
into account the correlations among all parameters. However, we are forced to 
use iterative solution methods, and it turns out that the introduction of $\delta C_{1,0}$ 
in AGIS effectively prevents the convergence of the block-iterative solution in its
original form. The non-convergence is however not caused by the strong correlation 
itself, but by the circumstance that the correlated parameters are in different blocks.
In AGIS the different blocks of source, attitude, calibration, and global parameters are 
treated as independent least-squares problems in a given iteration, thus ignoring 
correlations between, for example, the global block (containing $\delta C_{1,0}$)
and the source blocks (containing the parallaxes) when updates for the next iteration 
are computed.

We nevertheless found a way to obtain a converged solution including $\delta C_{1,0}$, by
using a special option in the global block of AGIS, called `consider parameters'. This device
was originally introduced for a different purpose,%
\footnote{The term `consider parameter' has various meanings in
the literature. Here we refer to a parameter that is included in the
estimation of updates to the current parameter values, but for which
no actual update is applied to the parameter. The consider parameter thus remains
at its original value (in this case zero), but the solution computes updates to the
other parameters, and uncertainties and correlations among all parameters,
exactly as if the consider parameter had been included in the fit. As the name suggests,
consider parameters are intended to help the researcher decide whether a particular
signal, modelled by the consider parameters, exists in the data, and how the covariance
of the solution would be affected if they were included in the fit.}
but here it is used to allow the AGIS iterations to converge in a reasonable time.
This particular use of consider parameters has been thoroughly tested in simulations,
and we are therefore confident in its fundamental correctness. Briefly, here is how it
works. In the global block, we introduce three more unknowns (consider parameters) 
that are strongly correlated with $\delta C_{1,0}$, namely one additive constant to all 
parallaxes, and two parameters for certain variations of the attitude (see Eq.~15 in 
\citeads{2017A&A...603A..45B}). All three consider parameters are fully degenerate with 
the parallaxes and the ordinary attitude parameters, but because their updates are ignored
in each iteration they remain at their initial zero values and do not affect the computation
of the right-hand side of the observation equations (the residuals). Their inclusion in
the left-hand side of the global block does however modify the updates to the regular
global parameters, including $\delta C_{1,0}$, and this is what allows the iterations to
converge. Because the three consider parameters remain at zero, it does not matter that
they are degenerate with other parameters, and the solution, after convergence, must be 
the same as a solution without them -- if such a solution could be obtained by some 
different algorithm. 
The role of the consider parameters in the block-iterative primary AGIS solution 
can formally be understood as a modification of the preconditioner of the adjustment 
scheme (e.g.\ \citealt{doi:10.1137/1.9780898718003}; \citeads{2012A&A...538A..77B}).

The condition number of the normal matrix for the fit of $\delta C_{1,0}$ 
is about $10^5$, so its inversion using normal 64-bit arithmetic is 
quite accurate. Although somewhat delicate, the fit works in practice and 
delivers a reasonably stable value of $\delta C_{1,0}$ after a number of 
AGIS iterations. The formal uncertainty of $\delta C_{1,0}$ from the fit 
is about 1~$\mu$as.
However, the fragile character of the fitting of $\delta C_{1,0}$ necessitates 
certain precautions: (i) $\delta C_{1,0}$ and the consider parameters 
should only be introduced at the very last stage of the AGIS iterations 
(see Table~\ref{tab:iterations}); (ii) only AL data should be used in the fit; 
(iii) in EDR3, the observations in window class WC0b have larger systematics
and were therefore omitted from the fit. 

In the primary AGIS solution for EDR3, $\delta C_{1,0}$ shifted the parallax zero 
point by about $+20~\mu$as compared with the same solution without $\delta C_{1,0}$, 
and by about $+10~\mu$as compared with DR2. The global parallax zero point of 
EDR3 is about $-17$~$\mu$as \citep{EDR3-DPACP-132}.
Although the inclusion of $\delta C_{1,0}$ in the global model for EDR3 did not bring
the global parallax zero point to zero, the partial success of the method is very 
encouraging and fosters the hope that the zero point issue can be resolved, at the level
of a few $\mu$as, in future releases that will benefit from much improved calibration models.


\section{Astrometric solutions}
\label{sec:agis}

\subsection{Main steps of the solutions}
\label{sec:steps}

The tasks labelled AGIS~3.1 and AGIS~3.2 in Fig.~\ref{fig:timeline} each consists of
several steps, the most important ones being:
\begin{enumerate} 
\item
Preprocess the input data (transits) from IDU: This includes filtering (removing transits
that are unmatched or of poor quality according to IPD flags, or outside the specified
time interval) and sorting the transits by position. Sorting uses the healpix index
\citepads{2005ApJ...622..759G} encoded in the \gacs{source\_id}.
\item
Select a set of primary sources to ensure a sufficient density of well-behaved sources with
a good coverage in magnitude and colours.
\item
Fit an initial attitude for the required time interval, using source parameters from a 
previous cycle or phase; also define data gaps where transits are missing or of poor quality.
\item
Generate the corrective attitude from rate data as described in Sect.~3.2 of 
\citetads{2018A&A...616A...2L}.
\item\label{step:prim}
Calculate a primary solution by simultaneously estimating source (S), attitude (A),
calibration (C), and global (G) parameters in an iterative least-squares solution involving
only the primary sources. See Sect.~\ref{sec:prim} for a brief explanation.
\item\label{step:Cprime}
Compute a separate set of calibration parameters (C$'$) for sources where
IPD used the default colour $\nu_\text{eff}^\text{\,def}=1.43~\mu$m$^{-1}$. This
calibration is based on a subset of the primary sources where image parameters 
were determined by IPD using both the actual colours and the default value.
\item
Calculate secondary solutions for all sources (Sect.~\ref{sec:sec}). 
The computation is equivalent to the S block in step~\ref{step:prim},
except that sources with default colour obtain six-parameter solutions
using calibration C$'$. In this step the acceptance criteria detailed in Sect.~\ref{sec:accept}
are checked and, if necessary, a fallback solution computed. 
\item
Postprocess the results: This includes calculating various statistics such as the
renormalised unit weight error (RUWE).
\item\label{step:CoMaRa}
Regenerate attitude and calibration data for use by downstream processes such as 
PhotPipe. This fills some time gaps and intervals (including the EPSL) that were excluded 
for the astrometry, but where the observations may still be useful for other processes.
The skymapper (SM) geometry is also calibrated at this point. Although the SM observations 
are not used in the astrometric solution, they are needed in downstream processes.
\item
Regenerate attitude and astrometric calibration data for the LSF and PSF calibrations in the  
CALIPD of the next processing cycle or phase. This uses the same calibration model as for 
the primary solution, but including only the purely geometric part of the model, that is the 
effects numbered 1, 2, 8, and 9 in Table~\ref{tab:cal}. This is known as the `NoCoMaRa' 
calibration: no dependency on colour, magnitude, or rate (as opposed to the normal,
`CoMaRa', calibration including all the effects). The rationale for this is that all 
dependencies on colour, magnitude, AC rate, saturation, subpixel phase, and CTI
effects should ultimately be accounted for by the LSF and PSF calibrations, so that 
AGIS can be a purely geometric solution. This goal will never be reached if the AGIS 
calibration used for the LSF and PSF calibrations already removes (part of) the dependencies. 
By using NoCoMaRa for the LSF and PSF calibrations, the latter processes see the full extent 
of the dependencies. (In principle the attitude generated in steps~\ref{step:prim} and
\ref{step:CoMaRa} is already purely geometric, but owing to the non-orthogonality 
of some CoMaRa and NoCoMaRa effects, the best-fitting geometric attitude is slightly 
different for the two calibrations.) The NoCoMaRa calibration and attitude are not used 
by any downstream processes, only by IDU.
\item
Export all results to the main database, making them available to other processes.
\end{enumerate}
The same steps were executed in AGIS~3.1 and 3.2, but with many differences
in the details. In particular, the selection of primary sources and the calibration models
were different in the two phases, and numerous improvements and bug fixes were 
implemented in between. The models described in Sect.~\ref{sec:models}, and all other 
details given hereafter, refer to AGIS~3.2.

\subsection{Primary solution for AGIS~3.2}
\label{sec:prim}

Although all steps listed in the previous section are needed for a successful 
astrometric solution, the primary solution (step~\ref{step:prim}) is by far the most 
important and difficult one. As described elsewhere
\citepads{2012A&A...538A..78L}, the primary solution iteratively updates the four
kinds of unknowns (source, attitude, calibration, and global parameters). The algorithm
can be described in terms of four separate blocks, designated S, A, C, and G.
In S the astrometric parameters of the primary sources are updated based on current 
values for the other unknowns; in A the attitude parameters are updated based on
current source, calibration, and global parameters; and so forth. The blocks are
normally executed in a cyclic manner, for example S-A-C-G-S-A-C-G-S-A-$\dots$, 
where S-A-C-G constitutes one iteration. (For specific purposes, some of the blocks 
may be left out, meaning that the corresponding unknowns are kept fixed.)
In the simple iteration (SI) algorithm, there is no memory of the updates in previous
iterations that can be used to optimise the next update; this algorithm 
reliably converges in all relevant cases
and is numerically very stable, but may require many iterations for complete convergence. 
The conjugate gradient (CG) algorithm \citepads{2012A&A...538A..77B} speeds up 
convergence considerably, but is less stable and (unlike SI) does not allow observation 
weights to be changed from one iteration to the next. Weight adjustment is necessary 
for a good treatment of outliers and for estimating the excess noise. Most often
AGIS employs a hybrid scheme consisting of three SI iterations during which the weights 
are adjusted, followed by three CG iterations with fixed weights; this sequence is 
then repeated as many times as required. A complete run typically ends with a 
sequence of simple iterations, confirming that the solution is sufficiently converged.

The primary solution for AGIS~3.2 processed about 6.5~billion ($6.5\times 10^9$) 
CCD observations for 14.3~million primary sources. The solution determined 
71.5~million source parameters together with 10.7~million attitude parameters, 
1.1~million calibration parameters, and 2.0~million global parameters; the redundancy 
factor (mean number of observations per unknown) is $\simeq\,$76.
(This does not count the corrective attitude, for which a largely different set of
observations was used to estimate some 17~million parameters.)

The sequence of iterations executed in the primary solution for AGIS~3.2 is 
detailed in Table~\ref{tab:iterations}. After a warm-up run to obtain good starting
values, a total of 165 iterations were made using the global model described in 
Sect.~\ref{sec:glob}. In the last 39 iterations the global parameter $\delta C_{1,0}$
was also adjusted, resulting in a significant reduction of the (negative)
parallax bias (Sect.~\ref{sec:C10}).



\begin{table}
\caption{Iteration sequences for the primary solutions of AGIS~3.2.  
\label{tab:iterations}}
\small
\begin{tabular}{lllll}
\hline\hline
\noalign{\smallskip}
Iterations & Blocks active & Algorithm & Remark\\
\noalign{\smallskip}
\hline
\noalign{\smallskip}
(1--27) & SAC & hybrid & warm up run\\
1--89 & SACG & hybrid & $\delta C_{1,0}=0$ \\
90--126 & SACG & simple & $\delta C_{1,0}=0$ \\
127--141 & SACG & hybrid & $\delta C_{1,0}$ free \\
142--165 & SAG & simple & $\delta C_{1,0}$ free, ad hoc corr. \\
166--181 & C & simple & step~6 (calibration C$'$) \\ 
\noalign{\smallskip}
\hline
\end{tabular}
\tablefoot{Columns~2 and 3 describe the AGIS configuration for the sequence of
iterations in the first column (see text).
$\delta C_{1,0}$ is the global parameter discussed in Sect.~\ref{sec:C10}. 
Iterations~1--165 correspond to step~5 in Sect.~\ref{sec:steps}.
Between iterations 141 and 142 the ad hoc correction to the calibration parameters 
for the bright sources (Sect.~\ref{sec:adHoc}) was applied, after which the calibration
was not updated. Iterations 166--181 correspond to step~6, computing the special 
calibration C$'$ needed for sources with default colour in IPD.
}
\end{table}



\subsection{Secondary solutions}
\label{sec:sec}

In this step the astrometric parameters were computed for all sources, including the primary sources.
Depending on the colour information used in the IPD for a particular source, it received a 
five- or six-parameter solution as described in Sect.~\ref{sec:nuEff}, but otherwise the treatment
was identical. The five-parameter solutions used calibration C obtained in step~5 of 
Sect.~\ref{sec:steps}, while the six-parameter solutions used calibration C$'$ obtained in step~6. 
For sources with an insufficient number of observations, or where the astrometric results failed to meet
the acceptance criteria for a five- or six-parameter solution (Sect.~\ref{sec:accept}), only 
the mean position at the reference epoch (J2016.0) is published. 

The secondary solutions processed nearly 78~billion FoV transits, generating converged
solutions for 2.495~billion sources (of which 585~million five-parameter, 883~million six-parameter, 
and 1027~million two-parameter solutions). Subsequently some of the five- and six-parameter solutions 
and most of the two-parameters solutions were removed because they failed to meet the acceptance 
criteria (Sect.~\ref{sec:accept}). The final number of sources and other statistics are given in 
Sect.~\ref{sec:results}.

  
\subsection{Acceptance criteria and fallback (two-parameter) solutions}
\label{sec:accept}  

The decision whether a converged, non-duplicated secondary solution is accepted as a 
five- or six-parameter solution, or at all retained for publication, depends on the four
quantities $N_\text{tr, astr}$, $N_\text{vpu}$, $\sigma_\text{pos, max}$, and $\sigma_\text{5d, max}$
calculated in the course of the source update process. Here,
$N_\text{tr, astr}$ is the number of FoV transits (detections) used in the AGIS solution; in
the \textit{Gaia} Archive it is given as \gacs{astrometric\_matched\_transits}.
$N_\text{vpu}$ (\gacs{visibility\_periods\_used}) is the number of distinct observation epochs 
(visibility periods) used in the solution, where a visibility period is a group of observations 
separated from other groups by a gap of at least four days.
$\sigma_\text{pos, max}$ (not in the \textit{Gaia} Archive) is the semi-major axis of the 
error ellipse in position at the reference epoch J2016.0 (Eq.~B1 in \citeads{2018A&A...616A...2L}).
Finally, $\sigma_\text{5d, max}$ (\gacs{astrometric\_sigma5d\_max}) is the five-dimensional 
equivalent to $\sigma_\text{pos, max}$, calculated as described in Sect.~4.3 of 
\citetads{2018A&A...616A...2L} but with $T=2.76383$~yr for the time coverage of the data 
used in the solution and ignoring the pseudocolour for six-parameter solutions.

For every source, a solution with five or six parameters (depending on the colour information used in
IPD) was first tried. This was accepted if it converged and satisfied the
criterion
\begin{equation}\label{eq:critP56}
G\le 21.0~~~  \& ~~~ N_\text{vpu}\ge 9 ~~~
\& ~~~ \sigma_\text{5d, max}<(1.2~\text{mas})\times\gamma(G)\, ,
\end{equation}
where $\gamma(G)=10^{0.2\max(6-G,~0,~G-18)}$.
This is similar to the DR2 criterion (Eq.~11 in \citeads{2018A&A...616A...2L}), except that the minimum
$N_\text{vpu}$ is higher and the upper limit on $\sigma_\text{5d, max}$ was increased for $G<6$ 
to accommodate the sharply rising uncertainty for the brightest sources (Fig.~\ref{fig:uncVsG}).
The present threshold on $N_\text{vpu}$ removes most cases where a lower threshold might produce 
spurious solutions, like the ones found in DR2 with very large (positive or negative) parallaxes.
We note that the $G$ used in Eq.~(\ref{eq:critP56}) is not the EDR3 value, which was unavailable 
at the time, but the value from DR2, or the real-time magnitude estimate from the on-board 
object detection if the source was not in DR2. In EDR3 there are 143\,546 sources with 
five- or six-parameter solutions and EDR3 magnitude $G > 21$, and conversely some sources with 
two-parameter solutions that would have passed Eq.~(\ref{eq:critP56}) if the EDR3 magnitude
had been used. (Elsewhere in this paper $G$ stands for the EDR3 value \gacs{phot\_g\_mean\_mag}.)

If the five- or six-parameter solution did not converge, 
or failed to satisfy Eq.~(\ref{eq:critP56}), prior information
on the parallax and proper motion was added, based on the Galactic model described in 
\citetads{2015A&A...583A..68M} and Sect.~4.3 of \citetads{2018A&A...616A...2L}. In such
cases only the position parameters ($\alpha$, $\delta$) at epoch J2016.0 and their covariances 
were retained out of the full
five- or six-parameter solution. As explained in \citetads{2015A&A...583A..68M}, the purpose of 
the Galactic prior is to provide more realistic uncertaintainties for the positions of sources with
a very small number of observations, by making some reasonable assumption about the
sizes of their parallaxes and proper motions. The resulting position is called a two-parameter 
solution ($\gacs{astrometric\_params\_solved}=3$), although 
in reality all five or six parameters are estimated. A two-parameter solution was accepted for 
publication if it satisfied the criterion
\begin{equation}\label{eq:crit}
N_\text{tr, astr}\ge 5 \quad \& \quad \sigma_\text{pos, max}<100~\text{mas} \, .
\end{equation}
It can be noted that any solution that satisfies Eq.~(\ref{eq:critP56}) also satisfies Eq.~(\ref{eq:crit}), 
which therefore holds for all sources in EDR3. In contrast to the corresponding criterion for 
\textit{Gaia} DR2 (Eq.~12 in \citeads{2018A&A...616A...2L}), Eq.~(\ref{eq:crit}) puts no upper limit 
on the \gacs{astrometric\_excess\_noise}, because it was found that such a limit rejects many 
partially resolved binaries that should be retained in the catalogue for completeness,
even though they do not have full astrometric data.

Finally, all sources, irrespective of the kind of solution, must be solitary in the sense that
there is no other source within a radius of 0.18~arcsec, as calculated from the position 
parameters ($\alpha$, $\delta$) at the reference epoch. If multiple sources are found at 
smaller separations, only one source is kept, namely, in order of precedence:
(i) a source previously identified as relevant for the extragalactic reference frame; 
(ii) the five- or six-parameter solution with the smallest $\sigma_\text{5d, max}$; or 
(iii) the two-parameter solution with the smallest $\sigma_\text{5d, max}$. 
In such cases the retained source has the flag \gacs{duplicated\_source} set in the
\textit{Gaia} Archive. 


\subsection{Ad hoc correction of WC0 calibration}
\label{sec:adHoc}   

Between iterations 141 and 142 an ad hoc correction was applied to the WC0
calibration parameters in order to mitigate a known problem with the bright 
($G\lesssim 13$) reference frame. The motivation and procedure for this 
correction, which should not be needed in future processing cycles, are briefly
as follows.

During the internal validation of the AGIS~2.2 solutions, carried out by the
astrometry team prior to the publication of \textit{Gaia} DR2, it was found that
the reference frame of the bright sources ($G\lesssim 12$--13) in DR2 was 
rotating, relative to the frame defined by the fainter quasars, at a rate of 
about 0.15~mas~yr$^{-1}$ (Sect.~5.1 of \citeads{2018A&A...616A...2L}).
The problem was confirmed by \citetads{2018ApJS..239...31B} in a 
comparison with proper motions calculated from the position differences 
between the DR2 and \textsc{Hipparcos} catalogues, and by
\citetads{2020A&A...633A...1L,2020A&A...637C...5L} in a comparison
with radio-interferometric (VLBI) observations of bright radio stars.
The likely cause of the effect is explained in Appendix~B 
of \citetads{2020A&A...633A...1L}.

A similar effect was seen during the production of AGIS~3.1. A major concern
then was that these systematics, if left uncorrected in AGIS~3.1, would propagate 
into the time-dependent LSF and PSF calibrations of WC0 sources in CALIPD~3.2,
only to appear again in AGIS~3.2. It was therefore decided to implement an 
ad hoc correction to the calibration parameters of WC0 in AGIS~3.1,
counteracting the effect. This procedure successfully fixed the bright reference 
frame in AGIS~3.1, but the problem nevertheless reappeared in AGIS~3.2,
albeit with different values. 
A similar ad hoc correction was therefore made after iteration 141 in the 
AGIS~3.2 iteration sequence (Table~\ref{tab:iterations}). Because the calibration
was not updated in the subsequent iterations, the correction remained effective
in the final results.

To explain the correction it is useful to consider how the AL astrometric measurements 
are affected by a change in the source positions corresponding to a small error in the 
celestial reference frame. The orientation error at a certain time is given by the 
(numerically small) rotation vector $\vec{\varepsilon}$, such that the change in the 
unit vector $\vec{u}$ towards a source is $\Delta\vec{u}=\vec{\varepsilon}\times\vec{u}$. 
Let $\vec{z}$ be the unit vector, at the same instant, along the nominal spin axis of 
\textit{Gaia}  (more precisely, $\vec{z}$ is the third axis of the scanning reference 
system SRS; e.g.\ Fig.~2 in \citeads{2012A&A...538A..78L}).
$\vec{z}$ and $\vec{u}$ must be nearly orthogonal for the source to be observed 
in one of the FoVs, and for simplicity we assume $\vec{z}'\vec{u}=0$. The 
tangent vector of the AL field angle $\eta$ at the source is then the unit vector 
$\vec{z}\times\vec{u}$, and the component of $\Delta\vec{u}$ in the AL direction is 
$\Delta\eta=(\vec{z}\times\vec{u})'\Delta\vec{u} = 
\vec{z}'\vec{\varepsilon}-(\vec{z}'\vec{u})(\vec{u}'\vec{\varepsilon})=\vec{z}'\vec{\varepsilon}$.
Both $\vec{\varepsilon}$ and $\vec{z}$ are functions of time, thus
$\Delta\eta(t) = \vec{z}(t)'\vec{\varepsilon}(t)$. Here $\vec{z}(t)$ is set by the 
scanning law, while the standard model of stellar motion (Sect.~\ref{sec:source}) 
requires that the frame orientation error is a linear function of time, 
$\vec{\varepsilon}(t)=\vec{\varepsilon}(t_\text{ep})+(t-t_\text{ep})\vec{\omega}$.
The function $\Delta\eta(t)$ therefore has six degrees of freedom corresponding to the 
components of the vectors $\vec{\varepsilon}(t_\text{ep})$ and $\vec{\omega}$. The
important conclusion from this brief discussion is that only very specific forms of
time-dependent AL displacements in the calibration model could be mistaken for  
a reference frame error. 

In the astrometric calibration
model (Sect.~\ref{sec:cal}), the AL large-scale calibration for WC1 ($G\simeq 13$ to 16) has 
a fixed origin, when averaged over both FoVs, but for WC0 and WC2 it is necessary
to permit time-dependent displacements relative to WC1. This means that each WC could 
in principle have its own reference frame, namely if the relative displacement between their
calibrations can be described in the form of the function $\Delta\eta(t)$ introduced above
for some vectors $\vec{\varepsilon}(t_\text{ep})$ and $\vec{\omega}$. 
In practice this should not be a problem, because many primary 
sources around magnitude 13 and 16 are not always observed in the same WC, and they
will only obtain consistent solutions if the reference frame is the same in all WC. This
mechanism apparently works as expected for the transition between WC1 and WC2 around 
$G=16$, but not for the transition between WC1 and WC0 around $G=13$. The probable
reason for this is the generally problematic calibrations in WC0, both in CALIPD and AGIS.

The ad hoc correction amounts to adding the time-dependent correction 
$\Delta\eta(t)=(t-t_\text{ep})\vec{z}(t)'\vec{\omega}$ to the WC0a and WC0b calibrations, where 
$\vec{\omega}=[-0.0166,~-0.0950,~+0.0283]'$~mas~yr$^{-1}$ was estimated from 
a comparison of the proper motions of \textsc{Hipparcos} stars, as obtained in iteration 
141 from the \textit{Gaia} observations, and as derived from the positional differences 
between \textit{Gaia} DR2 and the \textsc{Hipparcos} catalogue \citepads{2007ASSL..350.....V}.
For lack of better information, it was necessary to assume that the positional
systems of the different window classes agreed at the reference epoch, in other words that
$\vec{\varepsilon}(t_\text{ep})=\vec{0}$. In effect, the applied correction implies that the 
bright reference frame of EDR3, when extrapolated to the \textsc{Hipparcos} epoch 
J1991.25, agrees with the \textsc{Hipparcos} reference frame. The uncertainty in the 
alignment of the \textsc{Hipparcos} reference frame to the ICRS at epoch J1991.25 
was $\pm 0.6$~mas in each axis \citepads{1997A&A...323..620K}, which gives a systematic 
uncertainty of at least $(0.6~\text{mas})/(24.75~\text{yr})\simeq 0.024$~mas~yr$^{-1}$ 
per axis in the spin of the bright reference frame of \textit{Gaia}~EDR3.

It is important to note that the ad~hoc correction does not adjust the proper motions individually
for agreement with the \textsc{Hipparcos} positions (as was done for the \textit{Gaia} DR1 TGAS 
solution; \citeads{2016A&A...595A...2G}); only the reference frame is adjusted via the WC0
calibration parameters. Nevertheless, resorting to this procedure is very unsatisfactory
and hopefully exceptional: Improved calibrations in CALIPD and AGIS for the WC0 observations 
should eliminate the need for it in future releases. However, it highlights the need for 
independent means to verify the consistency of the \textit{Gaia} reference frame over the full 
range of magnitudes, for example by means of VLBI observations of radio stars 
\citepads{2020A&A...633A...1L}.

\begin{figure}
\center
  \resizebox{0.9\hsize}{!}{\includegraphics{fig05.pdf}}
    \caption{Magnitude distribution of sources in \textit{Gaia} EDR3. 
	The grey denotes all sources, the blue denotes sources with five-parameter solutions, 
	the green denotes sources with six-parameter solutions, and the red denotes sources 
	with two-parameter solutions.}
    \label{fig:histNsrc}
\end{figure}



\section{Results: Astrometric properties of EDR3}
\label{sec:results}

\subsection{Overview of the data}
\label{sec:overviewData}

The main table of \textit{Gaia} EDR3, \gacs{gaia\_source}, gives
astrometric data for more than 1.8~billion sources. The exact numbers are
585\,416\,709 sources with five-parameter solutions ($\gacs{astrom\_params\_solved}=31$),
882\,328\,109 with six-parameter solutions ($\gacs{astrom\_params\_solved}=95$), and
343\,964\,953 with two-parameter solutions ($\gacs{astrom\_params\_solved}=3$).
In total there are 1\,811\,709\,771 sources. 
Their distribution in $G$ magnitude (\gacs{photometric\_g\_mean\_mag})
is shown in Fig.~\ref{fig:histNsrc}.

In the following we give statistics related to the quantities 
listed below with their brief explanations.
\begin{itemize}
\item \gacs{ra\_error} = standard uncertainty in right ascension at epoch J2016.0, 
$\sigma_{\alpha*}=\sigma_\alpha\cos\delta$
\item  \gacs{dec\_error} = standard uncertainty in declination at epoch J2016.0, 
$\sigma_{\delta}$
\item  \gacs{parallax\_error} = standard uncertainty in parallax, 
$\sigma_{\varpi}$
\item \gacs{pmra\_error} = standard uncertainty of proper motion in right ascension, 
$\sigma_{\mu\alpha*}=\sigma_{\mu\alpha}\cos\delta$
\item \gacs{pmdec\_error} = standard uncertainty of proper motion in declination, 
$\sigma_{\mu\delta}$
\item \gacs{pseudocolour\_error} = standard uncertainty of the pseudocolour,
$\sigma_{\hat{\nu}\text{eff}}$ 
\item semi-major axis of error ellipse in position at epoch J2016.0, $\sigma_\text{pos,max}$
(Eq.~B.1 in \citeads{2018A&A...616A...2L})
\item semi-major axis of error ellipse in proper motion, $\sigma_\text{pm,max}$
(Eq.~B.2 in \citeads{2018A&A...616A...2L})
\item \gacs{ruwe} = renormalised unit weight error (RUWE). The unit weight error (UWE)
is the square root of the normalised chi-square of the astrometric fit to the AL observations, 
$\text{UWE}=[\,\chi^2/(n-n_\text{p})]^{1/2}$, where $n$ is the number of 
good CCD observations of the source (see below) and $n_p=5$ or 6 the number of
parameters fitted. $\text{UWE}\simeq 1.0$ is expected for a well-behaved source, but
that is often not the case owing to calibration errors. The RUWE is calculated by
empirical scaling of the UWE, depending on $G$ and $\nu_\text{eff}$ or 
$\hat{\nu}_\text{eff}$, such that $\text{RUWE}\simeq 1.0$ for well-behaved sources
(see also Sect.~\ref{sec:gof}). This statistic is not given for two-parameter solutions.
\item \gacs{astrometric\_excess\_noise} = excess source noise, $\epsilon_i$: This
is the extra noise per observation that must be postulated to explain the scatter
of residuals in the astrometric solution for the source (see also Sect.~\ref{sec:gof}). 
The excess source noise is considered to be statistically significant if 
$\gacs{astrometric\_excess\_noise\_sig}>2$
\item  \gacs{visibility\_periods\_used} = number of visibility periods of the source,
that is, groups of observations separated by at least four days
\item  \gacs{astrometric\_matched\_observations} = number of FoV
transits of the source used in the astrometric solution
\item \gacs{astrometric\_n\_good\_obs\_al} = number of good
CCD observations AL of the source used in the astrometric solution 
\item fraction of outliers (bad CCD observations AL of the source) = 
$\gacs{astrometric\_n\_bad\_obs\_al}/\gacs{astrometric\_n\_obs\_al}$
\item \gacs{ipd\_gof\_harmonic\_amplitude} = amplitude of the natural logarithm of the
goodness-of-fit obtained in the IPD versus position angle of scan
\item \gacs{ipd\_frac\_multi\_peak} = fraction of CCD observations where IPD detected 
more than one peak
\item \gacs{ipd\_frac\_odd\_win} = fraction of FoV transits with truncated windows 
or multiple gates
\end{itemize}
The last three statistics were generated at the image parameter determination (IPD) stage
prior to the astrometric solution, and may include transits that were not used for the
astrometry. They are listed here because they provide information on (potentially) 
problematic sources, complementary to what is obtained from the astrometric 
fit (see Sect.~\ref{sec:gof}). We refer to the \textit{Gaia} Archive 
on-line documentation for further explanation of these statistics.

Although both the five- and six-parameter solution provides estimates of the five 
astrometric parameters (position, parallax, and proper motion), the six-parameter 
solution (with pseudocolour as the sixth parameter) is intrinsically less accurate 
because the default colour had to be used for the IPD. Moreover, the six-parameter 
solution is normally only used for sources that are problematic in some respect, 
for example in very crowded areas, which tends to reduce its accuracy even more. 
It is therefore usually relevant to give separate statistics for the two kinds of 
solution. In Tables~\ref{tab:statP5}--\ref{tab:statP2} we report the mean or median 
values of most of the statistics listed above, as functions of magnitude and separated 
by the kind of solution. The median is used for quantities that have a long-tailed 
distribution, for which the mean value might be less representative.

\begin{figure}
\centering
  \includegraphics[width=0.90\hsize]{fig06.pdf} 
  \caption{Neighbourhood of 15th magnitude sources in \textit{Gaia} EDR3.
  The diagrams show all sources within 5~arcsec from any one of the 2.8~million sources
  in EDR3 with $14.95<G<15.05$~mag.
    \textit{Top:} Five-parameter solutions.
    \textit{Middle:} Six-parameter solutions.
    \textit{Bottom:} Two-parameter solutions.}
  \label{fig:resG15}
\end{figure}


\subsection{Angular resolution}
\label{sec:resolution}

Resolution here refers to the minimum angular separation between distinct sources
in \textit{Gaia} EDR3, that is between objects with different \gacs{source\_id}. As
explained in Sect.~\ref{sec:accept}, the separation will by construction never be
smaller than 0.18~arcsec, but relatively few sources are found with separations less
than about 
0.6~arcsec owing to other limitations. In a given situation, the effective resolution 
(however it is defined) depends on many different factors such as the magnitudes 
of both components in a pair, their relative orientation on the sky, and the 
kinds of solutions involved (with five, six, or two parameters). The complex situation
is illustrated in Fig.~\ref{fig:resG15}, which shows the neighbourhood of all
sources in EDR3 of magnitude $G\simeq 15$, subdivided by the kind of solution.
Clearly most neighbours at separation 0.18--0.6~arcsec only have two-parameter
solutions, while neighbours with five-parameter solutions are usually either the
brighter of the two sources or more distant than 2~arcsec. The six-parameter 
solutions partly fill the gap for separations between 0.6 and 2~arcsec. Additional
statistics on the small-scale completeness of EDR3 are given in \citet{EDR3-DPACP-126}.


\subsection{Goodness-of-fit statistics}
\label{sec:gof}

Several of the statistics listed in Sect.~\ref{sec:overviewData} quantify the goodness-of-fit
of the single-star model to the observations, either at the image parameter determination
(IPD), where a model LSF or PSF is fitted to the CCD samples, or in the subsequent astrometric 
solution, where the standard model of stellar motion (Sect.~\ref{sec:source}) is fitted to the 
resulting image locations. A few remarks should be made concerning the interpretation of
these statistics and their interrelations.

For the user, the most relevant goodness-of-fit statistics from the IPD are 
\gacs{ipd\_gof\_harmonic\_amplitude}, \gacs{ipd\_frac\_multi\_peak}, and 
\gacs{ipd\_frac\_odd\_win}; and from the astrometric fit, the RUWE and the 
excess source noise with its significance. (The fraction of outliers is probably 
less useful: The outlier detection is designed to remove occasional large deviations, 
caused by temporary perturbations that are usually unrelated to the source.)
All of them describe (real or spurious) deviations from the simplest possible point-source 
model, but they are sensitive to different kinds of modelling errors, and all of them are more
or less sensitive to calibration errors. The sensitivity is usually a strong function 
of the magnitude of the source, and may also depend on geometric factors such as the 
distribution of scans across the source. All of this complicates the interpretation of the 
statistics. For example, there is no simple way to convert them into $p$-values, using the 
single-star model as a null hypothesis; instead, the relevant distributions must be 
determined empirically, if it is at all possible.
 
The IPD statistics may be quite powerful for detecting certain kinds of binaries.
\gacs{ipd\_gof\_harmonic\_amplitude} is normally small but could become large
for sources that have elongated images, such as partially resolved binaries, provided
that their position angles are relatively fixed. In such cases \gacs{ipd\_gof\_harmonic\_phase} 
indicates the position angle of the major axis modulo $180^\circ$. 
\gacs{ipd\_frac\_multi\_peak} is sensitive to resolved binaries that in some scan
directions produce more than one peak in the window. Finally, \gacs{ipd\_frac\_odd\_win} 
is sensitive to the presence of another (usually brighter) source causing the window to be 
`odd', that is truncated or with multiple gating modes. The source causing the gating
could however be quite far away on the CCD, or even in the other FoV. Transits with 
odd windows were not used in the astrometric solution for EDR3.
From Tables~\ref{tab:statP5}--\ref{tab:statP2} it is seen that sources with 
five-parameter solutions are usually very clean, as indicated by the IPD statistics,
while the six-parameter solutions have higher fractions of observations with 
multiple peaks or odd windows (which partially explains why they did not
have sufficiently good BP and RP photometry for the calculation of $\nu_\text{eff}$), 
and the two-parameter solutions are even worse. Towards the faint magnitudes the 
\gacs{ipd\_frac\_multi\_peak} is always decreasing, because the diminishing 
signal-to-noise ratio makes the detection of secondary peaks increasingly difficult.

The astrometric goodness-of-fit measures RUWE and excess source noise 
(\gacs{astrometric\_excess\_noise}) quantify the same thing, namely how much 
the motion of the image centre (as determined by the IPD) deviates from 
the standard model of stellar motion fitted in the astrometric solution. However, while
the \gacs{astrometric\_excess\_noise} gives the discrepancy in angular measure (mas) 
per AL observation (ideally $=0$ for a good fit), the RUWE gives the discrepancy as a 
dimensionless factor (ideally $=1.0$). The RUWE was obtained from the unit weight 
error $\text{UWE}=[\,\chi^2/(n-n_\text{p})]^{1/2}$ by applying an empirical scaling 
factor to compensate for calibration errors, which tend to increase the UWE for bright, 
blue, and very red sources. A corresponding correction was not applied to the 
\gacs{astrometric\_excess\_noise}, which must therefore be interpreted with 
some caution for sources with $G\lesssim 13$, $\nu_\text{eff}\gtrsim 1.65$ 
($G_\text{BP}{-}G_\text{RP}\lesssim 0.4$), or $\nu_\text{eff}\lesssim 1.24$ 
($G_\text{BP}{-}G_\text{RP}\gtrsim 3.0$). The significance (signal-to-noise ratio) 
of the excess source noise is given by \gacs{astrometric\_excess\_noise\_sig}: 
\gacs{astrometric\_excess\_noise} should be regarded as insignificant (that is, 
effectively zero) if $\gacs{astrometric\_excess\_noise\_sig}\lesssim 2$. Alternatively, 
$\gacs{astrometric\_excess\_noise}/\gacs{astrometric\_excess\_noise\_sig}$ may
be taken as an estimate of the uncertainty of the excess source noise.
The RUWE and excess source noise are sensitive to the photocentric motions of 
unresolved objects, such as astrometric binaries, which are not revealed by the
IPD statistics, and therefore complement the latter in binary detection.

The standard uncertainties given in EDR3 have been adjusted to take into account 
the excess noise, whether it represents an astrometric mismatch or a calibration 
issue; they should not be further inflated based on the goodness-of-fit 
statistics.



\begin{table*}[t]
\caption{Summary statistics for the 585~million sources in \textit{Gaia} EDR3 
with five-parameter solutions.
\label{tab:statP5}}
\small   \setlength{\tabcolsep}{4.5pt}
\begin{tabular}{lrrrrrrrrrrl}
\hline\hline
\noalign{\smallskip}
& \multicolumn{10}{c}{Value at $G=$} \\
Quantity & $9$--$12$& $13$ & $14$ & $15$ & $16$ & $17$ & $18$ & $19$ & $20$ & $21$ & Unit \\
\noalign{\smallskip}\hline\noalign{\smallskip}
Fraction of sources with 5-param.\ solution & 92.9& 97.0& 97.1& 96.7& 93.8& 87.5& 76.7& 49.7& 14.5&  1.4& \% \\
Median standard uncertainty in $\alpha$ ($\sigma_{\alpha*}$) at J2016.0 & 0.014& 0.012& 0.015& 0.022& 0.034& 0.054& 0.094& 0.175& 0.374& 1.006& mas \\
Median standard uncertainty in $\delta$ ($\sigma_{\delta}$) at J2016.0 & 0.013& 0.011& 0.013& 0.019& 0.030& 0.049& 0.086& 0.161& 0.335& 0.977& mas \\
Median standard uncertainty in $\varpi$ ($\sigma_{\varpi}$) & 0.018& 0.015& 0.019& 0.027& 0.042& 0.069& 0.120& 0.221& 0.459& 1.320& mas \\
Median standard uncertainty in $\mu_{\alpha*}$ ($\sigma_{\mu\alpha*}$) & 0.018& 0.016& 0.019& 0.028& 0.044& 0.071& 0.123& 0.229& 0.487& 1.445& mas~yr$^{-1}$ \\
Median standard uncertainty in $\mu_{\delta}$ ($\sigma_{\mu\delta}$) & 0.016& 0.014& 0.017& 0.024& 0.038& 0.063& 0.111& 0.208& 0.428& 1.401& mas~yr$^{-1}$ \\
Median renormalised unit weight error (RUWE) & 1.039& 1.023& 1.016& 1.013& 1.012& 1.011& 1.010& 1.010& 1.010& 1.016&   \\
Fraction with significant excess source noise & 97.3& 59.4& 19.2& 16.8& 13.8& 10.6&  7.5&  5.0&  5.0& 11.7& \% \\
Median excess source noise (when significant) & 0.113& 0.089& 0.188& 0.224& 0.285& 0.409& 0.606& 0.976& 1.801& 3.668& mas \\
Mean number of visibility periods used & 20.9& 21.0& 20.9& 20.8& 20.7& 20.7& 20.7& 20.5& 19.2& 13.1&   \\
Mean number of FoV transits used & 43.7& 43.7& 43.6& 43.4& 43.2& 43.2& 43.4& 42.6& 38.7& 20.5&   \\
Mean number of good CCD observations AL & 378& 378& 377& 376& 374& 374& 376& 368& 334& 175&   \\
Mean fraction of bad CCD observations AL &  1.2&  0.4&  0.5&  0.5&  0.5&  0.4&  0.4&  0.4&  0.5&  0.6& \% \\
Median \gacs{ipd\_gof\_harmonic\_amplitude} & 0.018& 0.021& 0.020& 0.020& 0.020& 0.020& 0.020& 0.022& 0.029& 0.073&  \\
Mean \gacs{ipd\_frac\_multi\_peak} &  1.9&  1.5&  1.5&  1.7&  0.9&  0.6&  0.3&  0.2&  0.1&  0.9& \% \\
Mean \gacs{ipd\_frac\_odd\_win} &  0.0&  0.0&  0.0&  0.0&  0.0&  0.0&  0.1&  0.1&  0.2&  0.8& \% \\
\noalign{\smallskip}\hline
\end{tabular}
\end{table*}

\begin{table*}[t]
\caption{Summary statistics for the 882~million sources in \textit{Gaia} EDR3 
with six-parameter solutions.
\label{tab:statP6}}
\small   \setlength{\tabcolsep}{4.5pt}
\begin{tabular}{lrrrrrrrrrrl}
\hline\hline
\noalign{\smallskip}
& \multicolumn{10}{c}{Value at $G=$} \\
Quantity & $9$--$12$& $13$ & $14$ & $15$ & $16$ & $17$ & $18$ & $19$ & $20$ & $21$ & Unit \\
\noalign{\smallskip}\hline\noalign{\smallskip}
Fraction of sources with 6-param.\ solution &  6.2&  2.4&  2.2&  2.7&  5.6& 11.9& 22.4& 48.6& 77.7& 27.5& \% \\
Median standard uncertainty in $\alpha$ ($\sigma_{\alpha*}$) at J2016.0 & 0.017& 0.026& 0.025& 0.033& 0.052& 0.081& 0.131& 0.225& 0.430& 1.031& mas \\
Median standard uncertainty in $\delta$ ($\sigma_{\delta}$) at J2016.0 & 0.016& 0.024& 0.023& 0.030& 0.045& 0.069& 0.112& 0.194& 0.382& 1.025& mas \\
Median standard uncertainty in $\varpi$ ($\sigma_{\varpi}$) & 0.023& 0.033& 0.031& 0.040& 0.061& 0.093& 0.151& 0.266& 0.526& 1.390& mas \\
Median standard uncertainty in $\mu_{\alpha*}$ ($\sigma_{\mu\alpha*}$) & 0.023& 0.034& 0.033& 0.042& 0.065& 0.102& 0.167& 0.291& 0.565& 1.469& mas~yr$^{-1}$ \\
Median standard uncertainty in $\mu_{\delta}$ ($\sigma_{\mu\delta}$) & 0.021& 0.031& 0.029& 0.036& 0.051& 0.078& 0.132& 0.240& 0.490& 1.470& mas~yr$^{-1}$ \\
Median standard uncertainty in $\hat{\nu}_\text{eff}$ ($\sigma_{\hat{\nu}\text{eff}}$) & 0.004& 0.006& 0.006& 0.008& 0.012& 0.019& 0.031& 0.055& 0.109& 0.281& $\mu$m$^{-1}$ \\
Median renormalised unit weight error (RUWE) & 1.127& 1.362& 1.167& 1.104& 1.102& 1.085& 1.053& 1.033& 1.029& 1.048&   \\
Fraction with significant excess source noise & 99.6& 95.1& 75.4& 68.3& 64.0& 58.3& 45.5& 27.8& 17.7& 19.6& \% \\
Median excess source noise (when significant) & 0.169& 0.275& 0.387& 0.405& 0.471& 0.569& 0.809& 1.343& 2.153& 3.838& mas \\
Mean number of visibility periods used & 21.4& 21.1& 20.7& 20.1& 18.9& 18.2& 18.4& 18.5& 18.1& 13.3&   \\
Mean number of FoV transits used & 46.1& 44.6& 44.0& 42.6& 39.3& 38.2& 37.9& 37.8& 36.6& 21.1&   \\
Mean number of good CCD observations AL & 400& 384& 379& 366& 339& 331& 328& 327& 315& 181&   \\
Mean fraction of bad CCD observations AL &  1.3&  1.5&  1.9&  2.1&  1.5&  1.0&  0.7&  0.6&  0.6&  0.6& \% \\
Median \gacs{ipd\_gof\_harmonic\_amplitude} & 0.020& 0.033& 0.032& 0.030& 0.035& 0.034& 0.032& 0.030& 0.031& 0.071&  \\
Mean \gacs{ipd\_frac\_multi\_peak} &  9.4& 25.4& 21.6& 18.3& 10.7&  7.0&  5.5&  3.2&  1.0&  0.5& \% \\
Mean \gacs{ipd\_frac\_odd\_win} &  0.0&  4.8&  7.9&  8.1&  6.3&  5.2&  4.4&  4.2&  3.8&  2.1& \% \\
\noalign{\smallskip}\hline
\end{tabular}
\end{table*}


\begin{table*}[t]
\caption{Summary statistics for the 344~million sources in \textit{Gaia} EDR3 with 
two-parameter solutions.
\label{tab:statP2}}
\small   \setlength{\tabcolsep}{4.5pt}
\begin{tabular}{lrrrrrrrrrrl}
\hline\hline
\noalign{\smallskip}
& \multicolumn{10}{c}{Value at $G=$} \\
Quantity & $9$--$12$& $13$ & $14$ & $15$ & $16$ & $17$ & $18$ & $19$ & $20$ & $21$ & Unit \\
\noalign{\smallskip}\hline\noalign{\smallskip}
Fraction of sources with 2-param.\ solution &  0.8&  0.7&  0.7&  0.6&  0.6&  0.6&  0.9&  1.7&  7.8& 71.1& \% \\
Median standard uncertainty in $\alpha$ ($\sigma_{\alpha*}$) at J2016.0 & 1.771& 1.260& 1.127& 1.116& 1.060& 1.008& 0.979& 1.436& 2.579& 3.250& mas \\
Median standard uncertainty in $\delta$ ($\sigma_{\delta}$) at J2016.0 & 1.703& 1.211& 1.078& 1.067& 0.994& 0.925& 0.888& 1.252& 2.106& 2.847& mas \\
Mean number of visibility periods used & 17.9& 18.4& 18.2& 18.0& 17.7& 17.0& 15.7& 10.9&  8.0&  6.9&   \\
Mean number of FoV transits used & 35.2& 36.3& 35.9& 35.4& 34.7& 33.6& 30.7& 19.7& 13.3& 10.0&   \\
Mean number of good CCD observations AL & 299& 306& 304& 300& 294& 285& 261& 167& 113& 85&   \\
Mean fraction of bad CCD observations AL &  4.3&  3.5&  3.5&  3.5&  3.5&  3.4&  3.1&  2.5&  1.6&  0.9& \% \\
Median \gacs{ipd\_gof\_harmonic\_amplitude} & 0.217& 0.251& 0.258& 0.261& 0.279& 0.273& 0.256& 0.259& 0.175& 0.096&  \\
Mean \gacs{ipd\_frac\_multi\_peak} & 59.2& 53.0& 44.8& 43.6& 36.0& 27.8& 20.8& 14.2&  6.8&  0.3& \% \\
Mean \gacs{ipd\_frac\_odd\_win} &  0.0& 12.5& 19.2& 19.8& 27.2& 31.1& 33.3& 35.6& 25.2&  3.1& \% \\
\noalign{\smallskip}\hline
\end{tabular}
\end{table*}



\begin{figure*}
\centering
  \includegraphics[width=0.95\hsize]{fig07.pdf}
  \caption{Uncertainty in parallax versus magnitude. \textit{Left:} Five-parameter solutions.
  \textit{Right:} Six-parameter solutions. The plots include all sources
  with $G<11.5$ and a geometrically decreasing random fraction of the fainter sources,
  so as to give a roughly constant number of sources per magnitude interval. The colour
  scale from yellow to black indicates an increasing density of data points in the diagram. 
  The curves show the 10th, 50th, and 90th percentiles of the distribution at a
  given magnitude.}
  \label{fig:uncVsG}
\end{figure*}

\begin{figure}
\centering
  \includegraphics[width=0.90\hsize]{fig08.pdf}
  \caption{Formal uncertainties at $G\simeq 15$ for sources with a five-parameter 
  solution in EDR3.
    \textit{Top:} Semi-major axis of the error ellipse in position at epoch J2016.0.
    \textit{Middle:} Standard deviation in parallax.
    \textit{Bottom:} Semi-major axis of the error ellipse in proper motion.
    These and all other full-sky maps in the paper except Fig.~\ref{fig:mapsCorrEcl}
    use a Hammer--Aitoff projection in equatorial (ICRS) coordinates with 
    $\alpha=\delta=0$ at the centre, north up, and $\alpha$ increasing from right to left.}
  \label{fig:mapsUnc}
\end{figure}

\begin{figure}
\centering
  \includegraphics[width=0.90\hsize]{fig09.pdf}
  \caption{Selected observation statistics at $G\simeq 15$ for sources with
    a five-parameter solution in EDR3. These statistics are main factors
    governing the formal uncertainties of the astrometric data.
    \textit{Top:} Number of visibility periods used.
    \textit{Middle:} Number of good CCD observations AL. (The number of FoV 
    transits used in the solution looks very similar, only with a factor nine smaller numbers.)
    \textit{Bottom:} Mean of $\log_{10}[\max(0.001,\,\gacs{astrometric\_excess\_noise})]$.}
  \label{fig:mapsVpu}
\end{figure}

\subsection{Formal uncertainties}
\label{sec:form}

Tables~\ref{tab:statP5}--\ref{tab:statP2} give the median uncertainties of the
astrometric parameters at selected magnitudes for the different kinds of solutions. 
No statistics are given for $G<9$ owing to the relatively few sources and rapidly 
declining quality at this end (see Fig.~\ref{fig:uncVsG}).
Compared with DR2 (Tables~B.1 and B.2 in \citeads{2018A&A...616A...2L}), the gain 
in median uncertainty at $G=15$ is a factor 0.71 for the positions and parallaxes,
and 0.44 for the proper motions. This is slightly better than the factors 0.80 and 0.51
expected purely from the increased length of the data included in the solutions
(Sect.~\ref{sec:data}). The extra gain comes mainly from the improved robustness 
and homogeneity of results made possible by the higher redundancy of observations 
in the EDR3 solutions. 

For $G=9$--12 the gain in median uncertainty from DR2 to EDR3 is even more 
impressive thanks to the improved calibrations, which are relatively more important 
for the bright sources (cf.\  Sect.~\ref{sec:resid}): The factor is 0.43 for the positions 
and parallaxes, and 0.27 for the proper motions.

The comparison between the two releases is however complicated by the circumstance 
that in EDR3 there are three kinds of solutions (five, six, and two parameters), while 
in DR2 there are only five- and two-parameter solutions. At $G=15$ the five-parameter 
solutions comprised 99.0\% of the sources in DR2 and 96.7\% of the sources in EDR3, 
and the comparison above used the statistics for these subsets. At $G=15$ the median 
uncertainties in EDR3 are a factor 1.5 higher for the 2.7\% of the sources with 
six-parameter solutions than for the 96.7\% with five-parameter solutions. This
large ratio in the uncertainties reflects the generally more problematic nature 
of the sources receiving six-parameter solutions, also seen in the various
goodness-of-fit statistics discussed in Sect.~\ref{sec:gof}.

The fraction of sources that receive five-parameter solutions is higher than 90\%
down to $G\simeq 17$, but decreases rapidly for fainter sources. The fraction with
six-parameter solutions correspondingly increases down to $G\simeq 20$, after 
which there is instead a steep increase in the fraction of two-parameter solutions. 

At any magnitude there is a considerable spread in the uncertainties caused by
variations in the number of observations and the properties of the scanning law. For 
$G\le 12$ there are additional variations depending on the window classes and gates 
used for a particular source, and the onset of saturation for
the brightest sources. The spread is illustrated in Fig.~\ref{fig:uncVsG} 
for the parallaxes in the five- and six-parameter solutions. The uncertainties in position
and proper motion follow similar distributions. 

Figure~\ref{fig:mapsUnc} shows the median uncertainties in position, parallax, and 
proper motions at $G\simeq 15$ as functions of position. For the position and proper
motion data the semi-major axes of the error ellipses $\sigma_\text{pos,max}$,
$\sigma_\text{pm,max}$ are plotted.
The patterns are very similar at other magnitudes, only scaled according to the
general dependence on $G$ in Table~\ref{tab:statP5} or Fig.~\ref{fig:uncVsG}.
These patterns are mainly set by variations in the number, direction, and temporal
distribution of the scans across a given position, as governed by the scanning law.
In very crowded areas, such as along the Galactic plane and in the general direction
of the Galactic Centre, the increased level of excess source noise from background 
sources gives a local rise in the median uncertainties, which becomes more 
important at fainter magnitudes. Some relevant statistics are shown in 
Fig.~\ref{fig:mapsVpu}. A comparison with the corresponding maps for DR2
(Figs.~B.3 and B.4 in \citeads{2018A&A...616A...2L}) clearly shows an 
improved homogeneity in the uncertainties in the ecliptic belt, and the
smaller importance of the excess noise in EDR3.

\begin{figure*}
\centering
  \includegraphics[width=0.79\hsize]{fig10.pdf}\hspace{0.19\hsize}
  \caption{Median correlation coefficients among the astrometric parameters for
  five-parameter solutions in \textit{Gaia} EDR3. The plots were made from a random
  selection of five-parameter sources with $G$ magnitude between 13 and 16. The correlations
  at other magnitudes are very similar.}
  \label{fig:mapsCorr5p}
\end{figure*}

\begin{figure*}
\centering
  \includegraphics[width=0.98\hsize]{fig11.pdf}
  \caption{Median correlation coefficients among the astrometric parameters for
  six-parameter solutions in \textit{Gaia} EDR3. The plots were made using all 
  six-parameter sources with $G$ magnitude between 13 and 16. The correlations at other
  magnitudes are very similar.}
  \label{fig:mapsCorr6p}
\end{figure*}

\begin{figure*}
\sidecaption
  \includegraphics[width=120mm]{fig12.pdf}
  \caption{Correlations in ecliptic coordinates (five-parameter solutions, $G=13$--16). 
  \textit{Left:} Median correlation coefficient between errors in ecliptic longitude ($\lambda$)
  and latitude ($\beta$). 
    \textit{Right:} Same for the proper motion components ($\mu_{\lambda*}$, $\mu_\beta$).
  Contrary to other sky maps in the paper, these two use a projection in
  ecliptic coordinates with $\lambda=\beta=0$ in the centre.
  }
  \label{fig:mapsCorrEcl}
\end{figure*}

\subsection{Correlation coefficients}
\label{sec:corr}

\textit{Gaia} EDR3 gives the complete set of correlation coefficients $\rho$ between 
the astrometric parameters provided for a given source. For a source with
$n_\text{p}=5$, 6, or 2 parameters, we thus have  
$n_\text{p}(n_\text{p}-1)/2=10$, 15, or 1 non-redundant coefficients.
In the \textit{Gaia} Archive they are called \gacs{ra\_dec\_corr}, etc.; here we use the notation
$\rho(\alpha,\delta)$, etc. The correlations allow the elements of the 
$n_\text{p}\times n_\text{p}$ covariance matrix $\vec{K}$ to be reconstructed as
\begin{equation}\label{eq:cov}
\begin{split}
K_{00} = \sigma_{\alpha*}^2\,, \quad
&K_{11} = \sigma_\delta^2\,, \quad\dots,\\
&K_{01} = K_{10} = \rho(\alpha,\delta)\sigma_{\alpha*}\sigma_\delta\,, \quad\dots
\end{split} 
\end{equation}
where indices $0,\,1,\,\dots$ represent the parameters in the usual order, 
$\alpha$, $\delta$, $\varpi$, $\mu_{\alpha*}$, $\mu_\delta$, $\hat{\nu}_\text{eff}$.

The correlation coefficients for a given source are mainly determined by the distribution 
of scan directions and transit times among the observations of the source, which are
governed by the scanning law. The correlation coefficients are therefore practically
independent of magnitude, and we give here only statistics for sources with $G=13$
to 16~mag.
Figures~\ref{fig:mapsCorr5p} and \ref{fig:mapsCorr6p}
show the median correlation coefficients for five- and six-parameter solutions.
We note that the scanning law is (approximately) symmetric with respect to the ecliptic, 
which is reflected in many features depending on ecliptic latitude ($\beta$) rather than 
declination ($\delta$). Furthermore, the patterns are often distinctly different
for $|\,\beta\,|\lesssim 45^\circ$ (the ecliptic belt) and $|\,\beta\,|\gtrsim 45^\circ$ 
(the ecliptic caps).

Certain features of predominantly positive or negative correlations are caused by
the choice of ICRS (equatorial) coordinates for the position and proper motion 
parameters, and are much less pronounced if ecliptic coordinates are used. This 
is the case, for example, with the mainly positive correlations $\rho(\alpha,\delta)$
and $\rho(\mu_{\alpha*},\mu_\delta)$ in the ecliptic belt for $\alpha=270^\circ$ to 
$90^\circ$, and their mainly negative values for $\alpha=90^\circ$ to $270^\circ$.
Geometrically, this can be understood in terms of the orientation of the error ellipses
in position and proper motion: In the ecliptic belt their major axes are approximately 
aligned with the ecliptic, and consequently tilted by up to $\pm 23.5^\circ$ with 
respect to the equator, corresponding to non-zero correlations in the equatorial
components. As shown in Fig.~\ref{fig:mapsCorrEcl}, $\rho(\lambda,\beta)$ and 
$\rho(\mu_{\lambda*},\mu_\beta)$ (here shown on an ecliptic projection) are generally
smaller than $\rho(\alpha,\delta)$ and $\rho(\mu_{\alpha*},\mu_\delta)$. 

The correlations $\rho(\alpha,\mu_{\alpha*})$ and $\rho(\delta,\mu_\delta)$ are 
related to the mean epoch of observations contributing to the different parameters. 
A mean epoch later than the reference epoch J2016.0 gives 
a positive correlation between the position and proper motion, and vice versa. 
This is especially pronounced for $\rho(\delta,\mu_\delta)$, where regions of positive 
and negative correlations alternate along the ecliptic.   
The correlations among $\alpha$, $\delta$, $\varpi$, $\mu_{\alpha*}$, and $\mu_\delta$
are almost the same for the five- and six-parameter solutions. In the six-parameter
solutions, the correlations between $\hat{\nu}_\text{eff}$ and the other five parameters are 
generally small (RMS values around 0.1) compared with the correlations among the five
parameters (RMS values of 0.2--0.3). 

The generally small correlations between pseudocolour and the other five parameters
is a consequence of the variation in chromaticity along the path of the stellar image 
through the AF, and between successive observations in either FoV
in a few revolutions; by contrast, the scan directions and observation epochs relevant 
for the other correlations do not change much over several revolutions. The sizes of the 
correlations with pseudocolour are important for the potential to improve the six-parameter 
solutions by incorporating external colour information (Appendix~\ref{sec:col6p}).

  
\begin{figure*}
\centering
  \includegraphics[width=0.90\hsize]{fig13.pdf}
  \caption{Smoothed maps of quasar parallaxes and proper motions.
  \textit{Left column:} \textit{Gaia} EDR3, using data for about 1.2~million quasars.
  \textit{Right column:} \textit{Gaia} DR2, using data for the 94\% of the quasars in the
  left column that have full astrometric solutions also in DR2.
  From top to bottom the maps show parallax, proper motion in right ascension, and
  proper motion in declination. The maps were smoothed using a Gaussian kernel with 
  standard deviation $5^\circ$.
    No data are shown for $|\,b\,|<10^\circ$, where $b$ is Galactic latitude.}
  \label{fig:mapsQso}
\end{figure*}

\begin{figure*}
\centering
  \includegraphics[width=0.93\hsize]{fig14.pdf}
  \caption{Smoothed maps of parallaxes in the LMC area, visualising small-scale 
  systematics (the `checkered pattern') in \textit{Gaia} EDR3 and DR2. 
  \textit{Left:} Smoothed parallaxes in EDR3 for sources in the magnitude 
  range $G=16$--18 (median $G=17.4$), kinematically selected as probable members 
  of the system (see Appendix~B in \citeads{EDR3-DPACP-132} for details).
  \textit{Right:} Smoothed parallaxes in DR2 for the same sample of sources.
  Both maps were smoothed using a Gaussian kernel with standard deviation $0.1^\circ$.
  While the sample includes about 730\,000 sources within $5^\circ$ radius of the 
  adopted centre, only smoothed points within a radius of $4.5^\circ$ are shown to 
  avoid unwanted edge effects. Comparison between the two diagrams is facilitated 
  by the use of the same colour scale, only shifted by $10~\mu$as to compensate for the 
  mean difference in parallax between DR2 and EDR3.}
  \label{fig:mapsLmcPlx}
\end{figure*}

\subsection{Systematic errors}
\label{sec:syst}

Several aspects of the systematic errors (biases) in the astrometric data are 
examined in the EDR3 catalogue validation paper \citep{EDR3-DPACP-126}.
The bias in parallax, and its variation with magnitude, colour, and ecliptic latitude, 
is extensively discussed in a separate paper \citep{EDR3-DPACP-132}. The global
properties of the system of positions and proper motions are discussed in the
EDR3 celestial reference frame paper \citep{EDR3-DPACP-133}. 

In this section we focus on the statistical variations of parallax and proper motion
biases on various angular scales, as revealed by the quasars and (for parallaxes on
scales $\lesssim 1^\circ$) by sources in the direction of the Large Magellanic Cloud
(LMC). We also illustrate the improvements achieved since DR2. The results presented here 
are derived using relatively faint sources ($G\simeq 16$--20), and little is known about 
small- and medium-scale variations at brighter magnitudes, in particular for $G<13$, 
where the sources in many respects behave differently from the fainter sources. Furthermore, 
we only give results for the five-parameters solutions, which are used for most 
sources brighter than $G\simeq 19$ (Fig.~\ref{fig:histNsrc}). 
In general the six-parameter solutions are probably worse than the five-parameter 
solutions in terms of systematics, but it is difficult to know whether this is an intrinsic 
property of the solutions or a consequence of the faintness and more problematic 
nature of most of the sources getting a six-parameter solution (Sect.~\ref{sec:nuEff}). 

Figure~\ref{fig:mapsQso} (left) shows smoothed maps of the parallaxes and proper motion
components for a sample of 1\,215\,942 quasars, namely the subset of sources in 
\textit{Gaia} EDR3 Archive table \gacs{agn\_cross\_id} with five-parameter 
solutions in \gacs{gaia\_source} (median $G=19.9$).
The selection of quasars in \gacs{agn\_cross\_id} is discussed in \citet{EDR3-DPACP-133}.
Smoothed values were computed using a Gaussian kernel of $5^\circ$ standard deviation.%
\footnote{More precisely, the smoothed value at a given point is computed as
the weighted median of the individual values within a radius of $15^\circ$, using
weights proportional to $\exp[-\frac{1}{2}(\theta/5^\circ)^2]$, where $\theta$ is 
the angle between the quasar and the smoothed point.}
The smoothed points in the Galactic zone ($|\,b\,|<10^\circ$) are not displayed, 
as they are dominated by noise from small-number statistics. The standard 
deviations of the smoothed maps (for $|\,b\,|>10^\circ$) are 10.8~$\mu$as in 
$\varpi$, 11.2~$\mu$as~yr$^{-1}$ in $\mu_{\alpha*}$, and 10.7~$\mu$as~yr$^{-1}$
in $\mu_\delta$.

For comparison, we show in the right column of Fig.~\ref{fig:mapsQso} the corresponding 
maps for \textit{Gaia} DR2 astrometry, calculated in the same manner for the 1\,141\,470 of the
sources in the EDR3 quasar sample that have full astrometric data also in DR2. To facilitate 
comparison, the maps use the same colour scales as for the EDR3 data, only shifted by 
$10~\mu$as in parallax to compensate for the different mean
biases. The standard deviations in the DR2 maps are 15.5~$\mu$as,  
26.2~$\mu$as~yr$^{-1}$, and 23.5~$\mu$as~yr$^{-1}$. Thus, in EDR3 the systematics
are reduced by the factors 0.70 ($\varpi$), 0.41 ($\mu_{\alpha*}$), and 0.46 ($\mu_\delta$), 
that is very nearly the same factors as for the random uncertainties (Sect.~\ref{sec:form}).

On much smaller scales, down to $0.1^\circ$, Fig.~\ref{fig:mapsLmcPlx} shows the 
characteristic `checkered pattern' that was very prominent in the DR2 astrometry
for the LMC and in maps of the median parallax in the Galactic bulge 
area (Sect.~4.2 in \citeads{2018A&A...616A..17A}). In EDR3 there is a similar pattern,
but with a different structure and smaller amplitude as shown in Fig.~\ref{fig:mapsLmcPlx}.
The RMS amplitude of the smoothed variations in these plots is $7.7~\mu$as for EDR3 
and $14.3~\mu$as for DR2. 

The maps in Figs.~\ref{fig:mapsQso} and \ref{fig:mapsLmcPlx} were smoothed in order to
bring out clearly the pattern of systematic errors. Although the random errors are strongly
attenuated by the smoothing, they still contribute to the standard deviations quoted above,
which are therefore somewhat higher than the actual RMS systematics on the relevant 
angular scales. In order to correct for this bias, we randomly divided the sources into 
two subsets (A and B) of roughly equal size and computed separate smoothed maps 
$s_\text{A}(\alpha,\delta)$, $s_\text{B}(\alpha,\delta)$ for the subsets. 
Because the random errors are uncorrelated between A and B, while the systematics 
are the same, an unbiased estimate of the mean square systematics is obtained as
the sample covariance between the smoothed values, 
$\text{RMS}=\langle(s_\text{A}-\langle s_\text{A}\rangle)(s_\text{B}
-\langle s_\text{B}\rangle)\rangle^{1/2}$. 
Here $\langle\,\rangle$ denotes an average over the positions (for the quasars,
only $|\,b\,|>10^\circ$ was used; for the LMC, points within a radius of $4.5^\circ$).
Averaging over 50 different random divisions, we obtain the RMS values in 
Table~\ref{tab:rmsSyst}. Compared with DR2 (values in brackets), the RMS
systematics have improved by a factor 0.7 in the quasar parallaxes and 0.44 in the proper 
motions. For the small-scale parallax systematics in the LMC the improvement is a factor 0.53.

The RMS values in Table~\ref{tab:rmsSyst} for the quasars were computed using the full sample 
down to $G=21.0$ (median $G=19.9$), without taking into account that the individual 
uncertainties increase rapidly towards the faint end. This was done in order to benefit maximally 
from the large number of faint quasars in the sample. Unfortunately there are not enough 
of the brighter quasars to determine with any certainty how the systematics depend on 
magnitude, but it appears that the they improve marginally for brighter sources. For example, 
if the sample is restricted to the 16\% quasars brighter than $G=19$ (median $G=18.4$), the 
RMS systematics are 10--15\% smaller than in Table~\ref{tab:rmsSyst}.

For the LMC, the RMS values in Table~\ref{tab:rmsSyst} were computed after subtraction 
of the mean observed parallax in the area, which means that they do not include 
systematics on angular scales $\gtrsim 4.5^\circ$. This explains why the
RMS values in the last line break the increasing trends from the previous lines. 
The magnitude dependence mentioned above could also play a role here, the
LMC sources being on average brighter than the quasars, as well as the geometrically 
favourable location near the south ecliptic pole.


\begin{table}
\caption{RMS level of systematic errors in the quasar and LMC
samples at different angular scales for EDR3 and (in brackets) DR2.  
\label{tab:rmsSyst}}
\small
\begin{tabular}{cccrrrrrr}
\hline\hline
\noalign{\smallskip}
Sample && Angular && \multicolumn{2}{c}{Parallax} && \multicolumn{2}{c}{Proper motion}\\
&& scales &&  \multicolumn{2}{c}{[$\mu$as]} && \multicolumn{2}{c}{[$\mu$as~yr$^{-1}$]} \\
\noalign{\smallskip}
\hline
\noalign{\smallskip}
QSO && $>10^\circ$ && 8.1 & (11.3) && 7.7 & (17.7)\\
QSO && $>7^\circ$ && 8.7 & (12.6) &&  8.8 & (20.1)\\
QSO && $>5^\circ$ && 9.4 & (13.6) && 9.6 & (22.6)\\
QSO && $>3^\circ$ && 10.6 & (14.8) && 10.8 & (24.7)\\
QSO && $>2^\circ$ && 11.6 & (15.5) && 11.6 & (26.3)\\
QSO && $>1^\circ$ && 13.5 & (16.9) && 12.8 & (26.0)\\
QSO && $>0.5^\circ$ && 14.3 & (23.2) && 17.2 & (23.2)\\
LMC && $0.1{-}4.5^\circ$ && 6.9 & (13.1)\\
\noalign{\smallskip}
\hline
\end{tabular}
\tablefoot{The lower limit in the second column is the standard deviation of the Gaussian 
smoothing kernel; the upper limit (for LMC)
is the diameter of the area examined. Columns~3 and 4
give unbiased estimates of the RMS variations of the systematics, computed as 
described in the text. In brackets are the corresponding estimates 
for DR2, using as far as possible the same samples. For the proper motion, 
the values refer to a single component (the mean of the RMS in 
$\mu_{\alpha*}$ and $\mu_\delta$).}
\end{table}

\begin{figure*}
\centering
  \includegraphics[width=0.85\hsize]{fig15.pdf} 
    \caption{Angular covariances of the five-parameter quasar sample. 
    \textit{Left:} Covariance in parallax, $V_\varpi(\theta)$. 
    \textit{Right:} Covariance in proper motion, $V_\mu(\theta)$.
    The red circles are the individual estimates; the dashed black curves are
    fitted exponential functions. The bottom panels show the same data as in the top panels, 
    but for small separations only, with errors bars (68\% confidence intervals) and running 
    triangular mean values (blue curves).}
    \label{fig:qsoCov}
\end{figure*}

Similarly to what was done for DR2, angular covariance functions of the parallaxes and 
proper motions, $V_\varpi(\theta)$ and $V_\mu(\theta)$, have been computed for the
EDR3 quasar sample. See Sect.~5.4 in \citetads{2018A&A...616A...2L} for their definition.%
\footnote{In the DR2 paper these functions were called spatial covariance functions.
However, `angular' is a better qualifier, consistent with the established term
`angular power spectrum' and avoiding the association to (three-dimensional) spatial 
coordinates.} 
The results (Fig.~\ref{fig:qsoCov}) are qualitatively similar to the DR2 results, but
the covariances are smaller by a factor 2 to 4, consistent with other improvements. 
The black dashed curves in the upper panels are exponential fits for 
$0.5^\circ\lesssim\theta\lesssim 80^\circ$,
namely
\begin{align}\label{eq:covPlxExp}
V_\varpi(\theta) &\simeq (142~\mu\text{as}^2)\times\exp(-\theta/16^\circ) \,, \\
\label{eq:covPmExp}
V_\mu(\theta) &\simeq (292~\mu\text{as}^2\text{yr}^{-2})\times\exp(-\theta/12^\circ) \, .
\end{align}
The corresponding amplitudes for DR2 were 285~$\mu$as$^2$ and 
800~$\mu$as$^2$~yr$^{-2}$, with $e$-folding angles $14^\circ$ and $20^\circ$.
Taking the first bin ($0<\theta<0.125^\circ$) to represent the covariance of the 
systematic errors at zero separation, we have
\begin{align}\label{eq:covPlx0}
V_\varpi(0) &\simeq 700~\mu\text{as}^2 \,, \\
\label{eq:covPm0}
V_\mu(0) &\simeq 550~\mu\text{as}^2\text{yr}^{-2} \, .
\end{align}
Corresponding values for DR2 were, respectively, 1850~$\mu$as$^2$ and 
4400~$\mu$as$^2$~yr$^{-2}$. 

Both $V_\varpi(\theta)$ and $V_\mu(\theta)$
show oscillations with a period of the order of a degree, corresponding to
the checkered pattern. Consistent with the other findings, the oscillations in 
Fig.~\ref{fig:qsoCov} have significantly smaller amplitudes than their 
counterparts in DR2. 



\subsection{Angular power spectrum of parallax bias}
\label{sec:aps}

A comprehensive quantification of the positional variations of systematics on all 
angular scales can be given in the form of an angular power spectrum. This section
is an attempt to estimate the power spectrum of parallax bias from EDR3 quasar data.

In astrophysics, the angular power spectrum is 
perhaps best known in the context of the cosmic microwave background
(CMB; e.g.\ \citeads{2002ARA&A..40..171H}). Any scalar field $z(\alpha,\delta)$ 
defined on the full sky
(temperature for the CMB; or in our case, the quasar parallaxes) can be 
decomposed in spherical harmonics (SH) $Y_{\ell m}(\alpha,\delta)$ as
\begin{equation}\label{eq:SH1}
z(\alpha,\delta) = \sum_{\ell=0}^\infty \sum_{m=-\ell}^\ell a_{\ell m}Y_{\ell m}(\alpha,\delta)\, ,
\end{equation}
where $\ell$ is the degree of the SH (also known as multipole), $m$ is the order (or mode),
and $a_{\ell m}$ are the (complex) coefficients
\begin{equation}\label{eq:SH2}
a_{\ell m} = \int_\Omega z(\alpha,\delta) Y_{\ell m}^*\text{d}\Omega\, .
\end{equation}
Here $^*$ is the complex conjugate and $\int_\Omega$ denotes integration over the full
sphere with solid angle element $\text{d}\Omega=\cos\delta\,\text{d}\alpha\,\text{d}\delta$.
The equivalence of Eqs.~(\ref{eq:SH1}) and (\ref{eq:SH2}) can be verified using  
the orthonomality of the SH, 
$\int_\Omega Y_{\ell m}Y_{\ell' m'}^*\text{d}\Omega=\delta^{\ell\ell'}\delta^{mm'}$,
where $\delta^{ij}$ is the Kronecker symbol. By means of Eq.~(\ref{eq:SH1}) it is seen that 
the mean square value of $z$ on the sky is
\begin{equation}\label{eq:SH3}
\begin{split}
\overline{z^2}&=\frac{1}{4\pi}\int_\Omega |z|^2\text{d}\Omega =
\frac{1}{4\pi}\int_\Omega \left(\sum_{\ell,\, m} a_{\ell m}Y_{\ell m}\right)
\left(\sum_{\ell,\, m} a_{\ell m}^*Y_{\ell m}^*\right)\text{d}\Omega\\ 
&=\frac{1}{4\pi}\int_\Omega \sum_{\ell,\, m}|a_{\ell m}|^2 |Y_{\ell m}|^2\text{d}\Omega 
= \frac{1}{4\pi}\sum_{\ell=0}^\infty \sum_{m=-\ell}^\ell |a_{\ell m}|^2\, .
\end{split}
\end{equation}
The observed angular power spectrum, defined as 
\begin{equation}\label{eq:SH4}
C_\ell = \frac{1}{2\ell+1}\sum_{m=-\ell}^\ell |a_{\ell m}|^2\, 
\end{equation}
(\citeads{1973ApJ...185..413P}; \citeads{2003ApJS..148..135H}), thus measures the mean power
of the SH components of degree $\ell$, that is on angular scales 
\hbox{$\sim\! 180^\circ/\ell$}, integrated over the sphere. 
For our purpose it is convenient to plot the cumulative quantity
\begin{equation}\label{eq:SH5}
R(\ell_\text{max}) = \sqrt{\frac{1}{4\pi}\sum_{\ell=1}^{\ell_\text{max}} (2\ell+1)C_\ell}\, .
\end{equation}
The spherical harmonic of degree $\ell=0$, corresponding 
to the mean quasar parallax of about $-20~\mu$as, is not included in the sum. 
$R(\ell_\text{max})$ can therefore be interpreted as the RMS variation of the parallax 
systematics on angular scales $\gtrsim 180^\circ/\ell_\text{max}$.



\begin{figure}
\centering
  \includegraphics[width=\hsize]{fig16.pdf} 
    \caption{Cumulative angular power spectrum of systematics in quasar parallaxes.}
    \label{fig:caps}
\end{figure}



Figure~\ref{fig:caps} presents our estimates of the angular power spectrum of the
parallax bias in EDR3, derived from the five-parameter quasar sample using several different 
methods. Most straightforward is to determine the coefficients $a_{\ell m}$ by a weighted 
least-squares fit of Eq.~(\ref{eq:SH1}) directly to the quasar parallaxes, truncating the first 
sum at $\ell=L$ (that is, using $L^2$ unknowns). Owing to the lack of quasars at low Galactic 
latitudes ($|\sin b|<0.1$), this gives stable results only for $L\lesssim 15$ (angular scales 
$\gtrsim 12^\circ$), where the fit manages to bridge the no-data gap. 
Even so, the RMS values computed in this way (shown as open circles in
Fig.~\ref{fig:caps}) overestimate the true variations, as they include a contribution from
the random errors in the parallaxes. If the random errors are assumed to have standard
deviations equal to the formal parallax uncertainties, we can estimate their
contribution to $C_\ell$ by means of Monte Carlo simulations, and subtract from the 
power obtained in the fit. (Alternatively, the noise contributions could be estimated from
the formal variances obtained in the least-squares fits.) 
This gives the corrected RMS estimates shown as filled circles in Fig.~\ref{fig:caps}. 

The smoothed maps offer an alternative method to estimate the angular power spectrum
that is not restricted to $\ell_\text{max}\lesssim 15$, as was the case for the SH fit
owing to the no-data gap along the Galactic equator. The RMS values in Table~\ref{tab:rmsSyst}
were computed using only the smoothed points with $|\sin b|>0.1$, and are therefore 
not strongly affected by the gap.
However, the Gaussian smoothing does not correspond to a well-defined cut-off
in $\ell$, and the comparison of the smoothed RMS values with 
$R(\ell_\text{max})$ is not entirely straightforward. In Fig.~\ref{fig:caps} we have put the
RMS values from the table at the degree where the Gaussian beam transfer function,
$\exp\bigl[-(\ell+\frac{1}{2})^2\sigma^2\bigr]$ \citepads{1992PhRvD..46.4198W},
equals $\frac{1}{2}$, that is at $\ell_\text{max}\simeq 47.7^\circ/\sigma-0.5$.
The RMS values, shown as green squares, roughly continue the power-law trend of the 
corrected SH fit to higher $\ell_\text{max}$, but at RMS values that are 10--15\% higher. 
We have no explanation for this discrepancy, but conclude that the agreement is 
reasonable considering the approximations involved.  

The angular covariance function $V_\varpi(\theta)$ and the angular power spectrum $C_\ell$
contain equivalent information, and are related by the transformations \citepads{1973ApJ...185..413P}
\begin{align}
C_\ell &= 2\pi \int_{-1}^1 V_\varpi(\theta)\,P_\ell(\cos\theta)\,\text{d}(\cos\theta)\, ,\label{eq:SH6}\\
V_\varpi(\theta) &= \frac{1}{4\pi}\sum_\ell (2\ell+1)\,C_\ell\,P_\ell(\cos\theta)\, ,\label{eq:SH7}
\end{align}
where $P_\ell(x)$ are Legendre polynomials. Using the angular covariance function in
Fig.~\ref{fig:qsoCov} (left) and replacing the integral in Eq.~(\ref{eq:SH6}) by a sum over the
1440 covariance values, one can easily compute $C_\ell$ for arbitrary $\ell$. The result
is the red solid curve in Fig.~\ref{fig:caps}. Using the smoothed covariance values (blue
in Fig.~\ref{fig:qsoCov}, left) gives instead the black curve in Fig.~\ref{fig:caps}. Both curves
show some unphysical fluctuations: The cumulative RMS values cannot decrease for increasing
$\ell$. The fluctuations are caused by sampling noise and must be disregarded when 
interpreting the RMS values, although they give an impression of the statistical uncertainties
in $R(\ell_\text{max})$.

After accounting for estimation bias in the SH fit there is generally good agreement between
the different methods where they overlap. On the smallest scales, $\ell_\text{max}\gtrsim 100$, 
only the angular covariance function provides an estimate. As indicated by the dashed line 
in Fig.~\ref{fig:caps}, the overall trend can be described by a simple power law 
\begin{equation}\label{eq:SH8}
R(\ell_\text{max}) \simeq (5.8~\mu\text{as})\times \ell^{\,0.18}\, ,
\end{equation}
at least for for $\ell=3$ to ${\simeq\,}150$ (angular scales from 1.2$^\circ$ to 60$^\circ$).
For $\ell=1$ and 2 the RMS is significantly smaller than according to this relation. In particular, 
the power at $\ell=2$ (angle $\sim 90^\circ$) is remarkably small, which could be related to
the basic angle ($\Gamma=106.5^\circ$) providing a firm connection between areas
separated by angles of the order of $90^\circ$. For $\ell \gtrsim 150$ the RMS is higher 
than according to the
power law, corresponding to the `checkered pattern' in Fig.~\ref{fig:mapsLmcPlx} (left).
The value $V_\varpi(0) \simeq 700~\mu\text{as}^2$ in Eq.~(\ref{eq:covPlx0}) for
$\theta<0.125^\circ$ suggests that $R(\ell_\text{max})$ saturates at $\simeq 26~\mu$as 
for $\ell_\text{max}\gtrsim 1440$. Using a suitable (monotonic) model function
$R(\ell_\text{max})$ one can estimate the values $C_\ell$, which may be of interest 
in studies where the statistical variation of parallax bias with position is a concern.

No angular power spectrum is given for the systematics in quasar proper motions.
Large-scale systematics (for small $\ell$) are discussed in \citet{EDR3-DPACP-133} 
and \citet{EDR3-DPACP-134}.


\section{Improvements for \textit{Gaia} DR4 and beyond}
\label{sec:impr} 

Although EDR3 brings huge improvements over DR2 in terms of the overall quality of the astrometric
results, including systematics, it is obvious that the limits of \textit{Gaia}'s capability have not
been reached. At the time of writing (September 2020), \textit{Gaia} has already accumulated more 
than twice the amount of observations included in EDR3, and a solution using all these data, even 
without any improved modelling, will almost certainly bring down the random errors by a factor 0.7 for
the positions and parallaxes, and by a factor 0.35 for the proper motions. Some of the
systematics will also be reduced simply from the improved coverage in time and scan directions. 
But important advances will also come from dedicated efforts to improve and consolidate
the calibration models. In this context it is positive that many model deficiencies stand 
out clearly in the residuals (e.g.\ Appendix~\ref{sec:residT} and \ref{sec:residG}): It shows 
that the modelling errors are not degenerate with the astrometry and can be used to design 
better models. Below we list some areas where significant advances should 
be possible as a result of further model developments and analysis of the data.


\subsection{LSF and PSF modelling}
\label{sec:imprLSF}
 
The processing of CCD samples in CALIPD is intimately connected to the AGIS calibrations through 
the many instrumental effects that influence the shape and location of the image profiles. 
As detailed by \citet{EDR3-DPACP-73}, significant improvements to the LSF and PSF modelling are 
being implemented or planned for DR4 and subsequent cycles. These include improved basis 
functions for the LSF and PSF modelling, with an analytical representation of AL and AC smearing 
effects and a clear separation between parameters representing a shift of the profile and 
its shape; modelling of magnitude-dependent non-linear effects, caused for example by CTI; 
and a new bootstrapping of the attitude and geometric calibrations that will make the initial
CALIPD in a cycle more independent of the previous cycle, thus reducing the risk of 
propagating systematics from one cycle to the next. It is expected that these many developments
will further reduce systematics in the location and flux estimates per CCD coming from the IPD.
This will benefit the astrometry and photometry of all sources, but in particular the bright 
ones ($G\lesssim 13$). 


\subsection{CTI modelling}
 
As shown in Appendix~\ref{sec:residG}, the AL and AC residuals exhibit strong trends with
magnitude that can be interpreted as a manifestation of CTI. CTI effects are also diagnosed, 
but not corrected, by the CTI calibration parameters displayed in Fig.~\ref{sec:residG}.
The effects may already have some impact on the astrometry, and this will become more
important with time as the accumulated radiation dose increases while random errors
in the astrometric results continue to improve. Whereas ultimately a detailed, physics-based 
modelling of CTI effects is desirable (\citeads{2011MNRAS.414.2215P}; \citeads{2013MNRAS.430.3078S}),
much progress is still possible simply by better mapping of the empirical effect as a function 
of the main variables of the observation. On this macroscopic level, CTI is comparatively easy 
to separate from other effects, thanks to the periodic charge injection, and a wealth of additional 
calibration data are available to support the modelling \citepads{2016A&A...595A...6C}. 

\subsection{Time variations}

With the exception of the spin-related (quasi-periodic) distortion discussed in Sect.~\ref{sec:glob},
all instrument calibrations are assumed to be constant or at most linearly varying over a time interval
that could be as long as 63~days. For the AL large-scale geometry the time interval is at most 3~days.
The residual normal points in Fig.~\ref{fig:rnpAllTimeNoGc1} show that this resolution is adequate 
most of the time, but not always and especially not after thermal upsets such as the (partial) eclipses
of the Sun by the Moon. Clearly a quadratic model (or perhaps an exponential model with fixed 
time constant) 
would be a big improvement at these times. The spin-related periodicity of the residuals, visible 
along most of the time axis in the diagram, may be removed by the more extended VBAC/FOC 
modelling hinted at in Sect.~\ref{sec:glob}.

\subsection{AC rate dependency}
\label{sec:acRate}

The AC rate has a big impact on the PSF by smearing the images in the AC direction during a 
CCD exposure. The width of the PSF in the AC direction is increased by the smearing, and the 
number of photoelectrons per pixel is reduced, which affects saturation and CTI. The AC rate is 
therefore relevant for the calibration of both bright and faint sources. To first order, the smearing 
is proportional to $|\,\dot{\zeta}\,|^2$, which is the only dependence considered in the AGIS calibration 
model for EDR3 (effect~7 in Table~\ref{tab:cal}). Owing to the non-linearity of saturation and CTI 
effects, it is likely that the AL centroid biases they produce also have a component that is linear in
$\dot{\zeta}$. Such terms were not included in the AL calibration model for EDR3, as they might 
be difficult to disentangle from the parallax in view of the correlation between $\dot{\zeta}$ and 
the AL parallax factor described in Appendix~\ref{sec:plxAc}. This correlation is positive throughout 
the data segments used for EDR3, but negative in data segments DS6 and DS7 where reversed 
precession was used. The reversed precession during one year (2019.536--2020.576; 
cf.\ Fig.~\ref{fig:timeline}) was introduced precisely to address this and similar issues related 
to the non-symmetry of the nominal scanning law. Together with the much improved 
PSF modelling mentioned above, this should allow the main effects of AC rate variations to 
be resolved already in the data analysis for DR4, which includes data segment DS6 obtained 
in reversed precession mode.

\subsection{Use of colour information}

A rather unsatisfactory aspect of EDR3 is the division of sources with full astrometric information 
into two distinct subsets, namely the five- and six-parameter solutions. This was a consequence
of the unavailability of good colour information for some sources, which necessitated
their special treatment in IPD and AGIS. Although much more and better colour information
from the BP and RP spectra is available for DR4 through PhotPipe~3 (cf.\ Fig.~\ref{fig:timeline}), 
it is unavoidable that this information is missing or of poor quality for some sources. A more uniform 
treatment of all sources in IPD and AGIS can be achieved by consistently using the available
colour information, weighted according to its uncertainty. For AGIS this means that all 
sources obtain six-parameter solutions, but with available BP and RP data used as prior for the 
pseudocolour. 

\section{Conclusions}
\label{sec:concl}

Compared with \textit{Gaia} DR2, the number of sources in EDR3 that have
a parallax and proper motion is only 10\% higher. However, the average improvement 
on the standard uncertainties is roughly a factor 0.8 for the positions and parallaxes, 
and 0.5 for the proper motions. These factors reflect the higher number of observations
per source, by more than 50\% on average, and the longer time span of the data, which 
make the astrometric results 
considerably more robust and help to reduce systematic errors. The astrometric solution 
for EDR3 is also the
first one in the cyclic processing of DPAC to benefit from a full reprocessing of
the LSF and PSF calibrations and the image parameter determination. 
The next full-scale astrometric solution, for \textit{Gaia} DR4,
will be based on twice as many observations as EDR3. Considerable efforts are required
and planned to ensure a matching development of models and analysis methods. 


\input{acknowledgements}

\bibliographystyle{aa}
