Random Forests as a Viable Method to Select and Discover High-redshift Quasars

The data we are mining for quasars is a cross-match between the publicly available Pan-STARRS DR1 (described in Chambers et al. 2016) and ALLWISE (Cutri et al. 2021) catalogs. The ALLWISE survey is a release of the aggregated data from WISE and its extended mission NEOWISE (Mainzer et al. 2011) up until 2013. From Pan-STARRS, we use the five stacked PSF magnitudes (${g}_{{\rm{P}}{\rm{S}}{\rm{F}}},\,{r}_{{\rm{P}}{\rm{S}}{\rm{F}}},\,{i}_{{\rm{P}}{\rm{S}}{\rm{F}}},\,{z}_{{\rm{P}}{\rm{S}}{\rm{F}}}$, and y_PSF), stacked aperture magnitude in the z band (z_APERTURE), mean position, and the objectinfoFlag. From WISE, we use the 3.4 and 4.6 μm broadband magnitudes (${\rm{W}}1,{\rm{W}}2$), their signal-to-noise ratio (${\rm{W}}{1}_{{\rm{s}}/{\rm{n}}},{\rm{W}}{2}_{{\rm{s}}/{\rm{n}}}$), position, the active deblending flag (na)- and the number of PSF components used for the PSF fitting (nb).

We use the Python framework Large Survey Database (Juric 2012) to cross-match the two catalogs, applying the following selection criteria:

$$\begin{eqnarray}&&14\lt {z}_{\mathrm{PSF}}\leqslant 20.5\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{{\rm{y}}}_{\mathrm{PSF}}\ \mathrm{is}\ \mathrm{not}\ \mathrm{None}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&-0.3\leqslant {z}_{\mathrm{PSF}}-{z}_{\mathrm{APERTURE}}\leqslant 0.3\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\mathrm{galactic}\ \mathrm{latitude}\gt 20^\circ \ \mathrm{or}\lt -20^\circ \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\mathrm{objinfoFlag}\ \mathrm{has}\ {\mathtt{GOOD}}\ \mathrm{and}\ {\mathtt{GOOD}}\_{\mathtt{STACK}}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&2\buildrel{\prime\prime}\over{.} 0\ \mathrm{match}\ \mathrm{in}\ \mathrm{ALLWISE}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{\rm{W}}{1}_{{\rm{s}}/{\rm{n}}}\geqslant 5\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{\rm{W}}{2}_{{\rm{s}}/{\rm{n}}}\geqslant 3\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{\mathtt{na}}=0,{\mathtt{nb}}=1.\end{eqnarray} \tag{ 9 }$$

The resulting catalog has around 72 million objects. The z band is used to select the brightness range from 14 to 20.5 in magnitude. Since the brightest quasar at z ≥ 4.7 in our training data (see Section 2.2) has a z-band magnitude of z_PSF = 17.3, there is only a remote chance to miss quasar lenses by adopting a bright limit for our selection. We also choose a faint limiting magnitude on the z band that is well above the detection limit, to ensure the reliability of the photometry. The 5σ detection limit for the Pan-STARRS survey in the stacked z band for point sources is 22.3 mag in AB (see Table 11 of Chambers et al. 2016). They also showed that, in the z band, the 98% source completeness limit is fainter than 20.5 mag on most of the sky, especially away from the galactic plane (see Figure 17 of Chambers et al. 2016). The criteria on the z band automatically remove all objects with a missing z-band detection. We further require the y band not to be "None." However, the other bands of Pan-STARRS (g, r, and i) can be missing because we expect the targeted high-redshift quasars to have very little flux in these bluer bands.

We use the difference of the PSF and aperture magnitude in the z band to actively exclude sources with extended morphologies from our selection. The cutoff of ${z}_{\mathrm{PSF}}-{z}_{\mathrm{APERTURE}}=\pm 0.3$ is informed by Figure 3 in Bañados et al. (2016), where the magnitude difference is compared for stars, quasars, and galaxies. This cut is designed to effectively remove galaxies from our selection; however, it may also reduce our sensitivity to lensed quasars. In our final candidate selection (see Section 5), there are only a few remaining galaxies that were removed during visual selection, so this approach is sufficient for our purpose.

We furthermore restrict our selection to Galactic latitudes of $| b| \geqslant 20^\circ$, where the contamination by galactic stars decreases significantly. We require the photometry to fulfill the GOOD and GOOD_STACK flags in the objectinfoFlag from Pan-STARRS. These are quality flags provided by Pan-STARRS to indicate that the object has a good-quality measurement in the data and a good-quality object in the stack (>1 good stack measurement). We matched our objects with ALLWISE with a radius of 20 using only the closest match. We require that W1 and W2 are detected with respective signal-to-noise ratios of 5 and 3. We exclude obviously blended sources via the active deblending flag (na) and the number of PSF components used for the PSF fitting (nb) from WISE. These WISE flags ensure more reliable photometry, but we note that they reduce our sensitivity to lensed quasars and may remove some quasars with close-by sources.

