How Does the Polar Dust Affect the Correlation between Dust Covering Factor and Eddington Ratio in Type 1 Quasars Selected from the Sloan Digital Sky Survey Data Release 16?

A flowchart of our sample selection process is shown in Figure 1. In this work, we focus on spectroscopically confirmed type 1 quasars whose L_bol, M_BH, and λ_Edd can be securely estimated. We sampled quasars that were drawn from the DR16Q (v4 ¹¹ ), which contains 750,414 type 1 quasars up to z ∼ 6.0.

Zoom In Zoom Out Reset image size

Aside from the SDSS imaging data, 3.4 and 4.6 μm data taken from the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010) and information on source variability from the Palomar Transient Factory (PTF; Law et al. 2009; Rau et al. 2009) were used for target selection of quasars (see Sections 2.1 and 2.2 in Lyke et al. 2020, for more details). The targeted objects were automatically classified as QSO, STAR, or GALAXY based on the SDSS spectra. In addition, objects that were initially classified as QSO at z_pipe > 3.5 by the SDSS pipeline were reclassified for visual inspection (see Section 3.2 in Lyke et al. 2020).

Compared to previous SDSS quasar catalogs (e.g., DR7Q, Schneider et al. 2010; DR9Q, Pâris et al. 2012; DR14Q, Pâris et al. 2018), DR16Q extends coverage of luminosity to much less luminous quasars given a redshift (see Figures 7 in Pâris et al. 2018 and Lyke et al. 2020). Hence, our quasar sample is expected to cover a wide range of L_bol, M_BH, and λ_Edd compared with those derived from previous SDSS quasar catalogs employed in previous works on CF (e.g., Gu 2013; Ma & Wang 2013; Roseboom et al. 2013).

We first narrowed down the sample to objects with ZWARNING=0, which ensures a confident spectroscopic classification and redshift measurement for quasars (Bolton et al. 2012). We then extracted objects with 0 < z < 0.7 because we aim to estimate M_BH based on a recipe using Hβ (see Section 2.4) that is detectable for objects at z up to ∼0.7 given a spectral coverage of the SDSS spectrograph. Note that it is possible to measure M_BH for objects at z > 0.7 based on recipes using other emission lines (such as C iv and Mg ii). However, to avoid being affected by possible systematic uncertainties of M_BH because of using different recipes, we adopted the same recipe for all quasars and performed a correlation analysis (see Sections 2.4 and 4.1). Eventually, 37,181 quasars were selected as a sample. We note that the size of the quasar sample at z < 0.7 selected from the DR16Q is significantly increased by ∼15,000 from DR14Q.

For those 37,181 quasars, we compiled optical to mid-IR (MIR) photometry, most of which are already available in DR16Q. In addition to SDSS optical data (u, g, r, i, and z) corrected for Galactic extinction (Schlafly & Finkbeiner 2011), we utilized the IR data taken with the Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006), the UKIRT Infrared Deep Sky Survey (UKIDSS; Lawrence et al. 2007), and WISE (see Section 7 in Lyke et al. 2020 for a complete description of the above data).

Note that DR16Q newly employed unWISE (Lang 2014) instead of the ALLWISE catalog (Cutri et al. 2014). unWISE provides deep 3.4 and 4.6 μm data that are force-photometered at the locations of SDSS sources (Lang et al. 2016). On the other hand, WISE 12 and 22 μm data are still quite useful for this work because those wavelengths are expected to trace the emission from the dusty torus (e.g., Toba et al. 2014). Hence, we compiled those MIR data from ALLWISE with a search radius of 2'' in the same manner as Pâris et al. (2018).

Consequently, we have at maximum 13 photometric data (u, g, r, i, z, Y, J, H, K/K_s, and 3.4, 4.6, 12, and 22 μm) for the SED fitting (see Section 2.3). Among 37,181 quasars, 3345 (9.0%), 5515 (14.8%), 34,683 (93.3%), and 37,155 (99.9%) objects are detected by 2MASS, UKIDSS, ALLWISE, and unWISE, respectively. We used profile-fit photometry for 2MASS-detected sources with rd_flg = 2 ¹² and that for WISE 12 and/or 22 μm detected sources with cc_flag = 0 ¹³ (which ensures clean photometry without being affected by possible artifacts) at each band. If an object lies outside the UKIDSS footprint, its NIR flux densities are taken from 2MASS (if that object is detected). Otherwise, we always refer to the UKIDSS as NIR data.

