© The Authors 2024

1. Introduction

According to the unification schemes for active galactic nuclei (AGNs; Antonucci & Miller 1985; Urry & Padovani 1995; Netzer 2015; Hickox & Alexander 2018), dusty tori are a crucial component in our understanding of the AGN properties. Initial models (Antonucci & Miller 1985) describe such tori as axisymmetric structures made of dust that are found around the central super massive black holes (SMBHs), with the equivalent hydrogen column densities high enough to completely obscure the central engine in some directions. The inner radius of a torus can be constrained by the dust sublimation temperature beyond which dust grains evaporate. Given the UV emission of the innermost segments of the AGN accretion disks, this corresponds to scales on the order of parsecs and subparsecs (e.g., Suganuma et al. 2006; Koshida et al. 2014; Netzer 2015). Other unification models explore the connection between the circumnuclear dust in AGNs and the disk winds (Proga et al. 2000; Elvis 2000; Murray et al. 2005; Ricci et al. 2017). Recent studies highlight the potential for a more complex, notably highly inhomogeneous, or clumpy structure of a torus (Nenkova et al. 2008). Furthermore, these studies suggest a direct connection between the torus, the broad line region (BLR), and the accretion disk itself (Czerny & Hryniewicz 2011; Czerny et al. 2019). The study by Czerny & Hryniewicz (2011) specifically demonstrates that dust grains can survive in the accretion disk within a radial zone situated between the BLR and the inner surface of the torus. In this region the effective temperature of the disk remains below the dust sublimation temperature.

A covering factor (CF) describes the fraction of obscuration of a SMBH by a dusty torus. The initial studies defined a CF as the ratio of the solid angle Ω between the torus inner boundaries and a central SMBH to 4π, denoted as CF = Ω/4π (Hamann et al. 1993). This original definition, while accurate and straightforward, poses challenges in practical implementation using photometric observations. Currently, a widely adopted approach for CF calculation involves the ratio between the nuclear infrared luminosity, L_IR, and the bolometric AGN luminosity, L_agn, expressed simply as CF = L_IR/L_agn (Maiolino et al. 2007; Treister et al. 2008; Gu 2013; Toba et al. 2021). Three assumptions are necessary for this luminosity-based definition to correspond directly to the primary angle-based definition of the CF: (1) L_IR is dominated by the emission of the circumnuclear hot dust and depends on the amount of radiation captured from the accretion disk and reprocessed within the torus; (2) L_agn accounts for the bulk of the luminosity of the accretion disk in the active nucleus; (3) L_agn and L_IR are both intrinsically isotropic. Hence, the CF defined through the luminosity ratio should be directly proportional to the original definition of the factor. It should be noted that the values calculated from the angle-based definition are always in the range between 0 and 1. This cannot be said for the luminosity-based definition. The relation between the luminosity-based and the angle-based definitions of the CF, was addressed by Treister et al. (2008), among others.

The CF can also be calculated through various alternative methods employing different observables. One such approach involves modeling the observed X-ray continua of AGNs (Steffen et al. 2008; Brightman & Ueda 2012). Another option is to calculate the CF based as a population mean covering factor. In this approach, the CF value is derived as a ratio of type II and type I AGNs in a sample, expressed as CF = N(type II)/N(type I) (Hasinger 2008; Toba et al. 2014). The mean CF method is suitable for a sizable source sample, but is not applicable on an individual object basis. For an in-depth comparison of different methods, we refer to Netzer (2015).

The CF parameter has been widely discussed and analyzed in the literature (e.g., Maiolino et al. 2007; Treister et al. 2008; Hasinger 2008; Lawrence & Elvis 2010; Gu 2013; Toba et al. 2021). Several studies have addressed a possible evolution of the CF with redshift (La Franca et al. 2005; Treister et al. 2008; Hasinger 2008; Gu 2013; Netzer 2015; Toba et al. 2021). An anticorrelation between CF and L_agn have also been reported in several works (Gu 2013; Toba et al. 2021). In particular, Gu (2013), gathered a sample of nearly 6000 quasars from optical data in the Sloan Digital Sky Survey (SDSS) Data Releases 7 and 9 (hereafter DR7 and DR9, respectively), infrared data from the Wide-field Infrared Survey Explorer (WISE), ultraviolet data from Galaxy Evolution Explorer (GALEX), and near-infrared data from the UKIRT Infrared Deep Sky Survey (UKIDSS). They estimated L_IR and L_agn through spectral energy distribution (SED) integrals for targets in two redshift bins (2.0 ≤ z ≤ 2.4 and 0.7 ≤ z ≤ 1.1). Their analysis revealed significant anticorrelations between CF and L_agn in both redshift bins, a finding also supported by Toba et al. (2021). However, explaining the evolution of CF with redshift within the basic AGN unification scheme remains challenging (Netzer 2015). Other studies suggested no significant evolution of the CF with redshift (Gilli et al. 2007; Toba et al. 2014; Vito et al. 2018).

In this study we investigate the CF in a general population of quasars with available broadband coverage and explore the potential evolution of CF with redshift. The central question that we attempt to answer is whether the evolution of the CF is real or if the selection effects have a major influence on the CF estimates. In our study we use the largest up-to-date photometrical data sample covering the IR to UV spectral range. Based on the work of Kozłowski (2017), we estimate the masses of the SMBHs, M_BH, as well as the corresponding Eddington ratios, λ_Edd, for the objects in the sample. We explore correlations between CF and M_BH, or λ_Edd. In this context we thoroughly investigate multiple potential selection effects and signal-to-noise ratio (S/N) challenges in the infrared domain, as well as various approaches for calculating L_IR and L_agn.

The paper is organized as follows. Section 2 introduces the data and data reduction. Section 3 contains a discussion on luminosity estimates, and Sect. 4 reviews the methods used for the statistical analysis. Section 5 presents the main analysis results, which are further discussed in Sect. 6. Our conclusions are presented in Sect. 7. Cosmological parameters H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3, and Ω_Λ = 0.7 are used throughout the article. The spectral index α is defined as f_ν ∝ ν^−α, where f_ν is the flux spectral density at a frequency ν.

2. Data

The primary dataset for our analysis relies on photometric observations of quasars; specifically, we used the SDSS Quasar Catalog, Sixteenth Data Release (DR16Q¹; Lyke et al. 2020). The DR16Q contains ∼750 000 quasars, including ∼225 000 newly identified sources in comparison to previous data releases. Quasars in the SDSS DR16Q were classified based on the spectral fitting for the SDSS-IV/Extended Baryon Oscillation Spectroscopic Survey (eBOSS) spectra (Dawson et al. 2016). The final data catalog is estimated to have a completeness rate of 99.8% with contamination levels ranging from 0.3 to 1.3% (for further details, see Lyke et al. 2020). The SDSS photometric data are given in two ways, in asinh magnitudes and in nanomaggies² (Lupton et al. 1999). Measurements in nanomaggies can be converted to AB magnitudes after correcting for the shift in the u and z filters³.

The DR16Q data were split into two distinct redshift ranges. There are 110 563 targets in the 0.7 ≤ z ≤ 1.1 redshift range (hereafter the low-z sample) and 128 591 quasars in the 2.0 ≤ z ≤ 2.4 range (hereafter the high-z sample). By cross-matching with the WISE (Wright et al. 2010) All-Sky Survey, using a cross-matching radius of 2 arcsec, we identified 97 187 objects in the low-z sample and 60 056 in the high-z sample. The main selection criteria for the cross-match were significant detections in all four WISE bands W1 − W4. Further matching with the UKIDSS (Lawrence et al. 2007) DR8 returned 16 841 objects in the low-z sample and 7999 objects in the high-z sample. Finally, the low-z sample was cross-matched with the 3 arcsec matching radius with the GALEX (Martin et al. 2005) GR6/GR7, resulting in a low-z sample of 9024 objects. The main condition for the source selection here was that all the fluxes, including upper limits, were positive. In total, our final sample contains a set of 9024 low-z quasars (DR16Q × WISE × UKIDSS × GALEX), and 7999 high-z quasars (DR16Q × WISE × UKIDSS).

The WISE photometry was reported in previous work to have minor issues, such as redleak (see Sect. 2.1.1). To assess the reliability of WISE W4 (hereafter W4) photometry, we used the Spitzer Multiband Imaging Photometer for Spitzer (MIPS) 24 μm (hereafter M24) observations (Spitzer Science Center (SSC) & Infrared Science Archive (IRSA) 2021). The Spitzer M24 observations were available only for 290 low-z objects and 171 high-z sources. The entire data selection process is presented as a flow chart in Fig. 1.

Fig. 1. Flowchart showing the data selection process, as described in Sect. 2. The final samples are 1534 objects for low-z with S/NW3 & W4 > 5 and 309 objects for high-zS/NW3 & W4 > 5. The low-z and high-z dataset combined create the low-high-z. High-precision datasets with Spitzer observations have 198 objects for low-zSpitzer and 91 objects for high-zSpitzer, both with S/NW3 > 5. The datasets with different S/NW3 & W4 cuts are analyzed in Appendix C.

All the photometric data underwent correction for Galactic extinction using the reddening maps of Schlegel et al. (1998), and the extinction law of Cardelli et al. (1989). These data were transformed into rest-frame luminosity spectral densities $L_{\nu }=4\pi d_{\text{L}}^2{f}_{\nu }/(1+z)$, where f_ν represents the extinction-corrected flux, and d_L signifies the luminosity distance calculated using the chosen cosmological model based on the code of Wright (2006).

The mean SEDs for the low-z and high-z sources are presented in Fig. 2. As shown, the slopes of the low- and high-frequency segments of the SEDs are different between the low-z and high-z samples. In addition, the high-z quasars have higher monochromatic luminosities, especially in the IR range. Given the overall concave shape of the SEDs, hereafter we define the torus-related IR luminosity L_IR as that corresponding to the integrated continuum emission below the rest-frame 1 μm, and the accretion-related bolometric (optical–UV) luminosity L_agn corresponding to frequencies above 1 μm (see Fig. 2).