To determine which survey limits our quasar selection, we consider the set of additional high-redshift quasars, discussed in the next section. Of a total of 1001 quasars, 936 have a Pan-STARRS match, and 657 fulfill conditions (1–3) mostly limited by the brightness cut in the z band.

We contrast this with 710 objects that have an ALLWISE match and fulfill ALLWISE photometry conditions (6, 7), while 647 additionally fulfill condition (8). Both the Pan-STARRS and ALLWISE photometry requirements remove a similar fraction of known quasars, and the remaining objects have a large overlap: 565 objects fulfill conditions (1–3 and 6–8). This shows that our requirements on both surveys are well-balanced for our targeted class. To use fainter objects in the Pan-STARRS data, we would also require deeper infrared data. We note that, in Table 2, we only list additional quasars that are not already in the other set.

Random forests are a supervised machine-learning algorithm and therefore heavily rely on representative training sets. It is essential to assemble a training set consisting of spectroscopically identified objects representing the wide range of different objects in our catalog data. For a reliable selection of high-redshift quasars at z = 4.7–6, we need to make sure to construct a representative training set. This means that we need to include all potential contaminants that populate the same color space, like M, L, and T dwarfs. We do not include galaxies in our training set, because we already removed extended sources from our data set (see selection criterion 3). The training classes used in our random forest selection are listed in Table 1.

We do not include O- and B-type stars, as they are far from high-redshift quasars in color space. These classes are irrelevant since they likely get assigned the label of the most similar star class, but are not confused with our targeted quasars. We exclude the Y-type brown dwarfs since we do not have many objects of that class and they are also not relevant contaminants for quasars with z ∼ 5–6. Similarly, there are classes of objects that are underrepresented in our training set, like low-redshift BAL quasars, which are known contaminants for high-redshift quasars.

We built our training set with the spectroscopic training set from Schindler et al. (2019). It is based on the spectroscopically confirmed quasars from the Sloan Digital Sky Survey (SDSS) DR7 and DR12 quasar catalogs as well as the spectroscopically confirmed stars from SDSS DR 13. The SDSS data were matched to Pan-STARRS within 398, using only the closest match. For a full discussion of the data processing, see the referenced paper. We take only the Pan-STARRS position and classification as well as the redshift for the quasars, and reprocess the rest of the data for internal consistency. We expand this training set with more objects in the relevant color space region. To increase the number of red and brown dwarfs, the Dwarf Archive ¹⁰ is used. We match their positions to Pan-STARRS within 20. This should be sufficient for our purposes, but we acknowledge that a more careful cross-match considering proper motion would increase the number of dwarfs further. To supplement our training set with additional high-z quasars, we added a comprehensive list of quasars known as of mid-2018. Again, we cross-match the position to Pan-STARRS within 20. The major sources are Wang et al. (2016), Bañados et al. (2016), Jiang et al. (2016), McGreer et al. (2013), Matsuoka et al. (2018a), Yang et al. (2017), and the preliminary results of Yang et al. (2019) as of mid-2018.

For each of our three data sets, we download the Pan-STARRS data from MAST using the CasJobs ¹¹ interface. We cross-match with ALLWISE with a radius of 20 using IRSA ¹² and only using the closest match. To include as many training objects as possible, we omit some of our photometric selection criteria used for the full catalog data. The classification information added outweighs the downside of them not being fully representative of the catalog data we collected. We require conditions (2, 6, and 7) as well as a detection in the Pan-STARRS z and WISE W2 bands, which means a value entry in the catalog. Finally, we remove duplicates between our three data sets based on the Pan-STARRS ObjID.

In total, the Schindler et al. (2019) data set gives us 259,240 stars and 164,318 quasars that fulfill our photometric requirements. In addition to the SDSS DR7 and DR12 quasars, we add a total of 474 additional quasars from high-redshift surveys. We also add an additional 470 dwarfs from the Dwarf Archive. Table 2 lists the different sources and classifications for the training set, along with the number of usable objects they add. In Figure 1, we show a plot of the dust-corrected z-band magnitude from Pan-STARRS versus redshift for all training quasars. Outliers beyond magnitude 25 are cut off. At the high-redshift end, the additional quasars from recent surveys significantly extend the training set. We note that there is a visible underdensity of known quasars around redshift 5.5. At this redshift, quasars are challenging to find, as their colors in common bands are very similar to those of M stars (Yang et al.2019).

Zoom In Zoom Out Reset image size

In Figure 2, we show color–color plots for the objects in our training set. This highlights that the Pan-STARRS and WISE bands contain information that will allow us to differentiate the classes listed in Table 1. We note that this visualization emphasizes how the average color information of quasars differs from stars and evolves with redshift, but hides the complexity of applying this to real data. The 5% of stars with larger scatter than the ellipses shown far outnumber the high-z quasars, and the quasars themselves are also scattered. To reliably classify new objects, we need the full high-dimensional color information.

Zoom In Zoom Out Reset image size

