© The Authors 2026

1. Introduction

The unification schemes for active galactic nuclei (AGNs; Antonucci & Miller 1985; Urry & Padovani 1995; Netzer 2015; Hickox & Alexander 2018) were proposed to provide a simplified explanation of the structure of the AGN mechanism. Unification based on the inclination angle, where the dusty torus hides some of the internal structures, was questioned with superior observations with Gravity GRAVITY Collaboration (2020). Recent studies suggest that the structure of the torus might be more complex (consisting of clumps and clouds) Nenkova et al. (2008) and that it can be directly linked to the broad-line regions (BLRs) and the accretion disk (Czerny & Hryniewicz 2011; Czerny et al. 2019). There have also been indications that there is dust located in the polar regions of AGNs, so-called polar dust (PD; Toba et al. 2021; Gámez Rosas et al. 2022; Davies et al. 2022). The dusty torus can change the obscuration along the line of sight, its spectral properties, and the classification of an AGN. Additionally, some authors have found a significant variance in the amount of dust in AGNs (e.g., Roseboom et al. 2013).

Among AGN, quasars are particularly important due to their extreme brightness and wide redshift distribution. In recent years, significant progress has been made in calibrating quasars as standard candles for cosmological tests. The search for new standard candles is crucial, especially considering the pressing issue in modern observational cosmology known as the “Hubble tension” (Planck Collaboration XVI 2014; Verde et al. 2013; Planck Collaboration VI 2020; Verde et al. 2019; Riess et al. 2022). Recent studies have suggested that quasars observed with JWST at a high redshift (z > 2) may exhibit different properties or physical structures compared to local quasars (Maiolino et al. 2024).

Previous approaches utilized quasars for cosmology by making use of two major methodologies. The first methodology is based on the technique known as reverberation mapping (Kaspi et al. 2000; Bentz et al. 2013; Czerny & Hryniewicz 2011; Czerny et al. 2013; Zajaček et al. 2020; Czerny et al. 2019, 2021). This technique involves calculating the time lag between emissions from different regions within the AGN – specifically, between the accretion disk surrounding the supermassive black hole at the quasar’s center and the BLR around it. From the time lag, the radius of more distant regions can be calculated, along with the absolute luminosity. Using the absolute luminosity, the luminosity distance can be determined, and by comparing it to the redshift, cosmological models can be tested. The primary drawback of this method is its time consumption, as measuring the time lag requires years of observations for each object. Additionally the high dispersion of the R(BLR)-L method makes it problematic to perform reverberation mapping outside of the local Universe.

An alternative approach is based on the nonlinear luminosity relation proposed in the Risaliti-Lusso relation (Risaliti & Lusso 2015, 2019; Signorini et al. 2023). The Risaliti-Lusso relation is based on a widely accepted scheme that the ultraviolet emission of AGN is produced by an accretion disk, while the X-rays are produced by a hot corona around the accretion disk. Studies (Tananbaum et al. 1979; Zamorani et al. 1981; Avni & Tananbaum 1986) have shown that the relationship between X-ray and UV luminosities can be parametrized as log(L_X) = γlog(L_UV)+β, where L_X and L_UV are the rest-frame monochromatic luminosities at 2 keV and 2500 Å, respectively. Although L_X and L_UV are measured independently, they are physically correlated with each other and with the luminosity distance.

The optical (OPT) and infrared (IR) luminosity relation for quasars also shows promise. Several studies have indicated that this relation can be robust (Gu 2013; Rałowski et al. 2024; Trefoloni et al. 2024). The reason is that IR luminosity (L_IR) and optical luminosity (L_OPT) should be strictly connected. In bright quasars, the bulk of L_OPT emission can be explained by the emission from the accretion disk Netzer (2015). Part of this emission is absorbed and reemited by the dusty torus surrounding the supermassive black hole. Thus L_OPT and L_IR should be tightly correlated. The degree of this correlation is still not clear, but if it were possible to account for (1) potential selection effects and (2) any physical contaminations, it may then be possible to use it as a universal relation that is constant with redshift. In a previous study, Rałowski et al. (2024) we showed that the OPT–IR relation may be influenced by selection effects in the data, particularly involving the W3 and W4 filters in the IR, and that the covering factor, calculated from the L_IR/L_AGN luminosity ratio, does not evolve with redshift for quasars with similar physical properties. Our other conclusion was that the OPT-IR relation has to be calibrated, as it is influence by contamination, such as from polar dust.

The importance of such a calibration is even more significant because we believe that the OPT-IR nonlinear luminosity relation can be used for cosmological tests. In essence, it is similar to the X-UV Risaliti-Lusso relation that is widely used and debated today. For the OPT-IR photometrical relation to be useful for cosmology, it is necessary to exclude or/and understand the (1) systematics and the (2) additional contributions to the OPT or IR emission from components of the quasar host galaxy that are not connected to AGN emission. In the following, we analyze the photometrical observations, especially through the spectral energy distribution (SED) fitting, to find the most important factors and contaminations that can change the OPT-IR relation connected only to the central source AGN.

To calibrate the luminosity relation for objects as complex as quasars, dependencies on physical parameters must be carefully analyzed. Some works (e.g., Toba et al. 2021; Buat et al. 2021) have indicated that the L_IR − L_opt relation can depend on physical properties of AGN such as polar dust. It is also possible that the stellar light of the host galaxy influences the L_OPT. In Rałowski et al. (2024) we discuss possible explanations of how the low-luminosity part of the low-z sample breaks the trend set by the high-mass part in L_IR/L_OPT. This is in line with the work by Toba et al. (2021), who suggested strong polar dust contamination in L_IR in the lower-luminosity part of the AGN population. Reddening and dust obscuration at small redshift was reported to affect type-II AGNs more than type-I AGNs (Netzer 2015). The solution to this may be selecting the samples based on the IR, but this approach would be biased toward sources with a larger covering factor and against sources with large inclination angles (see Sect. 2 in Netzer 2015). Another important contamination may be emission from star-forming and inter regions in the galaxy as well as interstellar medium (ISM). The ISM also has an important effect on high redshift quasars (Decarli et al. 2023) and line emission (Kade et al. 2024).

With the work presented here, we aim to analyze the impact of such contaminations while focusing on two key aspects: (1) the role of contaminations in luminosity estimation and how they affect the OPT–IR relation, and (2) whether it is possible to calibrate our method after accounting for these contaminations. The following work is structured as follows: Section 2 presents the data selection process, in Sect. 3 we introduce the methodology we used, and Sect. 4 presents the most important results.

2. Data

The primary dataset used in this analysis is based on photometric observations of quasars. Specifically, the Sloan Digital Sky Survey Quasar Catalog, Sixteenth Data Release (SDSS DR16Q¹; Lyke et al. 2020), was employed. This catalog contains approximately 750 000 quasars, including around 225 000 newly identified sources relative to previous data releases. Quasars in SDSS DR16Q were classified using spectral fitting applied to 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 defilement levels ranging from 0.3 to 1.3% (for further details, see Lyke et al. 2020). The SDSS photometric data are provided in two formats: “asinh magnitudes” and “nanomaggies”² (Lupton et al. 1999). The fluxes in nanomaggies can be converted to AB magnitudes after applying corrections for filter-specific shifts in the u and z filters³.

There are 557 000 objects after the redshift selection within the range of z = 0.7 − 2.4. SDSS DR16Q was subsequently crossmatched using a 3 arcsec radius. After adding the IR from the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010) All-Sky Survey, the sample had 552 372 objects. Near-IR photometry was obtained from the UKIRT Infrared Deep Sky Survey (UKIDSS, Lawrence et al. 2007) resulting in just over 100 000 objects in the sample. For the ultraviolet wavelength range, we used the photometry from the Galaxy Evolution Explorer (GALEX GR6/GR7 data release Martin et al. 2005), which we restricted to low-redshift quasars, as the Lyman-α break becomes significant at higher redshifts. The UV observations were available for 17 463 objects with redshift z ≤ 1.1. Far-IR (FIR) data were gathered from the Herschel mission, specifically from the Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010) and Spectral and Photometric Imaging Receiver (SPIRE; Griffin et al. 2010) instruments. For the Hershel data, we used the IRSA database, with the main condition being for the quasars to have SPIRE 250 μm. The inclusion of Herschel data shrank the number of the whole sample massively, leaving us with only 351 objects. Additionally, SPITZER (Spitzer Science Center & Infrared Science Archive 2021) data were available for a few objects. In the procedure, we also included upper limits, which were accounted for in the SED fitting procedure (described in Sect. 3.1). SDSS, UKIDSS, WISE, and 250 μm SPIRE all have 100% coverage for all the selected objects. The data binning with respect to the redshift is presented below. For each dataset in brackets, the fraction of coverage in the remaining surveys is presented.

We split SDSS DR16Q into four samples based on the redshift ranges (see Fig. 1):

1. low-z: z ∈ [0.7, 1.1], containing 100 objects (PACS: blue 0%, green 2%, red 4%; SPIRE: 350 μm 30%, 500 μm 6%)