Fig. 2. Mean rest-frame SEDs for low-z quasars (top) and high-z quasars (bottom). The top and bottom panels depict the mean rest-frame SEDs for low-z and high-z quasars, respectively. The black and green shaded regions highlight the areas utilized for computing LIR and Lagn, respectively, in the all-points method. Each photometric observation is labeled with the filter name. The error bars denote the 1st and 3rd quartiles.

2.1. WISE data check

2.1.1. Redleak

Prior research has highlighted minor issues with the WISE data, such as redleak, which have been appropriately addressed in our analysis. In particular, following the methodology outlined in Brown et al. (2014), the effective wavelength of W4 was changed from 22.1 to 22.8 μm, and the AB magnitude of the Vega zero point of W4 was changed from m_W4 = 6.59 to m_W4 = 6.66. As follows, this correction is relatively minor.

2.1.2. Signal-to-noise ratio

For the faint objects in the sample, the S/N is notably low and the measured W3 and W4 fluxes have relatively large errors, making some of the observations unreliable. To address this issue, we imposed a stricter requirement for the S/N to exceed 5σ for the W3 and W4 filters. It is worth noting that the equivalent restriction on W1 and W2 does not alter the amount of data after the condition S/N_W3 & W4 > 5 is set, as the W3 and W4 data present more challenges in this context. In the end, there are 1534 high-S/N objects in the low-z sample and 309 in the high-z sample. Low-z and high-z combined together are called low-high-z. This dataset has 1843 objects. This associated subset of data was subsequently examined for possible alterations in the source parameter distributions, and a possible evolution of the CF (see Appendices A and C for further details). A similar restriction is applied to the associated Spitzer dataset, utilizing only W3 from the two problematic filters. When S/N_W3 > 5 is applied, the Spitzer low-z sample contains 198 objects, while Spitzer high-z contains 91 objects. These two sets, when combined, are collectively referred to as the low-high-zSpitzer dataset, consisting of 289 objects. Going forward, unless explicitly stated otherwise, the presented data adhere S/N_W3 & W4 > 5.

2.2. Value added catalogs

To estimate the physical parameters of the analyzed sources, such as the SMBH mass M_BH or Eddington ratios λ_Edd, high-quality spectroscopic data are essential. For this we utilized the value-added catalog (VAC) for the SDSS by Kozłowski (2017). This catalog is based on spectral line fitting using the data from SDSS DR12Q (Pâris et al. 2017). In particular, Kozłowski (2017) provided virial estimates for M_BH, using the BLR radius–disk luminosity relations, and the full width at half maximum (FWHM) of the Mg II and C IV emission lines. This approach relies on the phenomenological scaling relation $R_{\mathrm{BLR}}\propto L_{\mathrm{disk}}^{1/2}$ connecting the characteristic spatial scale of the BLR, R_BLR, with the continuum disk luminosity, L_disk. This relation is established through reverberation mapping (McLure & Jarvis 2002; Kaspi et al. 2007; Bentz et al. 2013). Utilizing this estimate, the virial theorem can be used with the relation M_BH ∝ R_BLR v², where v is the characteristic velocity of the BLR clouds inferred from the FWHM of the broad emission lines (Vestergaard & Peterson 2006; Shen 2013). We note that single-epoch viral M_BH estimates come with methodological biases, which are discussed in Shen et al. (2016) and Kozłowski (2017).

Based on the obtained values of M_BH, the Eddington luminosities were calculated as

$$\begin{array}{c}\hfill L_{\mathrm{Edd}}=\frac{4\pi GM{m}_{\mathrm{p}}c}{{\sigma }_{\mathrm{T}}}\simeq 1.26\times {10}^{38}\left(\frac{M_{\mathrm{BH}}}{M_{\odot }}\right)\mathrm{erg}{\mathrm{s}}^{-1},\end{array}$$(1)

and the corresponding logarithmic Eddington ratios as

$$\begin{array}{c}\hfill {\lambda }_{\mathrm{Edd}}=log\left(\frac{L_{\mathrm{agn}}}{L_{\mathrm{Edd}}}\right)·\end{array}$$(2)

We opted to exclusively utilize the Mg II-based estimations of M_BH from Kozłowski (2017). The cross-matching process with our S/N_W3 & W4 > 5 WISE sample returned 229 low-z and 125 high-z sources, and with the SpitzerS/N_W3 > 5 sample 41 and 33 targets, respectively. Subsequently, we derived the Eddington luminosities and Eddington ratios from these data.

3. Luminosity estimates

Using the comprehensive multiwavelength data we gathered, we precisely estimated the infrared and bolometric luminosities of quasars via the all-points method. This method offers an alternative to the commonly employed power-law method found in the literature (e.g., Gu 2013). Specifically, we utilized the trapezoid method to integrate the quasar SEDs. This involved fitting a line in log-log space between consecutive photometry data points and then integrating the area beneath it, as opposed to using a fixed wider frequency range. Furthermore, we incorporated interpolated luminosity densities at wavelengths of 0.11 μm (L_0.11 μm), 1 μm (L_1 μm), and 7 μm (L_7 μm). This allowed us to cover similar ranges of the rest-frame continua in both the low-z and high-z samples. For L_7 μm we interpolated values from the W4 and W3 WISE filter detections in both low-z and high-z samples. The value of L_1 μm was calculated based on interpolation between the W1 and K filters for low-z or W2 and W1 for high-z. The value of L_0.01 μm was calculated using interpolation between near-UV and far-UV for low-z, or g and u for high-z.

We defined L_IR as the integrated area between L_7 μm and L_1 μm, while L_agn was defined as the integrated area between L_1 μm and L_0.11 μm. The rationale for this definition is that in the AGN unification scheme, the dusty torus is illuminated by the AGN disk radiation and captures a fraction of the total energy proportional to the CF. The absorbed optical, UV, and soft X-ray radiation is then reemitted as a thermal IR radiation. In Appendix D.1 we discuss in more detail the applied luminosity integration procedure, in particular confronting it with the standard (power-law) method.

We reiterate our choice of estimating quasar infrared and bolometric luminosities using a model-independent method; we simply integrated the observed SEDs of the selected sources. An alternative approach could involve using templates to initially fit quasar SEDs and then deriving the corresponding IR and bolometric luminosities based on the best-fit values of the model parameters obtained. One clear advantage of the template fitting method is its consideration of intrinsic source extinction. However, it introduces systematic model uncertainties, such as those related to the specific extinction law assumed. A more detailed discussion of template fitting for the analyzed samples of quasars will be provided in a follow-up article.

The CF values were next calculated in the same manner for all the objects in the sample as

$$\begin{array}{c}\hfill \mathrm{CF}=\frac{L_{\mathrm{IR}}}{L_{\mathrm{agn}}}·\end{array}$$(3)

The resulting CF range is characterized by a wide spread, with the maximum values exceeding unity; we note again that such values are not allowed within the framework of the angle-based definition CF = Ω/4π.

4. Methods

4.1. Bayesian fitting

To quantify the differences between the low-z and high-z quasars in our sample, Bayesian regression analysis was performed. We used the Python library emcee⁴, which is the implementation of the Markov chain Monte Carlo (MCMC) ensemble sampler (Foreman-Mackey et al. 2013). The statistical model, M, in the form y = m × x + b was used in the likelihood function, with the model parameters Θ = (m, b). For the response variable we choose y = log L_IR or y = log L_IR/L_agn, and for the prediction variable x = log L_agn. With this approach we determined the posterior probability density function (PDF) via the Bayes theorem

$$\begin{array}{c}\hfill P(\mathrm{\Theta }|D,M)=\frac{P(D|\mathrm{\Theta },M)P(\mathrm{\Theta }|M)}{P(D|M)},\end{array}$$(4)

where D stands for data, P(D|Θ, M) is the likelihood, P(Θ|M) is the prior, and P(D|M) is the model evidence. The prior space was defined between −2 < m < 2, and −35 < b < 35. The starting values of the priors are based on the ordinary least-squares (OLS) regression fitting. We chose the log likelihood in the form

$$\begin{array}{cc}\hfill P(D|\mathrm{\Theta },M)=& lnp(y_n|x_n,s_n,m,b)=-{\displaystyle \sum _n[\frac{{(y_n-m{x}_n-b)}^2}{s_n^2}}\hfill \\ \hfill & \times ln\left(\frac{{(y_n-m{x}_n-b)}^2}{s_n^2}\right)],\hfill \end{array}$$(5)

where s_n is the data error (in this context, see Appendix B).

The resulting best-fit regression lines and parameters are given as the mean relations, at the 68% confidence levels, in Figs. 3a, 5, 9, and 10, along with the relations obtained in the literature for comparison.

Fig. 3. Relation between log LIR and log Lagn with S/NW3 & W4 > 5. The blue and orange circles indicate low-z quasars, whereas the green and red diamonds indicate high-z sources. The orange circles and red diamonds indicate low-z and high-z quasars with Spitzer M24 data. In panel a the blue line denotes the best-fit Bayesian regression for the low-high-z data with S/NW3 & W4 > 5, the orange line is the Bayesian regression for the low-high-zSpitzer data from both low-z and high-z, and the green and red dashed lines represent the best Bayesian fit with the low-z sample weighted and low-z with Spitzer weighted, respectively; the black line gives the 1:1 scaling relation between LIR and Lagn. In panel b the big black circles represent the medians calculated for the 0.5 dex in log Lagn for both redshift samples (low-high-z) and the orange big triangles stand for the medians calculated for both redshift samples with the Spitzer data (low-high-zSpitzer). For the black and the magenta points, the error bars represent MAD. The exact values of the fitted parameters are listed in Table 1.

