© The Authors 2025

1. Introduction

The cosmological concordance model, often called Lambda-CDM, or ΛCDM, where CDM stands for ‘cold dark matter’ and Λ represents the cosmological constant, is the best description of the current observations. The accelerated expansion of the Universe (see also Riess et al. 1998; Perlmutter et al. 1999; Weinberg et al. 2013 for a review), the cosmic microwave background (CMB; e.g. Bennett et al. 1996), and the measurement of baryon acoustic oscillations (BAOs; see Eisenstein 2005 for a review) are key predictions of ΛCDM, which led to its general consensus amongst the scientific community. Nonetheless, discrepancies with ΛCDM have also been observed, the prevailing ones being the ‘σ₈ tension’ and the ‘Hubble tension’ (see Abdalla et al. 2022 for a review of cosmological tensions and anomalies). Most recently, cosmological results from the first year of galaxy, quasar, and Lyman-α forest BAO measurements from the Dark Energy Spectroscopic Instrument (DESI Data Release 1) have shown a preference for w₀ > −1 and w_a < 0 when a time-varying energy density equation of state, w(z) = w₀ + w_a z/(1 + z), is considered, which departs from ΛCDM at about the 3σ statistical level, depending on the specific data set used (Adame et al. 2025). The same overabundance of very high-redshift galaxies that are being discovered by the James Webb Space Telescope (e.g. Finkelstein et al. 2023; Harikane et al. 2023) poses several challenges to the timescale of structure formation in a ΛCDM Universe (e.g. Boylan-Kolchin 2023; Melia 2023). In this landscape, our analysis of the distance-redshift relation (also referred to as the Hubble-Lemaître diagram) of luminous, unobscured active galactic nuclei (AGNs, where quasars represent the most powerful members of this population) has also revealed a deviation from the predictions of the flat ΛCDM model, which emerges at high redshift with a statistical significance of 3–4σ (Risaliti & Lusso 2019; Lusso et al. 2019, 2020, see also Bargiacchi et al. 2022). While complementary to the traditional resort to supernovae to estimate the cosmological parameters, quasars extend the distance-redshift diagram to a redshift range that is inaccessible to supernovae (z = 2–7.5). All the above-mentioned findings, if confirmed, would be important indications of the need of new physics beyond ΛCDM.

The determination of the distance versus redshift relation for quasars, and thus of the cosmological model parameters, is obtained from the empirical correlation between the X-ray and UV emission observed in AGNs (Tananbaum et al. 1979; Zamorani et al. 1981; Avni & Tananbaum 1982, 1986). Such a relation is non-linear, and has been observed over almost five orders of magnitude in UV luminosity with a dispersion of just ∼0.2 dex, implying that there must be a ubiquitous physical mechanism that regulates the energy transfer from the accretion disc to the X-ray emitting corona (e.g. Lusso & Risaliti 2017, and references therein). Several models have been proposed to explain the physics behind the X-ray to UV relation, involving, for example, reprocessing of radiation from a non-thermal electron-positron pair cascade (Svensson 1982, 1984; Guilbert et al. 1983; Zdziarski et al. 1990), a two-phase accretion disc model in which the entire gravitational power is dissipated via buoyancy and reconnection of magnetic fields in a uniform hot plasma in close vicinity of the cold opaque disc (Haardt & Maraschi 1991, 1993; Svensson & Zdziarski 1994; Di Matteo 1998), or a viscosity-heated corona (Meyer et al. 2000; Liu et al. 2002a). Magnetic field turbulence has also been recognised not only as a supplementary heating process in the formation of the corona itself (e.g. Galeev et al. 1979; Merloni & Fabian 2001; Liu et al. 2002b), but also as an efficient means for the transport of the disc angular momentum (e.g. Balbus 2003).

Although a comprehensive model is still missing, and the dependence of the coronal properties on black hole mass (M_BH) and accretion rate (M_⊙) remains an open question, there is little doubt that the X-ray to UV relation is inherent to the accretion process itself, and that a tight physical connection exists between these two bands in the spectral energy distribution (SED) of an AGN. We note that relations of some kind are known between various pairs of wavebands, virtually including the AGN emission in any part of the SED from radio to gamma-rays, yet these relations involve a combination of different processes and physical scales and do not automatically hold any cosmological value. Indeed, the use for cosmological purposes of any relation between the emission in two given bands must obey several conditions, besides the non-linearity of its slope. First and foremost, the dispersion should be small, implying a nearly consequential connection (albeit not necessarily known) between the two bands. Moreover, while the correlation has to be observed across a wide redshift range and for numerous sources, its key parameters should not show any evolution with redshift. Here, we aim to demonstrate that the correlation between X-ray and UV emission in AGN is rooted in black-hole accretion physics and boasts all the aforementioned features, and is therefore a valid cosmological tool irrespective of its largely empirical nature.

Given the apparent tension with ΛCDM of the distance-redshift relation of high-redshift quasars, it is mandatory to investigate whether hidden systematics, uncertainties in the data, or possible biases may account for such a discrepancy. Criticism was in fact raised by some authors regarding the reliability of this result. Khadka & Ratra (2020, 2021, 2022) suggested that the disagreement of the quasar distance-redshift relation with respect to the prediction of the flat ΛCDM model at high redshift could be due to the heterogeneity of the quasar sample considered in the analysis of Lusso et al. (2020), and specifically to some issues with the high-redshift subsamples. Their analysis was performed assuming different cosmological models, indicating that the values of the parameters involved in the determination of quasar distances (γ and β, see Section 3) are dependent upon the cosmological model assumed. This, together with the possible evolution with redshift of the same two parameters, would prevent the standardisation of quasars, and hence the applicability of the distances derived in this fashion to cosmological studies (see also Singal et al. 2022; Petrosian et al. 2022). In this paper, we address all the above issues through a quantitative comparison with supernova-derived distances in the common redshift range, complemented by simulations that prove that all the claimed shortcomings are not intrinsic to the data, but naturally arise when the cosmological model assumed in the analysis is incorrect. In other words, the purported non-standardisability of high-redshift quasars is a direct consequence of our ignorance of the true cosmological model.

The manuscript is structured as follows. In Section 2, we briefly describe our latest quasar sample used for cosmological analyses, whilst our method is summarised in Section 3. Section 4 presents a quantitative comparison of the quasar distances obtained through our technique with the ones of type Ia supernovae in the common redshift range, and in Section 5 we investigate any possible redshift evolution of the relation used to standardise quasars. Section 6 includes a set of simulations to demonstrate the intrinsic degeneracy between the determination of the cosmological parameters assuming different input values for Ω_M and Ω_Λ and a redshift evolution of the X-ray to UV relation. In Section 7, we discuss the issue of luminosity evolution. In Section 8, we provide general guidelines that should be followed to use quasars as standardisable candles for cosmology. Our conclusions are drawn in Section 9.

2. The quasar sample

The quasar sample published by Lusso et al. (2020, hereafter L20) contains 2421 sources over the redshift range 0.009 ≤ z ≤ 7.541. The interested reader should refer to Section 7 in L20 for a detailed discussion of the selection criteria. The sample was intended to be highly uniform, as radio-bright, broad-absorption-line, and significantly absorbed quasars at UV and/or X-ray energies were neglected. Nonetheless, some residual contamination from the host galaxy emission in the optical can still be present, especially for the low-redshift (z ≲ 0.5) and low-luminosity (≲10⁴⁵ erg s⁻¹) AGN, for which the contrast between the nuclear and the host-galaxy emission is limited. This additional contribution biases the determination of the continuum emission from the disc, which, for these sources, has to be extrapolated to the rest-frame 2500 Å. Therefore, unless ultraviolet data are available for the low-redshift AGN (i.e. covering the rest-frame 2500 Å, see Section 2.7 in L20), AGN at z < 0.7 should be conservatively excluded from the sample that is used for the calibration of quasars with supernovae, and also from the cosmological analysis. In this paper, we thus use 2036 quasar measurements; that is, the 2023 quasars at z > 0.7 plus the 13 AGN at very low redshift (0.009 < z < 0.087) with UV data from the International Ultraviolet Explorer (IUE) in the Mikulski Archive for Space Telescopes (MAST), for which L_UV could be determined directly from the spectra (see Section 2.7 in L20).