2. inter-z1: z ∈ (1.1, 1.5), containing 73 objects (PACS: blue 0%, green 1%, red 0%; SPIRE: 350 μm 22%, 500 μm 1%)

3. inter-z2: z ∈ [1.5, 2.0), containing 94 objects (PACS: blue 0%, green 5%, red 3%; SPIRE: 350 μm 34%, 500 μm 4%)

4. high-z: z ∈ [2.0, 2.4], containing 84 objects (PACS: blue 1%, green 0%, red 1%; SPIRE: 350 μm 37%, 500 μm 6%).

Fig. 1. Flowchart showing the data selection process as described in Sect. 2. The χ2 values were calculated for the best-fit SED in CIGALE. The Ldisk values are from the best-fit estimation of the SKIRTOR model, as described in Sect. 3.2.

The observations were corrected for extinction using the Cardelli extinction law (Cardelli et al. 1989) and dust maps from the work of Schlegel et al. (1998).

3. Methods

3.1. SED fitting

To properly assess the physical properties of a galaxy, it is important to carefully analyze its emitted spectrum. The SED fitting technique provides insight into the broad energy distribution of photons across the electromagnetic spectrum. Using broad-band photometry that spans a wide wavelength range combined with redshift information, we constructed an observed SED. Subsequently, synthetic SEDs were generated from a set of models based on parameter-grid templates, allowing for the estimation of the underlying physical properties.

The SED fitting was conducted using the Code Investigating GALaxy Emission (CIGALE; Boquien et al. 2019), version 2022.1⁴. The CIGALE methodology is based on the assumption of an energy balance principle, and thus the energy emitted by dust in the mid-IR (MIR) to FIR was estimated from the energy absorbed by dust in the UV-optical range. The dust emission can be even larger when taking the presence of the AGN and its dusty torus into account. This approach combines computational speed and a wide range of possible applications to galaxies (Buat et al. 2011, 2018; Boquien et al. 2013; Lo Faro et al. 2017; Salim et al. 2018; Małek et al. 2018), AGNs (Ciesla et al. 2015; Toba et al. 2021), ultraluminous IR galaxies (Małek et al. 2017; Dey et al. 2024). An important aspect of CIGALE is the use of the Bayesian-like approach. In addition, the physical properties are estimated by weighing all the models on the basis of their goodness of fit, with the best-fit models having the heaviest weight. (For more details we refer to Noll et al. 2009; Boquien et al. 2019).

The CIGALE modules include star formation history (SFH), stellar population synthesis, a dust attenuation model, cold dust emission, nebular emission, AGN emission, and X-ray and radio emission (these last two modules were not used). The SFH module calculates the evolution of different stellar populations based on different SFH models, taking into account different ages of stars and the τ-folding time. From these calculations, using the library of stellar population synthesis models, the spectra were generated. For the SFH model, we decided to use the one already available in the CIGALE code with delayed exponential burst, which can model both starburst and older stellar populations. For the population synthesis module, we used the Bruzual and Charlot stellar library (Bruzual & Charlot 2003), with the Salpeter (Salpeter & Bethe 1951) initial mass function. Dust attenuation and cold dust models estimate the fraction of the cold dust in the host galaxy, its absorption of UV light, and later its thermal emission in IR. There is a discussion on which dust attenuation model is best for different sources of objects (Małek et al. 2018). We decided to use the Charlot and Fall (CF00, Charlot & Fall 2000) attenuation law, which takes into account the dust clouds and HII regions around young stars, as it allowed us to precisely model each component.

To model cold dust emission, we chose the model proposed in Dale et al. (2014), which includes the polycyclic aromatic hydrocarbon and the radiation field. Thanks to the inclusion of these elements, the model is very flexible for the modeling of photometric data. To model the AGN components, we used the SKIRTOR model (Stalevski et al. 2016). It parametrizes AGN UV-OPT emission and physical parameters of the dusty torus responsible for IR emission from the nucleus. The half-opening angle of the torus (Δ) is measured between the equatorial plane and the edge of the dusty surface of the torus. The type I quasars in the SKIRTOR model usually have Δ values of about 30°, and for type II quasars, the value of Δ is about 70° (Yang et al. 2020). CIGALE uses a slightly modified version of SKIRTOR to account for the dust in the polar regions. The additional PD component was added in the work Yang et al. (2020), calculated as a homogeneous screen, with possible different extinction curves (Buat et al. 2018). The exact parameter values used for the SED fitting procedure can be found in Table 1. The observational upper limits were taken into account to reject nonphysical models and to better constrain the main physical properties, as was shown, for example, in Junais et al. (2023). An example SED plot for one of the quasars in our sample can be seen in Fig. E.1.

3.2. Estimates of luminosities

In Rałowski et al. (2024), we introduced a simple method for calculating observed SED luminosities based on multi-photometric datasets for each galaxy (L_Y^photo, where Y is the either IR or OPT wavelength range). This quantity is computed as the integral of the SED constructed from observed data points over a given wavelength range. The purpose of this method was to provide a practical estimate of luminosity using commonly available filter sets covering the near UV to the MIR. This wavelength range is typically covered by SDSS or other optical telescopes, UKIDSS or 2MASS, and WISE. However, this set may be insufficient for full SED modeling in the absence of UV coverage and reliable MIR and FIR data points. In this work, we perform SED modeling using an extended photometric set that includes GALEX and Herschel data. Our goal is to further calibrate the L_photo estimates by comparing them with luminosities derived using the CIGALE code to obtain more reliable results based on the simpler photometric method. The L_Y^photo values were computed via integration in two wavelength ranges: from 0.11 μm to 1 μm the IR luminosity is defined (hereafter $L_{\mathrm{IR}}^{\mathrm{photo}}$), and between 1 μm and 7 μm is the optical luminosity (hereafter $L_{\mathrm{OPT}}^{\mathrm{photo}}$). We assumed that for AGNs, the dominant contribution to L_IR originates from emission by the dusty torus, while the main component of L_OPT is produced by the accretion disk.

For the second method, we used best-fit results obtained from SED fitting with CIGALE. The AGN luminosities were estimated from the SKIRTOR model. We extracted the disk luminosity ($L_{\mathrm{IR}}^{\mathrm{disk}}$) and the torus luminosity estimate ($L_{\mathrm{IR}}^{\mathrm{torus}}$). All CIGALE luminosities were derived from the best-fit synthetic spectrum and integrated over the same wavelength ranges as in the photometric method. To evaluate contaminating contributions, we included additional components in the $L_{\mathrm{OPT}}^{\mathrm{disk}}$ and $L_{\mathrm{IR}}^{\mathrm{torus}}$ luminosities, specifically stellar emission from the host galaxy, cold interstellar dust, and dust located in the polar regions of the AGN (PD). This resulted in the following luminosity definitions:

• $L_{\mathrm{OPT}}^{\mathrm{disk}}$ is based on the AGN disk component, in the wavelength range of 1 μm to 7 μm.

• $L_{\mathrm{IR}}^{\mathrm{torus}}$ is calculated from the estimated luminosity of the dusty torus, in the wavelength range of 0.11 μm to 1 μm.

• $L_{\mathrm{OPT}}^{\mathrm{disk}+\mathrm{stellar}}$ is the sum of the disk and stellar components, in the wavelength range of 1 μm to 7 μm.

• $L_{\mathrm{IR}}^{\mathrm{torus}+\mathrm{PD}}$ is the sum of the dusty torus and polar dust components, in the wavelength range of 0.11 μm to 1 μm.

• $L_{\mathrm{IR}}^{\mathrm{torus}+\mathrm{PD}+\mathrm{disk}}$ is the sum of the dusty torus, polar dust, and accretion disk, in the wavelength range of 0.11 μm to 1 μm.

The two independent luminosity estimates enable comparison and calibration of our photometric method. This comparison is shown in Fig. A.1 and is discussed in detail in Sect. 4.1. It should be noted that the uncertainties provided by CIGALE for the luminosity estimates are likely overestimated. The CIGALE code computes uncertainties only for bolometric luminosities of each model component. As a result, the reported errors correspond to the full wavelength range of the model, whereas our integration of L_IR and L_OPT in IR and OPT luminosity ranges are consecutively 0.11 μm–1 μm for IR and 1 μm–7 μm. This was done for L_disk and L_torus but also for other components (stellar, cold dust, PD).

4. Results

The primary objective of this analysis is to identify the main sources of contamination that affect the luminosities L_IR and L_OPT. We then aim to calibrate the photometric luminosity estimation method so that it more accurately reflects the dusty torus and accretion disk luminosities of AGNs as determined through CIGALE fitting.

4.1. Comparison of luminosity estimations

