The probabilistic random forest applied to the QUBRICS survey: improving the selection of high-redshift quasars with synthetic data

ABSTRACT

Several recent works have focused on the search for bright, high-z quasars (QSOs) in the South. Among them, the QUasars as BRIght beacons for Cosmology in the Southern hemisphere (QUBRICS) survey has now delivered hundreds of new spectroscopically confirmed QSOs selected by means of machine learning algorithms. Building upon the results obtained by introducing the probabilistic random forest (PRF) for the QUBRICS selection, we explore in this work the feasibility of training the algorithm on synthetic data to improve the completeness in the higher redshift bins. We also compare the performances of the algorithm if colours are used as primary features instead of magnitudes. We generate synthetic data based on a composite QSO spectral energy distribution. We first train the PRF to identify QSOs among stars and galaxies, then separate high-z quasar from low-z contaminants. We apply the algorithm on an updated data set, based on SkyMapper DR3, combined with Gaia eDR3, 2MASS, and WISE magnitudes. We find that employing colours as features slightly improves the results with respect to the algorithm trained on magnitude data. Adding synthetic data to the training set provides significantly better results with respect to the PRF trained only on spectroscopically confirmed QSOs. We estimate, on a testing data set, a completeness of |$\sim 86{{\ \rm per\ cent}}$| and a contamination of |$\sim 36{{\ \rm per\ cent}}$|⁠. Finally, 206 PRF-selected candidates were observed: 149/206 turned out to be genuine QSOs with z > 2.5, 41 with z < 2.5, 3 galaxies and 13 stars. The result confirms the ability of the PRF to select high-z quasars in large data sets.

1 INTRODUCTION

Luminous, high-redshift quasars are of paramount importance for a wide range of extragalactic studies. These include, e.g. the number density of bright quasars at high-z (Schindler et al. 2019a; Boutsia et al. 2021; Onken et al. 2021; Grazian et al. 2022), their role in the re-ionization process (e.g. Fontanot, Cristiani & Vanzella 2012, Fontanot et al., in preparation), the early phases of galaxy formation and co-evolution with their central SMBHs (e.g. Fontanot et al. 2020), the inference of cosmological parameters from lensed QSOs or the characterisation of the accretion properties for SMBHs (e.g. Wu et al. 2015).

High-z, bright QSOs allow to investigate the properties of the intergalactic medium, which at z > 1.5 contains more than 80 per cent of all baryonic matter: light from these beacons carries information on the spatial distribution, chemical enrichment, and kinetic properties of the interposed gas. Example of relevant studies include: the determination of primordial element abundances (e.g. deuterium, Cooke et al. 2014; Cooke, Pettini & Steidel 2017), the variation of fundamental constants (e.g. the fine-structure constant; Milaković et al. 2021; Murphy et al. 2022) or the ability to directly probe the expansion history of the Universe by measuring the redshift drift (Liske et al. 2008).

The best candidates for these studies (i.e. bright, high-redshift QSOs) are, however, difficult to select in wide surveys. This is evident when comparing their surface densities with that of other objects at the same apparent magnitude: for instance, less than 1 per cent of the objects in the sample used in Boutsia et al. (2020) were quasars with i_psf < 18 and z > 2.5. Historically, the largest efforts in searching for QSO targets have been made in the Northern hemisphere: the SDSS, which has now delivered more than 7.5 × 10⁵ spectroscopically confirmed QSOs between redshift 0 and 7 (Lyke et al. 2020), is the most prominent example of these endeavours. Inferior investment in telescope time in the Southern hemisphere has led to a less favourable situation: undiscovered targets are waiting to be found, by exploiting recent photometric catalogues targeting the southern sky (SkyMapper, DESI, and PanSTARRS; Chambers et al. 2016; Onken et al. 2019; Sevilla-Noarbe et al. 2021).

Taking advantage of these large catalogues, the QUBRICS program started to search for bright and high-z QSOs by means of machine learning techniques applied on photometric data. The first catalogues of candidates have been derived using a canonical correlation analysis (CCA; Anderson 2003; Calderone et al. 2019) trained on known, spectroscopically confirmed QSOs from literature. We applied the CCA both as classifier (i.e. an algorithm to label individual sources as star, non-active galaxy, or QSO) and as regressor (i.e. as an estimator for the photometric redshift for each potential QSO). We obtained spectroscopical observations of the most promising quasar candidates with telescopes and confirmed the QSO nature of more than 250 sources (Fig. 1, Calderone et al. 2019; Boutsia et al. 2020). In Guarneri et al. (2021, hereafter Paper I), we extended the range of available ML techniques by considering a new selection method based on the PRF.

The main advantages of this second method are the non-linearity (and thus the ability to produce an improved predictive model), the ability to incorporate measurement uncertainties in the model, and to handle missing data. This paper builds upon the results presented in Paper I and reports our attempts in improving the sample selection of candidates (in particular at z > 3.5).

The paper is organized as follows: in Section 2, we describe the base catalogue and the selection algorithm; in Section 3, we review the approach followed in Paper I and possible improvements. In particular, in order to obtain a more effective training set we introduce synthetic QSOs. Section 4 describes the method for generating synthetic data; Sections 5 and 6 concerns the analysis of the performance of the algorithm and the characterization of the candidate sample; Section 7 reports the spectroscopic observations. Conclusions are finally drawn in Section 9.

Unless stated otherwise, magnitudes are given in the AB magnitude system; uncertainties represent 68 per cent confidence intervals. We adopt a flat ΛCDM cosmology, with |$\Omega _\mathrm{m}=0.3,\ \Omega _\Lambda = 0.7$|⁠, and H₀ = 67.7 kms⁻¹ Mpc⁻¹.

2 THE SELECTION METHODS AND THE CATALOGUE

2.1 The selection algorithm

Machine learning methods are becoming extremely common in astrophysics (e.g. Baron 2019) to mine the growing wealth of information produced by modern astronomical facilities. The algorithm used for this work is the Probabilistic Random Forest (Reis, Baron & Shahaf 2019), a generalization of the original random forest (RF; Breiman 2001) developed to account for measurements uncertainties. The PRF is applied as a supervised, classification algorithm, that maps features into discrete labels; in the current state, regression cannot be performed with the PRF. In machine learning terms, features are properties associated with a given object in a data set: in our case, features are magnitude estimates or colours (i.e. difference of magnitudes in two bands). Labels, instead, describe the class to which each object belongs to: in our case, we use star, non-active galaxies, low-z QSO, or high-z QSO (with the latter two classes being separated at z = 2.5).

