parameters are also reported in Gaia DR3. This is expected, because stellar
parameters were inferred for sources independently of their classification
status (Creevey et al., 2022). There are 255 948 sources in the qso_candidates
table and 7069 sources in the galaxy_candidates table that have effective
temperatures derived by the CU8 GSP-Phot module (Fousneau et al., 2022).
Checking a variety of metrics such as magnitude, sky distribution, and
effective temperature itself, there is nothing apparently peculiar with these
sources. Their presence is an inevitable consequence of the known stellar
contamination. It is also important to remember that DSC, which is the single
largest contributor to these integrated tables, did not filter out sources
simply because they were bright (only DSC-Allosmod classifies sources with
$G<14.5$ mag to be stars).

There are also 4027 sources with valid radial velocities in the qso_candidates
table, and 160 in the galaxy_candidates table. Considering that the
extragalactic tables are mostly populated with faint sources, these small
numbers are essentially due to the intrinsic magnitude limit of sources for
which radial velocities could be derived in Gaia DR3 (Katz et al., 2022).
Those featuring valid radial velocities have magnitudes that are usually
incompatible with extragalactic sources, so it is fair to assume that they are
stars.

###  5.4 Astrometric selection

Additional insight into the classified sources can be gained by analysing
their astrometric parameters. As has been demonstrated in Gaia Collaboration &
Klioner et al. (2022) and Gaia Collaboration et al. (2021), astrometry can be
used to improve the purity of a sample of quasar candidates. It is clear,
however, that this can only be achieved at the cost of reducing the
completeness.

The procedure here is similar to that used in the construction of Gaia-CRF3,
namely a two-step astrometric filtering of a sample of candidates (Gaia
Collaboration & Klioner et al., 2022). In the case of Gaia-CRF3, the sample
was obtained by cross-matching the Gaia EDR3 catalogue with several external
quasar catalogues. In the present study, each of the Gaia classifiers
contributing to the qso_candidates table is considered as an additional
catalogue, and the same procedure is applied to all the external and Gaia-own
selections of quasar candidates.

The first step of the astrometric filtering is to select individual sources
that have high-quality astrometric solutions in Gaia EDR3 and statistically
insignificant parallaxes and proper motions (see Gaia Collaboration & Klioner
et al., 2022, Sect. 2.1 for the exact mathematical formulations). This step
alone is insufficient to find genuine quasars (or extragalactic objects), as
about 214 million sources in Gaia EDR3, dubbed ’confusion sources’, satisfy
these astrometric criteria. These are mostly stars of our Galaxy (Gaia
Collaboration & Klioner et al., 2022, appendix C). At least at this stage of
the Gaia project, astrometry cannot be used as an independent quasar
classifier, although this may change in the future (see e.g. Heintz et al.,
2015, 2018).

A second step of filtering is therefore needed. In this step, only those
samples of sources are retained that show near-Gaussian distributions in the
uncertainty-normalized parallaxes and proper motions. Since extragalactic
sources are faint, the typical uncertainties of their astrometric parameters
in Gaia DR3 are about two orders of magnitude larger than either the known
level of systematic errors in Gaia DR3 (Lindegren et al., 2021a) or the known
physical systematic effects (Gaia Collaboration et al., 2021). Bearing in mind
that the true parallaxes and proper motions of genuine extragalactic sources
should be zero, one expects Gaussian distributions of the normalized
parameters. This requirement had proven to be very useful to distinguish
genuine quasars from the confusion sources.

Both steps of the astrometric filtering obviously reject some genuine quasars
that have considerable measured, but spurious, proper motions due to time-
varying source structure (see Sect. 2). A prominent example here is 3C273,
which is not part of Gaia-CRF3 for this reason. The samples considered at the
second step of the astrometric filtering can come from a particular external
catalogue or from one of the Gaia classifiers, but could also be selections
according to various criteria (e.g. avoiding the crowded areas on the sky) or
intersections of such selections (e.g. sources that were found to be quasars
by two classifiers).

An additional characteristic of a sample of genuine extragalactic objects is
that its sky distribution should not show overdensities in known stellar
structures in our Galaxy and its environments, such as clusters, although it
could still be influenced by such structures, for example variable Galactic
extinction. This can also be used to help decide whether a particular sample
of sources should be retained.

Figure 23: Distribution Galactic coordinates (Hammer--Aitoff projection) of
the 1 897 754 sources from the astrometric selection (i.e. sources in the
qso_candidates table with astrometric_selection_flag set). The plot shows the
density of sources per square degree computed from the source counts per pixel
at HEALPix level 7 (pixel size $\simeq 0.21$ deg2)

Using this two-step selection procedure we have identified a set of 1 897 754
quasar candidates, which we refer to as the ’astrometric selection’. They are
indicated by the astrometric_selection_flag in the qso_candidates table. The
purity of this sample is difficult to estimate, but we believe it to be 98% or
perhaps better. The vast majority of these sources were identified as quasars
by at least two independent external catalogues and/or Gaia classifiers. The
density distribution of these sources on the sky is shown on Fig. 23. This set
contains 1 406 729 sources (74%) with 5p astrometric solution and 491 025
sources (26%) with 6p solutions in Gaia DR3. The avoidance zone in the
Galactic plane as well as the lower density of sources around the LMC and SMC
result from the difficulty in reliably identifying quasars in those crowded
areas. This concerns both the external catalogues and the Gaia classifiers.

Figure 24: Distributions of the normalized parallaxes and proper motion
components for the sources in the astrometric selection with 5p (blue) and 6p
(green) solutions. The red curves show the corresponding best-fit Gaussian
distributions. The global parallax zero point of $-0.017$ mas of Gaia DR3 is
taken into account (Lindegren et al., 2021b, a). The standard deviations of
the best-fit Gaussian distributions for the sources with 5p (6p) solutions are
1.048 (1.068), 1.054 (1.092) and 1.063 (1.109) for the parallaxes, and proper
motions in right ascension and declination, respectively. As usual in Gaia,
the asterisk in $\alpha*$ in the middle panel indicates the implicit factor
$\cos\delta$, that is $\mu_{\alpha*}=\dot{\alpha}\,\cos\delta$.

Fig. 24 shows the distributions of the normalized parallaxes and proper
motions of the astrometric selection. They are close to Gaussian, which
suggests a reasonably low level of stellar contamination. The standard
deviations of the best-fit Gaussian distributions range from $1.05$ to $1.11$
and indicate by how much the formal uncertainties of the corresponding
astrometric parameters may be underestimated in Gaia DR3.

We attempted a similar astrometric selection for the galaxy_candidates table.
However, since most of the galaxies have only two-parameter astrometry in Gaia
DR3, and as more problems in Gaia astrometry can be expected for extended
sources, the astrometric selection for the galaxy table turned out to be less
useful, so we decided not to publish it. Nonetheless, this analysis did reveal
the properties of the population of sources in the astrometric selection that
were classified as both quasars and galaxies: the astrometric selection from
the qso_candidates table contains 54 892 sources that are also present in the
galaxy_candidates table (cf. overall overlap of these tables of 174 146
sources). 99% of those sources have 6p astrometric solutions. The normalized
parallaxes and proper motions of this set of sources also have near-Gaussian
distributions, but with standard deviations of 1.13--1.25, which is about 10%
larger than for the astrometric selection as a whole. This set of sources is
probably dominated by AGN for which source structure (i.e. the host galaxy)
notably affects the astrometric solution. Indeed, a host galaxy was detected
by Gaia for 23 805 of these sources (43%). Similar statistics of the
normalized astrometric parameters can also be found for the set of sources in
the astrometric selection for which a host galaxy was detected by Gaia (see
Sect. 4.3), which contains 51 586 sources.

Thus we encounter the problem in the optical that is well known in radio
astrometry (e.g. Charlot et al., 2020), namely the influence of the source
structure on the quality of the astrometry. This topic will need a special
attention in the future Gaia data releases.

###  5.5 Analysis of objects with lower probability classifications

Figure 25: OA class labels for extragalactic objects. HQN = high quality
neuron (quality 0--3), LQN = low quality neuron (quality 4--6).  Figure 26:
Distribution of OA neurons labelled as extragalactic for each quality
category. The number on each bar gives the number of sources in the
qso_candidates or galaxy_candidates tables. No extragalactic objects appear in
any of the best (0) or worst (6) quality neurons.

The unsupervised algorithm OA was used to analyse the sources with lower DSC
class probabilities (Sect. 2.4). Here we focus on those neurons that were
assigned to an extragalactic class label (QSO or GAL). These are shown in Fig.
25 for two different subsets: high quality neurons (HQN), that represent
quality categories 0 to 3, and low quality neurons (LQN), that represent
categories 4 to 6. We further limit our analyses to those sources that appear
in the integrated tables. Figure 26 shows the number of neurons and objects
assigned to each quality category. Approximately 80% of the sources are
assigned to a high quality neuron.

Table 6: Contingency table for OA classifications. Each entry gives the percentage of objects classified by DSC (using classlabel_dsc), or processed by QSOC or UGC, that are assigned to OA high-quality neurons (HQN) or low-quality neurons (LQN), for sources in the qso_candidates and galaxy_candidates tables.  | OA |
---|---|---
| QSO | GAL | other |
| HQN | LQN | HQN | LQN | HQN | LQN | Total
DSC | quasar |  $20$% |  $15$% |  $1$% |  $1$% |  $45$% |  $18$% | $2\,158\,916$
galaxy |  $7$% |  $2$% |  $77$% |  $6$% |  $3$% |  $5$% | $851\,127$
other |  $7$% |  $8$% |  $53$% |  $1$% |  $24$% |  $7$% | $1\,993\,592$
QSOC | quasar |  $19$% |  $15$% |  $2$% |  $1$% |  $45$% |  $17$% | $3\,069\,458$
other |  $3$% |  $1$% |  $87$% |  $3$% |  $4$% |  $2$% | $1\,934\,177$
UGC | galaxy |  $3$% |  $1$% |  $83$% |  $6$% |  $2$% |  $5$% | $199\,093$
other |  $13$% |  $10$% |  $33$% |  $2$% |  $30$% |  $12$% | $4\,804\,542$

All OA sources were processed by DSC, as well as QSOC or UGC depending on
their DSC probabilities. Table 6 is a contingency table showing the fraction
of objects in common between these classifiers and the various OA neurons. For
this we use classlabel_dsc. QSOC and UGC do not classify sources; for this
purpose we just look at the ones they provide redshifts for. Among the
galaxies identified by DSC, 83% of them were also found to be galaxies by the
OA module, of which 77% landed in a high quality neuron and 6% in a low
quality one. The coincidence increases for UGC, with 89% of its galaxies found
in a galaxy neuron, of which 83% have high quality. The coincidence for the
quasars is substantially lower, around 35% for both DSC and QSOC, with no
substantial difference between high and low quality neurons. We also see that
a large fraction of those objects that were not classified as a quasar or
galaxy by DSC, or that were not analysed by QSOC, are classified as galaxies
by OA: 54% and 90%, respectively, where most of them belong to a high quality
neuron.