The integrated luminosities were computed using two distinct approaches: (1) from the photometric SED, yielding $L_{\mathrm{IR}}^{\mathrm{photo}}$ and $L_{\mathrm{OPT}}^{\mathrm{photo}}$, and (2) based on the CIGALE estimations, $L_{\mathrm{IR}}^{\mathrm{torus}}$ and $L_{\mathrm{OPT}}^{\mathrm{disk}}$. The former estimations were subsequently combined with the dominant contaminating components. In Fig. A.1 the comparison between both estimates is shown. The fitter regression is the ordinary least squares regression. The spread (σ) was calculated as the standard deviation of residuals from the regression fit. Further, the mean uncertainty of observations (σ_err) was calculated as ${\sigma }_{\mathrm{err}}=\mathrm{mean}(\sqrt{x_{\mathrm{err}}^2+y_{\mathrm{err}}^2})$, where x_err, y_err are the observation error of X and Y variables, respectively. To test what the contribution of the spread from the observational uncertainties is within the total spread of the relation, the intrinsic scatter (σ_int) was calculated as ratio of σ/σ_err. The subsequent rows in Fig. A.1 show different sources of contamination that affect the CIGALE-based estimations. A tight correlation between the two luminosity estimates, in the IR and OPT ranges, can clearly be observed. As expected, the photometric estimates exhibit systematically higher values and effectively serve as upper limits to the CIGALE-based estimates. The observed offset from the 1:1 relation, along with the scatter in the $L_{\mathrm{IR}}^{\mathrm{torus}}$ versus $L_{\mathrm{IR}}^{\mathrm{photo}}$ relation (Fig. A.1, first row) indicates the presence of additional IR contributions that cannot be attributed solely to polar dust (Fig. A.1, second row). Our analysis shows that the accretion disk emission extends beyond the OPT range and contributes significantly to the IR regime as well. As shown in Table 2, the disk contribution to $L_{\mathrm{IR}}^{\mathrm{photo}}$ is significant. After inclusion of $L_{\mathrm{IR}}^{\mathrm{torus}+\mathrm{PD}+\mathrm{disk}}$, the relation is not only closer to 1:1, but it also shows a narrower spread (low-z 0.16−0.08, inter1-z 0.10−0.06, inter2-z 0.23−0.19, high-z 0.29−0.26; see the third row in Fig. A.1). The remaining discrepancy between the data points and the 1:1 relation is attributed primarily to stellar emission and cold dust. In general, for redshifts below 1.5, the intrinsic scatter drops (from ∼1.6 to 0.7 for low-z, ∼1.2 to 0.4 for inter1-z), while for higher redshifts, the σ_int slightly increases.

In the OPT range, the correlation between the modeled disk emission and the photometric estimates follows a well-defined main trend, as illustrated in Fig. A.1, fourth row. Additionally, the $L_{\mathrm{OPT}}^{\mathrm{disk}}$ versus $L_{\mathrm{OPT}}^{\mathrm{photo}}$ relation also exhibits a substantial number of outliers across all redshift bins. To investigate how these outliers differ from the main quasar population, we conducted two dedicated tests. The first test involved fitting two least-squares regression lines: one to the full sample of quasars (solid black line, shown in Fig. A.1, fourth row) and a second excluding outliers. Outliers were defined as objects lying beyond a 1σ threshold from the first regression line. The second regression, fit to the remaining sources, is represented by the dash-dot gray line and follows the main trend more closely, with a slope approaching unity (low-z from 1.53 ± 0.13 to 1.22 ± 0.04; inter1-z from 1.39 ± 0.09 to 1.13 ± 0.03; inter2-z from 1.35 ± 0.09 to 1.13 ± 0.04; high-z from 1.61 ± 0.28 to 1.24 ± 0.10). When the stellar component is included (Fig. A.1, fifth row), the relation between the photometric and CIGALE luminosities becomes consistent with a 1:1 relation, within the uncertainties, particularly for the previously identified outliers. A subset of these outliers can thus be attributed to a significant stellar contribution, which is absent in $L_{\mathrm{OPT}}^{\mathrm{disk}}$ but is present in the next row of panels with $L_{\mathrm{OPT}}^{\mathrm{disk}+\mathrm{stellar}}$. Most of the remaining sources exhibit high f_AGN values, indicating that stellar emission plays a less significant role in their total luminosity. The ratio of each contamination can be seen in Tables 2 and 3. The outliers were responsible for the majority of the σ_int, as it drops massively for each redshift range.

It should noted that the low-luminosity outliers with a higher AGN fraction are influenced by the addition of polar dust in L_IR (see inter-z1), visible in Fig. C.2. For outliers with a lower AGN fraction, the stellar component plays a significant role in estimating the L_OPT.

4.2. Luminosity cut

The quasars shown in the fourth row of Fig. A.1 exhibit increased scatter at the low-luminosity end of $L_{\mathrm{OPT}}^{\mathrm{disk}}$. To have a universal selection of objects, we applied a luminosity cut at $L_{\mathrm{OPT}}^{\mathrm{disk}}\sim {10}^{45}$ across all redshifts, similar to the threshold adopted by Trefoloni et al. (2024). The spread σ for the objects below $L_{\mathrm{OPT}}^{\mathrm{disk}}={10}^{45}$ is almost ten times wider compared to the quasars above the threshold (low-z: 0.63−0.07, inter1-z 0.36−0.08, inter2-z 0.51−0.06, high-z 1.06−0.11 for objects below and above the $L_{\mathrm{OPT}}^{\mathrm{disk}}={10}^{45}$ threshold). This separation divides the sample into two groups of objects: (1) lower-luminosity sources with greater scatter and stronger contamination and (2) brighter objects with a more robust trend. Our goal was to examine whether different types of contamination affect these two regimes in distinct ways As shown in Fig. A.1, a fixed threshold of $L_{\mathrm{OPT}}^{\mathrm{disk}}={10}^{45}$ [erg/s] does not perfectly follow the “break point” of the OPT-IR relation for all redshift due to luminosity evolution. Therefore, we propose a redshift-dependent cut, defined as $L_{\mathrm{OPT}}^{\mathrm{disk}}={10}^{45}+0.25\times |z|$, where |z| is the mean redshift in each bin. The selection of this particular threshold was mainly empirical, as proper selection should be based on fluxes. This issue is discussed in Sect. 6 and in Fig. 5. This adjusted threshold serves two main purposes: (1) to identify the luminosity regime in which both methods, photometric and CIGALE estimated, yield consistent OPT–IR relations and (2) to calibrate the photometric method in a region of reduced scatter. Although a similar effect can be observed in the IR part, we did not impose an additional cut for it.

4.3. Comparison of OPT-IR relation through different luminosity estimations

The OPT–IR relation between $L_{\mathrm{photo}}^{\mathrm{IR}}$ and $L_{\mathrm{photo}}^{\mathrm{OPT}}$ exhibits the strongest correlation (R² ≃ 0.8) and the narrowest spread (σ ≃ 0.11). Uncertainties for the photometric method were estimated using Monte Carlo Markov chain simulations. These uncertainties increase with redshift, reflecting the larger observational errors associated with higher redshift bins. The errors for the SED fitting estimations are overestimated, as they are computed over the entire model wavelength range, rather than being limited to the relevant IR or OPT intervals (as noted in Sect. 3.1). The second row of Figure A.2 shows the relation between $L_{\mathrm{OPT}}^{\mathrm{disk}}$ and $L_{\mathrm{IR}}^{\mathrm{torus}}$. A low-luminosity tail is clearly visible, characterized by increased scatter in $L_{\mathrm{OPT}}^{\mathrm{disk}}$, particularly among low-z, low-luminosity objects. The third row shows that the addition of polar dust luminosity in $L_{\mathrm{IR}}^{\mathrm{torus}+\mathrm{polar}\mathrm{dust}}$ does not significantly alter the distribution, but it does introduce large uncertainties for most objects. These uncertainties are again calculated over the full wavelength range of each model since CIGALE does not provide error estimates for individual wavelength intervals (e.g., IR or OPT) separately. Consequently, the resulting uncertainties are likely overestimated. The uncertainties for the sum of luminosities were computed via standard error propagation methods.

In contrast, the inclusion of the stellar luminosity component in $L_{\mathrm{OPT}}^{\mathrm{disk}+\mathrm{stellar}}$ substantially increases the total L_OPT, exerting a significant influence on the group of outliers previously observed in the panel showing only $L_{\mathrm{OPT}}^{\mathrm{disk}}$. Uncertainties for high-redshift objects remain the largest among the sample. The luminosity distributions that account for the contamination sources $L_{\mathrm{OPT}}^{\mathrm{disk}+\mathrm{stellar}}$, $L_{\mathrm{IR}}^{\mathrm{torus}+\mathrm{PD}+\mathrm{disk}}$ seem to be the most reminiscent of the photometric method, which is confirmed by the slopes of the fit regressions.

Surprisingly, the distribution of the AGN fraction is mixed even among the most luminous objects. No clear correlation was observed between the AGN fraction and the total luminosity. Notably, some of the highest-luminosity sources exhibit low AGN fractions, as illustrated in Fig. C.2. This can be explained by the imperfections in the SED fitting, as after checking the PDF and fit half-opening angles, we found that most of these objects were classified as type II quasars. However, for the majority of these cases, the type I quasar solutions yielded comparable χ² values and were accompanied by lower stellar luminosity components, higher AGN fractions, and smaller opening angles. We tested whether the low AGN fraction might be due to a higher starburst fraction (f_*). For the majority of objects, the average value of f_* is approximately f_* ∼ 0.15, which is similar to the main group of objects.