We correct both our catalog data and training set for galactic extinction based on the dust map of Schlegel et al. (1998), using the sfdmap ¹³ Python package. The Vega magnitudes in the ALLWISE catalog are converted to AB magnitudes using the constants ${\rm{W}}{1}_{\mathrm{AB}}-{\rm{W}}{1}_{\mathrm{Vega}}=2.699$ and ${\rm{W}}{2}_{\mathrm{AB}}-{\rm{W}}{2}_{\mathrm{Vega}}=3.339$ (Sec IV 4h, Cutri et al. 2012). All magnitudes are then converted to flux density in Jansky units. Our catalog restrictions allow objects with nondetections in the g, r, and i bands to be considered. However, the random forest method cannot handle null values. We work around this by replacing the missing values with a fixed value that is fainter than the detection limit of the catalog. This way, the resulting flux density ratios will be close to the true values. We choose 1e−10 Jy or a magnitude of 33.90 in the AB system as the replacement value. We replace all missing g, r, and i band measurements with this value.

For our analysis, we consider the flux densities F_g , F_r , F_i , F_z , F_y , F_W1, and F_W2, and the flux density ratios ${F}_{g}/{F}_{r}$, ${F}_{r}/{F}_{i}$, ${F}_{i}/{F}_{z}$, ${F}_{z}/{F}_{y}$, ${F}_{y}/{F}_{{\rm{W}}1}$, and ${F}_{{\rm{W}}1}/{F}_{{\rm{W}}2}$ as features. In Section 3.3, we choose a subset of these based on their individual information contribution.

The random forest algorithm trains a large set of binary decision trees using a training set with a set of features and known classes or redshift. Each binary tree makes a prediction for the probability distribution of classes or the expected redshift for our quasar candidates. In the scikit-learn implementation (Pedregosa et al. 2011) adopted here, the final pseudo-probability distribution for the classes or the expected redshift is the average of the prediction from each tree. We decide to use the photometric information in the form of flux density ratios as well as the two flux densities, resulting in k = 8 features.

The binary decision trees are built from the training set by determining the best cut along one feature axis via a minimization problem. For the classification, we minimize the Gini impurity $\left(G:= 1-{\sum }_{i}^{k}{p}_{i}^{2}\right)$, and for the regression, the sum of squared errors in redshift. This cut will split the sample in the current node of the tree into two subsamples—its children. The remaining objects in each child give the probabilities for each class as their percentage share of the leaf (p_i ) or the redshift estimate through the average redshift of the objects in the child. The tree will be developed until a stopping condition is reached (e.g., a specified minimum sample size per child node or the maximum depth of the tree). For a quasar candidate, the prediction is based on the p_i or average redshift of the deepest child node to which it belongs.

Single decision trees are prone to overfitting the training data. Hence, a random forest uses ensembles of randomly built decision trees to counteract overfitting. The source of randomization is twofold: (1) Individual bootstrap samples from the training set are drawn to build the trees. (2) Only a subset ($\lfloor \sqrt{k}\rfloor$, in our case) of all k features are considered to find the best split for each internal node. This decreases the running time and reduces the correlation between the individual trees further than just training on the bootstrap samples of the training data. Otherwise, features that are strong predictors would be cut very similar for most trees and thereby result in correlated trees. This was empirically demonstrated by Ho (1998). Correlated trees are undesired since the underlying assumption to be able to average the trees is that they are independent.

One of the main advantages of random forests is that their training is relatively fast. They can run in parallel since the different decision trees can be calculated independently and they scale well. A random forest with T trees and N training objects takes $O(T\,N\,\mathrm{log}\,N)$ time to build and can be applied in $O(T\,\mathrm{log}\,N)$ time.

More details about the random forest algorithm used can be found in Bishop (2006, Chapter 14) and Ivezić et al. (2014, Chapter 9). For this work, the implementation of the random forest classifier and regressor in scikit-learn (version 0.19.1) by Pedregosa et al. (2011) for Python is used. The hyperparameters min_sample_split, max_depth, and n_estimators were optimized for the training set using scikit-learn's GridSearch function. All unmentioned other hyperparameters are left at their default values.

To evaluate the performance of the classification, we use the two measures recall and precision. To evaluate them, we use cross validation. For this, we split our training set into two parts, train the random forest on one part, and then predict the classes for the other part. Considering the resulting true positives (T_p ), false positives (F_p ), and false negatives (F_n ), recall is defined as

$$\begin{eqnarray}&&R:= \displaystyle \frac{{T}_{p}}{{T}_{p}+{F}_{n}}\end{eqnarray} \tag{ 10 }$$

and precision is defined as

$$\begin{eqnarray}&&P:= \displaystyle \frac{{T}_{p}}{{T}_{p}+{F}_{p}}.\end{eqnarray} \tag{ 11 }$$

We are only interested in objects in the high-z class, so we consider the objects put into this class as the positives and all others as the negatives. Each candidate is given a probability to belong to the high-z class by the random forest classification. We set a cutoff probability for the high-z class to decide whether an object will be a valid candidate in our selection. Changing this cutoff probability allows us to increase the recall. However, in return, this will lower the precision of our classification. This is why we need both parameters in order to fully evaluate the performance of the random forest.