OA processes sources that tend to be faint with noisy ${\rm BP/RP}$ spectra,
some of which OA had to modify (e.g. remove negative fluxes) so that it could
process them. Table 6 suggests that the OA classification complements the
results from the other modules. OA coincides with DSC and UGC when identifying
galaxies in particular, and identifies objects rejected by those modules that
may be real galaxies. OA could also potentially help to identify extragalactic
candidates that are not in the qso_candidates or galaxy_candidates tables.

##  6 External comparison

###  6.1 WISE and proper motions

Figure 27: Gaia-catWISE colour--colour diagrams. Top: All sources in the
qso_candidates table. Bottom: All sources in the galaxy_candidates table. In
both cases regions around the LMC and SMC have been excluded. The colour scale
shows the density of sources on a log scale relative to the peak density
(densities 1000 times lower than the peak are not shown).  Figure 28: Gaia
proper motion vs catWISE W1$-$W2 colour for all sources in the qso_candidates
table (top) and the subset with classlabel_dsc_joint = quasar (bottom).
Regions around the LMC and SMC are excluded. The colour scale shows the
density of sources on a log scale relative to the peak density (densities 1000
times lower than the peak are not shown).

To investigate the infrared colours of the sources in the integrated tables,
we cross-matched them to the catWISE2020 catalogue (Marocco et al., 2021,
including the 2021 catalogue updates) using a 1″matching radius. We found 4.31
million sources (65%) matches in the qso_candidates table, and 4.59 million
(95%) in the galaxy_candidates table. Excluding the regions around the LMC and
SMC (defined in appendix B) left 2.99 million matches in the qso_candidates
table and 4.46 million in the galaxy_candidates table. Figure 27 shows the
distribution of these sources in a Gaia-catWISE colour--colour diagram. We see
that most galaxy candidates have W1$-$W2 colours between 0.0 and 0.5 mag. This
agrees with the range identified by Stern et al. (2012) for galaxies without
an active nucleus and redshifts below 0.6. The quasar candidates, in contrast,
show two overdensities in the catWISE colour. We explore this further by
looking at the quasar candidates in the proper motion space, as shown in Fig.
28. The upper panel is for all sources with 5p or 6p solutions (2.87 million
sources). We see that the bluer clump at around W1$-$W2 $\simeq$ 0 mag shows
the full range of proper motions. Recall that non-zero proper motions of true
quasars are spurious, either due to noise or to time-variable source
structure. Nonetheless, the larger proper motions in the bluer clump compared
to the redder clump is indicative of contamination by stars (and some
galaxies), and the W1$-$W2 colour would seem to confirm that. Indeed, the
lower panel of Fig. 28 is for the purer subset defined by classlabel_dsc_joint
= quasar (0.50 million sources), and this retains just the redder sources with
proper motions that are more consistent with zero (plus noise).

###  6.2 Quasars

Here we look in more detail at the properties of known quasars in Gaia. For
this purpose we cross-matched quasars from SDSS-DR14 (Pâris et al., 2018) that
have a visually confirmed redshifts (source_z = VI or source_z = DR7Q) to all
Gaia sources (those in gaia_source) using a 1″ matching radius. Such a match
is nominally identical to the quasars selected for training DSC (see
Delchambre et al. 2022). However, we further limited this set to those with
complete data in all photometric bands, at least five observations in BP and
RP, complete astrometry (i.e. 5p or 6p solutions), and with $G<20.75$ mag.
This gave 232 794 sources covering the redshift range 0.038--5.305.

A quasar in SDSS-DR14 is defined according to spectroscopic criteria.
Specifically, they are sources with: (a) either at least one broad emission
line with a full width at half maximum larger than 500 $\textrm{km\,s}^{-1}$,
or interesting or complex absorption features; and (b) sufficiently large
intrinsic luminosity ($M_{\rm i}[z=2]<-20.5$). Since only one broad emission
line is required, some objects may otherwise be classified as type 2 AGNs
(those with predominantly narrow emission lines). The second part of the first
condition aims to include Broad Absorption Line (BAL) quasars. This definition
is free of morphological criteria.

Figure 29: Completeness of classlabel_dsc = quasar (top) and
classlabel_dsc_joint = quasar (bottom) with respect to the SDSS-DR14Q-Gaia DR3
cross-match as a function of $G$ and SDSS redshift. Empty bins (white) have
fewer than 25 sources.

The sample defined above is similar to the superset from which the DSC-
Allosmod training set was drawn. However, DSC did not force classifications on
them, so we can use it to assess DSC’s completeness as a function of magnitude
and redshift (further assessments can be found in Bailer-Jones 2021). This is
shown in Fig. 29, using the two class labels from DSC. The dependence on
redshift is expected because of its weak correlation with $G_{\rm
BP}\\!-\\!G_{\rm RP}$, which increases the confusion with stars at high
redshifts and with galaxies at low redshifts. Nonetheless, the completeness is
above 80% for redshifts between 0.3 and 3.6 and $G$ $\leq$ 20.25 mag. The
lower completeness at fainter magnitudes is also expected, because lower
quality data are more likely to be classified by DSC as the majority class of
stars, according to the global prior (Sect. 2.1), especially for the more
conservative classlabel_dsc_joint label. This also explains why the overall
completeness is much lower for this label, although it is still above 60% from
$z=0$ to $z=2.5$ for $G$ $\leq$ 19.25 mag. The overall completeness of
classlabel_dsc is 215 721/232 794=93% and of classlabel_dsc_joint is 97
995/232 794=42%. However, given the non-uniform selection function of SDSS for
obtaining spectra, we should be careful not to over-interpret this specific
assessment of the DSC’s completeness.

Figure 30: $G_{\rm BP}\\!-\\!G_{\rm RP}$ vs redshift for the subsets
classlabel_dsc (top) and classlabel_dsc_joint (bottom) of the SDSS-DR14Q-Gaia
DR3 cross-match. The orange points are classlabel_dsc = galaxy (top) and
classlabel_dsc_joint = galaxy (bottom). The middle panel shows the quasars
with classprob_dsc_combmod$>$ 0.5 for any of the three stellar classes used in
DSC-Combmod. The red curves show the median $G_{\rm BP}\\!-\\!G_{\rm RP}$ and
the 16% and 84% quantiles for all sources in the cross-match, regardless of
the DSC class (so are the same in all panels).

The $G_{\rm BP}\\!-\\!G_{\rm RP}$ vs redshift relation for the sources
classified as quasars by classlabel_dsc and classlabel_dsc_joint is shown in
Fig. 30 (top and bottom panels). The pattern of undulations is expected, and
is due to the quasar emission lines moving across the bands with redshift. The
tail to the red for $z>3.5$ corresponds to the Ly$\alpha$ forest entering and
then filling the BP band. We see that some low redshift objects classified by
SDSS as quasars are classified as galaxies by DSC. These quasars likely have a
higher contribution of the host galaxy to the total emission, making them
redder. The quasars with classlabel_dsc_joint = quasar follow neatly the
bluest part of the colour-$z$ relation, avoiding the regions with $G_{\rm
BP}\\!-\\!G_{\rm RP}$ above the median (compare top and bottom panels). This
result and Fig. 29 indicate that this class selects mainly bright quasars, and
from these only the bluest ones. classlabel_dsc = quasar complements the
classlabel_dsc_joint = quasar class by covering the quasars that are redder
due to their intrinsic emission, Galactic extinction, or local absorption (as
in BALs, for example). The middle panel of Fig. 30 shows that the
incompleteness of classlabel_dsc at high-$z$ is due to the misclassification
of many of these quasars as stars (see also Fig. 29). Similarly, this plot
shows that the quasars in the envelope of the reddest colours over the range
$z=$ 0.5--3.5 are also classified as stars.

Figure 31: Colour--colour diagram for the set in Fig. 30 colour-coded by
redshift (top) and Gaia phot_bp_rp_excess_factor (bottom). The straight lines
correspond to fixed $G_{\rm BP}\\!-\\!G_{\rm RP}$ colour, with values 2, 1.5,
1, 0.5 and 0.

Figure 31 shows the colour--colour diagram for the sample colour-coded by SDSS
redshift (top panel) and by Gaia phot_bp_rp_excess_factor (bottom panel). In
the upper panel we see a clear trend of colour with redshift, with low
redshifts located in the upper left and redshift increasing as we descend, but
with the highest redshifts on the far right. The bottom panel shows that the
region of low-$z$ sources corresponds to sources with larger
phot_bp_rp_excess_factor, which indicates that the combined BP and RP bands
contain more flux than the $G$ band. This region overlaps with the location of
the galaxies (see plots of the DSC source densities in Delchambre et al.
2022). This suggests that the excess in $G_{\rm BP}$ and $G_{\rm RP}$ with
respect to $G$ is due to the wider photometric windows of $G_{\rm BP}$ and
$G_{\rm RP}$ compared to $G$, which allows detection of the quasar’s host
galaxy emission over a wider region than in the $G$ band. This, together with
the red $G_{\rm BP}\\!-\\!G_{\rm RP}$ colours, indicates a shift from quasars
dominated by the nucleus to quasars with an important contribution of the host
galaxy. This transition can be also appreciated in the colour--colour diagram
of the purer sub-samples from the qso_candidates and galaxy_candidates tables
(Fig. 37), followed by the transition to the general galaxy population.

###  6.3 Galaxies

Table 7: SDSS spectral classes for objects in the galaxy_candidates table. The first row gives the number of sources found in SDSS for the whole table (first column) and for each module (subsequent columns). The following rows give the percentage of sources of each SDSS class among these. The GALAXY row can be thought of as a measure of the purity against sources with SDSS spectral classifications, which by design is dominated by extragalactic sources, so is an overestimate of the purity of a sample selected at random from Gaia.  SDSS | galaxy_ | classlabel_ | classlabel_ | Vari- | OA | UGC | Surface |
---|---|---|---|---|---|---|---|---
| candidates | dsc | dsc_joint | Classification |  |  | brightness |
TOTAL | 534 154 | 477 939 | 48 460 | 360 217 | 105 296 | 248 196 | 96 918 | #
GALAXY | 98.0 | 97.9 | 95.1 | 99.7 | 95.6 | 97.9 | 99.7 | %
QSO | 1.6 | 1.7 | 4.9 | 0.3 | 3.9 | 2.0 | 0.3 | %
STAR | 0.4 | 0.4 | 0.1 | 0.0 | 0.5 | 0.0 | 0.0 | %

We compared the content of the galaxy_candidates table with the spectral
classes in SDSS DR16. We cross-matched the catalogues with a 1 arcsec radius
and removed duplicated sources or ones where zWarning was not equal to zero.
Of the 4 842 342 sources in the galaxy_candidates table, we found 534 154
matches in SDSS. 98.0% of these are classified by SDSS as GALAXY, 1.6% as QSO,
and 0.4% as STAR. Table 7 shows these percentages for each module that
contributed to the galaxy_candidates table. For example, 48 460 sources have
classlabel_dsc_joint = galaxy in the matched set, and 95.1% of these have SDSS
class GALAXY. This percentage is a measure of the purity of
classlabel_dsc_joint but only against those sources that have SDSS spectral
classifications. It is considerably higher than the purity of this DSC class
label reported in Sect. 2.1 (and detailed in Delchambre et al. 2022), even
though this purity estimate was also based on SDSS. The reason is that this
earlier purity estimate was computed for a set of sources selected at random
from Gaia: It includes the significant contamination from non-galaxies, which
outnumber galaxies by a factor of about one thousand in Gaia. The higher
figure reported in Table 7, in contrast, is just for those sources that have
SDSS spectral classifications. By design of SDSS this includes proportionally
very few stars and so very few potential contaminants. The numbers in Table 7
are a therefore a significant overestimate of the true purity and should be
treated with caution.