To derive IR luminosity contributed only from AGN dust torus (${L}_{\mathrm{IR}}^{\mathrm{torus}}$), we conducted the SED fitting by considering the energy balance between the UV/optical and IR. We employed a new version of Code Investigating GALaxy Emission (CIGALE; Burgarella et al. 2005; Noll et al. 2009; Boquien et al. 2019), so-called X-CIGALE ¹⁴ (Yang et al. 2020). This code enables us to handle many parameters, such as star formation history (SFH), single stellar population (SSP), attenuation law, AGN emission, dust emission, and radio synchrotron emission (see, e.g., Boquien et al. 2014, 2016; Buat et al. 2014, 2015; Ciesla et al. 2017; Lo Faro et al. 2017; Toba et al. 2019b; Burgarella et al. 2020). Aside from an implementation of the SED fitting even for X-ray data, X-CIGALE incorporates a clumpy two-phase torus model (SKIRTOR ¹⁵ ; Stalevski et al. 2012, 2016) and a polar dust emission as AGN templates. According to Yang et al. (2020), we consider the dust responsible for type 1 AGN obscuration as polar dust (i.e., the polar dust provides obscuration for the nucleus and the additional IR emission from the absorbed energy; see Figure 4 in Yang et al. 2020 for a schematic view of polar dust). Parameter ranges used in the SED fitting are tabulated in Table 1.

Table 1. Parameter Ranges Used in the SED Fitting with X-CIGALE

Parameter	Value
Delayed SFH
τmain (Gyr)	0.5, 1.0, 2.0, 4.0, 6.0, 8.0
age (Gyr)	5.0
SSP (Bruzual & Charlot 2003)
IMF	Chabrier 2003
Metallicity	0.02
Dust Attenuation (Charlot & Fall 2000)
${A}_{{\rm{V}}}^{\mathrm{ISM}}$	0.01, 0.1, 0.2, 0.3, 0.4, 0.6, 0.8, 1.0
slope_ISM	−0.7
slope_BC	−1.3
AGN Disk +Torus Emission (Stalevski et al. 2016)
τ9.7	3, 5, 9
p	0.0, 1.0, 1.5
q	0.0, 1.0, 1.5
Δ (deg)	10, 30, 50
${R}_{\max }/{R}_{\min }$	30
θ (°)	0
fAGN	0.5, 0.6, 0.7, 0.8, 0.9, 0.99
AGN Polar Dust Emission (Yang et al. 2020)
extinction low	SMC
E(B − V)	0.0, 0.05, 0.1, 0.15, 0.2, 0.3, 0.4
${T}_{\mathrm{dust}}^{\mathrm{polar}}$ (K)	100.0, 500.0, 1000.0
Emissivity β	1.6

Download table as:  ASCIITypeset image

We adopted a delayed SFH model, assuming a single starburst with an exponential decay (e.g., Ciesla et al. 2015, 2016), where we fixed the age of the main stellar population in the galaxy to be the same as what we used for the galaxy template when fitting to the SDSS spectra (see Section 2.4). We parameterized the e-folding time of the main stellar population (τ_main). The influence of fixing the age of the main stellar population on the CF is discussed in Section 4.3.

We chose the SSP model (Bruzual & Charlot 2003), assuming the IMF of Chabrier (2003), and the standard nebular emission model included in X-CIGALE (see Inoue 2011).

For attenuation of dust associated with the host galaxy, we utilized a model provided by Charlot & Fall (2000) that has two different power-law attenuation curves: the power-law slope of the attenuation in the interstellar medium (ISM), and birth clouds (BC) that were fixed in this work. We parameterized the V-band attenuation in the ISM (${A}_{{\rm{V}}}^{\mathrm{ISM}}$).

AGN emission from the accretion disk and dust torus was modeled by using SKIRTOR, a clumpy two-phase torus model produced in the framework of the 3D radiative transfer code, SKIRT (Baes et al. 2011; Camps & Baes 2015). This torus model consists of seven parameters: torus optical depth at 9.7 μm (τ_9.7), torus density radial parameter (p), torus density angular parameter (q), angle between the equatorial plane and edge of the torus (Δ), ratio of the maximum to minimum radii of the torus (${R}_{\max }/{R}_{\min }$), the viewing angle (θ), and the AGN fraction in total IR luminosity (f_AGN). In order to avoid a degeneracy of AGN templates (see Yang et al. 2020) and to estimate a reliable CF (see Section 3.3), we fixed ${R}_{\max }/{R}_{\min }$ and θ that are optimized for type 1 quasars. We note that ${R}_{\max }/{R}_{\min }$ and θ are insensitive to the CF at least for type 1 AGNs as Stalevski et al. (2012) reported.