The OPT-IR relation for photometric luminosities in the first row of Fig. A.2 shows significant similarities with those in the fifth row regarding all major contaminations in terms of both regression equation and the spread of the relation. The only exception is the high-z sample. This signifies that within the OPT-IR relation, there are indeed contaminations influencing the distribution. We address this issue in the next section.

5. Contamination analysis

In Sect. 4.2 we demonstrated that the photometric method is prone to contamination, especially for the objects with a low L_OPT. Therefore, to reliably investigate the relationship between the accretion disk and the torus in AGNs, it is essential to account for multiple sources of contamination, including polar dust within the AGN, stellar emission, cold dust from the host galaxy, and disk emission extending into the IR. In this section we analyze each of these major contamination sources in detail, taking into consideration their potential redshift-dependent effects. Tables 2 and 3 present a comparison between CIGALE-derived luminosity components (with additional contaminations) and the corresponding $L_{\mathrm{IR}}^{\mathrm{photo}}$ and $L_{\mathrm{OPT}}^{\mathrm{photo}}$ values. It is worth noting that $L_{\mathrm{IR}}^{\mathrm{torus}}$ and $L_{\mathrm{OPT}}^{\mathrm{disk}}$ respectively constitute the dominant contributions to $L_{\mathrm{IR}}^{\mathrm{photo}}$ and $L_{\mathrm{OPT}}^{\mathrm{photo}}$, with their fractional contributions increasing at higher redshifts (Table 2 rows 1 and 5).

5.1. Polar dust

One of the key sources of contamination analyzed in this work is dust located in the polar regions of the AGN, commonly referred to as PD. To investigate its impact, we performed SED fitting under two configurations: (1) with polar dust with parameters described in Table 1 and (2) without polar dust, i.e., parameter E(B − V)_PD was set to 0 (see also Table 1). The main results of this comparison are discussed in Sect. 5.1.2.

5.1.1. Fit with PD

Polar dust was modeled within the SED fitting framework using CIGALE, specifically through the SKIRTOR module. In this implementation, the SKIRTOR module was slightly modified by incorporating a polar dust component, represented as a blackbody emitter (for further details, see Yang et al. 2020). Polar dust has a relatively minor influence on L_IR. For the brighter objects, the polar dust luminosity is up to a few percent for all redshift ranges. For objects below the luminosity threshold, this contamination is higher (8% for low-z, 3% for inter1-z, 4% for inter2-z, and 9% for high-z).

Additionally, we examined the influence of polar dust by analyzing the ratio $L_{\mathrm{IR}}^{\mathrm{torus}}/L_{\mathrm{IR}}^{\mathrm{PD}}$. This metric quantifies the extent to which the torus luminosity dominates over the polar dust contamination. As shown in Fig. F.1, the torus is typically found to be four to six times more luminous (on average 79% across all redshift samples with $L_{\mathrm{IR}}^{\mathrm{torus}}/L_{\mathrm{IR}}^{\mathrm{PD}}$ greater than four and 40% when above five) than the polar dust component (indicated in brown and pink colors). Among the objects below the luminosity threshold, there are more objects with a lower L_torus/L_polardust ratio (only 7% objects above five). The only sample different from this trend is the inter2-z at Fig. F.1, where the majority of objects have a very high ratio. Objects above and below the luminosity threshold appear to exhibit similar values of the $L_{\mathrm{IR}}^{\mathrm{torus}}/L_{\mathrm{IR}}^{\mathrm{PD}}$ ratio. Although the most luminous sources with $L_{\mathrm{OPT}}^{\mathrm{photo}}>{10}^{45}$ [erg/s] tend to have the highest ratio, the torus remains the dominant IR component by a substantial margin. This trend seems to be consistent across redshift bins.

When comparing the reddening E(B − V)_PD caused by the presence of polar dust on the line of sight, for the objects below $L_{\mathrm{OPT}}^{\mathrm{disk}}\sim {10}^{45}$ across the full redshift range, the polar dust has a higher E(B − V)_polar dust value of 0.18, which is over two times greater than objects above the luminosity threshold. The high-z sample below the threshold has the lowest E(B − V)_polar dust of around 0.14. This is particularly evident in Fig. F.2 (lower panel), where objects below the threshold show increased dust attenuation. The estimated temperature of the polar dust in the low-z bin is, on average, the lowest among all redshift groups, with a median value of 170 K ± 70 K, and it strongly affects the FIR portion of the spectrum. inter1-z objects have a high E(B − V)_PD (Fig. F.2). These results confirm that the luminosity cut effectively excludes sources with significant polar dust attenuation, not only at low redshift but also in higher redshift bins. A similar trend is observed for inter2-z objects, where those with the highest E(B − V)_polar dust have $L_{\mathrm{OPT}}^{\mathrm{disk}}<{10}^{45}$ [erg/s] (see Fig. F.2). It should be noted that the E(B − V)_polar dust can be influenced not only by the amount of the polar dust but also by the inclination, which can significantly change its values.

5.1.2. Fit without polar dust

In the previous section, we described the characteristics of the SED fitting we performed with the inclusion of polar dust. As shown in Toba et al. (2021), the inclusion or exclusion of polar dust can significantly affect the derived model parameters. Motivated by this analysis, we conducted a second SED fitting run with polar dust excluded by setting E(B − V)_PD = 0 (see Table 1). The sample size after applying a reduced χ² cut was smaller: 50, 40, 50, and 38 sources respectively for the low-z, inter1-z, inter2-z, and high-z samples. In this configuration, the torus luminosity increased by up to 8% of the IR emission, while the disk contribution to the IR remained similar to the previous fit, at approximately 16%. (For the results of the fit without polar dust, see Tables 4 and 5). A visual inspection of the SEDs indicated that for objects previously fit with hot polar dust, the absence of this component leads to torus models with systematically higher temperatures. From the comparison between Figs. A.1 and A.4 in the $L_{\mathrm{IR}}^{\mathrm{torus}}$ versus $L_{\mathrm{IR}}^{\mathrm{photo}}$ (first row), one can see that (1) the torus luminosity is slightly higher in the model without polar dust, and (2) the overall scatter remains comparable across all redshift bins (on average σ = 0.03 compared to σ = 0.04 for the fit with polar dust). Additionally, the number of outliers based on the regression fit changed (13-4-11-4, for low-inter1-inter2-high datasets, compared to 8-8-9-10).

A detailed comparison between the two SED fitting approaches is presented in Fig. A.3. The best-fit SKIRTOR parameters differ significantly between the fits with and without polar dust. The fit without PD is characterized by the following parameter differences: (1) a steeper slope of the accretion disk emission, δ; (2) a significantly lower AGN fraction, f_AGN, with an average value of ∼70% (see Fig. A.3, second row); (3) a lower average inclination angle (i), approximately ∼30 [deg] compared to ∼40 [deg] in the fit with PD; (4) a higher opening angle, with a median of 44 [deg] compared to 30 [deg] when PD is included; and (5) the total accretion disk luminosity (fifth row in Fig. A.3) has higher errors for low luminous objects. The total AGN luminosity is slightly higher in the fit that includes polar dust, whereas the torus luminosity increases in the fit without polar dust.

Overall, considering both this result and the opening angle (panel 4, Fig. A.3), the majority of objects fit without polar dust can be interpreted as type II objects (low AGN fraction, high opening angle). In general, the inclusion of polar dust affects both the opening angle and inclination angle (i). For most objects, the fit opening angle values decrease when polar dust is included, as illustrated in Fig. 2. For redshifts in the range 0.7 < z < 2.0, a clear degeneracy is observed in the fit – manifested as a triangular distribution at an opening angle of around 40 [deg] and a diagonal trend. For high-z objects, the fit opening angle values tend to increase (typically reaching ∼50−60 [deg]), while the inclination angles decrease.

Fig. 2. Relation between opening angle on the X axis and inclination on the Y axis for the fit with polar dust (diamonds) and without (circles). The shift for each object for both fits is marked with gray arrows.

5.2. Cold dust

Cold dust from the ISM has a minor influence on $L_{\mathrm{IR}}^{\mathrm{photo}}$, as the majority of its emission lies beyond the upper wavelength limit of 7 μm used in our integration. This FIR component contributes approximately 3% to $L_{\mathrm{IR}}^{\mathrm{photo}}$ for objects above the luminosity threshold, across all redshift bins. For sources below the threshold, the impact is more pronounced: ∼8% for low-z, 12% for inter1-z, 7% for inter2-z, and 2% for high-z. In the fit without polar dust, the contribution of cold dust is reduced to around 2% across all redshifts. This is due to the dusty torus effectively absorbing the role of both the polar dust and cold dust components in the IR range.