We refer the interested reader to the original PRF paper (Reis et al. 2019) for a detailed description of the algorithm and the changes with respect to a classic random forest; we will provide below a very brief overview of the PRF algorithm.

Inheriting from the original random forest, the PRF is an ensemble algorithm that operates by creating a large number of decision trees. Each decision tree is a model described by a tree-like graph, built as a sequence of consecutive nodes. Intermediate nodes carry a condition, while the terminal node (i.e. a leaf of the tree) identifies a classification label; conditions and labels are determined during the training stage. During the prediction stage, objects are propagated along each decision tree and classified according to the leaves they reach. Each tree provides an independent classification: a majority vote among all trees provides the final class returned by the forest. In addition to a class label, both the PRF and the RF output the probability, for a given object, of belonging to a given class. Both the random forest and the PRF map features to labels: the PRF, however, improves this mapping by taking into account measured uncertainties as well. This is accomplished by treating each input feature as a probability distribution function, with variance equal to the provided squared error and expectation value equal to the given feature. This change provides the PRF advantages over the original random forest: the model is more robust, with improved generalization ability, while naturally handling missing data. This is achieved by propagating missing data to both sides of a node with equal probability. This last peculiarity is clearly desirable when working on photometric magnitudes, since missing data and non-detections are rather common.

Following the same two-step approach described in Paper I, we adopted the Probabilistic Random Forest in order to select the most promising high-redshift QSO candidate.

In the first step, we classify unlabelled sources as either QSOs, stars, or galaxies based on their photometric properties. In the second, objects predicted to be QSOs are then reclassified as high- (z > 2.5) or low- (z ≤ 2.5) redshift candidates. In both cases, we first train the PRF: as for the first step, we employ stars, galaxies and QSOs at all redshifts; we then use QSOs only as the training set for the second step. Once the algorithm is trained, it is applied on a validation/testing data set to evaluate its performances. Instead of using all stars, we undersample the available ones to uniformly cover the available i − z colour space and match the number of available QSOs and galaxies, reducing their number from ∼800 000 to ∼7500. For the second step, the algorithm is only trained on QSO data. Unfortunately, the number of high-redshift QSOs is ∼1/8 compared to the number of their low-z counterpart. To account for this, in Paper I high-redshift QSOs were oversampled to provide a more uniform redshift distribution; in this work instead synthetic data are introduced as an alternative to oversampling.

2.2 The new main sample

The multiwavelength photometric data base covers the observer’s frame from the UV to the far-IR and was derived from the following surveys:

• u_psf, v_psf, g_psf, r_psf, i_psf, z_psf magnitudes from the SkyMapper DR3 survey (Onken et al. 2019)

• G, G_BP, G_RP magnitudes from the Gaia eDR3 catalogue (Gaia Collaboration 2021)

• J, H, K from 2MASS (Skrutskie et al. 2006)

• W1, W2, W3, W4 from the AllWise (Wright et al. 2010) survey.

The choice of point spread function magnitudes is motivated by the point-like appearance of QSOs.

The main sample, hereafter MS, contains all the sources satisfying the following conditions:

• i_psf magnitudes with photometric flag i_flags = 0 (14 ≤ i_psf < 18). This is equivalent to require that the i_psf magnitude is free from warnings concerning saturation, close neighbours, edge-of-CCD effects, and other systematic that might affect the photometry¹;

• z_psf magnitudes with photometric flag z_flags = 0;

• Non nullG magnitudes from Gaia eDR3. Moreover, the associated source in the Gaia catalogue must be within 0.5 arcsec of the SkyMapper detection;

• Non null, SNR > 3 AllWISE magnitudes in the W1, W2, and W3 bands. Moreover, the associated source in the AllWise catalogue must be within 0.5 arcsec of the SkyMapper detection;

We restrict our selection to 14 ≤ i_psf < 18, excluding crowded regions with galactic latitude |b_gal| < 25 deg: this includes 1147 017 sources, spread on |$\sim 14,500 \rm ~deg^2$| mostly in the Southern hemisphere. The G_BP, G_RP (Gaia eDR3), J, H, K (2MASS), and W4 (AllWISE) magnitudes were added to the data base when available, but their existence was not a necessary condition for an object to be included in the MS. It should be noted that the sample considered for this work is rather similar to the Main Sample employed in Paper I. The most significant differences are the exclusion of near and far-UV magnitudes from GALEX (Bianchi & Shiao 2020) due to modelling uncertainties when generating synthetic data and a slightly larger area covered by the survey (Fig. 3), thanks to the third SkyMapper data release that includes regions around 0 < Dec. ≲ 15 deg.

Objects with parallax and proper motion significantly different from zero (>3σ) were classified as stars (79 per cent of the MS) while those with high likelihood of being extended were identified by comparing Petrosian and PSF magnitudes (17 per cent of the MS, see also Calderone et al. 2019, for more details). The remaining 4 per cent of sources were cross-matched against recent spectroscopic catalogues, e.g. the SDSSDR16Q (Lyke et al. 2020), the Veron-Cétty catalogue (Véron-Cetty & Véron 2010), the 2dF (Colless et al. 2001), the 6dF (Jones et al. 2009), in order to identify previously known QSOs and non-active galaxies. Additional spectroscopic identifications, when available, were drawn from the literature (e.g. Schindler et al. 2019a,b; Wolf et al. 2020; Onken et al. 2021) leaving 31 550 unclassified sources: these will automatically be classified by the PRF trained on the spectroscopically confirmed QSOs and galaxies, together with bonafide stars.

As described in Boutsia et al. (2020), the selection criteria adopted in this work might be biased against lensed sources, as positional matches were carried out with relatively stringent criteria to reduce contamination. This should not affect point-like sources, such as non-lensed QSOs, but might lead to the removal of extended objects (e.g. lensed sources with Einstein rings or crosses); moreover, the choice of some photometric bands (e.g. Gaia G) might remove very high redshift objects (z ≳ 5) from the initial sample. We decided against more relaxed constraints in order to maximize the success rate of telescope runs and maintain high success rate.

3 PERFORMANCE OF THE PRF IN PAPER I