In X-CIGALE, the polar dust emission is implemented as a graybody emission with an empirical extinction curve (see Section 2.4 in Yang et al. 2020, for more details). The graybody is formulated as 1 − e^τ(ν) ${\nu }^{\beta }\,{B}_{\nu }({T}_{\mathrm{dust}}^{\mathrm{polar}})$, where ν is the frequency, β is the emissivity index of the dust, and ${B}_{\nu }({T}_{\mathrm{dust}}^{\mathrm{polar}})$ is the Planck function. $\tau \equiv {(\nu /{\nu }_{0})}^{\beta }$ is the optical depth. ν₀ is the frequency where optical depth equals unity (Draine 2006). ν₀ = 1.5 THz (λ₀ = 200 μm) in X-CIGALE, and β is fixed to be 1.6 (e.g., Casey 2012), and thus we parameterized ${T}_{\mathrm{dust}}^{\mathrm{polar}}$ in the same manner as Toba et al. (2020b). We adopted the SMC extinction curve (Prevot et al. 1984) that may be reasonable for AGN dust (e.g., Pitman et al. 2000; Kuhn et al. 2001; Salvato et al. 2009), in which we parameterized E(B − V). E(B − V) = 0 means that the polar dust component is not necessary. We confirmed that our conclusion is not significantly affected if we assumed another extinction curve.

In the framework of our SED modeling, IR luminosity can be described by a combination of stellar component (${L}_{\mathrm{IR}}^{\mathrm{stellar}}$), AGN accretion disk component (${L}_{\mathrm{IR}}^{\mathrm{disk}}$), AGN dust torus component (${L}_{\mathrm{IR}}^{\mathrm{torus}}$), and AGN polar dust component (${L}_{\mathrm{IR}}^{\mathrm{polar}}$), which enables us to estimate ${L}_{\mathrm{IR}}^{\mathrm{torus}}$ that is relevant to the CF through the SED decomposition.

Under the parameter setting described in Table 1, we fit the stellar and AGN (accretion disk, dust torus, and polar dust) components to at most 13 photometric data as described in Section 2.2. Although the photometry employed in DR16Q is not identical, the profile-fit photometry is used for the SDSS and 2MASS (Lyke et al. 2020), while the UKIDSS and unWISE data are forced photometry at the SDSS centroids (Aihara et al. 2011; Lang 2014). Therefore, the flux densities in the optical–MIR bands in DR16Q are expected to trace the total flux densities, suggesting that the influence of different photometry is likely to be small (see also Section 4.5). We used flux density in a band when the signal-to-noise ratio (S/N) is greater than 3 at that band. Following Toba et al. (2019a), we put 3σ upper limits for objects with ph_qual = "U" (which means upper limit on magnitude) in the MIR bands.

In order to derive M_BH, L_bol, and λ_Edd of our quasar sample, we conducted a spectral fitting to SDSS spectra by using the Quasar Spectral Fitting package (QSFit v1.3.0; Calderone et al. 2017). ¹⁶ This code fits optical spectra by taking into account (i) AGN continuum with a single power-law, (ii) Balmer continuum modeled by Grandi (1982) and Dietrich et al. (2002), (iii) the host galaxy component with an empirical SED template, (iv) iron blended emission lines with UV–optical templates (Vestergaard & Wilkes 2001; Véron-Cetty et al. 2004), and (v) emissions line with Gaussian components. QSFit fits all the components simultaneously following a Levenberg–Marquardt least-squares minimization algorithm with the MPFIT (Markwardt 2009) procedure.

The main purpose of this spectral fitting is to measure the FWHM of the Hβ broad component (FWHM_Hβ ) and continuum luminosity at 5100 Å (L₅₁₀₀), which are ingredients for M_BH estimates. For the host galaxy contribution, ¹⁷ we employed a simulated 5 Gyr old elliptical galaxy template (Silva et al. 1998; Polletta et al. 2007), allowing the normalization to vary. After correcting the Galactic extinction provided by Schlafly & Finkbeiner (2011), which is the same as what we adopted for the SDSS photometry (see Section 2.2), spectral fitting was executed. Following Calderone et al. (2017), we fit the Hβ line with narrow and broad components in which the FWHM of narrow and broad components is constrained in the range of 100–1000 km s⁻¹ and 900–15,000 km s⁻¹, respectively, to allow a good decomposition of the line profile. ¹⁸

Based on outputs from QSFit, we estimated M_BH with a single-epoch method reported by Vestergaard & Peterson (2006):

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}\left(\displaystyle \frac{{M}_{\mathrm{BH}}}{{M}_{\odot }}\right) & = & \mathrm{log}\left[{\left(\displaystyle \frac{{\mathrm{FWHM}}_{{\rm{H}}\beta }}{1000\mathrm{km}{{\rm{s}}}^{-1}}\right)}^{2}{\left(\displaystyle \frac{{L}_{5100}}{{10}^{44}\mathrm{erg}{{\rm{s}}}^{-1}}\right)}^{0.50}\right]\\ & & +(6.91\pm 0.02).\end{array}\end{eqnarray} \tag{ 1 }$$

