© The Authors 2024

1. Introduction

Quasars are the most luminous persistent sources in the Universe and, as such, they enable its investigation up to redshifts of z > 6 (e.g. Mortlock et al. 2011; Bañados et al. 2018; Yang et al. 2020; Wang et al. 2021; Zappacosta et al. 2023). Thanks to their ubiquity and high luminosity, these sources were quickly recognized as powerful cosmological probes. Since the 1970s, correlations between quasars’ observables have been employed to determine their distances. The discovery of the Baldwin effect (i.e. the emission-line equivalent widths anti-correlate with the quasar luminosity; Baldwin 1977) envisaged that quasars could be utilised as cosmological probes. Techniques involving similar relations have been subsequently developed (see Section 6 in Czerny et al. 2018 for a detailed review). Yet the large dispersion as well as the scarcity of suitable sources, have in most cases limited their applicability for measuring the expansion parameters (e.g. H₀). Other attempts to extract cosmological information from quasars are based on the broad line region (BLR; Elvis & Karovska 2002) or dust (Yoshii et al. 2014) reverberation mapping in combination with spectro-astrometric measurements (SARM; Wang et al. 2020), the Eigenvector 1 main sequence (Marziani & Sulentic 2014), or the X-ray excess variance (La Franca et al. 2014).

In the last decade, the efforts of our group have been mostly devoted to developing and refining a novel technique based on the non-linear relation between the X-ray and UV luminosity (L_X − L_UV) in quasars (or equivalently, the α_OX − L_UV relation; e.g. Vignali et al. 2003; Lusso et al. 2010), which has long been known (Tananbaum et al. 1979). The L_X − L_UV relation, according to our current knowledge, reflects the interplay between the UV photons coming from the disc (Shakura & Sunyaev 1973; Novikov & Thorne 1973) and the hot ‘corona’ (Haardt & Maraschi 1991, 1993) that boosts their energy producing the X-ray emission. Despite the theoretical efforts (e.g. Nicastro 2000; Merloni 2003; Lusso & Risaliti 2017; Arcodia et al. 2019), a comprehensive description of the physical mechanism underlying the L_X − L_UV relation is still missing. Even the choice of the monochromatic luminosities at 2 keV and 2500 Å can be regarded as rather arbitrary, being mainly due to historical and empirical reasons. Indeed, we have argued that better proxies, such as the Mg II emission line (see Signorini et al. 2023 for further details), of the physical mechanisms responsible for the observed relation can be exploited.

In addition to important implications for accretion physics, the non-linearity of the L_X − L_UV relation, combined with the non-evolution with redshift of its parameters (Risaliti & Lusso 2019), makes it a useful tool to determine the luminosity distance of quasars and thus transforms them into a new class of ‘standardizable candles’. The L_X − L_UV relation indeed proved to be a fruitful cosmological tool after it was demonstrated that most of the observed scatter is not intrinsic but rather due to observational issues, as the dispersion can be narrowed down to 0.09 dex in case of samples with high quality data (Sacchi et al. 2022). This new category of standard candles becomes particularly convenient when it comes to populating the intermediate redshifts between the farthest supernovae at z ∼ 2 (Scolnic et al. 2018) and the cosmic microwave background (CMB). Despite a good agreement with the ΛCDM model in the low-redshift regime of the Hubble–Lemaître diagram, a 3-4σ tension with the model, driven by the quasar data at earlier cosmic epochs (z > 1.5), has been reported (Risaliti & Lusso 2019; Bargiacchi et al. 2021; Sacchi et al. 2022; Bargiacchi et al. 2022; Lenart et al. 2023).

In Lusso et al. (2020, henceforth L20), we discussed how we produced our latest compilation of 2421 quasars for cosmological purposes, describing how the measurements of X-ray and UV luminosities were performed. However, in order to validate the selection criteria adopted to build the sample, we still need to address two main issues that could be present to some extent: i) an evolution with redshift of the slope of the optical and UV continuum and ii) the contamination of host galaxy in the optical and flux attenuation due to dust and gas. While we discuss the validity of the cosmological techniques adopted to analyse the sample in a companion paper (Risaliti et al., in preparation), in this work we mainly focus on a spectral confirmation of the robustness of the selection criteria. In particular, we aim at demonstrating, by means of composite spectra, that the continuum of our sample (in which the UV proxy L_2500 Å is evaluated) does not show any appreciable evolution with redshift. At the same time, the spectral composites are a useful tool to reveal the presence of residual (i.e., still present despite the selection cuts) dust extinction and/or host galaxy contribution.

To the end of investigating any possible spectral diversity among the objects constituting our cosmological sample, we employed stacked templates. The lack of a significant evolution of the continuum among the stacked spectra, representative of different regions of the parameter space, ensures the universality of the accretion mechanism powering the optical–UV emission.

In this paper, we analyse the average spectral properties of the UV and optical data of the cosmological sample described in L20, with three main goals: i) directly verifying how close the average spectral properties of our sample are to those of typical blue quasars, by means of composite spectra in different regions of the M_BH − L_bol − z parameter space; ii) estimating the degree of host galaxy contamination; iii) quantifying the effect of intrinsic dust absorption in our sample of quasars. The last two effects would influence our estimates of fL_2500 Å and therefore undermine the cosmological application of the L_X − L_UV relation. Employing spectral data to understand – and if necessary, account for – such possible observational effects, that can alter both X-ray and UV measurements would definitely corroborate the evidence that the observed discrepancy with the concordance cosmological model is not due to observational issues.

The paper is structured as follows. A brief summary about the selection of the quasar sample is provided in Section 2, whilst the analysis performed on the sample spectra is described in Section 3. Results are discussed in Section 4, and conclusions are drawn in Section 5.

2. The data set

The quasar sample analysed here was published by L20 and is made of 2421 sources within the redshift range 0.13 < z < 5.42. The data set was produced by merging seven different subsamples available from the literature.

The Sloan Digital Sky Survey data release 14 (SDSS-DR14, Pâris et al. 2018) was cross-matched with the XMM-Newton 4XMM-DR9 (Webb et al. 2020) catalogue of serendipitous sources, the second release of the Chandra Source Catalogue (CSC2.0, Evans et al. 2010) and the AGN sample described in Menzel et al. (2016), belonging to the XMM-Newton XXL survey (XMM–XXL, PI: Pierre), producing respectively the SDSS-4XMM, the SDSS-Chandra, and the SDSS-XXL samples.

Local quasars (0.009 < z < 0.1) with available UV data from the International Ultraviolet Explorer (IUE) in the Minkulski Archive for Space Telescopes (MAST) and X-ray data were included in the sample to anchor the normalisation of the quasar Hubble diagram with Type Ia supernovae.

To improve the coverage at high redshifts, three samples of quasars with pointed observations were added. We added the sample of 29 luminous quasars at redshift 3.0 < z < 3.3 with dedicated XMM-Newton X-ray observations described in Nardini et al. (2019), the 53 high redshift (z > 4.01) quasars from Salvestrini et al. (2019), and the 25 quasars at z > 6 from Vito et al. (2019).

In order to produce such a compilation, several selection criteria were adopted. We thus redirect the interested reader to L20 (in particular, their Section 5) for a more comprehensive description, while here we briefly highlight the filters adopted to define the sample. All the selected sources are radio quiet¹ and not flagged as broad absorption line (BAL) quasars according to the SDSS-DR14 quasar catalogue (Pâris et al. 2018), which was the latest release at the time of the L20 analysis. A broadband photometric spectral energy distribution (SED) was build for each source, and only quasars with UV-optical colours close to the ones of typical blue quasars (Richards et al. 2006) were retained. In addition, steep X-ray spectra (1.7 < Γ < 2.8) and deep enough X-ray observations (to avoid the Eddington bias), were required. Such properties make the L20 sample a truly unique sample of blue quasars, sharing a high degree of homogeneity in terms of broadband properties and observational data. The present work is focused on the spectral properties of the sources with available SDSS data, which constitute the bulk (∼93%) of the sample.

3. Methods

3.1. The parameter space

Our aim is to study whether the continuum properties of our sample evolve with redshift or any of the governing parameters of the quasar spectrum, namely the bolometric luminosity (L_bol), black hole mass (M_BH), and accretion rate (parametrised through the Eddington ratio λ_Edd = L_bol/L_Edd, where L_Edd is the luminosity at the Eddington limit). To this end, we considered the latest quasar catalogue by Wu & Shen (2022, henceforth WS22), who compiled continuum and emission-line properties for 750,414 broad-line quasars from the SDSS DR16, including physical quantities such as single-epoch virial M_BH, L_bol, and λ_Edd. We also gathered the latest spectral data available for each source in our sample, hence the spectra of ∼4% of the sources belong to the DR17 release, where calibration files have been updated (Abdurro’uf et al. 2022).

Fig. 1 presents the log L_bol – log M_BH plane for our sample. The partition of the plane was defined by seeking the balance between a large number of bins to improve the sensitivity on the parameters (L_bol, M_BH) that control the average spectrum, and a statistically robust number of objects in each region. We further divided the sample into five redshift intervals, z = [0.13, 0.77, 1.12, 1.52, 2.03, 5.42], each one containing a similar number of objects (roughly 480 objects per bin). The resulting grid is shown in Fig. 1. The width of the log M_BH (log L_bol) bins is 0.7 dex (1 dex), large enough to encompass the large systematic uncertainties associated with the estimate of this quantity, typically 0.3 dex (but up to 0.4 dex in case of C IV-based M_BH, e.g. Shen et al. 2011; Kratzer & Richards 2015; Rakshit et al. 2020). The points of the grid are log M_BH = [7.0, 7.7, 8.4, 9.1, 9.8, 10.5] and log L_bol = [44.4, 45.4, 46.4, 47.4, 48.4]. This procedure defines a 3D grid made of 5 × 5 × 4 = 100 bins. The individual cells are labelled using the following nomenclature: increasing capital letters (A, B, C, D) denote decreasing bolometric luminosity across the bins, while increasing numbers (1, 2, 3, 4, 5) mark increasing average black hole masses, and increasing lower-case letters (a, b, c, d, e) represent increasing redshift intervals, which are also colour-coded according to the values in the legend.

Fig. 1. Parameter space of log MBH – log Lbol for the quasar sample published in L20 with colour-coded redshifts. Dot-dashed lines represent constant λEdd values. Typical uncertainties on log(MBH) and log(Lbol) are represented in the bottom-right corner.

Clearly, not all the bins are densely, nor uniformly populated. Lower redshift quasars preferentially occupy the lower luminosity and black hole mass region, whereas higher redshift sources are observed at higher values of both L_bol and M_BH. Moreover, there are cells in the L_bol – M_BH plane where the barycentre of the sources is close to the edge of a given interval (e.g. B2, B5, C1). In these cases, low sample statistics would produce not only low signal-to-noise (SNR) composite spectra, but also a stack not fully representative of its nominal locus in the parameter space.