In Paper I we carried out a preliminary analysis on the performance of the PRF algorithm applied to the selection of high-redshift (z > 2.5) QSOs. Extensive tests guided the choice of an approach aimed at maximizing the success rate of the observations while still maintaining a sufficiently high completeness for applications such as the estimate of luminosity functions. The algorithm showed remarkably good performance in separating QSOs candidates from stars or non-active galaxies (especially for z > 2.5 QSOs, as |$\sim 98{{\ \rm per\ cent}}$| of them were correctly classified as such). In a second stage, high-redshift QSOs were sieved from lower redshift QSOs: the estimated completeness, on a testing data set, was |$\sim 84{{\ \rm per\ cent}}$|⁠, with a relatively low contamination of |$\sim 22{{\ \rm per\ cent}}$| of z < 2.5 QSOs.

Spectroscopic follow-up and further tests allowed us to better characterize the properties of this PRF selection. A small number of candidates (16) turned out to be spectroscopically confirmed non-active galaxies or stars. In paper I we estimated their combined contribution to account for ∼25 of the 626 candidates found in the sample (4 per cent of the total); the current number of observed contaminants, with roughly half of the sample missing a spectroscopic confirmation, corroborates our estimate. However, further tests have shown that few z > 3 QSOs were not successfully classified by the PRF. Misclassified sources are either among the higher redshift targets (z ≳ 4.6) or the brighter, high-redshift QSOs, i.e. the incompleteness of the selection in Paper I affects the most precious QSOs for our scientific goals.

We attempt to interpret the latter shortcoming as a consequence of non-optimal sampling of high-redshift QSOs in the training set: a larger, controlled data set would improve the generalization capabilities of the algorithm and is expected to recover at least part of the misclassifications. As mentioned in the introduction, the number of actual, spectroscopically confirmed, QSOs available to the scientific community in the Southern hemisphere is still small, hence the training sample cannot be improved using the data from the literature. Besides, failures (such as those listed in Section 5.3) of spectroscopic pipelines (e.g. SDSS) might affect the quality of the training data set.

Synthetic data offer a solution to both problems: large samples can be generated rather easily, with complete control over the data set features, and fed to the selection algorithm to improve its performance.

4 SYNTHETIC DATA

Synthetic data are needed in particular at the highest redshifts and at brighter magnitudes, where only a few QSOs are present in the training set. To produce synthetic spectral energy distributions (SEDs) we rely on a relatively ‘standard’ approach, parametrizing a composite QSO SED and its expected variations, and extracting from it realizations of the QUBRICS photometry for objects at the desired redshifts.

4.1 The quasars’ composite SED

The first step for the generation of synthetic data consists in the parametrization of the spectral energy distribution (SED) of quasars. We compiled a data set with 6074 QSOs in the QUBRICS data base with Gaia G < 18.5 and spectroscopic redshift z > 0.5. For each object, when available, the photometry in the G_BP and G_RB bands, SkyMapper, 2MASS, WISE, SDSS-DR16, PanSTARRS1-DR2, VHS, GALEX was included. For each of the 6074 QSOs we have parametrized its photometric data points in terms of:

• an SED (shown in Fig. 4 with a continuos blue line), adapted from the work by Vanden Berk et al. (2001), Richards et al. (2006), and Krawczyk et al. (2013). We have empirically decomposed it in a ‘blue bump’ (cyan dotted line), representing the accretion disc (Laor & Netzer 1989; Sun & Malkan 1989), and an ‘IR bump’ (purple dotted line), representing a dusty torus (Pier & Krolik 1993; Mor, Netzer & Elitzur 2009), with a cross-over at 10⁴ Å. Different SEDs can be generated by chosing a different fractional contribution, f, of the ‘IR bump’ with respect to the ‘blue bump’;

• a multiplicative power law, PL ∝ (λ/2200 Å)^−α, applied to the ‘blue bump’ part of the SED to account for extinction in the rest-frame optical–UV and possibly different spectral slopes in the blue bump;

• the spectroscopic redshift, used to compute the IGM absorption shortwards of 1215 Å, according to Inoue et al. (2014), and shift the SED to the observer’s frame;

We obtain synthetic photometry by convolving the generated QSO SEDs, as would be measured in the observer’s frame, with the transmission curves of the filters adopted in the input photometric surveys.For each of the 6074 known QSOs, given the spectroscopic redshift, we have fitted the fractional contribution, f, and the slope, α, optimizing, chi-square-wise, the synthetic photometry with respect to the observed colours. An example of SED fit for a QSO with spectroscopic redshift z = 3.504 is shown in Fig. 5.

As a result, we obtain the distributions of f, whose median turns out to be <f > = 0.57, with a median absolute deviation (MAD) of 0.05 and of the slopes, α, with <α > = − 0.01 and MAD of 0.18. In the next section these distributions will be used to generate synthetic photometric data of quasars at different redshifts, aiming at a fair representation of the QSO population.

4.2 Generation of synthetic photometry

In order to produce synthetic QSOs to be added to our training set we proceeded as follows:

• realizations of QSO SEDs were produced on the basis of the procedure described in the previous section, in particular extracting random values for the fractional contribution, f, of the IR bump with respect to the blue bump, and for the slope α from the observed distributions;

• the SEDs were redshifted in the observer’s frame, and multiplied by the IGM transmission;

• the SEDs were multiplied with the surveys band-passes to obtain instrumental magnitudes;

• the instrumental magnitudes were re-scaled in order to reproduce the distribution of the real sample. To this end, we considered again the sample of spectroscopically confirmed QSOs. We binned the available quasars in redshift: bins were chosen to ensure a minimum number of objects per bin (∼150) and, when possible, a bin-width of ∼0.2. For each redshift bin, we considered the i_psf magnitude distribution and extracted a random number (using a Monte Carlo method) from the same distribution. We used the difference between the extracted number and the synthetic i magnitude as offset for the magnitudes in all other bands;