Figure 32: Difference between redshift_ugc and SDSS DR16 redshift, as a
function of the latter, for the 248 356 sources in common with SDSS spectral
class GALAXY.

Of the 1 367 153 sources in the galaxy_candidates table that have redshifts
provided by UGC, 248 356 match to sources in the SDSS-DR16 specObj table that
are classified as GALAXY. The small discrepancy with the number in Table 7 is
due to slightly different cross-match criteria. Figure 32 shows the difference
between redshift_ugc and the SDSS-DR16 redshift as a function of the latter.
The average of this difference is 0.06 with a standard deviation of 0.054
(which reduces to 0.029 when the 67 sources with redshifts above 0.6 are
excluded). Generally the agreement is good, although UGC seems to
systematically overestimate very low redshifts.

##  7 Composite quasar spectrum

Composite quasar spectra have many uses. First and foremost they are used as a
reference in cross-correlations of individual spectra in order to classify
these and determine their redshifts. Composite spectra are also used to
identify faint spectral features that would otherwise be undetectable, to
calibrate absolute magnitudes through the $k$-correction, as well as to
construct colour--colour relations for identifying and characterising quasars
based on photometry. Here we construct composite spectra to unveil the
capability offered by the Gaia ${\rm BP/RP}$ spectrophotometers to
characterise quasars.

Figure 33: Composite quasar spectra. The thick solid lines show composites made from the 42 944 ${\rm BP/RP}$ spectra with spectroscopically-confirmed redshifts from the Milliquas 7.2 quasar catalogue of (Flesch, 2021, i.e. type = Q in Milliquas). The different colours are for different redshift ranges. The thick dotted black line shows the composite made from 111 563 ${\rm BP/RP}$ spectra with reliable QSOC redshift estimates (identified using the query given in appendix B.2). The diagonal dotted line under each spectrum shows the quasar continuum, as described in Sect. 7 and defined in Table 9. Vertical dotted lines indicate common quasar emission lines. For comparison purposes we also show the median SDSS composite spectrum of Vanden Berk et al. (2001) (orange line). The flux densities are tabulated in Table 8.  Table 8: Composite quasar spectra shown in Figure 33. The columns are the arbitrarily scaled flux densities, $F$, and associated uncertainties, $F_{\rm err}$, computed at rest-frame wavelengths $\lambda$. The values below are for the Milliquas composite covering the redshift range $0.052$ to $4.358$. A full electronic version of this table is available with the online version of this article, along with the quasar composite spectra for the narrower redshift ranges shown in Figure 33 as well as the Gaia-only composite. $\lambda$ (nm) | $F$ | $F_{\rm err}$
---|---|---
75.669 | 0.00296 | 0.00292
81.074 | 0.00157 | 0.00133
81.318 | 0.00189 | 0.00128
81.562 | 0.00208 | 0.00124
81.807 | 0.00207 | 0.00119
… | … | …
980.799 | 0.00457 | 0.00095
983.746 | 0.00504 | 0.00142
986.701 | 0.00499 | 0.00188
989.666 | 0.00500 | 0.00222
992.639 | 0.00630 | 0.00297

In order to build composite ${\rm BP/RP}$ spectra we use the quasar sample
described in Sect. 4.2.1, which is based on the Milliquas 7.2 quasar catalogue
of Flesch (2021). Our sample comprises 42 944 sources for which we use the
Milliquas redshifts. We rely on an external catalogue instead of the Gaia DR3
internal classifications coupled with QSOC redshift determinations, because
the various instruments that contributed to the Milliquas catalogue have
better spectral resolutions and sensitivities than Gaia. We nonetheless also
compute a composite Gaia-only spectrum based on 111 563 sources coming from
the Gaia classifications and input lists together with the redshifts from
QSOC. The exact sample used for this is defined in appendix B.2. The method by
which we compute composite spectra is described in appendix C. The method
relies on a single parameter, the logarithmic wavelength sampling of the
composite spectra, chosen here to be $\log S=0.003$ (i.e. $S\approx 1.003$) as
a compromise between S/N and execution time. This sampling also applies to
transformation matrices -- $\boldsymbol{\mathsf{M_{i}}}$ in Eq. 1 \-- that
cover the observed wavelength region 309.5--1100.5 nm.

The resulting composite spectrum inferred using the Milliquas sample is shown
in Fig. 33 as the solid black line and the values listed in 8. It covers a
rest-frame wavelength range from $75.67$ nm to $992.64$ nm. The Gaia-only
composite spectrum is shown as the dotted line and covers the wavelength range
$47.08$ nm to $962.83$ nm. The composite spectra are trimmed in order to
discard wavelength regions having flux density S/N less than one. After a
multiplicative re-scaling of the flux densities so that their continua align,
we found absolute differences between these two composites -- relative to the
Ly$\alpha$ flux density in the Milliquas composite spectrum -- of less than
$4$% over the rest-frame wavelength region $100$\--$900$ nm, but up to 60% for
regions bluewards of this range. The cause of these larger deviations is
either contamination by sources with erroneous redshift estimates, as a
consequence of the low purity of QSOC in this very high redshift region (see
Delchambre et al. 2022), or border effects in the externally calibrated ${\rm
BP/RP}$ spectra that we describe below. Figure 33 also shows, for comparison,
a median composite spectrum from SDSS (Vanden Berk et al., 2001) that covers a
rest-frame wavelength range similar to that of our full composite spectrum.

Figure 34: Distribution of the literature redshift from the Milliquas 7.2
catalogue (Flesch, 2021) and Gaia $G$-band magnitudes for the 42 944 sources
used in the computation of the composite quasar spectra in Fig. 33.

Figure 34 shows the redshift and magnitude distributions of the sources used
to build the Milliquas-based composite spectrum. While the redshift
distribution is as expected, with imprints of the selection and observational
bias of each survey composing the Milliquas catalogue, the sharp drop at
$G>19$ mag is due to the filtering of the set of ${\rm BP/RP}$ spectra
published in Gaia DR3. The 17 sources with $G\geq 19$ mag are present because:
11 are associated with a best-matching node of the OA module (see Sect. 2.4),
four are used as external calibrators by CU5, and two are white dwarf
candidates used in Bellazzini at al. (2022).

Table 9: Physical quantities derived from the composite spectra of quasars from the Milliquas catalogue shown in Fig. 33. The reduced chi-square, $\chi_{\nu}^{2}$, is obtained from Eq. 1. The frequency continuum slope, $\alpha_{\nu}$, is computed from the wavelength continuum slope, $\alpha_{\lambda}=-\alpha_{\nu}-2$, where the continuum is modelled through a power law of the form $C_{\lambda}\propto\lambda^{\alpha_{\lambda}}$. In Fig. 33, the continuum is plotted as the lowest line joining the two most widely separated points in the range $121$\--$600$ nm without crossing the spectrum in this range. The emission line location, $\lambda_{\rm obs}$, and maximal line flux density, $F$, are retrieved from the fit of a quadratic polynomial in the vicinity of the laboratory wavelength, $\lambda_{\rm lab}$, after a local continuum has been subtracted that was computed in the same way as for $\alpha_{\lambda}$. The full-width at half maximum (FWHM) is retrieved from these continuum-subtracted emission lines using a linear interpolation of the spectral flux densities. All uncertainties are calculated, to first order, using the formal uncertainties on the composite spectra obtained from Eq. 4.  Composite spectrum |  Continuum slope, $\alpha_{\nu}$ |  Line+$\lambda_{\rm lab}$ (nm) |  $\lambda_{\rm obs}$ (nm) | $F/F_{\rm ref}^{a}$ | FWHM (nm)
---|---|---|---|---|---
$0.052\leq z<0.5$ | $-0.4506\pm 0.1453$ |  Mg ii$\lambda 279.875$ | $280.064\pm 0.234$ | $1.00000\pm 0.01036$ | $8.308\pm 0.214$
$2\,571$ sources |  |  H$\gamma$$\lambda 434.168$ | $435.221\pm 0.588$ | $0.13927\pm 0.00369$ | $10.432\pm 0.586$
$\chi_{\nu}^{2}=20.8$ |  |  H$\beta$$\lambda 486.268$ | $486.936\pm 0.081$ | $0.46937\pm 0.00233$ | $12.874\pm 0.131$
|  |  H$\alpha$$\lambda 656.461$ | $656.240\pm 0.010$ | $1.27766\pm 0.00095$ | $17.971\pm 0.035$
$0.5\leq z<1$ | $-0.4291\pm 0.0570$ |  C iii]$\lambda 190.873$ | $191.144\pm 0.134$ | $1.80831\pm 0.03931$ | $5.583\pm 0.262$
$8\,810$ sources |  |  Mg ii$\lambda 279.875$ | $279.644\pm 0.084$ | $1.00000\pm 0.00238$ | $10.613\pm 0.120$
$\chi_{\nu}^{2}=7.0$ |  |  H$\gamma$$\lambda 434.168$ | $433.932\pm 0.048$ | $0.25984\pm 0.00099$ | $10.355\pm 0.123$
|  |  H$\beta$$\lambda 486.268$ | $486.825\pm 0.042$ | $0.53950\pm 0.00085$ | $18.263\pm 0.082$
|  |  H$\alpha$$\lambda 656.461$ | $661.947\pm 0.100$ | $1.36875\pm 0.00735$ | $23.063\pm 0.684$
$1\leq z<2$ | $-0.6115\pm 0.0227$ |  Ly$\alpha$$\lambda 121.567$ | $122.264\pm 0.029$ | $1.00000\pm 0.00826$ | $4.411\pm 0.097$
$20\,755$ sources |  |  Si iv$\lambda 139.676$ | $139.925\pm 0.051$ | $0.14881\pm 0.00292$ | $3.637\pm 0.111$
$\chi_{\nu}^{2}=4.9$ |  |  C iv$\lambda 154.906$ | $154.636\pm 0.050$ | $0.29494\pm 0.00154$ | $4.477\pm 0.043$
|  |  C iii]$\lambda 190.873$ | $190.124\pm 0.093$ | $0.13300\pm 0.00039$ | $7.729\pm 0.068$
|  |  Mg ii$\lambda 279.875$ | $280.056\pm 0.053$ | $0.08408\pm 0.00026$ | $8.026\pm 0.047$
|  |  H$\gamma$$\lambda 434.168$ | $433.450\pm 0.101$ | $0.01752\pm 0.00013$ | $11.683\pm 0.292$
|  |  H$\beta$$\lambda 486.268$ | $490.168\pm 5.268$ | $0.03491\pm 0.00070$ | $23.559\pm 2.499$
$2\leq z<3$ | $-0.7752\pm 0.0113$ |  O vi$\lambda 103.383$ | $103.589\pm 0.152$ | $0.18338\pm 0.00686$ | $3.279\pm 0.342$
$9\,870$ sources |  |  Ly$\alpha$$\lambda 121.567$ | $122.156\pm 0.022$ | $1.00000\pm 0.00160$ | $4.908\pm 0.024$
$\chi_{\nu}^{2}=5.8$ |  |  Si iv$\lambda 139.676$ | $139.738\pm 0.086$ | $0.11467\pm 0.00068$ | $4.927\pm 0.129$
|  |  C iv$\lambda 154.906$ | $154.465\pm 0.049$ | $0.19586\pm 0.00041$ | $5.854\pm 0.114$
|  |  C iii]$\lambda 190.873$ | $190.589\pm 0.047$ | $0.11133\pm 0.00107$ | $6.742\pm 0.164$
|  |  Mg ii$\lambda 279.875$ | $279.941\pm 0.091$ | $0.07435\pm 0.00019$ | $9.471\pm 0.078$
$3\leq z<4$ | $-0.7155\pm 0.0548$ |  O vi$\lambda 103.383$ | $103.381\pm 0.112$ | $0.21926\pm 0.00277$ | $5.202\pm 0.935$
$929$ sources |  |  Ly$\alpha$$\lambda 121.567$ | $122.398\pm 0.086$ | $1.00000\pm 0.00489$ | $5.591\pm 0.083$
$\chi_{\nu}^{2}=6.0$ |  |  Si iv$\lambda 139.676$ | $139.980\pm 0.232$ | $0.12360\pm 0.00895$ | $4.588\pm 0.486$
|  |  C iv$\lambda 154.906$ | $154.581\pm 0.151$ | $0.24250\pm 0.00306$ | $3.928\pm 0.107$
|  |  C iii]$\lambda 190.873$ | $190.161\pm 0.118$ | $0.13617\pm 0.00097$ | $6.488\pm 0.094$
$0.052\leq z\leq 4.358$ | $-0.4639\pm 0.0057$ |  O vi$\lambda 103.383$ | $103.500\pm 0.103$ | $0.19285\pm 0.00584$ | $3.784\pm 0.517$
$42\,944$ sources |  |  Ly$\alpha$$\lambda 121.567$ | $122.202\pm 0.037$ | $1.00000\pm 0.00259$ | $5.088\pm 0.036$
$\chi_{\nu}^{2}=6.6$ |  |  Si iv$\lambda 139.676$ | $139.838\pm 0.141$ | $0.11467\pm 0.00146$ | $4.672\pm 0.139$
|  |  C iv$\lambda 154.906$ | $154.542\pm 0.052$ | $0.21615\pm 0.00132$ | $5.383\pm 0.066$
|  |  C iii]$\lambda 190.873$ | $190.232\pm 0.056$ | $0.11681\pm 0.00072$ | $7.938\pm 0.097$
|  |  Mg ii$\lambda 279.875$ | $280.006\pm 0.034$ | $0.08202\pm 0.00025$ | $8.891\pm 0.047$
|  |  H$\gamma$$\lambda 434.168$ | $434.164\pm 0.059$ | $0.02165\pm 0.00016$ | $10.706\pm 0.129$
|  |  H$\beta$$\lambda 486.268$ | $486.846\pm 0.109$ | $0.04867\pm 0.00011$ | $16.473\pm 0.179$
|  |  H$\alpha$$\lambda 656.461$ | $656.257\pm 0.014$ | $0.13924\pm 0.00011$ | $18.026\pm 0.032$
a If Ly$\alpha$ is covered by the composite spectrum, it is used as a
reference flux density, $F_{\rm ref}$, otherwise Mg ii is used.