5.3. Stellar emission

Stellar emission constitutes the primary source of contamination in $L_{\mathrm{OPT}}^{\mathrm{disk}}$ across all redshift bins. This effect is clearly visible in Fig. A.1, fourth row. All outliers (marked with stars) exhibiting a low $L_{\mathrm{OPT}}^{\mathrm{disk}}$ show a substantial stellar component. When stellar emission is added and the relation between $L_{\mathrm{OPT}}^{\mathrm{disk}+\mathrm{stellar}}$ and $L_{\mathrm{OPT}}^{\mathrm{photo}}$ is plotted, a tight correlation emerges, and the previously observed outliers disappear. The slope of the fit regression approaches unity for all redshift intervals. A lower stellar contribution is typically observed for sources with $L_{\mathrm{OPT}}^{\mathrm{disk}}>{10}^{45}$ [erg/s]. For these brighter objects, the average stellar contamination is approximately 5% (5% low-z, 4% inter1-z, 4% inter2-z, and 7% high-z). In contrast, for dimmer sources, the stellar contamination increases to around 10%, with even higher values for z > 1.5: 12% low-z, 7% inter1-z, 10% inter2-z, and 13% high-z. This trend is likely due to the reduced quality of UV photometric data at higher redshifts, which is critical for accurately separating AGN and stellar components in the UV/optical range. Overall, stellar emission remains the dominant source of contamination in L_OPT across the full sample. In the fit without polar dust, stellar contamination levels are similar for sources above the luminosity threshold: 6% for low-z, 5% for inter1-z, 6% for inter2-z, and 4% for high-z. However, for sources below the threshold, contamination is higher: 9% for low-z, 12% for inter1-z, 6% for inter2-z, and 9% for high-z.

5.4. Accretion disk emission in IR luminosity

The accretion disk emission modeled in SKIRTOR extends into the IR regime. This component of the disk continuum is visible in Fig. E.1 as an dashed orange line for the example SED. The calculated photometric luminosity is directly influenced by the definition of the wavelength ranges used for OPT and IR luminosities. Since the accretion disk emits across a broad spectral range, including the near-IR, its contribution to $L_{\mathrm{IR}}^{\mathrm{photo}}$ can be significant. Surprisingly, for sources above the luminosity threshold, the disk contributes approximately 17% to $L_{\mathrm{IR}}^{\mathrm{disk}}$. For dimmer objects, those below the threshold ($L_{\mathrm{OPT}}^{\mathrm{disk}}<{10}^{45}$ [erg/s]), the contribution from disk emission to IR decreases with redshift: 17% low-z, 21% inter1-z, 15% inter2-z, and 4% high-z. In the fit without polar dust, $L_{\mathrm{IR}}^{\mathrm{disk}}$ maintains a similar contribution of ∼17% to $L_{\mathrm{IR}}^{\mathrm{photo}}$ for brighter objects (see Table 4). For sources below the luminosity threshold, the disk contamination in IR shows a strong redshift dependence: 21% for low-z and decreasing to 1% for high-z (see Table 5).

5.5. AGN properties

To further investigate the differences between sources below and above the $L_{\mathrm{OPT}}^{\mathrm{disk}}\sim {10}^{45}$ [erg/s] threshold and SED fits performed with and without the polar dust component, we analyzed and compared the physical properties of the quasars. Our analysis focused primarily on AGN-related parameters as modeled by the SKIRTOR framework.

The vast majority of sources exhibit torus opening angle, Δ, values between 20 [deg] and 40 [deg], as can be seen in Fig. B.2, with an average of 30 [deg]. A systematic shift in the opening angle inclination, Δ − i, space can be seen in Fig. 2. Fits that include polar dust tend to yield smaller opening angles and larger inclination angles compared to those without polar dust. For fainter sources, the fit parameters are associated with significantly larger uncertainties, which can be attributed to two main factors: (1) lower flux levels and consequently higher observational errors and (2) an increased level of contamination, as discussed in the previous section.

We found that high-luminosity objects in high-z have a high opening angle (Δ ∼ 70 [deg]) and high L_IR > 10⁴⁷ [erg/s], suggesting a possible type II AGN and an ultraluminous IR galaxy. In general, high-z objects above the luminosity threshold have higher Δ values, as can be seen in Fig. F.3. Additionally, in Fig. F.4, the luminosity relation with the power-law of index δ is presented. The higher values of the δ slope imply a redder spectrum of the modeled AGN, especially in the high-luminous objects. Surprisingly, when we tried to quantify the AGN fraction (f_AGN), we found that the most luminous objects in the high-z bin have a quite low f_AGN ∼ 0.7 − 0.8. This can further suggest the hypothesis of type II sources. In Appendix B, we analyze the best-fit PDFs to check whether the second-best fit would change the type II to type I. We find that the fit for those objects is quite poor, as can be seen in Fig. B.1, where the PDF for the Δ parameter is shown. The PDF is flat, which results in a poor constraint on the Δ value.

5.6. Luminosity transformation

The main purpose of the luminosity cut at the threshold of $L_{\mathrm{OPT}}^{\mathrm{disk}}\sim {10}^{45}$ [erg/s] was to separate type I objects from outliers with significant contamination. Another goal was to check whether there is a simple transformation that changes the distribution of the photometric luminosity estimations to a distribution similar to the SKIRTOR $L_{\mathrm{disk}}^{\mathrm{OPT}}$ and $L_{\mathrm{torus}}^{\mathrm{IR}}$ estimations. This check should show whether the apparent photometry distribution can be easily calibrated to follow the AGN component’s luminosity. We used the following method:

1. Fit the ordinary least squares regression to the L_photo versus L_CIGALE for both L_IR and L_OPT in the form $L_{\mathrm{CIGALE}}^{\mathrm{X}}=a\times L_{\mathrm{photo}}^{\mathrm{X}}+b$, where X is either IR or OPT.

2. Use calculated coefficients to the L_TRANSFORM as follows: $L_{\mathrm{photo}-\mathrm{tranformed}}^{\mathrm{X}}=a\times L_{\mathrm{photo}}^{\mathrm{X}}+b$.

The results of the luminosity transformation can be seen in Figs. 3 and 4 (for objects with $L_{\mathrm{disk}}^{\mathrm{OPT}}>{10}^{45}$ [erg/s] and $L_{\mathrm{disk}}^{\mathrm{OPT}}<{10}^{45}$ [erg/s] respectively). The transformation shifts the whole distribution toward lower luminosities, which correspond to the luminosities after extraction of the major contaminations from both photometric luminosities. The change is visible in the increase of the slopes for all redshift bins. The slopes are closer to those fit to the SKIRTOR estimations. The mean luminosities for all redshift bins are also closer to those of SKIRTOR. It can be seen that the transformation based on the CIGALE-estimated luminosities does not work for objects below $L_{\mathrm{disk}}^{\mathrm{OPT}}<{10}^{45}$ [erg/s].

Fig. 3. Transformation of the luminosity relation across redshifts for objects with $L_{\mathrm{OPT}}^{\mathrm{disk}}>{10}^{45}$ [erg/s]. The first plot shows the relation between the photometric method estimated luminosities. The second plot shows the OPT-IR relation for CIGALE $L_{\mathrm{disk}}^{\mathrm{OPT}}$ and $L_{\mathrm{torus}}^{\mathrm{IR}}$. The third plot shows the $L_{\mathrm{OPT}}^{\mathrm{reduced}\mathrm{photo}}$ and $L_{\mathrm{IR}}^{\mathrm{reduced}\mathrm{photo}}$ with reduced contaminations, described in the Sect. 5.7. The fourth plot shows the transformed luminosities, described in Sect. 5.6.

Fig. 4. Transformation of luminosities relation across redshifts for objects with $L_{\mathrm{OPT}}^{\mathrm{disk}}<{10}^{45}$ [erg/s]. (For a description, see the caption of Fig. 3).

5.7. Reduction of contaminations

An alternative approach to the luminosity transformation proposed in Sect. 5.6 is to reduce the $L_{\mathrm{IR}}^{\mathrm{photo}}$ and $L_{\mathrm{OPT}}^{\mathrm{photo}}$ by a certain percentage. Based on our analysis and results shown in Table 2, we decided to reduce the luminosities for objects with $L_{\mathrm{OPT}}^{\mathrm{disk}}>{10}^{45}$ [erg/s] as follows: (1) $L_{\mathrm{IR}}^{\mathrm{Reduced}/\mathrm{photo}}=0.7\times L_{\mathrm{IR}}^{\mathrm{photo}}$ and (2) $L_{\mathrm{OPT}}^{\mathrm{Reduced}/\mathrm{photo}}=0.85\times L_{\mathrm{OPT}}^{\mathrm{photo}}$. The results of this approach can be seen in Fig. 3, third panel. As one can see, the described approach seems to yield distributions that are similar to the transformed luminosities, but compared to those estimations, it does not change the slopes of the relation. An additional advantage of this method of luminosity transformation is that it is easy to reproduce with other data. The transformations have almost no effect on the overall scatter of the relation.