To select the regions in the log M_BH – log L_bol plane populated enough to be included in the stacking analysis, we quantified the probability that a statistical fluctuation in terms of log M_BH and log L_bol of objects in adjacent bins pushed them into a given region. We simulated 1000 mock samples starting from the actual distribution of log M_BH and log L_bol in each bin within the grid. The values of log M_BH and log L_bol are then scattered by a value drawn from a Gaussian distribution whose amplitude was set by the total uncertainty on each given data point. In doing so, we conservatively assumed the largest possible systematic uncertainty on both parameters and summed it in quadrature to the statistical uncertainty reported in WS22. For the bolometric corrections, we assumed a systematic uncertainty of ∼50% on log L_bol for each data point in the plane, suitable in case the bolometric correction is performed on the rest-frame 1350 Å luminosity (see Shen et al. 2011), which translates into ≲0.22 dex. The systematic uncertainty on log M_BH was assumed to be 0.40 dex, a value based on the C IV calibration (i.e. the least reliable one). As a consequence of the random scatter, the points of the mock sample distribution could escape the bin where the original data were instead contained. We then computed the average log M_BH and log L_bol of the resulting scattered points and retained only the regions within the parameter space where in at least 90% of the mock samples the average log M_BH and log L_bol remained within the bin. We refer to bins retained by this procedure as bins meeting the ‘representativeness criterion’. In general, 44 bins out of 100 in the parameter space do not contain any object, while 25 bins are excluded by the aforementioned criterion. The whole selection procedure ultimately excluded 69 scarcely populated bins (with less than nine objects), and it ensures that the composite spectrum is truly representative of its associated region in the parameter space. At the end of this procedure, 2316 objects were retained in the parameter space.

3.2. Building the stacks

The spectra were collected according to their PLATE, MJD and FIBERID reported in WS22². We produced the stacked spectra and adopted a bootstrap recipe for the evaluation of the uncertainty on the stacks (more quantitative assessments on how different stacking prescriptions affect the results are explored in Appendix A) in each of the selected regions following the steps described below.

Each spectrum was corrected for Galactic absorption assuming the value for the colour excess E(B − V) available in the WS22 catalogue (Schlafly & Finkbeiner 2011) and a Fitzpatrick (1999) extinction curve. Then, the de-reddened spectra were shifted to the rest frame.

All the spectra in the same redshift bin were resampled, by means of linear interpolation, onto a fixed wavelength grid. The grid was defined in the interval λ₁ − λ₂ with a fixed dispersion Δλ so as to cover the shortest and the longest rest-frame wavelengths in each redshift interval. The values for λ₁, λ₂, Δλ are listed in Table 1. The value for Δλ was chosen in order to preserve the SDSS observed frame resolution also in the rest frame spectra, similarly to what described in Section 3.2 of Vanden Berk (2001, VB01), assuming an observed frame spectral resolution of ∼2000 at 6500 Å.

We also applied a few quality cuts. Objects with more than 10% of the pixels flagged as bad (AND_MASK > 0) in their SDSS spectra were excluded. Also objects where the normalisation wavelength $\overline{\lambda }$ (the wavelength chosen to scale all the spectra in the cell) was not covered in the rest-frame spectrum were excluded.

The average flux in each spectral channel and the relative uncertainty were evaluated adopting a bootstrap technique. In brief, a number of spectra equal to the total number of objects in the bin is drawn, allowing for replacement. These spectra are normalised by their flux value at the reference wavelength $\overline{\lambda }$, and an average stack is produced. This procedure is repeated with new draws until the desired number of stacks is reached. Subsequently, in each spectral channel the distribution of the normalised fluxes is computed, a 3σ clipping is applied, and the median value (less affected by possible strong outliers) is chosen. The uncertainty on the median value was evaluated as the semi-inter-percentile range between 1^st and 99^th percentile.

The information about the relevant quantities adopted in each redshift interval is shown in Table 1. We explored different possible ways of combining the spectra, averaging the flux in each spectral channel and estimating the uncertainty. How these assumptions affect the resulting composites and the subsequent results, is presented in Appendices A and B. All the composite spectra are presented in Appendix G.

3.3. Spectral fit of the composites

To estimate the characteristic spectral quantities of the continuum (chiefly the slope α_λ of the power law f_λ ∝ λ^−α_λ) and the emission lines, in particular the rest-frame equivalent width (EW), the full width at half-maximum (FWHM), and the peak offset velocity (v_off), we performed the spectral fits of all of our composite spectra. To this end we adopted a custom-made code, based on the IDL MPFIT package (Markwardt 2009), which takes advantage of the Levenberg-Marquardt technique (Moré 1978) to solve the least-squares problem. Our main goal is a faithful reproduction of the continuum, in order to investigate whether the continuum slope and consequently the 2500 Å flux exhibits any significant evolution with the accretion parameters.

The baseline model is given by a continuum power law pivoting at different emission-free wavelengths depending on the redshift interval (3000 Å in the three bins at z< 1.52, 2200 Å in the two bins at z≥ 1.52), a set of emission lines, the Balmer continuum, and a set of iron templates convolved with a Gaussian line-of-sight velocity distribution whose weights are free parameters in the fit. Each line required a different decomposition: the C IV line was generally fitted with two Lorentzian profiles (one broad, one narrow). A single broad Lorentzian profile combined with the iron templates is able to adequately reproduce the Mg II profile. The Hβ was modelled by considering a core narrow Gaussian component whose FWHM could not exceed ∼1000 km s⁻¹ and a broad component whose analytical shape (Gaussian or Lorentzian) was chosen each time to best reproduce the data. The [O III] doublet was modelled with Gaussian or Lorentzian profiles with a ratio of the [O III]λ4959Å to λ5007Å components fixed at one to three. In addition, other weaker forbidden ([Ne V] λ3426, [O II] λ3727, [Ne IIIλ3869]) and permitted (Hγ, Hδ) emission lines were included to improve the residuals. Above log(L_bol/(erg s⁻¹)) ≳ 46 the host galaxy contribution at optical wavelengths is expected to be negligible with respect to the quasar emission (Shen et al. 2011; Jalan et al. 2023). Therefore, we only added an elliptical galaxy template from Assef et al. (2010) to the model, to account for possible galactic emission, only in the two lowest redshift bins, where the bulk of the sample has L_bol close or below above value and the rest-frame wavelengths sampled in these redshift ranges could be imprinted by galactic emission. The fits were performed adopting the same wavelength interval within the same redshift bin to avoid mismatches between the adopted continuum windows. The fitting wavelengths are listed in Table 2 and examples of a low- and a high-redshift fits are shown in the left and right panels of Fig. 2, respectively.

Fig. 2. Examples of a low-redshift bin (left) and a high-redshift bin (right) fit. The continuum power law is marked in blue, Fe II templates in grey, emission lines in different colours (orange, purple, green, cyan). The total model is depicted in red.

The uncertainties on the relevant parameters were estimated by means of mock samples. For each composite we produced 100 mock samples by adding a random value to the actual flux, extracted from a Gaussian distribution whose amplitude was set by the uncertainty value in that spectral channel. We then fitted the mock samples and evaluated the 3σ clipped distribution of the relevant parameters. The final uncertainty on the reported parameters was estimated as the standard deviation of such a distribution.

4. Results

4.1. The full sample and literature composite spectra

As a reference for the whole sample we produced a stacked spectrum collecting all the spectra of the L20 sample matching our quality criteria. The result is shown in Fig. 3. The procedure adopted to produce our template is the same as that described in Sect. 3.2. The normalisation wavelength $\overline{\lambda }$ was set to 2200 Å, as this is the only relatively emission-line-free continuum window available for most of the objects in our sample. To place our sample in a broader context, we compared our new template with other literature average spectra, namely Francis et al. (1991), Brotherton et al. (2001), VB01, Selsing et al. (2016), and Harris et al. (2016).

Fig. 3. Spectral composite of the L20 sample. Top panel: composite spectrum of our cosmological sample described in this work. The Lyα emission line is cut for visualisation purposes. The inset zooms in on the Hβ–[O III] region. The main emission lines are labelled accordingly. Middle panel: relative residuals with respect to the VB01 template. Bottom panel: number of objects contributing to each spectral channel. The composite spectrum of the L20 sample is available in its entirety in a machine-readable form at the CDS.

The composite spectrum by Francis et al. (1991) was produced using 718 objects from the Large Bright Quasar Survey (LBQS) with absolute magnitude in the B-band M_B ≤ −20.5 and redshifts in the range 0.05  <  z  <  3.36; this sample is strongly biased towards luminous objects at high redshift. The FIRST Bright Quasar Survey (FBQS, Brotherton et al. 2001) took advantage of FIRST radio observations of known point-like sources with colour indices³O − E < 2.0 and E < 17.8 from the Palomar Observatory Sky Survey I (POSS-I). This match yielded 636 sources, here we show only the composite spectrum made of the radio-quiet sources as a comparison.

The first SDSS quasar template was produced by VB01, who selected 2204 sources between 0.004 ≤ z ≤ 4.789 based on the SDSS colour difference (Richards et al. 2001), a non-stellar continuum, and at least one broad (FWHM > 500 km s⁻¹) permitted emission line. The sample adopted by Harris et al. (2016) was instead based on the Baryon Oscillation Spectroscopic Survey (BOSS; Dawson et al. 2012), as part of SDSS-III (Eisenstein et al. 2011), aimed at characterizing the effects of the baryon acoustic oscillations (BAO) on the distribution of luminous red galaxies and studying the Lyα forest in high-redshift quasars. These purposes resulted in a sample of 102,150 objects between 2.1 ≤ z ≤ 3.5 (see Harris et al. 2016 for more details about the sample selection). Given the redshift range covered by such sources, the comparison with this composite spectrum was possible only in the UV up to ∼3000 Å. Lastly, Selsing et al. (2016) observed with XSHOOTER/VLT a sample of seven bright (r ≲ 17) quasars between 1.0 < z < 2.1, for which galactic contribution should be regarded as negligible.

The comparison between our template and the literature ones is shown in Fig. 4. Note that, the y-axis represents λf_λ, rather than the usual flux density. The agreement is generally good, especially in the region between 1100–4000 Å. At shorter wavelengths, the intergalactic medium (IGM) can affect the composites differently, according to the distribution in redshift of the individual spectra upon which the stacks are based (e.g. Madau 1995; Faucher-Giguère et al. 2008; Prochaska et al. 2009). On the red side, instead, some intrinsic differences start to become visible. Our composite is almost indistinguishable from the VB01 composite, but these two seem flatter than the others. There are several explanations for this slight difference in the long-wavelength end of the composites. First, the colour selections performed to select the quasar candidates is not homogeneous among the different samples, thus the SDSS objects that constitute the VB01 sample as well as ours could present slightly redder indices in the optical range below z = 0.8. Secondly, it is possible that the host galaxy, whose contribution is negligible in the UV (see Sect. 4.5) plays a minor role in flattening the optical side of the spectra. VB01 estimated the galaxy contribution to be roughly 7–15% at the Ca IIλ3934 wavelengths, and a similar fraction can be estimated in our global composite. On the other hand, the other surveys and the Selsing et al. (2016) sample are made, on average, of more luminous objects for the same redshift (see inset in Fig. 4). It is therefore likely that the galactic contribution in those templates is even lower than in ours.

