The Quasar Catalogue for S-PLUS DR4 (QuCatS) and the estimation of photometric redshifts

The advent of massive broad-band photometric surveys enabled photometric redshift estimates for unprecedented numbers of galaxies and quasars. These estimates can be improved using better algorithms or by obtaining complementary data such as narrow-band photometry, and broad-band photometry over an extended wavelength range. We investigate the impact of both approaches on photometric redshifts for quasars using data from Southern Photometric Local Universe Survey (S-PLUS) DR4, Galaxy Evolution Explorer (GALEX) DR6/7, and the unWISE catalog for the Wide-field Infrared Survey Explorer (WISE) in three machine learning methods: Random Forest, Flexible Conditional Density Estimation (FlexCoDE), and Bayesian Mixture Density Network (BMDN). Including narrow-band photometry improves the root-mean-square error by 11% in comparison to a model trained with only broad-band photometry. Narrow-band information only provided an improvement of 3.8% when GALEX and WISE colours were included. Thus narrow bands play a more important role for objects that do not have GALEX or WISE counterparts, which respectively makes 92% and 25% of S-PLUS data considered here. Nevertheless, the inclusion of narrow-band information provided better estimates of the probability density functions obtained with FlexCoDE and BMDN . We publicly release a value-added catalogue of photometrically selected quasars with the photo-z predictions from all methods studied here. The catalogue provided with this work covers the S-PLUS DR4 area ( ∼ 3000deg 2 ), containing 645 980, 244 912,


INTRODUCTION
Quasars are known to be the most energetic objects in the Universe (Kormendy & Ho 2013).Since the first quasar discovery by Schmidt (1963), the number of confirmed quasars increased to almost a million.Among many surveys that have contributed to the discovery of new quasars, the Sloan Digital Sky Survey (SDSS; York et al. 2000) has provided the largest number of confirmed quasars with ★ Email: lilianne.nakazono@usp.br the SDSS-I/II/III/IV phases, including the extended Baryon Oscillation Spectroscopic Survey (eBOSS; Dawson et al. 2016).The SDSS DR16 quasar catalogue from DR16 (DR16Q; Lyke et al. 2020) includes objects from all SDSS programmes and contains 920 110 spectroscopic observations of 750 414 quasars.The SDSS success in finding quasars is due to multiple efficient techniques for candidate quasar selection and redshift estimation based on optical and mid-infrared colours (see Myers et al. 2015 for an overview).
Information needed to select candidate quasars for follow-up spectroscopic observations is typically provided by photometric sky sur-veys.Furthermore, if adequate colours are available, quasars redshifts can be estimated using photometry alone, thus bypassing the need for expensive spectroscopic observations.Techniques to estimate redshifts using photometric observations (photometric redshifts or photo-zs) are well studied in the literature, but they remain an active research topic.There are two main approaches to photo-zs: SED (spectral energy distribution) template fitting (e.g.Salvato et al. 2009, Babbedge et al. 2004) and empirical methods based on training samples (e.g.Wu & Jia 2010, Brescia et al. 2013, Nakoneczny et al. 2021, Yang & Shen 2022), each with their own pros and cons.In general, machine learning as an empirical approach provides more accurate photo-z predictions than SED template fitting (e.g.Schmidt et al. 2020).However, the robust performance of such methods is limited to the parameter space of the training set, while SED template fitting has no such limitations (see Brescia et al. 2021 for further discussion).
In this work, we analyse the performance of photometric redshift estimates based on 12-band photometric data from the Southern Photometric Local Universe Survey (S-PLUS; Mendes de Oliveira et al. 2019).S-PLUS is a large-area sky survey that will cover ∼9300 deg 2 of the Southern Sky with the same optical system used in the Javalambre Photometric Local Universe Survey (J-PLUS; Cenarro et al. 2019).This optical system consists of seven narrow bands and five  broad bands ( are similar to the SDSS bands; Fukugita et al. 1996), providing effectively low-resolution spectra that are expected to improve photometric redshift estimates compared to those based on broad-band photometry alone.We also extend S-PLUS data with photometry from Galaxy Evolution Explorer (GALEX; Martin et al. 2005) and Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010) surveys to study how photo-z estimates improve due to an extended wavelength range compared to improvements due to the addition of narrow-band photometry.We focus here on photo-z techniques based on machine learning.
We tested three different and independent machine learning algorithms: Random Forest, Flexible Non-parametric Conditional Density Estimation (FlexCoDE), and Bayesian Mixture Density Network.We discuss the advantages and disadvantages of each method, and in particular assess the importance of the narrow-band photometry for the photo-z regression.Using three independent methods allows us to check for consistency in the interpretation of our results.We provide a value-added catalogue of quasar candidates for spectroscopic followup with their assigned probabilities of being a quasar (obtained from the star/quasar/galaxy classification by Nakazono et al. 2021) and their estimated photometric redshifts from this work.Estimates from all three methods are provided separately in the catalogue.
This paper is organised as follows: data sets used to estimate photozs are described in Section 2 and photo-z methods are described in Section 3. We present our results in Section 4 and how to select quasar candidates in Section 5. Finally, our main findings are discussed in Section 6.

