Modeling of the Quasar Main Sequence in the Optical Plane

We parameterize the incident continuum as two power laws with low- and high-energy cutoffs. The optical/UV power law represents the emission from the accretion disk, where most of the energy is dissipated. This component forms the big blue bump (Czerny & Elvis 1987; Richards et al. 2006). The slope of the power law α_UV is assumed to be 1/3 (in convention ${F}_{\nu }\propto {\nu }^{\alpha }$, consistent with the theory of accretion disks; Shakura & Sunyaev 1973), supported by polarimetric observations of quasars (Kishimoto et al. 2008) and broadband data fitting (Capellupo et al. 2015). The high-frequency cutoff, ν_*, is the basic parameter that we vary in our model, i.e., we assume that the accretion disk luminosity can be described by the formula

$$\begin{eqnarray}&&\nu {L}_{\nu }={\rm{A}}{\nu }^{4/3}\exp (-\nu /{\nu }_{* }).\end{eqnarray} \tag{ 1 }$$

This parameter is directly related to the maximum temperature in the accretion disk, T_max. If we use the Newtonian approximation of the disk from Shakura & Sunyaev (1973), there is a simple relation between the maximum temperature and global parameters of the accretion flow

$$\begin{eqnarray}&&{T}_{\max }=1.732\times {10}^{19}\ {\left(\displaystyle \frac{\dot{M}}{{M}^{2}}\right)}^{1/4}\,[{\rm{K}}],\end{eqnarray} \tag{ 2 }$$

where the black hole mass, M, and the accretion rate, $\dot{M}$, are expressed in cgs units (see, e.g., Panda et al. 2017). The disk maximum temperature is also related to the peak of the spectral energy distribution (SED) of a Shakura–Sunyaev disk on the νF_ν plot

$$\begin{eqnarray}&&\displaystyle \frac{h{\nu }_{\max }}{{{kT}}_{\max }}=2.464.\end{eqnarray} \tag{ 3 }$$

Here h is the Planck constant and k is the Boltzmann constant. The simplified shape of the spectrum given by Equation (1) peaks at the frequency 4ν_*/3. Therefore, combining these two factors, we finally obtain the convenient parameterization of the disk shape in the form

$$\begin{eqnarray}&&\nu {L}_{\nu }={\rm{A}}{\nu }^{4/3}\exp (-h\nu /1.853{{kT}}_{\max }).\end{eqnarray} \tag{ 4 }$$

In the case of an accretion disk around a rotating black hole, the relativistic corrections are very important, the position of the innermost stable circular orbit strongly depends on spin, and the relations between the T_max and the global parameters of the accretion flow are more complex, and in that case more advanced back-modeling is necessary, as done by Thomas et al. (2016). However, for the purpose of the current work we only consider nonrotating black holes, which are well parameterized by the proposed prescription (see Figure 1) and are characterized by only two free parameters: black hole mass (M) and accretion rate ($\dot{M}$).

Zoom In Zoom Out Reset image size

The normalization of the optical/UV power law in Equation (1) is given by the values of the black hole mass and accretion rate

$$\begin{eqnarray}&&A=1.33\times {10}^{-20}{(\dot{M}M)}^{2/3}\,[\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{Hz}}^{-4/3}]\end{eqnarray} \tag{ 5 }$$

under the assumption of the viewing angle 60°. Adopting such a viewing angle relates the integrated spectrum to the source bolometric luminosity. Small corrections due to an on average smaller viewing angle of ∼40°, due to the obscuration by the dusty torus (Lawrence & Elvis 2010), are not an important factor.

AGNs emit part of the energy in the form of the hard X-ray emission. We parameterize this component again as a power law with a low- and high-energy cutoff. The high-energy cutoff is set at 100 keV as in most fitting packages (e.g., OPTXAGN—standard AGN shape in CLOUDY; Done et al. 2012; Thomas et al. 2016). Observational data show a dispersion in this quantity from about 50 keV to over 200 keV, but the measurements are still rare (for a compilation of measurements from NuSTAR, see Fabian et al. 2015).

The relative normalization of the X-ray component with respect to UV is found from the universal scaling law recently discovered by Lusso & Risaliti (2017),

$$\begin{eqnarray}&&\mathrm{log}\ {L}_{{\rm{X}}}=0.610\ \mathrm{log}{L}_{\mathrm{UV}}+0.538\ \mathrm{log}\ {v}_{\mathrm{FWHM}}+3.40,\end{eqnarray} \tag{ 6 }$$