4.2. Dispersion statistical analysis

4.2.1. MCMC errors

The errors in the estimates of L_IR and L_agn (as discussed in Sect. 3) were calculated using Monte Carlo (MC) simulations based on the observed fluxes and their uncertainties. Due to the asymmetric distribution of simulated luminosities, two-sided errors were defined as the 16th and 84th percentiles. These MC errors were considered throughout the entire analysis, especially in the regression analysis.

4.2.2. Intrinsic scatter

Relations between physical quantities are hardly scatter-free. Inside each relation, there are perturbations involved (e.g., in our case dependence on other physical properties) that cause the scatter of this relation to some degree. That means that when one wants to measure the physical quantity θ, in reality $\varphi =\theta +\stackrel{\sim }{\theta }$ is measured, where $\stackrel{\sim }{\theta }$ is the perturbation.

To quantify the intrinsic scatter in the analyzed correlation, denoted as σ_i, we followed the procedure described in Stone et al. (2021, see Sect. 2 therein). Specifically, within the Bayesian framework, the PDF of σ_i can be calculated as

$$\begin{array}{c}\hfill P({\sigma }_i|\varphi ,{\sigma }_{\varphi })=\frac{P(\varphi |{\sigma }_i,{\sigma }_{\varphi })P({\sigma }_i)}{P(\varphi )}·\end{array}$$(6)

Here P(σ_i) represents the prior distribution of the intrinsic scatter, while P(ϕ) is the normalization corresponding to the measured quantity ϕ with its associated systematic and statistical error, σ_ϕ.

4.3. Efron & Petrosian test

When working with truncated data, one can never be sure whether the dependence between the analyzed parameters is real or results from selection biases in the sample. Our data sample is built from several surveys, each with different limitations and selection criteria. The most important of these effects is the flux-limitation of the sample. Handling such effects requires a statistical test to estimate how much the truncation (flux-limitation in our case) affects the final dependence between the luminosities and the redshift. Efron & Petrosian (1992, hereafter EP) developed a nonparametric statistical test to determine the existence of a correlation or independence of variables from the truncated dataset, with the special case being flux-limited data. A detailed description of each step can be found in Maloney & Petrosian (1999); here we only briefly describe the most important assumptions and steps leading to the calculation of the EP test.

If one deals with one-sided truncation of the data, as in the luminosity-redshift case in a flux-limited survey, the following steps are required by the EP test. First, for every object i in the data sample, one should create an associated set, $J_i=\{j:L_j>L_i,L_j^{-}<L_i\}$, where $L_j^{-}$ is the limiting luminosity (calculated for each object based on the chosen limiting flux), L_i is the luminosity of the considered object, and L_j is the luminosity of the object j to which we compare the object i. The set J_i contains all the objects that have luminosity above the chosen flux limit, and luminosity higher than the luminosity of the object i. Based on the set J_i, the rank R_i for each object i is calculated following the equation R_i = number{j ∈ J_i : z_j ≤ z_i}, where z_i is the redshift of the object i and z_j is the redshift of the object j from the set J_i. If luminosity and redshift are independent, then the rank R_i should have a uniform distribution between 1 and N_i, where N_i is the number of objects in J_i.

The test statistic is then defined as

$$\begin{array}{c}\hfill \tau =\frac{{\sum }_i(R_i-E_i)}{\sqrt{\sum _i{V}_i}},\end{array}$$(7)

where $E_i=\frac{1}{2}(N_i+1)$ is the mean of the uniform distribution and $V_i=\frac{1}{12}(N_i^2-1)$ is the variance. The hypothesis of independence is rejected between luminosity and redshift if |τ|≥nσ, where n can be chosen as the σ multiplier for a given confidence level.

5. Results

This section presents the results of our correlation analysis between AGN and IR luminosities in the constructed quasar samples. Additionally, we investigate the dependence of the CF on luminosities and present the results of the EP test. Finally, we assess the potential evolution of the CF with redshift.

5.1. Relations between luminosities

Figure 3 illustrates the relationship between L_IR and L_agn for the analyzed samples, where we compare the luminosities calculated with the WISE data to those derived based on Spitzer M24. The lines in the figure denote the results of the Bayesian regression analysis, for which the obtained model parameters are summarized in Table 1. The data were restricted to S/N_W3 & W4 > 5 (for the entire WISE sample, see Appendix C). Notably, we observe more WISE outliers despite the applied high-quality criterion, while the Spitzer sample exhibits greater condensation, particularly in the low-z range. The low-z and high-z samples tend to form a single correlation trend, and thus we decided to fit a single power law to the low-z and high-z sources simultaneously, which we call the low-high-z sample. The best-fit regression lines can be found in Fig. 3. The resulting power slope for the low-high-z data is 0.77 ± 0.01. For the low-high-zSpitzer data, where W4 was replaced with MIPS 24 μm, the index is 0.88 ± 0.01. These values exhibit only slight differences, and visually, the fits closely align. The major feature here is the lack of low-luminosity objects in the top panel of the figure (compared to the whole sample without the S/N restriction, see Fig. C.1a). This can be explained by the fact that the objects with the lowest luminosities have unreliable measurements in terms of the S/N in the W3 and W4 filters.

The S/N_W3 & W4 > 5 low-z sample is almost five times larger than the high-z sample. This can cause the regression fit to be largely determined by the low-z sources. To mitigate this effect, we conducted a simple test. From the low-z sample, 20% of the whole S/N_W3 & W4 > 5 data points were randomly selected (“low-z reduced”). A similar procedure was performed for the low-zSpitzer data, but only 50% of the low-zSpitzer sample with S/N_W3 & W4 > 5 was taken into account (“low-zSpitzer reduced”). The low-z reduced sample was then merged with the high-z sample with S/N_W3 & W4 > 5 (“low-high-z weighted”); likewise the low-zSpitzer reduced with the high-z sample with S/N_W3 & W4 > 5 (“low-high-zSpitzer weighted”). The Bayesian regression analysis was repeated for each of these combined samples, as outlined in Sect. 4. This procedure was executed 50 times and the final parameters were taken as the median of the 50 fits. This test is described as the weighted fit and is represented by the green line in Fig. 3a). The resulting slopes for the weighted data are 0.87 ± 0.01 (low-high-z), and 0.90 ± 0.01 (low-high-zSpitzer). This means that the larger low-z sample dominates the fit, in particular the tail below the dex 45.5 in log L_agn.

To address the issue of different sample sizes between the low-z and high-z sources, a second test was conducted using binning for the entire dataset. The bins were created as follows. Starting from log L_agn = 44.5, the subsequent bins were created every 0.5 dex of log L_agn. For each bin the median values and the median absolute deviations (MADs) were calculated both in log L_agn and log L_IR. This procedure was applied separately for the low-high-z and low-high-zSpitzer samples, and the results are given in the inset at the bottom of panel b of Fig. 3. The problematic low-luminosity segment of the low-z sample is now even more visible, exceeding the 1:1 scaling. A similar trend is observed in the Spitzer sample. However, due to a larger number of outliers toward higher L_IR, the low-luminosity bins in the WISE sample diverge more significantly from the overall trend, while in the Spitzer sample, the situation is far less severe. It is worth noting, however, that the lowest luminosity bin contains just a few sources in both the WISE and Spitzer samples, and this is reflected by a large error bar in the WISE sample, even though in the corresponding Spitzer bin the variance is much smaller.

To have some insight into how the restrictions on the S/N influence the sample, one can refer to the medians and their errors, as quantified in Table 2. It is evident that as we move toward higher S/N in the WISE filters, all the median luminosities increase. This is expected since more luminous sources usually have higher flux in a given redshift bin. The deviation in the low-z sample decreases with increasing S/N from 0.29 in log L_agn and 0.24 in log L_IR, to 0.19 and 0.12, respectively (i.e. it approaches the level observed in the high-z population). The high-z sample is significantly smaller, however, consisting of 309 high-S/N objects compared to 7999 without the S/N restriction, and so increasing the S/N does not result in this case in decreasing the variance. In Fig. 4 we summarize the best-fit slopes (panel a) and intercepts (panel b) for the five best Bayesian fits for the WISE data, S/N_W3 & W4 > 3, S/N_W3 & W4 > 5, Spitzer and the S/N_W3 & W4 > 5 samples with the weighted low-z sample. The general trend in the data with WISE W4 (blue squares) is that the higher the S/N, the lower the slope value, so the deviation from the 1:1 relation between L_IR and L_agn becomes larger.

Fig. 4. Comparison between the Bayesian regression analysis model parameters: (a) slopes, (b) intercepts, obtained for different analyzed samples of quasars (see Table 1 and Sect. 5.1).

As mentioned earlier, the low-z and high-z samples follow a similar correlation between L_agn and L_IR. In the case of the best-quality data with Spitzer, the obtained best-fit regression lines in logarithmic space only show slight differences, as illustrated in Fig. 5. For data without Spitzer the difference is more apparent, in particular the low-z fitted slope is 0.70 ± 0.01, whereas the high-z slope is 0.51 ± 0.01. A shallower high-z correlation line diverges at the high-luminosity tail of the sample. This is due to the majority of the members lying at bolometric luminosities below 46.5 dex.

Fig. 5. Regression analysis results for log Lagn and log LIR for low-z, low-zSpitzer, and high-z, high-zSpitzer sources separately. The model parameters of the Bayesian linear regression analysis are given in Table 1.

5.2. Relations between CF and luminosities