PHOTOMETRIC SURVEYS AND DATASETS
In this section, we describe three photometric surveys, S-PLUS, WISE, and GALEX, that provided photometry for photo-z estimation.In addition, a reliable sample of known quasars with spectroscopic redshifts is needed for training supervised machine learning methods and we used the SDSS DR16Q quasar catalogue (Lyke et al. 2020).In Table 1, we list the bands used in this work together with their effective wavelength, width and magnitude depth.

The Southern Photometric Local Universe Survey (S-PLUS)
The S-PLUS survey is an ongoing survey that, by the end, will cover ∼9300 deg 2 of the Southern Sky with the T80-South telescope located at Cerro Tololo Inter-American Observatory (CTIO), Chile.The S-PLUS optical filter system includes five broad-bands: , , , ,  (similar to the SDSS filter system, except for the  band) and seven narrow-band filters: 0378, 0395, 0410, 0430, 0515, 0660, and 0861, centred on [OII], Ca H+K, H, G-band, Mgb triplet, H and Ca triplet features, respectively (Marín-Franch et al. 2012).
We use aperture fluxes computed within a circular diameter of 3 arcsec and corrected for the fraction of the flux falling outside this limit (see Almeida-Fernandes et al. 2022 for more details).The procedure for data reduction and calibration of S-PLUS Data Release 4 (DR4), used in this work, is further detailed in Herpich et al. 2024.The S-PLUS DR4 comprises ∼3000 deg 2 , including the area from past S-PLUS data releases.There are 645 980, 244 912, 144 991 sources photometrically classified as quasars with 80%, 90%, and 95% probabilities in S-PLUS with  < 21.3 and good photometry flag in the detection image (SEX_FLAGS_DET = 0).We corrected the data for reddening using dust maps from Schlegel et al. (1998) and the extinction law from Cardelli et al. (1989) In Fig. 1 we show the SEDs of five examples of quasars from our sample with redshifts of 1. 639, 2.238, 2.516, 3.207, and 4.435.The figure shows how the S-PLUS SEDs are shaped when one or more emission lines fall in the region of a narrow band.For instance, for  ∼ 1.63, CIV is detected in 0410; for  ∼ 2.24, Ly is detected in 0395; for  ∼ 2.52, Ly is detected in 0430; for  ∼ 3.21, Ly is detected in 0515.At the red side of the optical wavelength range, S-PLUS has the 0660 narrow band with a magnitude limit of about 21.1 mag, which is deeper than all other narrow bands, and deeper than the  band (see Table 1).Note from the bottom panel of Fig. 1 that the Ly can be detected at the 0660 band in  ∼ 4.4.Considering the large area of the southern sky that will be surveyed by S-PLUS, it will play an important role in finding quasar candidates that are under-represented in the SDSS spectroscopic sample.639, 2.238, 2.516, 3.207, and 4.435 (top to bottom) showing an emission line being detected at an S-PLUS narrow band.In gray, we plot the spectra observed by the SDSS programs, with the identified emission lines being plotted in vertical dashed lines.The coloured triangles (squares) represent the flux measured in the S-PLUS narrow (broad) bands, from left to right: , 0378, 0395, 0410, 0430, , 0515, , 0660, , 0861, .The filter responses are shown accordingly.Empty triangles/squares represent the error in AB magnitude above 0.5 or non-detection at that particular band.

The Wide-field Infrared Survey Explorer (WISE)
The WISE (Wright et al. 2010) is an all-sky infrared survey that obtained photometry in four bands, in the wavelength range from 3.4 m to 22 m.In this work, we use the W1 and W2 bands from the unWISE catalogue (Schlafly et al. 2019), with effective wavelengths 3.4 m and 4.6 m, respectively.The addition of these two infrared bands is known to significantly improve the performance of photo-zs for quasars (e.g.Bovy et al. 2012, Brescia et al. 2013, DiPompeo et al. 2015, Yang et al. 2017).WISE magnitudes are reported on the Vega System and their conversion to AB system is given by Eq. 1: where Δ = 2.699 and Δ = 3.339 for W1 and W2, respectively.It is recommended to subtract a 4 mmag and a 32 mmag offset to the unWISE measurements (Schlafly et al. 2019).The Vega to AB offsets, as constants, do not affect machine learning photo-z methods as long as the feature space for the training set and for the unlabelled data are defined in the same way.Therefore, we do not convert these magnitudes to the AB system for this work.To apply our trained models in all sources observed in S-PLUS, we perform a CDS crossmatch with the unWISE catalogue within 2 ′′ .There are 75% out of 39 168 373 (14 <  < 21.3) sources with WISE counterpart in the S-PLUS DR4 area but only 54% have detection in the W2 band.