Fig. 4. Comparison between our quasar template and others from the literature. All the templates are scaled by their emission at 2200 Å. The agreement on the UV side, apart from minor differences in the Fe II pseudo-continuum, specifically the 2250–2650 Å and the 3090–3400 Å bumps, is apparent. The Lyα emission line is cut for visualisation purposes. The inset shows the distribution of the absolute R magnitude against the redshift for the different samples reported.

Throughout this work, we will refer to the full sample template as a reference to compare to the composite spectra of the individual bins in the M_BH – L_bol – z space, unless differently stated.

4.2. Non-evolution with redshift

A key observational feature bridging local and high-redshift quasars is the apparent similarity of their optical–UV continuum up to primordial times (e.g. Bañados et al. 2018; Yang et al. 2020; Wang et al. 2021). However, if the quasar spectra showed even a subtle evolution with redshift, the accretion physics underlying the X-ray to UV relation could also be evolving, making quasars unsuitable candidates as standard candles. Given the wide range in terms of redshift covered by our sample, we expect to be capable of detecting any such possible effect. This would appear as a systematic difference in the slope of the continuum produced by the accretion process at different cosmic times.

As a very basic step to check for any evolution with redshift of the spectral properties within our sample, in Fig. 5 we compare the average spectra in the redshift bins with the global composite (described in Sect. 4.1), obtained by stacking the spectra of all the sources in the sample. We find a very good agreement with the global composite at all redshifts, without any clear evidence of evolution, apart from the known trends of the emission lines, whose FWHM varies with M_BH (i.e. the virial relation; e.g. Peterson et al. 2004) and whose EW anti-correlates with L_bol (i.e. the Baldwin effect; Baldwin 1977).

Fig. 5. Composite spectra divided in redshift bins superimposed to the full sample composite. The dashed line denotes the number of objects contributing to the stack in each spectral channel. The composite spectra in each bin are scaled to match the template respectively at 4200 Å for the low redshift bin, 3000 Å for the z = [0.77 − 1.12), z = [1.12 − 1.52), z = [1.52 − 2.03) bins, and 2200 Å for the highest redshift bin. In the bottom-right panel we show the slopes of the continuum evaluated in different pairs of continuum windows, which are generally consistent.

The visual agreement is confirmed by a more quantitative assessment presented in the bottom-right panel in Fig. 5. There we show the slope of the continuum evaluated in continuum intervals sampled at different redshifts by the spectral composites. The slopes are generally consistent within 1σ, and no systematic trends are observed. The optimal agreement between the composite spectrum in every redshift bin and the global composite, as well as among the different redshift composites, demonstrates the non-evolution of the average continuum properties of the sample across cosmic time.

In addition, as discussed in Section 4.1, the global sample composite exhibits a remarkable agreement with the average SDSS spectrum described in VB01, meaning that the stacks in the different redshift bins are also similar to the spectrum of a typical blue quasar. Since even this direct comparison reveals that the average continuum properties of our sample do not suggest any redshift evolution, we now investigate and characterise in more detail how our sample behaves across the log M_BH – log L_bol parameter space.

4.3. Analysis of the evolution of the continuum slope of the stacks

The fits of the stacked spectra returned the slope of the continuum in each bin of the parameter space. We thus explored whether any possible evolution of the continuum slope exists throughout our sample. With the term slope (α_λ), we denote the slope of the continuum power law underlying the emission lines roughly between 1300 and 5500 Å. This wavelength range is sampled differently in the different redshift bins: for instance in the lowest redshift bin only the interval 2200–5500 Å is available, while in the highest redshift bin the composite spectra were fitted over the interval 1300–3200 Å. We also underline that, throughout this work, we never assign a physical meaning to the slope of the continuum, adopting it exclusively as a comparative parameter. The physical meaning of such a parameter would require the assumption of a physical continuum model, generally an accretion disc, but see also e.g. Antonucci et al. (1989) for some limitations, whose implications go beyond the scope of this work.

In Fig. 6, we show the values of the slope in the redshift bins as a function of log M_BH, log L_bol, and log λ_Edd. To avoid dust reddening or galaxy contamination, a selection in terms of photometric colours has been performed when constructing the cosmological sample (see L20, Section 5.1). In brief, for each object, the slope of the power law in the log ν – log(νf_ν) space was computed in the rest-frame ranges 3000–10 000 Å (Γ₁) and 1450–3000 Å (Γ₂). Objects whose photometric colours were compatible with those of typical blue quasars, Γ₁ = 0.82 and Γ₂ = 0.40 (Richards et al. 2006), only tolerating mild reddening, i.e. E(B − V)  <  0.1, were selected. Converting the photometric indices Γ in terms of spectral slopes as α_λ = −Γ − 1, the expected slopes in the 3000–10 000 Å and the 1450–3000 Å ranges, respectively covered by the low and high redshift bins, are α_λ, 1 = −1.82 and α_λ, 2 = −1.40. Although this process operates a sort of pre-selection on the possible continuum slopes, we note that only roughly 20% of the starting sample is excluded according to this colour selection.

Fig. 6. Values of the slope of the continuum produced by our fits. Bins not meeting the representativeness criterion are marked with a circle. The slope of the VB01 template (dot-dashed line) is shown as an external reference. Different pairs of line-free windows were adopted to define the continuum in each bin.

In Fig. 6 we show the slopes of the continuum evaluated in the bins of the parameter space. We mark with a hollow circle the values of the bins with a low N_obj which would be excluded as a result of the representativeness criterion described in Sect. 3.1. The values cluster around α_λ = −1.50 with a standard deviation of σ_α = 0.14, showing little or no evolution from low to high redshift. Considering also the excluded bins, the results are not significantly altered as the mean value negligibly shifts to α_λ = −1.51 and the standard deviation to σ_α = 0.17. The slope average value is close to those expected from the Richards et al. (2006) typical photometric index and from the templates of blue quasars described in Sect. 4.1, whose continuum slopes, evaluated using our fitting routine, are reported in the inset panel of Fig. 7. We recall that the slopes are evaluated using different continuum windows at different redshifts. For a fixed value of M_BH (or equivalently L_bol), the distribution in the observed slopes is therefore partly due to the spread along the other axis of the parameter space, which is anyway modest, and partly driven by the different continuum windows.

Fig. 7. De-reddened stacks, according to the same colour code as in Fig. 13, scaled by their mean emission between 1770 Å and 1810 Å. For comparison, the VB01 template is also shown in black as an external reference, rather than the usual L20 template. The rising continuum slopes with increasing redshift (and extinction) are clearly observed. Inset: UV continuum slopes of the literature samples described in Section 4.1 and plotted in Fig. 4 (squares), against the slopes of the de-reddened spectra (dashed lines).

We assessed the degree of correlation between the slope and the axes of the parameter space by means of the Spearman rank correlation index S (correlation increases in modulus from 0 to 1). We do not find any significant trend with any of the accretion parameters, (i.e. |S|≤0.3 and associated p − values  > 0.05). The inclusion of the bins excluded as a result of the representativeness criterion, does not change significantly the correlation index, as we show in Appendix E.

In addition to exploring the evolution of the continuum within our parameter space, we also explored how the number of sources in each bin affects the continuum properties with respect to the global template. In particular, in Fig. 8 we present the relative difference between the continuum slope for each bin (α_bin) and the full sample composite (α_FS) evaluated in the same wavelength range Δα = 1 − α_bin/α_FS against the number of objects in the bin. Regardless of the redshift, the difference decreases in modulus for increasingly populated bins, as if the samples were drawn from a single population. In this picture, the spread in the observed slopes results from sample variance, rather than from an actual evolution with the accretion parameters explored.

Fig. 8. Values of Δα (Δα = 1 − αbin/αFS) in the bins of our parameter space as a function of the objects in each bin. Bins not meeting the representativeness criterion are marked with a circle. The dashed lines define the $\sim 1/\sqrt{N_{\mathrm{obj}}}$ envelope characterising the standard error on the average slope assuming a standard deviation σα = 1.

4.4. Comparison of the 2500 Å monochromatic fluxes

We produced the normalised monochromatic 2500 Å fluxes (f_{2500 Å}) by interpolating the continuum resulting from the best fit slope and assuming a unitary monochromatic flux at 3000 Å. In Fig. 9 we present f_{2500 Å} against log M_BH, log L_bol, and log λ_Edd in order to study any possible evolution. In principle, a systematic change of the f_{2500 Å} with the accretion parameters would require a correction to ‘standardise’ the UV proxy. However, there is evidence against this possibility.

Fig. 9. Values of the normalised 2500 Å fluxes along the axes of the parameter space, assuming a common normalisation of the continuum at 3000 Å in all the bins. Bins not meeting the representativeness criterion are marked with a circle.

In Signorini et al. (2023) we showed that the monochromatic 2500 Å flux is a good UV proxy for the L_X − L_UV relation, although other wavelengths in the range 1350–5100 Å can perform equally well (see also Jin et al. 2024). Therefore, even assuming different slopes, the UV term of the L_X − L_UV relation would hardly be affected.

Here we aimed at quantifying the effect on the monochromatic 2500 Å fluxes of the differences in the spectroscopic slopes within our sample. The slopes are comprised between roughly [−2.0, −1.1], with the extremal values found in scarcely populated bins (respectively 9 and 6 objects), which would be excluded as a result of the procedure described in Sect. 3.1. Since our goal is to corroborate our sample for its cosmological application, we can estimate how much the spread in the observed slopes affects, for instance, the distance modulus (DM) estimates. Here we show that, even assuming a systematic evolution in the continuum slope with one of the explored parameters (an occurrence which is not supported by our analysis), the variation induced in the f_{2500 Å} would propagate into a negligible difference in the DM for our scopes. We can easily check this by expressing DM in terms of f_{2500 Å} and adopting the expression (see Risaliti & Lusso 2015, 2019; Lusso et al. 2020)

$$\begin{array}{cc}\hfill \mathrm{DM}& =5[\frac{log{f}_{\mathrm{X}}-\beta -\gamma (log{f}_{2500\mathrm{\AA }}+27.5)}{2(\gamma -1)}-\frac{1}{2}log(4\pi )+28.5]\hfill \\ \hfill & -5log(10\mathrm{pc}),\hfill \end{array}$$(1)

where f_X and f_2500 Å are, respectively, the X-ray and UV flux densities in units of erg s⁻¹ cm⁻² Hz⁻¹. The parameter log f_2500 Å is normalised to the value of 27.5, whereas the logarithm of the luminosity distance is normalised to 28.5. The γ = 0.600 ± 0.015 term is the slope of the L_X − L_UV relation reported in L20, and β is the intercept. Here we can express the 2500 Å monochromatic flux in a power-law form pivoting at 5100 Å (i.e. the farthest among the explored characteristic wavelengths). By differentiating Eq. (1) with respect to the slope of the continuum α_λ, we get

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathrm{DM}=|-\frac{5\gamma }{2(\gamma -1)}log\left(\frac{2500}{5100}\right)|\mathrm{\Delta }{\alpha }_{\lambda }.\end{array}$$(2)