• we associated to each synthetic magnitude an error estimate. In order to do so, available photometric data were used to estimate an empirical error function f_e(m_b). For each photometric band, we built a 2D histogram with magnitude on the x-axis and photometric error on the y-axis using all photometric data available in the MS (Fig. 6); the number of bins on both axis was chosen arbitrarily (250). For each magnitude bin (or, equivalently, each column in Fig. 6), we searched for the most populated error bin and interpolated, using a smoothing spline, each error point as a function of magnitude. This spline is the function f_e(m_b), which returns the typical uncertainty assigned to the magnitude m_b. The scatter in the uncertainties as a function of magnitude, s_e(m_b), was estimated as well using the same procedure. Rather than considering the most populated bin, we computed, as a first approximation for the scatter, the standard deviation across a column, and interpolated these points with a smoothing spline as a function of magnitude.

  Given f_e(m_b) and s_e(m_b), the actual photometric uncertainty |$e_{m_i}$| assigned to the magnitude m_i was estimated by drawing a random value from a Gaussian with:

  $$\begin{eqnarray} \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\mu = 0.1, \ \sigma = 0.01 & f_e(m_i) \le 0.1 \\ \mu = f_e(m_i), \ \sigma = s_e(m_i) & f_e(m_i) \gt 0.1 \end{array}\right. \end{eqnarray}$$

  The process was repeated for each band in the catalogue;

• finally, we recomputed a noisy photometric magnitude. A new value was drawn from a Gaussian distribution with μ equal to the noiseless magnitude, derived from the synthetic SED, and |$\sigma = e_{m_i}$|⁠. This is repeated for all magnitudes of a given object, and the newly found noisy magnitude, with the associated error, are added to the catalogue of synthetic objects.

In the catalogue of synthetic objects all sources have all the magnitude measurements and no missing data are present. Real photometric catalogues, however, often have missing data or non-detections, which are kept into account in the PRF selection, as described in Paper I and, briefly, in Section 2.

Once upper limits are introduced, the catalogue is ready to be fed to the PRF. The optimal number of synthetic data to employ and their redshift distribution will be discussed in Section 5. An example of SED of a synthetic object is shown in Fig. 7; the upper panel shows the SED of a spectroscopically confirmed QSO of comparable redshift and i_psf magnitude.

5 PERFORMANCE WITH SYNTHETIC DATA

In order to characterize the effect of synthetic data on the selection several tests have been performed. The selection of high-redshift QSOs proceeds in two steps: in the first the PRF is used to differentiate QSOs (at all redshifts) from stars, galaxies, and other contaminants. In the second task, instead, we aim to exclude z < 2.5 targets to obtain a sample of high-redshift QSOs ready for spectroscopic follow-up.

Unless stated otherwise, the same setup was considered: (i) the available, spectroscopically confirmed high-z QSOs (hereafter also referred to as ‘real’ QSOs) were split in training (75 per cent) and testing data set (25 per cent of the total); the same fractions (75/25) were used for low-z QSOs; additionally, 100 different realizations were drawn for each testing in order to not depend on the particular training-testing choice; (ii) synthetic QSOs were added only to the training data set, in addition to the spectroscopically confirmed QSOs. The optimal ratio synthetic/spectroscopically confirmed QSOs is not known a priori: several values were tested in order to identify the configuration that provides the best performance. This sample was used to understand the capabilities of the PRF in classifying QSOs as low- or high-z candidates.

If training the PRF to classify targets as star, galaxy or QSO at any redshift, stars and galaxies need to be added to the training and testing data set on top of the QSOs. Bonafide stars were down-sampled based on their i − z colours, i.e. a subset was randomly chosen among all available stars. The number of stars to consider is set by the number of QSOs (∼6000), in order to maintain a balanced data set [see e.g. Japkowicz & Stephen (2002) for a review on the class unbalance problem]. All spectroscopically confirmed galaxies are instead considered. Stars and galaxies are split according to the 75/25 per cent rule, and added to the appropriate QSO data set – training or testing.

Performance is tested solely on spectroscopically confirmed QSOs (and, when appropriate, known stars and galaxies – see Section 5.2.1): synthetic data are never considered in testing data sets when reporting performance values.

The algorithm was evaluated in terms of precision and recall (Fig. 8), defined as

• precision: the fraction of relevant instances among all those classified as high-z QSOs. The precision describes how good the algorithm is at excluding contaminants: in the context of this work it describes how good the PRF is at classifying stars, galaxies and low-redshift QSOs as such and not confuse these object as high-redshift QSOs;

• recall: the fraction of relevant instance (i.e. real high-z QSOs) that were correctly classified by the algorithm. The recall describes how good the algorithm is at retrieving relevant instances. In the context of this work, it describes how good the PRF is at identifying all high-redshift QSOs present in the MS.

Objects are classified as low- or high- redshift based on a threshold set at z = 2.5. However, only when evaluating the recall, two different redshift thresholds, 2.5 and 3.0, are considered. The inclusion of the second higher redshift threshold is justified since targets at z > 3.0 are the most interesting for our purposes. Moreover, we expect that most of the high-z misclassification to be in the 2.5 < z < 3.0 range, as they are closer to the classification threshold. A ‘global’ recall score would fail to highlight the performances on the most interesting objects (those with z > 3).

5.1 Colour selection

The selection of high-redshift quasar is often performed by employing colour cuts: optical and infrared colour criteria have been successfully employed in the past for the selection of high-redshift quasars (e.g. Wolf et al. 2020; Onken et al. 2021), and machine learning methods have been successfully employing colours as features (e.g. Nakoneczny et al. 2021; Wenzl et al. 2021). Low- and high- redshift QSOs, non active galaxies and stars occupy (mostly) distinct regions in a multidimensional colour space,² allowing an algorithm which employs colours to successfully differentiate between classes. Moreover, we expect the colour distributions to be independent from luminosity, thus allowing to effectively retrieve bright targets that may be under-represented in a training sample. With respect to Paper I, where only magnitudes were employed as features, in this work we tested colours as well.

5.1.1 Colour determination

Given a complete set of magnitudes, colours are trivial to generate, but, due to the introduction of photometric upper-limits, some care is required. In this work all colours are computed as difference between the b_th band and the SkyMapper i band, used as a reference. This choice is natural, given that the i band is required for all sources in the main sample. We distinguish three cases:

• if the magnitude in the b_th band is not an upper limit, the colour is simply computed as the difference between the measured magnitude and the corresponding i band magnitude. The associated error is the square root of the sum of the photometric errors squared;