The Galaxy Evolution Explorer (GALEX)
The GALEX (Martin et al. 2005) is an all-sky ultraviolet survey that obtained photometry in two bands: far-UV () and near-UV ().We retrieved data for both bands from GALEX DR6+7 (Bianchi et al. 2017) using CDS cross-match within 2 ′′ .The addition of these two UV bands is known to significantly improve the performance of photo-zs for quasars (e.g.Ball et al. 2008, Bovy et al. 2012, Brescia et al. 2013).There are 8% out of 39 168 373 (14 <  < 21.3) sources with GALEX counterpart in the S-PLUS DR4 area but only 2% have measurements in the FUV band.

Spectroscopic quasar sample
We cross-matched SDSS DR16Q with S-PLUS DR4 within 1 ′′ radius, using data from the SDSS Stripe 82 equatorial region for our analysis (∼300 deg 2 ) resulting in 33 151 matches.Stripe 82 region dominates the current overlap between SDSS DR16Q and S-PLUS DR4, and represents the best-studied region in the SDSS footprint.We selected a total of 33 151 quasars with  ≤ 22 and  ≤ 5 for our analyses, which were cross-matched with unWISE and GALEX within 2 ′′ .A fraction of 93.1% and 85.9% out of 33 151 have detection in 1 and 2 bands from unWISE, respectively.For GALEX, only 30.8% and 11.1% have detection in  and  bands, respectively.Note that GALEX does not provide full sky coverage and therefore some of the missing-band values are due to nonobservation.Ideally, non-detection and non-observation cases should be distinguished in the training process but determining the exact fraction of non-observation is not a trivial task.For the time being, we treat all non-observation missing values as they were non-detection (see §3.4 for further details).

METHODS
In this work, we assess the importance of narrow-band photometry in improving the quality of quasar photo-z's predictions.First, we discuss how much is the photo-z prediction improved when narrowband photometry is added into a feature space that only contains broad-band photometry.In order to answer this question, we compare single-point estimates from Random Forest experiments trained with fixed hyperparameters and data sets using two different feature colour spaces: • broad:  − ,  − ,  − ,  −  • broad+narrow: the above four broad-band colours extended with the following seven narrow-band colours: 0378 − , 0395 − , 0410 − , 0430 − , 0515 − ,  − 0660, and  − 0861.
We then discuss how the improvement factor changes when the GALEX and WISE broad-band photometry is also available.Similarly, as before, we compare two experiments trained with the following feature spaces: • broad+GALEX+WISE+narrow: the above eight broad-band colours extended with the same seven narrow-band colours as in the first case.
As GALEX and WISE bands improve photo-z's predictions even when narrow-band colours are not included (see Section 4), we only consider the broad+GALEX+WISE and broad+GALEX+WISE+narrow feature spaces for further analyses.Testing three independent methods allows us to check the consistency of our conclusions about the importance of narrow-band photometry.Moreover, different methods can generally learn different patterns from the data and thus we also evaluate an ensemble of the three methods.Lastly, we also estimate photo-z's probability density functions (PDFs) using FlexCoDE and Bayesian Mixture Density Network.
The quasar sample described in Section 2 is randomly split into training (75%) and testing (25%) sets.Their spectroscopic redshift distribution is shown in Fig. 2. Note that there are only few objects in the testing set for  > 4.2, which may affect the performance estimation in this range.Here, we introduce the methods used in this work: Random Forest ( §3.1), FlexCoDE ( §3.2), and Bayesian Mixture Density Network ( §3.3).In §3.4 we discuss how we handled missing-band values and outliers.Finally, in §3.5 we discuss the metrics considered in this work to evaluate our model performance.

Random Forest
Random Forest (RF) is a non-parametric method that is robust to outliers, noise and overfitting.Moreover, it is fast-learning, and, in some sense, interpretable due to the computation of feature importance.These characteristics and its known ability to provide high accuracy in diverse fields of study make RF an attractive algorithm.Since RF has been extensively used in astronomy, for more details we simply refer the reader to Breiman (2001).Regarding the specific task of estimating photometric redshifts, please see: Carliles et al. 2010 (galaxies), Henghes et al. 2021 (galaxies), andNakoneczny et al. 2021 (quasars).
We use the Python implementation of RF in scikit-learn1 (version 1.1.3)package.We utilized the complete colour space (broad+GALEX+WISE+narrow) to run a grid search, resulting in the following configuration: n_estimators=400, min_samples_split=2, min_samples_leaf=2, max_depth=None, bootstrap=True.The photometric redshift is estimated as the mean value of the 400 predictions.We set random_state to 47 for all processes, in order to enable reproducibility.The remaining hyperparameters were kept at default values of scikit-learn v1.3.0.In this work, we do not provide an estimated probability density function for the photometric redshift with RF.