We can then substitute the standard deviation of the observed slopes (σ_α = 0.17) as a measure of the spread of the Δα_λ term. The DM variation ΔDM is ∼0.2 in modulus, is not sufficient to explain the differences observed with the ΛCDM at high redshifts (see 13). If, as we believe, 2500 Å is a better approximation to the UV pivot point than 5100 Å, the average ΔDM is expected to be even smaller.

4.5. Host galaxy contamination

Although the broadband properties of the objects included in the sample are designed to minimise the additional contribution of the host galaxy in the quasar spectrum, here we perform an a posteriori test based on the spectral data to estimate the average degree of host galaxy emission in our sample. In Rakshit et al. (2020) a host galaxy–quasar decomposition was performed using a principal component analysis (PCA, Yip et al. 2004) in two windows close to 4200 Å and 5100 Å, obtaining reliable results for ∼12 700 objects below redshift z ≤ 0.8. We employed these measurements to assess the level of residual (i.e. not filtered out by our criteria) host galaxy contamination within the quasars of our sample. With this aim, we built a purposefully designed sample of sources below z = 0.8 to ensure the coverage of both the host fraction indices, from the latest WS22 catalogue by adopting the following criteria.

We considered only the objects for which the host galaxy fraction indices at 4200 Å (HOST_FR_4200) and 5100 Å (HOST_FR_5100) were both evaluated, with values ranging continuously from 0 (quasar-dominated emission) to 1 (galaxy-dominated emission). If any of the two indices was flagged as not measured (i.e. = − 999), the object was discarded.

We defined the host galaxy fraction index (HGF) and created 10 bins of width ΔHGF = 0.1. A source was a given HGF index when both the 4200 Å and the 5100 Å fractions fell within the same HGF bin in order to define an unambiguous index. The sample so constructed has 5748 objects.

We then stacked the spectra of the objects in each bin according to the same procedure described in Sect. 3.2. In Fig. 10 we show the stacked spectra as a function of their HGF index. To perform a fair comparison, we also produced a stack of all the 529 sources in our sample with z ≤ 0.8 (shown in black). On average, the quasars in the low redshift tail of our sample retain minimal levels (≲10%) of host galaxy contamination below 5400 Å, in good agreement with VB01 who report a galaxy flux fraction between 7% and 15% on the basis of the Ca IIλ3934 absorption feature. At redder wavelengths, the host contribution is expected to become increasingly important. In addition, we also note that the difference is terms of continuum slope and line strength is more evident in the UV region of the spectra, where our template exhibits the best match with the other samples of blue quasars. This test confirms that our UV proxy for the L_X − L_UV relation is not systematically affected by host galaxy contamination.

Fig. 10. Comparison between the stacks in bins of HGF and our reference spectrum (black) for the z ≤ 0.8 interval, as colour-coded in the legend. All the spectra are scaled to their 4200 Å flux. Increasingly host-dominated composites exhibit flatter UV continua.

4.6. Intrinsic extinction

Our photometric selection is designed to exclude dust reddened quasars, by retaining only sources near the locus in the Γ₁ – Γ₂ plane where blue unobscured quasars are expected to reside (Richards et al. 2006). We again checked, a posteriori, that the spectra composing our sample only accept minimal levels of extinction, and we bring three arguments in favour of this hypothesis.

First, as a general consideration, we find a very good agreement between our stacks and the VB01 SDSS template, which is made of sources whose colour-selection had been designed to select blue quasars (Richards et al. 2001). The agreement is excellent – especially on the bluer UV side – also with the other average spectra shown in Sect. 4.1. This implies that, if the quasars used to build such templates are not redder than the average, the same should hold for our sample. We also recall that, by selection, the VB01 template was built using objects with at least one broad emission line, thus the presence of Type 2 obscured objects should be limited.

As a second test, to get a more quantitative assessment of the effect on the spectral slope of increasing reddening, we reddened the VB01 template under the reasonable hypothesis that this is representative of the intrinsic spectrum of a typical blue quasar. We reddened the template with increasing E(B − V) from 0 to 0.1 with steps of 0.01, using the extinction curve described in Gordon et al. (2016) for an average Small Magellanic Cloud (SMC) extinction curve, which according to several authors is a good approximation of the AGN extinction curve (see e.g. Brotherton et al. 2001; Hopkins et al. 2004), and then checked the agreement with our full sample stack. A more detailed discussion on how a different extinction curve can affect this procedure is presented in Appendix C, where other prescriptions are explored. In brief, we found that the variation of the continuum slope for a given colour excess derived adopting the Gordon et al. (2016), used to account for intrinsic reddening throughout this work, is the closest to the average (see Fig. C.1). The result of the reddening procedure is shown in Fig. 11. The best match is obtained for minimal levels of intrinsic reddening, i.e. E(B − V) ≲ 0.01, an estimate similar to that provided by Reichard et al. (2003) on the basis of spectral data from the SDSS early data release. To quantify this effect, we performed a spectral fit of the reddened templates in the region between 1300 and 4100 Å, and compared the slope of the continuum with that obtained by the fit of our stack in the same region. This test confirmed the finding emerging from a crude comparison of the overall spectral shape, that is a continuum slope consistent with a low degree of intrinsic extinction. The slopes resulting from the spectral fit of the reddened templates are shown in the inset of Fig. 11. We note that, choosing a different template than VB01 would produce analogous results, given the similarity of all the explored templates in the near UV.

Fig. 11. Extinction of the VB01 template according to the colour-coded values. All the spectra are scaled by their 4200-Å flux, where extinction becomes less relevant. The black spectrum is our full sample stack and is in very good agreement with a low degree of extinction. The inset shows the best-fit slope in the range 1300–4100 Å of the VB01 template reddened as a function of the E(B − V) value. The black dot-dashed line represents the slope of the continuum fitted to our stack. Minimal levels of extinction are allowed to match the reddened template with our composite. We note that the visual match between the L20 and the VB01 reddened template depends slightly on the scaling wavelength. The slope derived from the fits, however, does not depend on such a wavelength.

Finally, we estimated the Balmer decrement of the low-redshift sources with available broad Hα and Hβ emission lines to assess possible levels of intrinsic extinction of the BLR. Dong et al. (2007) computed the distribution of the broad-line Balmer decrement Hα_br/Hβ_br ≡ f(Hα_br)/f(Hβ_br) for a sample of 446 low-redshift (z ≲ 0.35) blue AGN, finding that the distribution can be well approximated by a log-normal, with a mean value of ⟨log(Hα_br/Hβ_br)⟩ = 0.486 and a standard deviation of 0.03. A value of log(3.1) = 0.49 is generally adopted for the Balmer decrement in the narrow line region in AGN (Halpern & Steiner 1983; Gaskell & Ferland 1984), while it is found to vary widely according to different BLR conditions (e.g. Dong et al. 2007 and references therein).

In this context, large Balmer decrements are interpreted as due to local reddening, as also suggested by large ratios between infrared emission and broad lines (Dong et al. 2005) and significant X-ray absorption (Wang et al. 2009). In Fig. 12, we show the distribution of log(Hα_br/Hβ_br) ratios assuming the values reported in the WS22 catalogue for the 167 objects in our sample with reliable Hα and Hβ measurements (see Sect. 4 of WS22 for the detailed criteria). The mean value of the distribution is ⟨log(Hα_br/Hβ_br)⟩ = 0.48, and the standard deviation is 0.14. Under the assumption that the sample of blue quasars described in Dong et al. (2007) retains, on average, minimal levels of BLR extinction (see their Sect. 4.1), we find that the intrinsic extinction is negligible, at least for the low redshift sample.

Fig. 12. Distribution of log(Hαbr/Hβbr) for our sources with available data. The mean values from other literature samples (Osterbrock 1977, O77, Neugebauer et al. 1979, N79, Goodrich 1990, G90, Greene & Ho 2005, GH05, Zhou et al. 2006, Z06, Dong et al. 2007, D07, Mejía-Restrepo et al. 2022, MR22, Ma et al. 2023, M23) are also shown for reference according to the colour code. The mean of our values is shown in black. The BASS sample described in Mejía-Restrepo et al. (2022) has been split according to the column density NH inferred from the X-ray analysis.

Finally, the combination of the three tests proposed here strongly suggests that intrinsic reddening should be minimal accross the entire redshift range covered by our sample.

4.7. Is the ΛCDM tension driven by reddening?

One of the main consequences of reddening would be the underestimate of the 2500 Å monochromatic fluxes. A relatively small amount of extinction, such as E(B − V) ≃ 0.1, would diminish the 5100-Å flux by ∼25%, but the extincted radiation would be roughly double (depending on the assumed extinction curve) at 2500 Å, as a result of the differential cross Section of dust. The net effect of extinction is thus the underestimation of the DM, enhancing the tension between our cosmographic best-fit model and the high-redshift (z ≳ 2–3) extrapolation of the ΛCDM model (Risaliti & Lusso 2019; Bargiacchi et al. 2023). Here we focus on the possibility that the observed tension with the ΛCDM model is caused by local dust extinction affecting the UV measurements. To this end, we performed the following exercise, neglecting any possible effects on the X-ray fluxes. We assumed that the difference between the DM predicted by ΛCDM and that based on the L_X − L_UV relation is entirely driven by local (i.e. at the source redshift) extinction on f_{2500 Å}. Therefore, assuming a standard extinction curve, we can evaluate the average E(B − V) by which f_{2500 Å} should be de-extincted in order to match the predicted DM. With this quantity in hand, we can de-redden the UV spectra and check how much their intrinsic shapes would differ from those observed. In the following, the term “de-reddened” will be used to indicate the extinction-corrected spectral fluxes that produce a DM complying with the ΛCDM model.

Fig. 13 shows the samples used for this purpose. We gathered spectral data of objects at z ≥ 1, and divided them into 8 logarithmically equi-spaced redshift bins. For each of them, we evaluated the average redshift and DM_obs values. The DM can be written in terms of f_2500 Å by using Eq. (1) as:

$$\begin{array}{c}\hfill \mathrm{DM}=-\frac{5}{2}\frac{\gamma }{\gamma -1}log{f}_{2500\mathrm{\AA }}+C,\end{array}$$(3)

Fig. 13. Effect of the intrinsic extinction on the Hubble–Lemaître diagram. Top panel: Hubble–Lemaître diagram of the the sources used to investigate the possible effects of local reddening. The redshift bins used for this purpose are colour-coded in the legend together with the number of objects per bin and the corresponding E(B − V). Coloured stars represent the average values in each redshift bin, while black stars mark the ΛCDM model expectations, with red segments joining these two values. Increasingly high values of E(B − V) are needed in order to produce the observed tension with the concordance cosmological model. Bottom panel: difference between the predicted and observed distance modulus (ΔDM). Values are not extrapolated below z = 1.0 and above z = 4.8.

where C includes the X-ray flux contribution, the normalisation β, and the other constants. If we assume that the measured DM values are affected by reddening, we can evaluate the amount of extinction needed to place the intrinsic DM exactly onto the ΛCDM expectations and de-redden f_2500 Å accordingly. Such a correction can be expressed as (see Appendix F for the explicit derivation):

$$\begin{array}{c}\hfill E(B-V)=\mathrm{\Delta }\mathrm{DM}\frac{1}{R_V}\frac{A_V}{A_{2500\mathrm{\AA }}}\frac{1-\gamma }{\gamma },\end{array}$$(4)