3. The method: Computation of luminosity distances

Whilst the methodology we adopt to derive quasar distances has been presented in depth in a series of works (e.g. Risaliti & Lusso 2015, 2019; Lusso et al. 2020; Moresco et al. 2022), we briefly summarise the main points below. The method is based on the hypothesis of a universal, redshift-independent, non-linear relation between the X-ray and UV luminosities of quasars:

$$\begin{array}{c}\hfill log(L_{\mathrm{X}})=\gamma log(L_{\mathrm{UV}})+\beta ,\end{array}$$(1)

with γ ≃ 0.6 and β ≃ 8.0. We note that the exact values, at this point, are assumed to be an inherent attribute of the accretion process, and are therefore cosmology-independent. We justify this assumption at length in the following. The X-ray to UV relation has been observed since the 1980s (Tananbaum et al. 1979; Avni & Tananbaum 1982, 1986) and widely studied ever since by a number of independent groups and by using different sample selections. In this paper, we use the data in L20, in which L_X and L_UV are the monochromatic luminosities at the rest-frame 2 keV and 2500 Å, respectively, as derived from X-ray and UV photometric measurements. We extensively discuss the choice of the optimal indicators of the X-ray (coronal) and UV (accretion disc) emission in Signorini et al. (2023, see also Jin et al. 2024). While we refer to that work for more details, we have shown that the physical quantities that are more tightly linked to one another are the soft X-ray flux at ∼1 keV and the ionising UV flux blueward of the Lyman limit, and yet the usual monochromatic fluxes at 2 keV and 2500 Å estimated from photometric data provide an almost as tight X-ray to UV relation, which can thus be used to derive reliable quasar distances (Signorini et al. 2023).

If the relation in Equation (1) holds at all redshifts (we discuss this key point further in Section 5), then it is straightforward to derive the relation between X-ray and UV fluxes:

$$\begin{array}{c}\hfill log(F_{\mathrm{X}})=\gamma log(F_{\mathrm{UV}})+\widehat{\beta },\end{array}$$(2)

where

$$\begin{array}{c}\hfill \widehat{\beta }(z)=\beta -(1-\gamma )log(4\pi )-(2-2\gamma )log{D}_L(z),\end{array}$$(3)

and D_L is the luminosity distance¹. It is therefore possible to fit the data (i.e. F_X, F_UV, z) for each quasar, leaving the cosmological parameters of the model describing the D_L(z) expression and the variables γ and β of the L_X − L_UV relation free to vary.

We reiterate that for quasars to be standardisable candles, and thus for this technique to work, the correlation parameters γ and β must be intrinsically cosmology-independent. To validate (or to confute) the reliability of this method, the key question is whether this hypothesis holds. The rejection of this premise immediately leads to the circularity argument raised by Petrosian et al. (2022). In fact, by reading Equation (1) at face value, both γ and β are dependent upon the choice of the cosmological model. Different assumptions clearly result in slightly different values of the observed correlation parameters when analysing any real data set. This is not in contradiction with the underlying hypothesis, as the following analysis, entirely based on Equations (2) and (3), is devised to prove.

It is also worth noting that for any prospective standard candle, a direct observational validation is impossible in the redshift range in which no other established standard candle is present. For instance, one can observationally approve the distance measurements from type Ia supernovae only if these are found in galaxies close enough to have Cepheid-based distances available (up to 40 Mpc; Riess et al. 2016). For supernovae further away, the reliability of the distance measurements is based on two conditions: a deep physical understanding of the process leading to the distance estimates (i.e. the physical grounds of the Phillips relation; Phillips 1993), and the absence of observational and selection biases in the sample. In the case of type Ia supernovae, there is now a general consensus on both the above points, so their use as distance indicators is well established and the Hubble diagram of type Ia supernovae represents one of the pillars of observational cosmology. The situation is certainly different for quasars. In particular, there is not a satisfactory physical explanation of the X-ray to UV relation yet, which should necessarily involve the process of energy transfer from the accretion disc to the X-ray emitting corona to link the values of γ and β to the workings of accretion physics. To assess the reliability of the distance estimates based on the X-ray to UV relation, we need to: (1) rule out any possible redshift-dependent selection effect in the sample, including any evolution in the spectroscopic properties of the selected sources; (2) check for any possible systematic trend in the common redshift range with supernovae; and, based on the previous points, and (3) define a method to use quasars as standardisable candles in combination with other cosmological probes.

The sample selection (first point above) has already been discussed in great detail in our previous works (Lusso & Risaliti 2016; Risaliti & Lusso 2019; Lusso et al. 2020), while the analysis of the spectroscopic properties of the sample as a function of redshift, black-hole mass, and luminosity was presented in Trefoloni et al. (2024). Sacchi et al. (2022) specifically focused on the individual spectroscopic analysis of the z > 2.5 quasars. All these studies found that, when typical blue quasars are selected, the X-ray and UV continuum emission and the overall spectral properties are remarkably similar regardless of redshift. Moreover, when the X-ray and UV continuum luminosities are determined from a one-by-one analysis, the dispersion in the L_X − L_UV relation narrows down significantly, further hinting at an underlying causal connection.

Here, we concentrate on a quantitative discussion of the second point, and on the general guidelines for the third point. Our main goals are to uphold the use of quasars as standardisable candles in the light of the criticism recently advanced in some papers, and to lay down a set of conservative rules that should be followed for a safe use of quasars in cosmology, avoiding possible biases or, even worse, incorrect results. Specifically, the essential issues we want to address are the following:

• The Hubble diagram of quasars, if not analysed in combination with supernovae, prefers a value of the present total matter density parameter of Ω_M > 0.6, which is significantly higher than the concordance value of Ω_M ≃ 0.3 obtained through supernovae and BAO (see e.g. Bargiacchi et al. 2022). This could suggest a disagreement between the Hubble diagrams of quasars and supernovae (Khadka & Ratra 2020).

• When the X-ray to UV relation is studied in narrow redshift intervals to factor out any cosmological assumption, the finite width of the bins can introduce biases in the determination of the correlation parameters (Petrosian et al. 2022).

• A discrepancy on slope and intercept values of the X-ray to UV relation emerges when the Hubble diagram of quasars is studied assuming different cosmological models (including ΛCDM and several of its extensions). This is also observed when the L20 sample is split in two subsamples at low and high redshift (with a separation value of z ≃ 1.5; see Khadka & Ratra 2020, 2021). These results have been ascribed to the non-standardisability of the quasar sample, and may stem from a redshift evolution of the relation. Some groups have proposed a correction for the possible evolution of the luminosities in the cosmological fits (e.g. Singal et al. 2022; Dainotti et al. 2022, 2023; Wang et al. 2022, 2024).

While we thoroughly discuss all these concerns in the following sections, it is also worth noting that other groups have instead tried to improve the precision on the determination of the cosmological parameters resorting to our technique (e.g. Melia 2019; Colgáin et al. 2024). Indeed, despite the limitation represented by their large dispersion in the Hubble diagram (∼1.4 dex) with respect to supernovae (< 0.1 dex), quasars are complementary to the latter as cosmological probes: they come in much greater numbers and they offer a unique opportunity to explore the very first billion years of the Universe.

4. Comparison of the Hubble diagrams of quasars and type Ia supernovae

A quantitative comparison of quasars and supernovae as distance indicators is not straightforward. If one derives the luminosity distance, D_L, of quasars from Equations (2) and (3), the values of D_L will depend on the parameters γ and β (see Section 8.3). The parameter γ controls the shape of the Hubble diagram, while the value of β is degenerate with the absolute scale of D_L, which can be obtained only through the cross-calibration with other standard candles in the common redshift range. This is why quasars should always be cross-calibrated with supernovae to extend the distance ladder in a similar way to what is done for supernovae with Cepheids.

Crucially, the intrinsic value of the slope γ can be estimated independently from any cosmological fit by analysing the X-ray to UV relation in small redshift bins. We adopted this procedure in several previous papers (e.g. Risaliti & Lusso 2015, 2019; Lusso et al. 2020) to test whether any evolution of the slope of the relation with redshift can be detected. We performed this test by analysing the F_X − F_UV relation in narrow redshift intervals, such that the differences in distance are smaller than the dispersion of the relation within each bin. We demonstrate that this condition is met in Section 5.3; here, we just want to analyse the consistency between the Hubble diagram of supernovae and the one of quasars in the redshift range in which both populations are present. We can perform this test over the redshift interval 0.7–1.6, which contains 1157 quasars; that is, ∼50% of the whole sample of 2421 objects in L20.