FlexCoDE
FlexCoDE (Izbicki & Lee 2017;Dalmasso et al. 2020) is a nonparametric method that estimates conditional densities.FlexCoDE was applied to estimate photometric redshift distributions in the mock data produced for Vera Rubin Observatory's Legacy Survey of Space and Time (LSST; Schmidt et al. 2020) Let  (|x) be the probability density function of the redshift  of a quasar with photometric covariates x.The key idea of FlexCoDE is to project  (|x) on an orthonormal basis (  )  ∈N (we use the Fourier basis): where   (x) are the expansion coefficients.By construction, these coefficients are given by and thus each coefficient   can be estimated by regressing   () on x using the spectroscopic sample.In this paper, we use Random Forests to estimate each of these coefficients.The estimated density is then given by f (|x) =  =1 β (x)  (), where  is a hyperparameter that defines the total number of expansion coefficients in the model.In practice, we choose the value of  that minimizes the CDE loss function computed on a validation set (Equation 9).The single-point estimate is derived as the mode of the estimated probability density function.

Bayesian Mixture Density Network
We call as Bayesian Mixture Density Network (BMDN) the combination of two types of architectures: a Bayesian Neural Network (Bishop 1997) and a Mixture Density Network (Bishop 1994).BMDN is similar to a Dense network, where all neurons from one layer are connected to all neurons in the previous and following layers.The strength of these connections is given by weights, represented as single values.
In the Bayesian framework, these single-value weights are replaced by probability distributions.In our case, we assume a Gaussian distribution.For the implementation, we used the Tensorflow Probability (Dillon et al. 2017).The Mixture Density section of the network is what allows the estimation of distributions as output.This combination of architectures enables the network to represent arbitrary probability distributions while also allowing for an estimation of the epistemic (due to the model) and aleatoric (due to the data collection process) uncertainties, and provides a way to deal with the colour-redshift degeneracy (see Section 4.2).The architecture used in this work is shown in Fig. 3.
The MixtureNormal layer of the network outputs 7 Gaussian functions described by two parameters: the mean () and the standard deviation ().Each Gaussian has an associated weight (), so the photo-z's probability density function is a combination of these Gaussian functions.For each object, we also derive a single-point estimate by computing the mode of the estimated probability density function.

Dealing with missing-band values
Missing-band values are intrinsic to astronomical data.We assume that any missing data in a non-detection case is "missing not at random" (MNAR; see Rubin 1976 for more information).An example of MNAR is a high-redshift quasar that is not detected in the blue bands due to the Lyman Break (e.g.see the two bottom panels in Fig. 1).These missing values have physical meaning and cannot be ignored.For MNAR, the removal of objects from our training set (also known as the complete-case method), or the mean/median imputation would bias our models (Dong & Peng 2013).
For tree-based methods (such as Random Forest and FlexCoDE), we can handle the non-ignorable missing-band values by replacing them with an arbitrary out-of-the-range value (e.g.99).The colours are calculated after the missing-band replacements.
For the Bayesian Mixture Density Network, we first replace any missing-band values in S-PLUS by the minimum magnitude of a non-detection given in the error column for each band.Then, we standardize the data and re-scale them to the range [0,1], without taking into consideration the missing-band values in GALEX and unWISE.The last scaling process is motivated by the suggested procedure from Chollet (2017) of replacing missing features by zeros.Since this can only be done if the value zero is not meaningful for the problem, and since colours can have this value, we use this second scaling function to shift the distribution.After the scaling processes, we replace any missing-band values in GALEX and unWISE by zeros.

Model Evaluation
Here, we present the performance metrics considered in this work.Let  be the (true) spectroscopic redshift and ẑ, the estimated photometric redshift (sometimes, especially in figures, we also refer to them as z spec and z phot , respectively).
Defining  =  − ẑ,  norm = /(1 + ) and considering  being the number of evaluated objects, we define the following performance metrics for the single-point predictions: • Root-mean-squared error: • Normalized median absolute deviation: • Bias: • Fraction of objects with residual above a cutoff : The  NMAD is commonly used to evaluate redshift estimates as it is less sensitive to outliers (Salvato et al. 2019).Photo-z single-point estimates are better when these metrics are closer to zero.Now, let f be the estimated probability density function, and let F be the estimated cumulative density function.In order to evaluate the probability density functions estimated by FlexCoDE and Bayesian Neural Networks, we use the following metrics: • Conditional Density Estimate (CDE) loss function (Izbicki & Lee 2016;Izbicki, Lee & Freeman 2017), L(  , f ), which is an estimate of:

RESULTS
In this section, we analyse the importance of the narrow-band photometry for estimating quasar photometric redshifts using RF, Flex-CoDE, and BMDN, because different methods might learn different underlying relations between the colour space and the true redshift.
We first use RF with a 5-fold cross-validation method to evaluate all feature spaces described in Section 3.After narrowing the possible choices of feature spaces, we then evaluate the single-point estimates of all the three methods (see Section 4.1).In Section 4.2 we evaluate the estimated probability density functions from FlexCoDE and BMDN.Finally, in Section 4.3 we show the feature importances obtained from the tree-based models.

Assessing the importance of narrow-band photometry
Here, we compare the performance of Random Forest experiments trained on different feature spaces with fixed hyperparameters and fixed training/validation sets.With the choice of feature spaces described in Section 3, we aim to provide insights into the importance of narrow-band photometry for this regression problem.The quantitative comparisons were done with  RMSE ,  NMAD , ,  0.15 , and  0.3 under a 5-fold cross-validation scheme.These metrics are shown in Table 2.When narrow-band colours are added to the broad-band colour space, there is a decrease in the values of  RMSE and  NMAD by about 11% and 17%, respectively.As is discernible from Table 2, in all cases  RMSE is much larger than  NMAD due to the influence of outliers.The impact of outliers is also well-tracked through  metrics.The  0.15 metric decreased from 49% to 43%, while  0.3 also decreased from 23% to 19%.The biases are sufficiently small for all cases with an order of 10 −3 .
Table 3.This table shows the mean performance improvement in per cent due to the addition of narrow-band information for the subset of sources without GALEX ( and ) or WISE (1 and 2) counterparts.Here we only compare broad+narrow and broad colour spaces.The improvements are slightly below the 11% ( RMSE ) and 17% ( NMAD ) calculated from Table 2.For completeness, we also show the results for the complementary subsets.

Sample Subset
Decrease in  RMSE (%) These improvements are smaller than the improvements observed when GALEX and WISE colours are added to broad-band colours (i.e.without using narrow-band colours).As shown in Table 2, the extended wavelength range due to the addition of GALEX and WISE colours improve  NMAD and  0.15 metrics by as much as about a factor of two.When all the data are used, broad-band colours, GALEX and WISE colours and narrow-band colours, only slight additional improvements are observed.These improvements are visualized using the histograms of the residuals ( spec −  phot ) in Fig. 4; the residual distribution is more concentrated around zero when information from all bands is included in our model.Fig. 5 shows how  RMSE and  0.30 vary with the apparent magnitude .The addition of narrow-band and/or GALEX and WISE colours improves the regression performance, especially for brighter objects.
As discussed above, the addition of GALEX and WISE colours to broad-band colours improves the performance of photo-z more than the addition of narrow-band colours to broad-band colours.However, we emphasize that GALEX and WISE bands are effectively shallower than the S-PLUS bands and for faint objects, only narrowband photometry is available, besides broad-band photometry.To address these cases where narrow-band photometry will play a more important role, we calculate the metrics on subsets of the validation folds that do not have GALEX or WISE counterparts (see Table 3).For  RMSE the improvement is of an order of ∼9% for both cases.Regarding  NMAD , the improvement is 14.1% for sources without GALEX counterpart, and 16.5% for sources without WISE counterpart.
For further analyses that include FlexCoDE and BMDN, we only consider the broad+GALEX+WISE and broad+GALEX+WISE+narrow colour spaces.In Table 4 we show the metrics for all three methods; they all show small improvements from the addition of narrow bands regarding  RMSE ,  0.15 , and  0.30 .There is an interesting behaviour observed for  NMAD : without narrow bands, all three methods have similar performance but when narrow bands are added, FlexCoDE and BMDN outperform RF method by a factor of two, approximately.On the other hand, RF method achieves about an order of magnitude smaller bias in comparison to FlexCoDE and BMDN for all models, except for the BMDN with the narrow bands.When we average the predictions from these three models,  RMSE and  0.30 are improved compared to each individual model.Note that the first two rows in Table 4 are not identical to the last two rows in Table 2 because the former is based on the testing set while the latter is based on the validation sets.