In addition to the full composite spectrum computed over the whole redshift
range from $0.052$ to $4.358$, we also computed composite spectra over several
narrower redshift ranges, chosen such that their logarithmic rest-frame
wavelength coverage are approximately of equal size: $0.052\leq z<0.5$,
$0.5\leq z<1$, $1\leq z<2$, $2\leq z<3$ and $3\leq z\leq 4.358$. Composite
spectra associated with each redshift range are shown in Fig. 33. The number
of sources used in each redshift range, as well as the resulting reduced chi-
square derived from Eq. 1, are provided in Table 9 along with the frequency
($\nu$) continuum slope, $\alpha_{\nu}$ and some information on the strongest
observed emission lines: the rest-frame wavelengths, the relative flux density
at their peak compared to either Ly$\alpha$ or Mg ii, and the full width at
half maximum (FWHM). Unsurprisingly, all reduced chi-squares are larger than
unity: Each composite spectrum evidently does not completely model the
intrinsic variance seen in the observations. However, the moderate values that
are observed are indicative of reasonable fits to the observations, which is
corroborated by visual inspection and by the good agreement with the median
composite spectrum of Vanden Berk et al. (2001). The full composite spectrum
(with $\chi_{\nu}^{2}=6.6$) also explains 99.2% of the observed variance in
the entire set of spectra, although most of this variance comes from the
spectral continuum. The maximum S/N ranges from 224 per $\log S$ interval in
the composite spectrum with $3\leq z<4$, to 645 per $\log S$ interval in the
composite spectrum associated with the $0.5\leq z<1$ redshift range.

The highest S/N of 610 per $\log S$ interval of the full composite spectrum
allows us to identify many more emission lines than are otherwise visible in
the ${\rm BP/RP}$ spectra of individual sources. Whereas only the Ly$\alpha$,
C iv, C iii], Mg ii, H$\beta$, and H$\alpha$ emission lines are commonly seen
in the observed ${\rm BP/RP}$ spectra of quasars, a visual inspection of the
full composite spectrum reveals 22 emission lines, which are listed in Table
10. Despite the low resolution of ${\rm BP/RP}$ spectra (Sect. 4.2), all
common quasar emission lines were retrieved during this inspection procedure,
in addition to some weak or rarely-seen emission lines. Emission lines from
the wavelength region covering the Ly$\alpha$ forest are similarly recovered
in an unambiguous way. We achieve good agreement with laboratory wavelength
positions, with a maximum absolute difference of
$|\Delta\lambda|=|\lambda_{\rm lab}-\lambda_{\rm obs}|=0.951$ nm for the C
iii] emission line, where we consider only the nearest emission line in case
these are blended. The apparent blueshift of the C iii] emission line resides
in its asymmetry, which is due to the presence of the Si
iii]$\,\lambda\,189.203$ nm in its neighbourhood. The same rationale applies
to the shift of the Ly$\alpha$ and H$\beta$ emission lines
($\Delta\lambda=-0.683$ nm and $-0.433$ nm, respectively) due to the presence
of the N v$\,\lambda\,124.014$ nm and [O iii] doublets respectively.

Our composite spectra are highly coherent with one another, with little
variation depending on the redshift ranges that were used. After a
multiplicative re-scaling of the flux densities to align the continua, we
found differences relative to the Ly$\alpha$ flux density of less than $6$%
when compared to the full composite spectrum over several rest-frame
wavelength regions (230--900 nm for the $0.052\leq z<0.5$ composite spectrum;
180--680 nm for the $0.5\leq z<1.0$ composite spectrum; 115--500 nm for the
$1\leq z<2$ composite spectrum; 85--320 nm for the $2\leq z<3$ composite
spectrum and 85--240 nm for the $3\leq z<4$ composite spectrum). Some border
effects that were ignored in the previous comparison can still be seen in all
composite spectra. These are due to the low S/N as well as systematic errors
at the borders of the externally calibrated ${\rm BP/RP}$ spectra of quasars
that we used to build the composite spectra. Such an effect is particularly
noticeable in Fig. 33 redwards of $900$ nm, where a sharp rise is visible that
is also seen in about 90% of the individual spectra of quasars we used. This
artefact is consequently taken as a genuine signal by our method.

The differences noticed when comparing the full composite spectrum to the SDSS
one of Vanden Berk et al. (2001) redwards of the H$\beta$ emission line are
presumably due to different angular scales of the Gaia ${\rm BP/RP}$ footprint
and the SDSS fibres, which are $2.1$″(across scan direction; Carrasco et al.
2021) and $3$″respectively. This results in the contamination of the composite
spectrum by the light from the host galaxy being more suppressed in the Gaia
observations.

Table 10: Quasar emission lines found in the Milliquas-based composite spectrum and covering the redshift range $0.052\leq z\leq 4.35$ (Fig. 33). Each emission line is visually inspected before a quadratic polynomial is fit in the vicinity of its apparent peak using five samples of flux density. The maximum of the quadratic curve provides the observed rest-frame wavelength position, $\lambda_{\rm obs}$, of the line and its maximum flux density compared to Ly$\alpha$, $F/F_{{\rm Ly}\alpha}$ (where the number of significant digits is provided in parenthesis). Because of the intricacies inherent in the fit of a local continuum to faint, broad, and/or blended emission lines, such a continuum was not subtracted from the flux densities reported here. This explains the differences between this table and the values found in Table 9. Emission line | $\lambda_{\rm obs}$ | $F/F_{{\rm Ly}\alpha}$
---|---|---
+$\lambda_{\rm lab}$ (nm) | (nm) |
Ly$\epsilon$$\,\lambda\,93.780$ | $94.356\pm 1.327$ | $0.344(2)$
Ly$\delta$$\,\lambda\,94.974$ | $\cdots$ | $\cdots$
C iii$\,\lambda\,97.702$ | $98.632\pm 0.237$ | $0.393(2)$
N iii$\,\lambda\,99.069$ | $\cdots$ | $\cdots$
Ly$\beta$$\,\lambda\,102.572$ | $103.423\pm 0.139$ | $0.479(2)$
O vi$\,\lambda\,103.383$ | $\cdots$ | $\cdots$
Fe ii$\,\lambda\,111.208$a | $112.098\pm 0.500$ | $0.328(2)$
Ly$\alpha$$\,\lambda\,121.567$ | $122.250\pm 0.033$ | $1.000(2)$
O i$\,\lambda\,130.435$ | $130.495\pm 0.545$ | $0.402(2)$
Si ii$\,\lambda\,130.682$ | $\cdots$ | $\cdots$
Si iv$\,\lambda\,139.676$ | $139.568\pm 0.098$ | $0.427(2)$
O iv]$\,\lambda\,140.206$ | $\cdots$ | $\cdots$
C iv$\,\lambda\,154.906$ | $154.469\pm 0.061$ | $0.485(3)$
C iii$\,\lambda\,190.873$ | $189.922\pm 0.102$ | $0.327(3)$
Fe III$\,\lambda\,205.271$a | $205.600\pm 0.555$ | $0.230(3)$
C ii]$\,\lambda\,232.644$ | $232.390\pm 0.350$ | $0.202(3)$
[Ne iv]$\,\lambda\,242.383$ | $242.617\pm 0.324$ | $0.198(3)$
Mg ii$\,\lambda\,279.875$ | $279.888\pm 0.039$ | $0.212(3)$
He i$\,\lambda\,318.867$ | $318.612\pm 0.295$ | $0.136(3)$
H$\gamma$$\,\lambda\,434.168$ | $434.052\pm 0.067$ | $0.081(3)$
Fe ii$\,\lambda\,453.780$a | $454.281\pm 0.569$ | $0.070(4)$
H$\beta$$\,\lambda\,486.268$ | $486.701\pm 0.118$ | $0.085(4)$
[O iii]$\,\lambda\,496.030$ | $\cdots$ | $\cdots$
[O iii]$\,\lambda\,500.824$ | $\cdots$ | $\cdots$
[N i]$\,\lambda\,520.053$ | $522.270\pm 0.615$ | $0.054(4)$
He i$\,\lambda\,587.729$ | $588.441\pm 0.411$ | $0.044(4)$
H$\alpha$$\,\lambda\,656.461$ | $656.283\pm 0.014$ | $0.123(4)$
He i$\,\lambda\,706.720$ | $706.666\pm 0.945$ | $0.034(4)$
O i$\,\lambda\,844.868$ | $850.091\pm 0.413$ | $0.033(4)$
Ca ii$\,\lambda\,850.036$ | $\cdots$ | $\cdots$
a Broad feature composed of Fe multiplets.