We have fitted the F_X − F_UV relation in narrow redshift bins, and the best-fit slope and intrinsic dispersion from the regression analysis are plotted in Figure 1. As is evident from this figure, there is no redshift evolution, where this is intended as a systematic trend and not just as an empirical scatter of the slope: a linear fit as a function of redshift, of the form γ = m z + q, returns m = 0.032 ± 0.067 and q = 0.553 ± 0.075. The statistical analysis of the quasar sample at low redshift demonstrates that the γ parameter does not show any statistically significant evolution with redshift. This result does not change by using a smaller bin width (at the price of reducing the statistics within each bin) or luminosities instead.

Fig. 1. Best-fit slope and intrinsic dispersion from the regression analysis of the FX − FUV relation in narrow redshift intervals from z = 0.7 to z = 1.6. The solid red line represents γ = 0.6. A linear fit of γ as a function of redshift in the form γ = m z + q gives m = 0.032 ± 0.067 and q = 0.553 ± 0.075 (dash-dotted blue line).

We remind the reader that we have a decent number of sources at z < 0.7 (398 AGN in L20), but these low-redshift quasars are neglected as their UV flux measurements could still be biased by the presence of additional contribution from the host galaxy in the optical. The L_UV values for the L20 quasar sample are computed from the photometric SEDs (see their Section 3); in the case of low-luminosity AGNs whose contrast with the galaxy is modest, the overall continuum flattens, and thus L_UV for these galaxies is overestimated. The slight shift of the low-redshift (low-luminosity) AGNs to higher L_UV values for a given X-ray luminosity results in a steepening of the best-fit γ in the L_X − L_UV plane. This is shown in Figure 2, in which we performed the same analysis as above but including all the low-redshift AGNs in the sample. The slope begins to steepen at z < 0.7, mimicking a spurious evolution of this correlation parameter with redshift (see e.g. Li et al. 2024).

Fig. 2. Best-fit slope and intrinsic dispersion from the regression analysis of the FX − FUV relation in narrow redshift intervals from z = 0.2 to z = 1.6. A linear fit of γ as a function of redshift in the form γ = m z + q gives m = −0.052 ± 0.018 and q = 0.673 ± 0.026 (dash-dotted blue line). Host-galaxy contamination at z < 0.7 mimics a spurious redshift evolution of γ. For a comparison, we overplot the resulting fit (magenta point) of the 13 local AGN in the redshift range 0.009 < z < 0.087, whose flux at 2500 Å is directly evaluated from their UV spectra.

For completeness, we also overplot the resulting fit of the 13 local AGN (0.009 < z < 0.087) in the L20 sample. Their γ = 0.52 ± 0.27 is consistent within the (large) uncertainties with the expected value of 0.6. The determination of the rest-frame 2500-Å flux in the latter sources relies on the fit of the IUE spectra, which span a wavelength interval (1845–2980 Å) in which the additional contribution of the host galaxy is negligible. Unfortunately, the sample statistics of local AGNs with UV spectroscopy covering the rest-frame 2500 Å is not yet sufficient to obtain better constraints. Additional UV data of nearby AGNs are thus required to improve the determination of the correlation parameters in the redshift range that is key for the cross-calibration of quasars with supernovae.

Leveraging the above results, we can now test the consistency between the distances of supernovae and quasars. First, we can fit the quasar sample in the z = 0.7–1.6 interval with a standard flat ΛCDM model. To perform the fit, the slope, γ, and the intercept, β, of the relation are left as free parameters, together with Ω_M and an additional parameter, K, which is included in the definition of the distance modulus, DM:

$$\begin{array}{c}\hfill \mathrm{DM}=5log(D_L(z))+25+K,\end{array}$$(4)

with D_L in units of megaparsecs. In this way, we can test the cross-calibration of quasars with supernovae in the common redshift interval and the consistency of quasar data with the concordance model below z = 1.6. This test is also very useful to make a direct comparison between the best-fit γ parameter and the value obtained from the analysis of the F_X − F_UV relation in narrow redshift bins described above.

The results of the simultaneous fit of quasars and supernovae are Ω_M = 0.284 ± 0.013, γ = 0.605 ± 0.015, β = 8.11 ± 0.46, and K = 0.005 ± 0.007. The latter parameter is actually redundant and expected to be relatively small, as the cross-calibration was performed through β, but it was included to be fully conservative. The parameter β was determined by requiring a superposition of the Hubble diagram of quasars to that of supernovae in the common redshift range. The best-fit slope is consistent with the value obtained from the flux–flux analysis of the X-ray to UV relation discussed above, further demonstrating that our assumption of a constant γ is robust. Figure 3 presents the Hubble diagram of quasars and supernovae together with the resulting best cosmological fit to the data. This figure also shows the averages of the DM values in narrow redshift bins for quasars only. We emphasise that such averages are not used to perform any cosmological fit, and are presented only for ease of visualisation in the Hubble diagram plots. The comparison of the Hubble diagram of quasars with that of supernovae shows that the consistency between these two data sets is remarkable; in other words, the shape of the quasar Hubble diagram is in very good agreement with that of supernovae.

Fig. 3. Hubble diagram of Pantheon+ supernovae at z < 1.6 (magenta points) and quasars in the redshift range 0.7 < z < 1.6 (golden points; the average 1σ uncertainty on the DM measurements is also shown in the bottom right corner). The blue points are averages of the DM values in narrow redshift bins for quasars only, and are shown for an easier visual comparison. The solid line is the best fit of a flat ΛCDM model to supernovae and quasars.

A similar consistency test can be carried out by considering a cosmologically independent, cosmographic model. The fit of the Hubble diagram was again performed simultaneously for the quasar and supernova data sets in the redshift interval of overlap. We adopted a cosmographic third-order log-polynomial model as in Bargiacchi et al. (2021), leaving the γ and β parameters free to vary as usual. A third-order polynomial is a conservative choice, as any one of second or higher order would converge rather fast given the limited redshift range. Since the statistical weight of supernovae is largely dominant at z < 1.6, the shape of the cosmographic model (i.e. the best-fit coefficients of the third-order log-polynomial) will virtually adapt to that data set only with a negligible contribution of quasars, whilst the parameters β and γ will be predominantly set by the quasar measurements instead. Specifically, the intercept β will provide the global calibration between the two data sets and the slope γ will adapt its best-fit value to obtain the best possible correspondence between quasars and supernovae. The result of the fit is γ = 0.601 ± 0.015 and β = 8.23 ± 0.46. Again, the slope is in excellent statistical agreement with the cosmology-independent analysis in small redshift bins discussed above. We have also applied a sigma-clipping technique to the quasar data to test whether outliers (at > 3σ, thus excluding 12 quasars) could introduce any difference in the measurements. The results on γ and β are fully consistent within the uncertainties (i.e. γ = 0.599 ± 0.014 and β = 8.32 ± 0.45). This is a further proof of the excellent coincidence between the two data sets in the common redshift range.

5. Testing the evolution of the correlation parameters at high redshift

We now focus on the cosmological analysis of the Hubble diagram for quasars only, aiming to address the case of an evolution of the normalisation, β, and of the slope, γ, in the redshift domain inaccessible to supernovae. We do so in three steps. We first show that in our method a complete degeneracy is present between a redshift evolution, β(z), and the shape of the distance-redshift relation. As a consequence, contrary to the case of γ, the issue of a redshift evolution of β can only be evaluated in light of our current understanding of the X-ray to UV relation itself. We then provide several observational results supporting the hypothesis of no redshift evolution of the correlation parameters. Finally, we repeat the analysis of the F_X − F_UV relation extending to z > 1.6 and adopting different sizes of redshift bins.

5.1. Degeneracy between the distance-redshift and β-redshift relations