Probability density functions
Performance analyses from single-point estimates, however, do not account for the uncertainties in the photo-z predictions.It is expected that the photo-z's PDFs can be multimodal due to the colourredshift degeneracy, which can be seen in Figure 6 for a random selection of nine quasars from the test set at different ranges of spectroscopic redshift and magnitude.For instance, the colours of the quasar SDSS J232043.35-003049.3 are related to, approximately, two apart single-point estimates for FlexCoDE narrow-band model (i.e.broad+GALEX+WISE+narrow), for which the PDF's second highest peak is the one centred on the correct value.Nevertheless, the PDFs' primary peak from FlexCoDE with the narrow-band model covers well the "true" redshift for most of the examples shown.
In comparison with BMDN, the FlexCoDE PDFs are generally more accurate for these particular model and examples.It is notable, however, that the PDFs are less certain on the predicted photo-z (i.e.distributions present more peaks) for the cases in the tail of the spectroscopic redshift distribution of the training sample beyond the  magnitude depth of S-PLUS ( > 21.3).This behaviour is also seen with the broad-band model (i.e.broad+GALEX+WISE) and therefore might be related to the increasing photometric uncertainties at fainter magnitudes.Although the PDFs show more confusion at this range, the single-point prediction still agrees fairly well to the "true" redshift but only with the narrow-band model.
Other quantitative comparisons should also be considered in order to evaluate the whole PDF information, such as the PIT distribution and the CDE loss function.As we have seen in Fig. 6, narrow-band model PDFs are generally sharper than broad-band model PDFs for both FlexCoDE and BMDN, suggesting that the predicted distribution is more precise (but not necessarily more accurate) when narrow-band information is included in the model.On the other hand, this could mean that the predictive uncertainties are not being properly accounted for and the distributions are sharper than they should be.The high edges (U-shape) of the PIT distributions shown in Fig. 7 agree with the suggestion of underdispersed PDFs.From this figure, we see that the narrow-band PDFs are more underdispersed than the broad-band PDFs.It then requires a better calibration of these PDFs that can either be achieved with a different machine learning architecture or with a recalibration process (e.g.Dey et al. 2021), which is planned to be pursued in future work.
Nevertheless, we showed good evidence that the inclusion of narrow-band photometry is valuable for the precision of photometric redshift estimates.This is also supported by the calculated CDE losses, for which each individual PDF is evaluated (note that the loss can be estimated even if the true PDF is unknown).The CDE loss functions for FlexCoDE are −3.27 and −1.36 when trained with and without narrow bands, respectively.For BMDN, these values are −2.47 and −1.31.Therefore, the FlexCoDE narrow-band model provided the best photo-z predictions during the testing process.
Given that the feature space and the training sample are the same, the differences in PDF shapes are model intrinsic.Different algorithms can learn different underlying patterns from the data, hence these PDFs could be combined to provide more powerful information.Assessing the possible combination of these PDFs is beyond the  scope of this paper, but we provide the PDFs from both methods in our value-added catalogue (see Section 5).Based on the CDE loss, the FlexCoDE trained with narrow bands delivered the best photo-zs and is our general recommendation for broad usage.

Feature importances from tree-based models
In Fig. 8, we show the estimated feature importances for the Random Forest and the FlexCoDE models.In both models, the top two most important features for the quasar photo-z estimation are the same or strongly correlated.In first place, we have  − 1 for RF, and  − 2 for FlexCoDE.These two features are strongly correlated with Pearson correlation of 0.98.In the second place, the colour  −  appears for both methods.These features ( −1,  −2, and  − ) measure the overall shape of the SED over the wide wavelength range.Beyond the second place, there is no interpretation agreement between the two models.While the first narrow-band (0378 − ) appears in 8th place for Random Forest, this same colour is the third most influential feature for FlexCoDE.This difference could be explained by the correlations between some of the colours, which can mislead the interpretation of the feature importances and does not necessarily mean that less important features are not relevant.Interestingly, not only do narrow-band colours have higher importances for FlexCoDE but also the overall importances are more equally distributed than RF.The importances shown in this figure might also be closely related to the relative depths of each band (see Table 1) but further investigation would be needed.Note that the interpretation of these importances is limited to understanding what features had the most influence in these particular prediction models and no physical or general conclusions can be taken from those.Nevertheless, there are good indications of the narrow-band importance in predicting more accurate photo-zs, as suggested by the FlexCoDE importances estimates and, especially, considering the CDE loss estimates discussed in the last section.

QUCATS: THE QUASAR CATALOGUE FOR S-PLUS
The Quasar Catalogue for S-PLUS (QuCatS) provided with this paper contains 645 980 quasar candidates with classification probabilities above 80% up to the photometric depth of  band ( < 21.3) and good photometry quality in the detection image (SEX_FLAGS_DET = 0) that comes from SExtractor.For the selection of quasar candidates, we use the star/quasar/galaxy classification probabilities from Nakazono et al. 2021 that were obtained with a Random Forest fit on S-PLUS and WISE magnitudes, and morphological parameters extracted from the S-PLUS images.Sources observed in the CCD borders that were observed in multiple fields were removed from this catalogue via an internal cross-match in RA and Dec within 1 ′′ .Our catalogue covers a total of 1414 fields from S-PLUS DR4, totalling approximately 3 000 deg2 of the southern sky.This value-added catalogue is downloadable at https://splus.cloud/files/QuCatS_ Nakazono_and_Valenca_2024.csv and its columns are described in Table 5.We do not make any cuts in photometric errors for the catalogue we provide in this paper but we advise users to make proper cuts for their specific science cases.Below, we show an example of an ADQL query 2 .The selection Table 5. Description of the information provided in QuCatS.All columns here described are type float32.We also include useful information from the main survey (ID, Field, RA, DEC, PStotal magnitudes and their corresponding errors) and the classification value-added catalogues (PROB_QSO, PROB_STAR, PROB_GAL, model_flag).These extra columns are described in https://splus.cloud/documentation/DR4.criteria can be easily modified in this query to retrieve more (or less) sources, where [columns] 3 are defined by the user.We strongly recommend that a "WHERE det.Field = [field] 4 " condition is added to this query in order to run one query per field, as big queries can overload the server.Note that one can relax the selection criteria and retrieve more quasar candidates from the S-PLUS data base, as we derived photo-z estimations for all objects in each field regardless of their classification.Likewise, one can be more conservative with the selection criteria in order to improve the purity of the selected sample.