We interpret the recall as an estimate of the completeness of our selection, i.e., the fraction of all high-z quasars inherent in our photometric selection the random forest correctly classifies. Similarly, we interpret the precision as an upper limit of the efficiency, where the efficiency is the fraction of the final candidates that are high-z quasars. We will use the terms completeness and upper limit on efficiency in the following, and give a justification for our interpretation in Section 3.5.

To determine which features to use for our analysis, we run the random forest classification with all 13 available flux densities and flux density ratios. In Table 3, we show how much information gain each feature gives relative to the others. We calculate these importances as the (normalized) total reduction of the splitting criterion. As visualized in Figure 2, the different classes can be distinguished by their colors. We see this reflected in the feature importances: the flux density ratios lead to the most information gain for the random forest (see Table 3). The flux density ratios alone, however, remove the brightness information that would capture a luminosity evolution for different classes. To capture the brightness information, we choose to use one flux density per survey. Since each flux density has approximately equal importance, we choose the z and W1 bands, where our targeted quasars have the most reliable detections in Pan-STARRS and ALLWISE. Therefore, we choose to use all flux density ratios as well as F_z and F_W1.

Table 3. To Select the Features, We Run the Random Forest Classification with All Available Flux Densities and Ratios

Feature	Importance (%)
${F}_{g}/{F}_{r}$	17
${F}_{y}/{F}_{{\rm{W}}1}$	17
${F}_{r}/{F}_{i}$	14
${F}_{i}/{F}_{z}$	9
${F}_{z}/{F}_{y}$	9
${F}_{{\rm{W}}1}/{F}_{{\rm{W}}2}$	8
Fi	4
FW2	4
Fy	4
FW1	4
Fz	4
Fr	3
Fg	3

Note. The importance is calculated as the (normalized) total reduction of the splitting criterion brought by that feature. We decide to use all flux density ratios as well as F_z and F_W1.

Download table as:  ASCIITypeset image

Our training data is labeled with the classes given in Table 5. It is worth investigating if reformulating the problem as a binary problem (high-z quasar versus other) or three-class problem (high-z quasar versus other quasar versus star) would improve our results. For each case, we run a fivefold cross validation with our random forest. We calculate a range of statistics summarized in Table 4, with errors giving the standard deviation of the five runs. Precision and recall are calculated as defined in Section 3.2. The "macro" subcolumn is the average over all classes weighted equally, and the "high-z" subcolumn is the precision and recall just for our targeted class. As we reduce the complexity of the classification problem by using fewer classes, the macro metrics should get better because confusion between classes that we combine gets ignored. This is the case: for both the binary and three-class problems, the macro statistics are significantly better. However, the recall and precision for the high-z class are not getting better. In fact, the recall in the binary case is actually lower with 1.6σ. We conclude that, in our analysis, we do not see an improvement from reducing the number of classes, so we will use the full set of classes as defined in Table 1.

Under the assumption that the training sample represents the true distribution of objects on the sky and all sources are real, cross validation of our training sample can predict the performance of the algorithm. This means that, under the assumption that our training set contains a representative set of quasars, our definition of completeness is a reasonable measure for the fraction of all findable quasars that our random forest correctly identifies. However, we will overestimate the efficiency of our algorithm if we measure it with the precision, because the number of M stars is even more dominating on the sky than in our training data and we are neglecting artifacts in the Pan-STARRS+WISE data set. This is why we identify the precision as an upper limit for the efficiency and will use a different approach to get a more realistic estimate in Section 4.

Still, cross validation lets us evaluate the strengths and weaknesses of the algorithm. We train the random forest classification and regression on a random subsample of 80% of the training set and apply it to the remaining 20%. We can then compare the predicted class and redshift with the true ones. By separating the data for training and testing, cross validation avoids biased results from overfitting the data.

The random forest classification assigns each test object probabilities for each class. Later, we can simply select our high-z candidates by applying a threshold on the high-z probability. First, however, we are interested in a comparison of all classes, so the most logical choice is to assign each object to the class with the highest probability. In Figure 3(a), we show the results of the classification on the cross-validation set in the form of a confusion matrix. The matrix depicts how many objects of the true label class on the y-axis are classified to belong to the predicted label class on the x-axis. The majority of objects fall into the diagonal fields, demonstrating that our classifier assigned the correct label to them. The confusion between different types of stars is not concerning for our goal. Confusion between neighboring redshift bins is largely the result of objects right at the border between them, and therefore is also not concerning. As expected, the most relevant contaminants for high-z quasars are M, L, and T dwarfs, as the random forest classifies more than 15% of high-z quasars into those classes. In this case, there is only one star labeled as a high-z quasar, but the balance between completeness and efficiency depends on how we define the cutoff probability for a high-z classification. Therefore, the highest-probability approach chosen here maximizes efficiency for lower completeness. Since M, L, and T dwarfs far outnumber the high-z quasars on the sky, the contamination will become significant at the redshift region of the strongest overlap in color space at around z ≈ 5.4. Since the random forest regression will predict these objects around the same redshift, we will be able to exclude a large fraction of contaminants based on the regression by excluding highly contaminated redshift regions. A common approach to quantify the balance between completeness and efficiency is the ROC curve, which in our case is an almost perfect step function. The ROC score (area under the curve) for the high-z class versus the others is 0.99993.