L_bol was converted from BC₅₁₀₀ × L₅₁₀₀ where BC₅₁₀₀ = 8.1 ± 0.4 is the bolometric correction (Runnoe et al. 2012). The uncertainty in M_BH is calculated through error propagation of Equation (1), while the uncertainty in L_bol is propagated from 1σ errors of L₅₁₀₀ and BC₅₁₀₀.

Figure 2 shows examples of the SED fitting with X-CIGALE. We confirm that 34,541/37,181 (∼92.9%) objects have reduced ${\chi }_{\mathrm{CIGALE}}^{2}\lt 3.0$, meaning that the data are well fitted with the combination of the stellar, nebular, AGN accretion disk, dust torus, and polar dust components by X-CIGALE.

Zoom In Zoom Out Reset image size

On the other hand, it is worth investigating whether the polar dust component (which is the main topic in this work) is practically needed to improve the SED fitting. In order to test the requirement to add an AGN polar dust component to the SED fitting, we calculate the Bayesian information criterion (BIC; Schwarz 1978) for two fits that are derived with and without the polar dust component. The BIC is defined as BIC = χ² + k × ln(n), where χ² is the nonreduced chi-square, k is the number of degrees of freedom (dof), and n is the number of photometric data points used for the fitting. We then compare the results of two SED fittings without/with the polar dust component by using ΔBIC = BIC_wopolar—BIC_wpolar. The resultant ΔBIC tells whether or not the AGN polar dust component is required to give a better fit, taking into account the difference in dof (e.g., Ciesla et al. 2018; Buat et al. 2019; Aufort et al. 2020; Toba et al. 2020a).

Figure 3 shows the histogram of ΔBIC for our quasar sample. We adopt ΔBIC = 6 as a threshold to consider the differences in two fits in the same manner as Buat et al. (2019). If ΔBIC >6, this indicates that adding the AGN polar dust component significantly improves the fit. We find that about 30.0% of objects satisfy ΔBIC > 6. This suggests that for a faction of objects in our SDSS quasar sample at z < 0.7, an AGN polar dust component may not be necessary. We note that the threshold we adopted in this work is conservative compared to what was originally suggested, ΔBIC = 2 (Liddle 2004; Stanley et al. 2015). Nevertheless, the above result could suggest that our SED modeling may not be enough to constrain the AGN polar dust emission given a limited number of photometric data in the MIR regime, which we should keep in mind for the following analysis. We remind that the main purpose of this work is to see how the AGN polar dust would affect the CF of type 1 quasars. Hence, we consider objects with ${\chi }_{\mathrm{CIGALE}}^{2}\lt 3.0$ and ΔBIC > 6 for the correlation analysis (see Section 3.3).

Zoom In Zoom Out Reset image size

Figure 4 shows examples of the optical spectral fitting with QSFit. We confirm that 36,855/37,181 (∼99.1%) objects have reduced ${\chi }_{\mathrm{QSFit}}^{2}\lt 3.0$, while only 170 (∼0.5%) objects are failed to fit owing to low S/N. This indicates that the SDSS spectra of our sample are well fitted by QSFit. Note that objects shown in the top panel of Figure 4 are expected to have a large contribution from the host galaxy to optical spectrum according to the SED fitting (see top panel in Figure 2). We confirm that this is the case for their optical spectra, suggesting that our SED modeling and spectral fitting give consistent results.

Zoom In Zoom Out Reset image size

In addition to goodness of fitting (i.e., ${\chi }_{\mathrm{QSFit}}^{2}$), QSFit provides "quality flags" for measurements to assess the reliability of the results (see Appendix C in Calderone et al. 2017, for more details). We considered Cont_5100_Qual and HB_BR_Qual, which are relevant to the reliability of M_BH and L_bol (see also Table 4). We find that L₅₁₀₀ and FWHM_Hβ are securely estimated for 21,888/37,181 (∼58.9%) objects. To ensure the accuracy of fitting results, we focus on objects with ${\chi }_{\mathrm{QSFit}}^{2}\lt$ 3.0, Cont_5100_Qual = 0, and HB_BR_Qual = 0 for the following analysis (see Section 3.3).