Column name
SELECT [columns] from dr4_vacs.dr4_qso_photozas z JOIN dr4_vacs.dr4_star_galaxy_quasaras sqg on sqg.id = z.idJOIN dr4_dual.dr4_dual_ras r on r.id = z.idJOIN dr4_dual.dr4_detectionas det on det.id = z.idWHERE sqg.PROB_QSO >= 0.8 and r.r_PStotal < 21.3 and det.SEX_FLAGS_DET = 0 Density maps for the selected quasar candidates show a multimodal distribution of photo-zs for BMDN and FlexCoDE, while RF deliver 3 List of columns can be found at https://splus.cloud/catalogtools/ tap. 4 List of all unique field names and their central position can be obtained through https://splus.cloud/files/documentation/iDR4/tabelas/iDR4_pointings.csv a smoother single-peak distribution (Fig. 9).In Fig. 10 we show the photo-z distribution for the quasar candidates in the catalogue provided in this paper with probability of being a quasar down to 0.8, 0.9, and 0.95.We can see that as we restrict our sample to more high-confident quasars, there is an increasing convergence among the photo-z distributions estimated by the different methods considered in this work.For high-confident quasar candidates (PROB_QSO > 0.95), the distributions have medians at  phot ∼ 1.55.
Confirming which method is generalising better for unseen data is not trivial due to the intrinsic biases of using spectroscopic datasets to train the models.We advise users to use the information we provide in the catalogue with extra caution for the candidates that extrapolate the colour space distribution of the spectroscopic sample (for instance, see Fig. 11).Although we recommend FlexCoDE as it presented the lowest CDE loss, some science cases will demand certain requirements over specific magnitude/redshift ranges.For those, the performances per bin of magnitude , colour  −  and spectroscopic redshift allow for a more precise decision on which method to use (Fig. 12).Note that performances are poorer for lower ( ≲ 0.5) and higher ( ≳ 4.2) probably due to the lack of objects in this range (see Fig. 2).Performances for brighter ( ≲ 17.5) or bluer/redder (| −  | ≳ 0.5) objects are also affected for the same reason (see Fig. 11).The performance decrease as magnitude increases is most likely related to increasing photometric uncertainties.

