The Scavenger Hunt for Quasar Samples to Be Used as Cosmological Tools

Although the $\Lambda$ Cold Dark Matter model is the most accredited cosmological model, information at high redshifts ($z$) between type Ia supernovae ($z=2.26$) and the Cosmic Microwave Background ($z=1100$) is crucial to validate this model further. To this end, we have discovered a sample of 1132 quasars up to $z=7.54$ exhibiting a reduced intrinsic dispersion of the relation between ultraviolet and X-ray fluxes, $\delta_\mathrm{F}=0.22$ vs. $\delta_\mathrm{F}=0.29$ ($24\%$ less), than the original sample. This gold sample, once we correct the luminosities for selection biases and redshift evolution, enables us to determine the matter density parameter $\Omega_M$ with a precision of 0.09. Unprecedentedly, this quasar sample is the only one that, as a standalone cosmological probe, yields such tight constraints on $\Omega_M$ while being drawn from the same parent population of the initial sample.


Introduction
Recently, improved precision in measuring cosmological parameters has exposed tantalizing discrepancies within the widely accepted Λ Cold Dark Matter (CDM) model.This model describes the Universe relying on a CDM and dark energy components, where the dark energy is a cosmological constant (Λ), as required by the current accelerated expansion of the Universe [1,2].This model, with its enigmatic dark energy and cold dark matter, has served us well, explaining phenomena like the Cosmic Microwave Background (CMB) [3] and the accelerated expansion of the Universe proved by type Ia supernovae (SNe Ia).Despite its advantages, theoretical flaws still need to be understood.This is the case of the cosmological constant problem [4], which is the tension between the predicted and observed values of Λ, the nature of dark energy and its origin, and the fine-tuning problem, which derives from the fact that the current values of the matter density (Ω M ) and the dark energy density (Ω Λ ) are of the same order, whereas this is not expected due to their different evolution in time.In addition to these issues, recent measurements have highlighted the so-called Hubble constant (H 0 ) tension.This is the discrepancy between the value of H 0 measured locally from SNe Ia and Cepheids, which is H 0 = 73.04 ± 1.04 km s −1 Mpc −1 [5], and the value of H 0 extrapolated from the Planck data on the CMB within a flat ΛCDM model, H 0 = 67.4± 0.5 km s −1 Mpc −1 [3].The difference between these two measurements ranges between 4.4 and 6σ, according to the samples investigated [6][7][8].However, the maximum redshift reached by SNe Ia observations is z = 2.26 [9], while the CMB radiation is observed at z = 1100.Thus, it is crucial to probe the Universe in the intermediate epochs between these two to shed light on this tension, hence to confirm, alleviate, or even solve it.To this end, other probes rather than SNe Ia and CMB have already been investigated.These analyses have provided even a more complicated context: cosmic chronometers show a preference for the H 0 value derived from the CMB [10], time delay and strong lensing from Quasars (QSOs) favor the H 0 from SNe Ia [11], while QSOs [12], the Tip of the Red-Giant Branch [13], and Gamma-Ray Bursts (GRBs) [14,15] hint at an intermediate value of H 0 between the one of the CMB and the H 0 of SNe Ia.For the case of QSOs and GRBs, the value of H 0 obtained depends on several factors, such as if they are calibrated or not with SNe Ia, if they are fitted jointly with SNe Ia, and if their luminosities have been corrected for the redshift evolution and selection biases.Moreover, we here note that, when also SNe Ia are used, H 0 is a parameter degenerate with the absolute magnitude of SNe Ia.Nevertheless, the values of H 0 obtained, when QSOs or GRBs are considered, show a trend toward an intermediate value between the CMB and SNe Ia [12,14].To solve this intriguing puzzle and test if the flat ΛCDM model still represents the most suitable description of the Universe, reliable and powerful cosmological probes at redshift between z = 2.26 and z = 1100 are required.To crack this perplexing puzzle and put the ΛCDM model to the ultimate test, we need to voyage through uncharted territories: we are embarking on a mission to explore the Universe between the epochs of SNe Ia and the CMB.To date, the best candidates for this purpose are GRBs and QSOs.
In this framework, QSOs have recently attracted more and more interest in the cosmological community [12,[16][17][18][19][20][21][22][23][24], since they are observed up to z = 7.64 [25], at redshifts much higher than the maximum redshift of SNe Ia observations, z = 2.26 [9].The method to standardize QSOs as cosmological candles is based on the Risaliti-Lusso (RL) relation between the logarithms of the Ultraviolet (UV) luminosity at 2500 Å (L UV ) and the X-ray luminosity at 2 KeV (L X ).The RL relation reads as log 10 L X = γ log 10 L UV + β, where the slope γ and the intercept β usually have values of γ∼0.6 and β∼8.This empirical relation has been validated with several QSO samples [26][27][28][29][30][31][32][33][34][35] and it has been turned into a cosmological tool via a careful selection of the QSO sources aimed at removing observational biases [17,19,34,[36][37][38][39][40].To achieve this standard set starting from the initial QSO sample, several QSOs have been discarded through the investigation of different features.This procedure allows the sample to present well-defined properties, not to be hampered by a low signal-to-noise ratio, and not to be severely affected by extinction, UV reddening, and contamination of the host galaxy (see Section 2).The relation between UV and X-ray luminosities is also theoretically supported by the most accredited QSO model in which an accretion disk powers the central supermassive black hole converting mass into energy.In this scenario, the UV emission of the accretion disk is then reprocessed in X-rays from an external region of relativistic electrons via the inverse Compton effect.Nevertheless, this mechanism still needs to be fully understood to explain the stability of the X-ray emission, which is not expected since the electrons should cool down falling on the central region.Thus, this stability requires an efficient energy transfer between central and external regions, whose origin is yet to be unveiled [37,41].Moreover, Ref. [21] has confirmed the reliability of the RL relation in cosmology proving that this relation is completely intrinsic to the QSO physics and not induced by selection effects or redshift evolution.This is a crucial turning point for the reliable application of QSOs as cosmological tools.From now on, the notation "RL relation" refers to the relation obtained by the Risaliti-Lusso group, while the notation "RL correlation" is used when the relation is applied for fits, in bins, in the flux-flux space, and in the luminosity-luminosity space.
In common practice, QSOs are used jointly with other probes since the intrinsic dispersion, δ, of the RL relation (δ∼0.23 [17,40]) still limits their power in constraining cosmological parameters, compared to the precision of other probes, such as SNe Ia.We here stress that we use the term "intrinsic dispersion" to refer to the additional fit parameter of the RL relation which is implemented in the fitting likelihood.This allows for a spread around the ideal RL best fit.For this reason, we here focus on the scavenger hunt of a subsample of QSOs, a "gold" sample, which presents the optimal compromise between reduced intrinsic dispersion and a sufficient number of sources to be used as a standalone probe and constrain Ω M with unprecedented precision in the QSO realm.Similar efforts have already been made in the GRB domain [14,15,[42][43][44][45], leading to the definition of the "Platinum" GRB sample, which has been used in several cosmological analyses [14,45].Indeed, to choose a standard candle it is crucial to identify either the morphological properties of the lightcurves or the spectral features of the objects investigated.In the case of GRBs, the morphological feature that drives a standard candle is the plateau emission with peculiar characteristics, such as not be a steep plateau (with an angle < 41% and to have at least 5 data points at the beginning of the lightcurve and not to have flares or gaps inside the plateau region) (see, e.g., [46][47][48] for details).Regarding QSOs, the feature stressed here is the fulfillment of the RL relation.
We here also stress the reason why the hunt for a gold cosmological sample, not only for QSOs but also for GRBs, has assumed such a relevant role in the cosmological community.Indeed, as anticipated, due to the intrinsic dispersion of the RL relation, QSOs are not able alone to constrain cosmological parameters with precision.As a consequence, once they are combined with other probes, such as SNe Ia, they are not the "leading" probe, which means that the information that dominates is the one from the other probes.Besides that, QSOs could be extremely important to add information at z > 1.So, to make QSOs powerful cosmological probes, we need first to determine a QSO sample that is itself capable of constraining cosmological parameters when used alone.Only after that, we can join this gold sample with other probes.In this case, QSOs would contribute to the determination of cosmological parameters by adding their piece of information and changing, confirming, or eventually even improving the constraints given by the other probes.In view of the hunt for a gold QSO sample, we here present an unprecedented sample of QSOs reaching the maximum redshifts up to z = 7.54, which can be applied as a standalone cosmological probe to constrain Ω M .We here point out that our main aim is to develop a selection procedure that can be reliably applied to select a QSO sample that can be used alone to constrain cosmological parameters.Indeed, our purpose is not to investigate cosmological tensions, as for example the H 0 tension, or to test cosmological models alternative to the standard flat ΛCDM.This is the reason why we employ our final selected QSO sample to fit the specific case of a flat ΛCDM model with H 0 fixed and only Ω M free to vary.This analysis shows that this sample constrains Ω M with a precision unprecedentedly reached with only QSOs.
This manuscript is organized as follows.Section 2 describes the initial QSO sample.Section 3 details the correction of luminosities for selection biases and redshift evolution, the selection of the final QSO sub-samples, and the cosmological fitting method.In Section 4, we outline our results, and in Section 5, we draw conclusions.Appendix A discusses the different binning approaches investigated and Appendix B details the Efron and Petrosian method.

The Data Sample
The initial QSO data set for our analysis is the most recent one released for cosmological applications [40].It counts 2421 sources ranging between z = 0.009 and z = 7.54 [49] collected from eight different catalogs in the literature [38,50,51] and archives [52][53][54][55], with the addition of a sub-sample of low-redshift QSOs that present UV observations from the International Ultraviolet Explorer and X-ray data from archives.To obtain this QSO sample suitable for cosmological analyses, as many as possible observational biases have been carefully inspected and removed [17,34,36,40,50].We here briefly describe the steps of this selection.First, only measurements with a sufficient signal-to-noise ratio (S/N ≥ 1) are retained.Then, QSOs that manifest the presence of extinction (i.e., E(B − V) > 0.1) are removed to account for UV reddening and contamination of the host galaxy.The contribution of absorption in X-ray is also removed by imposing Γ X + ∆Γ X ≥ 1.7 and Γ X ≤ 2.8 if z < 4 and Γ X ≥ 1.7 if z ≥ 4, with Γ X and ∆Γ X being the photon index and its uncertainty, where the photon index is the coefficient of the power-law that describes the spectrum in X-ray.Eventually, the final sample is corrected for the Malmquist bias effect.Indeed, the Malmquist bias effect states that only larger fluxes can be observed at larger distances; thus, low fluxes are prevented from being seen due to the detector flux limitations, and this creates incomplete samples.This effect is overcome by requiring logF X,exp − logF min ≥ F , where F X,exp is the X-ray flux computed from the flux in UV by imposing the RL relation and assuming Ω M = 0.3 and H 0 = 70 km s −1 Mpc −1 in a flat ΛCDM model.Even though this correction requires the assumption of a specific cosmological model, simulations and mock samples of QSOs have been employed to prove that results are not affected by this choice, as explained in [40].F min is the minimum observable flux computed for each source from the time of observation of the chargecoupled device [34,56].F is the threshold value, which is fixed to F = 0.9 for QSOs derived from the cross-match of the Sloan Digital Sky Survey Data Release 14 (SDSS DR 14) with 4XMM Newton or with XXL, and to F = 0.5 for the ones with measurements from the SDSS DR 14 and Chandra.To reduce the effects of the X-ray variability, if a source has more than one X-ray observation after this selection, these observations are averaged.We here notice that the sources that have multiple X-ray observations in the QSO sample of [40] are only 289 [57], which is 12% of the total sample.Thus, for the remaining 88% of the sample, we cannot reduce the effect of the variability on the dispersion of the RL relation if we consider that the mechanisms responsible for this variability are still not completely understood.Nevertheless, it is possible to estimate the contribution of the X-ray variability to the observed scatter, as reported in [57] and discussed in Section 4.5.In this analysis we rely on a selection of a sample already studied in the literature and how the impact of these sources can affect the whole analysis is an interesting subject, but it goes outside the scope of the current paper.In our work, we start from this final sample of 2421 sources without any additional selection, such as the cut at redshift z = 0.7 previously used in some works [19], to avoid any possible induced bias due to the reduction in the redshift of the sample [12,21].

Methods
Since the 2421 QSOs sample is not yet ready for cosmological use, we here detail a method to find a golden sample, the best, optimal sample useful for cosmological studies.

Selection of the QSO Final Samples
To transform QSOs into powerful cosmological probes, we meticulously outline the procedure employed to define our final QSO sub-samples.The steps of this procedure are detailed below and are justified by physical and theoretical requirements.Indeed, these are necessary steps to define a suitable technique that can be reliably applied to slim the QSO sample aiming at turning them into standalone cosmological probes.To clearly detail and visually show our logical flow, we present the graphical representation of the following steps, regarding our selection methodology in Figure 1.We accept the bin as a final We move the whole bin to the untouched set We use the Huber regressor to obtain outlier free bin.
We combine all bins creating our final sample "without source".
We combine the sample "without sources" with the "untouched set" creating the sample "with sources".

1.
We have first divided the initial QSO sample into bins of redshift to fit a linear relation between the logarithms of fluxes in each redshift bin.The binning must be chosen to verify a specific condition that can be derived from the RL relation [17,35,36,40].The RL relation, log 10 L X = γ log 10 L UV + β, can be written in terms of fluxes as where a is the notation we use to refer to the slope in the flux-flux plane, and D L is the luminosity distance in units of cm.Equation ( 1) can be approximated as a linear relation between log 10 F X and log 10 F UV only if the contribution of log 10 D L is negligible compared to the other terms in Equation (1).Thus, in this case, the linear relation in fluxes represents a proxy of the RL relation in luminosities, with a different intercept, in the form log 10 F X = a log 10 F UV + b and with intrinsic dispersion δ F .More specifically, if we consider Equation (1) in a redshift bin, the contribution of the distance is negligible if the range of values of log 10 D L within the chosen redshift interval is smaller than the intrinsic dispersion of the relation in the same bin.This is the condition that must be fulfilled when choosing how to divide the QSO sample into redshift bins.Nevertheless, in addition to this requirement, we need to fine-tune our choice to ensure enough sources (we require a minimum of 10), at least in the majority of bins, to reliably perform the fit in each of them.The specific choice of 10 sources is arbitrary.Indeed, we could apply a different threshold, which allows sufficient statistics to perform the fit.We have actually performed our analysis also changing the minimum number to 4, 5, 6, and 10, without any change in our results.We here also notice that to reduce the scatter of the relation in each bin, the condition of narrow redshift bins is very relevant, as the smaller the bins, the smaller the difference in log 10 D L is.Hence, we can reduce the intrinsic dispersion up to the limit imposed by log 10 D L .Following all these prescriptions, we have defined our optimal division into redshift bins in terms of log 10 (1/(1 + ∆z)) with ∆z = 0.042 (see the yellow box in Figure 1).We have adopted the division as log 10 (1/(1 + z)), which is a natural choice for the division into redshift since it would retain the same division in volume.This way, we can keep the bin constant, and we do not need to derive arbitrary bins.This is an improvement of the method of bin division in [40].With this division in redshift, we obtain 32 bins with at least 10 sources (see Table 1), which is the threshold we require to guarantee sufficient statistics for the fit.We hereby stress that binning the data into redshift intervals is necessary to use fluxes instead of luminosities.This is a crucial point as it enables us to perform a circularity-free analysis.Indeed, fluxes are measured quantities that do not require any cosmological assumption, unlike luminosities.As a consequence, the use of fluxes in the selection of the sample guarantees that our cosmological results are not induced by any a priori cosmological assumption.It is true that binning leads to the reduction of the sample size (in each bin compared to the total sample size), and therefore, the estimates in each bin might be less accurate.However, in our case, the binning shows that the slope of the flux-flux RL correlation in each bin remains unchanged (see Figure 2), and is compatible with the slope reported in [40].In addition, binning is often used when it is necessary to highlight features that would otherwise be concealed when noisy data are combined altogether.In this analysis, the binning of the adopted to avoid the circularity problem (see also [58] for a discussion on the importance and reliability of the binning method).This is because the approximate RL correlation for fluxes, which does not depend on the cosmological parameters, holds only within bins of a limited length of redshift and hence of the distance luminosity.This relation within each bin allows us to highlight which QSO sources should be removed.Moreover, we have detailed in Section 5 that our analysis with the binning gives compatible results with the unbinned data (see [59] for comparison).We also further investigate different choices for the division into bins of the initial sample in Sections 3.1.2-3.1.4and their impact on the cosmological results in Appendix A.  1).Mean redshift < z > of each redshift bin with at least 10 sources along with the number of sources N, and the best-fit values of the slope a and the intrinsic dispersion δ F with their 1σ uncertainty after the removal of outliers.

2.
Once we have divided the redshift bins, we fit in each bin that presents at least 10 sources a linear relation between log 10 F X and log 10 F UV .This fit is performed using the Kelly method [60], which accounts for the uncertainties in both quantities and also for the intrinsic dispersion of the correlation.We have also imposed uniform priors in a wide range of values for the free parameters of the fit: the slope, the intercept, and the intrinsic dispersion.To verify that the condition described at Point (1) is satisfied, the best-fit value obtained with the Kelly method for the intrinsic dispersion is compared to the maximum difference of log 10 D L for the sources in the investigated redshift bin.This difference is computed by assuming a flat ΛCDM model.We here notice that the assumption of a specific cosmological model for this computation does not affect the result since we are considering a difference between two luminosity distances.We have retained unmodified sources in the redshift bins that do not provide enough statistics (less than 10 QSOs) to perform a reliable fit.From now on, we denote these sources with the notation "untouched".Also, we have distinguished two cases, one in which we do not include these sources and another one in which we have added them to the final selected sample obtained after Point (3).The strategy here is to balance and compromise among the smallest bin so that ∆log 10 D L < δ, but still sufficiently large so that the number of sources is at least 10 or more.

3.
At this stage, as the presence of outliers can decrease the performance and accuracy of least-squared-loss error-based regression, we have employed the consolidated statistical technique of the Huber algorithm [61][62][63] to reduce the intrinsic dispersion in each bin considered.The Huber regressor is indeed a method for estimating the parameters of a model, in this case the F X − F UV relation, to detect the outliers and weigh them less in the evaluation of the best-fit parameters of the fitted model.We are indeed aware that sources more scattered around this relation hamper significantly the finding of the most suitable sample with the smallest intrinsic dispersion.Thus, compared to traditional fitting procedures, such as the D'Agostini [64] or the Kelly [60] methods, the Huber regression identifies outliers, which can be caused, for example, by errors or problems in the measurements, and recognizes the actual best-fit based on the inliers.For these reasons, this technique is widely applied for robust regression problems.The Huber regressor has the advantage of not being heavily influenced by the outliers, while not completely ignoring them.This allows us to estimate the actual slope and intercept of the relation, not altered by outliers, and contemporaneously to identify the sources that are outliers of the model.Hence, we discard these sources from the QSO sample in each redshift bin.In order to quantitatively evaluate the Huber algorithm's numerical gain against the traditional fitting one, we have also compared the results obtained with the Huber regressor with those derived from the traditional sigma-clipping selection technique.This comparison is detailed in Section 4.3.
After this selection, we have also checked in each bin the following criteria: the null hypothesis that the populations of both UV and X-ray fluxes are drawn from the initial ones in the bin considered must not be rejected with p-value > 5% according to the Anderson two-sample test and if the distribution of the residuals about the best-fit line is Gaussian according to the Anderson-Darling normality test with an acceptance significance level of 5% (see, e.g., [65] for the Gaussianity discussion).The Anderson-Darling test for normality determines whether a data sample is drawn from the Gaussian distribution, and it is commonly applied in the literature (e.g., [43] in astrophysics and [66,67] in statistics).An important property of this test is that it can identify any small deviation from normality.We refer to [68] for a detailed description of the features of this test and its application to cosmological likelihoods.The Anderson-Darling two-sample test instead allows us to verify if the selected sample is still drawn from the original one.This guarantees that we are neither introducing biases nor significantly changing the physical properties of the initial sample when selecting the final sample.We here also stress that the Anderson twosample test is always fulfilled at a statistical level >25%.Table 1 reports the mean value of z (< z >) for each redshift bin with at least 10 initial sources, the number of sources retained, and the corresponding best-fit values for the slope and the intrinsic scatter of the linear relation.A visual representation of the trend of the best-fit values of the slope with the average redshift of each bin is also provided in Figure 2. To showcase the Huber regressor's advantage in each bin and how effectively it removes the outliers, in Figure 3, we present in green the selected sample and in red the sources identified as outliers.The two bins investigated on the left and right panels of this figure are the second most populated one and the second least populated one, respectively.4.
As anticipated, we have finally defined two ultimate samples: one with only the sources retained through the steps detailed above and another obtained by combining the sources retained in each bin with the unmodified sources of the bins without enough statistics.This way, we have generated the final selected QSO samples composed of 1065 and 1132 sources, respectively.We here anticipate that, among all the binning approaches investigated in this work, we choose as the best one the one that leads to the best precision on Ω M for both samples, with and without untouched sources.The sample obtained with this best method is the one referred to as the "gold sample".By using these final QSO samples, we have also fitted the linear relation in fluxes, as shown in the left panels of Figures 4 and 5.As a further step, we have transformed the relation from fluxes to luminosities (see Figures 6 and 7) to check if, for both samples, the slope of the luminosity-luminosity correlation obtained from our analysis is consistent with the slope of the RL relation corrected for the redshift evolution.More specifically, we have computed from the fluxes the corresponding luminosities, and we have fitted a linear relation among them with the following form:    The results of these analyses are shown in the right panels of Figures 4 and 5 and Figures 8 and 9.We also note that in the right panels of Figures 4 and 5 the RL correlation in the luminosity-luminosity space is presented by assuming a certain cosmological model and corrected for selection biases and redshift evolution.Therefore, these plots are only for the purpose of showing the reached tighten relation, and they are presented as an example assuming a flat ΛCDM model with Ω M = 0.3 and H 0 = 70 km s −1 Mpc −1 .The code and a comprehensive technical description of this method can be accessed from the Wolfram Mathematica Notebook Archive [69].Before performing any cosmological analysis with our selected samples, we have also demonstrated through simulations that we are able to retrieve any assumed input cosmology for a mock QSO sample with a similar redshift and flux distribution to our final sample.Indeed, we have generated 1065 and 1132 mock data, respectively, for the two final samples, with distributions of redshift (see Figure 10), fluxes, and uncertainties on fluxes drawn from the corresponding best-fit distributions of our final samples of observed data composed of 1065 and 1132 sources.Also, we have assumed a priori a cosmological model to compute the luminosities, and we have also fitted the cosmological parameters of the investigated cosmological model.Specifically, we have investigated different assumptions for the cosmological model: flat ΛCDM models with Ω M = 0.3 and H 0 = 70 km s −1 Mpc −1 , with Ω M = 0.1 and H 0 = 80 km s −1 Mpc −1 , with Ω M = 0.5 and H 0 = 65 km s −1 Mpc −1 , and with Ω M = 0.8 and H 0 = 60 km s −1 Mpc −1 .In all these cases, fixing the value of H 0 to the assumed one and applying the redshift correction with "varying evolution", we have recovered the assumed value of Ω M within 1σ.We here clarify that we use the notation of "outliers" not in a strict statistical sense but rather to refer to the QSOs that show more discrepancy from the RL relation line.With the above-described division of the sample into bins (see also Table 1), we can prove that there is not a particular trend or significant behavior of the slope a as a function of the redshift.We here show that the trend of the slope values corresponding to the average redshift fluctuates around an average value of a = 0.60 (see Figure 2).Some fluctuations are visible, but it is expected as the sample size is not equally divided according to the number of sources.

The Parameter ϵ in the Huber Procedure
We can also take advantage of the parameter ϵ, which is a free parameter in the Huber regression method in the range [1, inf) which can be arbitrarily fixed.This parameter controls the number of samples that should be classified as outliers: the smaller ϵ is, the more robust the Huber regression is to define outliers.We have indeed performed our selection by trying different values for ϵ, and we have obtained compatible results in all the attempts.Thus, we have identified the value of ϵ = 1.2 as the one that leads to the cosmological result with the smallest uncertainty on Ω M and yet a considerable statistical sample.The impact of the choice of ϵ on our analysis is also shown in Sections 3.1.2-3.1.4.

The Impact of the Binning on the Data Analysis: Bins Wider in Redshift
As anticipated, to further investigate the impact of our choice for the division in redshift bins (i.e., log 10 (1/(1 + ∆z)) with ∆z = 0.042) (see, e.g., [70] for a discussion on the binning), we have also selected the QSO sample by using three other choices for the division in bins, different from the one detailed above, which are described here and in Sections 3.1.3 and 3 We here start from the method of dividing the sources in bins whose width depends on the redshift (see the blue box in Figure 1).Specifically, we have divided the original sample according to the following prescription: and ∆log 10 z = 0.06 if z > 4. The choice of increasing the width of the bins for higher redshifts is justified by the trend of the luminosity distance D L (z) Indeed, the trend of D L (z) is much steeper at low z compared to the one at higher z, where the function D L (z) flattens; thus, we need to impose narrower ∆log 10 z at low redshifts to verify the condition that the difference ∆log 10 D L is smaller than the intrinsic dispersion in the corresponding bin.The specific values of ∆log 10 z in each redshift range have been chosen to guarantee the fulfillment of all the required criteria detailed above: a minimum number of 10 sources that fulfill the requirement conditions explained above, a range of values of log 10 z within the chosen redshift interval so that ∆log 10 D L is smaller than the intrinsic dispersion of the relation in the same bin, and the Anderson-Darling two-sample test that is passed at a statistical level of at least 5% to warrant that the selected sample is drawn from the original one in the investigated bin.By applying this division into bins, we generate 29 bins with at least 10 sources.As already discussed in the previous section, also in this case the value ϵ = 1.2 for the Huber regressor proved to be the one, among the several values of ϵ tested, that leads to the cosmological results with smaller uncertainties on Ω M .Furthermore, the results obtained by applying different values of ϵ are completely compatible with each other.Ultimately, with this procedure, we identify two final samples, as described in Point (4) of Section 3.1: one comprising 1084 QSOs, which does not contain the sources in the bins in which some of the conditions are not fulfilled, and one comprising 1125, which also includes these sources (see Table 2).This analysis, although a good alternative to the previous method, does not allow us to reach the golden sample if we consider the precision on Ω M obtained for both samples, with and without untouched sources.
Table 2. Comparison of results obtained from different selection approaches.The first column details the method applied to divide the initial QSO sample in bins, the second column specifies if the sample includes or not the sources in bins that do not fulfill the criteria required by our trimming analysis (see Sections 3.1 and 3.1.2-3.1.4for further details).The third and fourth columns report, respectively, the number of sources (N) in the considered sample after the removal of outliers and the estimated Ω M with its 1σ uncertainty.The symbol "*" identifies the cases in which Ω M is not constrained.The parameter z-score gold is a comparison between each result and the value with the smallest uncertainty obtained in this work, which is Ω M = 0.229 ± 0.061, while z-score SNe computes the compatibility of each result with Ω M = 0.338 ± 0.018 from [71] (see Appendix A for details).To avoid the arbitrary choice of the division into bins that is intrinsically needed by the two procedures outlined above, we have also developed an approach to optimize the binning.This method is based on the maximization of the number of sources (see the orange box in Figure 1).More specifically, we start from the source at the highest redshift in the initial sample or the source at the lowest redshift (we have investigated both cases).We define the edge of the bins so that the maximum number of sources is included in each bin and still the same criteria of the previous procedure are fulfilled.The maximum difference in D L (z) in the bin is smaller than the intrinsic dispersion in the same bin.Also, the selected sample is still drawn from the original one in the bin according to the Anderson-Darling two-sample test with a threshold for the p-value of 5%.With this procedure, we have found that the bin with fewer sources is the first, with 12 QSOs, while the most populated bin has 567 sources.After the first bin has been identified, the same procedure is repeated to create another bin adjacent to it in redshift, and this algorithm continues until the size of the initial sample is reached.Following this recursive approach, this algorithm automatically creates bins to divide the QSO sample without any arbitrary choice.

Method
When we arrive at the point where almost the whole sample is divided into bins, we could face an issue that the number of data points that still do not belong to any bin would be smaller than 10 (the assumed minimum size).Thus, it is impossible for those points to represent a reliable bin.We store those sources in a separate set called "untouched".Alternatively, this could also be the case in which the number of untouched sources is larger than 10, but the condition between the ∆log 10 (D L (z)) and the intrinsic scatter cannot be fulfilled for the set of these sources.These scenarios can occur whether we start from the source at the highest redshift or from the source at the lowest redshift.Thus, in these cases, we treat these sources that remain out of the binning division just like the untouched sources in bins without sufficient statistics.As already discussed in Section 3.1, we distinguish a sample in which these QSOs are not included and another one that instead includes them.With this approach, we obtain 9 bins and 15 untouched QSOs starting from the highest redshift and 9 bins and 2 untouched sources when we start our binning procedure from the lowest redshift.We here notice that the number of bins is significantly smaller compared to the ∼30 bins used in the two above-described methods.This can be ascribed to the fact that in this case, we are not a priori imposing a specific division into bins.Still, the binning is automatically generated while requiring the fulfillment of the necessary conditions for our analysis.This causes a smaller number of bins, which are more populated in the range of intermediate redshifts between z∼0.4 and z∼3.This innovative procedure of optimization of bins will be described with more details and in-depth analyses in [72].As in the previous approaches, we have also tested different values for the ϵ parameter of the Huber regressor.Since the results proved to be completely compatible, we have chosen the values ϵ = 1.5 and ϵ = 1.6, respectively, for the methods starting from the highest redshift and the lowest redshift, as these choices guarantee the cosmological results with the smallest uncertainty on Ω M .Using this approach, we have obtained the following final QSO samples (see Table 2): 1843 and 1858 sources, respectively, not including and including the sources that do not belong to any bin, for the method starting from the highest redshift, and 1965 and 1980 QSOs, for the corresponding cases when starting from the lowest redshift.This method allows us to reach a sample of low and high redshift, which is not covered by the other method.This allows a smaller number of untouched sources.Similarly to the previous method, it does not lead to the golden sample since it does not allow us to reach the same precision as the binning division in log 10 (1/(1 + ∆z)).

The Impact of the Binning on the Data Analysis: A Bin Centered on Each Source
Furthermore, we have also developed an additional method that allows us to completely free our analysis from the possible issues of the binning approach.As a matter of fact, in this procedure, we do not actually associate each QSO with a specific bin.Indeed, we consider each of the initial 2421 QSO sources as the center of an interval (see the purple box in Figure 1).This interval is symmetric in redshift; thus, it consists of the five sources next to the central one at lower redshifts and the five sources next to the central one at higher redshifts.In this way, each interval consists of a total of 11 sources, the minimum number we require to statistically perform a reliable fit, since we have a central source and 5 on the right-end side and 5 on the left-end side of the interval.Thus, this number must be odd by construction.Then, within each of these intervals generated around each of the initial QSOs, Steps (2), (3), and (4) outlined in Section 3.1 are performed.Specifically, in each of these intervals, the condition on D L (z) and the intrinsic scatter is checked, then the Huber regressor is applied, and the Anderson-Darling two-sample test is performed to verify that the selected sample is still drawn from the initial one and the same interval.
We define a source as an inlier only if the Huber regressor selects it as an inlier in all the intervals it belongs to.As in the other approaches detailed above, we still account for the QSOs in intervals that do not satisfy the conditions of our analysis by distinguishing the two final samples in which we do not include and include these sources.we have again checked that choosing different values of the parameter ϵ for the Huber regression does not impact our analysis, leading to compatible results.Thus, we have chosen ϵ = 1.3 as the value that leads to the cosmological result with the smallest uncertainty on Ω M .Hence, we obtain the following samples, summarized in Table 2: 811 and 825 sources, respectively, when we do not include and add the sources belonging to intervals that do not fulfill our criteria.Once again, this method does not lead to the golden sample if we consider the precision reached on Ω M .
In the following, we focus on the cosmological analysis performed with the final samples obtained from the first methodology detailed, which is the one with the binning in log 10 (1/(1 + ∆z)) with ∆z = 0.042.This is indeed the approach that, among all the methodologies investigated, leads to the cosmological result with the smallest uncertainty on Ω M .Nevertheless, we also discuss and compare the results obtained by applying the other selection procedures in Appendix A.

Treatment of Redshift Evolution and Selection Biases
Since QSOs are observed up to high redshifts, we need to correct their luminosities for selection biases and evolutionary effects [73], which could, in principle, distort or induce a correlation between luminosities, thus inducing an incorrect determination of cosmological parameters [74].To apply this correction, we employ the Efron and Petrosian (EP) statistical method [75] already used in several works [76,77] for GRBs [14,24,42,74,78,79] and in the QSO realm [12,21,24,79].We use our own package to better customize it for our own analysis (see the Mathematica notebook accessed in 18 May 2023) (https://notebookarchive.org/2023-05-8b2lbrh).In this study, we apply to the obtained QSO sub-samples the procedure outlined in the above-mentioned works.In this section, we summarize the method and outcomes.
In the EP technique, we are able to determine if there is evolution among the redshift and the variables at play.With the term "evolution" we refer to the trend of a given variable, in this case, the luminosity with the redshift as the variation of this variable with the redshift.The luminosities are assumed to evolve with z according to L ′ = L (1+z) k , where L is the observed luminosity, L ′ the corresponding corrected one without evolution, and k the parameter that mimics the evolution.The L ′ are the corrected luminosities, where both the intrinsic evolutionary effects and selection biases have been removed, and this is the main reason we should use these values and not the uncorrected luminosities in the final computation for cosmological use.From now on, with the symbol ′ , we indicate the new de-evolved quantities after the correction for the evolution, not only for the luminosities but also for the parameters of the RL relation (see Equation 2).Nonetheless, the choice of the functional form as a power-law does not affect the results [21,80,81] and hence we could also parameterize the dependence on the redshift through more complex functions.Then, Kendall's τ statistic is applied to identify the k value that eliminates the evolution with the redshift.In this procedure, τ is defined as where R i is the rank defined as the number of points in the associated set of the i-source; the associated set consists of all j-points for which z j ≤ z i and L z j ≥ L min,i , with L min,i being the minimum observable luminosity (L min,i ) at that redshift.The EP method defines subsamples of the data, which are called "associated sets" and defined as containing objects, denoted with j, which should have luminosity larger than the minimum luminosity pertinent to that object and still observable according to the satellite threshold limit and the redshift should always be smaller than the redshift of a given object.In Equation ( 3), E i = 1 2 (i + 1) and V i = 1  12 (i 2 + 1) are the expectation value and variance, respectively, when the evolution with redshift has been removed.As a consequence, the correlation with redshift disappears when τ = 0, which allows us to obtain the value of k that removes the dependence.The condition |τ| > n implies that the hypothesis of no correlation is rejected at nσ level.We provide the 1σ uncertainty on the k value by imposing |τ| ≤ 1.The found value of k can now be used to determine L ′ for the total sample.L min,i is computed by requiring a limiting flux.The value of this flux threshold is chosen such that the retained sample is composed at least of 90% of the total initial sources and that it resembles the overall original distribution according to the Kolmogorov-Smirnov test [21,42,45,82].Indeed, we here stress that the above-described procedure employed to correct the luminosities is applied to both the X-ray and UV luminosities separately.This means that we obtain two different evolutionary coefficients, which are k UV and k X for L UV and L X , respectively.The results of the application of the EP method to our two selected QSO samples, which will be described in this work, are provided in Appendix B. Indeed, Ref. [21] has already proved that the initial QSO sample of [40] suffers from redshift evolution (see Figure 2 of [21]) and thus the luminosities need to be corrected through the EP method.In [83], where the same evolutionary form of (1 + z) k is used, they found k opt = 3.0 ± 0.5 and corrected the luminosity function.Thus, the new luminosity function can be representative of the observed luminosity function, but it will be constructed with the local luminosities (de-evolved luminosities), and thus, they will be rescaled by the g(z) functions.
Nevertheless, from the description of the EP method, is clear that k is obtained assuming a specific cosmological model, needed to compute the luminosities from the fluxes.Usually, the assumed model is a flat ΛCDM model with Ω M = 0.3 and H 0 = 70 km s −1 Mpc −1 .This induces the so-called "circularity problem".This problem has been completely overcome for the first time by [14] for GRBs and [12,21] for QSOs, which have analyzed the trend of k as a function of the cosmological model assumed a priori.More precisely, in these studies, k is determined not by fixing the cosmological parameters of the assumed model, but over a grid of values of the cosmological parameters (i.e., Ω M , H 0 , and also other parameters for models different from the flat ΛCDM one), leading to the determination of the functions k(Ω M ) and k(H 0 ).Due to the invariance of τ under linear transformations of data, k does not depend on H 0 .However, it shows a dependence on Ω M , and thus k(Ω M ) can be applied in the cosmological fits while leaving k free to vary along with the free cosmological parameters.Hence, the RL relation log 10 L ′ X = γ log 10 L ′ UV + β can be written in terms of the evolutionary coefficients as log 10 L X = γ log 10 This overcomes the circularity problem since we do not fix any cosmology a priori.In all our computations, we employ this method, which we refer to as"varying evolution", since it allows us to avoid the assumption of a specific value of Ω M to correct the luminosities for this effect, for details, see the Appendix B. We here note that the "varying evolution" methodology allows us to elude any degeneracy between the evolutionary coefficients, k UV and k X , and the other fitted parameters, which are γ, β, and the cosmological parameters involved.On the other hand, letting k UV and k X to vary together inside a MCMC fitting, without knowing the dependence of cosmological parameters with these other parameters, would introduce degeneration.This is the reason why we determine k UV (Ω M ) and k X (Ω M ) as a step zero by applying the EP method in a completely cosmology-independent way, and then we use these functions in the cosmological fits.This issue was also already dealt with in [24].Since there is degeneracy among the k UV (Ω M ) and k X (Ω M ), we are able to determine precisely one variable if we know the other; thus, we prefer to leave the varying evolution approach with the two functions determined before we perform the cosmological fitting.We here also notice that the "varying evolution" method can be generalized to cosmological model other than the flat ΛCDM one.In this regard, [12] show how the evolutionary coefficients in UV and X-ray behave as a function of Ω M and w, the equation of state parameter, in a flat wCDM model and as a function of Ω M and Ω k , the curvature density parameter, in a non-flat ΛCDM model (see their Figures 3 and 4).Following the prescription of [12], one can fit any cosmological model, also more complex than the standard flat ΛCDM model, by applying the "varying evolution" approach.Thus, it will allow us to avoid the circularity problem.

Cosmological Fit
We have employed the final QSO samples to fit with the Kelly method a flat ΛCDM model, in which we fix H 0 = 70 km s −1 Mpc −1 , and we consider Ω M as a free parameter with a wide uniform prior between 0 and 1 (Figures 11 and 12).Under these assumptions, the formula for the luminosity distance D L reads (in units of Megaparsec) as Recalling Equation ( 1), we can notice that D L can be obtained from the observed quantities F UV and F X as a function of the parameters of the relation.Thus, combining Equations ( 1) and ( 4), we are able to fit the cosmological free parameter Ω M and the free parameter of the RL relation.Hence, we have also left γ ′ , β ′ , and δ ′ free to vary, we have imposed on them the uniform priors 0 < γ ′ ≤ 1, 0 < β ′ < 20, and 0 < δ ′ < 1, and we have applied the best-fit cosmological likelihoods: a Gaussian likelihood for the sample of 1065 sources and a logistic likelihood for the sample of 1132 QSOs.Indeed, Refs.[24,68,79] has proven that, as the Gaussian assumption is not satisfied, the commonly used Gaussian likelihood is not the appropriate likelihood to be applied for cosmological applications of SNe Ia of Pantheon and Pantheon + samples, the whole QSO sample of [40], and BAO, and that adopting the correct likelihood is crucial to reduce the uncertainties on cosmological parameters.Thus, following these works, we have checked the normality assumption for our final QSO samples uncovering that it is verified by the sample of 1065 QSOs, while it is not fulfilled by the sample of 1132 sources.Indeed, for this sample, the best-fit likelihood is a logistic one, as for the initial 2421 QSOs [24].We here point out that the fact that the best-fit likelihood for the sample of 1065 sources is Gaussian, differently from the one of the original sample of 2421 QSOs, does not contradict the fact that the initial sample is the parent population of the selected one, as verified through the Anderson-Darling two-sample test (see Section 3.1).Indeed, the two tests relate to different quantities: the Anderson-Darling two-sample test is applied to the distribution of fluxes, while the Anderson-Darling test for normality investigates the normalized residuals of luminosities.We refer to [24] for a detailed analysis of the non-Gaussianity of the full sample of QSOs.We notice that we always test the best-fit probability density function for the assumed cosmological model (e.g., Ω M = 0.3, H 0 = 70 km s −1 Mpc −1 ) once luminosities are corrected for the evolution.Hence, we have applied the best-fit likelihoods for each sample to fit the flat ΛCDM model.We here note that, since the best-fit distribution for the initial QSO sample is a logistic one, the tails of this distribution cannot be neglected and the standard deviation of this logistic distribution is rather large.In this fitting procedure, we have also accounted for the effects of the evolution in redshift of QSO luminosities.Indeed, to fit a cosmological model, and thus explicitly show the dependence on D L , we need to turn fluxes into luminosities according to L X,UV = 4 π D 2 L F X,UV , where D L is provided by Equation (4).As anticipated in Section 3.2, we here have applied the most general method for correcting for this evolution, the "varying evolution" [12], in which the correction varies as a function of Ω M .with best-fit values with 1σ uncertainties: γ ′ = 0.60 ± 0.01, β ′ = 7.9 ± 0.2, δ ′ = 0.077 ± 0.002, and Ω M = 0.231 ± 0.122.The dark region shows the 68% of probability of the parameters at play, while the lighter blue region the 95%.

The Gold Sample of 1132 QSOs
We here outline the main results obtained for the gold sample of 1132 QSOs, since, compared to the one of 1065 sources, this is the one that gives the best precision in terms of the cosmological results.This sample presents a dispersion in fluxes of δ F = 0.22, which is 24% less than the dispersion of the original sample (δ F = 0.29), and δ ′ = 0.07 vs. δ ′ = 0.09 for the luminosity relation.It still covers the whole redshift range from z = 0.009 to z = 7.54, and it is smaller in size compared to the initial sample (53% of sources are discarded).We here stress that the dispersion δ F cannot be derived from the values reported in Table 1.Indeed, Table 1 provides the best-fit values of the intrinsic dispersion in each redshift bin investigated, while the intrinsic dispersion of the final selected sample must be computed by fitting together the sources in the whole redshift range.More precisely, δ F = 0.22 is obtained by fitting the flux-flux linear relation log 10 F X = a log 10 F UV + b on the whole redshift range covered by the selected sample.It is also visible from Figure 5 that the dispersion is reduced compared to the one of the full sample.We note that our fitted RL parameters are re-estimated after removing the "outlying" observations.Thus, we obtain a new set of residuals.This means that the new set of residuals is not a truncated version of the original residuals.Moreover, our final sample still presents the same features as the parent sample.Indeed, we have applied the Anderson-Darling two-sample test to check that the two distributions in fluxes in each bin are drawn from the same parent population of the initial sample in the same bin.This allows us to statistically affirm compatibility with the null hypothesis that we are not introducing biases or significant changes in the initial QSO sample from a physical point of view.In addition, the reduction in the sample size is not surprisingly small, since the 1048 Pantheon SNe Ia have been slimmed down from an original sample of 3473 events, with a reduction in size of the 70% of the starting data set [84].In this regard, we acknowledge that we start from an already selected QSO sample [40], which has been determined from a much larger sample.Nevertheless, since we perform a cosmological analysis on QSOs, we do not start from all the original observed sources presented in the catalogs, but from the sample described in [40], which contains only the sources that are standardizable cosmological candles and thus can be used in cosmological studies.As a further step, we have also proved that our final sample still follows the RL relation after correcting for selection effects and evolution in redshift of luminosities.Indeed, as shown in Figure 9, we have obtained values of the slope and normalization consistent with those of the corresponding L ′ X − L ′ UV relation for the original sample [21].
Also, we have fitted a flat ΛCDM model, as detailed in Section 3.3, by fixing H 0 and applying the "varying evolution" method with the function k(Ω M ), in which the evolutionary coefficients of the EP method vary along with the free parameter Ω M to avoid any circularity problem.We have also obtained closed contours on Ω M along with a significant reduction in the uncertainty on Ω M compared to the one obtained with the whole QSO sample, which is 0.210 [12].Indeed, with the gold sample of 1132 QSOs, we have obtained Ω M = 0.256 ± 0.089, as shown in Figure 11, with a precision improved of 58% compared to 0.210.This precision is also slightly improved compared to the one reached with SNe Ia by [85], which is 0.10, even though in the case of SNe Ia Ω M is not the only free parameter of the fit.Since we are aware that QSOs at z > 3 (see, e.g., [86] for more details on the sources at z∼3) show a different distribution in the flux-flux plane, compared to the sources at lower redshift (see, e.g., Figure 3 of [36] and the left and right panels of Figure 6), we expect that the use of QSOs in bins at z > 3 could reduce the precision on the fitted Ω M .We refer to [87] for a theoretical analysis of the low-and high-redshift QSOs.In addition, the different cosmological roles of high and low-redshift QSOs can also be read in terms of the effects that should arise from the diversity of QSOs according to the Quasar Main Sequence [88][89][90][91][92]. Indeed, we expect a bias at higher redshift, where only extreme accretors are observed, while it may not be the case for closer QSOs, and that could affect the results.Hence, we have performed our analysis also considering only the bins at z < 3.
With this cut in redshift and including also the untouched sources, our selected sample is reduced from 1132 to 1062 QSOs.This sample constrains Ω M = 0.203 ± 0.073.As expected, the sample cut at z < 3 improves the sensitivity on Ω M by reducing the uncertainty by a factor of 18%.Nevertheless, we here stress that, for our study, it is important to consider QSOs at all redshifts because the restriction to a particular redshift range in the analysis may bias the results, as it will introduce further incompleteness in the sample.

The Comparison of the Two Final QSO Samples
As shown in Figure 4, it is also possible to achieve an increased reduction of the intrinsic dispersion, compared to that of the 1132 QSOs, if we discard the untouched sources (see Section 3.1).Following this approach, we have defined a QSO sample with 1065 sources with reduced redshift coverage between z = 0.2 and z = 3.4 and intrinsic dispersion of the flux relation δ F = 0.18, which is 18% less than the dispersion of the 1132 sources and 38% less than the dispersion of the original sample.Moreover, as for the sample of 1132 QSOs, the 1065 sources still present the same features as the parent sample, as tested through the Anderson-Darling two-sample test, and recover the RL relation, as shown in Figure 8.Compared to the gold sample of 1132 sources described above, which yields Ω M = 0.256 ± 0.089, the sample of 1065 QSOs, results in Ω M = 0.231 ± 0.122, as presented in Figure 12.
A remarkable difference between the gold samples of 1132 and 1065 is that in the smaller sample, the highest redshift is z = 3.435.The larger sample of 1132 contains 45 more QSOs, which are in the range of 3.435 ≤ z ≤ 7.5413 and are absent in the 1065 QSO sample.These additional high-z QSOs are distributed according to the redshifts as follows: 9 QSOs within 3.435 < z ≤ 4, 19 within 4 < z ≤ 5, 6 within 5 < z ≤ 6 and 11 within 6 < z ≤ 7.5.The left and right panels of Figures 6 and 7 show these differences since the color bar axes are different in the two panels: the left panel has a redshift ranging from 0.009 to 7. The right panel has a redshift ranging from 0.1948 to 3.435.This difference in the redshift coverage is also clearly visible looking at Figure 10 that shows the redshift distributions of the two samples and when comparing the top left and right panels of Figure 13.These figures show indeed the Hubble diagrams (i.e., distance modulus vs. redshift) for both our two final samples (yellow points) along with the 1σ uncertainty on the distance modulus.On the bottom panels of Figure 13, we compare if high-z data in our analysis show a significant deviation from the flat ΛCDM model with Ω M = 0.3.This analysis is similar to the one presented by [93], where a 4-σ incompatibility was claimed for a set of sources in the redshift range 3 < z < 3.3.Moreover, Ref. [93] obtained lower luminosity distance at high-z than predicted by flat ΛCDM with Ω M = 0.3.For the purpose of comparison with the aforementioned work, we averaged redshifts and the best-fit luminosity distances for the sources of our gold sample in the same redshift range.We computed the error bars as a simple standard deviation.The results are shown on the bottom-right panel of Figure 13.Noticeably, we obtain compatibility within < 2σ, and our averaged point (red) lays above the ΛCDM with Ω M = 0.3 (green line).The results presented here and in [93] are a reflection of obtained values of Ω M in treatments with and without correction for evolution.The high luminosity distance at high-z in our sample leads to a slightly smaller value of Ω M than 0.3.Similarly, when no correction for selection bias and redshift evolution is applied, one gets Ω M going towards the value of unity since the value of luminosity distance is much smaller than the one predicted by Ω M = 0.3 as obtained by [93].At the bottom left panel of Figure 13, we performed the same analysis as above.The only difference is that we averaged all the sources with z > 3. We deduct the same conclusion, but here, the averaged point is compatible with the ΛCDM (Ω M = 0.3) within < 1σ.
As already explained in Section 4.1, we have also tested to what extent the precision on Ω M is worsened by the inclusion of QSOs at z > 3. To this end, we have cut our sample, retaining only sources at z < 3 and discarding the untouched sources.This procedure trims the sample from 1065 to 1040 sources.This cut sample constraints Ω M = 0.180 ± 0.082, with an uncertainty reduced by a factor 33%, compared to the sample of 1065 QSOs in the whole redshift range.Nonetheless, as stressed above, it is important for our analysis to consider the QSO sample at all redshifts to avoid possible biases in our results.As just stressed, we can notice that the sample composed of 1065 sources presents a smaller intrinsic dispersion compared to the one of 1132 QSOs: δ F = 0.18 vs. δ F = 0.22 for the flux relation, and δ ′ = 0.07 vs. δ ′ = 0.09 for the luminosity relation.This is due to the fact that the smaller sample is obtained only by considering inliers determined through the Huber regression, while the larger sample includes sources at low (z < 0.2) and high (z > 3.4) redshift, for which, due to the insufficient statistics of the corresponding redshift bins, was not possible to perform the fit and thus remove outliers.The addition of these sources, which cannot be identified either as inliers or outliers, increases the dispersion of points, as clearly visible from the comparison of Figures 4 and 5. Despite this reduced intrinsic dispersion, the 1065 QSOs lead to a larger uncertainty on Ω M , 0.12, compared to the uncertainty of 0.09 reached by the sample of 1132 sources.The reason is that the latter sample covers the whole redshift range of the initial QSO sample, filling also the low-redshift interval that allows anchoring the Hubble diagram to the zero-point near z = 0, thus better constraining the value of the matter density today (at z = 0).The difference in redshift coverage of the two QSO samples is also reflected by the different bestfit cosmological likelihoods: the larger sample is better fitted with a logistic likelihood, as the original QSO sample, since it resembles the initial 2421 sources in the redshift coverage, while the smaller sample fulfills the Gaussianity assumption as it ranges over a reduced redshift interval.Based on all these results, we here claim that with more QSOs at low and high redshift (z < 0.2 and z > 3.4), that would allow us to remove outliers in these redshift ranges, we could reach a reduced intrinsic dispersion and smaller uncertainty on Ω M , which marks the relevance of future QSO surveys and observations.In Appendix A, we show how our main cosmological results are completely independent of the binning procedure.

The Advantage of Using the Huber Regression Technique versus the Standard Fitting Methods
As anticipated in Section 3.1, we have also estimated the gain achieved with the use of the Huber regressor by comparing the results obtained with the Huber technique with the ones obtained with a traditional fitting method.To this end, we have applied the sigma-clipping procedure (see [72] for details to select the final QSO samples), and also, we have performed the cosmological analyses of these new samples.More specifically, we have chosen the threshold of the sigma-clipping so that the size of the new final samples (with and without the addition of the untouched sources) is similar to the size of the corresponding samples obtained using the Huber algorithm.This allows us to be consistent with our analysis with the Huber approach.Specifically, by choosing a threshold for the sigma-clipping of 1.5, we have obtained a sample of 1078 QSOs when we do not include the untouched sources (compared to the 1065 of the Huber case) and a sample of 1145 QSOs when we include the untouched sources (compared in the 1132 of the Huber case).We have also fitted the flat ΛCDM model with Ω M free to vary, obtaining Ω M = 0.244 ± 0.132 for the sample of 1078 sources and Ω M = 0.307 ± 0.106 for the larger sample of 1145 QSOs.Hence, we can notice that the values of Ω M obtained with the traditional technique and with the Huber regressor are compatible within 0.07σ (for the sample without additional sources) and within 0.4σ (for the sample with additional sources), while the uncertainties on Ω M are reduced by 8% when we employ the Huber method (for the sample without untouched sources) and by 16% (for the sample with untouched sources).Thus, this comparison shows that the best-fit of Ω M is unchanged, as expected, but the application of the Huber selection technique allows us to improve the precision on Ω M since it is able to better detect outliers of the RL relation and better determine the "true" parameters of the relation itself, thus selecting a QSO sample that is more powerful to infer cosmological parameters.

The Comparison of These New Gold Samples (with and without the Untouched Sources) with the RL Relation for the Total Initial Sample
We here detail the comparison between the correlation in the de-evolved L ′ UV − L ′ X plane with the RL relation valid for the total sample and not corrected for selection biases.
In Figure 11, the best-fit parameters obtained when leaving also Ω M free to vary are γ ′ = 0.61 ± 0.01, β ′ = 7.8 ± 0.2, δ ′ = 0.084 ± 0.003, and Ω M = 0.256 ± 0.089 for the sample of 1132 QSOs.Without applying any correction for redshift evolution and selection biases, Ref. [40] found the following values: γ = 0.586 ± 0.061 and δ = 0.21 ± 0.06 with a simple forward fitting method.Hence, the slope obtained with our computation is compatible within 0.39σ with the one of the empirical RL relation, while the value of the intrinsic scatter is reduced by 60%.Interestingly, the smaller sample of 1065 QSOs yields the following bestfit parameters: γ ′ = 0.60 ± 0.01, β ′ = 7.9 ± 0.2, δ ′ = 0.077 ± 0.002, and Ω M = 0.231 ± 0.122 (see Figure 12).In this case, the slope of the RL relation and that derived with our approach are compatible within 0.23σ and the value of the intrinsic scatter is reduced by 63%.
Remarkably, if we consider a flat ΛCDM model in which both the value of Ω M and H 0 are fixed, the smaller sample composed of 1065 QSOs, shown in Figure 8, has similar values of the best-fit parameters: γ ′ = 0.60 ± 0.01, β ′ = 7.96 ± 0.20, and δ ′ = 0.069 ± 0.002.In this case, the slope of the RL relation and that obtained with our approach are compatible within 0.23σ and its intrinsic scatter is reduced by 67%.Similarly, for the first and larger sample, shown in Figure 9, the resulting best-fit parameters are: γ ′ = 0.61 ± 0.01, β ′ = 7.67 ± 0.21, and δ ′ = 0.088 ± 0.002.This shows that the slope is again within 0.39σ, similarly to the case in which instead Ω M is free to vary.In this case, the intrinsic scatter is reduced by 58%.We here note that the slope is degenerate with the normalization, but we have checked that, when we remove the degeneration by scaling the variables, the results remain compatible.We note that fixing Ω M or leaving it free to vary allows a difference of 2σ in the intrinsic scatter, where the scatter is smaller when all parameters are fixed.We expect indeed smaller values for fewer degrees of freedom.

The Need for This Analysis and the Interpretation of Results from a Physical Point of View
To better understand the origin of a larger intrinsic dispersion in the QSO sample of [40], we have compared the kurtosis of our sample to that in [40].The kurtosis is indeed the fourth standardized moment, which identifies extreme values in the tails of the distribution compared to Gaussian tails.As reported in [24], the full QSO sample shows a kurtosis of ∼0.8.On the other hand, for our final sample, we obtain a kurtosis of ∼−0.2 (with and without correction for evolution).Nevertheless, the larger kurtosis of the sample of [40] is not the only cause of the larger intrinsic scatter of this sample.Indeed, the kurtosis in the whole data set is not very large.The larger intrinsic scatter in the whole sample means that the standard deviation (not the kurtosis) is larger.It is quite straightforward that the standard deviation (intrinsic scatter) will be smaller in our reduced sample.This is the reason why we perform and build the full methodology of all this procedure.
We here would like to stress that, since the physics and the processes that induce the X-UV relation in QSOs are not known, we are not yet able to physically explain all the observational features of the QSO sample used in cosmology.Hence, we cannot provide a physical reason that would force an upper limit of the intrinsic dispersion.Nevertheless, recently, Refs.[57,93] has proven that several factors contribute with different percentages to the overall dispersion of the total sample and that the true intrinsic dispersion is much smaller than the one actually observed.Specifically, these factors are the following: the use of X-ray photometric measurements instead of the spectroscopic ones, the intrinsic variability of the sources, and the inclination of the torus of QSOs.Analyzing the contribution of each of these factors, Ref. [93] has shown that, for a sample of QSOs at z∼3 with high-quality observations, the intrinsic dispersion is only 0.09 dex, which is completely ascribed to the intrinsic variability of QSOs and geometry effects of the sources.The same analysis is performed for the full QSO sample of [40] in [57].In this case, the results are as follows: the intrinsic variability produces 0.08 dex of the intrinsic dispersion, with a larger contribution at low luminosities and a smaller one at high luminosities, while the inclination contributes to 0.08 dex assuming a torus with an opening angle, measured from the disc surface (see Figure 4 of [57] for a graphical representation), of 30°(this value of 0.08 dex is lowered for larger opening angles).The use of photometry instead is negligible.We refer to [57] for a detailed description of these contributions to the dispersion of the RL relation, but we here notice that the contribution of the inclination is computed from mock simulations, also taking into account the limb-darkening effect.Indeed, the authors start with the simple assumption of the absence of an absorbing torus.In this case, the inclination angle θ is randomly extracted from a distribution that is uniform in cos θ.Then, they improve the accuracy of the model by introducing an obscurer.This way, using a mock sample of 100000 QSOs, they derive a contribution of 0.08 dex to the observed dispersion from the inclination.In the end, the intrinsic dispersion proves to be 0.09 dex for the sub-sample of QSOs at z∼3 analyzed in [93] and ∼0.11 for the whole QSO sample studied in [57].The same value of 0.09 dex for the intrinsic dispersion is also obtained in [72], in which QSO sub-samples are selected in redshift bins through the sigma-clipping technique by retaining only the sources that better follow the RL correlation.Moreover, the value of ∼0.09 dex is recovered in [72] independently of the redshift interval investigated.This shows that trimming the sample by selecting only QSOs closer to the ideal RL relation line intrinsically removes low-quality data, which are thus outliers, reducing the dispersion towards its true intrinsic value.The value of 0.11 is much smaller than the actual observed dispersion of the full QSO sample, which is ∼0.2 dex.This means that the current quality of the X-ray and UV flux measurements is not yet sufficient to reveal the true (very small) intrinsic dispersion of the relationship leading to an observed dispersion that is larger than the real one.This proves that selecting the QSO sample to reduce the intrinsic dispersion, as in our procedure, allows us to build a sample much more similar and faithful to the one that properly follows the X-UV relation.
To further comment and interpret our results from the physical point of view of the RL relation, we have also investigated for each binning approach the compatibility between the best-fit values of the slope a in each bin and the mean slope, averaged from the obtained slopes in all bins.To this end, we have computed the corresponding "z-score" parameter defined as (a i − < a >)/ ∆ 2 a i + ∆ 2 < a >, where a i is the slope of the i-th bin and < a > is the mean value of the slope calculated over all bins.This analysis shows that in all the binning cases the slopes in each bin are consistent with the mean value of the slope in < 3σ.Indeed, the worst z-score is the one of the case with bins in log 10 (1/(1 + ∆z)), which has a minimum of −2.4 and a maximum of 2.6, while the best case is the one with the method of optimization of bins starting from the lowest redshift, for which the z-score is between −1.0 and 1.1.
In addition, possible improvements on the sample could depend on the diversity of QSOs according to the Quasar Main Sequence [88][89][90][91].One could expect that at higher redshift we have a bias seeing only extreme accretors, while it may not be the case for closer QSOs, and that could affect the results.The selection of a sub-sample that includes only high accretors, may improve this bias [90].This means that it would be interesting to perform an analysis on only extreme accretors, thus selecting only high accretors also at low z.Of course, one must be careful not to include the false candidates (see, for example, [92]) Although this is a very interesting topic, this analysis goes beyond the scope of this work.
In summary, our aim is to show that, if a proper sample with a reduced intrinsic dispersion is defined, then QSOs can be used as standalone probes with the precision of SNe Ia.Our gold QSO samples can help to reveal physical properties common to these sources in order to identify a QSO sample driven by fundamental physics.However, the investigation of these physical properties goes beyond the scope of the current paper.Overall, we are discarding outliers, identified this way by the robust Huber regressor, to define a QSO sample that is able to constrain Ω M with unprecedentedly high precision.The Huber algorithm is based on the identification of the sources that follow the true slope of the relation (inliers) and the ones that do not (outliers).In the end, the QSOs that are detected as outliers, even with measurement errors smaller than the intrinsic dispersion observed, could be affected by some observational problems as they do not follow the RL relation.An example of such an issue could be the host extinction for photometric measurements, where the lack of knowledge of the spectrum does not allow for a precise determination of the QSO's unabsorbed luminosity.

Summary, Discussion, and Conclusions
To uncover the ultimate QSO sample that can be used as a powerful cosmological tool, we started from the most comprehensive and up-to-date data set for cosmological studies [40], comprising 2421 sources spanning from z = 0.009 to z = 7.54, and we ushered in an original approach.We stress that, differently from other works [19,40], we use the full sample at all redshifts and we also correct for redshift evolution to ensure the utmost accuracy and reliability for cosmological applications.Our procedure to determine such a gold sample is general and versatile, and hence can also be used for other probes (e.g., GRBs) and larger samples.Additionally, this method is completely model-independent and thus, it avoids the circularity problem.Indeed, we apply our technique to the linear relation between the logarithms of QSO fluxes in UV and X-ray, that do not depend on the choice of a specific cosmological model, contrary to the case of luminosities.We employ the robust Huber regression which allows us to reduce the intrinsic scatter of the F X − F UV relation by removing the sources that are identified as "outliers" [62], a notation that we use not in its strict statistical sense but rather to refer to the QSOs that show more discrepancy from the best-fit relation line.The novelty in our approach lies in harnessing the remarkable power of the Huber fitting method to unearth the optimal QSO sample.Indeed, as detailed in Section 4.3, the employment of the Huber regressor in place of a standard fitting technique allows us to determine Ω M in a flat ΛCDM model with increased precision.Indeed, the Huber procedure better distinguishes the outliers of the RL relation, leading to a QSO final sample that better follows this relation and has more constraining power on the cosmological parameters.We here outline our main results and draw conclusions.

1.
The strategy in the selection of the QSO gold sample.Since our main challenge is to constrain cosmological parameters, such as Ω M , we strive not only to reach the smallest dispersion of the F X − F UV relation but also to keep a statistically sufficient number of sources in each redshift bin.This guarantees that the fitting is still possible from a statistical point of view.Hence, we have found a compromise between these two factors which are antagonistic.The optimal number of sources found is 1132 since these sources fulfill all the following required criteria: the minimum number of sources (i.e., 10 in our case), the Anderson-Darling two-sample tests in each bin, and the requirement on the distance luminosity that should be negligible compared to the dispersion of the relation in flux in each bin.This sample of 1132 QSOs presents δ F = 0.22 and still covers the whole redshift range of the original sample, from z = 0.009 to z = 7.54.The intrinsic dispersion of the flux relation can be even reduced if we discard the sources belonging to redshift bins with not enough statistics (i.e., <10), the "untouched" sources, to perform the Huber regression.By applying this choice, we have defined a QSO sample with 1065 sources with reduced redshift coverage between z = 0.2 and z = 3.4 and intrinsic dispersion of the flux relation δ F = 0.18.

2.
Comparison with the original RL relation and cosmological results.We have proven that the RL relation is still verified by our two final samples, once accounted for the redshift evolution.We have also used the obtained QSO gold samples for cosmological application to derive Ω M by fitting a flat ΛCDM model leaving contemporaneously free both Ω M and the RL relation parameters, γ ′ , β ′ , and δ ′ , while fixing H 0 = 70 km s −1 Mpc −1 .We have performed this fit by using the best-fit proper likelihoods for the samples, which proved to be a logistic one for the sub-sample of 1132 QSOs and a Gaussian one for the 1065 sources.We have also performed the fit by applying the correction for the redshift evolution as a function of Ω M .The gold sample of 1132 QSOs has provided Ω M = 0.256 ± 0.089, whereas the sample with 1065 sources has led to Ω M = 0.231 ± 0.122.Hence, we have reached a precision improvement of 58% compared to the one obtained with the whole QSO sample (i.e., 0.210, see [12]).
Moreover, these values are compatible with the current value of Ω M = 0.338 ± 0.018 [71], and hence in agreement with the expected value of the current matter density.Additionally, the obtained values of Ω M are compatible in ∼1σ with the one reported in [59], in which a non-binned analysis, independent from the one here presented, is performed to select the QSO sample and constrain Ω M .We here point out that our analysis is not biased or induced by any circular reasoning.The sample is trimmed by reducing the uncertainties in the flux-flux relationship, which does not depend on the cosmological parameters once we bin the data (see Section 3.1), but rather on the intrinsic scatter of this relation.Our results show that, after restricting our attention to the probes for which the intrinsic scatter is small, we obtain a substantially improved precision of the estimated cosmological parameter.

3.
The impact of bin division on cosmology.The above-detailed results have been obtained by trimming the initial QSO sample in bins of log 10 (1/(1 + ∆z)), as described in Section 3.1, since this is the approach that leads to the best cosmological results.Nevertheless, to further investigate the impact of this choice for the division into bins and to free our analysis from the arbitrariness and possible issues of the binning procedure, we have also performed our study by applying the three different selection methodologies outlined in Sections 3.1.2-3.1.4.We have thus proved that our results do not depend on the specific methodology employed to select the QSO sample.Indeed, as detailed in Appendix A and Table 2, the values of Ω M obtained in all the cases investigated are compatible within 1.2 σ and they are also consistent with Ω M = 0.338 ± 0.018 reported in [71].This compatibility with the most recent value of Ω M measured from SNe Ia ensures that our analyses recover the expected cosmology independently of the binning approach considered.Moreover, the comparison among the results obtained from different procedures has also shown that larger sample sizes shift Ω M towards values of ∼0.5, as expected from [12].Additionally, this extended analysis has also suggested that the low-z QSOs tend to lower the value of Ω M .Finally, the employment of these several independent approaches to select the final QSO samples has led to an extensive knowledge of the QSO selection and trimming procedure.Furthermore, leveraging the advantages of the methods investigated, we have avoided our analyses and results being biased or distorted by the arbitrary choice of a fixed binning.The outcomes of all these approaches have established the validity of our results and the robustness of our analysis.4.
The need for a larger sample and a physical interpretation.Based on all these considerations, we here point out that we would need a much larger sample with these properties to reach the current value of precision obtained with SNe Ia, which is δ Ω M = 0.018.This situation is similar to the one occurring for GRBs [45] where the precision of SNe Ia in [85] is reachable now, while for an improved precision as the one the SNe Ia have today, we would need to wait for two other decades.Nevertheless, if a common emission mechanism or properties in this sample will be driven by fundamental physics the waiting time for reaching this precision would be considerably improved.We here point out a parallel case for GRBs: if the plateau emission in the Platinum sample is driven by the magnetar, a given sample of GRBs with peculiar magnetic fields and spin period can drive the standard set.Here, the gold QSOs can help to reveal the physical meaning of these properties.Indeed, the current paper can allow us to identify a QSO sample, highlighted by our statistical procedure, that is driven by fundamental physics.However, this investigation goes beyond the scope of the current paper and will be analyed in a forthcoming paper.
In conclusion, we have shown that QSOs alone with the RL relation can now be promoted to reliable standard candles to measure cosmological parameters, such as Ω M , when a gold sample is defined.In this framework, we are able to constrain the cosmological parameter Ω M with significant precision at high redshift, up to 7.5.This ushers a new era of QSOs as effective standard candles, in which the efforts of the QSO community can be driven to understand and delve into the differences between this gold sample and the total one.
A scientific revolution is underway.QSOs, once enigmatic cosmic entities, now stand tall as standard candles illuminating the darkest corners of our Universe.The gold sample is our key to unlocking their true potential, propelling us towards an improved precision in cosmological studies.The journey has just begun, and the tantalizing mysteries that lie ahead beckon us to delve deeper.on the statistics and Takashi Hamana, Guido Risaliti, and Elisabeta Lusso for the discussion on the selection methodology.

Conflicts of Interest:
The authors declare no conflicts of interest.

Appendix A. The Impact of Different Binning Approaches on Cosmology
We here detail and discuss how the different binning approaches do not influence the cosmological results.In particular, we compare the results outlined above and derived from the binning in log 10 (1/(1 + ∆z)) with the ones obtained with the three different binning approaches detailed in Sections 3.1.2-3.1.4.By examining Table 2, we can compare results from corresponding cases: the sample in which the untouched sources are not included and the sample in which the untouched sources are added.Based on these samples, we can draw the following conclusions: the values of Ω M are always compatible within 0.8σ and its uncertainty is between 0.089 and 0.232 in all cases.We can also notice a specific trend: for the samples that include the untouched QSOs, the value of Ω M and its uncertainty are smaller compared to the samples in which these sources are discarded.This trend observed for the uncertainty on Ω M is valid for all the different selection approaches studied.We note a decrease of the central value of Ω M , from the sample without untouched sources to the one with untouched sources, in all methodologies except for the case with the binning in log 10 (1/(1 + ∆z)) (see Table 2).The reduction in the uncertainty on Ω M is due to the different redshift coverage of the two samples (with and without the untouched sources), as already detailed in Section 4.2.Considering instead the reduced value of Ω M , the trend suggests that the untouched sources, mainly located at very low z, prefer smaller values of Ω M .In this regard, we can also comment that the different trend of the case of the binning in log 10 (1/(1 + ∆z)) could be explained by considering that, compared to the other cases, the number of untouched sources is larger, and they are mainly located not at low redshifts but at z > 3.5.Hence, the contribution of these sources is expected to be different from the ones at low z, thus not leading to a reduction of Ω M .Nevertheless, this topic deserves to be further investigated in future analyses, e.g., [72].To verify the dependence of high-z QSOs on the precision on Ω M , we computed Ω M with correction for evolution for a whole sample with sources discarded at z > 3.435 (the maximum z of the 1065 sources sample).The posterior distribution is centered at Ω M ≈ 0.39 with a standard deviation of 0.21.With only 45 sources discarded from the whole sample, we obtain significantly less precise results for the determination of Ω M .
In addition, we can also notice from Table 2 that Ω M is not constrained in the case of bins centered on each QSO without untouched sources.This result clearly shows the crucial role of the compromise between a small intrinsic scatter of the flux-flux relation and the number of sources in view of using QSOs for precision cosmology.Indeed, on the one hand, large samples cannot constrain well Ω M due to their large intrinsic scatter.On the other hand, reducing the sample size leads to a reduced intrinsic dispersion and thus to better estimate Ω M but only until the number of sources does not become too small to determine Ω M with closed contours, as in the case of the 811 QSOs.Moreover, we observe that, when the sample size significantly increases, approaching the initial size of 2421 QSOs, the value of Ω M increases toward Ω M ∼0.5.This completely agrees with the value Ω M = 0.500 ± 0.210 reported by [12] and obtained with the full QSO sample.We here also notice that all the obtained values of Ω M are compatible with each other in < 1.2 σ.We show this consistency by reporting in Table 2 the "z-score gold " parameter, which computes the compatibility of each result with the result that presents the smallest uncertainty on Ω M .The z-score gold is computed as , where the index i indicates a given value from the table and the index "gold" refers to the measurement with the smallest uncertainty: Ω M, gold = 0.229 ± 0.061.Remarkably, the Ω M values obtained in this work are also compatible within 1.7σ with the most recent measurement of Ω M from Pantheon + SNe Ia [71], which is Ω M = 0.338 ± 0.018.More specifically, our values are compatible with the one reported in [71] within 0.9σ, with the exception of our case of 1125 selected QSOs, which manifests a discrepancy of 1.7σ.This is shown in the last column of Table 2, which reports the "z-score SNe " parameter.This parameter is calculated with the same formula of z-score gold detailed above, but replacing Ω M, gold with Ω M, SNe = 0.338 ± 0.018.The overall compatibility between the values of Ω M obtained with our selected QSO samples and the value of Ω M reported in [71] significantly proves that the cosmology is reliably recovered independently on the binning approach, and thus on the number of selected sources.

Appendix B. The Varying Evolution Method
Referring to Section 3.2, we here describe additional details of the EP statistical method applied to our selected samples of 1132 and 1065 QSOs considering the simple evolutionary form with g(z) = 1/(1 + z) k , where k is the slope of the power-law (see also Section 3.2).This treatment is similar to the application of this correction to the original sample of 2421 QSOs, detailed in [21].We here stress that we apply the EP method after, and not before, the selection of the QSO sample since we need to apply this correction only once, and we have chosen to apply to the luminosities since anyway the luminosities by definition carry redshift evolution, being dependent on the redshift.
As anticipated, from the measured flux we compute the luminosity for each QSO assuming a flat ΛCDM model with Ω M = 0.3 at the current time and H 0 = 70 km s −1 Mpc −1 .We also compute the flux limit F lim and the corresponding luminosity L min (z i ).Specifically, we have chosen F lim = 6 × 10 −29 erg s −1 cm −2 Hz −1 for the UV and F lim = 8 × 10 −33 erg s −1 cm −2 Hz −1 for the X-rays.We have also verified through the means of the Kolmogorov-Smirnov (KS) test that, for the full and the cut samples in both X-rays and UV, the probability of the null hypothesis that the two samples are drawn by the same distribution cannot be rejected at the p-value of p = 0.79 (for the sample of 1132 QSOs) and p = 0.40 (for the sample of 1065 QSOs) for the UV and p = 0.50 (for the sample of 1132 QSOs) and p = 0.38 (for the sample of 1065 QSOs) for X-rays.The limiting values for L UV and L X corresponding to the above-mentioned values of F lim are shown with a black continuous line in the left and right panels of Figures A1 and A2, respectively, over the whole set of data points represented by blue filled circles.We then apply the τ test to the data sets trimmed with values of the fluxes mentioned above and obtain the trend for τ(k) shown in the left and right panels of Figures A3 and A4 for the UV and X-rays, respectively.For the UV and X-rays, respectively, we obtain k = 4.39 ± 0.12 and k = 3.31 ± 0.08 for the sample of 1132 QSOs and k = 4.44 ± 0.13 and k = 3.39 ± 0.08 for the sample of 1065 QSOs.The fact that the values of the evolutionary parameter k are significantly different from 0 for both samples and both in UV and X-ray strongly proves that the redshift evolution plagues our final samples, as well as affects the original sample, as detailed in [21].
It is remarkable that the evolutionary function of the UV in our samples is compatible within 2.8σ with the optical evolutionary coefficient k opt obtained in [83], where the same form of g(z) is used.In their paper, they found k opt = 3.0 ± 0.5 and corrected the luminosity function.After we insert our values of k in g(z), we then compute the new de-evolved luminosities and the associated uncertainties for the whole QSO samples.The associated uncertainties are computed by propagating the uncertainties of the obtained k values, together with the uncertainties of the initial non-corrected luminosities, in the formula of L ′ (i.e., L ′ = L/(1 + z) k ).The comparison between these quantities and the initial ones is shown in Figure A5 in the (logL UV , logL X ) plane for both the sample of 1132 QSOs (left panel) and the sample of 1065 QSOs (right panel).Compared to the initial ones, the computed luminosities span a smaller region of the (logL UV , logL X ) plane and show a slightly greater dispersion.This fact is expected because the g(z) function, once the best-fit values for k are used, yields a greater correction (i.e., lower de-evolved values) for higher luminosities.In addition, we have accounted for the error on the determination of k by propagating the errors on the g(z) function.This naturally increases the associated uncertainties on the luminosities.The correction for g(z) affects the spread of the luminosities, hence the dispersion of the correlation, which is consequently larger.To summarize, the dispersion increases due to the larger spread of the luminosities, and it is minimally affected by the error propagation due to g(z).In other words, the dispersion yielded by the function g(z) is larger than the contribution given by the additional errors due to g(z).Larger errors on the variables may reduce the dispersion, but in this case not sufficiently enough to balance the increase of the dispersion due to the function g(z).Figure A5 clearly shows the effects of the application of the EP method on our data.Figure A6 shows the values of k obtained for a grid of values of Ω M for our gold sample of 1132 QSOs.We here stress that this figure is different from Figure 4 in [21] in which the full sample of 2421 QSOs is considered.

Figure 1 .
Figure 1.The flowchart of our selection procedure.The two violet boxes show the start and the end of the scheme, while purple, yellow, blue, and orange boxes identify different binning methods.

Figure 2 .
Figure 2. The slope, a, as a function of the average redshift for all bins.

Figure 3 .
Figure 3. Results of the Huber regression in redshift bins.Left panel: Outliers (in red) and inliers (in green) sources for the second most populated redshift bin along with the Huber best-fit shown with the black line.Right panel: Same as the left panel showing the second least populated redshift bin.

) 29 . 1 )Figure 4 . 1 )Figure 5 .
Figure 4.The left panel: The QSO sample of 1065 sources (in bright yellow) generated from the F X − F UV relation without adding the sources in redshift bins with insufficient statistics with the best-fit parameters being a = 0.704 ± 0.013, b = −12.01± 0.36, and δ F = 0.183 ± 0.004.The original parent sample is superimposed with gray points.Right panel: The same as the left panel but in the luminosity-luminosity space, once the correction for redshift evolution and selection biases has been applied.In this case, we assume a flat ΛCDM model with Ω M = 0.3 and H 0 = 70 km s −1 Mpc −1 .