The distributions of our sample with Cont_5100_Qual = 0 and HB_BR_Qual = 0 in the z − L_bol and M_BH − L_bol planes are shown in Figures 5 and 6, respectively. The mean and standard deviation of $\mathrm{log}{L}_{\mathrm{bol}}$ (erg s⁻¹), $\mathrm{log}({M}_{\mathrm{BH}}/{M}_{\odot })$, and $\mathrm{log}\,{\lambda }_{\mathrm{Edd}}$ are 45.0 ± 0.46, 8.29 ± 0.45, and −1.41 ± 0.38, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We investigate the dependences of the CF on L_bol, M_BH, and λ_Edd of our quasar sample. As mentioned in Sections 3.1 and 3.2, quasars with (i) ${\chi }_{\mathrm{CIGALE}}^{2}\lt 3.0$, (ii) ΔBIC > 6, (iii) ${\chi }_{\mathrm{QSFit}}^{2}\lt$ 3.0, (iv) Cont_5100_Qual = 0, and (v) HB_BR_Qual = 0 are used for the following correlation analysis, which yields 5752 objects (see also Figure 1). We discuss how the above selection cuts affect the dependences of the CF in Section 4.2. Their distributions of our subsample in the z − L_bol and M_BH − L_bol planes are also shown in Figures 5 and 6, respectively. We find that our quasar sample selected from the DR16Q covers 43.5 < $\mathrm{log}{L}_{\mathrm{bol}}$ (erg s⁻¹) < 46.4, 6.57 < $\mathrm{log}({M}_{\mathrm{BH}}/{M}_{\odot })$ < 9.89, and −3.16 < $\mathrm{log}{\lambda }_{\mathrm{Edd}}$ < −0.22. The mean and standard deviation of $\mathrm{log}{L}_{\mathrm{bol}}$ (erg s⁻¹), $\mathrm{log}({M}_{\mathrm{BH}}/{M}_{\odot })$, and $\mathrm{log}\,{\lambda }_{\mathrm{Edd}}$ are 45.0 ± 0.39, 8.29 ± 0.46, and −1.42 ± 0.38, respectively. The fact that the DR16Q contains fainter quasars than those in previous releases enables us to investigate the CF for less luminous quasars with less massive BH and smaller accretion rate compared to previous studies baaed on the SDSS quasar catalogs (e.g., Gu 2013; Ma & Wang 2013).

As we cautioned in Section 1, the IR-to-bolometric luminosity ratio would not always be a good tracer of the CF. We hence converted from ${L}_{\mathrm{IR}}^{\mathrm{torus}}/{L}_{\mathrm{bol}}$ to CF by using a nonlinear relation in Stalevski et al. (2016), who provides conversion formulae assuming SKIRTOR with ${R}_{\max }/{R}_{\min }$ = 30 and θ = 0° that are the same as in Table 1 (see also the middle panel of Figure 7 in Stalevski et al. 2016):

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{CF} & = & {a}_{4}{\left(\displaystyle \frac{{L}_{\mathrm{IR}}^{\mathrm{torus}}}{{L}_{\mathrm{bol}}}\right)}^{4}+{a}_{3}{\left(\displaystyle \frac{{L}_{\mathrm{IR}}^{\mathrm{torus}}}{{L}_{\mathrm{bol}}}\right)}^{3}\\ & \,+ & {a}_{2}{\left(\displaystyle \frac{{L}_{\mathrm{IR}}^{\mathrm{torus}}}{{L}_{\mathrm{bol}}}\right)}^{2}+{a}_{1}\left(\displaystyle \frac{{L}_{\mathrm{IR}}^{\mathrm{torus}}}{{L}_{\mathrm{bol}}}\right)+{a}_{0}.\end{array}\end{eqnarray} \tag{ 2 }$$

Since we parameterized τ_9.7 (see Table 1), we chose (a₀, a₁, a₂, a₃, a₄) = (0.192478, 1.40827, −1.48727, 0.875215, −0.177798), (0.195615, 1.20218, −1.04546, 0.47493, −0.0601471), and (0.196387, 1.02696, −0.782418, 0.299937, −0.0255416) for objects with τ_9.7 = 3, 5, and 9, respectively. The uncertainty of the CF was propagated from a relative error of ${L}_{\mathrm{IR}}^{\mathrm{torus}}/{L}_{\mathrm{bol}}$.

Figure 7 shows the relation between the CF of dust torus and AGN properties, i.e., CF–L_bol, CF–M_BH, and CF–λ_Edd for 5752 quasars. We confirm an anticorrelation of the above three quantities, which is consistent with recent works (Ezhikode et al. 2017; Ricci et al. 2017; Zhuang et al. 2018). We also investigate the dependence of the CF of polar dust (CF_PD) on L_bol, M_BH, and λ_Edd, which is shown in Figure 7. We find that CF_PD also clearly depends on AGN properties (see Section 4 for quantitative discussion).