Zoom In Zoom Out Reset image size

We further analyze how accurately the random forest regression predicts the redshift for a cross-validation set of 20%. We will differentiate two versions of the random forest regression. First, the full regression, where we train on quasars from the full range and can predict the redshift of any quasar. Second, the high-redshift regression, where we only train with z > 4.5 quasars and use it to predict the redshift of objects with class high-z. While the former covers a larger redshift range, the latter provides more accurate redshift estimates for the high-z class candidates, because low-redshift outliers in the full regression training set can skew the result to lower redshifts.

The top of Figure 3(b) shows the distribution of the difference between the predicted and true redshifts for the high-redshift regression when applied to the 20% test set of the training data. Ninety percent of objects have predictions within 0.2 of the true redshift. Since the algorithm can only find new quasars that look similar in color space to the training set, it is to be expected that this represents the performance of observations. As we will see in Section 6, this accuracy is consistent with test observations. It should be noted that the accuracy of our redshift estimate is a strong function of the redshift. This is highlighted in the bottom of Figure 3(b). Here, we show the absolute error of the prediction versus the true redshift. For this, we ran the prediction multiple times with different training set splits, and then averaged the error over bins with equal amounts of objects. One outlier around redshift 4.7 was removed. We expected this increase of error with redshift because (1) there are fewer training objects at higher redshift and (2) higher-redshift quasars appear fainter and thus have higher photometric uncertainties. We also note that, especially toward the high-redshift end, the number of training and test objects becomes very small, so overfitting and redshift gaps in the training set can lead to significant additional inaccuracy in the redshift estimate when applying it to new data that is not captured in our cross validation. The total training set for the high-redshift regression has 695 objects, only 50 of which are above redshift 6. The full regression applied to the same cross-validation set of z ≥ 4.5 quasars gives similar results, but with a bias toward lower redshifts. In this case, only 69% of cross-validation objects have predicted redshifts within 0.2 of the true redshift. In addition, the mean of the predicted redshifts is too low by δz = −0.24, because some objects are incorrectly predicted to be very low-redshift quasars. Therefore, we decided to use the high-redshift regression for our candidate selection.

We now apply the random forest classification and regression algorithms to our full Pan-STARRS+WISE photometric catalog data. To evaluate the completeness of our selection, we again split our training data into two parts: one for training and one for evaluating the completeness. We decide to use the objects within two stripes (R.A. ≤ 60° or R.A. ≥ 300) as well as (−1.26 ≤ decl. ≤ 1.26) for the evaluation. This includes the Stripe 82 area, which has been carefully surveyed for high-redshift quasars and thus makes our completeness estimate more reliable (McGreer et al. 2013). Overall, we use about ∼22% of the training set for evaluation and the rest to train the algorithm.

Our selection picks up larger numbers of candidates in regions of high Galactic extinction and near Andromeda. Therefore, we decide to apply additional restrictions on our photometric catalog data. We require a separation of at least 30° from the galactic center and a separation of 5° from Andromeda (0^h 42^m 44^s +41^d 16^m 9^s), and apply a dust extinction cut of E_B–V  < 0.1.

After we removed the described areas, our final photometric catalog includes around 59 million objects, covering 45% of the sky. ¹⁴ Table 5 shows the results of our random forest classification when assigning the class of the highest probability to the source. As expected, the predicted M stars far outnumber our predicted quasars in the high-z class. However, our training set overrepresents high-z quasars, therefore we expect the number of good high-z quasar candidates to be even lower. Similarly, since our training set underrepresents L and T dwarfs in comparison to high-z quasars, the number of predicted brown dwarfs is much less than high-z quasars, even though we know from observations that it is the other way around.

The random forest classification algorithm provides us with a pseudo-probability for each class. So far, we have simply assigned the class of highest probability, but now we instead look directly at the high-z quasar class probability. Putting a cutoff on this pseudo-probability lets us make a candidate selection where the cutoff can be tuned to our choice of efficiency versus completeness. While these pseudo-probabilities provided by the random forest depend on the input training set and cannot be trusted to represent an absolute measure, a greater high-z class probability makes for a better quasar candidate. Therefore, we can improve the efficiency of our selection by increasing the cutoff on the high-z class probability. At the same time, an increase in the cutoff will reduce the completeness because we are excluding more objects.