Figure 33 and Table 9 reveal differences in the continuum slopes,
$\alpha_{\nu}$, over the various redshift ranges. These are mostly a result of
the different rest-frame wavelength coverage of each composite spectrum.
Indeed, the presence of broad Fe multiplets and the Balmer continuum in the
wavelength region $200$\--$550$ nm -- the so-called $300$ nm bump --
complicates the continuum fitting in this range due to the lack of pure
continuum. Consequently the values of $\alpha_{\nu}$ decrease once we consider
quasars with redshifts $z>1$. The value we found on the full composite
spectrum, $\alpha_{\nu}=-0.464\pm 0.005$, agrees reasonably well with
literature values. The median composite spectrum from Vanden Berk et al.
(2001) shown in Fig. 33, for example, has $\alpha_{\nu}=-0.46$.

We inspected the measurements associated with the dominant quasar emission
lines given in Table 10. We did not find any meaningful correlation between
their values and the redshift range that was used for computing each composite
spectrum.

##  8 How to select purer sub-samples

The qso_candidates and galaxy_candidates tables collate together results for
most extragalactic candidates in Gaia DR3 from a number of modules, as
described in Sect. 3.1. Overall these tables aim for high completeness, rather
than high purity. We have done this intentionally to allow the user to select
their own sub-samples using the quality indicators and class probabilities.
Here we describe how to select purer sub-samples from these tables.

Table 11: ADQL query to select the purer quasar sub-sample.

    SELECT * FROM gaiadr3.qso_candidates
    WHERE (gaia_crf_source=’true’ OR
           host_galaxy_flag<6 OR
           classlabel_dsc_joint=’quasar’ OR
           vari_best_class_name=’AGN’)

Table 12: ADQL query to select the purer galaxy sub-sample.

    SELECT * FROM gaiadr3.galaxy_candidates
    WHERE (radius_sersic IS NOT NULL OR
           classlabel_dsc_joint=’galaxy’ OR
           vari_best_class_name=’GALAXY’)

For the quasars, EO and CRF3 used input lists of quasars identified by other
surveys, so their samples are believed to be quite pure, above 90%. From the
EO sample we exclude those with close neighbours (host_galaxy_flag = 6). For
DSC, the joint subset has a purity of 62%, increasing to 79% when the Galactic
plane ($|b|<11.54\deg$) is avoided (Delchambre et al., 2022). The Vari
classifier results for AGN, which already exclude the Galactic plane, has been
assessed to have a purity of over 90% (Rimoldini et al., 2022). We therefore
recommend the query in Table 11 to select a purer sub-sample of quasars. This
returns 1 942 825 sources, which is 29% of the original table. Using the same
approach at the end of section 3.1, and assuming a 96% purity for the non-DSC
modules, we estimate the overall purity of this sub-sample to be 95%. Of
these, 1.7 million have published redshifts from QSOC.

We use similar criteria to define a purer galaxy sub-sample, except that here
there is no contribution from CRF3. Again we take all of the sources from EO
(provided by an input list), the purer subset of DSC (64% pure; 82% outside
the Galactic plane), and all of Vari. This query, in Table 12, returns 2 891
132 sources, 60% of the original table. We estimate the purity of this sub-
sample to be 94%. Of these, 1.1 million have published redshifts from UGC.

There are 14 471 sources in common between these two purer sub-samples, and
their union contains 4.8 million sources.

Figure 35: Galactic sky distribution of all the purer sub-sample of sources in
the qso_candidates table (left) and galaxy_candidates table (right). The plot
is shown at HEALpixel level 7 (0.210 sq. deg.) in Hammer--Aitoff projection.
The colour scale, which is logarithmic, covers the full range for each panel,
so is different for each panel. Compare to Fig. 5 for the full tables.

The sky distributions for these purer sub-samples are shown in Fig. 35 and can
be compared with those for the full tables in Fig. 5. We immediately see how
the purer sub-samples have lower densities in the Galactic plane. There are
still artefacts from the Gaia scanning law, which is an indication of the less
than perfect purity. We also still see overdensities at the LMC and SMC. This
comes mostly from DSC, because unlike the other modules DSC did not do any
sky-position-dependent filtering. Such sources can be easily removed by the
user, as shown in appendix B.3.

Figure 36: G-band magnitude distribution of the purer sub-sample of objects in
the qso_candidates (blue) and galaxy_candidates (orange) table on a
logarithmic scale. The dotted lines show the distributions for the full
tables.

The magnitude distributions of the purer sub-samples are shown in Fig. 36.
Compared to the distribution for the full set (dotted lines), we see that the
purer sub-sample has excluded the brightest sources (the presence of which
appears exaggerated, however, due to the logarithmic number scale). The
faintest quasars have also been removed.

Figure 37: Colour--colour diagram (top) and colour--magnitude diagram (bottom)
for the purer sub-sample of sources in the qso_candidates table (blue) and
galaxy_candidates table (orange). The contours show density on a linear scale.
The points are a random selection of 10 000 sources for each class. Compare to
Fig. 7 for the full tables.

Colour--magnitude and colour--colour diagrams of the purer sub-samples are
shown in Fig. 37 and can be compared with the same diagrams for the full
tables in Fig. 7. This shows that the purer sub-samples have a tighter colour
distribution, and remove many of the fainter galaxies.

There are other ways to obtain purer sub-samples of the integrated tables. One
could, for example, select on a higher probability threshold of the DSC
probabilities (the probabilities for all three DSC classifiers are listed in
the astrophysical_parameters table). An example is shown in appendix B.3. The
variation of purity with threshold is explored in Bailer-Jones (2021). The
joint flag used in the purer sub-samples (Tables 11 and Tables 12) corresponds
to a threshold of 0.5. One could also use a higher threshold on the
vari_best_class_score and vari_agn_membership_score rankings from the Vari-
Classification and Vari-AGN modules listed in the qso_candidates table.

In Sect. 5.4, we identified a purer sub-sample of the qso_candidates table via
an analysis of astrometric distributions that are consistent with a uniform
sky distribution of infinitely distant objects. The sources in this
astrometric selection are indicated by the boolean flag
astrometric_selection_flag in the qso_candidates table. Although this
certainly excludes genuine quasars, it should be a reasonably pure list. Of
the 1 897 754 sources in the astrometric sample, 1 801 255 are in common with
the purer quasar sub-sample defined in Table 11, and the union of these two
sets contains 2 039 324 sources. An equivalent flag is not available for the
galaxy_candidates table, for reasons discussed in Sect. 5.4.

##  9 Conclusions

We have described the data products released in Gaia DR3 for 11.3 million
candidate quasars and galaxies. This set arises from both a classification
using the Gaia data and from an analysis of sources identified by external
surveys cross-matched to Gaia. The information on these sources is presented
in the qso_candidates and galaxy_candidates integrated tables. Further
information, also on additional lower probability candidates, is provided in
several other tables (see Sect. 3). Our integrated tables are completeness
driven, so many sources in them will not be true extragalactic objects. We
therefore also provide a purer sub-sample of 4.8 million quasars and galaxies
(see Sect. 8).

We foresee a number of use cases for our results, including: aiding
confirmation of candidates found in other surveys; identifying unusual or rare
objects; providing targets for spectroscopic follow-up; providing input data
for more focused classifications or characterisations. It is our hope and
expectation that the community can build on and improve our results also by
combining them with data from other surveys.

As with previous data releases, Gaia DR3 is an intermediate data release, this
time based on 34 months of mission data. The next data release will be based
on 66 months of data, with correspondingly higher S/N and lower systematic
errors. We plan to use those data, along with improvements in our algorithms,
to update both our classifications and our characterisations of extragalactic
objects.

## Acknowledgements

This work presents results from the European Space Agency (ESA) space mission
Gaia. Gaia data are processed by the Data Processing and Analysis Consortium
(DPAC). Funding for the DPAC is provided by national institutions, in
particular the institutions participating in the Gaia MultiLateral Agreement
(MLA). The Gaia mission website is https://www.cosmos.esa.int/gaia and the
archive website is https://archives.esac.esa.int/gaia. Further acknowledgments
are given in appendix A.