To quantify the anticorrection shown in Figure 7, we conducted a correlation analysis for the subsample of 5752 quasars by using a Bayesian maximum likelihood method provided by Kelly (2007), providing a correlation coefficient (r) that takes into account uncertainty on both the x and y values (see, e.g., Mateos et al. 2015; Toba et al. 2019a). The same analysis was done for the same subsample of 5752 quasars whose CF was derived by the SED fitting without adding the polar dust component in order to check how the presence or absence of polar dust would affect the correlation strength. The resultant correlation coefficients are summarized in Table 2.

Zoom In Zoom Out Reset image size

Figure 8 shows the correlation coefficient (r) of CF–L_bol, CF–M_BH, and CF–λ_Edd with and without adding the AGN polar dust component to the SED fitting. We find that the ∣r∣ values of objects when considering polar dust emission are larger than those without considering the polar dust. In particular, the influence of the polar dust component on the correlation strength may be significant for CF–λ_Edd rather than CF–L_bol and CF–M_BH. This indicates that polar dust wind, probably driven by radiation pressure from the AGN, may be crucial for regulating obscuring dusty structure surrounding SMBHs. Table 2 also provides insight into the properties of polar dust; CF_PD seems to strongly depend on L_bol and λ_Edd rather than M_BH, which is consistent with what we discussed above. Strong radiation pressure from luminous quasars with high λ_Edd may blow out dust in the polar direction, which would be associated with an AGN-driven outflow (e.g., Schartmann et al. 2014). This suggests that CF_PD is regulated by L_bol and possibly λ_Edd.

Zoom In Zoom Out Reset image size

4.2.1. Parent SDSS Quasar Sample

4.2.2. Malmquist Bias

4.2.3. Subsample of Quasars

As described in Sections 3.1 and 3.2, we created a subsample for the correlation analysis by considering the reduced ${\chi }_{\mathrm{CIGALE}}^{2}$, ${\chi }_{\mathrm{QSFit}}^{2}$, ΔBIC, Cont_5100_Qual, and HB_BR_Qual to understand how AGN polar dust affects dust CF. In particular, the threshold of ΔBIC contributed to a decrease in the sample size of quasars from 37,181 to 5752, although the size of the sample is still large enough for statistical investigation. We discuss how the selection criteria for the correlation analysis would be biased toward the resultant correlation coefficient. To address this issue, we performed a Monte Carlo simulation in which we randomly chose a set of threshold values (reduced ${\chi }_{\mathrm{CIGALE}}^{2}$, reduced ${\chi }_{\mathrm{QSFit}}^{2}$, and ΔBIC) and conducted the correlation analysis. We iterated the above process 1000 times where calculation was executed only if the sample size exceeds 1000 in order to ensure reliability of the resultant correlation coefficient.

Figure 9 shows the distribution of r. The mean and standard deviations of r for CF–L_bol, CF–M_BH, and CF–λ_Edd are −0.84 ± 0.10, −0.67 ± 0.13 and −0.26 ± 0.13, respectively, which are roughly consistent with what we obtained in this work (see Table 2). This suggests that selection bias caused by adopting threshold values in this work does not significantly affect the result of the correlation analysis.

Zoom In Zoom Out Reset image size

For the SED and spectral fitting in a consistent manner, we fixed the age of the main stellar population to be 5.0 Gyr for the SED fitting, while we used the template of an elliptical galaxy with a stellar age of 5.0 Gyr for the spectral fitting (see Sections 2.3 and 2.4). Although this assumption (i.e., low-z quasars are hosted by elliptical galaxies) is expected to apply to the majority of our quasar sample (e.g., Bahcall et al. 1997; Dunlop et al. 2003; Kauffmann et al. 2003; Floyd et al. 2004), a fraction of quasars could not be the case (e.g., McLure et al. 1999; Schawinski et al. 2010; Ishino et al. 2020). How would these possible variations of quasar host properties affect resulting correlation coefficients?

To address this issue, we modified (i) the parameter ranges in the age of the main stellar population used for the SED fitting and (ii) AGN host templates for the spectral fitting. For the SED fitting, we allowed the age of the main stellar population to vary from 1.0 to 11.0 Gyr at intervals of 2.0 Gyr. For the spectral fitting, we utilized host galaxy templates of a simulated 2.0 Gyr old elliptical galaxy, S0 galaxy, and three types of spiral galaxies (i.e., Sa, Sb, and Sc; Polletta et al. 2007), in addition to the 5.0 Gyr old elliptical galaxy template we originally used. We randomly chose a host galaxy template from the above options. We then performed the SED fitting and spectral fitting in the same manner as was described in Sections 2.3 and 2.4. We iterated this process 1000 times and estimated the distribution of r.

The resulting distributions of r are shown in Figure 10. We find that the mean and standard deviations of r for CF–L_bol, CF–M_BH, and CF–λ_Edd are −0.79 ± 0.00, −0.62 ± 0.01, and −0.24 ± 0.01, respectively, which are in good agreement with what we obtained in this work (see Table 2). Therefore, we conclude that the influence of a limited number of stellar templates on the correlation coefficients is likely to be small.

