Machine Learning-based Search of High-redshift Quasars
Guangping Ye Huanian Zhang Qingwen Wu







Department of Astronomy, Huazhong University of Science and Technology, Wuhan, Hubei 430074, China
guangping@hust.edu.cn


Department of Astronomy, Huazhong University of Science and Technology, Wuhan, Hubei 430074, China
Steward Observatory, University of Arizona, Tucson, AZ 85719, USA
huanian@hust.edu.cn


Department of Astronomy, Huazhong University of Science and Technology, Wuhan, Hubei 430074, China


§ ABSTRACT 

We present a machine learning search for high-redshift (5.0 < z < 6.5) quasars using the combined photometric data from the DESI Imaging Legacy Surveys and the WISE survey. We explore the imputation of missing values for high-redshift quasars, discuss the feature selections, compare different machine learning algorithms, and investigate the selections of class ensemble for the training sample, then we find that the random forest model is very effective in separating the high-redshift quasars from various contaminators. The 11-class random forest model can achieve a precision of 96.43% and a recall of 91.53% for high-redshift quasars for the test set. We demonstrate that the completeness of the high-redshift quasars can reach as high as 82.20%. The final catalog consists of 216,949 high-redshift quasar candidates with 476 high probable ones in the entire Legacy Surveys DR9 footprint, and we make the catalog publicly available. Using MUSE and DESI-EDR public spectra, we find that 14 true high-redshift quasars (11 in the training sample) out of 21 candidates are correctly identified for MUSE, and 20 true high-redshift quasars (11 in the training sample) out of 21 candidates are correctly identified for DESI-EDR. Additionally, we estimate photometric redshift for the high-redshift quasar candidates using random forest regression model with a high precision. 





























§ INTRODUCTION
 

Quasars in general are driven by a supermassive black hole (SMBH) at the centre of the host galaxy through a process of accretion, and are the brightest non-transient sources of light in the Universe. SMBH activities are a key ingredient of galaxy formation, and are critical to subsequent galaxy evolution. Quasars at z>5 are often referred to as high-redshift quasars since they are at the end of the reionization epoch of the Universe, when the vast majority of the Universe's neutral hydrogen has been reionised <cit.>. High-redshift quasars provide effective probes for the study of galaxy evolution and cosmology, including the evolution of the intergalactic medium <cit.> and the circumgalactic medium <cit.>, the formation of early supermassive black holes, the co-evolution of SMBH and their host galaxies <cit.>. 

There exist many challenges in searching and verifying high-redshift quasars. On one hand, the number density of high-redshift quasars is very low in the Universe <cit.>, resulting to a small number of observable high-redshift quasars, which means that the single-object spectroscopic observations on large-aperture telescope of quasar candidates are needed and expensive; on the other hand, the contamination is overwhelming and could be a few orders of magnitude more than signals <cit.>. The major contaminations include cool galactic dwarfs with spectral types of M, L, and T <cit.>, which are similar to those of high-redshift quasars in the color space constructed from optical and near-infrared broad-band photometry. 

The traditional method for searching the high-redshift quasars is the color-cut selection <cit.> based on color drop-out, which is caused
by the Lyα break in quasar spectra due to significant IGM
absorption at the wavelength blueward of the Lyα emission
line at the rest frame of the quasar. This method successfully discovers the majority of the currently known high-redshift quasars <cit.>. The big advantage of this method is that it leads to well-defined selections that are easily reproducible and can be justified with known physics (e.g., the redshift evolution of the Lyα emission through the broadband filters and the drop-out due to neutral IGM absorption). However, color-cut selections might not make use of all the available information in the high-dimensional color space and the two-dimensional color-color diagram might be misleading or biased. And the strict color cuts might result in missing quasars due to the scattering out of the selection regions based on color-color cuts. In contrast to the color-cut selection method, machine learning-based methods can make full use of the color information and construct correlations in the high-dimensional space. Moreover, the large imaging surveys collect data far more than what can be handled manually, machine learning-based automatic methods can easily process those data and can be easily applied to new data. 

So far, there have been many successful examples of machine learning algorithm-based searches for quasars up to redshift z ∼ 6 quasars <cit.>. A variety of supervised machine-learning algorithms have been successfully applied for quasar selections, such as random forest <cit.>, Support Vector Machines (SVM), XGBoost, and Artificial Neural Networks <cit.>. <cit.> successfully applied the random forest algorithm to the high-redshift quasars (4.8 < z < 6.3) search using the combined dataset of Pan-STARRS DR1 <cit.> and ALLWISE <cit.>. Although the y band from Pan-STARRS is essential for high-redshift quasars at z > 6.5, it is comparably shallow (3π stack 5σ depth is ∼ 21.4), which might not be efficient in capturing faint high-redshift quasars. Here in this study we will use data from Dark Energy Spectroscopic Instrument (DESI) Legacy Imaging Surveys <cit.>, which is roughly 2 magnitudes deeper than the 3π stack 5σ depth of y band of Pan-STARRS. We will demonstrate that the constructed features other than the simple photometry provided in the catalog are critical to a successful machine learning model for high-redshift quasars search. 

The aim of this paper is to search new high redshift quasars from the Legacy Survey data <cit.> combined with Wide-field Infrared Survey Explorer (WISE) all-sky survey data <cit.> based on machine learning algorithms. 
The structure of the paper is organized as follows. In Sec. <ref> we introduce the catalog data we will use in this study and the training sample for the machine learning algorithm. In Sec. <ref> we briefly discuss the various aspects of the machine learning algorithms. In Sec. <ref> we estimate the photo-z of the high-redshift quasar candidates. In Sec. <ref> we discuss the procedure to obtain the high-redshift quasar candidates, the spectral verification of the candidates and present the final catalog. Finally we summarize. The magnitudes used in this paper are generated based on the AB system after applying Galactic extinction corrections. We use a ΛCDM cosmology with Ω_Λ = 0.7, Ω_m = 0.3 and H_0 = 70  km s^-1Mpc^-1.




§ DATA
 


 §.§ Legacy Survey Data
 

We use the photometric data from the Data Release 9 (DR9) of the Legacy Survey, which includes observations obtained by the DECam at the CTIO 4 m (DECaLS), an upgraded MOSAIC camera at the KPNO 4 m telescope (MzLS, Mayall z-band Legacy Survey), and the 90Prime camera <cit.> at the Steward Observatory 2.3 m telescope (BASS, Beijing-Arizona Sky Survey). Briefly, the Legacy Survey <cit.> was initiated to provide targets for the DESI survey drawn from deep, three-band (g=24.7, r=23.9, and z=23.0 AB mag, 5σ point-source limits) images, roughly two magnitudes deeper than the Sloan Digital Sky Survey (SDSS) data. The survey covers about 14,000 deg^2 of sky visible from the northern hemisphere between declinations approximately bounded by -18^∘ and +84^∘. The depth of SDSS <cit.> and Pan-STARRS <cit.>) is insufficient to provide reliable DESI targets. The footprint of DR9 (Figure <ref>) also includes an additional 6,000 deg^2 extending down to -68^∘ imaged at the CTIO by the Dark Energy Survey <cit.>. Furthermore, the Legacy Survey also includes g, r, z band images in the North Galactic Cap (NGC) region that overlaps with DES, extending the coverage of the Legacy Survey to over 20,000 square degrees. 








 §.§ The Wide-field Infrared Survey Explorer Data
 

In addition to the Legacy Survey data, we also take advantage of the WISE data release, providing infrared photometry of four bands at central wavelength of 3.4, 4.6, 12, and 22 μm (W1, W2, W3, W4) over the entire sky <cit.>. The W1, W2, W3 and W4 data in the Legacy Survey catalog are based on forced photometry in the unWISE stacked images (including all imaging through year 7 of NEOWISE-Reactivation) at the locations of the Legacy Survey optical sources. For our selection process we only focus on the W1 (3.4 μm) and W2 (4.6 μm) photometry with limiting magnitudes of 19.8, 19.0 in AB magnitude (Vega: 17.1, 15.7), respectively. The conversion between AB and Vega magnitude for W1 and W2 is W1_ AB = W1_ Vega + 2.699 and W2_ AB = W2_ Vega + 3.339, respectively.



 §.§ Training sample
 


We use supervised machine learning algorithms to perform the classification, therefore the training sample is critical to construct a good and reliable model. Any bias in the training sample will propagate to the new data, resulting in highly unreliable outputs. In order to construct a representative and unbiased training sample, we collect the spectra-confirmed signal (high-redshift quasars at 5.0 < z < 6.5) and all possible contaminators which populate the same or similar color space (eg. r-z, z-W1, z-W2) because of similar absorption features in the optical and near-infrared bands for both the signal and the contaminators. For quasars at z ∼ 5, Lyman series absorption systems begin to dominate in the r band and Lyα emission moves to i or z band, resulting in redder r-z color. However, as the redshift increases, most z > 5 quasars enter the MLT dwarfs locus in g-r/r-z color-color diagram (as shown in the upper-left subplot of Figure <ref>). At z ≳ 5.7, Lyα emission line moves into the z band and the g, r, i bands become the drop-out bands, resulting in very red r-z color and comparably red z-W1 color. And the infrared photometry W1 and W2 (constructed colors such as z-W1, z-W2, W1-W2) are useful for all redshifts. The medium-redshift quasars with redshift just below 5, the quasars at lower redshift, the MLT dwarfs, the different types of stars also share similar color-space with the high-redshift quasars, consisting of the major contaminants for the signal of high-redshift quasars.