where L_UV is a monochromatic luminosity νL_ν measured at 2500 Å and L_X is measured at 2 keV. The importance of the X-ray component decreases with the increase of the UV flux, as illustrated in Figure 1. Equation (6) gives us the value of the broadband index α_ox measured between 2500 Å and 2 keV. In Figure 2 we show the range of this index covered in our computations.

Zoom In Zoom Out Reset image size

However, in order to use this formula, we actually need the line v_FWHM, which depends on the cloud location. We need this distance also for calculation of the line emissivity.

To determine the distance to the BLR clouds, we use the observationally established relation from Bentz et al. (2013). We choose the version Clean from their Table 14,

$$\begin{eqnarray}&&\mathrm{log}\ {R}_{\mathrm{BLR}}=1.555+0.542\ \mathrm{log}\ {L}_{\mathrm{44,5100}}\ \ [\mathrm{light}\,\mathrm{days}],\end{eqnarray} \tag{ 7 }$$

where the luminosity at 5100 Å is measured in units of 10⁴⁴ ergs s⁻¹.

Thus, the shape of the incident radiation with its normalization provides us with the distance to the BLR. The value of the black hole mass gives the required line width

$$\begin{eqnarray}&&{v}_{\mathrm{FWHM}}={\left(\displaystyle \frac{\mathrm{GM}}{{{\rm{R}}}_{\mathrm{BLR}}}\right)}^{1/2},\end{eqnarray} \tag{ 8 }$$

where we adopt the virial factor 1 for the BLR clouds. This is indeed a simplified approach since the virial factor seems to be a function of the measured line width but the average value is actually close to 1 (Mejía-Restrepo et al. 2018), which is enough for our current purpose.

Thus, for the two parameters characterizing the SED (e.g., normalization and the maximum disk temperature, or, equivalently, M and $\dot{M}$) all the other parameters required to do the radiative transfer calculation are uniquely determined.

We calculate Hβ and optical Fe ii emission performing the computations with the CLOUDY code (Ferland et al. 2017), version 17. We use the incident radiation and the distance to the BLR as described above. We adopt a traditional single-cloud approximation (e.g., Mushotzky & Ferland 1984). This is a reasonable approach since we are interested in a single line ratio, so going to the more complex locally optimized cloud (LOC) model of Baldwin et al. (1995a) is not necessary. In addition, the BLR extension is not large, and the outer-to-inner ratio is estimated to be of order of 4–5 (Koshida et al. 2014). In CLOUDY computations we assume the plane-parallel geometry, which provides emission both from the illuminated side of the cloud and from the dark side of the cloud, in equal proportions. We assume that the density inside the cloud is constant. This is again an approximation since the clouds in the vicinity of the nucleus are expected to be in pressure equilibrium in order to survive at least in the dynamical timescale (e.g., Różańska et al. 2006). Some models based on radiation pressure confinement suggested that the density gradient in the cloud must be steep (Stern et al. 2014; Baskin & Laor 2018). However, at least in a significant fraction of sources that show the presence of the intermediate-line region instead of a well-separated BLR and NLR, the local density at the cloud surface must be high, and the density gradient is rather shallow (Adhikari et al. 2016, 2018a). We assume two additional free parameters of the cloud: density, n, and hydrogen column density, N_H.

Hβ flux is taken from the code output, and Fe ii emission is calculated by summing up all the Fe ii lines in the range from 4434 to 4684 Å, as defined in Shen & Ho (2014). The parameter R_Fe is calculated as the ratio of these two numbers.

For convenience, we decided to use the following parameters to present our results: T_max, L/L_Edd, n, and N_H. Here L/L_Edd is defined using the Eddington accretion rate value 1.26 ×10³⁸ (M/M_⊙), and to calculate the Eddington accretion rate, we used the Newtonian accretion disk efficiency 1/12 (Shakura & Sunyaev 1973). For the Eddington ratio, we choose values between 0.01 and 1. We cannot consider higher values since then the slim disk effects would be important (Abramowicz et al. 1988; Wang et al. 2014). At lower accretion rates inner optically thin flow likely replaces the cold geometrically thin optically thick disk, which leads to modification of both the SED (Narayan & Yi 1994; Nemmen et al. 2014) and the BLR itself (e.g., Czerny et al. 2004; Balmaverde & Capetti 2015). In general, we consider the range of the T_max between 1.6 × 10⁴ and 5 × 10⁵ K, appropriate for the AGN big blue bump. The corresponding range of the black hole masses depends on the adopted Eddington ratio. We give example values in Table 1. We assume that the viable range of masses is between 10⁶ M_⊙ and 3 × 10⁹ M_⊙ for low-redshift sources with a detected Hβ line. Thus, not all of the T_max range is realistic for a given choice of the Eddington ratio, and we include this effect in our plots. Therefore, the allowed range of the T_max for a given Eddington ratio is further constrained by the limits on the black hole mass to be between 10⁶ M_⊙ and 5 × 10⁹ M_⊙. For the cloud density n, we choose values from 10¹⁰ to 10¹² cm⁻³, appropriate for the LIL part of the BLR. The column density was assumed to vary from 10²² to 10²⁴ cm⁻². We assume a constant-density cloud profile, but we address the issue of the constant-pressure clouds in the discussion. We also allow for the turbulent velocity since the need for a velocity of order of 10–20 km s⁻¹ is evident from the previous studies of the Fe ii production (Bruhweiler & Verner 2008).