• if the magnitude in the b_th band is instead an upper-limit, we first compute a new m_{t, n} factor (equation 3), specific for each source and photometric band, as the difference between the previously determined value of m_{t, n} and the corresponding i magnitude. The associated error is |$\sqrt{\sigma _i^2+\sigma _{t,n}^2}$|⁠, where σ_i is the error associated with the i magnitude.

• if, finally, the magnitude in the b_th band is a proper missing data (i.e. it was not observed by a given survey), the colour is marked as a proper missing data as well.

Following these three prescriptions, colours are computed both for synthetic and spectroscopically confirmed QSOs.

5.2 Performance evaluation

5.2.1 QSO, stars, and galaxies

The first task for the selection consists in separating QSOs at any redshift from stars and galaxies. A pre-selection against stars is performed by applying cuts on measured proper motion and parallaxes (i.e. objects with parallax and proper motion significantly different from zero (>3σ) are considered stars), and against non-active low-z galaxies by excluding extended sources (Section 2.2).

The same procedure applied in paper I, and briefly summarized in the previous sections, was repeated. The training set is composed of stars, galaxies, low-z and high-z real QSOs and a varying number of synthetic QSOs in roughly equal parts. Stars were randomly under-sampled in order to uniformly cover the entire i − z colour space, while all available QSOs and galaxies were considered. We test the performance of the trained algorithm by applying it on a testing data set (stars, non-active galaxies and QSOs that were not used during the training) and evaluating the recall and the precision.

The effect of synthetic data was tested, finding no significant improvements in recall and precision over the approach adopted in Paper I. Synthetic QSOs are introduced with z > z_th and in varying number. Several redshift threshold z_th were tried, from z_th = 1.5 to z_th = 3.0 with Δz = 0.25, as well as different numbers of QSOs, ranging from N_synth = 3000 to N_synth = 21 000, with ΔN_synth = 2000. In all cases, using precision and recall as metrics, the performance measured on the testing data set are compatible with one another, and with those reported in Paper I. Although non relevant here, synthetic data will become precious to distinguish low-z from high-z QSOs (see Section 5.2.2). To avoid biases related to a particular testing set, 100 runs were performed, each with a different, random, training/testing split. On average, the PRF achieves a recall of 95 per cent for QSOs at all redshifts (99 per cent for QSOs at z > 2.5, missing 0 or 1 QSO per test run). We measure a precision of 96 per cent: most of the contaminants are non-active galaxies, accounting for |$\sim 75{{\ \rm per\ cent}}$| of all misclassifications.

The reported results do not change if training the PRF on colours or on magnitudes data: in this context the PRF performs equally well.

5.2.2 Distinguishing low from high-redshift QSOs

The second step of the selection aims at selecting high-z QSOs: we once again introduce synthetic data in the training data set in order to improve the generalization abilities of the PRF at high redshift. As reported in Section 5.2.1, the performance of the PRF in distinguishing among QSOs, stars and galaxies was comparable either if the algorithm was trained only on spectroscopically confirmed QSOs or on a data set also including synthetic data. It is thus worth investigating if training the algorithm using the combined data set provides improvements in terms of predictive capabilities of the PRF for classifying low- and high-redshift QSOs. To address this question, we considered a training data set composed of both spectroscopically confirmed QSOs and synthetic data in equal parts. All available real QSOs at z > 0.3 were considered; objects with 0 < z ≤ 0.3 were instead undersampled to produce a more uniform redshift distribution. The down-sampling strategy was simple: low-z sources were chosen at random among all available targets. After down-sampling, ∼350 QSOs with 0 ≤ z ≤ 0.3 were kept.

The PRF was then trained on the combined data set (i.e. a combination of synthetic and spectroscopically confirmed QSOs) or, alternatively, on the subset containing only spectroscopically confirmed or synthetic QSOs. The predictive performance, measured on a testing data set of only spectroscopically confirmed QSOs, were compared one with the other. Results are shown in Table 3. It is clear that a sample of both real and synthetic QSOs produces, in term of recall, the best possible result; improvements in the recall negatively impact the purity of the sample: none the less, since the main interest is to improve the completeness of the final sample, both synthetic and spectroscopically confirmed QSOs will be used in Section 6. We consider the trade-off acceptable, since the number of candidates is still relatively low (≲1000): a different approach might be needed when dealing with a larger number of targets.

Once the usefulness of synthetic data is established, the optimal number of generated QSOs to employ and the redshift threshold z_th need to be identified. In order to do so, multiple runs of the PRF were considered; for each of these, the number n_synth and redshift threshold over which synthetic data are introduced were varied, while keeping the same testing data set. Results for each pair are illustrated in Fig. 9.

Improvements in the recall are negligible when introducing more than ∼6000 synthetic QSOs in the training data set, independently from the chosen redshift threshold. As a matter of fact, employing more than ∼6000 synthetic QSOs degrades the predictive abilities of the PRF in terms of precision. With respect to the number of real QSOs available in the training data set, ∼6000 synthetic QSOs correspond to twice the number of real QSOs at all redshifts: in the following this is the proportion employed for choosing the number of synthetic data to add to the training data set. This result is coherent with the premises of this work: introducing too many synthetic QSOs, mainly at z > 2.5, produces a training data set skewed towards high redshift. As such, the performance of the algorithm will degrade when classifying low redshift targets, thus lowering the overall precision.

We tested the performance of the PRF using both magnitudes and colours in the training and testing data sets. We found comparable results: the main difference is, on average, a slightly higher precision measured when using colours (⁠|$\sim 2{{\ \rm per\ cent}}$|⁠) with respect to the selection performed with magnitudes. Given the higher purity, the preferred candidate list will be generated with the PRF trained on colour data.

5.2.3 Leave-one-out test

Recall and precision are useful to estimate the global performance of the PRF over the entire data set; due to the low number of very high-redshift QSOs in the data set they may however paint a biased picture, underestimating misclassification of these targets. To account for this effect a leave-one-out test was performed. A leave-one-out test consists in training an algorithm on all available data but one object. The latter, left out from the training set, is considered as test. The process is then repeated iteratively for all available objects. To reduce the computation time needed we restrict the QSOs used for the leave-one-out to the 320 available targets at z > 3.0. The PRF was thus trained 320 times; in each of these 320 times, the training data set consists of synthetic data on top of all real QSOs, except one; the testing data set is just the one of the 320 QSO left out. This procedure allows to test the PRF on all QSOs at z > 3, one by one. The exercise was repeated with the oversampling approach employed in Paper I as well (see Paper I, end of Section 2).