where ΔDM is the difference between the average distance modulus in a given redshift bin according to the ΛCDM model and that observed based on the L_X − L_UV relation. Assuming a selective to total extinction ratio R_V = 3.1 and the SMC extinction curve from Gordon et al. (2016) used above, we retrieve the values listed in the legend of Fig. 13.

Next, we produced the composite spectrum in each bin, de-reddened it by applying the E(B − V) values found in the previous step, and performed the spectral fits to each of them to evaluate the slope of the continuum. We show the so derived spectra and the slopes of the continuum obtained by means of spectral fits in Fig. 7. The slopes characterising the continua of the de-reddened spectra are not compatible with the expectations based on literature samples, as shown in the inset of Fig. 7. Indeed, moving towards higher redshift, under the hypothesis of local reddening, we would select intrinsically bluer objects with increasing amounts of extinction, whose UV slopes are anomalously steeper with respect to those typically found in blue quasars, while still delivering similar spectral shapes at all redshifts.

With the aim to explore the consequences of de-reddening the spectra by such amounts in terms of quasar physics, we built the semi-analytic (i.e. using our bin-averaged data and adopting analytic approximations from the literature in the other bands) intrinsic quasar SED in each redshift bin and computed the hard X-ray (2–10 keV) to bolometric corrections (k_bol, X) to be compared with other literature samples. Anomalously high k_bol, X would further testify the selection of a very peculiar class of objects at high redshift, seemingly at odds with the typical spectral properties upon which the sample is built, and which we have verified so far.

In order to produce the average SEDs we built a semi-analytical parametrisation, similar to that discussed in Marconi et al. 2004 (M04). For wavelengths λ > 1 μm we adopted a power-law spectrum L_ν ∝ ν^α_ν, with slope α_ν,IR = 2 as in the Rayleigh-Jeans tail of the disc. Between 1 μm > λ > 1200 Å we adopted the spectral index α_ν,UV found in the fits of the de-reddened spectra, under the assumption that at high redshift and high luminosity the continuum can be approximated as a single power law. In the range 1200 Å > λ > 600 Å we chose the average far UV spectral slope α_ν,FUV = −1.70 described in Lusso et al. (2015), who took into account a statistical correction on the neutral hydrogen column density of the inter-galactic medium. On the X-ray side, we used the pexrav model (Magdziarz & Zdziarski 1995) in the XSPEC software (Arnaud 1996) to simulate a power-law spectrum with a photon index Γ equal to the average photometric photon index derived from XMM-Newton photometry (see Sect. 4.8) with a cutoff at 500 keV, and a reflection component with solid angle of 2π, inclination angle of cos i = 0.5, and solar metal abundances⁴. Each X-ray spectrum was then scaled to the average 2 keV luminosity in each bin.

Lastly, we connected the 600 Å and the 0.5 keV luminosities with a power law. We defined 8 template SEDs according to the 8 values of L_{2500 Å,int} and α_ν,UV, and used them to compute the hard X-ray (2–10 keV) bolometric corrections k_bol, X. We estimated the uncertainty on k_bol, X by building 100 mock SEDs in each bin assuming Gaussian distributions for the adopted parameters with standard deviations σ(α_ν,IR) = 0.1, σ(α_ν,UV) equal to the uncertainty on the continuum slope evaluated from the fit of each de-reddened spectrum, σ(L_{2500 Å,int}) equal to the standard deviation on the mean of the L_{2500 Å,int} values in each redshift bin, σ(α_ν,FUV) = 0.6, σ(L_2 keV) and σ(Γ) equal to the standard deviation on the mean of the L_2 keV and Γ distributions in each redshift bin.

The average SEDs resulting from this procedure are shown in Fig. 14, together with other literature samples as considered in Duras et al. 2020 (D20). Our k_bol, X values from the de-reddened SEDs clearly detach from the high-luminosity extrapolation of the M04 best-fit relation. At the same time, they are all located on the upper envelope of the D20 best fit for luminous type-1 sources. This implies that, if such sources existed, we would be selecting an exotic population made of the most luminous, yet dust-extincted quasars. These sources would exhibit the highest known bolometric conversion factors, while still showing the UV-optical colours and spectral properties of typical blue quasars. Moreover, if our SED templates are consistent by construction with M04, which then represents a proper comparison, the same does not hold for the empirical SEDs in D20. Indeed, we did not include any physically motivated emission in the spectral region between the far UV and the soft X-rays, as we only linked the L_6000 Å and L_0.5 keV with a straight power law, while in the latter work the authors could also rely on far-UV photometry from GALEX. In order to quantify the possible inconsistency of our average estimates with those from the other literature samples, we performed a basic fit of the k_bol, X – L_bol relation in our de-reddened sources. We adopted a simple power-law functional form, which is a good approximation for the high luminosity (L_bol ≥ 10⁴⁶ erg s⁻¹) end of the D20 distribution:

$$\begin{array}{c}\hfill log(k_{\mathrm{bol},\mathrm{X}})=a+b[log\left(\frac{L_{\mathrm{bol}}}{{10}^{47.5}{\mathrm{erg}\mathrm{s}}^{-1}}\right)].\end{array}$$(5)

Fig. 14. Average hard X-ray bolometric corrections in each redshift bin. As a term of comparison, we report the kbol, X for the semi-analytical SED templates from M04 (blue squares) and the empirical ones from D20 (pentagons). The blue dashed line represents the best-fit reported in M04. Black and red continuous lines show, respectively, the best fit of our sample and the D20 objects in the high Lbol (≥1046 erg s−1) regime. The red dashed lines represent the 95% confidence interval of the best-fit regression line for the selected high-luminosity objects from D20 (red pentagons). Our values of kbol, X from the high-redshift, de-reddened SEDs are inconsistent with the extrapolation of the lower luminosity trend from M04, and steadily offset with respect to our best fit of the D20 high-luminosity type-1 quasars. Inset: average SEDs in each redshift bin, according to the same colour code as in Fig. 13. The shaded area is given by the individual mock SEDs in each bin.

The best fit parameters and their respective uncertainties are listed in Table 3. In Fig. 14 we show the k_bol, X – L_bol relation for our sample and the D20 compilation of type-1 objects. The bolometric corrections of our higher redshift bins are located outside the 95% confidence interval for the best fit of the high-luminosity tail of the D20 distribution (red pentagons). The offset a of the best fit relation for our sources is significantly higher (3.7σ) than for those from D20. We checked that this result is not biased by the choice of the luminosity threshold for the objects included in the fit, by performing again the fit increasing the threshold by steps of Δlog(L_bol/erg s⁻¹) = 0.5 from 45.5 to 47. In each case, the high-redshift (z ≳ 2.5) bolometric correction factors remain located above the 95% confidence interval. In conclusion, the bolometric correction estimates derived from our average de-reddened SEDs are significantly offset with respect to the expectations for high-luminosity type-1 quasars, hinting at a seemingly exotic – and yet unknown – population.

As a result of all the tests performed in this Section, we conclude that a dust extinction scenario as an explanation for the tension with the ΛCDM model emerging in the quasar Hubble–Lemaître diagram seems disfavoured for several reasons. Firstly, significant amounts of extinction at high redshift, in combination with steeper than the average UV spectra, are required in order to reconcile the tension with the ΛCDM model, unless we assume a cosmological evolution from high-redshift extremely blue sources to low-redshift typical quasars. Secondly, the reddening should increase with redshift, in opposition to the current paradigm of increasing dust with time (see, e.g. Aoyama et al. 2018, and references therein). In addition, the dust content in the high-redshift quasars should increase by the exact amount needed to compensate for the steepening of the intrinsic spectrum, thus producing the typical spectrum of blue quasars at all redshifts. Lastly, the hard X-ray bolometric correction factors derived from the average de-reddened SEDs are in tension with the high luminosity extrapolation from known sources.

All the above considerations would imply the selection of exotic objects at high redshift whose emission properties are intrinsically different from those at lower redshift and luminosity, while their apparent broadband properties are consistent with those of typical blue quasars.

4.8. The X-ray stacks

The availability of X-ray data, which underlies the creation of the L20 sample, allowed us to address the question on the evolution of the shape of X-ray emission with respect to the accretion parameters and redshift. However, at this stage, a full stacking of the X-ray spectra, taking into account the largely different redshift intervals sampled, as well as the different response functions is still unfeasible. For this reason, we adopted a slightly different approach and built ‘photometric X-ray stacks’. The intrinsic X-ray spectrum of an unobscured quasar above 0.2 keV can be roughly described as a power law with spectral index α = 1 − Γ, where Γ is the photon index. For our sample the best estimate for the Γ parameter is the photometric photon index (Γ_ph) derived following the procedure described in Section 4 of L20 (where also a more detailed description of photometric spectral indices and of their limitations can be found). Although the photometric Γ values are somewhat less reliable than those spectroscopically derived (Γ_sp), the lack of a systematic offset in the ΔΓ distribution (defined as ΔΓ = Γ_sp − Γ_ph) and its non-evolution with redshift allow us to use these estimates in a statistical sense. For each source we built a synthetic spectrum between 0.2 and 10 keV using Γ_ph and the 2 keV flux reported in L20. Then, the procedure to build the X-ray stacks is analogous to that explained in Sect. 3.2.

We note, however, that the X-ray stacks are less informative than those produced in the UV-optical for at least two reasons: 1) they are not computed using the actual X-ray spectra, but rather photometric proxies, which are not sensitive to the peculiar features of each spectrum other than the steepness of the spectrum itself; 2) the quality cuts performed in the X-ray selection (especially that concering the lower limit on the photon index, i.e. Γ_ph − δΓ_ph > 1.7, with δΓ_ph being the uncertainty) are more conservative, given the generally lower quality of the X-ray data, and thus have a much greater impact on obliterating the spectral diversity than the optical cuts. Nonetheless, the analysis of the slope of the stacked X-ray spectra provides valuable information about the goodness of our selection as well as the possible evolution of the X-ray properties within our parameter space. In Fig. 15 we recover the well-known correlation between Γ – log λ_Edd (Shemmer et al. 2008; Risaliti et al. 2009; Jin et al. 2012; Brightman et al. 2013; Trakhtenbrot et al. 2017), with a correlation index S = 0.38. Interestingly, the correlation index is marginally larger for the Γ – log M_BH (S = −0.41), but this could be due to the way the photometric data were grouped to produce the stacks, rather than to an intrinsically tighter correlation.

Fig. 15. Values of Γ for our X-ray stacks across the parameter space. Bins not meeting the representativeness criterion are marked with a circle.

In theory, the dependence of the 2 keV fluxes on M_BH or λ_Edd implies that adding one of these parameters, or both, to the L_X − L_UV space, the dispersion could be further narrowed down. For instance, the dependence on M_BH, proportional to the Mg II FWHM, has been explored in Signorini et al. (2023), where it has been shown that the inclusion of this parameter has a negligible effect in improving the intrinsic dispersion of the relation. At the same time, the zero-covariance between the X-ray flux estimated at the pivot point (i.e., the effective-area-weighted average energy of any given band of interest) and Γ (see Section 4 in L20) assures that even a systematic steepening of the photon index would not affect our estimates of the X-ray luminosity.