Table 1.  Examples of the Correspondence between the Black Hole Masses and the Assumed T_max and λ_Edd (see Equation (2))

λEdd	0.01	0.01	0.1	0.1	1.0	1.0
$\mathrm{log}\ {T}_{\max }$	5.145	4.270	5.395	4.520	5.645	4.770
$\mathrm{log}\ M/{M}_{\odot }$	6.0	9.5	6.0	9.5	6.0	9.5

Download table as:  ASCIITypeset image

Finally, AGNs can have a range of metallicities, and most studies have found that the abundances are actually at least solar and mostly supersolar (factor of 1–10; Hamann & Ferland 1999; Tortosa et al. 2018), even for high-redshift quasars (by a factor of 5–10; e.g., Simon & Hamann 2010). In our computations we either assumed the standard chemical composition, which is the default in CLOUDY, or we allowed for an increase in Fe abundance. In the first option the default values are provided by CLOUDY, and in this case (see Table 7.1 of CLOUDY manual) C and O abundances come from photospheric abundances of Allende Prieto et al. (2001, 2002), while N, Ne, Mg, Si, and Fe are from Holweger (2001), and the Fe is taken from Holweger (2001). If we now specify a solar abundance, Fe is taken from the GASS model (Grevesse et al. 2010), and in this option we also consider supersolar abundance.

The emitting medium is most likely turbulent, and the turbulence decreases the optical depth of the clouds for lines. Fe ii emissivity is quite sensitive to this value (e.g., Bruhweiler & Verner 2008). Therefore, we performed tests of the influence of the turbulence velocity on the calculated R_Fe, varying it between 0 and 100 km s⁻¹. Overall, the Fe ii emissivity increases, but the trend is not monotonic. The emissivity rises with the rise of the turbulent velocity up to 20 km s⁻¹, and a further increase in the velocity leads to a decrease of Fe ii production apart from the low-temperature tail (see Figure 5). Thus, the turbulence in the range 10–20 km s⁻¹ is generally favored for more efficient Fe ii production, as such values are actually favored in detailed fitting of specific objects (e.g., Bruhweiler & Verner 2008; Marziani et al. 2013; Hryniewicz et al. 2014; Średzińska et al. 2017).

In the computations above we assumed a constant-density model, traditionally adopted in the computations of the BLR clouds (Davidson 1977). However, physically it is not justified, particularly for such thick clouds since they are irradiated from the side exposed to the radiation flux from the central region, while the other side of the cloud is relatively cold. Such a cloud cannot be in hydrostatic equilibrium and would be rapidly destroyed. Therefore, a more appropriate description of the cloud structure is to assume a constant pressure throughout the cloud, as discussed for narrow-line region clouds (e.g., Davidson 1972), warm absorbers (e.g., Różańska et al. 2006; Adhikari et al. 2015), and BLR clouds (e.g., Baskin et al. 2014a). The effect is particularly strong for low local density clouds, but for high-density clouds at BLR distances the compression is relatively less important, by a factor of a few (Adhikari et al. 2018a). Nevertheless, the effect may be noticeable. We thus calculated an example of the constant-pressure cloud corresponding to the most typical values for the observational sample: $\mathrm{log}\ {{\rm{T}}}_{\max }=4.8$, λ_Edd = 0.1, $\mathrm{log}n=12$ at the cloud illuminated surface, and $\mathrm{log}\ {{\rm{N}}}_{{\rm{H}}}=24$. The results are shown in Table 2. The value of R_Fe calculated for such a cloud was somewhat higher than for a constant-density cloud, if the turbulent velocity was neglected, and the effect decreased with the rise of the turbulence velocity. Thus, for such dense clouds, constant-density and constant-pressure models give very similar results.