Out of the 320 available QSOs at z > 3, we find that the PRF misclassifies 4 targets when trained on synthetic data, 13 with the method described in Paper I. Considering the results at face value, both approaches work rather well, with at most 13 (⁠|$\sim 4{{\ \rm per\ cent}}$| of the total) misclassifications. A more careful analysis, however, shows that the algorithm trained with synthetic data produces better results:

• the number of misclassified targets is lower (4, compared to the 13 misclassified registered with the method outlined in Paper I). 3 of the misclassifications are in common: two of them have peculiar SED, which are not similar to the typical SED of high-z QSOs (Fig. 10). IR spectra would be required to understand whether this peculiarity is intrinsic or, e.g. due to effect of variability. The SED of the third object, on the other hand, is rather typical; the most likely explanation for its misclassification as low-z QSO is a detection in the SkyMapper v band, extremely uncommon in the sample for quasars with z ∼ 3, possibly an indication that the line of sight towards this target is particularly unabsorbed.

• the method outlined in paper I misclassifies 4 out of the 13 available targets at z > 4.5, which are instead recovered by the PRF trained on synthetic data;

• a bright, high-z target, recently identified with a novel selection approach (Calderone et al., in preparation) is not selected by the method outlined in paper I, while it is correctly classified when training the PRF with synthetic data. Misclassifying this target is especially undesirable, as it is both bright and at high redshift: it is thus an ideal candidate for the redshift drift Sandage test (e.g. Liske et al. 2008; Boutsia et al. 2020), the search of variation for fundamental constants (e.g. Milaković et al. 2021; Murphy et al. 2022), or any other high resolution, precision, and stability experiment to be carried out in the future.

The Leave-One-Out test was repeated with the PRF trained both on colours and on magnitudes, finding the same results reported above.

5.3 Literature misclassifications

As a byproduct of the tests described in the previous sections we identified some misclassifications in published catalogues. We report them in Table 4, together with a revised redshift and/or class. An example of a spectrum for a misclassified target is shown in Fig. 11.

These targets were identified by visual inspection of the misclassified objects in a given testing/validation data set. After each evaluation run, the spectral energy distribution and, if available, the spectrum of misclassified targets at z > 2.5 is visually inspected. This is done to try to understand the reason for the misclassification and to check whether the target is genuinely a high-redshift QSO or was assigned an incorrect redshift estimate or class in the original catalogue. Following this approach, 4 targets with an available spectrum were re-classified or given a new redshift estimate.

Additionally, outliers identified as objects with a high reduced χ² of the synthetic photometry with respect to the observed magnitudes while building the composite SED (Section 4.1) were visually inspected, resulting in the discovery of 5 additional errors, in the redshift determination and/or in the object classification; they are also reported in Table 4.

6 THE CANDIDATE SAMPLES

The previous sections described the tests performed on known data in order to understand how well the algorithm performs in the selection of high-redshift QSOs. We now turn to the application of the algorithm on unlabelled data, in order to extract a new candidate list.

Two lists were produced: one obtained from a training sample where magnitudes (both of real and synthetic QSOs) are used as features; another where, instead, colours are employed. Based on the results presented in the previous sections, the two lists are expected to be similar.

After removing likely extended targets (Calderone et al. 2019), 31 550 out of the 1147 017 sources in the original main sample are still lacking a classification label. As a first step, the PRF is trained to separate QSOs, stars and non-active galaxies (Section 5.2.1). When applied on the data set of unclassified sources, the algorithm produces a list of 11 612 (11 696) QSO candidates at all redshifts if trained on magnitudes (colours), 19 198 (19 191) star candidates and 732 (663) non-active galaxy candidates. The lists of stars and non-active galaxies are discarded: they are, for our purposes, contaminants. During the second step, and given the sample of QSO candidates at all redshift, the algorithm is trained on QSO data in order to distinguish low- and high-redshift targets (Section 5.2.2). The list consists of 1028 candidates if trained with magnitudes, 898 when the algorithm is trained with colours; the two lists share ∼85 per cent of the targets.

Both selections were performed using the optimal hyperparameters identified in Section 5. Consistently with the results discussed in previous sections we find similar completeness (⁠|$\sim 87{{\ \rm per\ cent}}$|⁠) for both lists, but slightly lower contamination from the candidate list generated by training the PRF on colour data. This is consistent with the number of candidates recovered, assuming that a large part of the non-overlapping candidates are contaminants (either galaxies, stars or low-redshift targets).

As a consistency check, we compared the colour distribution of the high-redshift QSO candidates with that of the training sample. We expect the distribution of the training set, which includes synthetic data, to extend to larger values of a given colour: synthetic data were introduced to provide a more extensive coverage in the colour space. The i − z colour was chosen, as all sources in the main sample are required to have a i magnitude measurement, and all but 32 targets, mostly stars, have a z magnitude measurement. Fig. 12 compares the distribution of colours for the training set and the high-z candidate list, obtained when training the PRF on magnitudes: as expected, the two are consistent with each other but with a wider distribution for the training set. Similar results are found when training the PRF on colours, as most of the candidates are in common between the two lists.

7 SPECTROSCOPIC VALIDATION AND CATALOGUE

New spectra of quasar candidates are continually taken in the framework of the QUBRICS program. To maximize the success rate and reduce the contamination, targets are selected following an ‘iterative’ approach: after each spectroscopic run the training sample is updated to include the new identifications. In this way the training set and the list of candidates is evolving after each run: for the purposes of this paper they have been frozen at the situation at the end of February 2022.

In Table A1, we list the new spectroscopic identifications – not reported in the previous QUBRICS papers – of PRF candidates from the list of Paper I and updated lists based on the SkyMapper DR3. They are part of the training sample described in this work.

A total of 206 PRF-selected candidates have been observed at Las Campanas Observatory, Telescopio Nazionale Galileo (TNG, La Palma) and ESO-NTT using, respectively, the LDSS-3 (Clay Telescope), DOLORES, and EFOSC-2 spectrographs. Table 5 summarizes the observing setups and significant information about each observing run. Fig. 13 shows some selected spectra drawn from the catalogue.