Zoom In Zoom Out Reset image size

We reported in Section 4.1 that CF depends on λ_Edd by taking into account polar dust emission. It should be noted that λ_Edd is directly related to L_bol that strongly correlates with the CF, which would induce an artificial correlation of CF–λ_Edd.

To disentangle the dependence and see how r (CF–λ_Edd) would be affected by strong correlation of CF–L_bol, we employed a diagnostic method presented in Hasinger (2008), who investigated redshift dependence on CF based on the relation between the CF and X-ray luminosity (see also Toba et al. 2014). First, we estimated the average value of the slope of the relation between the CF and L_bol in the range of $-2.0\lt \mathrm{log}\,{\lambda }_{\mathrm{Edd}}\lt 0.0$. We then estimated the normalization value at a luminosity of $\mathrm{log}\,{L}_{\mathrm{bol}}$ = 45.0 erg s⁻¹, in the middle of the observed range, as a function of λ_Edd by keeping the slope fixed to the average value. The λ_Edd bins were $-2.0\,\lt \mathrm{log}{\lambda }_{\mathrm{Edd}}\lt -1.5$, $-1.5\,\lt \mathrm{log}{\lambda }_{\mathrm{Edd}}\lt -1.0$, $-1.0\lt \mathrm{log}\,{\lambda }_{\mathrm{Edd}}\,\lt -0,5$, and $-0.5\,\lt \mathrm{log}\,{\lambda }_{\mathrm{Edd}}\lt 0.0$ in the range of $\mathrm{log}\,{L}_{\mathrm{bol}}$ = 44–46 erg s⁻¹. We find that the resulting correlation coefficient is r(CF–λ_Edd) ∼ −0.31, which is consistent with what we obtained in Section 4.1. Therefore, we conclude that a strong CF–L_bol correlation may not significantly affect the correlation coefficient of CF–λ_Edd.

We then check whether or not the derived CF can actually be estimated reliably, given the limited number of photometry and its uncertainty. We execute a mock analysis that is a procedure provided by X-CIGALE (see, e.g., Buat et al. 2012, 2014; Ciesla et al. 2015; Lo Faro et al. 2017; Boquien et al. 2019; Toba et al. 2019b, 2020a, for more details) for 5752 quasars. This analysis is done by a mock catalog in which a value is taken from a Gaussian distribution with the same standard deviation as the observation is added to each photometry originally used, which would also enable us to test how the difference in photometry influences CF (see Section 2.3).

Figure 11 shows the histogram of differences in CFs that are derived from X-CIGALE in this work and from the mock catalog. The mean and standard deviations of ΔCF = 0.04 ± 0.07, suggesting that the CF of the majority of our quasar sample is insensitive to the limited number of photometric points and difference in photometry.

Zoom In Zoom Out Reset image size

In Section 4.1, we report that adding AGN polar dust emission to SED fitting makes anticorrelations of CF–L_bol, CF–M_BH, and CF–λ_Edd stronger than those without adding it. On the other hand, L_bol, M_BH, λ_Edd, and also redshift should be correlated with each other as mentioned in Sections 4.2.2 and 4.4 (see also Figures 5 and 6). With all the caveats discussed in Sections 4.2–4.5 in mind, we discuss how the CF–λ_Edd correlation depends on AGN properties (i.e., L_bol and M_BH) and redshift. The resultant correlation coefficients are summarized in Table 3.

Table 3. Dependence of Correlation Coefficient of CF–λ_Edd on Redshift, BH Mass, and Bolometric Luminosity

	Parameter Range	r (CF–λEdd)
	z < 0.12	−0.49 ± 0.32
Redshift	0.12 < z < 0.2	−0.24 ± 0.14
	z > 0.2	−0.22 ± 0.05
	${\rm{log}}({M}_{{\rm{BH}}}/{M}_{\odot })\lt 8.0$	−0.56 ± 0.07
BH mass	$8.0\lt {\rm{log}}({M}_{{\rm{BH}}}/{M}_{\odot })\lt 9.0$	−0.50 ± 0.05
	${\rm{log}}({M}_{{\rm{BH}}}/{M}_{\odot })\gt 9.0$	−0.52 ± 0.22
	${\rm{log}}{L}_{{\rm{bol}}}\gt 44.8$	−0.10 ± 0.09
Lbol (erg s−1)	$44.8\lt {\rm{log}}{L}_{{\rm{bol}}}\lt 45.0$	−0.12 ± 0.13
	$\mathrm{log}{L}_{\mathrm{bol}}\gt 45.0$	0.06 ± 0.07