DISCUSSION AND CONCLUSIONS
Up to this point, S-PLUS is the Southern hemisphere survey that provided data within the largest area (∼3000 square degrees in DR4)  and the highest number of narrow bands.This makes S-PLUS a good laboratory for empirically evaluating the narrow-band importance in tasks such as photometric redshift (photo-z) regression.In this work, we compared three independent methods (Random Forest, FlexCoDE, Bayesian Mixture Density Network) trained with S-PLUS, GALEX, and WISE colours to obtain quasar photo-zs.
A reasonable improvement in quasar photo-z predictions was observed when narrow-band colours were included in the model over training with only broad-band colours.On the other hand, narrowband information provided a small impact on the single-point estimates when GALEX and WISE colours were available.Thus narrow bands play a more important role for objects that do not have GALEX or WISE counterparts, which respectively makes 92% and 25% of the S-PLUS data.This conclusion is consistent with all tested algorithms.
As multimodal photo-z distributions are expected due to the colourredshift degeneracy, it is more appropriate to use the whole PDF information instead of single-point estimate summaries when evaluating performances.The models trained with narrow-band provide sharper distributions than the broad-band models, suggesting that the prediction uncertainties are not being properly accounted.This is confirmed through the PIT distributions, which show that the estimated PDFs are underdispersed for both FlexCoDE and BMDN (in this work, we do not analyse RF's PDFs) and a recalibration procedure is still needed.Nevertheless, for the scope of this paper, we have good indications that narrow-band photometry is helping to provide more accurate photo-zs, based on the CDE loss (which accounts for the PDF and not only the single-point estimates) and feature importance estimates for the tree-based models (especially from FlexCoDE).
Within this paper, we provide a value-added catalogue of photometric redshifts from all tested methods for 1414 fields (∼3 000 deg 2 ) of S-PLUS DR4.We broadly recommend using the narrowband FlexCoDE model as it presents the best CDE loss.
We conclude with a comment about the next-generation survey, the Rubin Observatory's LSST (Ivezić et al. 2019).LSST will obtain broad-band photometry for billions of galaxies and millions of quasars.Given our analysis here, photometric redshift performance for these objects could be improved in two principal ways.First, the wavelength coverage could be extended over that accessible to Rubin Observatory using near-infrared space-based survey facilities such as Euclid (Laureĳs et al. 2011) and Roman Space Telescope (Spergel et al. 2015).Alternatively, narrow-band photometry could be obtained after the initial 10-year broad-band survey is completed using the same telescope and camera as for LSST, but with an updated filter set as discussed in e.g.Yoachim et al. (2019) and Kahn et al. (2019).The S-PLUS project, including the T80-South robotic telescope and the S-PLUS scientific survey, was founded as a partnership between FAPESP, the Observatório Nacional (ON), the Federal University of Sergipe (UFS), and the Federal University of Santa Catarina (UFSC), with important financial and practical contributions from

Figure 1 .
Figure 1.Spectral energy distribution of five quasars with spectroscopic redshifts of 1.639, 2.238, 2.516, 3.207, and 4.435 (top to bottom)  showing an emission line being detected at an S-PLUS narrow band.In gray, we plot the spectra observed by the SDSS programs, with the identified emission lines being plotted in vertical dashed lines.The coloured triangles (squares) represent the flux measured in the S-PLUS narrow (broad) bands, from left to right: , 0378, 0395, 0410, 0430, , 0515, , 0660, , 0861, .The filter responses are shown accordingly.Empty triangles/squares represent the error in AB magnitude above 0.5 or non-detection at that particular band.

Figure 2 .
Figure 2. Distributions of spectroscopic redshift for the training and testing samples.
, providing the best Conditional Density Estimation (CDE) loss (see Section 3.5 for details) among several other methods.It is implemented in R (Izbicki & Pospisil 2019, the implementation used in this work), and Python (Pospisil, Dalmasso & Inacio 2019).

Figure 3 .
Figure 3. Architecture of the Bayesian Mixture Density Network.The input layer is followed by three blocks of DenseVariational and BachNormalization layers and a Dense layer.The numbers represent the number of neurons in that layer.The output layer (a MixtureNormal, in red) returns 7 weights (), means () and standard deviations ().

Figure 6 .
Figure 6.Density estimations  () and single-point estimates from FlexCoDE and BMDN trained with broad+GALEX+WISE (green) and broad+GALEX+WISE+narrow (pink) for the photometric redshift of 9 quasars sampled from the testing set.The spectroscopic redshift is pointed out with a triangle marker.The single-point estimates are shown with dashed vertical lines.

Figure 7 .
Figure 7. PIT distribution for (a) FlexCoDE trained without narrow bands (broad+GALEX+WISE) and with narrow bands (broad+GALEX+WISE+narrow); (b) BMDN trained without narrow bands and with narrow bands.The PDFs are considered well-calibrated if the distribution of PIT is described by a Uniform distribution.

Figure 8 .
Figure 8.Average estimated importances and standard deviation for (a) Random Forest and (b) FlexCoDE.

Figure 9 .
Figure 9. Density maps of photometrically selected quasar in S-PLUS DR4 with classification probabilities above 80% (Nakazono et al. 2021) and photometric redshift estimated with RF, BMDN, FlexCoDE, and the average of these three methods.Only sources with good photometry flag (SEX_FLAGS_DET = 0) and  < 21.3 are plotted.The colour map indicates the number of quasar candidates in bins of 0.1 for both apparent magnitude  and photo-z.In the upper panels, we show the number of quasar candidates per squared degree for 0.1 bins of photo-z.

Figure 10 .Figure 11 .
Figure10.Distribution of photometric redshift for the sources that are photometrically classified as quasars in S-PLUS with a probability above 0.8, 0.9, and 0.95.Only sources with good photometry flag (SEX_FLAGS_DET = 0) and  < 21.3 are plotted.We can note that different estimation methods are leading to more alike photo-z distributions for high-confident quasar candidates (PROB_QSO > 0.95).The median of each distribution is given in the legend.

Table 1 .
Schlafly et al. 20197 characterization of GALEX (5;Morrissey et al. 2007), S-PLUS (peak of the magnitude distribution at S/N > 3; Almeida-Fernandes et al. 2022) and WISE (50% completeness at S/N > 5;Schlafly et al. 2019).The effective wavelengths and widths are given in Angstroms, and the depths are given in AB magnitudes.

Table 2 .
Performances from cross-validation with Random Forest algorithm in terms of the root-mean-squared error ( RMSE ), normalized median absolute deviation ( NMAD ), bias (), and the fraction of sources with residual above 0.15 ( 0.15 ) and 0.3 ( 0.3 ).Standard deviations are typically of order 10 −3 .

Table 4 .
Final performances calculated on the test set for Random Forest, FlexCoDE and Bayesian Mixture Density Network.The last two rows correspond to the performance when we average the single-point estimates from these three models.The best metric values are highlighted in bold.