## References

  * Ahumada et al. (2020) Ahumada, R., Prieto, C. A., Almeida, A., et al. 2020, ApJS, 249, 3
  * Álvarez et al. (2021) Álvarez, M. A., Dafonte, C., Manteiga, M., Garabato, D., & Santoveña, R. 2021, Neural Computing and Applications
  * Assef et al. (2018) Assef, R. J., Stern, D., Noirot, G., et al. 2018, ApJS, 234, 23
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Bailer-Jones (2021) Bailer-Jones, C. A. L. 2021, GAIA-C8-TN-MPIA-CBJ-094
  * Bailer-Jones et al. (2019) Bailer-Jones, C. A. L., Fouesneau, M., & Andrae, R. 2019, MNRAS, 490, 5615
  * Bellazzini at al. (2022) Bellazzini at al. 2022, A&A in prep.
  * Blanton et al. (2017) Blanton, M. R., Bershady, M. A., Abolfathi, B., et al. 2017, AJ, 154, 28
  * Butler & Bloom (2011) Butler, N. R. & Bloom, J. S. 2011, AJ, 141, 93
  * Carnerero et al. (2022) Carnerero et al. 2022, A&A
  * Carrasco et al. (2021) Carrasco, J. M., Weiler, M., Jordi, C., et al. 2021, A&A, 652, A86
  * Chang & Lin (2011) Chang, C.-C. & Lin, C.-J. 2011, ACM Transactions on Intelligent Systems and Technology, 2, 27:1
  * Charlot et al. (2020) Charlot, P., Jacobs, C. S., Gordon, D., et al. 2020, A&A, 644, A159
  * Cortes & Vapnik (1995) Cortes, C. & Vapnik, V. 1995, Machine Learning, 20, 273
  * Creevey et al. (2022) Creevey, O., Sordo, R., Pailler, F., & et al. 2022, A&A in prep.
  * Cutri & et al. (2012) Cutri, R. M. & et al. 2012, VizieR Online Data Catalog, II/311
  * De Angeli et al. (2022) De Angeli et al. 2022, A&A
  * de Bruijne et al. (2015) de Bruijne, J. H. J., Allen, M., Azaz, S., et al. 2015, A&A, 576, A74
  * de Souza et al. (2014) de Souza, R. E., Krone-Martins, A., dos Anjos, S., Ducourant, C., & Teixeira, R. 2014, A&A, 568, A124
  * Delchambre et al. (2019) Delchambre, L., Krone-Martins, A., Wertz, O., et al. 2019, A&A, 622, A165
  * Delchambre et al. (2022) Delchambre et al. 2022, A&A
  * Ducourant et al. (2022) Ducourant et al. 2022, A&A in prep.
  * Evans et al. (2022) Evans et al. 2022, A&A in prep.
  * Flesch (2015) Flesch, E. W. 2015, PASA, 32, e010
  * Flesch (2021) Flesch, E. W. 2021, arXiv e-prints, arXiv:2105.12985
  * Fousneau et al. (2022) Fousneau et al. 2022, A&A
  * Fu et al. (2021) Fu, Y., Wu, X.-B., Yang, Q., et al. 2021, ApJS, 254, 6
  * Gaia Collaboration (2018) Gaia Collaboration. 2018, A&A, 616, A14
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Klioner, S. A., Mignard, F., et al. 2021, A&A, 649, A9
  * Gaia Collaboration & Klioner et al. (2022) Gaia Collaboration & Klioner et al. 2022, A&A
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1
  * Heintz et al. (2015) Heintz, K. E., Fynbo, J. P. U., & Høg, E. 2015, A&A, 578, A91
  * Heintz et al. (2018) Heintz, K. E., Fynbo, J. P. U., Høg, E., et al. 2018, A&A, 615, L8
  * Holl et al. (2022) Holl et al. 2022, A&A in prep.
  * Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90
  * IAU (2021) IAU. 2021, Resolution B3, https://iau.org/administration/resolutions/general_assemblies/
  * Katz et al. (2022) Katz et al. 2022, A&A
  * Kohonen (1982) Kohonen, T. 1982, Biological Cybernetics, 43, 59
  * Krone-Martins et al. (2018) Krone-Martins, A., Delchambre, L., Wertz, O., et al. 2018, A&A, 616, L11
  * Krone-Martins et al. (2022) Krone-Martins, A., Gavras, P., Ducourant, C., et al. 2022, A&A
  * Li et al. (2021) Li, J., Silverman, J. D., Ding, X., et al. 2021, ApJ, 918, 22
  * Lindegren et al. (2021a) Lindegren, L., Bastian, U., Biermann, M., et al. 2021a, A&A, 649, A4
  * Lindegren et al. (2021b) Lindegren, L., Klioner, S. A., Hernández, J., et al. 2021b, A&A, 649, A2
  * Ma et al. (2009) Ma, C., Arias, E. F., Bianco, G., et al. 2009, IERS Technical Note, 35, 1
  * Marocco et al. (2021) Marocco, F., Eisenhardt, P. R. M., Fowler, J. W., et al. 2021, ApJS, 253, 8
  * Ochsenbein et al. (2000) Ochsenbein, F., Bauer, P., & Marcout, J. 2000, A&AS, 143, 23
  * Padilla & Strauss (2008) Padilla, N. D. & Strauss, M. A. 2008, MNRAS, 388, 1321
  * Paine et al. (2018) Paine, J., Darling, J., & Truebenbach, A. 2018, ApJS, 236, 37
  * Pâris et al. (2018) Pâris, I., Petitjean, P., Aubourg, É., et al. 2018, A&A, 613, A51
  * Pâris et al. (2017) Pâris, I., Petitjean, P., Ross, N. P., et al. 2017, A&A, 597, A79
  * Pérez & Granger (2007) Pérez, F. & Granger, B. E. 2007, Computing in Science and Engineering, 9, 21
  * R Core Team (2020) R Core Team. 2020, R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria
  * Raiteri et al. (2017) Raiteri, C. M., Villata, M., Acosta-Pulido, J. A., et al. 2017, Nature, 552, 374
  * Rimoldini et al. (2019) Rimoldini, L., Holl, B., Audard, M., et al. 2019, A&A, 625, A97
  * Rimoldini et al. (2022) Rimoldini et al. 2022, A&A
  * Secrest et al. (2015) Secrest, N. J., Dudik, R. P., Dorland, B. N., et al. 2015, ApJS, 221, 12
  * Simonetti et al. (1985) Simonetti, J. H., Cordes, J. M., & Heeschen, D. S. 1985, ApJ, 296, 46
  * Souchay et al. (2015) Souchay, J., Andrei, A. H., Barache, C., et al. 2015, A&A, 583, A75
  * Souchay et al. (2019) Souchay, J., Gattano, C., Andrei, A. H., et al. 2019, A&A, 624, A145
  * Stern et al. (2012) Stern, D., Assef, R. J., Benford, D. J., et al. 2012, ApJ, 753, 30
  * Taylor (2005) Taylor, M. B. 2005, in Astronomical Society of the Pacific Conference Series, Vol. 347, Astronomical Data Analysis Software and Systems XIV, ed. P. Shopbell, M. Britton, & R. Ebert, 29
  * Vanden Berk et al. (2001) Vanden Berk, D. E., Richards, G. T., Bauer, A., et al. 2001, AJ, 122, 549
  * Wenger et al. (2000) Wenger, M., Ochsenbein, F., Egret, D., et al. 2000, A&AS, 143, 9
  * Wu et al. (2021) Wu, Q., Liao, S., Qi, Z., et al. 2021, arXiv e-prints, arXiv:2111.02131

##  Appendix A Acknowledgements

The Gaia mission and data processing have been financially supported by, in
alphabetical order by country:

  * •

the Algerian Centre de Recherche en Astronomie, Astrophysique et Géophysique
of Bouzareah Observatory;

  * •

the Austrian Fonds zur Förderung der wissenschaftlichen Forschung (FWF) Hertha
Firnberg Programme through grants T359, P20046, and P23737;

  * •

the BELgian federal Science Policy Office (BELSPO) through various PROgramme
de Développement d’Expériences scientifiques (PRODEX) grants and the Polish
Academy of Sciences - Fonds Wetenschappelijk Onderzoek through grant
VS.091.16N, and the Fonds de la Recherche Scientifique (FNRS), and the
Research Council of Katholieke Universiteit (KU) Leuven through grant
C16/18/005 (Pushing AsteRoseismology to the next level with TESS, GaiA, and
the Sloan DIgital Sky SurvEy -- PARADISE);

  * •

the Brazil-France exchange programmes Fundação de Amparo à Pesquisa do Estado
de São Paulo (FAPESP) and Coordenação de Aperfeicoamento de Pessoal de Nível
Superior (CAPES) - Comité Français d’Evaluation de la Coopération
Universitaire et Scientifique avec le Brésil (COFECUB);

  * •

the Chilean Agencia Nacional de Investigación y Desarrollo (ANID) through
Fondo Nacional de Desarrollo Científico y Tecnológico (FONDECYT) Regular
Project 1210992 (L. Chemin);

  * •

the National Natural Science Foundation of China (NSFC) through grants
11573054, 11703065, and 12173069, the China Scholarship Council through grant
201806040200, and the Natural Science Foundation of Shanghai through grant
21ZR1474100;

  * •

the Tenure Track Pilot Programme of the Croatian Science Foundation and the
École Polytechnique Fédérale de Lausanne and the project TTP-2018-07-1171
‘Mining the Variable Sky’, with the funds of the Croatian-Swiss Research
Programme;

  * •

the Czech-Republic Ministry of Education, Youth, and Sports through grant LG
15010 and INTER-EXCELLENCE grant LTAUSA18093, and the Czech Space Office
through ESA PECS contract 98058;

  * •

the Danish Ministry of Science;

  * •

the Estonian Ministry of Education and Research through grant IUT40-1;

  * •

the European Commission’s Sixth Framework Programme through the European
Leadership in Space Astrometry (ELSA) Marie Curie Research Training Network
(MRTN-CT-2006-033481), through Marie Curie project PIOF-GA-2009-255267 (Space
AsteroSeismology & RR Lyrae stars, SAS-RRL), and through a Marie Curie
Transfer-of-Knowledge (ToK) fellowship (MTKD-CT-2004-014188); the European
Commission’s Seventh Framework Programme through grant FP7-606740
(FP7-SPACE-2013-1) for the Gaia European Network for Improved data User
Services (GENIUS) and through grant 264895 for the Gaia Research for European
Astronomy Training (GREAT-ITN) network;

  * •

the European Cooperation in Science and Technology (COST) through COST Action
CA18104 ‘Revealing the Milky Way with Gaia (MW-Gaia)’;

  * •

the European Research Council (ERC) through grants 320360, 647208, and 834148
and through the European Union’s Horizon 2020 research and innovation and
excellent science programmes through Marie Skłodowska-Curie grant 745617 (Our
Galaxy at full HD -- Gal-HD) and 895174 (The build-up and fate of self-
gravitating systems in the Universe) as well as grants 687378 (Small Bodies:
Near and Far), 682115 (Using the Magellanic Clouds to Understand the
Interaction of Galaxies), 695099 (A sub-percent distance scale from binaries
and Cepheids -- CepBin), 716155 (Structured ACCREtion Disks -- SACCRED),
951549 (Sub-percent calibration of the extragalactic distance scale in the era
of big surveys -- UniverScale), and 101004214 (Innovative Scientific Data
Exploration and Exploitation Applications for Space Sciences -- EXPLORE);

  * •

the European Science Foundation (ESF), in the framework of the Gaia Research
for European Astronomy Training Research Network Programme (GREAT-ESF);

  * •

the European Space Agency (ESA) in the framework of the Gaia project, through
the Plan for European Cooperating States (PECS) programme through contracts
C98090 and 4000106398/12/NL/KML for Hungary, through contract
4000115263/15/NL/IB for Germany, and through PROgramme de Développement
d’Expériences scientifiques (PRODEX) grant 4000127986 for Slovenia;

  * •

the Academy of Finland through grants 299543, 307157, 325805, 328654, 336546,
and 345115 and the Magnus Ehrnrooth Foundation;

  * •

the French Centre National d’Études Spatiales (CNES), the Agence Nationale de
la Recherche (ANR) through grant ANR-10-IDEX-0001-02 for the ‘Investissements
d’avenir’ programme, through grant ANR-15-CE31-0007 for project ‘Modelling the
Milky Way in the Gaia era’ (MOD4Gaia), through grant ANR-14-CE33-0014-01 for
project ‘The Milky Way disc formation in the Gaia era’ (ARCHEOGAL), through
grant ANR-15-CE31-0012-01 for project ‘Unlocking the potential of Cepheids as
primary distance calibrators’ (UnlockCepheids), through grant ANR-19-CE31-0017
for project ‘Secular evolution of galxies’ (SEGAL), and through grant
ANR-18-CE31-0006 for project ‘Galactic Dark Matter’ (GaDaMa), the Centre
National de la Recherche Scientifique (CNRS) and its SNO Gaia of the Institut
des Sciences de l’Univers (INSU), its Programmes Nationaux: Cosmologie et
Galaxies (PNCG), Gravitation Références Astronomie Métrologie (PNGRAM),
Planétologie (PNP), Physique et Chimie du Milieu Interstellaire (PCMI), and
Physique Stellaire (PNPS), the ‘Action Fédératrice Gaia’ of the Observatoire
de Paris, the Région de Franche-Comté, the Institut National Polytechnique
(INP) and the Institut National de Physique nucléaire et de Physique des
Particules (IN2P3) co-funded by CNES;

  * •