We start by arguing that an intrinsic degeneracy exists between the distance-redshift relation and a possible evolution with redshift of the intercept, β, and that this degeneracy is present irrespective of the cosmological model. Based on Equation (3), it is straightforward to demonstrate that, if we assume a ‘true’ cosmological model and a non-evolving X-ray to UV relation as in Equation (1), such a relation can be reproduced within any alternative cosmological model, provided that the parameter β evolves as

$$\begin{array}{c}\hfill \beta (z)=\beta -(2\gamma -2)log\frac{D_L^{\ast }(z)}{D_L(z)},\end{array}$$(5)

where D_L(z) and D_L^*(z) are the luminosity distances computed with the true model and with the alternative one, respectively. This formal equivalence is an intrinsic feature of the standard-candle method. Since it holds precisely at all redshifts, it implies that if we assume an analytic function, β(z), with enough flexibility to reproduce the redshift-dependent term in Equation (4), it is impossible to disentangle the correct cosmological model from a fitting procedure of observational data, even for perfect standard candles: all the fits with any possible cosmological model will always produce the same likelihood value. Besides the formal correctness of the above statement, it is easy to verify that, even considering cosmological models very different from each other, the redshift-dependent term in Equation (4) typically has a simple shape that can be readily reproduced by an ordinary analytical function containing only a couple of free parameters.

To further elucidate this point, we show the effects of the degeneracy using the Pantheon+ supernova sample, which is widely accepted as a reliable distance indicator. We remind the reader that within our method, the quantity that contains the cosmologically relevant information is the parameter $\widehat{\beta }$ in Equation (3); that is, the absolute calibration of the relation in flux units at a given redshift. In this respect, assuming a redshift evolution of β in our analysis is formally equivalent to assuming a redshift dependence of the zero-point reference magnitude of type Ia supernovae; discussing whether these assumptions are justified is beyond the scope of this work, but we want to elaborate on their formal consequences. We first fitted the Pantheon+ sample in the standard way; that is, assuming the distance moduli as estimated by Scolnic et al. (2022) and adopting a flat ΛCDM model. We obtained the well-known results presented in Brout et al. (2022), and we computed the dispersion from the best fit as a function of redshift. The results are shown in Figure 4. We then repeated the fit assuming a fixed ΛCDM model with Ω_M = 1 and Ω_Λ = 0, and an evolution of the supernova zero-point magnitude of the form a₁exp(a₂z), with a₁ and a₂ as free parameters. Despite the rather extreme values chosen for the cosmological model, we were able to obtain a fit perfectly equivalent to the first one, which is also shown in Figure 4. The best-fit values of the free parameters are fully reasonable (a₁ ∼ −0.3; a₂ ∼ −0.5).

Fig. 4. Top panel: Hubble diagram of the Pantheon+ supernovae sample, fitted with an open ΛCDM model with ΩM and ΩΛ as free parameters, and with a ΛCDM model with fixed ΩM = 1 and ΩΛ = 0 and allowing for a redshift evolution of the absolute zero-point magnitude of supernovae. Bottom panels: Dispersion with respect to the best fit for the two models.

We conclude that there is no way to rule out the latter scenario based on the cosmological analysis, although we had purposely chosen an unreasonable cosmological model. It is obviously possible to dismiss the matter-only model based on other cosmological probes, yet we can safely discard this solution in the context of an analysis solely based on supernovae simply because we believe we have enough observational evidence and physical understanding about supernovae Ia to rule out the luminosity evolution implied by the best-fit values of the parameters a₁ and a₂. This exercise shows that any departure of the assumed cosmological model from the correct (and unknown) one can be compensated for by a suitable, albeit artificial evolution of some physical property of the standard candle at issue. For this reason, one cannot rely on a cosmological analysis to determine the ultimate standardisability of a given data set. The analogy with quasars is straightforward: it is not possible to either reject or confirm any evolution of the absolute calibration of the X-ray to UV relation based on a cosmological analysis. Therefore, the reliability of the method entirely depends on the observational evidence supporting the non-evolution hypothesis.

5.2. Observational evidence for a non-evolving X-ray to UV relation

While the analysis of the F_X − F_UV relation in narrow redshift intervals confirms the redshift-independence of γ, we do not have any such direct test for β. However, a first, fundamental argument against the possible evolution of β with redshift is the universal nature of black-hole accretion. It is well established that the accretion process is one and the same irrespective of black-hole mass, and that AGN are a scaled-up version of stellar-mass black holes (e.g. McHardy et al. 2006). Likewise, once a black hole starts to steadily accrete matter, and no other physical mechanisms are involved (e.g. strong outflows, or jets), there is no reason to believe that accretion physics should be any different at low and high redshifts. Indeed, no theoretical model has yet proposed, or even postulated, a variation in the accretion process only based on redshift. In this wake, there is ample indirect observational evidence that β must be independent of redshift:

1. We observe an excellent agreement between the Hubble diagrams of quasars and supernovae in the common redshift range. Despite the large scatter of quasars, there is neither a systematic offset nor any trend that may suggest a possible evolution of β, at least up to z ≃ 1.5. An evolving β at higher redshifts thus appears like a sort of fine-tuning required to reconcile the discrepancy observed in the Hubble diagram of high-redshift quasars with the ΛCDM predictions.

2. Δα_OX (i.e. the difference between the observed and the predicted α_OX) is independent of redshift. For any given luminosity, the parameter α_OX, defined as 0.384log(L_X/L_UV) after Tananbaum et al. (1979), has been found to be constant across the entire redshift range probed so far (see, for instance, Figure 4 in Vito et al. 2019).

3. The similarity between the average optical/UV spectrum of standard blue quasars at z ≲ 1.6 (e.g. Vanden Berk 2001) and that observed in the highest-redshift quasars (z ≳ 7; see e.g. Mortlock et al. 2011; Bañados et al. 2018) is striking. Such a spectral homogeneity implies that the underlying physical process is the same across a wide redshift range. This is quantitatively confirmed by the analysis presented in Trefoloni et al. (2024), in which we reveal no evolution of the quasar properties (e.g. the nuclear continuum slope) when the sources are selected to be blue, with little radio emission and no broad absorption lines, ensuring the absence of strong jets and/or outflows.

4. We individually analysed the X-ray and UV spectra of the 130 quasars at z > 2.5 in L20 in Sacchi et al. (2022), finding that there is no difference in either band between the nuclear continuum and overall spectral properties of high-redshift quasars and those of their low-redshift analogues.

5. The correlation between L_X and L_UV remains tight, if not slightly tighter, even when different proxies for the X-ray and UV emission are adopted (Signorini et al. 2023). When the residual contributions to the observed dispersion are investigated (i.e. those that cannot be removed through a careful sample selection, like variability, or accretion disc inclination), the intrinsic dispersion in the X-ray to UV relation turns out to be < 0.06 dex (Signorini et al. 2024).

All these pieces of observational evidence support the non-evolution of the γ and β parameters with redshift in typical quasars, in keeping with the notion that accretion physics is the same at low and high redshifts. Even so, we always need to check whether the F_X − F_UV relation for the adopted quasar sample shows any trend with redshift before performing any cosmological analysis. This issue was addressed in depth in previous works from our group (see Section 8 in L20); yet, below we further demonstrate that our method to test for the redshift evolution (or lack thereof) of the correlation parameters is robust.

5.3. Non-evolution of γ at z > 1.6

Following from Section 4, we divided the sample of quasars in narrow redshift bins within the interval 0.7–3.5, with a logarithmic step whose width is set by two criteria: each redshift interval must be characterised by enough statistics (the number of objects is > 15 in each bin), and by a dispersion in distances smaller than the dispersion in the F_X − F_UV relation. As we have already mentioned, this test serves not only to verify that γ does not evolve with redshift, but also to ascertain its intrinsic value, independent from any cosmological model. Figure 5 shows the slope resulting from the regression analysis of the F_X − F_UV relation computed in narrow redshift bins as a function of redshift. The bottom panel shows the values for the dispersion (δ) along the best fit of the F_X − F_UV relation compared to the dispersion of the luminosity distances; that is, δ_DL = σ(D_L²)/⟨D_L²⟩, where σ(D_L²) and ⟨D_L²⟩ are, respectively, the standard deviation and the mean of the distribution of the squared distances within each bin, which have been considered to be more conservative. The dispersion associated with the distance distribution is indeed much smaller than the overall dispersion of the F_X − F_UV relation in each bin².