109 of the 206 presented candidates have been observed with LDSS-3 at the Clay Telescope. Observations were obtained in several nights with varied conditions (e.g. bright time, variable weather conditions). The VPH-all grism with the 1 arcsec central slit and no blocking filter was employed, covering a wavelength range between 4000 and 10 000 Å with a low resolution of R ∼ 800. Exposure times ranging between 800 and 1800 s were used, depending on the candidate magnitude. Data obtained from the LDSS-3 instrument were reduced with a custom pipeline based on MIDAS (European Southern Observatory 2013) scripts. Each spectrum has been processed to subtract the bias and normalized by the flat; wavelength calibration is achieved using helium, neon, and argon lamps, finding an rms of ∼0.5 Å. Observing conditions have not always been photometric: the flux calibration is thus relative to a spectrophotometric standard star, observed each night. Eighty candidates were observed with the DOLORES instrument mounted on the Telescopio Nazionale Galileo. Exposures have been taken during the AOT43 and AOT44 periods, from April to November, under two proposals (PI: F. Guarneri); the LR-B grism (resolution ∼600) with a 1 arcsec slit aperture was used with an exposure time between 300 and 600 s.

In November 2021, we were awarded four nights (PI. F. Guarneri, proposal 108.22L1.001) at NTT, employing the EFOSC-2 instrument Grism #13 (wavelength range λ ∼ 3700–9300 Å), with typical exposure times ranging between 400 and 800 s; 17 candidates from the list obtained with Paper I and this work were observed and given a robust identification and redshift. The same data reduction procedure applied for LDSS-3 data was used for the EFOSC-2 spectra.

Out of the 206 targets, 149 turned out to be genuine QSOs with z > 2.5, 41 QSOs with z < 2.5 (<z> = 2.18), 3 galaxies and 13 stars. Based on the results of these runs, the achieved success rate is of 72 per cent, with the most significant contaminant being low-redshift QSOs. Among the 41 low-redshift QSOs, 7 show broad absorption features in their spectra, and have been marked as * in Table A1; BAL features mimic the colours of high-z QSOs, making it hard to select against these particular targets (Cupani et al. 2022).

8 ESTIMATED SUCCESS RATE OF THE REMAINING CANDIDATES

The precision of our algorithm for high-z QSOs selection has been estimated in Section 5. However, if we were to observe all the remaining candidates, we expect a lower precision since (i) during candidate prioritization we already observed the most promising ones and (ii) the adoption of synthetic spectra for training may bias our estimate. In this section, we will estimate a more reliable precision for the remaining candidates; this estimate is equivalent to the success rate expected on the candidate sample.

As mentioned in Section 2, the PRF computes, for each object, both a class label and the probability of belonging to the assigned class. This probability can be interpreted as a reliability score: the higher the probability, the most likely a candidate is an actual high-z QSO. Comparing the probabilities distribution of high-z QSO candidates with those of a testing data set (Fig. 14) it is clear from the blue histogram that a significant part of the objects in the candidate list have a lower probability of being high-z QSOs with respect to the testing data set. In order to keep this into account we use the QSO luminosity function to estimate the expected number of high-z QSOs in the candidate list and estimate the success rate of a telescope run observing all the 898 targets selected in the previous section with the PRF trained on colour data,

As in paper I, we assume that the PRF is able to exactly select QSOs, stars and non-active galaxies.

We first try to estimate how many non-QSOs are expected to contaminate the candidate list. The PRF is thus trained to distinguish QSOs, stars, and non-active galaxies. Applied on a testing data set (the same used in Section 5.2.1), a small fraction of the stars (⁠|$\lesssim 1{{\ \rm per\ cent}}$|⁠) and galaxies (∼6.5 per cent) into it are misclassified as QSOs. Following the steps described in the previous sections these sources, actually contaminants, would then be re-processed to identify low- and high-redshift targets. As such, the algorithm is trained to distinguish low- and high-z QSOs, and used to label the stars and non-active galaxies just misclassified as QSOs. This leads to |$\sim 1{{\ \rm per\ cent}}$| of the galaxies and |$\lesssim 0.1{{\ \rm per\ cent}}$| of the stars initially present in the testing data set to be labelled as high-redshift QSO. We then multiply these estimates with the predicted number of stars and galaxies among the unclassified sources. This should provide the number of stars and galaxies expected to contaminate the candidate list, which turns out to be ∼10 galaxies and ∼20 stars.

The most significant contaminant is however represented by low-redshift QSOs (Calderone et al. 2019; Boutsia et al. 2020): in Paper I, these were estimated to be |$\sim 85{{\ \rm per\ cent}}$| of all contaminants. In order to assess their significance for the current candidate list, we multiply the expected number of QSOs as a function of redshift, derived from Kulkarni, Worseck & Hennawi (2019) and corrected at high-redshift with Grazian et al. (2022) and Boutsia et al. (2021) with the misclassification rate. We define the latter as the ratio between the number of misclassified QSOs in a given redshift bin and the total number of QSOs in the same redshift bin. In order to calculate the misclassification rate, we train the PRF to separate low- and high-redshift QSOs; to avoid biases due to a particular choice of training/testing data sets, 100 different iterations were considered, each with a different training/testing split. We calculate the misclassification rate separately for each of these runs, and a global misclassification rate as the median across all 100 iterations; an uncertainty is estimated as well by using the median absolute deviation. The results are shown in Fig. 15. Most of the misclassifications happen at 2 ≤ z ≤ 2.5, close to the threshold separating the two classes.

Keeping into account the misclassification rate shown in Fig. 15 and that, as suggested by Boutsia et al. (2021), SkyMapper has an inherent incompleteness of |$\sim 15{{\ \rm per\ cent}}$| with respect to QSO data, we expect to find ∼590 low-redshift QSOs in the candidate sample, mostly in the 2 < z < 2.5 redshift range. In this way, the success rate for the candidates still to be observed is expected to be around |$31{{\ \rm per\ cent}}$| with about 278 new QSOs at z > 2.5 and i < 18 to be found.

9 DISCUSSION AND CONCLUSION

The probabilistic random forest is an effective tool in selecting bright (i < 18), and high-redshift (z > 2.5) QSOs: in Guarneri et al. (2021) we presented a first approach to the selection of these targets based on photometric data. In this paper, we aimed at improving the selection method outlined in Paper I by exploring the effect of synthetic data used as training set. Synthetic data are commonly used when training machine learning algorithms, and their role is expected to become even more prominent in future years (e.g. Nikolenko 2021). Diversified training examples can make an ML algorithm more precise, including corner cases and under-represented samples to make up for the lack of examples found in real-world data sets. Bright, high-z QSOs are a clear example of under-represented targets: training an ML algorithm for identifying them can prove difficult; introducing synthetic data significantly eases the task.