We then obtain the photometry data from both Legacy Survey (g, r, z bands) and WISE (W1, W2 bands) for the spectra-confirmed high-redshift quasars (5.0 < z < 6.5), and the spectra-confirmed contaminators, which include the MLT brown dwarfs, AFGK stars, medium-redshift quasars (3.5 < z < 5.0), low-redshift quasars (1.5 < z < 3.5), and very low-redshift quasars (0 < z < 1.5). Although there are a number of quasars at z > 6.5, we will not include those since the Lyα emission of those quasars is already beyond the z band coverage of the Legacy Survey. We also do not include O- and B-type stars, because they are far from the high-redshift quasars in color space <cit.>. The distant compact early-type galaxies at intermediate or higher redshift, exhibiting very red r-z color, might be a possible type of contaminants, which will be further discussed in Sec. <ref>.







The sample of high-redshift quasars is extracted from the Supplementary Database[https://www.annualreviews.org/content/journals/10.1146/annurev-astro-052920-102455#supplementary_dataSupplementaryDatabase.org/] <cit.>, which are constructed from previous observations <cit.>, and the DESI spectra-confirmed quasars at z > 5[https://cdsarc.cds.unistra.fr/viz-bin/cat/J/ApJS/269/27#/browseDESIHighzQSO.org/] <cit.>. 
The MLT brown dwarfs are extracted from both UltracoolSheet[http://bit.ly/UltracoolSheetUltracoolSheet.org/] <cit.>
and SDSS DR16 Database[https://www.sdss.org/dr16sdss.org/dr16] <cit.>. 
The quasars at z < 5 are extracted from SDSS quasars catalog DR16[https://data.sdss.org/sas/dr16/eboss/qso/DR16Q/data.sdss.org/sas/dr16/eboss/qso/DR16Q]<cit.>. 
The AFGK stars are obtained from the SDSS DR16 Database <cit.>. 
To obtain a training sample as representative as possible, we remove the entries whose z band magnitude is either too bright (top 1%) or too faint (last 1%) for each category, which are not representative among the whole training sample. The number of instances for each category in the training sample is shown in Table <ref> and the redshift distribution of the quasars is shown in Figure <ref>. Note that the data size for each class in the training sample is not equal, which means that we are dealing with the imbalanced classification. This issue will be further discussed in detail in Sec. <ref>.




The Legacy Survey catalog data also provides fluxes at different aperture radius for each band, calling `apflux', which is useful because it provides more features for the machine learning model and essential to discriminate the signal from extended but compact sources like distant compact galaxies. The radius of the apertures are [0.5, 0.75, 1.0, 1.5, 2.0, 3.5, 5.0, 7.0] arcseconds for the g, r, z bands and [3, 5, 7, 9, 11] arcseconds for the W1, W2 bands. We will construct apflux ratios as the machine learning features, with a detailed description  in Sec. <ref>.


Moreover, the combined fluxes grz and W, which are constructed from g, r, z fluxes and W1, W2 fluxes, respectively, are useful in quasars search according to a recent study <cit.>. The definitions of the combined fluxes of grz and W and the conversion between flux and magnitude are:

    
    flux_grz = (flux_g + 0.8*flux_r + 0.5*flux_z)/2.3
	

    
    flux_W = 0.75*flux_W1 + 0.25*flux_W2
	

    
    mag = 22.5 - 2.5*lg(flux)
	





The measurements we will use in the machine learning model include: g, r, z, grz, W1, W2, W fluxes and magnitudes, as well as apfluxes. Some selected color-color distributions of the training sample are shown in Figure <ref>. The top panels are based on the data set with no missing values of all measurements for high-redshift quasars (252 sources), and for all background sources. The bottom panels are based on the entire dataset, where the missing values of high-redshift quasars have been filled in. The process of filling in the missing values is described in detail in the following Sec. <ref>. As is shown clearly that some colors such as g-r, g-W2, r-grz are effective in separating high-redshift quasars from the contaminators. We would expect that the machine learning algorithm would significantly enhance the performance in separating the high-redshift quasars from the contaminators when making full use of all color information. It is also obvious that the MLT brown dwarfs and the mid-z quasars are the major contaminators since they overlap with the high-redshift quasars to a greater extent in the color space.
 


 §.§ Imputation of Missing Values




Currently, the number of spectra-confirmed quasars at z > 5 is 727, of which 602 quasars are in the Legacy Survey footprint. Because of the color drop-out, a large number of quasars (335) have nearly zero or even negative fluxes in the g, r bands, and a small fraction of quasars (29) also have nearly zero or negative flux measurements in the W1 and W2 bands due to the low sensitivity of WISE survey. There are no well-defined magnitudes when converting from zero/negative fluxes into magnitudes, which we call “missing values". The missing value issue highly limits the high-redshift quasars sample since the machine learning algorithm could not deal with missing values. Therefore, special considerations and techniques may be required to address this issue to ensure accurate analysis and interpretation of the available data.

To fully utilize the limited photometric information available for the known high-redshift quasars, we first employ common imputation methods to handle missing values in the training dataset. One common imputation method is to fill in zeros for the missing values, but it will highly bias the classification model because of the extreme color constructed from the filled-in zeros. Moreover, a zero AB mag has a physical meaning, so we will not discuss it here. The other common imputation methods include: 1) mean imputation, replacing missing values with the mean value of the respective feature; 2) LS limit imputation, substituting missing values with the limiting magnitudes of the Legacy Survey; 3) future imaging survey limit imputation(LR limit), adopting the 5σ point-source depth of 27.5 AB mag from the Large Synoptic Survey Telescope <cit.> for g, r, z bands <cit.> and of 24 from
the Roman Space Telescope <cit.> for W1, W2 bands <cit.> to replace the corresponding missing values, which will be further discussed in Sec. <ref>; 4) random forest (RF) imputation, training a random forest regression model to generate reasonable values for missing values; 5) Multiple Imputation by Chained Equations <cit.> forest, generating imputed values by iterative random forest model on multiple independent divisions of dataset. To validate the feasibility of the five methods mentioned above, we compare the distributions of the imputed values with the complete dataset without missing values and assess the differences.

Here we briefly introduce MICE. MICE is based on the random forest algorithm, and it constructs multiple different imputed datasets randomly to examine the uncertainty and other effects caused by missing values. Each dataset undergoes independent imputation processes, where the imputed values for each feature in the dataset are predicted using a random forest regression model trained on the non-missing dataset. The above process is repeated to update the imputed values. This iterative process does not stop until the average of the imputed values for each feature tends to converge. One big advantage of MICE is its ability to automatically select features related to the imputation variable without requiring manual specification. It also calculates weights for each variable to optimize the variable selection and model fitting during the imputation process. By combining the predictions from multiple random forest models, MICE provides accurate and reliable imputation results. 

We determine the best imputation method by comparing the density distributions of the magnitudes between the imputed dataset and the complete dataset. We use R^2 to represent the differences between the imputed datasets and the complete dataset for these features. R^2 assesses how well the model explains variable changes (0 to 1, closer to 1 for better fit). Figure <ref> illustrates the density distributions of the four magnitudes with missing values for the complete dataset as well as the five imputed datasets. For the g magnitude with more than half missing values (55.6%), the RF and MICE methods can achieve a R^2 value of > 0.90, while the results of the other methods are poor. It is easily seen that MICE performs significantly better than the RF model. For the r, W1, W2 magnitudes with fractions of missing values of 12.29%, 2.33%, 3.82%, respectively, the performance of RF and MICE are similar and are much better than the other methods. Note that the R^2 criteria could not reflect the effectiveness of the imputation method using future imaging survey limits, and we will further discuss it in Sec. <ref>.

Overall, MICE is an effective approach that leverages the power of random forest regression models for missing value imputation. It is very effective in that it automates feature selection and enhances the imputation process through weighted variables. The combination of multiple random forest models enables MICE to provide accurate and reliable imputation results. Therefore, we will adopt the MICE method to perform imputation in the following discussion. The robustness of the MICE imputation will be discussed in the following discussions.





§ MACHINE LEARNING
 

In the past decades, many large surveys such as SDSS <cit.>, Legacy Survey <cit.>, DESI <cit.> collected tons of data that is challenging to analyze. And the approach using automatic methods like machine learning is both efficient and easily reproducible. Here we briefly introduce the various aspects of the machine learning technique.




 §.§ Metrics
 

We already introduce the training sample in Sec. <ref>. Here we will introduce the other aspects of the machine learning technique. To evaluate the performance of a classification model, we typically consider the evaluation metrics including precision, recall, accuracy, f1 score, and the area under the Receiver Operating Characteristic curve (ROC AUC). In the evaluation process, the comparison between the predicted labels from the machine learning algorithms and the true labels yields four outcomes: true positives (TP), true negatives (TN), false positives (FP), and false negatives (FN), from which we can obtain the desired evaluation metrics.

Precision refers to the fraction of true positives in total positives predicted, defined as:


    precision = TP/TP + FP
	

Recall refers to the ratio of positive class instances correctly identified by a classifier, defined as:


    recall = TP/TP + FN
	

The f1 score metrics combines precision and recall. In fact, the f1 score is the harmonic mean of the two. A high f1 score symbolizes a high precision as well as a high recall. It is defined as:


    f1 = 2 ×precision×recall/precision + recall
	



The f1 score belongs to a unified formula f_β where β = 1. f_β is defined as:


    f_β = (1 + β^2)(precision×recall)/(β^2 ×precision) + recall
	

On datasets with balanced class quantities, the f1 score performs well. For highly imbalanced dataset, the Adjusted f-score (AGF) is an improvement upon the f_β score <cit.>. The AGF considers all elements of the original confusion matrix, making it a more equitable evaluation metrics for the classification of the minority class. It is defined as:


    AGF = √(f_β=2× inv f_β=0.5)

where the f_β=2 score weights recall more than precision, the inv f_β=0.5 refers to the f_β=0.5 score calculated using a confusion matrix obtained by inverting the labels, where TP is replaced by TN, and then substituted into f_β. 









In principle, the precision, recall, f1 as well as AGF is least affected by the imbalanced data sample. And it is almost impossible to tune the desired metrics for all classes for a multi-class classification task. Therefore, we will focus on the metrics of precision, recall, f1 and AGF for the signal of high-redshift quasars in the following discussions, because the correct identification of the high-redshift quasars from enumerous backgrounds is our focus in this study. Although there will be confusion between the classes of contaminants, such as the confusion between the different types of stars, they are not the focus of this study. But we will discuss the weighted metrics of all other classes except the high-redshift quasars for completeness. 






 §.§ Feature Selection
 


The feature selection, the appropriate algorithm and the choices of class ensemble is essential to a good machine learning model, all of which will be discussed in the following. Many studies <cit.> have demonstrated that random forest is very effective in separating signals from contaminators, therefore we will adopt the 11-class random forest classifier as our fiducial model to discuss the feature selection, and the 11 classes are shown in Table <ref>. And then we will compare the different machine learning algorithms and finally we will discuss whether the different class ensembles will help improve the model performance.

The measurements we will use in this study include dereddened magnitudes, dereddened fluxes, and aperture fluxes in the g, r, z, grz, W1, W2, and W bands. The aperture radius are [0.5, 0.75, 1.0, 1.5, 2.0, 3.5, 5.0, 7.0] arcseconds for the g, r, z bands, and are [3, 5, 7, 9, 11] arcseconds for the W1 and W2 bands. The apfluxes are essential to remove extended but compact sources like galaxies, which will be further discussed in Sec. <ref>. 
The features used in the machine learning algorithm are constructed from magnitudes (colors), dereddened fluxes (flux ratios) and aperture fluxes (flux ratios). The magnitude features include g, r, z, W1, W2, grz, W. The color features are constructed from the subtraction of two magnitudes. The flux ratios are constructed from the division of two fluxes, which is the same as the construction of color features. The apflux ratios are constructed as follows: 1) for apfluxes in the same band, we construct the apflux ratios by taking the ratio of aperture fluxes at adjacent radius indices; 2) for apfluxes across different bands, we construct the apflux ratios  by taking the ratio of aperture fluxes at the same radius index across adjacent bands for the first 5 radius for both optical (g, r, z) and the infrared (W1, W2) bands. For the extra three radius in the optical bands (g, r, z), the apflux ratios are constructed by the ratio of aperture fluxes at the same radius across different bands.

The apflux ratios can be represented by a general formula: ap_[a]_[i]/ap_[b]_[j], where [a] and [b] represent indices for the photometric bands (g, r, z, W1, W2). And the indices [i] and [j] representing the aperture radius, range from 1 to 8 for the optical bands and range from 1 to 5 for the infrared bands. The conditions are: 


    
  * [a] = [b] (the same band) and [i] - [j] = 1 (adjacent radius indices).
    
  * [a] - [b] = -1 (adjacent bands) and [i] = [j] (the same radius index).
















There are 7 magnitudes, 21 colors or flux ratios and 55 apflux ratios in total. To fully understand the contribution to improving the performance of the machine learning algorithm, we investigate the following subsamples with different feature subsets:


	
  * FeatureSet-A (7): g, r, z, W1, W2, grz, W magnitude features. 

	
  * FeatureSet-B (28): magnitude (7) and color (21) features.
	
	
  * FeatureSet-C (83): magnitude (7), color (21) and apflux ratio (55) features.
	
	
  * FeatureSet-D (83): flux (7), flux ratio (21) and apflux ratio (55) features.
	



FeatureSet-D is the mirror of FeatureSet-C and the big advantage of FeatureSet-D is that there is no missing value issue. We train the random forest classification model on the above four different datasets separately. In the model training, we use RandomizedSearchCV to find the best set of hyperparameters and we split the entire training sample into “training set" and “test set". Table <ref> presents the evaluation scores of the four feature sets for the class of high-redshift quasars. We also use cross-validation (CV) to avoid overfitting on the training sample, and the fold number is set to be 10. In Table <ref>, the “Val" column represents the average results of ten equally folds of the cross-validation datasets. We adopt the standard deviation within the ten folds as the error. The “Test" column shows the results of the models for the test set. 





The comparison between FeatureSet-A and B tells us that the addition of 21 color features significantly improve the performance for the validation data set. This also indicates that in high-dimension color space, the high-redshift quasars are separable from the contaminators. When incorporating the additional 55 apflux ratio into the models, we find that the performance is also enhanced for the test set. Additionally, we use pure flux values from different bands and their ratios as features in FeatureSet-D. Comparing the performance of FeatureSet-C and FeatureSet-D in the table, it can be observed that their precision are very similar, with FeatureSet-C having slightly higher recall, f1 and AGF. The consistency of all evaluation metrics between FeatureSet-C and FeatureSet-D as shown in Table <ref> demonstrates the robustness of our imputation method. 

Table <ref> presents top 20 features ranked by the importance for FeatureSet-C. We can see that the color features align with the trends shown in Figure <ref>. As can be seen clearly, the newly introduced features of grz and W magnitudes have comparably high importance, demonstrating that these additional features indeed play an important role in improving the model performance. We notice that the most important feature, z-W2, is not effective in separating the high-redshift quasars from the contamination shown in Figure <ref>, which will be further discussed in detail in Sec. <ref>.







 §.§ Comparison of Different Algorithms
 

There are many machine learning algorithms with their own strength and weakness, here we compare the performance of different machine learning algorithms on separating high-redshift quasars from contaminators. The classification algorithms used in this study are the following:




  * : A non-parametric algorithm that forms a model by placing the data from the training set in the high-dimensional feature space and makes predictions based on the similarity between the data point and its k nearest neighbours of the training set in the high-dimension feature space. The hyper-parameters in this algorithm include number of neighbors (n_neighbors), weights (w) and distance metric (p). The distance metric can be either Manhattan distance or Euclidean distance, used to measure the similarity.


  * : A non-parametric supervised learning that constructs a tree-like classification or regression structure by continuously dichotomising features based on a discrete set of values. In these tree structures, leaves represent class labels and branches represent conjunctions of features that lead to those class labels. The hyper-parameters in this algorithm include criterion for split quality, maximum depth of trees, minimum number of samples required to split an internal node, minimum number of samples required at each leaf node, maximum number of features to consider for each split, etc.


  * : An integrated approach that performs classification or regression by combining multiple decision tree models, where each tree is constructed using a random subset of the training data and a random subset of features. This randomness helps to reduce the overfitting issue and to improve generalisation. The hyper-parameters of this algorithm include number of estimators and the ones of the decision tree algorithm. 


  * : A gradient boosting algorithm that makes predictions by iteratively training multiple decision trees. Each decision tree is trained based on the residuals of the previous tree to gradually reduce the prediction error. The learning rate in LGBM algorithm is critical to the model performance. Other hyper-parameters are similar to the ones in random forest classifier.


  * : A plain Bayesian probabilistic classification algorithm assuming that features are independent of each other and that their conditional probabilities follow a Gaussian distribution. The prior probability (p) is the most important hyper-parameter of this algorithm. By default, p is estimated from the proportion of samples in each category in the training data.



We apply the above algorithms to FeatureSet-C and obtain the corresponding precision, recall, f1 and AGF scores of the high-redshift quasars for these models, as shown in Table <ref>. The maximum value for each evaluation metrics is highlighted. It can be seen clearly that the random forest model is the best one for all four evaluation metrics. Then we will adopt the random forest algorithm in the following discussion. 

Note that FeatureSet-C is the best feature set selected via the random forest model, but it might not be the best feature set for other algorithms. Here we investigate whether FeatureSet-A or FeatureSet-B on other algorithms have better performance than FeatureSet-C. We find that all metrics on FeatureSet-A for other algorithms are worse than those on FeatureSet-B. We present the best metrics using FeatureSet-B for the other four algorithms, we have the best precision (0.88 ± 0.05, from KNN), best recall (0.84 ± 0.06, from LGBM), best f1 score (0.83 ± 0.05, from LGBM) and best AGF (0.83± 0.05, from LGBM). It is obvious that the results obtained on FeatureSet-B for these four algorithms are all lower than the corresponding metrics on Featureset-C. Therefore, we conclude that Featureset-C is the best feature set not only for the random forest algorithm, but also for other classification algorithms considered in this work.






 §.§ Class Selection
 

We adopt the 11-class classification model in the above discussion, and the 11 classes are shown in Table <ref>. In order to investigate whether further reducing the number of classes in the training sample can improve the model's performance, we will now discuss four scenarios in which certain classes from FeatureSet-C are merged. The four scenarios are the following:


	
  * 11-class classification (P_11): 
		
		
  * P_11: vlowz, lowz, midz, highz quasars and M, L, T, A, F, G, K dwarfs.
		
	
	
  * 4-class classification (P_4): 
		
			
  * P_4^0: vlowz, lowz, midz quasars
			
  * P_4^1: highz quasars
			
  * P_4^2: M, L, T dwarfs
			
  * P_4^3: A, F, G, K stars 
				
	
  * 3-class classification (P_3):
 		
			
  * P_3^0: vlowz, lowz, midz quasars
			
  * P_3^1: highz quasars
			
  * P_3^2: M, L, T, A, F, G, K dwarfs
   				
	
  * 2-class classification (P_2):
 		
			
  * P_2^0: vlowz, lowz, midz quasars and M, L, T, A, F, G, K dwarfs
			
  * P_2^1: highz quasar
				



We apply the above four scenarios to FeatureSet-C and obtained four new datasets, called P_11, P_4, P_3, P_2. In P_4, P_3, P_2, the reason why we can combine vlowz, lowz, midz quasar samples is that they are in the period when the cosmic reionisation epoch is already over. In P_4 scenario, we separate the M, L, T dwarfs from other stars because they are the main contaminants of high-redshift quasars and they are quite different from other types of stars in the optical and near-infrared bands. M, L, T dwarfs are grouped together because of the similarity between them. A, F, G, K dwarfs are also grouped together because they are relatively far away from the high-redshift quasars on the color-color diagrams, as can be seen from Figure <ref>. 



We train the random forest classification model on the above four datasets. The precision, recall, f1 and AGF score of the high-redshift quasars for the four scenarios are presented in Table <ref>. Meanwhile, we also present the  average metrics weighted by the number of instances of each class  for all other classes except the high-z class for the test set. It is clearly seen that the precision and recall are consistent within error for the different scenarios identified by the number of classes in the dataset. This might indicate that there is certain correlation among the different contaminators, which is demonstrated by the weighted metrics of all other classes. The relatively stable precision indicates that the sample of high-redshift quasars is clean and possesses inherent features that is distinguishable from other contaminants, ensuring high precision and recall regardless of the number of classes in the classification models. As we reduce the complexity of the classification model by fewer classes via combining certain classes, the weighted metrics should get better
because the confusion between classes that we combine gets ignored, which is consistent with the results shown in Table <ref>.

Based on the statistics in Table <ref>, we decide to pick the 11-class model for two reasons. First of all, the 11-class model would provide us the full picture of contamination within the different classes, expecially the contamination for the high-redshift quasars. Secondly, we find that the strict boundary between “mid-z" and “high-z" redshift bins introduces confusion between these two classes, which should not be a concern. The inspection on those “mid-z" quasars, which are classified as “high-z" quasars, tells us that the redshift for those “mid-z" quasars are between 4.84 and 4.98, right around the redsfhit boundary of 5. Further inspection on the spectra of these “mid-z” quasars indicates that they bear a striking resemblance to the typical spectra of high-redshift quasars. This implies that these “mid-z" quasars which are predicted as the “high-z" class are not actually contaminants of the “high-z" quasars. Therefore, the precision of the “high-z" class for the 11-class model has been underestimated. If we exclude these “mid-z" quasars, the precision of of the “high-z" class for the 11-class model would reach 0.99. 




 §.§ Important feature: i-band photometry
 

For quasars with redshifts between 5 and 6.5, the Lyα emission line is shifted to the wavelength range of 7296 Å - 9120 Å. However, the Legacy Survey g, r bands does not cover this wavelength range, and the z band only has a small coverage in this wavelength range. Then the i band, whose wavelength coverage mostly overlaps with this wavelength range, is essential for high-redshift quasar (5 < z < 6.5) searching . Fortunately, the Legacy Survey DR10 not only includes images in the g, r, and z bands from DECaLS, but it also contains DECam observations for g, r, i, and z bands  from several non-DECaLS surveys, primarily the Dark Energy Survey, the DELVE Survey <cit.>, and the DeROSITA Survey[https://noirlab.edu/science/programs/ctio/instruments/Dark-Energy-Camera/DeROSITASDeROSITAS.org/]. These surveys mainly cover the southern sky (Declination ≤ 32.375^∘). The area covered in the i band for number of passes greater than 3 is 13,024 square degrees, while the total area covered in the g, r, i, and z bands jointly is 9,923 square degrees.

Because the combined sky coverage of the g, r, i, and z bands is approximately half of that for the g, r, z bands, we divide the entire dataset into the one that includes the i band data (FeatureSet-i) and the one that does not include the i band data (FeatureSet-non-i). The feature construction is the same as above except that the FeatureSet-i dataset includes all i band related measurements. We apply MICE to fill in missing values as discussed in Sec.<ref> and we train models using the 11-class 
random forest classification algorithm on these two datasets. We evaluate the models on the validation set and test set to assess the i band importance.





Table <ref> shows the evaluation scores of validation and test sets for the two datasets (FeatureSet-i v.s. FeatureSet-non-i), as well as the first few most important features of the models trained on FeatureSet-i. Compared to FeatureSet-non-i, which does not include features related to the i band, the models trained on FeatureSet-i achieve higher precision on the test set, while the other evaluation scores are similar. This indicates that the i band can improve the performance of our model in searching for high-redshift quasar to some extent. The columns of “Feature" and “Importance" display the top-ranked features and the corresponding importance in the models trained on FeatureSet-i. It can be observed that the feature g-i has high importance, which reflects the significance of i band related features in our model.


In this section, we already demonstrate that the i band and its corresponding features are indeed playing an important role in the classification model for high-redshift quasars. However, due to the limited amount of i band data, our model is ultimately trained on the data provided by Legacy Survey DR9. In future work, we will consider incorporating photometric data from more bands to obtain additional features and further improve the model's performance.




 §.§ Discussion of Random Forest performance
 


After the discussion on feature selections, the comparison of different classification algorithms and class selections of the dataset, we find that the 11-class random forest classification model on the FeatureSet-C and D yields the best performance for separating the high-redshift quasars from the contaminators. Figure <ref> shows the confusion matrix of the model using FeatureSet-C for the test set, where each row represents the true class and each column represents the predicted class. The percentages in each square indicate the fraction of the corresponding class in each column. Squares with a percentage of less than 1% are not displayed. For the diagonal elements, these percentages represent the precision of the corresponding class. 


 
For the class of high-redshift quasars, the precision achieves as high as 96.43% and the recall can also reach 91.53%. The precision of the high-redshift quasars for this model is significantly enhanced compared to the previous work <cit.>. The high precision ensures that the final high-redshift quasar candidates have much larger probability to be true high-redshift quasars, which significantly improves the efficiency of the future spectra verification. The recall is also significantly enhanced compared to the previous work <cit.>, which means much less true high-redshift quasars will be missed. The average precision and recall weighted by the number of instances of each class for all other classes except the signal can also reach 86% and 87%, respectively, which are also significantly higher than those presented in <cit.>.

As can be seen that the main contaminators to the high-redshift quasars are M, L, and T dwarfs, which agrees with the studies using traditional color-cuts, because M, L, and T dwarfs exhibit many absorption features in the optical and near-infrared bands, making them to have similar colors as the high-redshift quasars in those bands. To further improve the classification model, on one hand, we need to expand the training sample size, especially the high-redshift quasars, L and T dwarfs. On the other hand, considering the inclusion of photometric data from additional bands, such as i, and y bands, which are not available in Legacy Survey DR9 but will be available in future imaging surveys, may provide us with useful features to distinguish the high-redshift quasars from the contaminators. 

For quasars with z < 5, they all have ∼ 90% precision, indicating that each class has its own intrinsic characteristics that make them distinguishable. However, there is also significant contamination among them. On one hand, this is because they are in the post-reionization era of the universe, where the IGM is already mostly ionized, resulting in color similarities among these quasars. On the other hand, our rough setting of the redshift boundaries might lead to misclassification of quasars around the redshift boundary into other classes.

For M, L, and T dwarfs, the contamination among them is significant. This is because they all belong to subtypes of brown dwarfs, which are differentiated based on the spectra <cit.>, and photometric data lacks the precision of spectroscopic data. We find that a significant portion of L and T dwarfs are misclassified as other classes, which is due to the high similarity in the optical spectra and comparably small sample size. Moreover, because of the limitation of the sample size, the color space for L and T dwarfs are not fully spanned and not representative enough to separate them. The issues for both M dwarfs and quasars with z < 5 are highly mitigated. As for other type stars, although there is serious contamination among them, we have not find any significant impact on the predictions of the high-redshift quasars. 

We notice that the top-ranked feature of z-W2 in Table <ref> contradicts to the color-color diagram presented in Figure <ref>. Here we investigate this issue in detail. It is known that the color-color diagram presented in Figure <ref> is the projection of the high-dimensional color space. Then the projection on certain direction might diminish the separation between the high-redshift quasars and the contaminators. In order to fully resolve this issue, we construct a three-dimensional (3-D) color space from the top ranked features presented in Table <ref>. Figure <ref> illustrates the color diagram of high-redshift quasars and the main contaminant sources, MLT dwarfs, in the 3-D color space of z-W2, g-z, and r-grz. It can be observed that while there was significant overlap in the 2-D color space of z-W2 and r-grz, these sources exhibit clear separation in the 3-D color space. This demonstrates that the signal and the contaminators can be better distinguished in high-dimensional color space, which is not evident in the projected 2-D color space.
 





 §.§ Discussion of the imbalanced issue




We notice that the number of instances in each class in the training sample is not balanced, as shown in Table <ref>. Here we are discussing the impact of the imbalanced data sample on the random forest classification model. The ratio of the most numerous class (low-z quasars) to the least abundant class (high-z quasars) is as high as 700:1. 
Typically, classification models trained on the extremely imbalanced datasets tend to favor the majority class, making it easier for the minority classes to be incorrectly classified as the majority. This skews the model's classification performance for the minority class, more specifically, it results in a lower recall for the minority class, which also 
indirectly affects the precision of the minority class.


There are multiple approaches to deal with the imbalanced datasets. For example, the random forest classification algorithm offers a class_weight parameter to adjust the weights of each class based on the number of instances in the input dataset, with the goal of giving the classes that are less frequent a higher weight. This approach aims at balancing the dataset by ensuring that the minority classes are not overwhelmed by the majority classes. There are two modes (“balanced” and “balanced_subsample”) to adjust the class weights. The “balanced” mode uses the values of labels to automatically adjust weights inversely proportional to class frequencies in the input data. The “balanced_subsample” mode is the same as “balanced” except that weights are computed based on the bootstrap sample for every tree grown. 




The two methods mentioned above do not involve any reprocessing of the training samples but rather assigning different weights to different classes to achieve balanced datasets. There are methods that directly alter the number of instances for each class in the training sample to balance the dataset, and the methods can be further divided into two categories. One approach is to reduce the entire training sample according to the minority class, which is called “under-sampling". And the other approach is generating more training sample according to the majority class, which is called “over-sampling".






The commonly used under-sampling algorithms include: Balanced Random Forest <cit.>, Near Miss <cit.>, and One Sided Selection <cit.>. The main difference among the algorithms lies in the way they undersample the training sample. BRF draws a bootstrap sample from the minority class and sample with replacement the same number of samples from the majority class. NM uses the KNN algorithm to remove majority class samples near the minority class. OSS rejects the majority class samples around the minority class and near the decision boundary, and those that have minimal impact on the model. We adopt these three under-sampling methods to reprocess the training sample. 



The commonly used over-sampling algorithms include: Random Over Sampler (ROS), Synthetic Minority Oversampling Technique <cit.>, and Adaptive Synthetic Sampling <cit.>. The approaches to over sample the training sample are different for these algorithms. ROS randomly samples and duplicates the minority class with replacement. SMOTE synthesizes new samples between the nearest neighbors of minority class samples in the feature space, while ADASYN focuses on synthesizing new samples between the nearest neighbors of minority class samples misclassified by the KNN classifier. We adopt these three algorithms to reprocess the training sample. 


Table <ref> shows the performance for the RF model without resampling and the RF models on the balanced samples for the test set. 
For the class_weight balanced methods, it can be observed that setting the class_weight parameter has little impact on recall but reduces precision compared to the results on imbalanced samples. 
For the under-sampling balanced methods, the precision of both BRF and NM methods decreases significantly, while both recalls show a significant improvement. This is because the under-sampling methods reduces the number of the majority samples significantly, leading to the loss of important information from the contaminators and shrinking in the color space. This issue is particularly obvious for the BRF method which is based on random under-sampling. 
Surprisingly, we find that the OSS method returns comparably good precision and excellent recall. 
Despite such a small training set, we are able to obtain comparably good models with balanced dataset. This also reflects the fact that the high-redshift quasars possess inherent features such that they are separable from the contaminators in the high-dimensional color space. For the over-sampling balanced methods, the metrics obtained by these three methods are quite consistent.
Compared to the metrics on the imbalanced dataset, they show little change in recall but a decrease in precision. Although the performance of the over-sampling balanced methods is excellent, one issue about these methods is that the over-sampled instances  might even not exist or might be biased in the color space. 


In conclusion, most of the methods discussed for handling imbalanced datasets do not show significant differences in all four metrics (precision, recall, f1, AGF) compared to those on the imbalanced dataset. Additionally, these balanced methods often prioritize on enhancing the recall of the high-z class while sacrificing its precision. However, what we are more concerned about is the precision of the high-z class, which is more crucial for correctly classifying a vast amount of data with unknown label using our model. Therefore, we will still adopt the 11-class random forest classification model on our original training sample.









 §.§ Future Imaging Surveys


There are multiple complementary imaging surveys on the way, which will be the frontier for high-redshift quasars searching. On one hand, there will be more confirmed high-redshift quasars in the near future, expanding the size of training sample significantly. On the other hand, those deep imaging surveys not only provide much better images with high-quality, but also do they provide many more features for both high-redshift quasars searching and photo-z estimation. For example, the Chinese Space Station Telescope <cit.> will cover ∼ 17000 square degrees of the sky and has a wavelength coverage of 2600-10000 Å with broad-band filters of near-ultraviolet (NUV), u, g, r, i, z, y. The Roman Space Telescope 
<cit.> has a wide wavelength coverage of 0.48-2.3 microns with 8 broad-band filters, whose near-infrared images will be a lot better than the WISE W1, W2 images. The 5σ AB magnitude limits could reach ∼ 24[https://roman.gsfc.nasa.gov/science/WFI_technical.htmlromanWFItechnical.org/] <cit.> 
for filter F184 and F213 (central wavelength are 1.84 and 2.31 μm, respectively) even with only 55 seconds exposure. The Euclid Space Telescope <cit.> will deliver data over 15,000 square degrees of the sky with a wavelength range of 0.95-2.02 microns. The Large Synoptic Survey Telescope <cit.> will deliver the best images in 6 filters for over 18,000 square degrees sky coverage, and the 5σ point-source depth of the coadded maps can reach 27.5 AB magnitude for r-band <cit.>. The near- or mid-infrared photometric data from those wide-field imaging surveys will significantly improve the performance of high-redshift quasars searching.

Furthermore, those future imaging surveys could reach much deeper detection limit, which might resolve the missing values issue for some of the high-redshift quasars. Here we investigate the imputation method in which the missing values are replaced with the 5σ point-source depth 
from future imaging surveys (refer to Sec. <ref> for details). The evaluation scores (precision: 0.96, recall: 0.91, f1: 0.93, and AGF: 0.95) of the imputation method using future imaging survey limits on the test set are very close to those using the MICE imputation method (precision: 0.96, recall: 0.92, f1: 0.94, and AGF: 0.96). This indicates that the future imaging surveys are effective to identify the high-redshift quasars. When incorporating additional features other than the g, r, z, W1, W2 bands from future imaging surveys, one could imagine how effective our random forest classification model will be.





§ PHOTOMETRIC REDSHIFT ESTIMATION
 

Generally, spectroscopic redshift is challenging to obtain for a large number of high-redshift quasar candidates, then the photometric redshift could provide us a good estimation to study various distance-related aspects of the high-redshift quasars. The photometric redshift estimation algorithms can be divided into two categories: template fitting and machine learning methods. Template fitting involves establishing a relationship between photometric magnitudes or fluxes and spectroscopic redshifts using a series of spectral energy distribution (SED) templates. For example, EAZY <cit.> is an algorithm based on template fitting, in which the default template set, as well as the default functional forms of the priors are from semianalytical models. However, machine learning methods seek to find the relationship between photometric information and spectroscopic redshifts based on spectra-confirmed samples. Algorithms that are commonly used include k-nearest neighbors <cit.>, random forest regression <cit.>, and CatBoost <cit.>. Here we are using the random forest regression algorithm to establish a model to estimate the photometric redshift of the high-redshift quasars and introducing the evaluation metrics of the photometric redshift estimation.

The evaluation of photometric redshift is often based on the ratio of the number where the absolute value of Δ z = z_ spec - z_ photo is less than some threshold, e, to the total number, as shown in Equation <ref>. <cit.> conducted statistical analyses on the spectroscopic redshifts of previously verified quasars and the photometric redshifts derived from their optical and near-infrared photometric data. 


    ϕ_e = 1/N∑_z∈ Z |Δ z| < e

here Z is the set of spectroscopic redshifts, N is the total number of objects with spectroscopic redshift.  e is typically set to be 0.1, 0.2, 0.3 <cit.>. 
 
 The evaluation criteria mentioned earlier directly utilize the difference between the spectroscopic redshift and the estimated redshift, which is referred to as “non-normalized" and remains unaffected by the magnitude of the spectroscopic redshift. Another evaluation criterion is based on using the normalized difference scaled by the spectroscopic redshift, as shown in Equation <ref>. 



	
    Δ z_ norm = Δ z/1 + z_ spec ,
	
 here z_ spec is the spectroscopic redshift.

The fraction of entries where Δ z_ norm exceeds a certain threshold compared to the total number is commonly referred to as the outlier rate, η. As shown in Equation <ref>, we typically set this threshold to 0.1. A smaller outlier rate indicates that more of the predicted values are within the predetermined margin of error. This outlier rate evaluation criterion is more suitable for lower photometric redshift evaluations. The outlier rate is defined as:

    η_0.1 = 1/N∑_z∈ Z|Δ z_ norm| > 0.1 ,

where  represents an eigenfunction, defined as:

    x =
        1,true
    0,false
        ,





We utilize three widely employed machine learning algorithms (KNN, RF, CatBoost) to construct regression models across two distinct datasets, FeatureSet-mag and FeatureSet-flux, respectively. FeatureSet-mag comprises magnitudes of the g, r, z, W1, W2, W, and grz bands, alongside the apfluxes in the g, r, z, W1, and W2 bands. FeatureSet-flux encompasses fluxes of the g, r, z, W1, W2, W, and grz bands, in addition to the apfluxes in the g, r, z, W1, and W2 bands. There are only 602 high-redshift quasar samples (5 < z_ spec < 6.5), in order to enhance the accuracy of redshift predictions, we include the mid-z quasars with redshifts greater than 4.5 in the training sample. Similar to the classification models above, we divide the dataset into training, testing, and validation sets and use RandomizedSearchCV to find the best set of hyperparameters in the model's hyperparameter space.

We have ultimately developed six regression models based on the quasars  with redshifts ranging from 4.5 to 6.5. Table <ref> presents the R^2 (R-squared) and MSE (Mean Squared Error) scores of these six models. R^2 and MSE are common regression model evaluation metrics. R^2 assesses how well the model explains variable changes (0 to 1, closer to 1 for better fit), while MSE directly reflects prediction accuracy (smaller values for higher accuracy). Among the different training sets, the models trained on FeatureSet-mag outperform those trained on FeatureSet-flux, which also demonstrates the effectiveness of the MICE imputation method we use earlier. When considering different regression algorithms, the KNN model performs significantly worse compared to the other two models. Both the CatBoost and RF models exhibit very close and relatively good scores, with the RF model outperforming the CatBoost model.










Figure <ref> illustrates the distribution of spectroscopic redshifts and photometric redshifts for all quasars with redshifts from 4.5 to 6.5 for the three machine learning regression algorithms (KNN, RF, CatBoost) based on FeatureSet-mag, along with the distribution of Δ z_ norm. According to the traditional evaluation criterion ϕ_ e, the KNN, RF, and CatBoost models achieved ϕ_0.1 values of 42.30%, 66.98%, and 54.30%, respectively, representing the proportion of Δ z < 0.1 to the total count. The ϕ_0.2 values for these models are 65.62%, 90.64%, 86.72%. The ϕ_0.3 values are 78.38%, 96.60%, and 95.40%. Judging based on the normalized evaluation criterion η_0.1, with results of 2.89%, 0.68%, and 0.51% for the KNN, RF, and CatBoost models. In the lower panel of the three top subplots, the horizontal cyan lines represent Δ z_ norm = ± 0.1. It is clearly seen that there are almost no points falling outside the cyan parallel lines.

All three models show exceptional performance, although the results of the KNN model are slightly inferior to the other two models. In conclusion, it is evident that the RF regression model predicts redshifts with higher precision and accuracy. Therefore, we adopt the random forest regression model to predict the photometric redshifts for the high-redshift quasar candidates which will be discussed in the following. 









§ HIGH-Z CANDIDATES AND VERIFICATION
 




 §.§ Selection criteria



Based on the discussions above, we decide to adopt the 11-class random forest classification model on FeatureSet-C and FeatureSet-D. These models are named the “mag model" and the “flux model", respectively. Subsequently, we apply these two models separately to the entire Legacy Survey DR9 dataset to obtain the high-redshift quasar candidates. Before that, we have applied several selection criteria to reduce the size of the catalog data from LS DR9 (more than 1 billion entries) without missing too many signals. The selection criteria followed by the reasoning are: 

    
  * dered_mag_g,r,z,W1,W2 is not null if the features include colors.
    
  * brick_primary=1 and maskbits!=[1,10,12,13].
    
  * The type is  `PSF'.
    
  * snr_z>5, snr_W1>3 and snr_W2>2.
    
  * dered_mag_z > 15 and dered_mag_z<21.5.

1) We apply the first criterion to remove sources with missing values because the trained random forest model does not allow missing values, and imputation is challenging due to the large data size and unknown object labels. 2) The second criterion is used to ensure that the quality of the source is minimally affected by uncontrollable factors (such as being on the boundary, near bright sources, etc). 3) Setting the source type to Point Spread Function (PSF) helps eliminate contamination from extended sources like galaxies. 4) The fourth criterion sets the signal-to-noise ratio constraints, as referenced in <cit.>, helping reduce the catalog size without missing too many signals. 5) The last criterion sets the magnitude constraints to help reduce either too bright or too faint sources in the z-band. 