Figure 6 presents the redshift evolution of the slope of the F_X − F_UV relation adopting a finer binning. As was expected, the average slope is fully consistent with the one shown in Figure 5, although the uncertainties are larger as the number of sources within each bin is much smaller, especially at z > 2.5. Therefore, the average γ value is not dependent upon the binning choice, as long as the distance dispersion within the bin is such that δ_DL < 0.1; that is, lower than the minimum observed dispersion reasonably expected for the X-ray to UV relation (cf. Sacchi et al. 2022). Here, we also show the resulting $\widehat{\beta }$ as a function of redshift, which displays the expected evolution of the luminosity distance with redshift expressed in Equation (3). We fitted the data points with a second-order cosmographic fit of the form

Fig. 5. Redshift evolution of the slope of the FX − FUV relation in narrow redshift bins (top panel). To perform the regression fit, X-ray and UV fluxes are normalised to 10−31 and 10−27 erg s−1 cm−2 Hz−1, respectively. The data points are plotted at the average redshift within the interval. Error bars represent the 1σ uncertainty on the mean in each bin. The solid and dashed grey lines are the mean and 1σ uncertainty range on the slope (γ). The red line marks γ = 0.6. The bottom panel shows the dispersion (δ, blue filled circles) along the best fit of the FX − FUV relation and the dispersion of the distance distribution (δDL ≃ 0.06, square symbols) in each bin. See Section 5.3 for details.

Fig. 6. Same as for Figure 5 but with a much finer binning within the redshift range 0.7–3.5. The bottom panel presents $\widehat{\beta }$ from Equation 3 as a function of redshift. The solid black line is a second-order cosmographic fit of the data, for which a Gaussian prior is considered on both γ and β (i.e. γ = 0.600 ± 0.005 and β = 8.00 ± 0.25). The dash-dotted line is a flat ΛCDM model with ΩM fixed to 0.3, γ and β fixed to 0.6 and 8.0, and normalised to the cosmographic fit at z = 1.

$$\begin{array}{c}\hfill D_L=\frac{cln(10)}{H_0}(x+a_1{x}^2),\end{array}$$(6)

where x = log(1 + z), a₁ is the polynomial coefficient, c is the speed of light, and H₀ is fixed to 70 km s⁻¹ Mpc⁻¹. The second-order cosmographic fit is presented for visualisation purposes only, as $\widehat{\beta }$ is not calibrated; hence, the best-fit values of the model have no cosmological meaning. A Gaussian prior was considered for both γ and β in the cosmographic fit (i.e. γ = 0.600 ± 0.005 and β = 8.00 ± 0.25). We also plot a flat ΛCDM model with Ω_M fixed to 0.3, and γ and β fixed to 0.6 and 8.0, respectively. This model is normalised to the cosmographic fit at z = 1. This (qualitative) comparison shows that, whilst there is consistency of the distances below z ≈ 1.7, the data points deviate from the predictions of a flat ΛCDM model at higher redshifts.

Again, such a discrepancy can be interpreted in two alternative ways: since we have demonstrated that γ does not show any systematic trend, either there is an evolution of β with redshift, or a flat ΛCDM model is not a good representation of the data. However, the Hubble diagram of quasars is in excellent agreement with that of supernovae in the common redshift range (see Figure 3), hence also β is constant at least up to z ∼ 1.5. Moreover, all the physical and spectral properties of the quasars selected at z > 2.5 do not differ from those of lower redshift AGN (i.e. Sacchi et al. 2022; Trefoloni et al. 2024). We therefore conclude that this deviation is most likely a limitation of the flat ΛCDM model.

6. The degeneracy between evolutionary and cosmological models

We now expand on the analysis presented in Section 5.1 and run a series of simulations to illustrate that it is not possible to draw any conclusions regarding the reliability of our technique purely based on the comparison of cosmological models (e.g. Khadka & Ratra 2020, 2021, 2022), even when accounting for the possible evolution of the correlation parameters γ and β with redshift. For such an approach to work, one would have to know the ‘true’ cosmological model, whereas even in the case of the widely accepted ΛCDM parameterisation, it has been demonstrated by several groups that this model has pitfalls. The logical fallacy lies in the fact that a cosmological model (or a suite of very different ones) cannot be employed to validate the standardisability of a specific cosmological probe, or data set. The determination of the correlation parameters γ and β can thus be obtained only from observational evidence, independently of any cosmological model, as we have discussed in Section 5.

6.1. The mock quasar sample

We simulated a quasar sample starting from the redshift and UV flux distributions of the L20 data set, where only z > 0.7 quasars are considered (i.e. 2023 quasars within the redshift interval 0.703 < z < 7.540, with an average redshift of ≃ 1.6). With respect to the observed distributions, we introduced a small random offset to the value of both redshifts and UV fluxes (in log units), assuming a Gaussian distribution centred at zero with dispersion σ_z = 0.001 and σ_FUV = 0.01. From the simulated redshift distribution we computed the luminosity distances adopting Ω_M = 0.1 and Ω_Λ = 0.9, whose values were chosen to be very different from the ones of the ΛCDM model or its simplest extensions. By combining the resulting distances and the simulated UV fluxes, we obtained mock UV luminosities, from which we subsequently derived the corresponding X-ray luminosities by assuming γ = 0.6, β = 8.0, and a dispersion δ of the L_X − L_UV relation of 0.24 dex (see e.g. Lusso & Risaliti 2016). We then converted X-ray luminosities back into fluxes, still in the framework of the above-mentioned cosmology. Uncertainties were assigned to both X-ray and UV fluxes by applying, again, to the observed ones a small random offset drawn from a Gaussian distribution. In order to evaluate the effects of sample statistics on the results, we considered two samples: one with exactly the same number of quasars as in the real data set (2023, sample 1), and another where we increased the statistics by a factor of 10 (i.e. 20230 sources, sample 2).

Before embarking upon the cosmological analysis, we have also examined whether the determination of the average correlation parameter γ in this clearly exotic cosmology presents any anomaly or bias. For extra caution, we adopted two different widths of the redshift bins. Figure 7 presents the evolution of the slope of the F_X − F_UV relation with redshift for the simulated quasar sample 1 (2023 sources). The plot in the left panel is obtained with the same number of intervals as in the real data (i.e. 24 intervals, δ_DL ≃ 0.06), while the plot in the right panel displays a coarser binning (i.e. 10 intervals, δ_DL ≃ 0.1). In both cases, we can retrieve the input γ value within uncertainties, thus confirming that we can safely assume the same γ for both the L_X − L_UV and the F_X − F_UV relations (i.e. γ_L = γ_F in the formalism of Petrosian et al. 2022) when applying this technique to the analysis of real data.

Fig. 7. Same as Figure 5 but for the simulated quasar sample 1 (2023 sources) with two different binning criteria: the one on the left is the same as for the real data (24 intervals, δDL ≃ 0.06), while the one on the right considers a coarser binning (10 intervals, δDL ≃ 0.1). In both cases, we can retrieve the input value of γ within uncertainties, supporting the assumption that γ is constant in the analysis of real data independently from any cosmological model (here ΩM = 0.1 and ΩΛ = 0.9 were assumed).

6.2. Results for spatially flat ΛCDM models

We first fitted both sample 1 and sample 2 with a ΛCDM model in which Ω_M and Ω_Λ were free to vary (model 1). The curvature Ω_k was set to zero (i.e. spatially flat case) and the radiation density is fixed to the observed value Ω_r = 2.47 × 10⁻⁵ h⁻², with a Hubble parameter h = 0.7 km s⁻¹ Mpc⁻¹. We also imposed the following two conditions: Ω_r = 1 − Ω_M − Ω_Λ and Ω_Λ < 0.99 + 1.8 Ω_M, the latter aimed at excluding solutions with no Big Bang. The results of the fit are shown in Figure 8. The marginalised one-dimensional best-fitting parameters with 1σ confidence intervals are summarised in Tables 1 and 2 for samples 1 and 2, respectively. As was expected, we can retrieve the input values of both the cosmological model and the correlation parameters within the uncertainties, with just a marginal offset of Ω_Λ for sample 1. This effect is due to the low sample statistics, especially at z < 1, where the energy density dominates. When sample 2 is considered, the value of Ω_Λ is better constrained thanks to the larger sample statistics in the redshift range 0.7–1.6.