the German Aerospace Agency (Deutsches Zentrum für Luft- und Raumfahrt e.V.,
DLR) through grants 50QG0501, 50QG0601, 50QG0602, 50QG0701, 50QG0901,
50QG1001, 50QG1101, 50QG1401, 50QG1402, 50QG1403, 50QG1404, 50QG1904,
50QG2101, 50QG2102, and 50QG2202, and the Centre for Information Services and
High Performance Computing (ZIH) at the Technische Universität Dresden for
generous allocations of computer time;

  * •

the Hungarian Academy of Sciences through the Lendület Programme grants
LP2014-17 and LP2018-7 and the Hungarian National Research, Development, and
Innovation Office (NKFIH) through grant KKP-137523 (‘SeismoLab’);

  * •

the Science Foundation Ireland (SFI) through a Royal Society - SFI University
Research Fellowship (M. Fraser);

  * •

the Israel Ministry of Science and Technology through grant 3-18143 and the
Tel Aviv University Center for Artificial Intelligence and Data Science (TAD)
through a grant;

  * •

the Agenzia Spaziale Italiana (ASI) through contracts I/037/08/0, I/058/10/0,
2014-025-R.0, 2014-025-R.1.2015, and 2018-24-HH.0 to the Italian Istituto
Nazionale di Astrofisica (INAF), contract 2014-049-R.0/1/2 to INAF for the
Space Science Data Centre (SSDC, formerly known as the ASI Science Data
Center, ASDC), contracts I/008/10/0, 2013/030/I.0, 2013-030-I.0.1-2015, and
2016-17-I.0 to the Aerospace Logistics Technology Engineering Company (ALTEC
S.p.A.), INAF, and the Italian Ministry of Education, University, and Research
(Ministero dell’Istruzione, dell’Università e della Ricerca) through the
Premiale project ‘MIning The Cosmos Big Data and Innovative Italian Technology
for Frontier Astrophysics and Cosmology’ (MITiC);

  * •

the Netherlands Organisation for Scientific Research (NWO) through grant
NWO-M-614.061.414, through a VICI grant (A. Helmi), and through a Spinoza
prize (A. Helmi), and the Netherlands Research School for Astronomy (NOVA);

  * •

the Polish National Science Centre through HARMONIA grant 2018/30/M/ST9/00311
and DAINA grant 2017/27/L/ST9/03221 and the Ministry of Science and Higher
Education (MNiSW) through grant DIR/WK/2018/12;

  * •

the Portuguese Fundação para a Ciência e a Tecnologia (FCT) through national
funds, grants SFRH/BD/128840/2017 and PTDC/FIS-AST/30389/2017, and work
contract DL 57/2016/CP1364/CT0006, the Fundo Europeu de Desenvolvimento
Regional (FEDER) through grant POCI-01-0145-FEDER-030389 and its Programa
Operacional Competitividade e Internacionalização (COMPETE2020) through grants
UIDB/04434/2020 and UIDP/04434/2020, and the Strategic Programme
UIDB/00099/2020 for the Centro de Astrofísica e Gravitação (CENTRA);

  * •

the Slovenian Research Agency through grant P1-0188;

  * •

the Spanish Ministry of Economy (MINECO/FEDER, UE), the Spanish Ministry of
Science and Innovation (MICIN), the Spanish Ministry of Education, Culture,
and Sports, and the Spanish Government through grants BES-2016-078499,
BES-2017-083126, BES-C-2017-0085, ESP2016-80079-C2-1-R, ESP2016-80079-C2-2-R,
FPU16/03827, PDC2021-121059-C22, RTI2018-095076-B-C22, and TIN2015-65316-P
(‘Computación de Altas Prestaciones VII’), the Juan de la Cierva Incorporación
Programme (FJCI-2015-2671 and IJC2019-04862-I for F. Anders), the Severo Ochoa
Centre of Excellence Programme (SEV2015-0493), and
MICIN/AEI/10.13039/501100011033 (and the European Union through European
Regional Development Fund ‘A way of making Europe’) through grant
RTI2018-095076-B-C21, the Institute of Cosmos Sciences University of Barcelona
(ICCUB, Unidad de Excelencia ‘María de Maeztu’) through grant
CEX2019-000918-M, the University of Barcelona’s official doctoral programme
for the development of an R+D+i project through an Ajuts de Personal
Investigador en Formació (APIF) grant, the Spanish Virtual Observatory through
project AyA2017-84089, the Galician Regional Government, Xunta de Galicia,
through grants ED431B-2021/36, ED481A-2019/155, and ED481A-2021/296, the
Centro de Investigación en Tecnologías de la Información y las Comunicaciones
(CITIC), funded by the Xunta de Galicia and the European Union (European
Regional Development Fund -- Galicia 2014-2020 Programme), through grant
ED431G-2019/01, the Red Española de Supercomputación (RES) computer resources
at MareNostrum, the Barcelona Supercomputing Centre - Centro Nacional de
Supercomputación (BSC-CNS) through activities AECT-2017-2-0002,
AECT-2017-3-0006, AECT-2018-1-0017, AECT-2018-2-0013, AECT-2018-3-0011,
AECT-2019-1-0010, AECT-2019-2-0014, AECT-2019-3-0003, AECT-2020-1-0004, and
DATA-2020-1-0010, the Departament d’Innovació, Universitats i Empresa de la
Generalitat de Catalunya through grant 2014-SGR-1051 for project ‘Models de
Programació i Entorns d’Execució Parallels’ (MPEXPAR), and Ramon y Cajal
Fellowship RYC2018-025968-I funded by MICIN/AEI/10.13039/501100011033 and the
European Science Foundation (‘Investing in your future’);

  * •

the Swedish National Space Agency (SNSA/Rymdstyrelsen);

  * •

the Swiss State Secretariat for Education, Research, and Innovation through
the Swiss Activités Nationales Complémentaires and the Swiss National Science
Foundation through an Eccellenza Professorial Fellowship (award PCEFP2_194638
for R. Anderson);

  * •

the United Kingdom Particle Physics and Astronomy Research Council (PPARC),
the United Kingdom Science and Technology Facilities Council (STFC), and the
United Kingdom Space Agency (UKSA) through the following grants to the
University of Bristol, the University of Cambridge, the University of
Edinburgh, the University of Leicester, the Mullard Space Sciences Laboratory
of University College London, and the United Kingdom Rutherford Appleton
Laboratory (RAL): PP/D006511/1, PP/D006546/1, PP/D006570/1, ST/I000852/1,
ST/J005045/1, ST/K00056X/1, ST/K000209/1, ST/K000756/1, ST/L006561/1,
ST/N000595/1, ST/N000641/1, ST/N000978/1, ST/N001117/1, ST/S000089/1,
ST/S000976/1, ST/S000984/1, ST/S001123/1, ST/S001948/1, ST/S001980/1,
ST/S002103/1, ST/V000969/1, ST/W002469/1, ST/W002493/1, ST/W002671/1,
ST/W002809/1, and EP/V520342/1.

We made use of the following tools in the preparation of this paper: (SIMBAD,
Wenger et al. 2000) and VizieR (Ochsenbein et al. 2000) operated at (CDS)
Strasbourg; NASA ADS; TOPCAT (Taylor 2005); Matplotlib (Hunter 2007); IPython
(Pérez & Granger 2007); Astropy, a community-developed core Python package for
Astronomy (Astropy Collaboration et al. 2018); R (R Core Team 2020).

Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the
Participating Institutions, the National Science Foundation, and the U.S.
Department of Energy Office of Science. The SDSS-III web site is
http://www.sdss3.org/. SDSS-III is managed by the Astrophysical Research
Consortium for the Participating Institutions of the SDSS-III Collaboration
including the University of Arizona, the Brazilian Participation Group,
Brookhaven National Laboratory, Carnegie Mellon University, University of
Florida, the French Participation Group, the German Participation Group,
Harvard University, the Instituto de Astrofisica de Canarias, the Michigan
State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence
Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max
Planck Institute for Extraterrestrial Physics, New Mexico State University,
New York University, Ohio State University, Pennsylvania State University,
University of Portsmouth, Princeton University, the Spanish Participation
Group, University of Tokyo, University of Utah, Vanderbilt University,
University of Virginia, University of Washington, and Yale University.

Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High Performance Computing at the University of Utah. The SDSS website is www.sdss.org. SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics | Harvard & Smithsonian, the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

##  Appendix B ADQL queries

###  B.1 Purer sub-samples

The queries that provide the recommended purer sub-samples of quasars and
galaxies are shown in Tables 11 and 12. These sub-samples are discussed in
Sect. 8.

###  B.2 Sources used for the construction of a Gaia-only composite spectrum
of quasars

The following query returns the source_id of the set of sources used to
construct the Gaia-only composite spectra of quasars. It selects 111 563
sources from the purer quasar sub-sample (Table 11) that have both a reliable
redshift estimate from QSOC (flags_qsoc = 0 or flags_qsoc = 16, as explained
in Delchambre et al. 2022) and ${\rm BP/RP}$ coefficients that are published
in Gaia DR3.

    SELECT source_id
    FROM gaiadr3.qso_candidates
    JOIN gaiadr3.xp_summary USING (source_id)
    WHERE (classlabel_dsc_joint=’quasar’
           OR vari_best_class_name=’AGN’
           OR host_galaxy_detected=’true’
           OR gaia_crf_source=’true’)
      AND (flags_qsoc = 0 OR flags_qsoc = 16)

###  B.3 Using DSC probability thresholds

The following query uses thresholds on DSC-Specmod and DSC-Allosmod, which are
in astrophysical_parameters table, to select sources from the qso_candidates
table. This example uses thresholds of 0.9 and returns 371 708 sources. This
compares to 547 201 returned when using thresholds of 0.5 (which is equivalent
to, but slower than, selecting on classlabel_dsc_joint = quasar supplied in
the qso_candidates table itself). An analogous query can be used for the
galaxy_candidates table.

    SELECT source_id
    FROM gaiadr3.qso_candidates
    JOIN gaiadr3.astrophysical_parameters USING
         (source_id)
    WHERE classprob_dsc_specmod_quasar>0.9 AND
          classprob_dsc_allosmod_quasar>0.9

If we want to exclude generous regions around the LMC and SMC from the DSC
purer subset, we can use the following query (change the initial line to
select the fields you want). It removes 22 705 sources in the LMC and 7456 in
the SMC, to leave 517 040 sources.

    SELECT source_id
    FROM (SELECT *
    Ψ  FROM gaiadr3.qso_candidates
          WHERE classlabel_dsc_joint=’quasar’) AS temp
    JOIN gaiadr3.gaia_source USING (source_id)
    WHERE 1!=CONTAINS(
      POINT(’ICRS’, 81.3, -68.7),
      CIRCLE(’ICRS’, ra, dec, 9)) AND
      1!=CONTAINS(
      POINT(’ICRS’, 16.0, -72.8),
      CIRCLE(’ICRS’, ra, dec, 6))