Figure 6 shows the relations CF versus L_IR and CF versus L_agn. Despite precise estimates of both luminosities, a substantial spread is still evident, particularly in the log CF − log L_IR relation. The corresponding medians, calculated in the same way as in Fig. 3, are also displayed in Fig. 6. In panel a, an anticorrelation between log CF and log L_agn is apparent, consistent with the regression fit presented in Fig. 3. What is particularly evident here is the influence of the low-luminosity segment of the low-z sample, below ∼log L_agn/erg s⁻¹ = 46. All luminosities are in cgs units, unless stated otherwise. The low-z quasars (blue dots in the figure) in general dominate up to the bin at log L_agn ∼ 46.2; the high-z quasars (green points) dominate above 46.5 dex in log L_agn. This is especially visible for the median calculated for both the low-z and high-z samples together (low-high-z: big black dots). A similar trend is present in the low-z and high-zSpitzer subsamples. On the other hand, no anticorrelation between CF and L_IR is present in panel b of Fig. 6. The medians here, calculated with the bins every 0.5 dex log L_IR, form a horizontal sequence, with the only outliers being the low-luminosity tail of the low-zSpitzer sample. The differences between the low-z and high-z samples in Fig. 6a are likely caused by instrumental effects. In particular, the low-luminosity bin in the high-z sample (first green diamond), which has the highest CF, would simply be absent with the increased S/N threshold. Meanwhile, the first bin in high-zSpitzer (red diamond) has the highest variance, with many outliers. The second point in high-z also has higher values than the trend for both samples combined, which is not the case for the second bin in the corresponding Spitzer sample.

Fig. 6. Relationships between CF and LIR (upper panel a) and Lagn (lower panel b), including only the sources with S/NW3 & W4 > 5. The blue circles indicate low-z quasars, and the green and red diamonds indicate high-z sources. The orange circles and red diamonds indicate low-z and high-z quasars with the Spitzer M24 data. Additionally, the median statistics were calculated for the samples, as described in the legend in the lower left corner. The error bars represent the MAD errors.

5.3. Testing for the luminosity evolution

The challenges in addressing the luminosity evolution of quasars in the analyzed samples arise from different flux limits in various survey releases and the integration of source luminosities across broad frequency ranges. The SDSS DR16Q data were folded from DR7Q and DR12Q, and each of these releases had a different truncation and magnitude limit. In the end, we chose r = 22 mag as the limit for L_agn, which was the most restrictive limit among the SDSS DR16Q data, and W4 = 6 mJy for L_IR, the most restrictive limit of the IR photometry. The fluxes were then converted to monochromatic luminosities in erg s⁻¹ units. Next we followed the steps of the EP test described in Sect. 4.3 for the joint low-z and high-z population, working on (a) the low-high-z sample without S/N cuts, (b) the sample with S/N_W3 & W4 > 3, (c) the low-high-z sample with S/N_W3 & W4 > 5, and (d) the low-high-zSpitzer sample. For each of these subsamples, the EP test was performed separately for L_IR and L_agn. The resulting test values of τ are (a) τ_IR = 91.6, τ_agn = 103.9; (b) τ_IR = 54.4, τ_agn = 53.0; (c) τ_IR = 33.1, τ_agn = 30.4; and (d) τ_IR = 13.8, τ_agn = 15.6. In all these cases the hypothesis of independence between L_IR and z, and between L_agn and z should therefore be rejected, and so a strong redshift evolution of the luminosities in both bands is observed.

5.4. Testing for the CF evolution

The EP test results described earlier indicate the evolution of both L_IR and L_agn with redshift. Although relevant, this result still does not answer the question of a possible evolution of the CF. To address this issue, we selected quasars from our sample of the objects with comparable M_BH values. Quasars selected in this way should have corresponding physical properties regardless of their redshifts. Moreover, in addition to the spectroscopic-based virial M_BH estimates, we also required good-quality Spitzer M24 data for the targets. Finally, only the objects with the M_BH values above the threshold of either log M_BH/M_⊙ = 8.5 (which is the lowest value for the high-z population) were selected within the subsample with S/N_W3 & W4 > 5. The final data sample for this test consists of 59 objects, including 26 low-z objects, and 33 high-z sources. Figure 7, panel a shows the resulting distribution of the CF values with M_BH. As can be seen, the analyzed dataset reveals no significant differences between the low-z and high-z targets in the projected log CF distribution. To quantify this statement, a one-dimensional two-sample Kolmogorov–Smirnov test was performed. The resulting p-value = 0.80 is high, meaning that the hypothesis that the two samples were drawn from the same distribution (the null hypothesis) should not be rejected, at a very high confidence level. In the final step, we performed the two-dimensional generalization of the two-sample Kolmogorov–Smirnov test (Fasano & Franceschini 1987) for the log M_BH − log CF distribution. The p-value in this case turns out to be relatively small, below 0.01, implying that the hypothesis that the two-dimensional sample is from the same distribution can be rejected at the > 99% confidence level. Additional tests, however, revealed that this was due to a shift in the mean M_BH values between our low-z and high-z populations, consistent with a trend observed earlier and discussed in McLure & Dunlop (2004).

Fig. 7. Relations of the CF with redshift. The top panel shows the relation between log CF and log MBH for quasars with the Spitzer M24 and VAC data. The red line represents the cut on log MBH/M⊙ = 8.5. The bottom panel shows the relation between log CF and redshift for quasars with log MBH/M⊙ > 8.5. The orange and red points in both panels represent the low-z and high-z quasars, respectively. The large symbols with MAD error bars stand for the medians calculated for each data sample: low-z with no cut in MBH (red triangle), high-z (orange triangle), low-z with the cut in log MBH/M⊙ > 8.2 (pink triangle), low-zSpitzer with no cut in MBH (black square), high-zSpitzer (green square), low-zSpitzer with the cut in log MBH/M⊙ > 8.5 (blue square). For the luminosity and Eddington ratio selection, see Appendix E.

Figure 7, panel b shows the distribution of the CF with respect to the redshift for the low-zSpitzer sample (after the cut in M_BH) and the high-zSpitzer sample. Additionally, the medians in both log CF and z for each data sample were calculated, as represented in the figure by large symbols with error bars. As follows, all the medians are consistent within the errors with the scenario of no significant evolution of the CF with redshift. As discussed in more detail in Appendix E, similar results also hold for quasars selected based on luminosities or Eddington ratios.

5.5. Dispersion analysis

The intrinsic scatters σ_i in the log L_IR versus log L_agn relations in the primary datasets with S/N cuts were analyzed, as described in Sect. 4.2. The results are summarized below:

• (i)

  low-zS/N_W3 & W4 > 5: σ_i = 0.157;

• (ii)

  high-zS/N_W3 & W4 > 5: σ_i = 0.165;

• (iii)

  low-high-zS/N_W3 & W4 > 5: σ_i = 0.208.

For the Spitzer data, on the other hand, we obtained

• (iv)

  low-zSpitzerS/N_W3 > 5: σ_i = 0.150;

• (v)

  high-zSpitzerS/N_W3 > 5: σ_i = 0.222;

• (vi)

  low-high-zSpitzerS/N_W3 > 5: σ_i = 0.148.

The σ_i values are very similar in both low-z and low-zSpitzer. The difference is visible for the high-zSpitzer sample, where the variance is large as the relative number of outliers is increased, especially at the low-luminosity end of the high-z sample (Fig. 5).

The second test aimed to assess variances in both fitted regression and observational errors. This test was adapted from a similar approach presented in the literature (Risaliti & Lusso 2015), with the necessary modifications. To ensure a cosmology-free comparison, fluxes were used rather than luminosities. Specifically, we utilized the SDSS r (the most limiting filter in L_agn present in both low-z and high-z with and without Spitzer) and WISE W2 filters (the most limiting filter in L_IR). Figure 8a shows the relation between the logarithms of the r-band and W2 extinction-corrected fluxes, log f_r and log f_W2, respectively.

Fig. 8. Analysis of the intrinsic scatter. Top panel: scatter plot for the extinction corrected log fr and log fW2, for low-z, high-z, low-zSpitzer, and high-zSpitzer. Bottom panel: standard deviation for binned fluxes in the bands r and W2, and the error of log r − W2 relations σr W2 scatter. The bins were calculated as equally distanced in log r space. The sizes of the symbols correspond to the log fr flux values in each bin.

The data were spread into ten bins in log f_r, with equal distances in logarithmic space. For each bin the standard deviations σ were calculated considering observational (flux measurement) uncertainties, namely ${\mathrm{\Delta }}_{r-W2-\mathrm{error}}=\sqrt{\mathrm{\Delta }{(log{f}_r)}^2+\mathrm{\Delta }{(log{f}_{W2})}^2}$. Moreover, for the given flux distributions (as given in Fig. 8a) we also calculated ${\sigma }_{rW2\mathrm{scatter}}=\sqrt{\sigma {(log{f}_r)}^2+\sigma {(log{f}_{W2})}^2}$. We repeated this process for all four datasets (low-z, high-z, low-zSpitzer, high-zSpitzer). The results of the test are presented in Fig. 8b, where the size of the markers is coded to represent the mean r-band flux of each bin; the larger the points, the higher the flux values.

From this test one can see first that the trend with σ_{r W2 scatter} is only visible for the lower luminous bins (with higher flux the σ_{r W2 scatter} is lower); the higher luminosity bins have no significant trend. Second, on average, the dispersion of the relation is higher than the observational errors for all datasets. The dispersion of relation σ_{r W2 scatter} is comparable with the values of σ_i. In our opinion, the two proposed tests give similar results.

6. Discussion

Our detailed analysis has confirmed that calculating quasar infrared and bolometric luminosities by integrating directly all available photometric data points is more accurate than when using a power-law approximation (Figs. D.1 and D.2). This finding leads to an improved correlation between L_agn and L_IR, as presented in Fig. D.1, especially in the high-z sample, for which the power-law method underestimates the L_agn values and overestimates L_IR. The integral method also reduces the L_IR dispersion in the low-luminosity segment of the high-z distribution, and improves the consistency between the low-z and high-z samples in the overlapping luminosity range.