Fig. 8. Results of the cosmological fit for both sample 1 (green) and sample 2 (black) for a spatially flat ΛCDM with ΩM and ΩΛ free to vary. The correlation parameters γ and β (as well as the dispersion δ) of the X-ray to UV relation are also free to vary (see Section 6.2 for details). Contours are at 68% and 95% confidence levels. The input values of both the cosmological and correlation parameters are marked by the red square. All the input parameters considered to build the mock quasar samples are retrieved within the uncertainties.

We then fitted both samples with the same cosmological model but with the matter and energy density parameters fixed to the values Ω_M = 0.3 and Ω_Λ = 0.7 (model 2), and Ω_M = 0.9 and Ω_Λ = 0.1 (model 5). The results of the relative cosmological fits are shown in Figures 9 and 10 for samples 1 and 2, respectively. The correlation parameters γ and β deviate from the input values at more than the 5σ statistical level, and this effect is even more significant when Ω_M and Ω_Λ are fixed to values that substantially differ from the ‘true’ ones, as is the case for Ω_M = 0.9 and Ω_Λ = 0.1.

Fig. 9. Results of the cosmological fit for sample 1 in a spatially flat ΛCDM model with matter and energy density parameters fixed to the values ΩM = 0.3 and ΩΛ = 0.7 (red, model 2 in Table 1) and ΩM = 0.9 and ΩΛ = 0.1 (blue, model 5 in Table 1). Contour levels are at 68% and 95% confidence. The input values of correlation parameters are marked by the black square. See Section 6.2 for details.

Fig. 10. Same as in Figure 9 for sample 2.

We finally fitted samples 1 and 2 by allowing also for an evolution with redshift of both slope and intercept of the X-ray to UV relation, of the form: γ(z) = γ₀ + γ₁(1 + z) and β(z) = β₀ + β₁(1 + z). Figures 11 and 12 show the results of this test. For sample 1, all the correlation parameters are in broad statistical agreement with the true values irrespective of the assumed cosmological parameters, which are fixed to the values Ω_M = 0.3 and Ω_Λ = 0.7 (model 4) or Ω_M = 0.9 and Ω_Λ = 0.1 (model 6). The evolution parameters γ₁ and β₁, in fact, are consistent with zero within ∼ 1.5σ at most. However, when we appreciably increase the sample statistics, the farther the cosmological model from the true values, the more significant the emergence of an evolution of the correlation parameters (Table 2, Figure 12), even though the input mock sample was built without assuming any evolution. This effect is clearly due to the fact that the correlation parameters try to compensate for the deviation of the cosmological model at issue from the input one we had assumed as ‘true’.

Fig. 11. Results of the cosmological fit for sample 1 in a spatially flat ΛCDM with matter and energy density parameters fixed to the values ΩM = 0.3 and ΩΛ = 0.7 (red, model 4 in Table 1), and ΩM = 0.9 and ΩΛ = 0.1 (blue, model 6 in Table 1). The slope and intercept are assumed to evolve with redshift as γ(z) = γ0 + γ1(1 + z) and β(z) = β0 + β1(1 + z). Contours are at 68% and 95% confidence levels. The input values of all the correlation parameters are marked by the black square. See Section 6.2 for details.

Fig. 12. Same as in Figure 11 but for sample 2.

Supposing that the input cosmological model considered to build the mock quasar sample is unknown, and that the adopted working model is very different from the ‘true’ one, any inconsistency between the data and the current model may lead to the wrong conclusion that there is a standardisability problem with the data, rather than a defect in the model that precludes a good representation of the data. This is especially true when a potential evolution of the correlation parameters is taken into account: it is not possible to distinguish between the non-standardisability of the data and a limitation of the model itself.

All these tests support the argument that cosmological models cannot be employed to draw any conclusion on the validity of any given observed data set for cosmological purposes. Along these lines, the now widely accepted accelerated expansion of the Universe could have been seen as a standardisability issue with supernovae. In the case of blue quasars, a mere comparison of the values obtained for the correlation parameters, γ and β, by fitting any quasar sample (or subsample) within different cosmological models is intrinsically flawed. The results of such an experiment do not (and cannot) entail any cosmology dependence of γ and β.

7. The evolution of X-ray and UV luminosities

Some authors (e.g. Dainotti et al. 2022; Wang et al. 2022, 2024) have applied a correction for the evolution of both X-ray and UV luminosities with redshift, with the aim of retrieving the ‘intrinsic’, de-evolved correlation between the two luminosities. For instance, by analysing a sample of AGN from SDSS-DR7 (with an i-band magnitude ≤19.1) cross-matched with Chandra and XMM-Newton, Singal et al. (2022) determine that the best-fit intrinsic power-law correlation between the X-ray and UV luminosities should be much flatter than what has ever been observed in any AGN sample; that is, γ = 0.28 ± 0.03. This value results from the application of the Efron–Petrosian non-parametric method (Efron & Petrosian 1992), where the evolution of the luminosities with redshift is assumed in Singal et al. (Singal et al., see their Section 3.2 for details) to have the following form:

$$\begin{array}{c}\hfill L^{\prime }=L/g(z),\end{array}$$(7)

where the evolution function g(z) is defined as

$$\begin{array}{c}\hfill g(z)=\frac{{(1+z)}^{\eta }}{1+{\left(\frac{1+z}{Z_{\mathrm{cr}}}\right)}^{\eta }},\end{array}$$(8)

with Z_cr = 3.7, while η takes the values of 3.3 and 0.55 for the optical and the X-ray data, respectively.

In practical terms, this correction is largely alternative to allowing for a redshift evolution of β. The results are strongly dependent upon both the function and the values of the exponents assumed. Indeed, by applying the same technique but with a simpler evolution function, specifically g(z) = (1 + z)^η, with the exponent taking the values of 4.36 and 3.36 for the optical and the X-ray data, respectively, Lenart et al. (2023, see also Dainotti et al. 2022) and references therein) found that the intrinsic slope of the X-ray to UV luminosity relation is statistically consistent with 0.6 instead. The Efron–Petrosian test statistics applies to heavily truncated data, as those delivered by blind surveys. Neither the UV data from the SDSS, nor the X-ray ones after the Eddington-bias correction (see Section 8.2) can be plainly considered as such. Bryant et al. (2021) point out that the Efron–Petrosian test also depends in a critical way on the choice of the detection threshold. An inaccurate implementation of this method can lead to unreliable or biased conclusions, specifically favouring a stronger, when not spurious, redshift evolution. For these reasons, we believe that this kind of analysis, at this stage, is not a compelling argument against the standardisability of quasars.

8. General guidelines on how to use quasars for cosmology

In this section, we outline some general guidelines that should be followed to safely use quasars for cosmology in the framework of the X-ray to UV correlation, through the application of Equations (2) and (3).

8.1. Sample selection

To build a quasar sample that can be utilised for cosmological purposes, multi-wavelength data from optical to X-rays are required to derive the rest-frame 2-keV and 2500-Å fluxes. Radio-bright quasars and broad-absorption-line (BAL) sources should be removed to reliably compute X-ray and UV nuclear fluxes that are the direct output of the accretion process and are not affected by jets and/or outflows. Further care should be taken to select low-redshift quasars with minimal contribution from the host galaxy (otherwise a redshift cut at z ∼ 0.7 is a conservative, yet reasonable choice), and sources with negligible flux attenuation from gas and dust at both X-ray and UV energies. Starting from any parent sample, the main aim is to pinpoint the intrinsic, nuclear emission from the accretion disc and the X-ray corona, thus taking advantage of the universal nature of the black-hole accretion process.

Even so, the large dispersion characterising the Hubble diagram of quasars (∼1.4 dex) compared to that of type Ia supernovae (< 0.1 dex) still represents a major limitation for the use of quasars as cosmological probes. To overcome this large dispersion, some authors suggest to select for cosmological fits quasar subsamples that have a smaller scatter along the X-ray to UV relation, through the application of severe sigma-clipping (down to a threshold of 1.5; Dainotti et al. 2024a,b). However, this procedure is not based on any physical argument, as the parent sample should already be the outcome of a homogeneous selection. Moreover, also the statistical grounds are questionable, since the choice of a tight sigma-clipping threshold has the only effect of cutting the wings of both the X-ray and UV flux distributions. We believe instead that the best way forward is to reduce the overall dispersion (e.g. Signorini et al. 2024), which requires more effort to better understand the physics behind the X-ray to UV relation and to improve the data quality of the observations.