Note. The polar dust component is taken into account to estimate r (CF–λ_Edd).

Download table as:  ASCIITypeset image

4.6.1. Redshift Dependence

Although we confirm the anticorrelations of CF and L_bol, M_BH, and λ_Edd, the correlation strength differs among them. We find that the CF–L_bol correlation is stronger than other two correlations as shown in Table 2 and Figure 8, which seems inconsistent with the X-ray-based work of Ricci et al. (2017), who reported that the CF–λ_Edd correlation is stronger than the CF–L_bol correlation. What causes this discrepancy?

There are some possibilities, such as the difference in (i) redshift, (ii) sample selection (i.e., our sample is limited to type 1 quasars, while the sample in Ricci et al. 2017 is not), and (iii) definition of CF (i.e., ${L}_{\mathrm{IR}}^{\mathrm{torus}}/{L}_{\mathrm{bol}}$ is employed in this work, while a fraction of X-ray-obscured sources is employed in Ricci et al. 2017). It is hard to quantify, however, the influence of the difference in the sample selection and definition of the CF on the correlation strength, and it is beyond the scope of this work. We thus argue for the possibility that the difference in redshift could cause the discrepancy. Ricci et al. (2017) targeted nearby AGNs with a median redshift of z ∼ 0.037, whereas we focus on AGNs with much higher redshift up to z = 0.7.

Figure 12 shows the redshift dependence of the CF–λ_Edd relation. We find that the CF of AGNs with z < 0.12 (which is similar to that of the sample in Ricci et al. 2017) tends to strongly depend on λ_Edd, with a correlation coefficient of r = − 0.49 ± 0.32, while CF–λ_Edd for AGNs at z > 0.2 shows weaker correlations than that at z < 0.12. This indicates that the aforementioned discrepancy can be explained by the difference in redshift between the samples in Ricci et al. (2017) and this work. We also remind that our correlation analysis would be affected by the Malmquist bias, which could make the correlation strength weaker, as discussed in Section 4.2.2.

Zoom In Zoom Out Reset image size

This result also might suggest that the CF decreases with increasing redshift for at least z < 0.7, a trend that is consistent with what is reported in Toba et al. (2014), although there is large uncertainty. On the other hand, the fraction of obscured X-ray AGNs increases with redshift for the same luminosity (e.g., Ueda et al. 2014). This discrepancy might be caused by the difference in sample selection and definition of CF.

4.6.2. BH Mass Dependence

Recently, Kawakatu et al. (2020) constructed a model of a nuclear starburst disk supported by the turbulent pressure from Type II supernovae and investigated how the CF depends on AGN properties by taking account of anisotropic radiation pressure from AGNs. They found that the CF strongly depends on λ_Edd for AGNs with M_BH < 10⁸ M_⊙, while the dependence of the CF on λ_Edd is much weaker for AGNs with M_BH > 10⁹−10¹⁰ M_⊙. Although they did not incorporate the AGN polar dust outflow, it is worth checking the dependence of the CF–λ_Edd correlation on M_BH.

Figure 13 shows how the difference in M_BH would affect the correlation strength of CF–λ_Edd. The resultant r of CF–λ_Edd for quasars with ${\rm{log}}({M}_{{\rm{BH}}}/{M}_{\odot })\lt 8.0$, $8.0\lt {\rm{log}}({M}_{{\rm{BH}}}/{M}_{\odot })\lt 9.0$, and ${\rm{log}}({M}_{{\rm{BH}}}/{M}_{\odot })\gt 9.0$ are −0.56 ± 0.07, −0.50 ± 0.05, and −0.52 ± 0.22, respectively. This result indicates that the CF of quasars with < 10⁸ M_⊙ strongly depends on λ_Edd. The overall statistical trend is in good agreement with what was reported in Kawakatu et al. (2020).

Zoom In Zoom Out Reset image size

4.6.3. Bolometric Luminosity Dependence

Figure 14 shows the dependence of L_bol on CF–λ_Edd. We find that the anticorrelation of CF–λ_Edd may disappear for AGNs with $\mathrm{log}\,{L}_{\mathrm{bol}}\gt 45.0$. In this L_bol range, the correlation coefficients are even positive with large uncertainty (see Table 3). Zhuang et al. (2018) reported that the CF of luminous quasars with L_bol ≳ 10^45.0 erg s⁻¹ increases with increasing λ_Edd, although this positive correlation would appear when $\mathrm{log}{\lambda }_{\mathrm{Edd}}$ ranges from −0.25 and 0.5. This result could indicate that the overall trend of CF–λ_Edd reported in Section 4.1 is determined by AGNs with log L_bol < 45.0.

Zoom In Zoom Out Reset image size