The bias that remains is caused by the accuracy of the IR photometry. This was tested for the least accurate filter W4, in particular by comparing it with the Spitzer MIPS 24 μm whenever possible. In many cases the Spitzer photometry gave systematically lower fluxes than those measured in WISE W4. Furthermore, the Spitzer data have a smaller overall spread (Fig. C.1, top panel). Imposing a data quality cut on the W3 and W4 filters by setting S/N > 5 improves the situation in this respect by tightening the luminosity-luminosity correlations (see Figs. 3, A.1, and A.2). Although the much reduced number of sources with the available M24 data limits the scope of the analysis, the medians of the primary parameters in the Spitzer-detected quasars are more representative for the entire sample than medians calculated with the high-quality WISE data (S/N_W3 & W4 > 5), as given in Table 2.

We also checked the variance versus error for L_agn and L_IR following the method by Risaliti & Lusso (2015). It seems that the lowest variance is in the medium error bin. This means that the lowest spread is present near the middle of the distribution. The medium luminosities also have the highest density of members in the bin.

Regarding the analyzed samples including either WISE or Spitzer data, the major differences follow from WISE W4 versus Spitzer MIPS 24 μm photometry, as discussed in detail in Appendix A.2. We believe that the differences are caused by instrumental effects rather than source variability. In particular, there is a systematic difference in the distribution of data points on the log L_agn − log L_IR diagram between the W4 and MIPS24 samples, with the higher-luminosity sources being underrepresented in the Spitzer data.

Our main point of interest was to investigate the distribution of the covering factor of the dusty torus within the quasar population. For that purpose, we adopted a simple definition of the CF parameter as the L_IR-to-L_agn ratio. When we consider the entire sample of sources binned in log CF and log L_agn, the bins above 45.5 in log L_agn form a horizontal sequence with log CF ≃ −0.2 (see Fig. 6, panel a). The difference between the WISE and Spitzer data emerges here only at the lowest luminosity bins, for which the median CF values are in general higher and noticeably divergent (although still within the errors). On the other hand, in the log CF − log L_IR binned distribution, we see a constant value of log CF ≃ −0.2 throughout the entire luminosity range.

When considering separately the low-z and high-z samples, in both cases we observe an anticorrelation between log CF and log L_agn. The low-luminosity tail of the low-z sample is however sensitive to the adopted instrument at 24 μm: the slope coefficient of the line fitted to the W4 data is almost twice smaller than that obtained with the Spitzer data. In Fig. 9, we compare our regression results in this respect with the results by Toba et al. (2021), who performed the broadband quasar SED fitting including the dusty torus and polar dust. In their analysis, the CF of the torus exclusively turns out weakly dependent on L_agn, but when the polar dust is added as the SED model component, the anticorrelation becomes steeper and, in fact, basically the same as observed in our low-zSpitzer subsample. This may suggest that the integrated IR quasar luminosities, as considered in our work, are dominated by the torus in the high-luminosity regime, and are affected by a significant contamination from the polar dust in the low-luminosity low-z regime.

Fig. 9. Regression analysis results for the covering factor and bolometric luminosity for the low-z quasars with S/NW3 & W4 > 5. The blue and orange lines are the same as in Fig. 5. The green line corresponds to the regression from Toba et al. (2021) fitted to the model with torus dust only, while the red line is the sum of regressions fitted to the torus and polar dust models separately.

Somewhat puzzling, in this context, is that the CF in the high-z sample alone also decreases with L_agn. This is clearly visible in Figs. 6 and 10. An interesting finding here is that, with the Spitzer data, we obtain different slopes of the regression lines for the low-z and the high-z subsamples (−0.17 ± 0.01 and 0.12 ± 0.01, respectively; see Fig. 10). The low-z relation can be accounted for, as argued above, by the polar dust contamination, with the log CF dropping from 0 to −0.25 within the log L_agn range from 45 to 46, while in the high-z sources the same drop occurs on log L_agn from 46 to 47. It is important to note that, unlike low-z, the high-z fits for luminosity ratios using different methods give inconsistent results. It is clearly visible in Figs. 10 and B.1, where W4 versus M24 data are used or unweighted versus error-weighted regression is used. This can suggest that the data quality and distribution do not allow for a good constraint on the slope parameter. Thus, it is uncertain what physical behavior is present in the high-z sample CF.

Fig. 10. Relation between log CF and log Lagn for the low-z and high-z sources with S/NW3 & W4 > 5 (upper and lower panels, respectively). The regression lines in both panels fitted to our various subsamples, with and without the Spitzer data, along with the Gu (2013) regression lines for comparison, are described in the legends of the two panels.

Ricci et al. (2017, 2022) pointed out that the column density of material in the torus or material above the torus anticorrelates with the Eddington ratio. In their radiation-regulated unification scenario in low-accretion rate AGNs (λ_Edd ≥ −2) the solid angle is larger, which suggests thicker tori. The study focused on constraining X-ray absorption in the local low-luminosity AGN population which is well represented in the Swift Burst Alert Telescope (BAT) catalog. This finding is generally compatible with studies by other authors based on SED fitting. For example, Toba et al. (2021) found an anticorrelation of the covering factor with AGN luminosity and pointed out that polar dust emission adds up to the torus emission leading to elevated CF values. On the other hand, the study by Yamada et al. (2023) for a specific sample of local ultraluminous infrared galaxies (ULIRGs) shows that luminosities of the torus and polar dust components are correlated with each other and with the Eddington ratio. However, the ratio of L_torus/L_agn and L_polar/L_agn are not correlated with the Eddington ratio. Those studies were done for local AGNs. How they scale up to high-z sources and higher luminosity in that matter is an open question. Ricci et al. (2022) noted that, due to a higher amount of ISM or host dust, or obscuring material and mergers, the bulk of the dust-related absorption in high-z sources may take place outside of the SMBH sphere of influence, and as such should depend on L_agn rather than λ_Edd.

As we demonstrated when investigating the difference between W4 and MIPS24 flux measurements, the IR data accuracy is a particularly important factor, which may affect the main conclusions of the analysis. As discussed in detail in Appendix A, the low-luminosity tails of both the low-z and high-z samples diverge from the overall L_agn − L_IR correlation.

In the gathered sample of quasars we do see significant bolometric and infrared luminosity evolution with redshift (in this context see, e.g., Lusso et al. 2013). This evolution manifests as a systematic difference between the low-z and high-z luminosity values. The luminosity evolution is also formally confirmed by the EP test. However, if we limit the Spitzer sample to the higher SMBH masses (over 10^8.5 M_⊙) in both low-z and high-z ranges, then we do not see any CF ≡ L_IR/L_agn evolution with redshift (see Fig. 7), and this observation is confirmed with the Kolmogorov–Smirnov test. The scenario with non-evolving CF means that the dusty torus properties are universal across a wide redshift range, at least in the high SMBH mass regime. This result stands, within the errors, with the selection of quasars with similar L_agn and similar λ_Edd, as described in Appendix E.

The CF distribution was compared with two previous works (Gu 2013; Toba et al. 2021). In Fig. 9 we show the relationship between the CF and L_agn for the low-z subset exclusively, with the regression line fitted to our data, and along with the Toba et al. (2021) regression lines. Toba et al. (2021) used the Bayesian approach, based on the likelihood function from the work by Kelly (2007), which uses error weights. The two different regression lines stand for different models used by Toba et al. (2021) to calculate the luminosities: the torus dust (green line) is the basic model, while the torus and polar dust model (red line) also takes into account the dust placed in polar areas of an AGN.

In Fig. 10 we compare our results regarding the log CF − log L_agn anticorrelation with those of Gu (2013). It should be noted that Gu (2013) used the power-law method of integrating the luminosities. The differences between the estimations were described in Appendix D.1. The overall trend is similar for the low-z quasars, with the slope in both regressions having only a minor difference, namely −0.29 ± 0.01 in our work and −0.19 ± 0.01 in Gu (2013). For the high-z objects, the slopes emerging from our work and Gu (2013) are again similar (−0.46 ± 0.01 and −0.42 ± 0.01, respectively). The difference becomes more pronounced, however, when only the high-zSpitzer subsample is taken into account (slope +0.12). In Fig. 5 the Spitzer subsample is in better agreement with the Gu (2013) correlations, with the low-z slope of 0.82 ± 0.01, and with the high-z slope of 0.74 ± 0.01 (for a more in-depth comparison with the work by Gu 2013, see Appendix D).

The CF value translates to the torus viewing angle in the type I AGN population. Possible effects of the viewing angle evolution are crucial for the quasar-based cosmology (Prince et al. 2021). Our “no CF evolution” result is in agreement with the studies of the torus viewing angle evolution (Prince et al. 2022), although observational biases and systematics induced by the applied methods must also be addressed. This result is, on the other hand, in disagreement with Gu (2013), who argued for an evolution of the CF with redshift. However, as we pointed out in this study, our improvement in deriving L_IR and L_agn leads to systematic differences in luminosity estimates in comparison to the Gu (2013) method, and as a result also in the CF estimates.

To further improve the accuracy of our simple methodology we must address several pending issues. As we have shown, using the W4 data, especially without a restricting quality cut, introduces a strong bias. Spitzer MIPS is a better choice in terms of data accuracy. Altogether, error-weighting is a good practice. We describe these points further in the Appendices of this article. Additionally, the low-luminosity tail of the low-z sample can be affected by an integration limit in L_agn, due to hotter standard accretion disks around lower mass SMBHs, and by an unknown intrinsic extinction. Both effects may cause an underestimation of the bolometric quasar luminosity. Further investigation of the data limitation in the IR and in the optical/UV ranges may also be needed; they will be the subjects of future studies.

7. Conclusions