8.2. Correction for the Eddington bias

Any flux-limited sample is biased towards brighter sources at high redshifts, and this effect is more relevant to the X-rays since the relative observed flux interval is narrower than in the UV. To minimise this bias, one possible approach is to neglect all the X-ray detections below a certain threshold, defined as κ times the intrinsic dispersion of the F_X − F_UV relation (δ) computed in narrow redshift intervals (Lusso & Risaliti 2016; Risaliti & Lusso 2019). Specifically:

$$\begin{array}{c}\hfill log{F}_{2\mathrm{keV},\exp }-log{F}_{2\mathrm{keV},\min }<\kappa \delta ,\end{array}$$(9)

where F_{2 keV,  min} represents the flux limit of a given observation or survey at the rest-frame 2 keV, whilst the product κδ takes a value that should be estimated for the considered quasar sample. The parameter F_{2 keV,  exp} is the monochromatic flux at 2 keV expected from the observed rest-frame quasar flux at 2500 Å, with the assumption of an intrinsic value for γ of 0.6, and it is calculated as follows:

$$\begin{array}{c}\hfill log{F}_{2\mathrm{keV},\exp }=(\gamma -1)log(4\pi )+(2\gamma -2)log{D}_L+\gamma log{F}_{\mathrm{UV}}+\beta ,\end{array}$$(10)

where D_L is the luminosity distance calculated for each redshift with a fixed cosmology (the results do not depend upon the choice of any reasonable model), and the parameter β represents the pivot point of the non-linear relation in luminosities, β = 26.5 − 30.5γ ≃ 8.0. The Eddington bias can be then reduced by including only X-ray detections for which the minimum detectable flux F_{2 keV,  min} in a given observation is lower than the expected X-ray flux F_{2 keV,  exp} by a factor proportional to the dispersion in the F_X − F_UV relation in narrow redshift bins (see Appendix A in Lusso & Risaliti 2016 and Risaliti & Lusso 2019).

8.3. Cosmological analysis: The likelihood

Fitting fluxes with a Gaussian prior on both γ and β should be the best option under the assumption that the intrinsic values for the slope and intercept do not vary; that is, the L_X − L_UV is universal and redshift-independent (see Section 5). An additional parameter should always be considered to perform the cross-calibration with supernovae, as in Equation (4). Nonetheless, it is good practice to fit the quasar sample with variable correlation parameters to verify that their values are in agreement with the assumed ones (i.e. γ ≃ 0.6 and β ≃ 8.0). Moreover, by fitting the data with a variable slope, the resulting uncertainties on the output parameters are more conservative than those obtained in the case of fixed correlation parameters. The likelihood, when fluxes are fitted, can be written as follows:

$$\begin{array}{c}\hfill ln\mathcal{L}=-\frac{1}{2}{\displaystyle \sum _{i=1}^n}[\frac{{(log{F}_{\mathrm{X},i}^{\mathrm{obs}}-log{F}_{\mathrm{X},i}^{\mathrm{mod}})}^2}{{(dlog{F}_{\mathrm{X},i}^{\mathrm{obs}})}^2+e^{2ln\delta }}+ln[{(dlog{F}_{\mathrm{X},i}^{\mathrm{obs}})}^2+e^{2ln\delta }]],\end{array}$$(11)

where

$$\begin{array}{c}\hfill log{F}_{\mathrm{X},i}^{\mathrm{mod}}=\beta +\gamma (log{F}_{\mathrm{UV},i}^{\mathrm{obs}}+27.5)-2(1-\gamma )(\mathrm{DM}-28.5)+(\gamma -1)log(4\pi ),\end{array}$$(12)

in which DM = DM(Ω_M, Ω_Λ, z, K) as given in Equation (4), and all the data units are in the cgs system.

Another possibility is to directly fit the distance moduli, in which case the likelihood in Equation (11) is modified as follows:

$$\begin{array}{c}\hfill ln\mathcal{L}=-\frac{1}{2}{\displaystyle \sum _{i=1}^n}[\frac{{({\mathrm{DM}}_i^{\mathrm{obs}}-{\mathrm{DM}}_i^{\mathrm{mod}})}^2}{{\left(d{\mathrm{DM}}_i^{\mathrm{obs}}\right)}^2+e^{2ln\delta }}+ln[{\left(d{\mathrm{DM}}_i^{\mathrm{obs}}\right)}^2+e^{2ln\delta }]],\end{array}$$(13)

where the luminosity distance is

$$\begin{array}{c}\hfill D_L=5\{\frac{1}{2(1-\gamma )}[\beta +\gamma (log{F}_{\mathrm{UV}}+27.5)-log{F}_{\mathrm{X}}]+28.1\},\end{array}$$(14)

in centimetres. To quantify the difference between these two approaches, we fit the L20 sample together with the Pantheon+ sample (Scolnic et al. 2022) in a spatially flat ΛCDM with Ω_M and Ω_Λ free to vary. The cosmological fit of quasars and supernovae considering fluxes has γ and β free to vary, whilst the cosmological fit considering luminosity distances assumed fixed values for γ and β (i.e. γ = 0.61 ± 0.02 and β = −31.16 ± 0.01, where the uncertainties are considered for the error propagation on the distances). As β is fixed, a cross-calibration parameter, K, is included in both cases to allow for extra flexibility. As was expected, the value of K is very close to zero, implying that the input value for β is almost correct.

In Figure 13, we present the comparison of the cosmological parameters, Ω_M, Ω_Λ, and K, which clearly shows that they are in statistical agreement within 1σ, with slightly larger uncertainties for the case of variable γ and β, as expected. The summary of the marginalised one-dimensional best-fitting parameters with 1σ confidence intervals is provided in Table 3.

Fig. 13. Results of the cosmological fit for the L20 quasar sample combined with Pantheon+ supernovae in a spatially flat ΛCDM with ΩM and ΩΛ free to vary. Orange contours: cosmological fit of quasars and supernovae considering fluxes, with γ and β free to vary. Violet contours: cosmological fit of quasars and supernovae considering luminosity distances, with γ and β fixed (i.e. γ = 0.61 and β = −31.16). A cross-calibration parameter K is also considered in both cases.

9. Conclusions

The aim of this work is to address the question of whether the non-linear relation between the X-ray and UV emission of quasars can be employed to derive quasar distances, and whether typical blue quasars can therefore be regarded as standardisable candles for cosmology. For such a technique to be reliable, the correlation parameters, γ (slope, in log space) and β (intercept), must be cosmology- and redshift-independent. As has already been shown in previous works of our group, the distances estimated with this method agree with the standard flat ΛCDM model up to z ∼ 1.5, but start to deviate significantly from the model predictions at higher redshifts. For this reason, some authors have suggested that the observed discrepancy could be due to heterogeneities between the low- and high-redshift subsamples considered in Lusso et al. (2020), or to an evolution of the relation with redshift, ultimately implying the non-standardisability of quasars. To refute the latter arguments, here we argue that the values of γ ≃ 0.6 and β ≃ 8.0 are inherent to the black-hole accretion process and consequential to its universal nature, and thus intrinsically cosmology-independent. Indeed, we have demonstrated that the relation between X-ray to UV fluxes does not show any redshift evolution in its slope when studied in narrow redshift bins, confirming that γ is both redshift- and cosmology-independent. Moreover, we performed a set of simulations showing that all the claimed inconsistencies naturally arise from any limitations of the cosmological model(s) used in the analysis, and do not necessarily indicate any problem with the data set itself. Specifically, the claimed non-standardisability of high-redshift quasars is an obvious manifestation of our ignorance of the true cosmological model. Our main results can be summarised as follows:

• The fit the F_X − F_UV relation in narrow redshift bins and the comparison with supernova-derived distances in the common redshift range show that quasars are fully consistent with supernovae and that there is no redshift evolution of the slope, γ, at z < 1.6. The statistical analysis of the quasar sample is then extended up to z ≃ 3.5, confirming that the γ parameter does not exhibit any systematic trend across the entire redshift range for which we have enough statistics to perform the regression fit. This result does not depend on the width of the adopted redshift bins, provided the dispersion of the luminosity distance within each bin is smaller than the one expected for the X-ray to UV relation. The deviation from the flat ΛCDM that we start to detect beyond z ≈ 1.6 cannot thus be attributed to an evolution of the slope, γ. While no direct test can be performed to also verify the non-evolution of β, its sudden change at z > 1.6 is not supported by observational evidence.

• Simulations of a mock quasar sample demonstrate that, when different cosmological models are compared, there is an unavoidable degeneracy between the detection of a possible evolution of γ and β (or a departure from their ‘intrinsic’ values), and the determination of the cosmological parameters. This is a general feature of the standard-candle method: any departure from the correct, unknown cosmology can be neutralised through a suitable modification or evolution of some parameter intrinsic to the physical sources in use, whether these are quasars or supernovae. Hence, no cosmological models or comparison thereof can conclusively validate or reject any given probe or data set.

We conclude that the reliability of our method is fully corroborated by the comparison with supernovae up to z ≃ 1.5, while outside the common redshift range it can only derive from a cosmology-independent evaluation of the hypothesis of non-evolution, especially for β, subject to a meticulous check of the sample selection and of the flux measurements for possible redshift-dependent systematic effects. Multi-wavelength observations of quasars across a wide redshift range show a high level of homogeneity in their spectral properties. As a matter of fact, when a supermassive black hole steadily accretes matter, the physics underlying this process is the same regardless of redshift; and so is its observational appearance in the X-ray and UV bands, as long as no additional energy-extraction mechanism (e.g. the launch of radio jets and/or fast uncollimated outflows) or contamination effects (e.g. host-galaxy emission, dust and/or gas attenuation) are present. Since there are no clear systematics in either the sample selection or in the flux measurements for all the quasars that we have employed in our works, we are confident the application of the X-ray to UV relation to cosmology is solid and well motivated. A better understanding of the physical process behind the L_X − L_UV relation could definitely strengthen the method, and more effort is required along this line. Future measurements of supernovae at redshifts of z > 2, possibly confirming the discrepancy with ΛCDM found with quasars, will provide independent observational proof that new physics is involved.

Acknowledgments

We thank the anonymous reviewer for their thorough reading and for useful comments. This work was performed in part at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-2210452. All the authors acknowledge partial support from the Simons Foundation (1161654, Troyer).

References

All Tables

All Figures

	Fig. 1. Best-fit slope and intrinsic dispersion from the regression analysis of the FX − FUV relation in narrow redshift intervals from z = 0.7 to z = 1.6. The solid red line represents γ = 0.6. A linear fit of γ as a function of redshift in the form γ = m z + q gives m = 0.032 ± 0.067 and q = 0.553 ± 0.075 (dash-dotted blue line).
In the text

Fig. 2. Best-fit slope and intrinsic dispersion from the regression analysis of the FX − FUV relation in narrow redshift intervals from z = 0.2 to z = 1.6. A linear fit of γ as a function of redshift in the form γ = m z + q gives m = −0.052 ± 0.018 and q = 0.673 ± 0.026 (dash-dotted blue line). Host-galaxy contamination at z < 0.7 mimics a spurious redshift evolution of γ. For a comparison, we overplot the resulting fit (magenta point) of the 13 local AGN in the redshift range 0.009 < z < 0.087, whose flux at 2500 Å is directly evaluated from their UV spectra.

In the text

Fig. 3. Hubble diagram of Pantheon+ supernovae at z < 1.6 (magenta points) and quasars in the redshift range 0.7 < z < 1.6 (golden points; the average 1σ uncertainty on the DM measurements is also shown in the bottom right corner). The blue points are averages of the DM values in narrow redshift bins for quasars only, and are shown for an easier visual comparison. The solid line is the best fit of a flat ΛCDM model to supernovae and quasars.

In the text

	Fig. 4. Top panel: Hubble diagram of the Pantheon+ supernovae sample, fitted with an open ΛCDM model with ΩM and ΩΛ as free parameters, and with a ΛCDM model with fixed ΩM = 1 and ΩΛ = 0 and allowing for a redshift evolution of the absolute zero-point magnitude of supernovae. Bottom panels: Dispersion with respect to the best fit for the two models.
In the text

Fig. 5. Redshift evolution of the slope of the FX − FUV relation in narrow redshift bins (top panel). To perform the regression fit, X-ray and UV fluxes are normalised to 10−31 and 10−27 erg s−1 cm−2 Hz−1, respectively. The data points are plotted at the average redshift within the interval. Error bars represent the 1σ uncertainty on the mean in each bin. The solid and dashed grey lines are the mean and 1σ uncertainty range on the slope (γ). The red line marks γ = 0.6. The bottom panel shows the dispersion (δ, blue filled circles) along the best fit of the FX − FUV relation and the dispersion of the distance distribution (δDL ≃ 0.06, square symbols) in each bin. See Section 5.3 for details.

In the text

Fig. 6. Same as for Figure 5 but with a much finer binning within the redshift range 0.7–3.5. The bottom panel presents $\widehat{\beta }$ from Equation 3 as a function of redshift. The solid black line is a second-order cosmographic fit of the data, for which a Gaussian prior is considered on both γ and β (i.e. γ = 0.600 ± 0.005 and β = 8.00 ± 0.25). The dash-dotted line is a flat ΛCDM model with ΩM fixed to 0.3, γ and β fixed to 0.6 and 8.0, and normalised to the cosmographic fit at z = 1.

In the text

Fig. 7. Same as Figure 5 but for the simulated quasar sample 1 (2023 sources) with two different binning criteria: the one on the left is the same as for the real data (24 intervals, δDL ≃ 0.06), while the one on the right considers a coarser binning (10 intervals, δDL ≃ 0.1). In both cases, we can retrieve the input value of γ within uncertainties, supporting the assumption that γ is constant in the analysis of real data independently from any cosmological model (here ΩM = 0.1 and ΩΛ = 0.9 were assumed).

In the text

Fig. 8. Results of the cosmological fit for both sample 1 (green) and sample 2 (black) for a spatially flat ΛCDM with ΩM and ΩΛ free to vary. The correlation parameters γ and β (as well as the dispersion δ) of the X-ray to UV relation are also free to vary (see Section 6.2 for details). Contours are at 68% and 95% confidence levels. The input values of both the cosmological and correlation parameters are marked by the red square. All the input parameters considered to build the mock quasar samples are retrieved within the uncertainties.

In the text

	Fig. 9. Results of the cosmological fit for sample 1 in a spatially flat ΛCDM model with matter and energy density parameters fixed to the values ΩM = 0.3 and ΩΛ = 0.7 (red, model 2 in Table 1) and ΩM = 0.9 and ΩΛ = 0.1 (blue, model 5 in Table 1). Contour levels are at 68% and 95% confidence. The input values of correlation parameters are marked by the black square. See Section 6.2 for details.
In the text

	Fig. 10. Same as in Figure 9 for sample 2.
In the text

Fig. 11. Results of the cosmological fit for sample 1 in a spatially flat ΛCDM with matter and energy density parameters fixed to the values ΩM = 0.3 and ΩΛ = 0.7 (red, model 4 in Table 1), and ΩM = 0.9 and ΩΛ = 0.1 (blue, model 6 in Table 1). The slope and intercept are assumed to evolve with redshift as γ(z) = γ0 + γ1(1 + z) and β(z) = β0 + β1(1 + z). Contours are at 68% and 95% confidence levels. The input values of all the correlation parameters are marked by the black square. See Section 6.2 for details.

In the text

	Fig. 12. Same as in Figure 11 but for sample 2.
In the text

Fig. 13. Results of the cosmological fit for the L20 quasar sample combined with Pantheon+ supernovae in a spatially flat ΛCDM with ΩM and ΩΛ free to vary. Orange contours: cosmological fit of quasars and supernovae considering fluxes, with γ and β free to vary. Violet contours: cosmological fit of quasars and supernovae considering luminosity distances, with γ and β fixed (i.e. γ = 0.61 and β = −31.16). A cross-calibration parameter K is also considered in both cases.

In the text