Combining the two queries above to use higher thresholds on the DSC
probabilities and to exclude the LMC and SMC, we get the following query,
which returns 366 574 sources.

    SELECT source_id
    FROM (SELECT *
    Ψ  FROM gaiadr3.qso_candidates
    Ψ  JOIN gaiadr3.astrophysical_parameters USING
    Ψ  (source_id)
          WHERE classprob_dsc_specmod_quasar>0.9 AND
          classprob_dsc_allosmod_quasar>0.9) AS temp
    JOIN gaiadr3.gaia_source USING (source_id)
    WHERE 1!=CONTAINS(
      POINT(’ICRS’, 81.3, -68.7),
      CIRCLE(’ICRS’, ra, dec, 9)) AND
      1!=CONTAINS(
      POINT(’ICRS’, 16.0, -72.8),
      CIRCLE(’ICRS’, ra, dec, 6))

##  Appendix C Computation of the composite spectra of quasars

The computation of a composite spectrum of quasars from individual spectra may
appear to be a straightforward task, but it turns out not to be for a number
of reasons.

  1. 1.

Quasars cover a large range of redshifts, equivalently a large range of rest-
frame wavelengths, whereas each spectrum covers only a limited fraction of
these rest-frame wavelengths.

  2. 2.

Spectra have very different apparent luminosities, either because some are
intrinsically brighter or fainter, or because of their difference in redshift,
or because they are gravitationally lensed. Accordingly, each spectrum
contributing to the composite spectrum must be scaled in order to reduce the
dispersion of flux in rest-frame wavelength.

  3. 3.

Spectra often have correlated noise in their fluxes that should be taken into
account.

  4. 4.

Quasars can have different continuum slopes (or any other background signal)
that we may want to subtract in order to produce a pure emission line
composite spectrum.

  5. 5.

Spectra used to build the composite spectrum may have different resolutions,
sampling and line spread function, as in ${\rm BP/RP}$ spectra, that should be
first homogenized so as to model the sole signal of interest.

With all of these arguments in mind, we developed a new method for computing a
composite ${\rm BP/RP}$ spectrum based on maximum likelihood estimation
through the minimization of333In our notation,
$\left\|\@vec{x}\right\|_{\boldsymbol{\mathsf{W}}}^{2}=\left\|\boldsymbol{\mathsf{W}}\@vec{x}\right\|^{2}=\@vec{x}^{T}\boldsymbol{\mathsf{W}}^{T}\boldsymbol{\mathsf{W}}\@vec{x}$.

$\chi^{2}=\sum_{i=1}^{N}\;\left\|\;\@vec{x}_{i}-\boldsymbol{\mathsf{M_{i}}}\;\left[\boldsymbol{\mathsf{P}}\@vec{f}_{i}+\@vec{m}s_{i}\right]\;\right\|_{\boldsymbol{\mathsf{W_{i}}}}^{2},$
(1)

where

  * •

$N$ is the number of observations (spectra).

  * •

$\@vec{x}_{i}$ is the $i$-th observation vector, here taken as a concatenation
of the BP and RP spectral coefficients, which are coefficients associated with
a linear spline basis functions that represent the ${\rm BP/RP}$ spectra
Carrasco et al. (2021).

  * •

$\boldsymbol{\mathsf{W_{i}}}$ is the weight matrix associated with
$\@vec{x}_{i}$. If $\boldsymbol{\mathsf{L_{i}}}$ is the Cholesky decomposition
of the covariance matrix associated with $\@vec{x}_{i}$,
$\boldsymbol{\mathsf{C_{i}}}=\boldsymbol{\mathsf{L_{i}}}\boldsymbol{\mathsf{L_{i}}}^{T}$,
then $\boldsymbol{\mathsf{W_{i}}}=\boldsymbol{\mathsf{L_{i}}}^{-1}$ such that
$\boldsymbol{\mathsf{C_{i}}}^{-1}=\boldsymbol{\mathsf{W_{i}}}^{T}\boldsymbol{\mathsf{W_{i}}}$.

  * •

$\@vec{m}$ is the composite spectrum we are inferring. We additionally infer
the scaling factors, $s_{i}$, one associated with each observation $i$.

  * •

$\boldsymbol{\mathsf{P}}$ is a matrix composed of a set of basis functions
used to model a background signal to be subtracted from the spectra. The
linear coefficients associated with $\boldsymbol{\mathsf{P}}$ for the $i$-th
observation are given by the column vector $\@vec{f}_{i}$. An example use of
$\boldsymbol{\mathsf{P}}$ would be to model the quasar continua as a low order
polynomial (whose coefficients are computed in $\@vec{f}_{i}$) and to subtract
these continua from the spectra in Eq. 1. This produces a pure emission line
composite spectrum in $\@vec{m}$. The matrix $\boldsymbol{\mathsf{P}}$ could
be a set of vectors resulting from a previous minimization of Eq. 1, namely
$\boldsymbol{\mathsf{P}}=\boldsymbol{\mathsf{P_{t}}}=\begin{pmatrix}\boldsymbol{\mathsf{P_{t-1}}}&\@vec{m}_{t}\end{pmatrix}$
where $\@vec{m}_{t}$ minimizes
$\chi_{t}^{2}=\sum_{i}\;\left\|\;\@vec{x}_{i}-\boldsymbol{\mathsf{M_{i}}}\;\left[\boldsymbol{\mathsf{P_{t-1}}}\@vec{f}_{i}^{(t)}+\@vec{m}_{t}s_{i}^{(t)}\right]\;\right\|_{\boldsymbol{\mathsf{W_{i}}}}^{2}$.
This method can be seen as a weighted principal component analysis, in the
sense that $\boldsymbol{\mathsf{P_{t}}}$ are the minimal set of $t$ components
that minimize $\chi_{t}^{2}$. Although mentioned here for completeness, we
decided not to subtract the quasar continua in the present study, so we set
$\boldsymbol{\mathsf{P}}=0$ and $\@vec{f}_{i}=0$.

  * •

$\boldsymbol{\mathsf{M_{i}}}$ is a transformation matrix, associated with
$\@vec{x}_{i}$, that projects $\boldsymbol{\mathsf{P}}$ and $\@vec{m}$ into
the space of $\@vec{x}_{i}$. In the present application, the goal of
$\boldsymbol{\mathsf{M_{i}}}$ is twofold. First it isolates the rest-frame
wavelength regions from $\boldsymbol{\mathsf{P}}$ and $\@vec{m}$ that
correspond to the observed wavelength region in $\@vec{x}_{i}$ (the source
redshift must therefore be taken into account). Second, the shifted and
resampled spectrum is converted into BP/RP spectral coefficients through the
use of the GaiaXPy simulator. See the documentation on simulate_continuous for
more information on the calibration procedure.

In the minimization of Eq. 1, $\@vec{m}$ and
$\@vec{a}_{i}=\begin{pmatrix}\@vec{f}_{i}\\\ s_{i}\end{pmatrix}$ are free to
vary. However, for a given value of $\@vec{m}$, we can differentiate Eq. 1
with respect to $\@vec{a}_{i}$ and set the resulting gradient to zero to
get444Direct inversion of the normal equations is known to be numerically
unstable and should be avoided.

$\@vec{a}_{i}=\left(\boldsymbol{\mathsf{T_{i}}}^{T}\boldsymbol{\mathsf{W_{i}}}^{T}\boldsymbol{\mathsf{W_{i}}}\boldsymbol{\mathsf{T_{i}}}\right)^{-1}\boldsymbol{\mathsf{T_{i}}}^{T}\boldsymbol{\mathsf{W_{i}}}^{T}\boldsymbol{\mathsf{W_{i}}}\@vec{x}_{i}\hskip
14.22636pt\text{where}\hskip
<EMAIL_ADDRESS>(2)

Substituting this last equation into Eq. 1 makes it depend only on $\@vec{m}$,
such that any (global) optimization algorithm can be used with a number of
unknowns given by the number of fluxes in $\@vec{m}$. In the present study, we
use an expectation-maximization algorithm with momentum in batch mode. The
steps of the expectation-maximization algorithm are: (i) compute
$\@vec{a}_{i}$ using Eq. 2; (ii) fit $\@vec{m}$ to all $\@vec{x}_{i}$ given
the previously computed $\@vec{a}_{i}$,

$\@vec{m}=\left(\sum_{i=1}^{N}s_{i}^{2}\;\boldsymbol{\mathsf{M_{i}}}^{T}\boldsymbol{\mathsf{W_{i}}}^{T}\boldsymbol{\mathsf{W_{i}}}\boldsymbol{\mathsf{M_{i}}}\right)^{-1}\left(\sum_{i=1}^{N}s_{i}\;\boldsymbol{\mathsf{M_{i}}}^{T}\boldsymbol{\mathsf{W_{i}}}^{T}\boldsymbol{\mathsf{W_{i}}}\left[\@vec{x}_{i}-\boldsymbol{\mathsf{M_{i}}}\boldsymbol{\mathsf{P}}\@vec{f}_{i}\right]\right).$
(3)

Steps (i) and (ii) are then iterated until the reduced chi-square improves by
no more that 0.001 for 16 consecutive iterations.

To first order, the covariance matrix associated with the formal uncertainties
of the computed composite spectrum can be approximated through the asymptotic
normality property of the maximum likelihood estimator as

$\boldsymbol{\mathsf{\Sigma_{m}}}\approx\left(\;\sum_{i=1}^{N}\boldsymbol{\mathsf{J_{i}}}^{T}\boldsymbol{\mathsf{W_{i}}}^{T}\boldsymbol{\mathsf{W_{i}}}\boldsymbol{\mathsf{J_{i}}}\;\right)^{-1}$
(4)

where

$\boldsymbol{\mathsf{J_{i}}}=\frac{\partial\boldsymbol{\mathsf{T_{i}}}\@vec{a}_{i}}{\partial\@vec{m}}=\boldsymbol{\mathsf{T_{i}}}\frac{\partial\@vec{a}_{i}}{\partial\@vec{m}}+\boldsymbol{\mathsf{M_{i}}}\;s_{i}$

and

$\frac{\partial\@vec{a}_{i}}{\partial\@vec{m}}=\left(\boldsymbol{\mathsf{T_{i}}}^{T}\boldsymbol{\mathsf{W_{i}}}^{T}\boldsymbol{\mathsf{W_{i}}}\boldsymbol{\mathsf{T_{i}}}\right)^{-1}\begin{pmatrix}-s_{i}\left[\boldsymbol{\mathsf{M_{i}}}\boldsymbol{\mathsf{P}}\right]^{T}\\\
\left[\@vec{x}_{i}-\boldsymbol{\mathsf{M_{i}}}\boldsymbol{\mathsf{P}}\<EMAIL_ADDRESS>
Equation 4 is only valid for large values of $N$. How large $N$ should be is
problem dependent. We performed simulations using 65 536 noisy realisations of
a problem with $N=64$ and eight variables in $\@vec{m}$, where all matrices,
except $\boldsymbol{\mathsf{W_{i}}}$, are uniformly distributed in $[-1,1]$
and $\boldsymbol{\mathsf{W_{i}}}$ is an orthogonal transformation of matrices
whose eigenvalues are uniformly drawn in $[0.001,1]$. This led to a maximum
absolute error in the median correlation coefficients of $0.007$.

The Octave/Matlab source code for minimizing Equation 1 and for computing the
approximate covariance matrix from Equation 4 is available at
https://github.com/ldelchambre/gls_mean/.