5. Discussion

The suitability of quasars as ‘standardized candles’ relies upon the consistency of the method (i.e. the absence of selection/analysis biases, or correction thereof) and the intrinsic non-evolution of quasar spectra with redshift. Although increasingly larger samples are desirable, our group has been lately devoting a substantial effort to investigate any possible issues affecting the L20 sample with SDSS spectral coverage. In this context, this work has been chiefly focused on the analysis of the average spectral properties of 2316 of the sources described in L20. The main results obtained through the various analyses performed in the previous sections are the following: 1) The overall shape of the composite spectra (apart from emission-line properties) does not show any clear evolution across the explored parameter space. 2) Host galaxy contamination is not a major issue in our sources below z ≤ 0.8, and we therefore expect it to be completely negligible for the bulk of our sample, at higher redshfit and luminosity. 3) Our sample is compatible, on average, with a low degree of intrinsic extinction, E(B − V) ≲ 0.01. 4). The deviation of the cosmographic best fit of the quasar–SN Ia Hubble–Lemaître diagram from the ΛCDM model reported by our group is unlikely to be caused by dust reddening, based on the spectral properties of our sample. These findings have far-reaching consequences on the physical properties of the sources described here, as well as on the solidity of the criteria adopted to select our cosmological sample.

The similarity of quasar spectra has led through the years to the production of several composite spectra (e.g. Francis et al. 1991, VB01, Zheng et al. 1997; Brotherton et al. 2001; Harris et al. 2016), which have found large employment as templates for the typical quasar emission, in terms of line and continuum properties. In our case, building high signal-to-noise composite spectra allowed us to explore the spectral diversity within the parameter space underlying SMBH accretion. We do not report any clear evidence for a redshift-dependent evolution of the average continua in our sample, as the agreement with either the global L20 or the independent VB01 composite spectra does not deteriorate with redshift. This appears to hold up to very high redshift, where luminous quasars have been observed to show spectral properties very akin to those seen at lower redshift (e.g. Fig. 1 in Mortlock et al. 2011). Under the paradigm that such emission is driven by SMBH accretion, these results would imply that the emission mechanism underlying the observed spectra does not vary significantly in the blue quasars spanning the last 12 Gyr. Although we do observe the expected correlations between line equivalent width and continuum luminosity (i.e. the Baldwin effect, Baldwin 1977) and between line width and BH mass (virial relations; e.g. Peterson et al. 2004), the non-evolution of the continuum shape of our composite with any of the quantities in our parameter space is remarkable (see Fig. G.1).

Contamination from the host galaxy light is a known issue when trying to isolate the AGN component (e.g. Baldwin et al. 1981; Kauffmann 2003). This effect is capable, in principle, of inducing the overestimate of the actual 2500 Å fluxes, by introducing a spurious UV contribution and eventually biasing the estimate of the cosmological parameters based on the quasar L_X − L_UV relation. For this reason, our selection criteria require steep photometric continua to single out sources where the contribution from the galactic stellar component is negligible with respect to the AGN emission (Elvis et al. 1994; Richards et al. 2006; Elvis et al. 2012; Krawczyk et al. 2013; Hao et al. 2013). As shown in Section 4.5, host galaxy contamination does not play a major role in our sample, as we generally observe steep continua, coupled with the lack of strong stellar absorption features and weak narrow emission lines from stellar photoionisation. These findings are particularly evident on the UV side of the spectrum, where our composites are dominated by higher redshift, luminous objects, in which the host contribution is expected to be negligible above L_bol ∼ 10⁴⁶ erg s⁻¹ (Shen et al. 2011). A detailed disentanglement of the host galaxy and quasar emission by means, for instance, of a principal component analysis (e.g. Yip et al. 2004) is beyond the scope of this work. Nonetheless, we can roughly estimate the average incidence of host galaxy light by comparing the absorption features between our full composite and VB01. The EW of the Ca IIλ3934 absorption in our composite (0.7 ± 0.4 Å ) is similar to that observed in the VB01 template (0.9 ± 0.1 Å ), yet a fair comparison is hampered by the differences in the resolution and the signal to noise between the two templates. This goes in the direction of our sample having an average host contribution comparable to that of VB01, deemed 7%–15% at such wavelengths, with an incidence decreasing bluewards.

Recently, Zajaček et al. (2023) have investigated the possibility that luminosity distance measurements based on the L_X − L_UV relation are biased by extinction on a sample of 58 quasars at z ≲ 1.7 with available reverberation mapping data (Khadka et al. 2023). In the latter work, the authors used Mg II time-delays to infer distances by means of the radius–luminosity relation (Mg IIR–L relation, with intrinsic dispersion σ ∼ 0.3 dex; Khadka et al. 2022). Assuming that the difference between the distances based on the Mg IIR–L and L_X − L_UV relations were attributable to local reddening, and adopting different cosmological models, they derive an intrinsic colour excess E(B − V) ≃ 0.002–0.003, an amount of extinction that would be intrinsic to all type-1 quasars. However, such an effect, even if present, would still not be enough to explain the reported high-redshift tension with the current cosmological model.

In this framework, the shape of the quasar continuum can be interpreted as ultimately resulting by the combination of two main parameters: the intrinsic spectral slope of the UV spectrum and the column density of gas and dust causing the extinction. The degeneracy between these two parameters could be capable, in principle, of producing similar spectra by combining increasingly larger amounts of dust in high redshift, luminous sources with intrinsically steeper spectra (e.g. Xie et al. 2015).

There are several reasons to believe instead that dust extinction, even if present in small amounts in our sample, is not systematically biasing our results. First, the photometric selection, besides excluding sources where the host galaxy significantly reddens the UV-optical SED, also singles out objects where the amount of dust extinction is expected to be negligible. Indeed, assuming that the photometric slopes of the average quasar SED described in Richards et al. (2006) are truly representative of blue quasars, by selection, the colour excess of our sources should not exceed E(B − V) ≃ 0.1. We also note that the incidence of significantly reddened objects with E(B − V) > 0.1, at least in early SDSS data (whose average properties are encapsulated in the VB01 template) was estimated to be modest, below < 2% according to Hopkins et al. (2004). In this work, we took a step further, searching a posteriori for spectral signatures of dust reddening, but we found none. Moreover, the distribution of the broad-line Balmer decrements from the Hα and Hβ reported in the WS22 catalogue for the low-redshift tail (z  <  0.57) of our sample is in excellent agreement with other samples made of allegedly unobscured sources (e.g. Zhou et al. 2006; Dong et al. 2007). In Section 4.6 we showed that the spectral shape of our global template is consistent with the VB01 template reddened by minimal levels of local extinction, i.e. E(B − V) ≲ 0.01. Similar conclusions have also been drawn in Reichard et al. (2003) who estimated the intrinsic reddening in non-BAL AGN from the early SDSS release to be E(B − V)∼0.004. This implies, under the assumption that the VB01 template is made of mostly blue unobscured quasars, a low incidence of reddening (dust/gas extinction) and galactic contamination. Finally, the similarity of the spectral shape between our sample and other quasar composite spectra from several surveys designed to select luminous blue quasars (Francis et al. 1991; Brotherton et al. 2001; VB01; Selsing et al. 2016; Harris et al. 2016) can be interpreted as a proof of the lack of significant reddening, and a clue of the non-evolution through cosmic time of the intrinsic emission properties of quasars. The other possible explanations to interpret the overall agreement among composite spectra from different samples involve a certain degree of fine-tuning. Several samples of quasars on different scales of luminosity and redshift should have intrinsically steep spectral shapes and significant amounts of dust, but still producing similar spectral shapes. As a corollary to this consideration, we should also conclude that all the quasars used to build the average SEDs available in the literature (Richards et al. 2006; Krawczyk et al. 2013; Gu 2013), whose colours have been adopted to select our sample, are also affected by the same issue.

If we were to attribute the tension with the ΛCDM in the Hubble–Lemaître diagram to dust reddening, a constant extinction with redshift would not be enough. This would require a differential extinction (increasing with redshift) combined with increasingly steeper intrinsic spectra, resulting in observed spectral shapes virtually independent of redshift.

Yet, this introduces additional questions, as larger intrinsic luminosities would also imply larger SMBH masses, enlarging the discrepancy between the largest possible SMBH masses in the Eddington-limited regime and those observed, unless the accretion efficiency varies systematically at high redshift (see, for example, Inayoshi et al. 2020 and Lusso et al. 2022 for recent reviews on this topic). In addition, significant local extinction E(B − V)∼0.1 would imply an intrinsic UV flux larger by a factor whose exact value depends on extinction curve, but on average between 3 and 4 times at 1450 Å. In turn, assuming as a rough approximation a first order constant bolometric correction, the intrinsic bolometric luminosity should be accordingly larger than the observed one. Therefore, the total mass in black holes in the Universe as evaluated by Soltan (1982), would need a significant increase with increasing redshift, unless the M_BH values are already underestimated by similar amounts. The same argument applies when considering greybody-like extinction curves (e.g. Gaskell et al. 2004), which do not modify the continuum slope, but underestimate the bolometric luminosities.

The discrepancy between the observed DM and those derived in the ΛCDM framework cannot even be solved invoking for cosmic opacity on the 2500 Å photons which would cause the observed DM to diminish. In addition, an opposite effect on DM would arise when considering cosmic opacity acting on X-ray photons. We can thus consider the decrease in the DM caused by the dimming of the UV emission as an upper limit to the detachment from the ΛCDM model caused by free electrons extinction along the line of sight. In quantitative terms, the extrapolation of the cosmic dimming values derived for SNe Ia in B band in Zhang (2008), Hu et al. (2017) lie below 0.02 at z < 3, while our Hubble–Lemaître diagram (Fig. 13) shows much larger values, as ΔDM approaches 1 at z = 3. Alternatively, some ubiquitous – yet unknown – material in the IGM should be responsible for a grey-body like absorption, though this option seems unlikely.

On the X-ray side, gas absorption plays an opposite role with respect to the UV counterpart, as the correction for absorption would increase the departure from the ΛCDM model. Generally speaking, the incidence of gas absorption among our sources is minimised by the Γ_ph  >  1.7 criterion. In addition, this effect is further mitigated at high redshift by the shift within the observed 0.5–2 keV band of much higher rest-frame energies, for which the effect of gas absorption becomes increasingly weaker. Indeed, in all the high-redshift objects for which a dedicated X-ray spectral analysis was performed (Vito et al. 2019; Nardini et al. 2019; Salvestrini et al. 2019; Sacchi et al. 2022), column densities in excess of 10²³ cm⁻² would be required to modify the spectral slope and thus affect the extrapolated rest-frame 2 keV flux, but such values were generally excluded.

As a final consideration, the results of the spectroscopic analysis described in this work ultimately prove the reliability of the photometric selection criteria adopted in our works to define a sample of quasars with typical properties. Using the photometric cuts, rather than the more time-consuming spectroscopic analysis, will allow us to considerably enlarge the cosmological sample with the forthcoming surveys (e.g. LSST, DESI), while retaining a high degree of consistency with the spectroscopic properties.