6. Discussion

Here, our goal is to calibrate a simple L_OPT − L_IR relation based on more precise method of fitting the SED. By doing so, we can estimate the contributions of various components to the apparent photometric approach. As a result, it improves the accuracy of deriving AGN emission in the optical and IR bands from broadband photometry. We calculated the intrinsic luminosities of key components, including stellar emission, the accretion disk, the dusty torus, polar dust, and cold dust. Subsequently, we evaluated their respective contributions within the typical photometric coverage in the optical and near-IR bands.

Our interest was in identifying the parameter range in which a simple estimation of the disk and torus luminosities is feasible. In particular, we investigated whether integrating the SED in two fixed wavelength intervals, 0.11−1 μm for optical ($L_{\mathrm{OPT}}^{\mathrm{photo}}$) and 1−7 μm for IR bands ($L_{\mathrm{IR}}^{\mathrm{photo}}$), can be reliably used, following Rałowski et al. (2024). Figure 5 displays the threshold luminosities for $L_{\mathrm{OPT}}^{\mathrm{photo}}$ in blue and $L_{\mathrm{IR}}^{\mathrm{photo}}$ in red, below which the trend in correlation with the AGN components breaks in a given redshift bin. The approximate limits on optical log luminosity are 45, 45.2, 45.5, and 46 for the consecutive redshift bins, as shown in the first row of Fig. A.1. Above these values, $L_{\mathrm{OPT}}^{\mathrm{photo}}$ strongly correlates with $L_{\mathrm{OPT}}^{\mathrm{disk}}$, and the accretion disk dominates the SED in the 0.11−1 μm range. For the IR, the respective log luminosity limits are 45.2, 45.5, 45.8, and 46.3, as visible in the fourth row of Fig. A.1. Those limits denote the lower end of the good behaving correlation of $L_{\mathrm{IR}}^{\mathrm{photo}}$ with $L_{\mathrm{IR}}^{\mathrm{torus}}$, so they give luminosities where the dusty torus dominates in the 1−7 μm band. Within these luminosity ranges, the contribution from the host galaxy components can be considered negligible, and thus the luminosity cut allows for robust photometric estimation of AGN emission. The final size of the sample is 171 objects.

Fig. 5. Transition luminosity above which the SED fit is dominated by the AGN components in the consecutive redshift bins. The red points show $L_{\mathrm{IR}}^{\mathrm{photo}}$, which approximates the dusty torus, and the green points reveal the corresponding fluxes. The blue points visualize $L_{\mathrm{OPT}}^{\mathrm{photo}}$, which correspond to the accretion disk.

The luminosity cut, which divides the sample into sources with a well-defined correlation between observed and modeled AGN luminosity and those with a broken correlation, can be translated into fluxes, as shown in Fig. 5. Our analysis indicated that these cuts in both the optical and IR bands are primarily driven by data quality. In the 0.11−1 μm range, the most important factor is the short-wavelength point, which in most cases is the NUV from GALEX (for low-z) and u from SDSS. This point determines the CIGALE solution, as it influences the classification between UV dominated by attenuated starlight (type 2 AGN) and accretion disk dominated UV (type 1 AGN). As a result, fainter sources often exhibit prominent outliers with the fit low disk luminosity. In the IR part, long wavelength points impose their limits. The W3 and W4 WISE filters are import for dusty torus and polar dust fit accuracy. The W4 is especially known to be flaky (i.e. Rałowski et al. 2024). Additionally, Herschel data, especially for high-z sources, in the dimmer regime are challenging. For instance, their blending issue compromises the accuracy of constraining cold dust and polar dust contributions. Beyond observational issues, physical effects may also contribute to the dispersion. For example, in low-z less luminous sources are expected to be more affected by polar dust contamination (i.e. Toba et al. 2021; Rałowski et al. 2024). On the other side, z > 2 galaxies tend to exhibit higher dust content overall.

Based on Fig. 2, we find that the majority of sources in the inclination versus opening angle plane are clustered within a triangular region. This distribution suggests that many objects fall into the intermediate category, where the edge of the dusty torus is expected to be close to the line of sight. Removing the polar dust component increases the opening angle by at least 20°, while the inclination typically decreases.

The primary motivation for incorporating a polar dust component into the SED fitting procedure was to enhance the precision of addressing photometric points between 20 and 300 μm. In a subset of sources, the MIR and FIR points were insufficient. We find that polar dust can strongly influence the estimated physical properties of quasars. Although its luminosity is typically four orders of magnitude lower than that of the torus, the polar dust constitutes an important contaminant when estimating the parameters of the torus. In the presence of polar dust, the opening angle is on average smaller compared to the fit without it. This leads to a model-dependent degeneracy: either a solution with polar dust and a low opening angle or one without polar dust and a higher opening angle. This finding is similar to Toba et al. (2021), where the PD was also shown to influence the estimated covering factor, and the authors stated that the use of PD improves fit quality, based on the Bayesian information criterion. However, the authors of that study did not account for contamination from the ISM, specifically cold dust. When fitting the IR photometry up to W4, the cold dust influence in the longer wavelength starts to be significant (see Fig. E.1). It is therefore expected that including a model component dedicated to FIR emission (e.g., polar dust) leads to improved fits compared to models that exclude it. Nevertheless, even in our models without polar dust, the cold dust component alone does not fully reproduce the observed emission. Ciesla et al. (2015) reported that the AGN contribution estimates with CIGALE are often overestimated for type 1 AGN and underestimated for type 2 AGN, with the threshold being an AGN fraction of around 50%.

It has been reported that the SKIRTOR model does not produce deep enough silicate absorption features to fit the SEDs of ultraluminous IR galaxies (Efstathiou et al. 2022). In CIGALE, slight modifications to the SKIRTOR model were implemented, in particular, the inclusion of the polar dust component. However, this addition appears insufficient to resolve the issue, as the polar dust component does not contribute significantly in the wavelength range where silicate features are prominent. This limitation was evident when comparing fits with polar dust and without polar dust.

However, the photometric method of luminosity estimation is not perfect, as it contains 15−30% contaminations ($L_{\mathrm{IR}}^{\mathrm{photo}}$ and $L_{\mathrm{OPT}}^{\mathrm{photo}}$ respectively). Despite these uncertainties, the transformed luminosities agree well with SKIRTOR-based estimates of $L_{\mathrm{IR}}^{\mathrm{torus}}$ and $L_{\mathrm{OPT}}^{\mathrm{disk}}$. It should be noted that although the transformed luminosity values are close to SKIRTOR, they may still inherit the same model-dependent limitations discussed earlier. However, the transformed OPT-IR relation, visible in Fig. 3, remains robust across the sample, with a similar spread before and after the transformations, around 0.15 for joined redshifts 0.7−2.4. The slope of the transformed OPT-IR relation for all redshift bins is in the range of 0.8−0.9, which suggests its potential applicability in cosmological studies (Lusso et al. 2025). Analogous to the X-UV relation described in the introduction of this paper.

7. Conclusions

In the presented work, the OPT-IR relation based on the large sample of OPT-UV quasars was analyzed. We outline our major findings as follows:

• The major contaminations for $L_{\mathrm{IR}}^{\mathrm{photo}}$ were disk emission and cold dust, while polar dust have a minor influence. For $L_{\mathrm{OPT}}^{\mathrm{photo}}$ the major contamination was stellar emission.

• The accretion disk, starlight, and cold dust respectively contribute up to 17%, 15%, and 12% to the $L_{\mathrm{IR}}^{\mathrm{photo}}$.

• The photometric luminosity estimation method works well for objects with a luminosity above $L_{\mathrm{OPT}}^{\mathrm{disk}}\sim {10}^{45}$ [erg/s] compared to the luminosity estimations from CIGALE.

• Instead of the cut on $L_{\mathrm{OPT}}^{\mathrm{disk}}$ we can minimize the influence of contaminations on optical and IR luminosity by using redshift-dependent cuts on $L_{\mathrm{OPT}}^{\mathrm{photo}}$ and $L_{\mathrm{IR}}^{\mathrm{photo}}$.

• The population of type 1 AGNs can be contaminated by misclassified sources of type 2 – L_disk outliers. This type of misclassification can result in a drop in the accuracy of the SED fitting below the luminosity threshold. This effect can be strongly suppressed by imposing luminosity cuts.

• The proposed transformation for photometric luminosities produces the OPT-IR relation with a distribution similar to $L_{\mathrm{IR}}^{\mathrm{torus}}$, $L_{\mathrm{OPT}}^{\mathrm{disk}}$ estimated with the SKIRTOR model.

• The inclusion of the polar dust model changes the estimation of the physical properties of the quasars, such as the opening angle, inclination, and AGN fraction. This can be seen in the comparison between the fit with polar dust and the fit without polar dust (see Fig. A.3). However, it should be noted that the effect on luminosities is minor.

• The polar dust, which can be thought of as the extended atmosphere of the torus, is important in the lower luminosity end of the population (especially low-z, where it contributes up to 10% L_IR).