Although setting the source type to PSF can help us eliminate extended sources significantly, it cannot exclude all galaxies, such as the compact early-type galaxies at intermediate or higher redshift. These galaxies, due to their comparably high redshift and compactness, exhibit the morphology of PSF and have a redder color. We will investigate whether those distant compact galaxies are contaminants for the high-redshift quasars by training a binary RF classifier. For the training sample, we first obtain the galaxy sample from the SDSS DR17 Database[https://www.sdss4.org/dr17/sdss.org/dr17] <cit.> with spectral type of “Galaxy" and with z ≥ 0.2 (2,207,611 objects in total). Then, we cross-match the galaxy sample to the LS DR9 database  with the above selection criteria to obtain the sample of distant compact galaxies (19,374 objects in total). The sample of high-redshift quasars is the same as the one presented in Table <ref>. 

Similar to the 11-class classification model, we also obtain two parallel feature sets (FeatureSet-C and FeatureSet-D) with the best hyper-parameters determined by RandomizedSearchCV. The evaluation scores (precision, recall, f1, AGF) of the “high-z" class for the test set are 0.98, 0.99, 0.98, 0.99, respectively, for the “binary mag model", and the results for the “binary flux model" are consistent with those of the “binary mag model". The corresponding metrics (precision, recall, f1, AGF) of the distant compact galaxies for the test set are 1.00, 1.00, 1.00, 0.99, respectively, for both the “binary mag model" and the “binary flux model". These evaluation scores indicate that the distant compact galaxies almost do not contaminate the high-redshift quasars. Moreover, we find that the apflux ratios are indeed very effective to distinguish the high-redshift quasars from the distant compact galaxies.

Furthermore, we apply the trained 11-class classification models to the distant compact galaxy sample, and we found that only 4 galaxies (0.021%) are predicted as the “high-z" class by the “mag model", and only 1 galaxy (0.005%) is predicted as the “high-z" class by the “flux model". In conclusion, the contamination of the distant compact galaxies to the signal of high-redshift quasars is  minimal and negligible. 



 §.§ High-z candidates


After applying the above selection criteria, we have obtained a total of 140 million sources from LS DR9, which is passed to the two trained random forest classification models to obtain two parallel sets of high-redshift quasar candidates. The random forest classification algorithm provides the predicted probability of each category for each source object, we take the class with the highest predicted probability as the predicted class for each source object. The number of objects for each predicted class using the “mag model" is shown in the Table <ref>. The “flux model" returns similar results, with 420,208 mid-z quasar candidates and 568,188 high-z quasar candidates, respectively. And there are 216,949 overlapping high-redshift quasar candidates which are both identified by the “mag model" and the “flux model". 






To obtain more reliable high-redshift quasar candidates, further screening is necessary. The random forest classification model provides us with the predicted probabilities for each source object. By setting reasonable thresholds on these probabilities, we can select much more reliable high-redshift quasar candidates. Figure <ref> displays the distribution of photometric redshifts for high-redshift quasar candidates and known quasars with the probability of being classified as the high-z class. The orange and cyan dots depict our high-redshift quasar candidates, which are the overlapping candidates obtained from both the “mag model" and the “flux model" for the entire LS DR9 footprint. The red and green dots represent spectra-confirmed high-redshift quasars. The photometric redshift for the high-z class is provided by the random forest regression model trained on FeatureSet-mag.




We try first to set a probability threshold such that the precision of the known high-redshift quasars in our predicted results can reach 100%. To do that, we identify the highest predicted probability among the known high-redshift quasars that are incorrectly predicted as non-high-redshift quasars, and we call it as the first probability threshold, and the values are p_ thre1,mag = 0.41 and p_ thre1,flux = 0.40 for the “mag model" and the “flux model", respectively. They are represented by the lower blue and black horizontal lines in Figure <ref>, with the green dots representing known high-redshift quasars that are misclassified to other categories. This threshold helps us remove approximately one-third of the candidates.

The second probability threshold (p_ thre2) is related to the objects which are not “high-z" class but predicted as the “high-z" class. Among the objects predicted as the “high-z" class, we identify the highest probability among those that is incorrectly predicted as the “high-z" class. They are p_ thre2,mag = 0.83 for the “mag model" and p_ thre2,flux = 0.75 for the “flux model", and are represented by the upper blue and black dashed horizontal lines in Figure <ref>. The blue points in Figure <ref> represent the “mid-z" quasars that have been incorrectly predicted as the “high-z" class. As has been discussed earlier, the confusion is not a concern because their redshifts are around the boundary. This means that the precision of our trained random forest regression model is actually higher than what is currently being observed.

At this stage, we have obtained a relatively pure candidate set. To further ensure the plausibility of the high-redshift quasar candidates, we also quote the results for higher probability thresholds of 0.90 and 0.95. The number left after each of the probability threshold for the entire LS DR9 footprint is presented in Table <ref>. For the “mag model", we have p_ quasar≥ p_ thre1,mag: 198,339 candidates, p_ quasar≥ p_ thre2,mag: 2,984 candidates, p_ quasar≥ 0.90: 476 candidates, p_ quasar≥ 0.95: 32 candidates. For the “flux model", we have p_ quasar≥ p_ thre1,flux: 473,401 candidates, p_ quasar≥ p_ thre2, flux: 69,736 candidates, p_ quasar≥ 0.90: 11,012, p_ quasar≥ 0.95: 740. And for the overlapping candidates, we adopt the probability threshold derived from the “mag model", and we have p_ quasar≥ p_ thre1: 165,734 candidates, p_ quasar≥ p_ thre2: 2,984 candidates, p_ quasar≥ 0.90: 476, p_ quasar≥ 0.95: 32. Our catalog includes information of all candidates, and Table <ref> describes the columns of the data file, which is available in its entirety in machine-readable form. We present the catalogs obtained from the “mag model", the “flux model" and the overlapping results.





We notice that there is a big difference between the “mag candidate" and the “flux candidate" in Table <ref>. This difference becomes more pronounced as the probability threshold increases. Remember that the big advantage of the “flux model" (FeatureSet-D) is that there is no missing value issue. We suspect that the big difference between the “mag candidate" and the “flux candidate" is due to the missing value issue in the LS catalog. To fully understand this issue, we present the probability distributions of the high-redshift quasar candidates from the two models in Figure <ref>. The green area corresponds to the probability distribution of the candidates for the “mag model". The red area represents the probability distribution of the candidates for the “flux model".

It is clearly seen that there are many more candidates of the “flux model" than the “mag model", especially at the high probability region. It is obvious that the objects with missing values tend to have larger probability to be high-redshift quasars due to the color drop-out, which can be inferred from the training sample (more than half of the current known high-redshift quasars have missing value). If we remove the objects with missing values for the “flux model", which is represented by the sky blue area. It is obvious that the majority candidates of the “flux model" at high probability region have missing values, resulting in the big difference between the two models as shown in Table <ref>. The phenomenon observed in the high-redshift quasar candidates is consistent with that of the training sample.  










Furthermore, by calculating the completeness of known high-redshift quasars in the test set at different thresholds, we provide a reference for the completeness of our high-redshift quasar candidates. The calculation of completeness is similar to recall. At different thresholds, if the probability of the high-z class obtained from the random forest is higher than the threshold, it is considered a correctly classified high-redshift quasar; otherwise, it is deemed a misclassified high-redshift quasar. The completeness of the high-z class is the ratio of correctly classified high-redshift quasars to the total number of high-redshift quasars in the test set. For the “mag model", the completeness of the known high-redshift quasars in the test set at different thresholds is as follows: 
p_ quasar≥ p_ thre1,mag: 82.20%; p_ quasar≥ p_ thre2,mag: 33.05%; p_ quasar≥ 0.90: 18.64%; p_ quasar≥ 0.95: 5.09%. Similarly, for the “flux model", the completeness is as follows: p_ quasar≥ p_ thre1,flux: 79.49%; p_ quasar≥ p_ thre2,flux: 41.88%; p_ quasar≥ 0.90: 15.39%; p_ quasar≥ 0.95: 0.00%.







 §.§ Verification using MUSE


The Multi Unit Spectroscopic Explorer <cit.>, an integral field unit (IFU) mounted on the Very Large Telescope (VLT-Yepun, UT4), has a field of view (FOV) of 1 arcmin^2. The MUSE instrument provides high throughput (35% end-to-end, including the telescope, at 7000 Å), moderate spectral resolution (R ≃ 3000 at ∼7000 Å), full optical coverage 4650 - 9300 Å) spectroscopy
at a spatial sampling scale of 0.2 arcsec across a 1 × 1 arcmin^2 FOV. The wide wavelength coverage and moderate spectral resolution enable us to identify high-redshift quasars up to redshift of ∼ 6.6 and suitable for our purpose of high-redshift quasar candidates verification.