Clouds are irradiated from one side, and the dark sides of the clouds have a different proportion in Hβ and Fe ii emission. In the computations shown throughout this paper we used a plane-parallel geometry, and the line intensity was calculated from the Intrinsic line intensities section of the main CLOUDY output. These include the combined emission from the dark and the bright sides of the cloud. This approach is justified if an observer is not highly inclined, clouds do not shield each other, the reprocessing of the emission of one cloud by the other cloud is negligible, and the geometrical thickness of the BLR is relatively small. With these approximations, we see on average the same total illuminated and dark surfaces of all clouds. However, if any of these assumptions are violated, the obtained R_Fe ratio will be different. Such extreme setups might be responsible for the extreme Fe ii emitters that were recovered in Section 3 only for supersolar metallicity.

Thus, another possibility is that the abundances are always solar but the BLR is so geometrically thick that it covers most of the quasar sky, and the number of clouds so large that the reprocessing is important. To check this option, we calculated one cloud model using the closed geometry. We assumed the same cloud parameters as in Section 2. In this case the increase of the R_Fe is not very large since two effects counteract. One is that we now see basically the dark sides of the clouds, but the other is that multiple scattering increases the local incident flux and the cloud ionization. In addition, in order not to heavily absorb the continuum, we need a gap in the cloud distribution just along the line of sight to the innermost part of the accretion disk.

The second possibility is similar to the one above, but with the cloud number not as high, so the cross-illumination of clouds can be neglected. Again, in this case we see more of the dark sides of the clouds than of the bright sides. Such a picture has been used by Ferland et al. (2009), where the Fe ii emission has been modeled as coming from infalling clouds. In this case we calculated separately the R_Fe values for the bright side and for the dark side of the cloud, by using a plane-parallel approximation, but with the sphere command to store just the outward line emission, and the inward emission has been calculated as a difference between total (intrinsic) and outward line flux:

$$\begin{eqnarray}&&{R}_{\mathrm{Fe}{\rm{}}\mathrm{II}}(\mathrm{dark})=\displaystyle \frac{\mathrm{Fe}\,{\rm{}}{\mathrm{II}}_{\mathrm{Intrinsic}}-\mathrm{Fe}\,{\rm{}}{\mathrm{II}}_{\mathrm{Bright}}}{{\rm{H}}{\beta }_{\mathrm{Intrinsic}}-{\rm{H}}{\beta }_{\mathrm{Bright}}}.\end{eqnarray} \tag{ 9 }$$

The results are given in Table 3. In the case of the same cloud (log T = 5, log n = 12, $\mathrm{log}\ {N}_{{\rm{H}}}$ = 24), we obtain R_Fe = 0.2883 for the bright side and R_Fe = 0.2782 for the dark side. Thus, if the clouds cover the nucleus densely, we do not see an enhancement in the Fe ii production since the intercloud scattering increases the overall ionization. The situation is different if we calculate emission from the dark sides of the clouds in a standard plane geometry, when no such intercloud scattering is present. In such a geometry, the Fe ii emission from the dark side of the cloud measured with respect to Hβ is enhanced by a factor of 6. So, in rare cases, when we see predominantly the dark sides of the clouds, we can reproduce R_Fe values as high as a few, required to explain the extreme data points without postulating supersolar metallicity.

The proper coverage of the optical plane, including the right corner occupied by the extreme Fe ii emitters, required allowing for supersolar abundances, although the average quasar parameters were well represented without (see Section 3). The increase of the abundances simply enhances the Fe ii emissivity, shifting the theoretical curves to the right in the optical plane (see Figure 4). The same effect is seen if the parameter R_Fe is studied directly as a function of the maximum disk temperature. We show the corresponding plots in Figure 7.

Zoom In Zoom Out Reset image size

Upon using the GASS model (at Z_⊙), we find an increase in the Fe ii strength by 7%–9% in all three cases of changing microturbulence compared to the default case for solar metallicity. The Fe ii strength increases by 95%–118% and by 237%–406% when the metallicity (Z) is increased to 3 and 10 Z_⊙, respectively (Table 4). In this case only the most extreme values of R_Fe (above 4) are not accounted for and might require that we see predominantly dark sides of the clouds (see Figure 7).