In the gathered sample of quasars with the available IR–UV data, we see significant luminosity evolutions with redshift, but no significant scaling of the covering factor of the dusty torus (defined as the IR-to-bolometric luminosity ratio) with luminosities. That is, the low-z and high-z sources follow a similar correlation between the IR and bolometric luminosities. We also do not observe any covering factor evolution with redshift, at least in the regime of high SMBH masses, M_BH > 10^8.5 M_⊙, and when using more accurate Spitzer data. On the other hand, the low-luminosity tail of the low-z sample breaks the trend set by the high-mass segment. This is in line with the work by Toba et al. (2021), who proved strong polar dust contamination in the covering factor in the lower luminosity range of the low-z quasar population. Surprisingly, the high-z sample alone also shows a similar kind of dependence, but in a higher-luminosity regime, although here it is caused by limited data quality. To verify whether this effect is physical, we need to further test our methodology and data accuracy, while keeping in mind how the WISE W4 fluxes differ from the Spitzer MIPS 24 μm data. Any potential scaling of the covering factor with luminosity should be robustly identified and calibrated for accurate quasar-based cosmology.

Acknowledgments

We thank the anonymous referee for their remarks. We thank Bożena Czerny for helpful discussions. This research was supported by the Polish National Science Centre grant UMO-2018/30/M/ST9/00757 and by Polish Ministry of Science and Higher Education grant DIR/WK/2018/12. Ł.S. was supported by the Polish NSC grant 2016/22/E/ST9/00061. This research was funded by the Priority Research Area Digiworld under the program Excellence Initiative – Research University at the Jagiellonian University in Kraków.

References

Appendix A: Comparison of WISE W4 and MIPS M24 datasets

A.1. Data quality comparison

Here we address a calibration issue known as the “redleak”, using the correction method described by Brown et al. (2014). This problem arises from the dependence of the filter correction factors on the slope of an SED. While most of the WISE calibration was conducted with the use of stars, which appear blue in the IR band, galaxies due to the cold dust emission exhibit different SED slopes and appear red. Consequently, applying a standard correction to extragalactic sources may underestimate W3 and overestimate W4 by about 9% (Wright et al. 2010). However, in our case, accounting for the redleak correction had a minimal impact.

To validate the accuracy of the WISE W4 observations, the Spitzer M24 data were examined. We conducted a cross-match using data from the Spitzer Science Center (SSC) & Infrared Science Archive (IRSA) (2021) with a cross-match radius of 3 arcsec. The data were available for only 290 objects in the Low-z sample, and 157 objects in the High-z sample. The Spitzer data are shown in Figures C.1 and C.2. The relation L_IR versus L_agn shown in Figure C.1 top panel is characterized by a much reduced dispersion, but follows a similar roughly linear trend. Figure C.1 middle panel shows a slightly different relation for both redshift samples, with the Spitzer data generally having a lower variance in log CF. Notably, the log CF values of Low-z and High-z sources are much closer to the M24-based L_IR, as can be seen in Table 2. The change in the amount of Spitzer data points between Figure C.1 and Figure C.2 is due to the S/N_W3 > 3 criterion. The resultant datasets from Spitzer include 239 Low-z sources and 122 High-z targets. In our analysis we carefully examined the difference in luminosity estimates based on the W4 and M24. Even with S/N_{W3& W4} > 5 there is still some discrepancy between each monochromatic luminosity; W4 tends to overestimate the luminosity values compared to M24 (see Figure A.2, lower panel).

Figure A.2 top panel shows the comparison between the monochromatic luminosities of W4 and M24, both calculated with the all-points method. Notably, there is a clear deviation from the 1:1 line at lower luminosities, in both redshift ranges. This discrepancy between W4 and M24 diminishes as luminosity increases. The major deviation from the 1:1 line signals significant issues with the photometry measurements in one of the bands. This issue remains even in the case of the best-quality W4 measurements with the S/N_W4 > 3, as shown in the bottom panel of Figure A.2. We conducted an additional test on the quality of both W4 and M24 by directly comparing photometry observations from maps with the determination of measured flux, and comparing it with the values from the archive. In general, W4 tends to extract fluxes from larger areas on the sky when compared with the higher angular resolution of M24 (2.55″ vs 12″).

The Spitzer and WISE instruments differ in spatial resolution and point spread functions (PSFs). As we show in Figure A.3, in some cases W4 blends neighboring sources. Additionally, the wide-field character of WISE increases the probability of a light pollution caused by a bright star in the field of view, or nonsidereal objects such as asteroids. Lower-luminosity sources are affected in this respect to a larger extent, as shown directly in Figs. A.1 and A.2.

Fig. A.1. Comparison of data quality for the M24 and W4 with different restrictions on S/NW4. The panels present (from top to bottom) log LW4 − log LM24 vs log LM24 with S/NW4 > 0, > 3, and > 5. The orange circles indicate the Low-z quasars; the red diamonds are the High-z sources.

Fig. A.2. Comparison between monochromatic WISE W4 and Spitzer M24 luminosities for Low-z and for High-z sources. The black line in the panel shows the 1:1 relation. The orange circles and red diamonds indicate Low-z and High-z quasars with Spitzer M24 data. The panels (from top to bottom) correspond to S/NW4 > 0, > 3, and > 5.

Fig. A.3. Sample raw images of the Spitzer MIPS 24 μm and WISE W4 frames. The source position is shown as a magenta point at the center of each frame. The colors represent flux levels in each frame: violet–dark blue for low flux, and red for the highest flux level.

A.2. Spectral energy distribution

To quantify the differences between WISE and Spitzer photometry even further, we compared the mean SEDs for Low-z, High-z (S/N_{W3& W4} > 5) and Low-zSpitzer, High-zSpitzer (S/N_W3 > 5). The comparison is shown in Fig. A.4. As described in the data selection in Sect. 2, the main difference between data with and without Spitzer is the replacement of the WISE W4 filter with MIPS 24 μm. It is worth noting that for both redshift bins Spitzer-based datasets have lower L_IR and L_agn luminosities (approximately 0.1 dex luminosity difference for Low-z and 0.2 dex for High-z samples). The difference in luminosity is even larger between W4 and M24 filters and is visible for both redshift bins. Generally, despite the systematic differences in luminosity, the average SEDs for W4 and M24 data exhibit very similar shapes. Similar SED shapes support the idea that average CF should be comparable and the apparent difference should be caused by the systematic difference between WISE and Spitzer photometry. It is worth noting that the Spitzer sample mostly contains sources of lower value of L_IR distribution of the WISE sample (with S/N_W3 > 5).

Fig. A.4. Mean rest-frame SEDs for the Low-z and High-z quasars (top and bottom panels, respectively). The red points represent data obtained solely from WISE IR (Low-z top panel, and High-z bottom panel) and the blue points from Spitzer M24 (Low-zSpitzer top panel and High-zSpitzer bottom panel). The Low-z and High-z data shown have S/NW3& W4 > 5, while the Spitzer datasets have S/NW3 > 5. The black and green shading indicate the areas used for the calculations of LIR and Lagn, respectively for both Low-z and High-z. The red and yellow shading indicate the areas used for the calculations of LIR and Lagn, respectively for both Low-zSpitzer and High-zSpitzer. Each photometric observation is labeled with the filter name. The error bars are defined as the 1st and 3rd quartiles.

Appendix B: Different likelihood for Bayesian fitting

One of the tested likelihoods is similar to the main one used in the analysis, but without the error weights:

$$\begin{array}{c}\hfill P(D|\mathrm{\Theta },M)=-{\displaystyle \sum _n[{(y_n-m{x}_n-b)}^2\times ln\left({(y_n-m{x}_n-b)}^2\right)].}\end{array}$$(B.1)

The fit is shown in Figure B.1. The black lines represent Bayesian fitting with the likelihood function B.1. The high spread is visible in both panels, a) and b). Although the fitted lines for Low-z sources are similar to weighted likelihood, for the High-z targets the difference is significant (slopes 0.12 and −0.50 for likelihood with and without error weighted accordingly). The reason for this discrepancy is most likely the errors of low-luminosity objects which are present in both Low-z and High-z samples, but affect mostly the High-z fit.

Fig. B.1. Comparison of two different likelihoods in Bayesian fitting methods for Spitzer data between log CF vs log Lagn for a) Low-z sources and b) High-z sources. The orange and red lines represent the same Bayesian fitting as in Figure 10, while the black lines represent the fitting based on the likelihood without error weights.

Appendix C: Luminosity scalings

C.1. The entire WISE sample

Figure C.1 illustrates the relations between the CF, L_IR, and L_agn, for the complete WISE sample. Additional Spitzer data, which are also shown in the figures, are discussed further in Appendix A. Figure C.1 top panel shows, in particular, the relation between log L_IR and log L_agn. The Pearson correlation coefficient values for this relationship are 0.9 for the Low-z sample and 0.70 for the High-z sample, as determined by the best-fitting least-squares linear regression. Notably, both the Low-z and High-z quasars exhibit similar trends.

Fig. C.1. Luminosity scalings within the complete WISE datasets for Low-z and High-z quasars. The top panel shows log LIR vs log Lagn, the middle panel illustrates log CF vs log Lagn, and the bottom panel displays log CF vs log LIR. The blue circles indicate Low-z quasars, whereas the green and red diamonds indicate High-z sources. The yellow circles and red diamonds indicate Low-z and High-z quasars with the Spitzer M24 data.

Figure C.1 middle panel shows the relation between log CF and log L_agn. Both the Low-z and High-z quasars exhibit similar scaling, but the Low-z quasars have a greater spread in L_agn. On the other hand, the High-z quasars have higher values of log CF, and the general scaling relation is more robust. An anticorrelation between log CF and log L_agn is present, albeit with a high variance. However, no such anticorrelation is evident in the bottom panel of Figure C.1, which shows the relation between log CF and log L_IR.