28 FFigure 6 .
Figure 6.Left panel:The logarithmic flux-flux space with the gold sample of 1132 QSOs color coded according to the redshift.Right panel: Same as the left panel but for the smaller sample of 1065 QSOs.

Figure 7 .
Figure 7. Left panel: The logarithmic luminosity-luminosity space corrected for the evolution on both luminosities with the gold sample of 1132 QSOs color coded according to the redshift.In this case, we assume a flat ΛCDM model with Ω M = 0.3 and H 0 = 70 km s −1 Mpc −1 .Right panel: Same as the left panel but for the smaller sample of 1065 QSOs.

Figure 8 .
Figure 8. Corner plot of the L ′X − L ′ UV relation, corresponding to the F X − F UV relation, once corrected for the effects of selection and the evolution in redshift for the sample of 1065 QSOs.The resulting best-fit parameters are: γ ′ = 0.60 ± 0.01, β ′ = 7.96 ± 0.20, and δ ′ = 0.069 ± 0.002.Here, the fiducial cosmology assumed is a flat ΛCDM model with Ω M = 0.3 and H 0 = 70 km s −1 Mpc −1 .The brighter color indicates the probability of the occurrence of the parameters at 95% and the darker color at 68%.

Figure 9 .
Figure 9. Corner plot of the L ′X − L ′ UV relation, corresponding to the F X − F UV relation, once corrected for the effects of selection and the evolution in redshift for the sample of 1132 QSOs.The resulting best-fit parameters are: γ ′ = 0.61 ± 0.01, β ′ = 7.67 ± 0.21, and δ ′ = 0.088 ± 0.002.Here, the fiducial cosmology assumed is a flat ΛCDM model with Ω M = 0.3 and H 0 = 70 km s −1 Mpc −1 .The dark region shows the 68% of probability of the parameters at play, while the lighter blue region the 95%.