In this work synthetic data have been generated based on the mean properties of previously spectroscopically known QSOs, including new QUBRICS identifications. A composite spectral energy distribution was computed, and used to obtain synthetic magnitude measurements by multiplying the estimated fluxes by instrumental pass-bands. The candidate selection was split in a two-step process: in the first step we aim at sifting QSOs, at all redshifts, from stars and galaxies. In the second, we re-classify QSO candidates in order to distinguish low- and high-redshift targets. Synthetic data are used for training the algorithm during this second step. In all cases, we use precision and recall as evaluation metrics for the algorithm.

• The comparison of the performances of the new approach with respect to Paper I shows significant improvements in recovering a large part of the bright and/or very high redshift targets. This is especially significant, as bright and z > 3 QSOs are the most precious targets for future experiments, such as the search for varying fundamental constants (Milaković et al. 2021; Murphy et al. 2022) or the redshift drift experiment (Liske et al. 2008; Boutsia et al. 2020). The improvements in recall cost something in terms of precision, lowering the overall score from 78 to 64 per cent. Since we are mostly interested in improving the recall and the number of QSO candidates remains manageable, we consider this trade-off acceptable;

• the effect of colours on the selection was tested, with and without synthetic data. We found slightly better results when colours are adopted, mainly in the precision metric, with respect to those obtained when using magnitudes as features. Similarly to the results reported above, synthetic data significantly improve the measured recall also when colours are adopted.

We stress that precision and recall are only used to evaluate the performance of the algorithm and should not be intended as, respectively, success rate of the survey and completeness of the final sample. In particular:

• we use precision as the evaluation metric, but note that the expected success rate on the remaining candidates may be lower (Section 8). The two are not expected to be equal to each other: due to the iterative nature of the survey, as the survey progresses new spectroscopically confirmed QSOs are included in the training set and are thus removed from the candidate list. The success rate measured on the sample will thus drop with time, as more spectroscopically confirmations are carried out;

• we chose the recall, rather than the QSO completeness, as metric to characterize the algorithm since the former depends only on the selection performances, while the latter depends also on the MS selection. More specifically, in order to estimate the completeness we would also needs to keep into account the number of quasars excluded a priori from the MS due to the requirements mentioned in Section 2.2. According to Boutsia et al. (2021), the MS has an inherent incompleteness of 15 per cent with respect to QSO data, we expect the completeness of the selection to be ∼70 per cent, i.e. 0.85 times the recall, estimated to be ∼85 per cent.

The selection method outlined in this work was applied on unlabelled data from a new main sample, based on the third data release of the SkyMapper survey. The data include photometric measurements from the SkyMapper Dr3, Gaia eDR3, AllWISE, and 2MASS surveys. When trained on colour data, the algorithm produces a list of 898 candidates.

A large number of spectroscopic follow-ups have been carried out since Paper I, and 206 targets have been observed. These 206 have been selected in the past by the PRF, and are now part of the training sample: in order to maximise the success rate, we follow an iterative approach, updating the QSOs candidate list after each observing run. Of the 206 targets, 149 turned out to be genuine high-z QSOs. About 278 more QSOs at z > 2.5 and i < 18 are expected to be contained in our present candidate list.

We are further refining the selection method and continuing the spectroscopic campaign, in order to expand the spectroscopic sample and provide more training to our machine learning algorithms. Besides, novel selection techniques are being developed and currently tested, while new photometric data sets are also being employed (Calderone et al., in preparation).

SUPPORTING INFORMATION

suppl_data

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

ACKNOWLEDGEMENTS

We thank M. Fumagalli, J. X. Prochanska, and J. Hennawi for the SpecDB data base (Prochaska 2017), and for kindly providing spectra of previously observed targets. We thank an anonymous referee for a timely and thorough report that helped to clarify various aspects of the paper. We thank Società Astronomica Italiana (SAIt), Ennio Poretti, Gloria Andreuzzi, Marco Pedani, Vittoria Altomonte and Andrea Cama for the observation support at TNG. Part of the observations discussed in this work are based on observations made with the Italian Telescopio Nazionale Galileo (TNG) operated on the island of La Palma by the Fundacion Galileo Galilei of the INAF (Istituto Nazionale di Astrofisica) at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias.

The national facility capability for SkyMapper has been funded through ARC LIEF grant LE130100104 from the Australian Research Council, awarded to the University of Sydney, the Australian National University, Swinburne University of Technology, the University of Queensland, the University of Western Australia, the University of Melbourne, Curtin University of Technology, Monash University, and the Australian Astronomical Observatory. SkyMapper is owned and operated by the Australian National University’s Research School of Astronomy and Astrophysics. The survey data have been processed and provided by the SkyMapper Team at ANU. The SkyMapper node of the All-Sky Virtual Observatory (ASVO) is hosted at the National Computational Infrastructure (NCI). Development and support the SkyMapper node of the ASVO has been funded in part by Astronomy Australia Limited (AAL) and the Australian Government through the Commonwealth’s Education Investment Fund (EIF) and National Collaborative Research Infrastructure Strategy (NCRIS), particularly the National eResearch Collaboration Toolsand Resources (NeCTAR) and the Australian National Data Service Projects (ANDS).

This work is based on data products from observations made with European Southern Observatory (ESO) Telescopes at La Silla Paranal Observatory under ESO program ID 108.22L1.

This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation.

This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration.

This paper includes data gathered with the 6.5 meter Magellan Telescopes located at Las Campanas Observatory, Chile.

This work made use of Astrocook (Cupani et al. 2020), astropy (Astropy Collaboration 2013, 2018, 2022), IPython (Pérez & Granger 2007), the Julia programming language (Bezanson et al. 2017), Jupyter Notebooks (Kluyver et al. 2016), matplotlib (Hunter 2007), Numpy (Harris et al. 2020), Sklearn (Pedregosa et al. 2011), and tqdm (da Costa-Luis et al. 2022).

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

REFERENCES

APPENDIX A: SPECTROSCOPIC CATALOGUE