C.2. S/N > 3

The Low-z sample contains 3,229 objects, while the High-z sample includes 1,279 objects, all of which meet the S/N_{W3& W4} > 3 criterion. Figure C.2 shows the relations between the L_IR, L_agn, and CF parameters in direct analogy to those given in Figure C.1, but with the S/N_{W3& W4} > 3 cut. One notable difference is the absence of low-luminosity objects in the top panel of the figure, when compared to Figure C.1 top panel.

Fig. C.2. Luminosity scalings for the WISE samples of the Low-z and High-z quasars with S/NW3& W4 > 3. The top panel shows log LIR vs log Lagn, the middle panel log CF vs log Lagn, and the bottom panel log CF vs log LIR. The blue circles indicate Low-z quasars, whereas the green and red diamonds indicate High-z sources. The yellow circles and red diamonds indicate Low-z and High-z quasars with the Spitzer M24 data. In the top panel the blue line denotes the Bayesian regression for the data with S/NW3& W4 > 3, the orange line is the Bayesian regression for the Spitzer data from both Low-z and High-z subsamples, and the black line is the 1:1 relation between the LIR and Lagn.

Appendix D: Comparison with Gu 2013

The power-law luminosity estimation method proposed by Gu (2013) involved fitting two power-law functions to the IR and optical–UV segments of the SED for each source. Here we followed this method for a direct comparison with our all-points method, and therefore interpolating luminosities at specific wavelengths of 0.01 μm (L_0.11 μm), 1 μm (L_1 μm), and 7 μm (L_7 μm) from the same filters. The calculation of L_IR involved integrating the area under the power-law fit between L_7 μm and L_1 μm, while L_agn was calculated as the integral of the fitted power law between the luminosities L_1 μm and L_0.11 μm.

D.1. Comparison of luminosity estimates

In this section we compare the results obtained using the all-points method for luminosity estimation with those obtained using the power-law method. As previously described, in the power-law method both L_IR and L_agn are calculated as the areas under the power-law fitted curves between the two points of the interpolated monochromatic luminosities. The power-law method (red and blue areas in the figure) agrees relatively well with the photometry for Low-z quasars, but deviates from the monochromatic luminosities for High-z quasars. A notable limitation of the power-law method is the underestimation of L_agn and the overestimation of L_IR for High-z quasars. Consequently, this method may lead to a significant overestimation of the CF.

Figure D.1 presents the comparison of L_IR versus L_agn calculated with each method (upper and lower panels for the all-points and power-law methods, respectively). As shown, the luminosities of the High-z quasars following from the power-law method have a greater dispersion compared to the all-points method in the log L_IR versus log L_agn representation. However, the number of outliers in the Low-z sample increases with the all-points method. The High-z quasars feature a major difference in the log L_IR versus log L_agn relation, having a slight tail at lower values of log L_agn.

Fig. D.1. Relation between log Lagn and log LIR calculated with the all-points method and the power-law method (upper and lower panels, respectively), for the Low-z and High-z quasars (blue circles and green diamonds, respectively).

Figure D.2 shows the direct comparison between the all-points and power-law methods for L_IR (upper panel) and L_agn (lower panel). For the majority of the Low-z objects, the difference becomes significant above 45.5 for both log L_IR and log L_agn, with the power-law method underestimating the values compared to the all-points method. For the High-z sample, the power-law method tends to slightly overestimate IR luminosities up to log L_IR = 46.5, above which the trend reverses. In comparison, the L_agn for High-z quasars is overall significantly underestimated with the power-law method.

Fig. D.2. Comparison between luminosities based on the power-law (Y-axis) and all-points (X-axis) methods of integration. The upper panel shows the comparison between LIR, and the lower panel panel between Lagn. The black line stands for the 1:1 relation. In both panels, the blue circles and green diamonds denote the Low-z and High-z quasars, respectively. The red and orange lines represent the best Bayesian regression fit, with the best-fit relation given in each panel.

D.2. Recreation of Gu 2013

To recreate the results of the Gu (2013) work, the same data were analyzed. While there may be minor discrepancies in the data due to potential updates in the database, we aimed for a direct comparison, and therefore we utilized the same SDSS data releases as in Gu (2013), and the two exact redshift ranges were defined: the Low-z range within 0.7 ≤ z ≤ 1.1 and the High-z range within 2.0 ≤ z ≤ 2.4. Optical data for the Low-z objects, SDSS DR7Q, and DR9Q, were taken from Schneider et al. (2010) and Pâris et al. (2012), and for the High-z quasars only the SDSS DR9Q was used. The DR7Q contains 16,695 Low-z quasars, while the DR9Q contains 10,553 Low-z and 25,934 High-z quasars. The SDSS data were then cross-matched within a 2 arcsec radius with the WISE All-Sky Data Release. This returned 11,943 Low-z quasars from DR7Q, and 4,149 Low-z and 16,228 High-z quasars from DR9Q. The following data cross-matching was done with the UKIDSS DR8 (also within 2 arcsec radius), with 4,855 objects in DR7Q Low-z, and 2,517 in Low-z and 2,630 in High-z in DR9Q. Finally, only the Low-z sample was cross-matched with the GALEX GR6/GR7 data release with a 3 arcsec radius, resulting in 3,136 DR7Q and 1,134 DR9Q quasars. The final Low-z sample contains 4,218 quasars (after removal of duplicates), while the final High-z sample contains 2,630 quasars. In comparison, Gu (2013) claimed 3,940 and 2,056 objects in the Low-z and High-z samples, respectively.

Despite the differences in sample size, our results for the integrated luminosities and the CF are consistent within the error margins. The median integrated luminosities for the Low-z sample are log L_IR = 45.91 ± 0.44, log L_agn = 46.08 ± 0.42, while for the High-z sample they are log L_IR = 46.13 ± 0.11, and log L_agn = 46.01 ± 0.17. The median values of the CF for the Low-z and High-z data samples are CF_{low − z} = 0.70 ± 0.11 and CF_{high − z} = 1.23 ± 0.37. These error estimates are the median absolute deviations.

The OLS regressions were fitted to the data, with the following results:

$$\begin{array}{c}\hfill log{L}_{\mathrm{IR}}=(0.96\pm 0.01)log{L}_{\mathrm{agn}}+(1.48\pm 0.12)\end{array}$$

for the High-z quasars, and

$$\begin{array}{c}\hfill log{L}_{\mathrm{IR}}=(0.42\pm 0.01)log{L}_{\mathrm{agn}}+(22.84\pm 0.33)\end{array}$$

for the Low-z quasars;

$$\begin{array}{c}\hfill log{L}_{\mathrm{IR}}/L_{\mathrm{agn}}=(-0.04\pm 0.01)log{L}_{\mathrm{agn}}+(1.48\pm 0.12)\end{array}$$

for the High-z sample, and

$$\begin{array}{c}\hfill log{L}_{\mathrm{IR}}/L_{\mathrm{agn}}=(-0.58\pm 0.008)log{L}_{\mathrm{agn}}+(22.84\pm 0.33)\end{array}$$

for the Low-z sample.

All in all, the exact recreation of the data sample from Gu (2013) turned out to be quite problematic. We carefully checked, for example, for one of the quasars for which the SEDs were shown (SDSS J001600.60-003859.2), that the current version of the SDSS DR12Q catalog gives negative flux measurements, so it is unclear how the target was selected in Gu (2013).

Appendix E: Other physical properties cuts

In this article we primarily focus on selecting objects with similar SMBH masses, as shown in Figure 7. However, this approach has a limitation because objects with similar M_BH can exhibit very different luminosities, depending on the accretion rate. To double-check the different possible cuts, we performed a similar luminosity selection and a similar Eddington ratio cut. We also checked how the M_SMBH cut selected the Low-zSpitzer sources in L_agn and L_IR. In Figure E.1 we present the objects before and after the cut at log M_BH/M_⊙ > 10.5. As shown, the M_BH cut excluded approximately half of the Low-zSpitzer objects with luminosities below log L_IR & agn < 46.0.

Fig. E.1. LIR vs Lagn relation for Spitzer data. The left panel shows the objects without the log MBH = 8.5 cut; the right panel shows the objects above the log MBH = 8.5 threshold.

The L_agn luminosity cut was also performed with the selection of both Low-zSpitzer and High-zSpitzer at log L_agn > 46.5. The lack of consistently wide coverage of the L_IR − L_agn space in both datasets is problematic for this method, which is shown in Figure E.2a) (left panel). Consequently, the samples after the cut were unevenly sized, and the CF values differed, though they were comparable within the margins of error (Figure E.2b). In the Low-z sample the most luminous and the least obscured objects are present after the cut. For comparison, we also included the M_BH cut, which showed consistent results with the luminosity cut within the errors.

Fig. E.2. Spitzer data after the S/NW3 > 5, but without the cross-match with VAC. We performed the cut in log Lagn > 46.5 for both Low-z and High-z. The right panel presents the CF distribution for the truncated sample. The black and blue points represents medians for the MBH cut with MAD errors. The green and red points represent the medians and MAD errors for newly cut data.

Finally, the Eddington ratio cut was performed. The resulting scaling between the Eddington ratio and L_agn is shown in Figure E.3a (left panel). As one can see, the spread in the Eddington ratio is quite large, which can influence this method. We decided to apply a cut at an Eddington ratio value of −1.1. Figure E.3b shows datasets with a value above −1.1 for both Low-zSpitzer (blue and orange squares) and High-zSpitzer (green and red squares). The distribution of CF is comparable within the errors, with the Eddington ratio cut.