Figure 10 .
Figure 10.Left panel: The redshift distribution of the sample with 1132 sources.Right panel: Same as the left panel but for the smaller sample of 1065 QSOs.3.1.3.The Impact of the Binning on the Data Analysis: Optimization of the Width of Bins

Figure 11 .Figure 12 .
Figure 11.Results obtained from the gold sample of 1132 QSOs from the cosmological fit of the flat ΛCDM model with γ ′ , β ′ , and δ ′ of the RL relation, corrected for redshift evolution in the luminosities and Ω M left as free parameter together with the ones of the relation.H 0 is fixed to 70 km s −1 Mpc −1with best-fit values with 1σ uncertainties: γ ′ = 0.61 ± 0.01, β ′ = 7.8 ± 0.2, δ ′ = 0.084 ± 0.003, and Ω M = 0.256 ± 0.089.The dark region shows the 68% of probability of the parameters at play, while the lighter blue region the 95%.

Figure 13 .
Figure 13.Upper left panel: The Hubble diagram of the gold sample of 1132 QSOs (in bright yellow) derived by assuming a flat ΛCDM model with H 0 fixed to 70 km s −1 Mpc −1 and Ω M = 0.256 ± 0.089, as obtained by fitting this sample.The error bars represent the statistical 1σ uncertainties.Upper right panel: Same as the left panel but for the sample of 1065 QSOs, for which the best-fit value of Ω M is 0.231 ± 0.122.Bottom left panel: Same as above, but the QSOs at z > 3 are averaged into one data point (shown in red).Bottom right panel: Same as above, but the QSOs at 3 < z < 3.3 are averaged into one data point (shown in red).Plots in the bottom part mimics the analysis performed by[93].The error bars on the averaged data points are computed as simple standard deviation.The green line is a plotted theoretical luminosity distance for a flat ΛCDM model with Ω M = 0.3.