6. Conclusions

In this sixth paper of the series, we presented a thorough analysis of the average spectral properties of 2316 sources with available SDSS data from the sample described in Lusso et al. (2020). Our main aims were to investigate the possible evolution of the UV-optical properties with redshift, and to study the possible systematic effects (e.g. host galaxy contamination, dust extinction) that could undermine the validity of quasars as standard candles. To this end, we built spectral composites in the parameter space defined by M_BH, L_bol, and redshift to investigate how the spectral properties of the cosmological sample vary as a function of the key accretion parameters and cosmic time. We also carried out specific tests to determine the presence of residual galactic contamination and dust reddening based on their spectral signatures.

Our main findings are listed below:

• Apart from the line emission, which varies according to well-known relations with L_bol and M_BH, we report a strong similarity in terms of continuum shapes among the composite spectra across the explored parameter space. The slopes of the continuum of our composite spectra cluster around ⟨α_λ⟩= − 1.50 ± 0.14 without statistically significant correlations (i.e. |S|≤0.3) with the explored parameters (M_BH, L_bol, λ_Edd, z). This points to a non-evolution of the broadband continuum emission within our sample with respect to the considered parameters.

• Our global composite spectrum shows a good agreement with other templates of typical blue quasars as well as with the individual composite spectra in the various regions of the parameter space. This confirms that the L20 catalogue is, on average, made of typical objects.

• By comparing the composite spectrum of our low-redshift (z ≤ 0.8) sources with SDSS objects affected by heavy host galaxy contamination, we were able to verify that host galaxy emission is negligible in our sample. We estimate the host galaxy fraction at optical wavelengths to be ≲10%, and likely even lower at UV wavelengths.

• After performing several tests, as listed in Section 4.6, we demonstrated that the objects constituting the L20 sample are consistent with a low degree of intrinsic extinction, with an average colour excess value E(B − V)≲0.01.

• We showed that the tension between our cosmographic best fit of the joint quasar-SN Ia Hubble–Lemaître diagram, and the ΛCDM model cannot be simply accounted for by allowing for cosmic opacity or differential dust extinction. Otherwise, a colour excess E(B − V)≳0.1 increasing with z, accompanied by anomalously steep UV spectra, should ubiquitously affect high-z sources while delivering similar spectral shapes at all redshifts.

• The spectroscopic analysis carried out in this work corroborates the photometric selection criteria adopted to define a sample of quasars with typical properties. This in turn allows for the continued use of the more efficient photometric cuts rather than a more time-consuming individual spectroscopic analysis.

The search for possible systematics performed in this work has broad implications for the cosmological applications of our compilation of quasars. First, the spectral confirmation that dust or galaxy contamination do not significantly affect, without invoking peculiar or anomalous explanations, our cosmological sample further corroborates the reliability of our methods and results. In this context, quasars emerge as consistent standardisable candles, and the high-redshift deviation from the current cosmological model remains outstanding. Secondly, as our method has proven fruitful in producing large samples of homogeneous sources, we can extend and refine its application to the forthcoming data from new facilities. As a matter of fact, the extended ROentgen Survey with an Imaging Telescope Array (eROSITA, Merloni et al. 2012) began populating the pivotal z ∼ 1 region in the X-ray domain, allowing for a more robust cross-calibration between quasars and SNe Ia and a more precise determination of Ω_Λ (Lusso 2020). In the mid- to long-term future, the cross-match between Euclid (Laureijs et al. 2011) and LSST (Ivezić et al. 2019) observations in the optical–UV and Athena (and/or other proposed facilities) in the X-rays will significantly enlarge the sample, helping to fill up the region at intermediate to high redshift. New large and reliable samples will be key to testing the cosmological models in synergy with other well-established probes with the intent to precisely constrain the cosmological parameters.

Data availability

The spectrum in Fig. 3 and Table D.1 are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/689/A109

Acknowledgments

We gratefully acknowledge the anonymous referee for the thoughtful comments which resulted in a significantly improved paper. For catalogue cross analysis we have used the Virtual Observatory software (Taylor 2005), available online at http://www.star.bris.ac.uk/mbt/topcat/. The figures were generated with (Hunter 2007), a PYTHON library for publication-quality graphics. The reddening/unreddening of the spectral templates was carried out through the PYTHON package ‘dust_extinction’ (Gordon et al. 2021). For catalogue cross analysis we have used the Virtual Observatory software (Taylor 2005), available online at http://www.star.bris.ac.uk/mbt/topcat/. The figures were generated with (Hunter 2007), a PYTHON library for publication-quality graphics. The reddening/unreddening of the spectral templates was carried out through the PYTHON package ‘dust_extinction’ (Gordon et al. 2021).

References

Appendix A: Stacking methods

To build the stacks, we considered different techniques to resample the flux values, average them, and estimate the uncertainty in each spectral channel. The abbreviations in italics are used to denote these options in Table B.1.

A.1. Stacking algorithm

Here we list the stacking algorithms adopted to combine the single spectra into the composites.

Direct stacking (avg): each spectrum was rebinned on a fixed wavelength grid designed to contain all the possible wavelengths in the same redshift bin. Then, each spectrum was normalised (see below). In each spectral channel the final flux was computed as the mean, median, or geometric mean of the normalised non-zero fluxes. The uncertainty on the flux on a spectral channel was computed as the standard deviation or the inter-percentile deviation divided by the number of objects contributing to the spectral channel.

MC resampling (mc): a target number was decided before the stacking (generally N_target = 10,000). Then, we randomly extracted N_target spectra allowing for replacement and created a mock spectrum by randomizing the flux according to the uncertainty vector, in the same way as described in Section 3.3 to generate the mock spectra. The mock spectra are subject to the same procedure as in the case of direct stacking.

Bootstrap resampling (bs): a target number was decided before the stacking (generally N_target = 200). Then we randomly extracted N_draw spectra allowing for replacement and stacked them together to create the i-th composite spectrum. N_draw was chosen to be equal to the total number of good spectra in the bin, according to the quality criteria defined in Sec. 3.2. This procedure was performed N_target times, building the distribution of the mean fluxes in each spectral channel. The final spectrum in each spectral channel was obtained as the mean of the mean fluxes and the uncertainty as the standard deviation.

A.2. Normalisation

Here we list the normalisation options adopted to scale the single spectra during the combination process.

Monochromatic (mcr): the spectrum is normalised by its value at a certain wavelength $\overline{\lambda }$, generally chosen in a relatively emission-free region. Actually, the normalisation is not performed using the exact monochromatic interpolated flux, but rather the mean flux values in a wavelength range (10-100 Å) around $\overline{\lambda }$, in order to avoid anomalously high or low fluxes due to bad pixels.

Integral (itg): the spectrum is normalised by its integral over the whole wavelength range.

A.3. Average estimator

Here we list the several estimators adopted to get the average value in each spectral channel when building the spectral composites.

Mean: the direct mean is the most straightforward way to average the flux in each spectral channel. However, especially close to the wings of the spectra where the SNR lowers, composites are more prone to anomalously high fluxes, likely due to bad pixels.

Median: using the median value in each spectral channel is safer when there are strong outliers, such as bad pixels, which can significantly alter the average flux. Moreover, the median spectrum preserves the relative fluxes of the emission features (see Sec. 3 in VB01).

Geometric mean (gmean): the geometric mean spectrum preserves instead the global continuum shape. Indeed, it can be shown (Sec. 3 in VB01) that the geometric mean of a distribution of spectra describable as power laws is equivalent to a power law with a spectral index given by the arithmetic mean of the spectral indices ⟨f_λ⟩_gmean ∝ λ^{−⟨α_λ⟩}.

A.4. Uncertainty estimator

Here we list the estimators adopted to evaluate the uncertainty on the average flux in each spectral channel.

Standard deviation (std): the uncertainty is estimated as the standard deviation of the non-zero normalised fluxes divided by the square root of the objects contributing to a spectral channel.

Semi-interpercentile range (pct): the uncertainty is estimated as half the difference between the 99^th and 1^st percentile values divided by the square root of the objects contributing to a spectral channel.

Appendix B: The test bin

To test how the different stacking options could affect the final result we performed all the 36 possible combinations of the stacking options using the averagely populated (N_obj = 77) bin ‘cB4’. For the MC (bootstrap) stacks 10,000 (200) re-samples were performed.

We performed the spectral fits of all the test-bin stacks, in order to assess whether the slope of the continuum and the monochromatic fluxes were affected by the choices adopted for building the stack. In Fig. B.1 we show that α_λ, f_λ, and Γ are barely sensitive to the different stacking assumptions, as their relative variations are of the order of 2%.

Fig. B.1. Stack-to-stack variations of αλ (green), 2500-Å monochromatic flux (red), and X-ray Γ (violet). The solid black line represents the mean value, and dashed lines enclose values within a standard deviation. Stacks are numbered according to Table B.1.

Appendix C: Impact of the extinction curve

With the end to explore the effect of the choice of the extinction curve adopted to redden or de-redden the composites throughout this work, we tested different extinction curves on the composite spectrum in the ‘cB4’ bin. We included SMC-like curves (Prevot et al. 1984; Gordon et al. 2016), as well as AGN-derived extinction curves (Gaskell et al. 2004; Gallerani et al. 2010, excluding BALs; Wild et al. 2011; Zafar et al. 2015). After de-reddening the composite we performed the spectral fit to infer the slope of the continuum α_λ and compare it to that of the actual composite. We show the results of this procedure in Fig. C.1. The slope of the de-reddened spectrum assuming the Gordon et al. (2016) SMC-like curve (adopted in this work) is close to the average of the distribution. The strongest outlier in the extinction distribution is the Gaskell et al. (2004) extinction curve, which is akin to a grey-body for wavelengths below ≲3300 Å, and therefore does not alter the spectral shape at short wavelengths. However, their sample is based on 72 radio-selected quasars with the addition of mostly radio-quiet objects from the LBQS survey, and it is possible, as already suggested by Hopkins et al. (2004), that the flattening of their extinction curve is driven by the spectral properties of the radio-loud quasars in their sample. Our choice of a SMC-like curve as representative of the intrinsic quasar extinction curve seems also justified by the findings of the latter work. There, the authors found a red tail in the distribution of SDSS photometric colours according to different redshifts. This trend proved to be consistent with SMC-like reddening, but not with Large Magellanic Cloud- or Milky Way-like.

Fig. C.1. Test-bin composite spectrum de-reddened according to different extinction curves reported in the inset and assuming E(B − V) = 0.1 and RV = 3.1. The slopes produced from the spectral fits according to the different extinction curves are represented in the inset. The spectrum and the slope of the continuum before the de-reddening (αλ,  obs) are reported in magenta as references. The slope derived assuming the extinction curve from Gordon et al. (2016) is the closest to the average (black line).

Appendix D: Table of results

In Table D.1 we report all the average accretion parameters of the bins in which the parameter space is divided, as well as the spectral properties of the composites.

Appendix E: Correlation index tables

In Table E.1 we report the correlation index S between the slope α_λ and the accretion parameters, including/excluding the bins ruled out by the representativeness criterion.

Appendix F: Reddening expression