• Accounting for the polar dust improves fit quality and affects the inclination and the torus thickness estimations.

We believe that we identified the major contaminations within the $L_{\mathrm{OPT}}^{\mathrm{photo}}$ and $L_{\mathrm{IR}}^{\mathrm{photo}}$. The calibration procedure for the OPT-IR relation is crucial to using it in the future as a cosmological probe.

Data availability

Tables with photometrical data are available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/708/A31

Acknowledgments

MR has been supported by the Polish National Agency for Academic Exchange (Bekker grant BPN/BEK/2024/1/00298). This research was supported by the Polish National Science Center grants 2018/30/M/ST9/00757 and 2023/50/A/ST9/00579 and by Polish Ministry of Science and Higher Education grant DIR/WK/2018/12. This research was funded by the Priority Research Area Digiworld under the program Excellence Initiative – Research University at the Jagiellonian University in Kraków. We gratefully acknowledge Polish high-performance computing infrastructure PLGrid (HPC Center: ACK Cyfronet AGH) for providing computer facilities and support within computational grant no. PLG/2024/017491. This article is based upon work from COST Action CA21136 – “Addressing observational tensions in cosmology with systematics and fundamental physics (CosmoVerse)”, supported by COST (European Cooperation in Science and Technology). This research has made use of the NASA/IPAC Infrared Science Archive, which is funded by the National Aeronautics and Space Administration and operated by the California Institute of Technology. Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA.

References

Appendix A: Plots of luminosity estimations and OPT-IR relations

In this appendix we present the plots referenced in the main body of the article. Figure A.1 illustrates the relationship between photometrical luminosity estimates and their estimates from CIGALE, with and without components from the host-galaxy. Starting from the L${}_{\mathit{IR}}^{\mathrm{torus}}$ versus L${}_{\mathit{IR}}^{\mathrm{photo}}$, we add contamination components estimated through the SED fitting are added sequentially, culminating in L${}_{\mathrm{IR}}^{\mathrm{torus}+\mathrm{PD}+\mathrm{disk}}$, where PD and disk stand for polar dust and accretion disk emission in the IR range. After including these components, the relation between modeled and photometric IR luminosities becomes tighter, with a slope approaching unity (see rows 1 to 3 of Fig. A.1). A similar comparison is shown for the optical range between L${}_{\mathrm{OPT}}^{\mathrm{disk}}$ and L${}_{\mathrm{OPT}}^{\mathrm{photo}}$ is shown (row 4 in Fig. A.1). This time a significant number of outliers is visible across a redshift bins, being the reason for the luminosity cut at L${}_{\mathrm{OPT}}^{\mathrm{disk}}={10}^{45}+|z|$, where |z| is the mean redshift for each bin. Only objects above this threshold are considered in the main analysis. Stellar emission from the host galaxy is identified as the dominant source of contamination in the optical range. Further details are provided in Sects. 4.1, 4.2, and 4.3.

Fig. A.1. Comparison of luminosities estimated from 1) our photometric method (presented on X-axis) and 2) CIGALE estimates without and with contaminations for four redshift bins. The dashed line represents the 1:1 relation, while the solid black line corresponds to the best-fit ordinary least squares regression. In the fourth row, outliers were identified as objects lying below the 1σ threshold (indicated by green lines) and are marked with star symbols. These outliers were excluded, and a second ordinary least squares regression was fit to the remaining data points (dash-dot gray line). Additionally, in the fourth row, the dashed red line indicates the luminosity cut described in Sect. 4.2. The σ is the spread of the relation.

Figure A.2 shows relation between photometrical estimations of the L${}_{\mathrm{IR}}^{\mathrm{photo}}$ versus L${}_{\mathrm{OPT}}^{\mathrm{photo}}$ and the luminosity-luminosity relation for the CIGALE based estimations. Only after taking into account all contamination components the relation L${}_{\mathrm{IR}}^{\mathrm{torus}+\mathrm{PD}+\mathrm{disk}}$ versus L${}_{\mathrm{OPT}}^{\mathrm{disk}+\mathrm{stellar}}$ is similar to photometrical in terms of slope and scatter. For a detailed description, please see Sect. 4.3.

Fig. A.2. OPT-IR relation for all redshift bins through different estimates of luminosities. The solid black lines are the ordinary least squares regressions. In the first row the LIR and LOPT are estimated with the photometric method. The legend is explained under previous plot.

Although polar dust is not the dominant source of contamination it significantly affects the fit AGN model parameters, particularly inclination and accretion disk slope. The impact of polar dust was assessed through two separate CIGALE runs, performed with and without this component. A comparison of the resulting model parameters and AGN luminosities is shown in Fig. A.3. Figure A.4 presents an analog of Fig. A.2 for fits excluding polar dust. detailed discussion is provided in Sect. 5.1.2.

Fig. A.3. Comparison of SKIRTOR model parameters derived from two independent SED fitting runs: 1) With polar dust, and 2) without polar dust (E(B − V)PD = 0). The X-axis shows values from the fit including polar dust, while the Y-axis corresponds to the fit without polar dust. The plotted parameter is indicated in the panel title. The dash-dotted black line represents the 1:1 relation.

Fig. A.4. No PD: Comparison of luminosities estimated from 1) our photometric method (presented on X-axis) and 2) CIGALE estimates without and with contaminations for 4 redshift bins. The dashed lines is the 1:1 relation, while the solid black line is the best-fit regression. In the fourth row the outliers were also accounted for, defined as objects below 1σ error (green lines), marked with stars. Then the outliers were extracted and another regression was fit (dashed gray line).

Appendix B: Testing Cigale PDFs

After running an additional test for outliers after the luminosity cut described in Sect. 4.2, we considered whether this group of objects is wrongly fit by the Cigale. Our hypothesis was, that some of those objects were wrongly classified as type II source, and the solution for type I had only slightly smaller probability distribution function (PDF) for values of parameters described as type I. We performed additional SED fitting to calculate PDFs for SKIRTOR and stellar models. The example PDFs are shown in Fig. B.1. As on can see the χ² values are extremely similar for this object, and it is unclear whether the chosen value 60 [deg] is truly the best one.

Fig. B.1. Example PDF plot for Δ parameter for one of the problematic galaxies at high-z.

Fig. B.2. Δ error analysis. On Y axis the Δ is presented. On X axis the difference between Δerr and median of Δerr is shown. It is visible that 1) for high-z the general error is higher, 2) the errors increase, when CIGALE estimates are further from given a priory values.

Appendix C: OPT-IR relation versus AGN fraction

Figures C.1 and C.2 are modifications of Fig. A.1 and A.2, the color representing the f_AGN. Figure C.1 shows that the relation with f_AGN is not obvious. On one side the low values of f_AGN are present with low luminosity outliers (see the third row, which would mean that those objects could be recognized as type II AGN). On the other hand, few of high luminous objects also have a low f_AGN of < 0.8.

Fig. C.1. Comparison of luminosities estimated from 1) our photometric method (presented on X-axis) and 2) CIGALE estimates without and with contaminations for 4 redshift bins.

Fig. C.2. OPT-IR relation for all redshift bins through different estimates of luminosities. In the first row the LIR and LOPT are estimated with photometric method. For the relation to be robust the cut in the luminosities may be necessary.

Appendix D: Star formation rate versus stellar mass

We checked the estimated stellar masses (M_*) and star formation rate (SFR) on Fig. D.1. In recent work by Buat et al. (2021) the authors also analyzed SED for quasars with CIGALE. We wanted to compare the received distributions of M_* and SFR to the results therein. Buat et al. (2021) calculated the SED with, and without polar dust component, and then compare SFR and M_* histograms (see Fig. 9 therein). Distributions for our samples in Fig. D.1 resemble values calculated with polar dust.

Fig. D.1. SFR vs. M* for our sample of quasars.

Appendix E: SED example

Below we present the SED plot generated with CIGALE code in Fig. E.1. The plot was slightly modified to calculate the reduced χ² in IR and OPT domains. The solid black line represents the model spectrum which is the sum of stellar, AGN, cold dust models.

Fig. E.1. Example SED for one of quasars in observer frame. The upper panel shows the best-fit SED model, which is showed with thick black line. The components of model spectrum are stellar emission (both attenuated - solid yellow line, and unattenuated - dashed blue line; Bruzual & Charlot 2003), AGN emission with orange color (Stalevski et al. 2016), cold dust emission with solid red line (Dale et al. 2014). The red dots mark the modeled flux observations, based on the broadband filters transmissions. The pink circles mark the photometrical observations in the observed-frame of the objects. The border of integration between OPT and IR luminosities is at 1μm, in rest-frame and 2.04μm in observer frame.

Appendix F: Diagnostic plots