Fig. E.3. Spitzer data after the S/NW3 > 5, but without the cross-match with VAC. The cut was performed in log LIR/Lagn = −1.1 and two subsamples were created, for Low-z and High-z. The right panel presents the CF distribution for the truncated sample. The black and pink points represent medians for the MBH cut with MAD errors. The orange and red points represent the medians and MAD errors for the low and high sample with the Eddington ratio exceeding −1.1, while the blue and green points represents the samples with the Eddington ratio below −1.1. The samples are comparable within the errors.

All Tables

All Figures

Fig. 1. Flowchart showing the data selection process, as described in Sect. 2. The final samples are 1534 objects for low-z with S/NW3 & W4 > 5 and 309 objects for high-zS/NW3 & W4 > 5. The low-z and high-z dataset combined create the low-high-z. High-precision datasets with Spitzer observations have 198 objects for low-zSpitzer and 91 objects for high-zSpitzer, both with S/NW3 > 5. The datasets with different S/NW3 & W4 cuts are analyzed in Appendix C.

In the text

Fig. 2. Mean rest-frame SEDs for low-z quasars (top) and high-z quasars (bottom). The top and bottom panels depict the mean rest-frame SEDs for low-z and high-z quasars, respectively. The black and green shaded regions highlight the areas utilized for computing LIR and Lagn, respectively, in the all-points method. Each photometric observation is labeled with the filter name. The error bars denote the 1st and 3rd quartiles.

In the text

Fig. 3. Relation between log LIR and log Lagn with S/NW3 & W4 > 5. The blue and orange circles indicate low-z quasars, whereas the green and red diamonds indicate high-z sources. The orange circles and red diamonds indicate low-z and high-z quasars with Spitzer M24 data. In panel a the blue line denotes the best-fit Bayesian regression for the low-high-z data with S/NW3 & W4 > 5, the orange line is the Bayesian regression for the low-high-zSpitzer data from both low-z and high-z, and the green and red dashed lines represent the best Bayesian fit with the low-z sample weighted and low-z with Spitzer weighted, respectively; the black line gives the 1:1 scaling relation between LIR and Lagn. In panel b the big black circles represent the medians calculated for the 0.5 dex in log Lagn for both redshift samples (low-high-z) and the orange big triangles stand for the medians calculated for both redshift samples with the Spitzer data (low-high-zSpitzer). For the black and the magenta points, the error bars represent MAD. The exact values of the fitted parameters are listed in Table 1.

In the text

	Fig. 4. Comparison between the Bayesian regression analysis model parameters: (a) slopes, (b) intercepts, obtained for different analyzed samples of quasars (see Table 1 and Sect. 5.1).
In the text

	Fig. 5. Regression analysis results for log Lagn and log LIR for low-z, low-zSpitzer, and high-z, high-zSpitzer sources separately. The model parameters of the Bayesian linear regression analysis are given in Table 1.
In the text

Fig. 6. Relationships between CF and LIR (upper panel a) and Lagn (lower panel b), including only the sources with S/NW3 & W4 > 5. The blue circles indicate low-z quasars, and the green and red diamonds indicate high-z sources. The orange circles and red diamonds indicate low-z and high-z quasars with the Spitzer M24 data. Additionally, the median statistics were calculated for the samples, as described in the legend in the lower left corner. The error bars represent the MAD errors.

In the text

Fig. 7. Relations of the CF with redshift. The top panel shows the relation between log CF and log MBH for quasars with the Spitzer M24 and VAC data. The red line represents the cut on log MBH/M⊙ = 8.5. The bottom panel shows the relation between log CF and redshift for quasars with log MBH/M⊙ > 8.5. The orange and red points in both panels represent the low-z and high-z quasars, respectively. The large symbols with MAD error bars stand for the medians calculated for each data sample: low-z with no cut in MBH (red triangle), high-z (orange triangle), low-z with the cut in log MBH/M⊙ > 8.2 (pink triangle), low-zSpitzer with no cut in MBH (black square), high-zSpitzer (green square), low-zSpitzer with the cut in log MBH/M⊙ > 8.5 (blue square). For the luminosity and Eddington ratio selection, see Appendix E.

In the text

Fig. 8. Analysis of the intrinsic scatter. Top panel: scatter plot for the extinction corrected log fr and log fW2, for low-z, high-z, low-zSpitzer, and high-zSpitzer. Bottom panel: standard deviation for binned fluxes in the bands r and W2, and the error of log r − W2 relations σr W2 scatter. The bins were calculated as equally distanced in log r space. The sizes of the symbols correspond to the log fr flux values in each bin.

In the text

	Fig. 9. Regression analysis results for the covering factor and bolometric luminosity for the low-z quasars with S/NW3 & W4 > 5. The blue and orange lines are the same as in Fig. 5. The green line corresponds to the regression from Toba et al. (2021) fitted to the model with torus dust only, while the red line is the sum of regressions fitted to the torus and polar dust models separately.
In the text

	Fig. 10. Relation between log CF and log Lagn for the low-z and high-z sources with S/NW3 & W4 > 5 (upper and lower panels, respectively). The regression lines in both panels fitted to our various subsamples, with and without the Spitzer data, along with the Gu (2013) regression lines for comparison, are described in the legends of the two panels.
In the text

	Fig. A.1. Comparison of data quality for the M24 and W4 with different restrictions on S/NW4. The panels present (from top to bottom) log LW4 − log LM24 vs log LM24 with S/NW4 > 0, > 3, and > 5. The orange circles indicate the Low-z quasars; the red diamonds are the High-z sources.
In the text

	Fig. A.2. Comparison between monochromatic WISE W4 and Spitzer M24 luminosities for Low-z and for High-z sources. The black line in the panel shows the 1:1 relation. The orange circles and red diamonds indicate Low-z and High-z quasars with Spitzer M24 data. The panels (from top to bottom) correspond to S/NW4 > 0, > 3, and > 5.
In the text

	Fig. A.3. Sample raw images of the Spitzer MIPS 24 μm and WISE W4 frames. The source position is shown as a magenta point at the center of each frame. The colors represent flux levels in each frame: violet–dark blue for low flux, and red for the highest flux level.
In the text

Fig. A.4. Mean rest-frame SEDs for the Low-z and High-z quasars (top and bottom panels, respectively). The red points represent data obtained solely from WISE IR (Low-z top panel, and High-z bottom panel) and the blue points from Spitzer M24 (Low-zSpitzer top panel and High-zSpitzer bottom panel). The Low-z and High-z data shown have S/NW3& W4 > 5, while the Spitzer datasets have S/NW3 > 5. The black and green shading indicate the areas used for the calculations of LIR and Lagn, respectively for both Low-z and High-z. The red and yellow shading indicate the areas used for the calculations of LIR and Lagn, respectively for both Low-zSpitzer and High-zSpitzer. Each photometric observation is labeled with the filter name. The error bars are defined as the 1st and 3rd quartiles.

In the text

	Fig. B.1. Comparison of two different likelihoods in Bayesian fitting methods for Spitzer data between log CF vs log Lagn for a) Low-z sources and b) High-z sources. The orange and red lines represent the same Bayesian fitting as in Figure 10, while the black lines represent the fitting based on the likelihood without error weights.
In the text

Fig. C.1. Luminosity scalings within the complete WISE datasets for Low-z and High-z quasars. The top panel shows log LIR vs log Lagn, the middle panel illustrates log CF vs log Lagn, and the bottom panel displays log CF vs log LIR. The blue circles indicate Low-z quasars, whereas the green and red diamonds indicate High-z sources. The yellow circles and red diamonds indicate Low-z and High-z quasars with the Spitzer M24 data.

In the text

Fig. C.2. Luminosity scalings for the WISE samples of the Low-z and High-z quasars with S/NW3& W4 > 3. The top panel shows log LIR vs log Lagn, the middle panel log CF vs log Lagn, and the bottom panel log CF vs log LIR. The blue circles indicate Low-z quasars, whereas the green and red diamonds indicate High-z sources. The yellow circles and red diamonds indicate Low-z and High-z quasars with the Spitzer M24 data. In the top panel the blue line denotes the Bayesian regression for the data with S/NW3& W4 > 3, the orange line is the Bayesian regression for the Spitzer data from both Low-z and High-z subsamples, and the black line is the 1:1 relation between the LIR and Lagn.

In the text

	Fig. D.1. Relation between log Lagn and log LIR calculated with the all-points method and the power-law method (upper and lower panels, respectively), for the Low-z and High-z quasars (blue circles and green diamonds, respectively).
In the text

Fig. D.2. Comparison between luminosities based on the power-law (Y-axis) and all-points (X-axis) methods of integration. The upper panel shows the comparison between LIR, and the lower panel panel between Lagn. The black line stands for the 1:1 relation. In both panels, the blue circles and green diamonds denote the Low-z and High-z quasars, respectively. The red and orange lines represent the best Bayesian regression fit, with the best-fit relation given in each panel.

In the text

	Fig. E.1. LIR vs Lagn relation for Spitzer data. The left panel shows the objects without the log MBH = 8.5 cut; the right panel shows the objects above the log MBH = 8.5 threshold.
In the text

	Fig. E.2. Spitzer data after the S/NW3 > 5, but without the cross-match with VAC. We performed the cut in log Lagn > 46.5 for both Low-z and High-z. The right panel presents the CF distribution for the truncated sample. The black and blue points represents medians for the MBH cut with MAD errors. The green and red points represent the medians and MAD errors for newly cut data.
In the text

Fig. E.3. Spitzer data after the S/NW3 > 5, but without the cross-match with VAC. The cut was performed in log LIR/Lagn = −1.1 and two subsamples were created, for Low-z and High-z. The right panel presents the CF distribution for the truncated sample. The black and pink points represent medians for the MBH cut with MAD errors. The orange and red points represent the medians and MAD errors for the low and high sample with the Eddington ratio exceeding −1.1, while the blue and green points represents the samples with the Eddington ratio below −1.1. The samples are comparable within the errors.

In the text