Figure 6 shows the probability for the high-z class versus the predicted redshift using the high-redshift regression. We show the majority of candidates with a contour plot to visualize where the density of candidates is highest. To estimate the probability density, we use a Gaussian kernel density estimation applied to all candidates with high-z probabilities above 15%. We then show the probability density contours for three arbitrary density levels that are increasing by factors of 10. This way, we can directly see that there is a large overdensity of candidates around redshift 5.4 and high-z probabilities around 20%. All candidates outside of the lowest probability contour are directly plotted as black dots. At the low-redshift edge, the number of objects with large high-z probability drops off due to the transition from the high-z to the mid-z class. We also find only a few high-z candidates beyond redshift 6.2. This is expected, as only one of our known z ≥ 6.3 quasars passes our photometric requirements on the Pan-STARRS+WISE catalog data. In general, we expect a monotonous decrease in candidates with redshift, since they are fainter. We identify an overdensity of high-z quasar candidates at z ∼ 5.4. There, the trend of monotonous decrease in candidates is interrupted, and toward low high-z probability, the number of candidates goes up much faster than at lower or higher redshift. There is no physical reason to expect many more quasars at that redshift, therefore we are likely seeing significant contamination from stars with similar colors. This is consistent with our evaluation of the cross validation of the random forest classification (Section 3.5): around z ∼ 5.4, the contamination fraction rises because of the photometric similarity between M stars and quasars at this redshift.

Zoom In Zoom Out Reset image size

For the quasar candidate selection presented here, we have chosen to divide all candidates into two separate redshift ranges and treat them separately: 4.8 < z ≤ 5.6 and 5.6 < z ≤ 6.3. Our choice is motivated by the sharp drop of candidate density around z = 5.6. In our final selection, no cross-validation quasar is predicted to be in the wrong redshift range, allowing us to calculate the completeness for both sections separately.

We decide on the cutoff for the high-z probability for each range by evaluating the efficiencies for a range of cutoffs. For our method of estimating the efficiency based on the sky distribution discussed in Section 4, we need to first remove remaining artifacts and blended or extended sources. For this, we visually inspect image cutouts of the Pan-STARRS and WISE photometry for a manageable amount of objects. We inspect the objects with high-z probability above 0.6/0.4 for the lower/higher redshift range. We remove an object if one of its Pan-STARRS images has an artifact interfering with the observation. We consider a nearby Pan-STARRS source detected in the z or y band as blended if it is within the 1σ radius of the PSF fit to the WISE source. We also remove all objects that are clearly extended in multiple bands of the Pan-STARRS imaging.

Then, we calculate the efficiency for a range of values to choose an optimum probability cutoff. The efficiency calculation follows Section 4 and uses the sky area of our Pan-STARRS+WISE catalog data. It covers about 45% of the sky and is defined by:

• 1.  
  Decl. > −30 deg;
• 2.  
  $| b| \geqslant 20\,\deg$;
• 3.  
  >30° angular distance from the galactic center;
• 4.  
  >5° angular distance from Andromeda;
• 5.  
  E_B–V  < 0.1.

The efficiency estimates for different cutoffs on the high-z class probability are shown in Figure 7. The uncertainties on the efficiency reflect the 50% confidence interval. We expect lower values for the probability cutoff to include more contaminants, resulting in a lower selection efficiency.

Zoom In Zoom Out Reset image size

This is exactly what we see in the lower redshift range (4.8 < z ≤ 5.6) of our selection. The lower the probability cutoff is, the lower is our estimated efficiency. In this redshift range, the efficiency declines steeply for efficiency cutoffs below ∼80%. Therefore, we choose 80% as our lower cutoff, as indicated by the blue line in Figure 7.

In the higher redshift range (5.6 < z ≤ 6.3), much lower cutoffs still are predicted to have high efficiency. We choose the minimum lower limit tested: 40%. We note that, at first glance, it seems like the efficiency drops for higher cutoffs, which is not expected. However, we argue that this just reflects the increase in uncertainty for higher cutoffs, since very few objects remain. The full interval for each efficiency prediction in the higher redshift range is consistent with 1.

The lower limit on the high-z class probability is also indicated by the blue line in Figure 6. We retrieve a total of 617 candidates, of which we removed 102 during the visual inspection above (35 image artifacts, 42 blended sources, 25 extended sources). A total of seven known quasars are removed in that process. However, we do not relax our criteria on the visual inspection process, to only select candidates with highly reliable WISE photometry. Of the seven known quasars we removed during visual inspection, six were removed because they are blended in WISE. The seven known quasars removed during visual inspection represent 3.0% of the known quasars in the candidate set, while overall 16% of the candidate set is removed. We interpret this as an indication that the processing step is indeed reducing the fraction of contaminants in the final selection.

In the end, we select a total of 515 promising quasar candidates, which we call the high-z candidate set. Of these, 226 (or ∼43%) are already known quasars, demonstrating the success of our selection method.

For this sample of 515 quasar candidates, we now estimate the completeness and the efficiency. We estimate the completeness with the known quasars that we withheld from the training set. In particular, we define the completeness as the fraction of these known quasars that are in our final candidate set. The quoted uncertainties represent the 1σ confidence interval.

Since we are dealing with small sample sizes, we use the Wilson interval to estimate the confidence interval for this binomial distribution, following the recommendation of Brown et al. (2001). For large enough data sets, this converges back to the usual standard deviation of a Gaussian. For small data sets, it better captures the asymmetry in the error while retaining that the 1σ range captures 68.27% of data points. We calculate a completeness of 66% ± 7% for the redshift range of 4.8 < z ≤ 5.6, and a completeness of ${83}_{-9}^{+6} \%$ for the higher redshift range ($5.6\lt {\text{}}z\leqslant 6.3$).