The additional diagnostic plots for the L${}_{\mathit{OPT}}^{\mathit{photo}}$ versus L${}_{\mathit{IR}}^{\mathit{photo}}$ relation, showing with the color the E(B-V)_PD in Fig. F.2 and the ratio of L_torus/L_PD visible with color in Fig. F.1. Figure F.2 show that a number of objects below the luminosity threshold have the highest values in E(B-V)_PD. There are still some pbjects with the lowest luminosities above the threshold with similarly high E(B-V)_PD, which suggests that not all of the contaminated objects were excluded with the luminosity cut, but for the majority of objects the luminosity cut work well, as on average the objects with lower values of E(B-V)_PD are left.

Fig. F.1. The LIR vs. LOPT relation for all data samples. Upper panel shows the objects with L${}_{\mathit{OPT}}^{\mathit{disk}}>{10}^{45}$ [erg/s], while lower shows dimmed objects with L${}_{\mathit{OPT}}^{\mathit{disk}}<{10}^{45}$ [erg/s]. With the color the ratio of L${}_{\mathit{IR}}^{\mathit{polardust}}$/L${}_{\mathit{IR}}^{\mathit{torus}}$ is shown.

Figure F.1 shows that the objects above the luminosity threshold on average have torus luminosity 5 or more times higher than the PD luminosity. The objects below the threshold have lower ratio, though still few times higher than PD luminosity.

Fig. F.2. LIR vs. LOPT relation for all data samples. Upper panel shows the objects with L${}_{\mathit{OPT}}^{\mathit{disk}}>{10}^{45}$ [erg/s], while lower shows dimmed objects with L${}_{\mathit{OPT}}^{\mathit{disk}}<{10}^{45}$ [erg/s]. With the color the E(B-V) for polar dust is shown.

Fig. F.3. LIR vs. LOPT relation for high-z sample. Upper panel shows the objects with L${}_{\mathit{OPT}}^{\mathit{disk}}>{10}^{45}$ [erg/s], while lower shows dimmed objects with L${}_{\mathit{OPT}}^{\mathit{disk}}<{10}^{45}$ [erg/s]. With the color the half opening angle of the dusty torus is shown.

Fig. F.4. LIR vs. LOPT relation for high-z sample. Upper panel shows the objects with L${}_{\mathit{OPT}}^{\mathit{disk}}>{10}^{45}$ [erg/s], while lower shows dimmed objects with L${}_{\mathit{OPT}}^{\mathit{disk}}<{10}^{45}$ [erg/s]. With the color the ’delta’ slope of SKIRTOR is shown.

All Tables

All Figures

	Fig. 1. Flowchart showing the data selection process as described in Sect. 2. The χ2 values were calculated for the best-fit SED in CIGALE. The Ldisk values are from the best-fit estimation of the SKIRTOR model, as described in Sect. 3.2.
In the text

	Fig. 2. Relation between opening angle on the X axis and inclination on the Y axis for the fit with polar dust (diamonds) and without (circles). The shift for each object for both fits is marked with gray arrows.
In the text

Fig. 3. Transformation of the luminosity relation across redshifts for objects with $L_{\mathrm{OPT}}^{\mathrm{disk}}>{10}^{45}$ [erg/s]. The first plot shows the relation between the photometric method estimated luminosities. The second plot shows the OPT-IR relation for CIGALE $L_{\mathrm{disk}}^{\mathrm{OPT}}$ and $L_{\mathrm{torus}}^{\mathrm{IR}}$. The third plot shows the $L_{\mathrm{OPT}}^{\mathrm{reduced}\mathrm{photo}}$ and $L_{\mathrm{IR}}^{\mathrm{reduced}\mathrm{photo}}$ with reduced contaminations, described in the Sect. 5.7. The fourth plot shows the transformed luminosities, described in Sect. 5.6.

In the text

	Fig. 4. Transformation of luminosities relation across redshifts for objects with $L_{\mathrm{OPT}}^{\mathrm{disk}}<{10}^{45}$ [erg/s]. (For a description, see the caption of Fig. 3).
In the text

	Fig. 5. Transition luminosity above which the SED fit is dominated by the AGN components in the consecutive redshift bins. The red points show $L_{\mathrm{IR}}^{\mathrm{photo}}$, which approximates the dusty torus, and the green points reveal the corresponding fluxes. The blue points visualize $L_{\mathrm{OPT}}^{\mathrm{photo}}$, which correspond to the accretion disk.
In the text

Fig. A.1. Comparison of luminosities estimated from 1) our photometric method (presented on X-axis) and 2) CIGALE estimates without and with contaminations for four redshift bins. The dashed line represents the 1:1 relation, while the solid black line corresponds to the best-fit ordinary least squares regression. In the fourth row, outliers were identified as objects lying below the 1σ threshold (indicated by green lines) and are marked with star symbols. These outliers were excluded, and a second ordinary least squares regression was fit to the remaining data points (dash-dot gray line). Additionally, in the fourth row, the dashed red line indicates the luminosity cut described in Sect. 4.2. The σ is the spread of the relation.

In the text

	Fig. A.2. OPT-IR relation for all redshift bins through different estimates of luminosities. The solid black lines are the ordinary least squares regressions. In the first row the LIR and LOPT are estimated with the photometric method. The legend is explained under previous plot.
In the text

	Fig. A.3. Comparison of SKIRTOR model parameters derived from two independent SED fitting runs: 1) With polar dust, and 2) without polar dust (E(B − V)PD = 0). The X-axis shows values from the fit including polar dust, while the Y-axis corresponds to the fit without polar dust. The plotted parameter is indicated in the panel title. The dash-dotted black line represents the 1:1 relation.
In the text

Fig. A.4. No PD: Comparison of luminosities estimated from 1) our photometric method (presented on X-axis) and 2) CIGALE estimates without and with contaminations for 4 redshift bins. The dashed lines is the 1:1 relation, while the solid black line is the best-fit regression. In the fourth row the outliers were also accounted for, defined as objects below 1σ error (green lines), marked with stars. Then the outliers were extracted and another regression was fit (dashed gray line).

In the text

	Fig. B.1. Example PDF plot for Δ parameter for one of the problematic galaxies at high-z.
In the text

	Fig. B.2. Δ error analysis. On Y axis the Δ is presented. On X axis the difference between Δerr and median of Δerr is shown. It is visible that 1) for high-z the general error is higher, 2) the errors increase, when CIGALE estimates are further from given a priory values.
In the text

	Fig. C.1. Comparison of luminosities estimated from 1) our photometric method (presented on X-axis) and 2) CIGALE estimates without and with contaminations for 4 redshift bins.
In the text

	Fig. C.2. OPT-IR relation for all redshift bins through different estimates of luminosities. In the first row the LIR and LOPT are estimated with photometric method. For the relation to be robust the cut in the luminosities may be necessary.
In the text

	Fig. D.1. SFR vs. M* for our sample of quasars.
In the text

Fig. E.1. Example SED for one of quasars in observer frame. The upper panel shows the best-fit SED model, which is showed with thick black line. The components of model spectrum are stellar emission (both attenuated - solid yellow line, and unattenuated - dashed blue line; Bruzual & Charlot 2003), AGN emission with orange color (Stalevski et al. 2016), cold dust emission with solid red line (Dale et al. 2014). The red dots mark the modeled flux observations, based on the broadband filters transmissions. The pink circles mark the photometrical observations in the observed-frame of the objects. The border of integration between OPT and IR luminosities is at 1μm, in rest-frame and 2.04μm in observer frame.

In the text

	Fig. F.1. The LIR vs. LOPT relation for all data samples. Upper panel shows the objects with L${}_{\mathit{OPT}}^{\mathit{disk}}>{10}^{45}$ [erg/s], while lower shows dimmed objects with L${}_{\mathit{OPT}}^{\mathit{disk}}<{10}^{45}$ [erg/s]. With the color the ratio of L${}_{\mathit{IR}}^{\mathit{polardust}}$/L${}_{\mathit{IR}}^{\mathit{torus}}$ is shown.
In the text

	Fig. F.2. LIR vs. LOPT relation for all data samples. Upper panel shows the objects with L${}_{\mathit{OPT}}^{\mathit{disk}}>{10}^{45}$ [erg/s], while lower shows dimmed objects with L${}_{\mathit{OPT}}^{\mathit{disk}}<{10}^{45}$ [erg/s]. With the color the E(B-V) for polar dust is shown.
In the text

	Fig. F.3. LIR vs. LOPT relation for high-z sample. Upper panel shows the objects with L${}_{\mathit{OPT}}^{\mathit{disk}}>{10}^{45}$ [erg/s], while lower shows dimmed objects with L${}_{\mathit{OPT}}^{\mathit{disk}}<{10}^{45}$ [erg/s]. With the color the half opening angle of the dusty torus is shown.
In the text

	Fig. F.4. LIR vs. LOPT relation for high-z sample. Upper panel shows the objects with L${}_{\mathit{OPT}}^{\mathit{disk}}>{10}^{45}$ [erg/s], while lower shows dimmed objects with L${}_{\mathit{OPT}}^{\mathit{disk}}<{10}^{45}$ [erg/s]. With the color the ’delta’ slope of SKIRTOR is shown.
In the text