In this section, we explicitly derive Eq. 4. We begin from Eq. 7 in L20 in order to express the distance modulus DM starting from the luminosity distance d_L computed from the L_X − L_UV relation:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathrm{DM}& =5[\frac{log{f}_{\mathrm{X}}-\beta -\gamma (log{f}_{2500\mathrm{\AA }}+27.5)}{2\gamma -1}-\frac{1}{2}log(4\pi )+28.5]\hfill \\ \hfill & -5log(10\mathrm{pc}).\hfill \end{array}\end{array}$$

Keeping only the dependence on f_2500 Å, DM can be written as

$$\begin{array}{c}\hfill \mathrm{DM}=5\frac{\gamma }{2(\gamma -1)}log{f}_{2500\mathrm{\AA }}+\theta ,\end{array}$$

where θ encapsulates all the terms where f_2500 Å does not appear. Next, we assume, by hypothesis, that the observed UV fluxes are significantly extincted to the point that the discrepancy from the ΛCDM model that we report is driven by such extinction. We can then denote the observed DM as

$$\begin{array}{c}\hfill {\mathrm{DM}}_{\mathrm{ext}}=-5\frac{\gamma }{2(\gamma -1)}log{f}_{2500\mathrm{\AA },\mathrm{e}xt}+\theta ,\end{array}$$

and the intrinsic (unextincted) DM as

$$\begin{array}{c}\hfill {\mathrm{DM}}_{\mathrm{int}}=-5\frac{\gamma }{2(\gamma -1)}log{f}_{2500\mathrm{\AA },\mathrm{i}nt}+\theta .\end{array}$$

The intrinsic UV flux is related to the extincted one by

$$\begin{array}{c}\hfill f_{2500\mathrm{\AA },\mathrm{i}nt}=f_{2500\mathrm{\AA },\mathrm{e}xt}{10}^{0.4e_{\lambda }{R}_{V}E(B-V)}\end{array}$$

where e_λ is the assumed extinction curve, R_V is the total to selective extinction ratio, and E(B − V) is the colour excess. We can then write DM_int in terms of the extincted UV flux:

$$\begin{array}{c}\hfill {\mathrm{DM}}_{\mathrm{int}}=-5\frac{\gamma }{2(\gamma -1)}[log{f}_{2500\mathrm{\AA },\mathrm{e}xt}+0.4e_{\lambda }{R}_{V}E(B-V)]+\theta .\end{array}$$

The difference between the intrinsic and the observed distance moduli (ΔDM) thus becomes

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathrm{DM}={\mathrm{DM}}_{\mathrm{int}}-{\mathrm{DM}}_{\mathrm{ext}}=-5\frac{\gamma }{2(\gamma -1)}0.4e_{\lambda }{R}_{V}E(B-V).\end{array}$$

Supposing that the DM_int is the DM corresponding to the ΛCDM model, ΔDM can be evaluated by imposing DM_int = DM_ΛCDM. By isolating the E(B − V) term, Eq. 4 is found:

$$\begin{array}{c}\hfill E(B-V)=\frac{\mathrm{\Delta }\mathrm{D}M}{-\frac{5}{2}\frac{\gamma }{\gamma -1}0.4R_V\frac{A_{2500\mathrm{\AA }}}{A_V}}.\end{array}$$

Appendix G: Synopsis of the UV stacks

In Fig. G.1, we show all the UV spectral stacks across the parameter space in different colours according to their redshift bin, along with the average spectrum of our whole sample in black as a reference. Regions without any spectrum are those discarded in Sec. 3.1. Only the UV side for the spectra is shown, as the optical one is only available in the lower redshift bins.

Fig. G.1. Similarity of the stacked spectra across the parameter space. All the spectra are normalised by their average emission between 2490–2510 Å. The black spectrum represents the average spectrum of our full sample. Only spectra of the bins meeting the representativeness criterion are included in the plot.

All Tables

All Figures

	Fig. 1. Parameter space of log MBH – log Lbol for the quasar sample published in L20 with colour-coded redshifts. Dot-dashed lines represent constant λEdd values. Typical uncertainties on log(MBH) and log(Lbol) are represented in the bottom-right corner.
In the text

	Fig. 2. Examples of a low-redshift bin (left) and a high-redshift bin (right) fit. The continuum power law is marked in blue, Fe II templates in grey, emission lines in different colours (orange, purple, green, cyan). The total model is depicted in red.
In the text

Fig. 3. Spectral composite of the L20 sample. Top panel: composite spectrum of our cosmological sample described in this work. The Lyα emission line is cut for visualisation purposes. The inset zooms in on the Hβ–[O III] region. The main emission lines are labelled accordingly. Middle panel: relative residuals with respect to the VB01 template. Bottom panel: number of objects contributing to each spectral channel. The composite spectrum of the L20 sample is available in its entirety in a machine-readable form at the CDS.

In the text

Fig. 4. Comparison between our quasar template and others from the literature. All the templates are scaled by their emission at 2200 Å. The agreement on the UV side, apart from minor differences in the Fe II pseudo-continuum, specifically the 2250–2650 Å and the 3090–3400 Å bumps, is apparent. The Lyα emission line is cut for visualisation purposes. The inset shows the distribution of the absolute R magnitude against the redshift for the different samples reported.

In the text

Fig. 5. Composite spectra divided in redshift bins superimposed to the full sample composite. The dashed line denotes the number of objects contributing to the stack in each spectral channel. The composite spectra in each bin are scaled to match the template respectively at 4200 Å for the low redshift bin, 3000 Å for the z = [0.77 − 1.12), z = [1.12 − 1.52), z = [1.52 − 2.03) bins, and 2200 Å for the highest redshift bin. In the bottom-right panel we show the slopes of the continuum evaluated in different pairs of continuum windows, which are generally consistent.

In the text

	Fig. 6. Values of the slope of the continuum produced by our fits. Bins not meeting the representativeness criterion are marked with a circle. The slope of the VB01 template (dot-dashed line) is shown as an external reference. Different pairs of line-free windows were adopted to define the continuum in each bin.
In the text

Fig. 7. De-reddened stacks, according to the same colour code as in Fig. 13, scaled by their mean emission between 1770 Å and 1810 Å. For comparison, the VB01 template is also shown in black as an external reference, rather than the usual L20 template. The rising continuum slopes with increasing redshift (and extinction) are clearly observed. Inset: UV continuum slopes of the literature samples described in Section 4.1 and plotted in Fig. 4 (squares), against the slopes of the de-reddened spectra (dashed lines).

In the text

	Fig. 8. Values of Δα (Δα = 1 − αbin/αFS) in the bins of our parameter space as a function of the objects in each bin. Bins not meeting the representativeness criterion are marked with a circle. The dashed lines define the $\sim 1/\sqrt{N_{\mathrm{obj}}}$ envelope characterising the standard error on the average slope assuming a standard deviation σα = 1.
In the text

	Fig. 9. Values of the normalised 2500 Å fluxes along the axes of the parameter space, assuming a common normalisation of the continuum at 3000 Å in all the bins. Bins not meeting the representativeness criterion are marked with a circle.
In the text

	Fig. 10. Comparison between the stacks in bins of HGF and our reference spectrum (black) for the z ≤ 0.8 interval, as colour-coded in the legend. All the spectra are scaled to their 4200 Å flux. Increasingly host-dominated composites exhibit flatter UV continua.
In the text

Fig. 11. Extinction of the VB01 template according to the colour-coded values. All the spectra are scaled by their 4200-Å flux, where extinction becomes less relevant. The black spectrum is our full sample stack and is in very good agreement with a low degree of extinction. The inset shows the best-fit slope in the range 1300–4100 Å of the VB01 template reddened as a function of the E(B − V) value. The black dot-dashed line represents the slope of the continuum fitted to our stack. Minimal levels of extinction are allowed to match the reddened template with our composite. We note that the visual match between the L20 and the VB01 reddened template depends slightly on the scaling wavelength. The slope derived from the fits, however, does not depend on such a wavelength.

In the text

Fig. 12. Distribution of log(Hαbr/Hβbr) for our sources with available data. The mean values from other literature samples (Osterbrock 1977, O77, Neugebauer et al. 1979, N79, Goodrich 1990, G90, Greene & Ho 2005, GH05, Zhou et al. 2006, Z06, Dong et al. 2007, D07, Mejía-Restrepo et al. 2022, MR22, Ma et al. 2023, M23) are also shown for reference according to the colour code. The mean of our values is shown in black. The BASS sample described in Mejía-Restrepo et al. (2022) has been split according to the column density NH inferred from the X-ray analysis.

In the text

Fig. 13. Effect of the intrinsic extinction on the Hubble–Lemaître diagram. Top panel: Hubble–Lemaître diagram of the the sources used to investigate the possible effects of local reddening. The redshift bins used for this purpose are colour-coded in the legend together with the number of objects per bin and the corresponding E(B − V). Coloured stars represent the average values in each redshift bin, while black stars mark the ΛCDM model expectations, with red segments joining these two values. Increasingly high values of E(B − V) are needed in order to produce the observed tension with the concordance cosmological model. Bottom panel: difference between the predicted and observed distance modulus (ΔDM). Values are not extrapolated below z = 1.0 and above z = 4.8.

In the text

Fig. 14. Average hard X-ray bolometric corrections in each redshift bin. As a term of comparison, we report the kbol, X for the semi-analytical SED templates from M04 (blue squares) and the empirical ones from D20 (pentagons). The blue dashed line represents the best-fit reported in M04. Black and red continuous lines show, respectively, the best fit of our sample and the D20 objects in the high Lbol (≥1046 erg s−1) regime. The red dashed lines represent the 95% confidence interval of the best-fit regression line for the selected high-luminosity objects from D20 (red pentagons). Our values of kbol, X from the high-redshift, de-reddened SEDs are inconsistent with the extrapolation of the lower luminosity trend from M04, and steadily offset with respect to our best fit of the D20 high-luminosity type-1 quasars. Inset: average SEDs in each redshift bin, according to the same colour code as in Fig. 13. The shaded area is given by the individual mock SEDs in each bin.

In the text

	Fig. 15. Values of Γ for our X-ray stacks across the parameter space. Bins not meeting the representativeness criterion are marked with a circle.
In the text

	Fig. B.1. Stack-to-stack variations of αλ (green), 2500-Å monochromatic flux (red), and X-ray Γ (violet). The solid black line represents the mean value, and dashed lines enclose values within a standard deviation. Stacks are numbered according to Table B.1.
In the text

Fig. C.1. Test-bin composite spectrum de-reddened according to different extinction curves reported in the inset and assuming E(B − V) = 0.1 and RV = 3.1. The slopes produced from the spectral fits according to the different extinction curves are represented in the inset. The spectrum and the slope of the continuum before the de-reddening (αλ,  obs) are reported in magenta as references. The slope derived assuming the extinction curve from Gordon et al. (2016) is the closest to the average (black line).

In the text

	Fig. G.1. Similarity of the stacked spectra across the parameter space. All the spectra are normalised by their average emission between 2490–2510 Å. The black spectrum represents the average spectrum of our full sample. Only spectra of the bins meeting the representativeness criterion are included in the plot.
In the text