We show the completeness as a function of redshift in Figure 8. To calculate this, we applied a kernel-density estimate (kde) to both the targeted known cross-validation quasars remaining in the selection and in total. The ratio then gives our completeness estimate. For the kde, we used Gaussian kernels with equal weights for all points and bandwidths chosen with Scott's rule (Scott 1992). Below redshifts of z ≈ 5.6, the completeness is nearly constant around a value of ∼67%. Above z ≈ 5.6, it rises to peak around 88% at z ≈ 5.9. This behavior simply reflects that, above predicted redshift z = 5.6, we accept candidates with a lower high-z class probability.

Zoom In Zoom Out Reset image size

Above redshift 6, however, the completeness declines sharply. While this behavior is estimated based on only three cross-validation quasars with z ≥ 6, it signals that our method stops being effective at z ≥ 6. Potentially, the small number of z ≥ 6 quasars in our training set (52 total) might not allow for proper classification using the random forest method.

We have fine-tuned our selection, in particular the lower limit on the high-z class probability, to ensure high selection efficiencies. Indeed, our selection efficiency in the lower redshift range (4.8 < z ≤ 5.6) is on average ${78}_{-8}^{+10} \%$, and for redshift of $5.6\lt {\text{}}z\leqslant 6.3$ we reach ${94}_{-8}^{+5} \%$ (Figure 7). The quoted uncertainties correspond to the 68% confidence interval. It is a good check for consistency that our method predicts efficiencies that are at least as high as the fraction of known quasars in our selection: The fraction of known quasars in the final set is 42% in the lower redshift range and 50% in the higher redshift range, respectively.

Next, we estimate the redshift dependence of the efficiency of our candidate set. The efficiency is the fraction of our candidates that are actually quasars. To estimate this, we again use our method from Section 4. We calculate the efficiency for bins with a width of 0.1 in redshift. Figure 9 shows our selection efficiency as a function of redshift. Part of our candidate set consists of known quasars, which tells us that the efficiency of the selection is at least as large as the fraction of known quasars in that bin. We show this minimum efficiency in orange. Whenever the lower efficiency limit and our estimated efficiency agree, we do not expect to find new quasars in that redshift bin. When the estimated efficiency is larger, we do expect the candidate set to contain quasars that are not yet known. In the redshift bin of z = 5.4–5.5, we see the selection efficiency drop to the lowest value in our entire redshift range. This is likely the result of the significant overlap in color space of quasars with M stars at z ≈ 5.4, as we discussed in Section 3.5. Based on our efficiency estimate, we expect to find the highest-redshift quasars with this selection around 5.5 < z < 5.8, where our predicted selection efficiency is above the lower limit. At z ≥ 5.8, the efficiency prediction and the lower limit are consistent with each other. Therefore, we expect to find few new quasars at z ≥ 5.8. Finally, at the low end of our targeted redshift range, we also expect new quasars, since the efficiency estimate is well above the lower limit.

Zoom In Zoom Out Reset image size

Based on the estimate of the efficiency we can predict that our high-z candidate sample contains ${319}_{-33}^{+41}$ quasars at 4.8 < z ≤ 5.6 and ${100}_{-8}^{+5}$ at 5.6 < z ≤ 6.3. Subtracting the known quasars, we expect that our candidate set contains ${148}_{-33}^{+41}$ and ${45}_{-8}^{+5}$ new quasars in the lower and higher redshift ranges, respectively, where the error is a 68% confidence interval. Table 6 summarizes our predictions for the selection. We deliver the paper with a data file containing the full high-z candidate set. Table 7 describes the columns of the data file.

Table 6. Summary for the High-z Candidate Set

Redshift range	4.8 < z ≤ 5.6	5.6 < z ≤ 6.3
Number of candidates	409	106
Completeness	66 ± 7%	${83}_{-9}^{+6} \%$
Efficiency	${78}_{-8}^{+10} \%$	${94}_{-8}^{+5} \%$
Known quasars	171	55
Predicted new quasars	${148}_{-33}^{+41}$	${45}_{-8}^{+5}$

Note. The calculation of these properties is discussed in Section 5. All errors give 68% confidence intervals.

Download table as:  ASCIITypeset image

We select a subset of our final high-z candidate catalog and require all candidates to have prediction redshifts of ${z}_{\mathrm{RF}}\gt 5.6$ in both the full regression, using training quasars at all redshifts, and the high-redshift regression, using training quasars at z > 4.5. We retain 59 promising quasar candidates, which nominally have a selection efficiency of ${86}_{-34}^{+11} \%$ as estimated by our new method (see Section 4). At the time of our selection, 32 of the candidates were already known quasars, leaving 27 unknown objects. From our efficiency estimate, we expect through error propagation that, for these 27 objects, our success rate to find quasars should be ${69}_{-52}^{+24} \%$, with the 68% confidence interval as the error. We can compare this to a naive estimate based on Equation (11), which would lead us to predict 100% for the efficiency because there is no known star from our test set in our final selection (see Figure 6).

One of these candidates, J112143.62-071839.4, was recently identified by Yang et al. (2019) as a z = 5.71 quasar.