So far there are nearly 20,000 publicly available MUSE datacubes, of which the sky coverage is approximately 6 square degrees. Among the entire high-redshift candidates, 21 of them are located within the MUSE field. After spectroscopic verification, 11 of these candidates are high-redshift quasars included in our training set, 3 are known high-redshift quasars but not present in our training sample. 
5 of these candidates are M dwarfs, while the rest two have spectra that are much noisier and are difficult to identify emission/absorption lines, making it hard to determine the type of those sources. By estimating their redshifts using Lyα, N V or C IV emission lines, we find that these three “new" high-redshift quasars precisely align with the redshift boundaries we set for high-redshift quasars. Two of them are having redshifts of z ∼ 5 , and one is near the boundary of 6.5. Because they fall on the redshift boundaries,  we initially miss them when collecting the high-redshift quasars training sample. Since they are not part of our training sample, they provide us a valuable opportunity to validate our model.

Overall, 14 out of 21 high-redshift quasar candidates from our random forest model has been confirmed to be true high-redshift quasars, reaching a success rate of 66.7%. If we only focus on the high-redshift quasars not present in the training sample, the success rate is 30% (3 “new" high-redshift quasars out of 10 candidates). 
Figure <ref> shows the MUSE spectra of these three “new" high-redshift quasar obtained from our random forest model. The details of these three “new" high-redshift quasars are the following:

J014132.4-542749.9. <cit.> discovered this quasar by cross-matching the first data release of the Dark Energy Survey with the Sidney University Molonglo Survey radio catalog at 0.843 GHz. The z band magnitude of this quasar is 21.07, 
and the photometric redshift provided by our trained regression model is 5.087. <cit.> utilized Gaussian profiles based on the Lyα, O VI, N V, and C IV emissions to estimate the redshift of this quasar, with the value of 5.000 ± 0.002.

J103418.65+203300.2. <cit.> discovered this quasar in the SDSS Quasar Catalog. Its magnitude in the z band is 19.61. 
The photometric redshift provided by our trained regression model is 5.214. Based on the emission lines of Hα, Hβ, Mg II, C III, C IV, and Lyα,  the redshift of this quasar is estimated to be 5.0150 ± 0.0005 <cit.> .

J022426.54–471129.4. <cit.> discovered this quasar using the European Southern Observatory New Technology Telescope and Gemini South telescopes. The z band magnitude of this quasar is 20.23, 
and the photometric redshift provided by our trained regression model is 6.139. The significant difference between the photometric and spectroscopic redshifts is due to the scarcity of samples with redshift around 6.5 in our training set. <cit.> calculated the redshift of this quasar by fitting its spectrum with a known quasar model, obtaining a value of 6.50 ± 0.01.






 §.§ Verification using DESI

DESI will collect millions of qusars in its 5-year operation <cit.>, of which thousands of high-redshift quasars will be identified. Currently there are 306 objects at z > 5 in the DESI-EDR catalog <cit.>. We then apply our trained RF classification model to these 306 objects at z > 5 in the DESI-EDR catalog to validate and evaluate our model. Before that, we first inspect the spectra of those objects visually. After spectra inspection, we find that only 22 of them are true high-redshift quasars at z > 5, 12 of which are in our training sample and 10 of which are new high-redshift quasars. The rest 284 objects are not high-redshift quasars after spectra inspection, which might be due to either the incorrect classification or the incorrect redshift determination by the DESI team.

This ensemble of objects at z >5 from DESI-EDR provides a perfect labratory to test and validate our random forest classification model. For the “mag model", 188 out of 306 sources have no missing values, with only 6 out of the 22 true high-redshift quasars found among these 188 sources. 4 out of 188 objects are predicted to be high-redshift quasars, and these 4 objects are among the 6 true high-redshift quasars. The rest 2 true high-redshift quasars are not correctly classified: one is classified as a M dwarf, and the other is classified as a mid-z quasar. This implies that the precision is 1.00 and the recall is 0.67 for the “mag model". Using the MICE imputation to fill in the missing values for the 118 sources, we find that  17 out of 118 sources are predicted to be high-redshift quasars. 16 out of the 17 candidates are true high-redshift quasars, and the type of the rest 1 source is suspicious because it is too faint without distinct emission or absorption lines. Overall, 21 objects are predicted to be high-redshift quasars, 20 of which are true high-redshift quasars. Therefore, the precision is 0.95 and the recall is 0.91 for the “mag model".

For the “flux model", all 306 sources have well-defined features. Among these sources, 21 objects are predicted to be high-redshift quasars, 20 of which are true high-redshift quasars. The rest 1 object is the same as the above one whose type is suspicious. The precision is 0.95, and the recall is 0.91 for the “flux model". As discussed earlier, there exist a confusion between “mid-z" and “high-z" quasars because of the set of arbitrary redshift boundary, which should not be a concern. In this sample, one “high-z" quasar with true redshift of 5.036 is predicted to be a “mid-z" quasar, and we should not consider this as an error. This means that the precision and recall will be slightly higher, with a precision and recall of 0.96, 0.96 for both the “mag model" and the “flux model". Even if we only focus on the 10 new high-redshift quasars, 9 of which are successfully recovered. It turns out that our random forest classification model can achieve much higher success rate than the DESI target selection. The consistency between the “mag model" and the “flux model" further demonstrates the robustness and unbiasness of the imputation method we adopt.

Although the estimation of the precision and recall for the high-redshift quasars using the DESI-EDR database is biased because these objects are not randomly picked from our high-redshift quasar candidates and they are pre-selected by the DESI team, they can still shed a light on the robustness of our machine learning model. On the other hand, the DESI target selection of quasars is mainly for quasars at 0.9 < z < 2.1 and seldom of these quasars are at z > 2.1 <cit.>. Moreover, the DESI target selection of high-redshift quasars is still based on the traditional color-cut method that has been adopted in previous high-redshift quasar surveys <cit.>. In this sense, these 306 objects extracted from  DESI-EDR database can be viewed as randomly picked objects from our high-redshift quasar candidates. Therefore, if our estimation of the precision and recall for the high-redshift quasars from DESI-EDR database is biased, they are least biased. The high success rate of our high-redshift quasar candidates using the objects at z>5 from DESI-EDR further validates our random forest classification model. Therefore, if we pick the high-redshift quasar candidates with high probability (p_ quasar > 0.9 for instance) from our classification model for follow-up spectroscopy observations, we will discover many more high-redshift quasars than expected.




§ SUMMARY
 

In this work, we obtain a representative training sample which are all spectra-confirmed and it consists of 588 high-redshift quasars and millions of contaminators in the DESI Imaging Legacy Survey DR9 footprint. We then obtain the photometric data for the training sample from Legacy Survey DR9 and WISE survey and constructed various features to train a machine learning model. In addition to the photometric measurements in the g, r, z, W1, and W2 bands, we also construct the combined photometric measurements of grz and W. We  implement the fluxes at different aperture radius provided by the LS DR9, the apfluxes are very effective in distinguishing high-redshift quasars from the extended sources like distant compact galaxies. And we demonstrate that the contamination from distant compact galaxies can be negligible via a binary classification for the first time. In addition, we address the issue of missing values in the known high-redshift quasar dataset. Among the various imputation methods, we find that the predictive results of the model trained with the MICE method are best consistent with the complete dataset.

We compare the performance of several different commonly used machine learning classification algorithms, including KNN, decision tree, random forest, LGBM, GaussianNB. The random forest classification model emerged as the most optimal algorithm. To further enhance the performance of the random forest classification model, we perform feature selections on the training sample. Additionally, we discuss whether the number of classes in the training sample affects the performance of the classification model. Considering the intrinsic characteristics of each class in the training sample, we decide to adopt the 11-class classification model as our final model. We discuss the issue of imbalance sample and find that different algorithms handling the imbalanced data sample do not show significant differences in all metrics compared to those with the imbalanced dataset.

The performance of the random forest with 11-class can reach a high precision of 96.43%, a high recall of 91.53% for the high-redshift quasars for the test set, both of which are significantly higher than the previous studies. In addition to the classification model, we also train a regression model on a dataset of quasars with redshift between 4.5 and 6.5 to predict the photometric redshift of the high-redshift candidates. Among the different machine learning regression algorithms (KNN, random forest, CatBoost), we find that the random forest regression model has the best performance, with 99.3% data points within the range of Δ z/1 + z_ spec < 0.1 and MSE ≤ 0.025. 

We apply several selection criteria to reduce the catalog size of the entire Legacy Survey DR9 catalog. Applying the final 11-class models using random forest algorithm on FeatureSet-C (FeatureSet-D) to the 140 millions sources from LS DR9, we obtain 272,424 (568,188) high-redshift candidates. And there are 216,949 candidates which are classified as high-redshift quasars by both the “mag model" and the “flux model". We further narrow down the high-redhisft quasar candidates by setting two cutoff probabilities from the random forest model: the maximum probability of a known high-redshift quasar being misclassified as another class (40.82%, 40.03% for the “mag model" and “flux model", respectively) and the maximum probability of a contaminator being misclassified as a high-redshift quasar (82.49%, 75.23% for the “mag model" and “flux model", respectively). By applying these cutoff probabilities to the common high-redhisft quasar candidates identified by both the “mag model" and the “flux model", we obtain 165,734 and 2,984 high-redhisft quasar candidates, respectively. There are 476 candidates with probability greater than 90%, which could be the targets of the highest priority for future spectroscopy followup.


Using MUSE spectra to do the spectra verification of the high-redshift quasar candidates, we find that there are 21 high-redshift quasar candidates in the total ∼ 20000 MUSE observations footprint. After inspection, we find that 11 of these candidates are high-redshift quasars already included in our training sample, 3 are known high-redshift quasars but not present in our training sample. Using DESI-EDR spectra, we confirm that 21 out of 22 true high-redshift quasars with correct redshift can be successfully identified for the “mag model" with missing values imputed and 21 out of 22 true high-redshift quasars can also be successfully classified by the “flux model", reaching much higher success rate than the DESI target selection. 


Our current model still has  room for improvement. The future DESI survey will discover thousands of new high-redshift quasars, which will significantly expand the training sample size. The current photometric measurements available are very limited, therefore introducing new photometric measurements in the near- or mid-infrared could enhance the model performance significantly. The Legacy Survey DR10 currently includes the i band for ∼ half sky coverage of DR9, and the future data release including the y band is on the way as well. We find that the inclusion of i band photometric measurement improves the classification performance significantly, which demonstrates that more photometric measurements are critical to improving the machine learning performance. The future imaging surveys such as CSST, RST, EST and LSST, would deliver much deeper images with many more photometric measurements such as the y, J, H and K bands, which are beyond the current surveys and would be the frontier for high-redshift quasars searching with redshift even up to 7 or 8. 



§ ACKNOWLEDGMENTS


GY and HZ acknowledge financial support from the start-up funding of the Huazhong University of Science and Technology and the National Science Foundation of China grant (No. 12303007). QW acknowledges financial support from the National Science Foundation of China grant (No. 12233007). The authors thank the referees for comments that helped us improve the manuscript. This research was supported by the Munich Institute for Astro-, Particle and BioPhysics (MIAPbP) which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy EXC-2094 390783311.

The Legacy Surveys consist of three individual and complementary projects: the Dark Energy Camera Legacy Survey (DECaLS; Proposal ID #2014B-0404; PIs: David Schlegel and Arjun Dey), the Beijing-Arizona Sky Survey (BASS; NOAO Prop. ID #2015A-0801; PIs: Zhou Xu and Xiaohui Fan), and the Mayall z-band Legacy Survey (MzLS; Prop. ID #2016A-0453; PI: Arjun Dey). DECaLS, BASS and MzLS together include data obtained, respectively, at the Blanco telescope, Cerro Tololo Inter-American Observatory, NSF's NOIRLab; the Bok telescope, Steward Observatory, University of Arizona; and the Mayall telescope, Kitt Peak National Observatory, NOIRLab. Pipeline processing and analyses of the data are supported by NOIRLab and the Lawrence Berkeley National Laboratory (LBNL). The Legacy Surveys project is honored to be permitted to conduct astronomical research on Iolkam Du'ag (Kitt Peak), a mountain with particular significance to the Tohono O'odham Nation.

NOIRLab is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. LBNL is managed by the Regents of the University of California under contract to the U.S. Department of Energy.

This project used data obtained with the Dark Energy Camera (DECam), which was constructed by the Dark Energy Survey (DES) collaboration. Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundacao Carlos Chagas Filho de Amparo, Financiadora de Estudos e Projetos, Fundacao Carlos Chagas Filho de Amparo a Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Cientifico e Tecnologico and the Ministerio da Ciencia, Tecnologia e Inovacao, the Deutsche Forschungsgemeinschaft and the Collaborating Institutions in the Dark Energy Survey. The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energeticas, Medioambientales y Tecnologicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenossische Technische Hochschule (ETH) Zurich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciencies de l'Espai (IEEC/CSIC), the Institut de Fisica d’Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig Maximilians Universitat Munchen and the associated Excellence Cluster Universe, the University of Michigan, NSF’s NOIRLab, the University of Nottingham, the Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, and Texas A&M University.

BASS is a key project of the Telescope Access Program (TAP), which has been funded by the National Astronomical Observatories of China, the Chinese Academy of Sciences (the Strategic Priority Research Program “The Emergence of Cosmological Structures” Grant #XDB09000000), and the Special Fund for Astronomy from the Ministry of Finance. The BASS is also supported by the External Cooperation Program of Chinese Academy of Sciences (Grant #4A11KYSB20160057), and Chinese National Natural Science Foundation (Grant #12120101003, #11433005).

The Legacy Survey team makes use of data products from the Near-Earth Object Wide-field Infrared Survey Explorer (NEOWISE), which is a project of the Jet Propulsion Laboratory/California Institute of Technology. NEOWISE is funded by the National Aeronautics and Space Administration.

The Legacy Surveys imaging of the DESI footprint is supported by the Director, Office of Science, Office of High Energy Physics of the U.S. Department of Energy under Contract No. DE-AC02-05CH1123, by the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility under the same contract; and by the U.S. National Science Foundation, Division of Astronomical Sciences under Contract No. AST-0950945 to NOAO.

This research used data obtained with the Dark Energy Spectroscopic Instrument (DESI). DESI construction and operations is managed by the Lawrence Berkeley National Laboratory. This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of High-Energy Physics, under Contract No. DE–AC02–05CH11231, and by the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility under the same contract. Additional support for DESI was provided by the U.S. National Science Foundation (NSF), Division of Astronomical Sciences under Contract No. AST-0950945 to the NSF’s National Optical-Infrared Astronomy Research Laboratory; the Science and Technology Facilities Council of the United Kingdom; the Gordon and Betty Moore Foundation; the Heising-Simons Foundation; the French Alternative Energies and Atomic Energy Commission (CEA); the National Council of Science and Technology of Mexico (CONACYT); the Ministry of Science and Innovation of Spain (MICINN), and by the DESI Member Institutions: www.desi.lbl.gov/collaborating-institutions. The DESI collaboration is honored to be permitted to conduct scientific research on Iolkam Du’ag (Kitt Peak), a mountain with particular significance to the Tohono O’odham Nation. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the U.S. National Science Foundation, the U.S. Department of Energy, or any of the listed funding agencies.

Based on observations collected at the European Southern Observatory under ESO programme(s) 109.238W.003, 105.208F.001, and 0104.A-0812(A), and/or data obtained from the ESO Science Archive Facility with DOI(s) under https://doi.org/10.18727/archive/41.