- 1 Hz - 1 )Figure A1 .- 1 Hz - 1 )Figure A2 .
Figure A1.Redshift evolution of L UV (left panel) and L X (right panel) in units of erg s −1 Hz −1 for the QSO sample of 1132 sources.The black line in both panels shows the limiting luminosity chosen according to the prescription here described.

Figure
Figure A3.τ(k) function (dashed red line) for both the UV (left panel) and X-ray (right panel) analyses for the 1132 QSOs.The point τ = 0 gives us the k parameter for the redshift evolution of L UV and L X , while |τ| ≤ 1 (gray lines) the 1σ uncertainty on it (dashed purple lines).
Figure A4.τ(k) function (dashed red line) for both the UV (left panel) and X-ray (right panel) analyses for the 1065 QSOs.The point τ = 0 gives us the k parameter for the redshift evolution of L UV and L X , while |τ| ≤ 1 (gray lines) the 1σ uncertainty on it (dashed purple lines).

Figure A6 .
Figure A6.The results of the EP method obtained for a grid of values of Ω M for the sample of 1132 sources.k L UV (Ω M ) is shown on the left, while k L X (Ω M ) on the right.The error bars correspond to the 1σ confidence intervals.

Table 1 .
Results of the selection procedure in redshift bins with our fiducial binning in log 10 (1/(1 + ∆z)) (see yellow box in Figure