The J-band is very effective in differentiating the classes because red and brown dwarfs tend to have more flux in their spectrum for this band than otherwise similar-looking quasars. We followed up on six candidates with J-band photometry using NOTCam. The results are summarized in Table 8.

Figure 10 shows the zJ–JW2 color–color diagram combing the Pan-STARRS magnitudes with the J band from our NOTCam follow-up observations (blue dots). As a point of comparison, we also plot our known objects for which J-band information from 2MASS is available (Skrutskie et al. 2006). ¹⁵ We note that 2MASS is shallower than our follow-up observations. The known quasars and stars have mean J-band AB magnitudes of 17.6 and 16.5, while our six follow-up observations have a mean of 19.6. However, to first order, quasars at the same redshift have similar colors with only minor luminosity evolution. Promising quasar candidates can be separated from likely dwarf stars with a color cut shown in black ($z-J\lt 1.9$ and $J-{\rm{W}}2\gt -0.1$, both in AB mag). Five out of our six observed objects make the color cut. A close-by source is evident in the J-band photometry of our only nonpromising candidate, J145950.96-181251.7. Therefore, it is likely blended in WISE, which could explain the false classification. Overall, our NOT photometric follow-up observations indicate that our candidate set does contain promising candidates.

Zoom In Zoom Out Reset image size

Furthermore, we were able to obtain five follow-up spectra of our selection. These objects were not prioritized by the high-z probability; we observed the objects in the candidate set with the best visibility at the observatories. We identified three objects as contaminants and two as quasars at z ∼ 5.7.

Counting the Yang et al. (2019) quasar and the likely quasar at lower redshift, the selection efficiency would be three out of six (or 50%). Since we did not observe the sixth object, a more conservative counting would be two new quasars in the targeted redshift range out of five observed, or a selection efficiency of 40%. From a naive approach, we would have expected an efficiency of 100%, while our method for estimating efficiency predicted ${69}_{-52}^{+24} \%$. While our very small sample size does not allow conclusions about the accuracy of our approach, we do argue that our method gives more realistic results that are consistent with our small test observation.

In the following, we discuss the two newly discovered quasars individually. We note that we observed a sixth object: J032615.68-061358.2. The continuum looks like a power law; however, we do not see a Lyα break in our spectral range. This indicates that it likely is a quasar but at z < 5.4 where the Lyα line moves out of our spectral range. The predicted redshift of z = 5.62 was too high. We do not consider it for our efficiency test here, because we cannot confirm the classification. In Figure 6, we show the object with a cross.

6.1.1. J152330.66+293539.1—z = 5.73

J152330.66+293539.1 is a newly discovered quasar at redshift 5.73 based on Lyα emission. The predicted redshift was 5.72 and the high-z probability was 0.83. The object was part of our NOT photometric follow-up observations, where we measured the J-band magnitude. The obtained colors where $J-{\rm{W}}2=0.80$ and $z-J=-0.32$ in AB, which are consistent with quasars at this redshift as seen in Figure 10. We observed this object with the MODS spectrograph on the LBT, and we present the spectrum in Figure 11 together with the other discovered quasars.

Zoom In Zoom Out Reset image size

This object has a relatively blue color of $i-z=1.84$, so it would not be part of a typical color cut selection like Bañados et al. (2016), where candidates were cut at $i-z\gt 2$. This indicates that our method can find quasars that are missed by traditional color cuts even if our random forest is trained with objects largely from these selections. Our full high-z candidate set that we publish with this work contains a further 37 candidates with predicted redshift above 5.6 that do not fulfill this color cut. There are only seven known quasars above redshift 5.6 that do not fulfill the color cut.

6.1.2. J163752.18+024158.1—z = 5.76

We have tested our random forest selection method with pilot observations during the development process. We used a preliminary version of the algorithm to make a selection targeted at redshift 4.8–5.4 and obtained 31 optical spectra, out of which we successfully identified 17 new quasars. Eight of these quasars are retained within our final high-z candidate set, but none of the contaminant objects are selected anymore. This indicates that our selection improved in robustness. The newly discovered quasars, which did not make it into our final selection (a total of nine), are either at lower redshift than our targeted selection (three have observed z < = 4.8) or are very close to the redshift boundary. A total of eight have observed redshifts below z = 4.92 and hence their classification shifted toward the mid-z class. Another newly discovered quasar at z = 5.03 just barely missed the 80% cutoff on the high-z probability (79.6%) and thus was excluded from our final candidate list.

These preliminary observations were already very successful, with 55% of observed candidates being newly discovered quasars. With the improvements to our selection discussed above, the final selection is expected to be even better in this redshift range around $z\approx 5$.

Furthermore, we were able to obtain follow-up observations for one additional object in the lower redshift range of our final high-z candidate set using FORS2. We identify the object, J110942.97-285521.0, as a quasar at z = 5.01. The predicted redshift was 5.09. In the following, we discuss the discovered quasars.

6.2.1. J001150.03-244400.1—z = 5.41

6.2.2. Seventeen New Quasars at 4.6 ≤ z ≤ 5.1