In Śniegowska et al. (2018) 27 objects were selected for study from the Shen et al. (2011a) catalog, but after careful analysis only six objects were confirmed as strong Fe ii emitters. Three of these sources had broad Hβ lines, with FWHM above 4500 km s⁻¹, mean black hole mass $\mathrm{log}M=8.9$, and Eddington ratio of about 0.01, while the other three had very narrow Hβ (below 2100 km s⁻¹), mean black hole mass $\mathrm{log}M=7.5$, and Eddington ratio above 0.3. The first family of quasars is consistent with high expected values of R_Fe since the typical maximum temperature in this case is about 20,000 K. The second group has the temperatures of the order of 2 × 10⁵ K, and from the model computations the expected values of R_Fe are low, particularly for high Eddington ratio sources. The model predicts only the further rise of R_Fe if the temperatures are well above 10⁶ K. This cannot happen within the frame of the blackbody representation of the big blue bump.

The inner radius of the BLR cloud that has been used in this paper follows Bentz et al. (2013) (Section 2.2). However, the BLR is actually extended. Here we test the effect of changing the radius from 0.3R_Bentz to 5R_Bentz. The lower limit (0.3R_Bentz) used has been set corresponding to maximum disk temperature (∼2000 K) expected from the BLR model based on dust presence in the disk atmosphere (Czerny & Hryniewicz 2011b). The upper limit is set assuming that R_out ∼ 5R_Bentz comes from the dust reverberation studies of the torus (Koshida et al. 2014). The results for one such case are shown in Figure 8. There is a monotonic behavior of the Fe ii strength with respect to changing BLR radius. In a single-zone approximation Fe ii emission is relatively more efficient if the BLR is located closer in, with all the other parameters fixed. This suggests that future studies should include the radial stratification of the BLR, but it is not simple since the results would depend on the weighted emissivity as a function of radius. Additionally, continuum is variable and BLR responds to it after a delay. One needs to take present-day continuum luminosity and line width that traces radius corresponding to continuum luminosity from the past. Single-epoch spectra, as in the Sloan Digital Sky Survey, do not trace this effect, which introduces some bias. Study of the full BLR structure is beyond the scope of the current work.

Zoom In Zoom Out Reset image size

In our model we did not include the range of viewing angles toward the nucleus. As pointed out by many authors (Zhang & Wu 2002; Collin et al. 2006; Shen & Ho 2014), BLR is not flat. If the emission of Fe ii and Hβ comes roughly from the same region, as assumed in the current paper, the ratio R_Fe is not affected, but the measured Hβ width is expected to depend on the viewing angle, i. In type 1 sources this viewing angle is never very large; otherwise, the torus shields the view of the nucleus. The frequently adopted range of the viewing angles is thus between 0° and 45°. Marziani et al. (2001) and others have tried to connect the observational plane with the source orientation and Eddington ratio. The line width, in turn, depends on both $\sin i$ and the turbulent (random) velocity field. Thus, the FWHM can be affected by a factor of 2 for the BLR thickness of order of 0.3 (see Equation (8) in Collin et al. 2006), and a similar effect was determined in the studies of the virial factor trends (Mejía-Restrepo et al. 2018). Thus, the vertical extension of the covered area can be increased by a factor of less than 2 if this effect is included. On the other hand, we see from our modeling that the spread in the vertical direction is mostly caused by the coupling between the model parameters and the Fe ii production efficiency. We intend to look into the effects of the orientation on the main sequence more carefully in our subsequent work.

In our approach we have not yet explicitly connected the Eddington ratios to the black hole masses in the observed sample. If we plot the quality-controlled sample of 4989 objects from the original Shen et al. (2011b) study, the Eddington ratio and the black hole masses in the current optical plane can be constrained with

$$\begin{eqnarray}&&\mathrm{log}{\lambda }_{\mathrm{Edd}}=-1.05\,\mathrm{log}{M}_{\mathrm{BH}}+7.15,\end{eqnarray} \tag{ 10 }$$

where M_BH is considered in M_⊙. This is shown in the left panel of Figure 9. If the analysis to model the optical plane is to be contained within the limits of this boundary defined between the Eddington ratio and the black hole mass, then we cover a much lesser portion of the optical plane than before (see right panel of Figure 9). Still, with this additional constraint, we cover 4903 out of the 4989 objects, i.e., over 98% of the total subsample. Indeed, as shown in Figure 9, we require deeper observations to exploit the lower regime of the $\mathrm{log}\ {\lambda }_{\mathrm{Edd}}-\mathrm{log}\ {M}_{\mathrm{BH}}$ in the context of the quasar optical plane.

Zoom In Zoom Out Reset image size