Probing the rapid formation of black holes and their galaxy hosts in QSOs

TL;DR

This study models the spectral energy distributions of 1,359 QSOs to investigate their host galaxy formation, star formation rates, and black hole growth, revealing rapid galaxy formation and a phase of accelerated SMBH growth.

Contribution

It provides new insights into the rapid formation of QSO host galaxies and the evolutionary phase where SMBH growth outpaces galaxy growth.

Findings

Star formation histories are short, indicating rapid galaxy formation.

No evidence of AGN quenching star formation.

QSOs show a phase of SMBH growth exceeding galaxy growth.

Abstract

Using the modelling code X-CIGALE, we reproduced the SEDs of 1,359 SDSS QSOs within the redshift range 0 < z < 4, for which we have NIR/MIR fluxes with the highest quality and spectral data characterizing their SMBHs. Consistent with a rapid formation of the host galaxies, the star formation histories (SFHs) have small e-folding, at most 750 Myrs using an SFH function for Spiral or 1000 Myrs using one for Elliptical. Above z \sim 1.6, the two solutions are degenerate, the SEDs being dominated by the AGN continuum and high star formation rates (SFRs), typical of starburst galaxies, while at lower redshifts the starburst nature of the host, independent from its morphology, is better reproduced by an Spiral SFH. In general, the SFR increases with the redshift, the mass of the bulge, the AGN luminosity and Eddington ratio, suggesting there is no evidence of AGN quenching of star formation. Comparing the specific BHAR with specific SFR, all the QSOs at any redshift trace a linear sequence below the Eddington luminosity, in parallel and above the one-to-one relation, implying that QSOs are in a special phase of evolution during which the growth in mass of their SMBH is more rapid than the growth in mass of their galaxy hosts. This particular phase is consistent with a scenario where the galaxy hosts of QSOs in the past grew in mass more rapidly than their SMBHs, suggesting that a high star formation efficiency during their formation was responsible in limiting their masses

Tables

Table 1. Table 1: Median properties of HQWISE QSOs in 14 redshift bins. The numbers in col. 3 give the number of QSOs with BH masses measured using the Mg  ii or C  iv emission lines, compared with the total number of HQWISE QSOs in each bin. Note that when both lines are present in a spectrum, we only used the mass determined using the Mg  ii line.

(1)	(2)	(3)	(4)	(5)	(6)
#	z bin	N	Log MBH	Log LAGN	$nEdd$
			(M⊙)	(erg s-1)
0	$\left(0.00-0.25\right]$	0/103	–	–	–
1	$\left(0.25-0.50\right]$	103/188	7.53	44.12	-1.17
2	$\left(0.50-0.75\right]$	120/135	8.36	45.02	-0.75
3	$\left(0.75-1.00\right]$	211/226	8.97	45.71	-0.70
4	$\left(1.00-1.25\right]$	171/177	9.25	45.98	-0.72
5	$\left(1.25-1.50\right]$	121/124	9.39	46.22	-0.58
6	$\left(1.50-1.75\right]$	109/112	9.48	46.38	-0.52
7	$\left(1.75-2.00\right]$	113/118	9.60	46.50	-0.58
8	$\left(2.00-2.25\right]$	69/70	9.58	46.59	-0.44
9	$\left(2.25-2.50\right]$	21/31	9.43	46.66	-0.56
10	$\left(2.50-2.75\right]$	19/23	9.89	47.00	-0.47
11	$\left(2.75-3.00\right]$	21/27	9.51	46.83	-0.53
12	$\left(3.00-3.25\right]$	11/13	9.06	46.94	-0.61
13	$\left(3.25-4.00\right]$	9/12	9.90	47.07	-0.67

Table 2. Table 2: Summary of inputs used in X-Cigale for ultimate run with polar emission.

Parameters	Values	Descriptions
Star formation history (SFH) functions: sfhdelayed and sfh2exp
${\tau }_{\text{main}}$	50, 100, 500, 1000, 2000, 4000, 7000	e-folding time in Myr.
Age	CT	Age of the oldest stars (Myr); we used the cosmic time (CT) at the redshift of the QSOs
fritz06
R${}_{\text{max}}$/R${}_{\text{min}}$	10.0, 30.0, 60.0, 100.0, 150.0	Ratio of the maximum to minimum radii of the dust torus.
$\tau$	0.1, 0.6, 1.0, 6.0, 10.0	Optical depth at 9.7 $\mathrm{\mu m}$.
$\beta$	$-$1.00, $-$0.75, $-$0.50, $-$0.25, 0.00	Beta describing the dust density distribution with r the radius.
$\gamma$	0.0, 2.0	Gamma from the power-law density distribution for the polar component of the dust torus.
Opening Angle ($\theta$)	60.0, 100.0, 140.0	Full opening angle of the dust torus.
$\psi$	(20.000, 25.000, 30.000), 89.990	Angle between equatorial axis and line of sight; values within parenthesis were tried in different run omitting value out of parenthesis.
$\delta$	-1, -0.5, -0.7, 0	Power-law index $\delta$ modifying optical slope of AGN continuum. Negative values make slope steeper/positive values make slope shallower.
$f_{AGN}$	0.1, 0.5, 0.9	Fraction of AGN contribution to the IR luminosity.
Extinction Law	1	Extinction law of the polar dust: 0 (SMC), 1 (Calzetti 2000), or 2 (Gaskell et al. 2004).
$E\left(B-V\right)$	0, 0.05, 0.1, 0.15, 0.2, 0.3, 0.4	E(B-V) for the extinction in the polar direction in magnitudes.
Temperature	100.0, 500.0, 1000.0	Temperature of the polar dust in K.
Dust emission: dl2014
$q_{\mathrm{PHA}}$	0.45 (0.47, 2.50, 4.58, 5.95, 7.32)	Mass fraction of PAH (in percentage); values within parenthesis were tried in different run omitting value out of parenthesis.
umin	50	Minimum radiation field.
$\alpha$	1.0	Powerlaw slope $dU/dM\propto U^{\alpha }$
$\gamma$	0.1	Fraction illuminated from Umin to Umax.
Single-age stellar population (SSP): bc03
$Z$	(0.0004, 0.004, 0.008), 0.02, 0.05	Metallicity; values within parenthesis were tried in different runs omitting values out of parenthesis.
Separation Age	10	Age difference between the youngest and the oldest stellar populations (Myr).
Dust attenuation: dustatt_modified_starburst
E${\left(B-V\right)}_{\mathrm{line}}$	0.005, 0.05, 0.5, 0.10, 0.70	The colour excess of the nebular lines light for both the young and old population.
E${\left(B-V\right)}_{\mathrm{factor}}$	0.25, 0.5, 0.75	Reduction factor E${\left(B-V\right)}_{line}$ of stellar continuum attenuation reduction.
Power-law slope	0.0, 0.25, 0.5	Slope of the power law modifying the attenuation curve.

Table 3. Table 3: Best fits for sfhdelayed function; the values reported are based on the medians of the models with χ 2 < 1 superscript 𝜒 2 1 \chi^{2}<1 (61% of the sample).

(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)	(15)
bin	z	CT	${\chi }^2$	f(${\chi }^2>3$)	W1-W2	W2-W3	O.A.	$\delta$	fAGN	EBV	EBV	Log(SFR)	${\tau }_{main}$	Tdust
		Gyrs		%			(deg)			Lines	Factor	(M${}_{\odot }{}^{}\mathrm{yr}_{}^{-1}$)	(Myr)	(K)
0	0.18	11.2	0.64	19	0.96	2.96	60	-0.7	0.9	0.70	0.75	-1.22	1000	100
1	0.38	9.3	0.57	26	1.02	3.07	60	-1.0	0.9	0.70	0.75	-0.82	1000	100
2	0.62	7.3	0.72	27	1.04	2.94	60	-1.0	0.9	0.70	0.75	0.17	1000	100
3	0.89	6.3	0.42	10	1.17	3.43	140	-0.7	0.9	0.70	0.75	0.85	500	100
4	1.11	5.3	0.33	8	1.31	3.00	100	-0.7	0.9	0.70	0.75	1.89	500	100
5	1.38	4.5	0.29	7	1.38	3.08	140	-0.7	0.9	0.70	0.75	1.69	500	100
6	1.62	3.9	0.18	6	1.37	3.25	140	-0.5	0.9	0.70	0.75	2.08	1000	100
7	1.86	3.4	0.24	5	1.30	3.38	140	-0.7	0.9	0.50	0.50	2.52	1000	500
8	2.12	3.0	0.28	6	1.23	3.53	100	-0.5	0.9	0.50	0.50	2.76	1000	500
9	2.38	2.7	0.14	10	1.18	3.61	140	-0.6	0.9	0.50	0.50	3.02	1000	500
10	2.63	2.4	0.16	9	1.14	3.65	140	-0.5	0.9	0.50	0.50	3.33	1000	500
11	2.85	2.2	0.42	11	1.08	3.58	140	-0.7	0.9	0.10	0.25	3.39	500	500
12	3.09	2.0	0.32	15	0.98	3.66	140	-0.7	0.9	0.10	0.50	3.53	500	100
13	3.41	1.7	0.44	8	1.03	3.73	140	-0.5	0.9	0.10	0.50	3.54	500	100

Table 4. Table 4: Best fits for sfh2exp function; the values reported are based on the medians of the models with χ 2 < 1 superscript 𝜒 2 1 \chi^{2}<1 (61% of the sample).

(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)	(15)
bin	z	CT	${\chi }^2$	f(${\chi }^2>3$)	W1-W2	W2-W3	O.A.	$\delta$	fAGN	EBV	EBV	Log(SFR)	${\tau }_{main}$	Tdust
		Gyrs		%			(deg)			Lines	Factor	(M${}_{\odot }{}^{}\mathrm{yr}_{}^{-1}$)	(Myr)	(K)
0	0.18	11.2	0.56	27	0.98	3.01	80	-0.7	0.9	0.50	0.50	0.22	750	100
1	0.40	9.3	0.62	32	1.02	3.05	60	-1.0	0.9	0.70	0.50	0.91	50	100
2	0.66	7.3	0.70	28	1.07	2.94	60	-1.0	0.9	0.70	0.75	1.48	50	100
3	0.88	6.3	0.41	14	1.22	2.97	140	-1.0	0.9	0.70	0.75	1.72	50	100
4	1.11	5.3	0.33	12	1.34	2.99	140	-0.7	0.9	0.70	0.75	1.86	50	100
5	1.38	4.5	0.30	9	1.39	3.09	140	-0.7	0.9	0.70	0.75	2.15	100	100
6	1.62	3.9	0.17	8	1.37	3.25	140	-0.5	0.9	0.70	0.50	2.35	500	100
7	1.86	3.4	0.17	6	1.27	3.41	140	-0.5	0.9	0.50	0.50	2.70	500	100
8	2.12	3.0	0.23	7	1.20	3.50	140	-0.3	0.9	0.50	0.50	2.84	500	500
9	2.38	2.7	0.10	16	1.17	3.64	100	-0.5	0.9	0.50	0.50	2.98	500	500
10	2.63	2.4	0.12	9	1.12	3.66	140	-0.5	0.9	0.10	0.50	3.19	50	500
11	2.85	2.2	0.39	30	1.08	3.60	140	-0.5	0.9	0.50	0.25	3.26	500	500
12	3.10	2.0	0.75	31	1.04	2.82	100	0.0	0.9	0.10	0.50	3.59	500	100
13	3.41	1.7	0.58	8	1.07	3.71	140	0.0	0.9	0.10	0.75	3.23	500	100

Table 5. Table 5: Best fitted results for sfh2exp with χ 2 > 3 superscript 𝜒 2 3 \chi^{2}>3 varying the psy angle

(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)	(15)	(16)
bin	z	N	${\chi }^2$	${\chi }^2$	$\psi$	W1-W2	W2-W3	O.A.	$\delta$	fAGN	EBV	EBV	Log(SFR)	${\tau }_{\mathrm{main}}$	Tdust
			previous	new	(deg)			(deg)			Lines	Factor	(M${}_{\odot }{}^{}\mathrm{yr}_{}^{-1}$)	(Myr)	(K)
0	0.19	28	8.0	1.1	20	1.10	3.29	60	-0.7	0.9	0.70	0.50	0.24	100	100
1	0.39	60	6.4	1.4	20	1.11	3.19	60	-0.7	0.9	0.50	0.75	0.99	1000	100
2	0.66	38	5.5	2.1	20	1.24	3.22	60	-1.0	0.9	0.70	0.75	1.58	1000	100
3	0.86	32	6.7	3.4	30	1.25	3.25	100	-0.5	0.9	0.70	0.75	1.75	100	100
4	1.11	21	6.2	5.4	30	1.27	3.38	60	0.0	0.9	0.70	0.75	2.67	500	100
5	1.40	11	7.3	6.0	30	1.16	3.23	60	0.0	0.9	0.70	0.75	2.36	1000	100
6	1.60	9	7.8	5.8	30	1.34	3.31	60	0.0	0.9	0.70	0.75	2.50	500	100
7	1.88	7	15.3	13.3	30	1.23	3.43	60	0.0	0.5	0.70	0.75	2.72	300	100
8	2.12	5	18.8	18.5	25	1.26	3.61	60	0.0	0.5	0.70	0.75	3.03	100	300
9	2.38	5	16.5	18.0	30	1.21	3.52	60	0.0	0.5	0.70	0.75	3.10	550	500
10	2.51	2	22.5	17.6	25	1.47	3.69	100	0.0	0.5	0.70	0.75	3.11	550	750
11	2.85	8	5.2	7.7	10	1.27	5.26	140	-0.5	0.9	0.70	0.50	2.70	550	500
12	3.08	4	5.8	6.9	10	1.31	4.30	80	0.0	0.9	0.30	0.50	2.84	300	750
13	3.31	1	6.5	8.2	10	1.29	3.74	60	0.0	0.9	0.70	0.50	2.83	100	100

Equations

$$m=\frac{-2.5}{\ln 10}\cdot [\mathrm{asinh}(\frac{f/f_0}{2b})+\ln b]$$

$$\left(\frac{{\mathrm{M}}_{BH}}{10^9{\mathrm{M}}_{\odot }}\right)=(0.22){\left(\frac{{\mathrm{M}}_{bulge}}{10^{11}{\mathrm{M}}_{\odot }}\right)}^{1.24}$$

Full text

††pubyear: 2015††pagerange: Probing the rapid formation of black holes and their galaxy hosts in QSOs–Probing the rapid formation of black holes and their galaxy hosts in QSOs

1 Introduction

Since their discovery in the 1960s, theories about quasars have experienced a dramatic paradigm shift. Because they were discovered as strong radio sources, at first it was thought that quasars are examples of a rare and short phase in the evolution of radio galaxies (RGs), with a typical duration time of the order of a few $10^{6}$ years (Sandage, 1965; Gold, 1967). However, based on a photometric survey of blue quasi-stellar objects (QSOs), Sandage (1969) made clear that radio loud quasars only constitute a few percent of a much larger population of QSOs that are intrinsically radio quiet (see Coziol et al., 2017, and references therein), suggesting that this phenomenon must represent a more common phase of galaxy evolution. This last conclusion was readily accepted once the cosmological nature of the redshift was recognized and different surveys revealed QSOs were more common in the past, their density and luminosity significantly increasing at high redshifts (e.g., Croom et al., 2009).

The other breakthrough in QSO studies consisted in identifying the physical source of their activity: a super-massive black hole (SMBH) accreting matter at the centre of their galaxy hosts (Rees, 1978; Soltan, 1982). This implies that QSOs represent a primordial phase in the formation of galaxies during which SMBHs are actively forming at the centre of their nuclei (Cavaliere & Szalay, 1986; Lapi et al., 2006; Letawe, Letawe, & Magain, 2010). Adopting this new paradigm transformed the study of QSOs into an investigation about galaxy formation and the role the formation of SMBHs plays during this process (e.g., Silk & Rees, 1998).

There are three key observations to consider in search for such a connection. The first one is that SMBHs with masses of the order of $10^{9}$ M⊙ are found to exist well above a redshift $z\sim 6$ (Fan et al., 2003; Mortlock et al., 2011; Wu et al., 2015; Bañados et al., 2018; Yang et al., 2020; Wang et al., 2021). This observational fact emphasizes how fast must be the formation process of SMBHs (Woods et al., 2019; Inayoshi, Visbal, & Haiman, 2020; Pacucci & Loeb, 2022). For example, consider J0313$-$1806, which is the farthest QSOs discovered so far at $z=7.642$ (Wang et al., 2021). Based on its emission lines, and adopting a standard $\Lambda$CDM cosmology (with $H_{0}=70$ km s*-1* Mpc*-1*, $\Omega_{DM}=0.30$, and $\Omega_{\Lambda}=0.70$), this QSO, appearing at a relatively young age of the universe, that is, a cosmological time (CT) of only 0.667 Gyrs, would already have reached a mass of $(1.6\ \pm\ 0.4)\times 10^{9}$ M⊙. Now, assuming the formation of the SMBH started just after recombination at $z=1,100$, this implies a formation time $\Delta t\sim 7\times 10^{8}$ years. How can that be possible? For instance, assuming the black hole (BH) mass increases solely by accretion, at a rate $\dot{M}_{BH}=(1-\eta)\dot{m}_{acc}$, adopting a typical radiation efficiency $\eta=0.1$ would imply a constant rate of the order $\dot{m}_{acc}\sim 3$ M${}_{\odot}/{\rm yr}$, which is slightly lower than necessary to produce its observed bolometric luminosity, $L_{bol}=3.6\times 10^{13}$ L⊙, or $\sim 0.7$ $L_{\rm Edd}$. However, assuming a constant accretion rate just after recombination is highly improbable, considering that $\dot{m}$ varies as the mass of the BH grows, and its growth itself is entangled with convoluted initial conditions related to the formation of its host galaxy, which, within the context of hierarchical structures formation (massive galaxies forming by the mergers of massive, gas-rich protogalaxies), is a highly complex, non-linear physical phenomenon; this explains why some researchers have considered super massive BHs (SMBHs) form by accreting matter at super Eddington rates, that is, well above the Eddington limit.

The second clue about the formation of SMBHs is that their masses are correlated with the masses of the bulges of their host galaxies (the so-called BH mass-stellar velocity dispersion relation, or $M_{\rm BH}-\sigma$; Magorrian et al., 1998; Ferrarese & Merritt, 2000; Häring & Rix, 2004; Gültekin et al., 2009; Graham et al., 2011). For J0313$-$1806, assuming the relation deduced based on local AGNs applies at high redshifts (since the formation of SMBH is possibly fixed early on; Fan, 2006; Shen et al., 2019), would yield a velocity dispersion of stars in the bulge equal to $\sigma\sim 323$ km s*-1* (using the relation in Graham et al., 2011), which is typical of elliptical galaxies with minimal mass $\sim 1.5\times 10^{11}$ M⊙ (according to Reines & Volonteri, 2015). For QSOs with SMBHs, this not only suggests their host galaxies formed at the same time as their SMBHs, but also that they form the bulk of their stellar populations extremely rapidly, a formation process that is typical of massive, bulge-dominated, early-type galaxies (Sandage, 1986; Warner, Hamann, & Dietrich, 2003; Bischetti et al., 2021).

A rapid formation of galaxy host (or galaxy bulge) is also implied by the third observational clue, which is that in any QSO at any redshift, the gas in the broad line region surrounding the accretion disk has a solar or higher than solar metallicity (Hamann & Ferland, 1993; Jiang et al., 2007; Juarez et al., 2009; Śniegowska et al., 2021). Since metals are formed by stars, the high metallicity of the gas accreting on a SMBH could only mean higher star formation rates in the past accompanied the rapid formation of the SMBHs (the process possibly favouring massive stars; a starburst with flat initial mass function, IMF). A fast formation process for galaxies might also explain the various mass-metallicity relations observed (Dietrich et al., 2003a; Neri-Larios et al., 2011; Matsuoka et al., 2018).

But probably the most constraining facts about the abundances is that the line ratio Fe ii/Mg ii observed in any QSO at any redshift is relatively high (between 2 and 4; De Rosa et al., 2014; Wang et al., 2022), which, taken as a proxy for the ratio Fe /$\alpha$-elements, suggests a rapid enrichment in Fe  at an epoch much earlier than the CT (Dietrich et al., 2003b; Barth et al., 2003; Kurk et al., 2007; Shin et al., 2019). That characteristic could represent a problem for QSOs at high redshifts, because, according to the standard interpretation, while $\alpha$-elements are produced by SNe II and Ib, and thus their abundances in the interstellar medium are expected to increase rapidly during a massive burst of star formation, Fe  is produced by SNe Ia and thus its abundance in the interstellar medium should be delayed compared to the $\alpha$-elements by a at least 1 Gyr, which is much longer than the typical time-scale for the growth of SMBH by accretion (a few $\times 10^{8}$ years Inayoshi, Visbal, & Haiman, 2020). However, this apparent difficulty might also translates into an intriguing alternative, which is that the protogalaxies that merged to form the bulges of galaxies (Tinsley & Larson, 1979; Silk & Norman, 1981) were populated by very massive and metal-poor stars, consistent with the elusive Pop III (De Rosa et al., 2014): since these stars are metal-poor and massive they would eject in the ISM very few Mg ii compared to Fe  as they evolve, explaining the typical high ratio observed in QSOs (Thomas, 1999; Heger & Woosley, 2002).

The hypothesis of protogalaxies formed by Pop III stars would also readily fit within a model for the formation of SMBH seeds based on the very popular hierarchical galaxy formation model (Devecchi & Volonteri, 2009; Yajima & Khochfar, 2016; Sakurai et al., 2017; Reinoso et al., 2018): 1- a huge quantity of pristine gas fall into dark matter (DM) mini-haloes, with masses of the order $10^{5-6}$ M⊙, to form protogalaxies where the IMF is top-heavy, 2- the most massive stars ($m_{*}\geq 10^{2}$ M⊙) at the centre of these protogalaxies rapidly evolve and merge to form SMBH seeds with masses of the order $M_{BH}\sim 10^{3-4}$  M⊙, 3- the mass of these seeds rise rapidly up to $M_{BH}\sim 10^{5-6}$ M⊙ through haloes/protogalaxies/BHs mergers, forming galaxies with massive bulges where a central SMBH rapidly increases its mass by gas accretion up to $M_{BH}\sim 10^{8-9}$ M⊙.

However, this is the naive scenario, and details how it really happened are still open to discussion (see Inayoshi, Visbal, & Haiman, 2020). Some important questions that need to be answered are the following. Does the formation process of SMBH in high density regions (at high z) differ from what happens in lower density regions (intermediate and low z)? One possibility is that SMBHs observed at high redshifts are rare objects forming only in the most massive DM haloes ($>10^{8-9}$ M⊙), suggesting that the majority of SMBHs have much lower masses and thus their formation process (and that of their host galaxies) might have been different. Some researcher also claimed that the most massive SMBHs at high redshifts reside in slightly less massive galaxies than observed at lower redshifts, which would imply their formation pre-dated the formation of their hosts (Shimasaku & Izumi, 2019; Vayner et al., 2021). However, this would be difficult to explain, considering the spectral similarities and metallicities of QSOs at any redshift. Related topics consist in searching for correlations between the BH accretion rate (BHAR) and star formation rate (SFR) in AGNs at different redshifts (e.g., McDonald et al., 2021) and determining what role (if any; e.g., Peng, 2007; Jahnke & Macciò, 2011) AGN feedback (outflows or winds) play in establishing the $M_{\rm BH}-\sigma$ relation (e.g., Harrison et al., 2018; Torres-Papaqui et al., 2020).

In this article, we present the results of a study that aims to better constrain the formation galaxies hosting SMBHs by tracing the change of SFRs with redshift in a sample of QSOs with different BHAR. Although the UV-Opt part of the spectrum of a QSO is dominated by the AGN continuum and intense broad line emission, information about star formation in their host galaxies still appears in their SEDs through their dominant stellar populations and dust emission in MIR and FIR. In principle, therefore, the host galaxy SFR in an AGN can be estimated by reproducing its SED synthetically using sophisticated tools like X-CIGALE (Boquien et al., 2019; Yang et al., 2022),111https://cigale.lam.fr which allows distinguishing in a SED the stellar population components (young and old) from the SBMH component. The organization of this paper is the following: in Section 2, our sample of QSOs is described, and their characteristics in terms of BH mass and BHAR are discussed in Subsection 2.1. Our method of analysis using X-CIGALE is explained in Section 3 and the results are presented in Section 4, followed by a brief discussion in Section 5. Our conclusions can be found in Section 6.

2 Sample and data

To better constrain the SEDs of QSOs using X-CIGALE, we need the largest sample possible with the largest range in high quality fluxes measured. Starting with the QSOs compilation based on SDSS DR12 built by Pâris et al. (2017) and the AllWISE catalog (Cutri et al., 2021), both available through CDS VizieR,222https://vizier.cds.unistra.fr/viz-bin/VizieR we used the X-Match application tool to retrieve a preliminary list of detected QSOs in the MIR by cross-correlating the position of each entry, obtaining 190,415 candidates (64% of the QSOs in Pâris et al., 2017). However, because the fluxes in WISE do not have all the same quality, we restricted our selection to those QSOs having the highest level: Flag A on all four magnitudes, W1 (3.368 $\mu$m), W2 (4.618 $\mu$m), W3 (12.082 $\mu$m), and W4 (22.194 $\mu$m). This reduced our final list of high-quality WISE QSOs (hereafter HQWISE QSOs) to 1,359 QSOs covering a range in redshift from $0<z<4$. Initially we contemplated adding more fluxes in NIR or even FIR. However, due to the small number of HQWISE QSOs and scarcity of information in other bands, cross-correlation with other data banks (in particular, IRAS for data in the FIR) turned out to be fruitless.

Note that we favor the highest quality fluxes rather than high number of objects in order to fully understand the physical consequences related with the variations of the SED of QSOs at different redshifts. To explain better this last point, we trace in Figure 1 the variations with redshift of the IR colours of the HQWISE QSOs. At any redshift, we observe a relatively high variance, which, considering our strict selection criterion for the fluxes, suggests its physical cause must be intrinsic. For the color W2-W3 (middle panel) an increase in variance is most obvious at low redshift (below $z\sim 1$), where many sources have colours typical of high redshift QSOs, suggesting they might share similar characteristics. One possibility is higher than normal star formation (SF) enshrouded in dust, which could be evidence at low redshifts of ultra luminous IR galaxies (ULIRGs), that is, starburst galaxies where SMBHs hide behind thick veils of dust. In 1988, Sanders et al. (1988) proposed that the ULIRG state is a normal phase in the formation/evolution of QSOs and confirming their presence in our sample would thus be a decisive step. This justifies using high quality IR fluxes, since using lower quality fluxes would increase the variance (creating false results) and impede us to have a clear view about this phenomenon.

One important consequences of the ULIRGs model is the connection implied between Starbursts and AGNs. According to this model, an important part of the dust in QSOs must be heated by star bursts, and as SF in galaxies increases at high redshifts (Madau & Dickinson, 2014), an increase of SFR in QSOs would explain why the W2-W3 get redder. At low redshifts ($z<0.25)$, such a reddening sequence, Liners$\rightarrow$Seyfert 2 (Sy2)$\rightarrow$Star Forming Galaxies (SFGs), was found by our group to be empirically connected with a genuine increase of SF (Torres-Papaqui et al., 2013; Coziol et al., 2014; Coziol, Torres-Papaqui, & Andernach, 2015). What happens in QSOs however is more complicated. Comparing the colors produced by different SED models for QSOs at low redshifts (Figure 27 in Coziol, Torres-Papaqui, & Andernach, 2015) suggested mild star formation rates, between 0.3 M⊙ yr*-1* and 10 M⊙ yr*-1* (consistent with SFRs as measured by Diamond-Stanic & Rieke, 2012; Xu et al., 2015). However, this would imply as W2-W3 gets redder in QSOs at higher redshifts, an increase in SFR by a factor 10 or even 100 (e.g., Drouart et al., 2014; Dong & Wu, 2016). Note that the variation of the W3-W4 colour in the lower panel of Figure 1 seems consistent with such scenario, the colours at low redshifts being comparable to those of Sy2s (local AGNs with high SFRs), becoming bluer at high redshifts; the stellar component of the SED becoming flatter as the SFR increase with the redshift.

However, the behaviour of the W1-W2 colour with the redshift in the upper panel of Figure 1 is more complicated to explain. At first, W1-W2 becomes redder as the redshift increases, but after $z\sim 1.4$ it becomes bluer. Note that to locate this inflection precisely, we need high quality fluxes. Adopting the standard interpretation of WISE colours (Figure 12 in Wright et al., 2010, also based on high quality fluxes), a change in AGN type could be suggested, since Seyfert 1 (Sy1s) are bluer than QSOs in W1-W2. However, the shift in colour happening at high redshift contradicts this explanation, since Sy1s (mostly spiral galaxies) appear only at low redshifts (cf. Torres-Papaqui et al., 2020). Also, the fact that the W2 fluxes show no evidence of a similar shift in W2-W3 at this redshift suggests this particularity of the W1-W2 colour must be related to a change in the SED of the QSOs. But why is that change happening at $z=1.4$ is perplexing. This cannot be due to a spectral line entering or getting out of one of the WISE filter pass bands, since at this redshift the most intense emission line in the red in AGNs, H$\alpha$, is still in the blue, at 1.6 $\mu$m, while PAH emission pass bands (e.g. Cesarsky et al., 1996) fall way farther to the red than W1 and W2; besides, PAH emission is usually weak or absent due to the strong AGN continuum component in QSOs (Draine & Li, 2007).

Interestingly, Richards et al. (2009) observed the same trend for the colours 3.6-4.5 $\mu$m (consistent with W1-W2) for a sample 5,546 quasars observed with Spitzer (IRAC instrument). To explain this feature, the authors alluded to the inflection at 1 $\mu$m, separating the blue and red part of the SED, as was detected previously by Elvis et al. (1994) using a sample of 47 quasars at low redshifts. However, the covering range in redshifts does not match. This is better seen in the study made by Labita, Treves, & Falomo (2008, see their Figure 3): the inflection at low redshift is at a frequency of $10^{14.5}$ Hz, and the blue bump at $10^{15}$ Hz, while at $z=1.4$, the part of the SED affecting W1 is only at $10^{14.3}$ Hz, the 1 $\mu$m only reaching the W1 bandpass at $z=2.5$, far from the inflection of W1-W2 colours. The 1 $\mu$m inflection consequently cannot explain the W1-W2 change in colours towards the blue. What is needed is a component of the SED that increases the flux in W1 relative to W2, making this part of the SED flatter. Since the UV bump is still far to the blue at $z=1.4$, this might suggest an extra component, like a star burst in the circumnuclear region (e.g., Xie et al., 2021). This once again legitimizes using high quality fluxes to distinguish what component of the SED can have this property.

Researching ADS for the literature and using NED and CDS VizieR tool to compare our data with available catalogues, we found no evidence of detection in radio or X-ray for any of the HQWISE QSO, which makes this sample a purely UV-Opt/MIR selected sample. For the present study, our data are thus restricted to the five SDSS fluxes in u, g, r, i and z filters and the four high quality fluxes in WISE. In SDSS the magnitudes are expressed as inverse hyperbolic sine function (or asinh), as described by (Lupton, Gunn, & Szalay, 1999).333https://www.sdss.org/dr12/algorithms/magnitudes/ The transformation from linear flux measurements to asinh magnitudes was designed to be identical to the standard astronomical magnitude at high signal-to-noise ratio (S/N) and to behave reasonably at low S/N. The magnitudes are characterized by a softening parameter $b$, which is the typical 1-sigma noise of the sky in a PSF aperture for a 1 arcsecond seeing. The detected flux, $f$, is then obtained by inverting the following relation for the asinh magnitude:

[TABLE]

Where $f_{0}$ is the flux of an object with conventional magnitude of zero. The quantities $b$ as measured relative to $f_{0}$ are given for each filter in Table 21 of Lupton, Gunn, & Szalay (1999), along with the asinh magnitude associated with a zero flux calibrator. The table also lists the flux corresponding to $10\times f_{0}$, above which the asinh magnitude and the traditional logarithmic magnitude differ by less than 1% in flux.

The instruction how to transform the WISE magnitudes in fluxes can be found in Wright et al. (2010).444https://wise2.ipac.caltech.edu/docs/release/allsky/expsup/sec4_4h.html According to the colour distributions shown in Figure 1, we opted to use the standard calibration for the power law $F_{\nu}\propto\nu^{-2}$.

2.1 Properties of HQWISE QSOs.

Information gleaned from two spectral analysis of SDSS spectra allows us to characterize these QSOs in terms of AGN luminosity, BH mass and Eddington ratio. The first analysis is from Kozłowski (2017). Using QSOs from the Pâris et al. (2017)’s sample, within redshift range $0.1\leq z\leq 5.5$, Kozłowski measured in the extinction-corrected SDSS spectra the fluxes and FWHM of the broad Mg ii and C iv emission lines, and their adjacent fluxes in the continuum, respectively at $\lambda=3000$ and 1350 Å, using these measurements to estimate the bolometric luminosities, Lbol, AGN monochromatic luminosities, L${}_{AGN}=\lambda{\rm L}_{\lambda}$, virial masses, MBH and Eddington ratios, $nEdd={\rm Log}({\rm L}_{bol}/{\rm L}_{Edd})$ (see details about measurements and calculations in Kozłowski, 2017). From his compilation, we retrieved the characteristics of the BHs for 1,098 of the HQWISE QSOs (only those with Mg ii and C iv emission lines; Kozłowski only used these two lines). For our analysis, we separated the HQWISE QSOs in 14 redshift bins with width 0.25, except for the last bin which contains only 12 QSOs within the range $3.25<z\leq 4$. We report the medians of the parameters characterizing the QSOs in each bin in Table 1. Note that there is no data for the QSOs in bin #0 because at this redshift the Mg ii line does not appear in the SDSS spectra. The number of QSOs with measurements in each bin is indicated in col. 3, compared to the number of HQWISE QSOs in the bin considered. For the BH mass and luminosity we utilized the values based on the Mg ii line, only using the C iv line at high redshifts where the Mg ii line is not visible anymore.

In figure 2, we present the box-whisker plots for LAGN, MBH and $nEdd$ for our sample of HQWISE QSOs. As expected, LAGN (upper panel) continuously rises with the redshift, by almost three orders, increasing rapidly up to $z\sim 1.62$ (bin #6) and more slowly after that, until $z\sim 2.38$ (bin #9), where the statistics becomes less secured as the number of QSOs in the bins decrease. Because the BH mass estimation depends on both the luminosity and FWHM (cf. equation 1 in Kozłowski, 2017), its uptrend with the redshift (middle panel) naturally follows the increase in luminosity. However, this is only until bin #6 where MBH is almost constant. Due to the lower number of QSOs above $z=2.38$ in our sample, we cannot be sure the trend towards higher BH masses is physical or due to a bias in luminosity. Note that the redshift range $1.63\leq z\leq 2.12$ (bin #6 to bin #8), which is well covered in our sample, includes the epoch in the evolution of the universe when the SFR in galaxies peaks and start to decrease at low redshifts (cf. Madau & Dickinson, 2014). The trend for the BH mass in QSOs to peak at almost the same epoch is usually interpreted as evidence for a connection between the formation of SMBH and their host galaxies (Boyle & Terlevich, 1998; Madau & Dickinson, 2014).

Because the Eddington ratio (lower panel) corresponds to the ratio of accretion rate (BHAR) to BH mass ($nEdd\propto{\rm BHAR}/{\rm M}_{BH}$) and the BH mass grows with the luminosity, the fact that we observe a lower amplitude variation of this parameter with the redshift is not unexpected. However, there seems to be a weak trend for $nEdd$ to increase above $z=1.37$ (bin #5), suggesting higher BHAR at high redshifts than at low redshifts in more massive BH and luminous QSOs.

The second spectral analysis useful for our study was done by Rakshit, Stalin, & Kotilainen (2020), who contrary to Kozłowski (2017) used spectra in SDSS DR14 instead of DR12. Comparing MBH as calculated in these two different sources in Figure 3, we found the results to be almost the same. The linear fit shows that in general the BH masses in Rakshit, Stalin, & Kotilainen (2020) tend to be slightly lower than in Kozłowski (2017), the difference most probably being due to the subtraction by the latter of Fe templates. Indeed, to determine the contribution of Fe in each spectrum, Rakshit, Stalin, & Kotilainen (2020) did a robust fitting analysis of different templates, which allowed them not only to correct the spectra but also estimate the Fe  and Mg ii abundances. In Figure 4, we show the ratios Fe /Mg ii for these two elements in our HQWISE QSO sample based on the abundances calculated by Rakshit, Stalin, & Kotilainen (2020). This result shows the typical high ratios found in QSOs independently of the redshift (compare with Figure 11 in De Rosa et al., 2014). Note that to preserve the homogeneity of our analysis, we will not use the results of Rakshit, Stalin, & Kotilainen (2020) for MBH based on measurements of the broad component of H$\beta$ at low redshifts.

3 Method of analysis

Our analysis consists in reproducing the SEDs of the 1,359 HQWISE QSOs using X-CIGALE (Boquien et al., 2019; Yang et al., 2022).555X-CIGALE manual can be found here: https://cigale.lam.fr/documentation/ The key parameter to estimate how fast the galaxy hosts of quasars formed the bulk of their stellar populations is the star formation history (SFH) function of the host galaxies. In our universe, galaxies form their stars following only two different paths, which in great part also explains their morphology (Sandage, 1986; Charlot & Bruzual A, 1991): they either form most of their stars early on, like elliptical galaxies, or form their stars more gradually like in spiral galaxies, the time-scale difference being due to the high angular momentum of gas collapsing in spirals to form large-scale disks. In X-CIGALE the SFH functions that reproduce these two patterns are called sfhdelayed, for Elliptical-like or early-type, massive bugle galaxies, and sfh2exp, an exponential SFH typical of smaller bulge, spiral galaxies (see Fig. 1 in Boquien et al., 2019). Because we cannot assume a priori the morphology of the QSO hosts in our sample both SFH functions were tried using X-CIGALE.

How fast galaxies with different SFH function form their stars is determined by the e-folding time, $\tau_{\rm main}$, which corresponds to the time it takes a galaxy to form 69% of its stars. The smaller $\tau_{\rm main}$ and the faster the formation of the galaxy. Theoretically (Sandage, 1986) elliptical galaxies have smaller e-folding time than spiral galaxies and very large e-folding times (a few Gyrs) are rare, since they imply SF in these galaxies have been delayed over a long period of time (the cases of small-mass irregular galaxies). However, assuming the host galaxies of QSOs are massive and the $M_{\rm BH}-\sigma$ relation is universal, star formation cannot be delayed over arbitrarily long periods of time compared to the formation of their SMBHs, and we would thus expect the $\tau_{\rm main}$ produced by X-CIGALE to be relatively small. How small can it be, this is the condition that was tested (as shown in Table 2) by trying different values for $\tau_{\rm main}$, varying between 50 Myrs to 7000 Myrs, that is, all the values allowed in X-CIGALE (see Table 2).

One important module in X-CIGALE tackles the AGN contribution. Because QSOs are Type-1 AGNs, we expect the AGN continuum to dominate the SED in the UV-Opt and NIR (most specifically the W1 and W2 passbands). Having experimented with the two main options, skirtor2016 (Stalevski et al., 2016) and fritz2006 (Fritz, Franceschini, & Hatziminaoglou, 2006), and having obtained slightly better fits (lower $\chi^{2}$) with the latter, we decided to use only this module for our final testing. The main difference between the two models lies in the isotropy of re-emission by dust, the torus in fritz2006 model being less clumpy than in the fritz2006 model, producing more isotropic re-emission. Consistent with Type-1 AGN in the Fritz module, we fixed the line of sight to $89.99^{\circ}$ (implying we see the torus of dust face on) but left the program free to try different inputs: most specifically, varying the fraction of AGN, fAGN, and opening angle, OA, of the dust torus (see manual of X-CIGALE and original article by Fritz, Franceschini, & Hatziminaoglou, 2006, for detailed descriptions).

One crucial improvement to X-CIGALE made in 2022 by Yang et al. (2022) is the addition of a $\delta$ parameter, which fixes the slope of the AGN continuum between $0.125\ \micron<\lambda\leq 10\ \micron$ ($\lambda{\rm L}_{\lambda}=\lambda^{-0.5+\delta}$). Without this parameter, the contribution of the AGN to the SED would be underestimated, increasing artificially the stellar contribution in the UV-Opt, which then translates into higher SFRs. Having experimented with the previous model of X-CIGALE without the slope correction, we did observed such SFR overestimation in the host galaxies of our QSOs, although not as high as we could have expected. In our modelling, different values of $\delta$ were tried, varying between -1 and 0 (see Table 2).

Another module that is important for the MIR range (W3 and W4 in WISE) is the module related to re-emission in IR by dust of UV light due to SF. Experimenting with the different options offered by the program (see Figure 5 in Boquien et al., 2019), we have adopted the dl2014 module (Draine et al., 2014). After experimenting with different values for the PAH mass fraction, we used in our final modelling the minimal value $q_{\rm PAH}\sim 0.45$% (Draine & Li, 2007). Indeed, even with high SFR in the QSO hosts, the PAH features should not play a great role due to the dominant AGN component (e.g., Marshall et al., 2018); the fine dust producing PAH emission can also be destroyed by the intense AGN luminosity (Ramos Almeida et al., 2023).

For the single stellar population (SSP), we used the module bc03 based on Bruzual & Charlot (2003), adopting a Salpeter’s IMF. After experimenting with different values for the metallicity, from $Z=0.004$ to $Z=0.05$, we determined that our best fits were those where the metallicity was solar or higher at any redshift (consistent with observations; cf. Figure 4). For the attenuation of UV-Opt light due to dust extinction, we used the module dustatt_modified_starburst (based on Calzetti et al., 2000; Leitherer, Calzetti, & Martins, 2002), which includes a bump at 2200 Å with varying amplitudes.

Finally, after experimenting with different SED models without polar emission, we realized that we can greatly improve the fits by including such component in all the HQWISE QSOs, at all redshift. This was done using the special module introduced in X-CIGALE by Yang et al. (2020). In the literature, the existence of polar dust components (excess of dust emission emitted in a direction perpendicular to the obscuring torus; Tristram et al., 2014) was already noted and thoroughly studied by various authors in different AGNs (e.g., Toba et al., 2021; Isbell et al., 2022; Lyu & Rieke, 2022), suggesting this characteristic is more common than previously thought. The inclusion of polar dust in our models systematically produced best fits with $\chi^{2}<3$ for most of the QSOs in our sample (87% using sfhdelayed and 83% using sfh2exp).

In Table 2, we summarize the inputs used in X-CIGALE for our final solutions, indicating the range of values for the parameters that we allowed to varied to get the best SED models (i.e., the models that minimize the $\chi^{2}$). Note that for some parameters, we tried more values (indicated in parenthesis) in separated runs, to compare with our best models. The mock option was also used to verify the robustness of our final solutions (see full decription in Mountrichas et al., 2021).

4 Results

In this section, we compare in parallel in Figure 5 the best solutions in each bin, as obtained using the two SFH functions sfh2exp (left panels) and sfhdelayed (right panels). In Tables 3 and  4, we summarize the results for both functions by giving the median values for important parameters in each bin, calculated considering only those models with $\chi^{2}<1$ (61% of our sample for both solutions). Comparing col. 3 in the two tables, we conclude that both SFH functions yield equally reasonable solutions. In general, the AGN component dominates the SEDs (col. 10) and the SFR is high, increasing at high redshift (col. 13). The opening angle is huge (col. 8) and extinction relatively large (col. 11). These SEDs are consistent with obscured starbursts with dominant AGNs (Marshall et al., 2018; Xie et al., 2021).

Examining in detail col. 3 in Table 3 for the sfhdelayed function, the first three bins at low redshifts (#0, #1 and #2) have slightly higher $\chi^{2}$, the fits getting better after that, up to bin #6 ($z=1.62$), the $\chi^{2}$ then fluctuating but still staying relatively low until the last three bins (#11, #12, #13) where the numbers of QSOs is falling drastically. Disregarding the last three bins, we assert that the sfhdelayed function yields excellent solutions over a range of redshift $0.8<z\leq 2.8$. This conclusion is also supported by the fraction of models with $\chi^{2}>3$ as indicated in col. 4, which does not pass 10% within this redshift range.

Doing the same analysis for the sfh2exp function in Table 4, we also see higher $\chi^{2}$ values and higher fraction of models with $\chi^{2}>3$ in the first and last three bins, but comparably lower values, as observed for the sfhdelayed function, within the range $0.8<z\leq 2.8$. This suggests that our solutions are degenerate, since we cannot differentiate between the two SFH functions. The reason for this degeneracy is easy to understand: the common traits of the SEDs at any redshift are the dominance of the AGN component (with some exceptions at high redshifts, for examples in bins #9 to #12) and the high SFR (col. 13). The high SFR is consistent in both models with relatively small e-folding (col. 14): not passing 1,000 Myrs in the sfhdelayed models and 750 Myr in the sfh2exp, suggesting that, independently of their SFH function (morphology), the host galaxies of HQWISE QSOs generally built the bulk of their stellar populations rapidly.

In our sample the dominance of the AGN and the high SFRs explain the degeneracy. However, examining in detail the best fits in Figure 5 some important nuances must be added. More specifically, we observe significant differences in the first five bins at low redshifts in terms of the stellar populations (unattenuated curves in blue in the figures), suggesting they are systematically older using the sfhdelayed function than the sfh2exp function. This characteristic is consistent with slightly larger e-folding and much lower SFRs in the sfhdelayed SEDs than in the sfh2exp SEDs. On the other hand, starting at bin #6 ($z\sim 1.62$) the stellar components in both models become similar. This happens as the cosmic time (CT in col. 3) becomes comparable to the e-folding time (CT falling from 3.9 to 1.7 Gyrs).

Because the SEDs of the stellar populations in the sfhdelayed models at low redshifts are typical of early-type galaxies, the above difference could be interpreted as two solutions for two different morphological types for the host galaxies, Elliptical vs. Spiral. Considering the evolution in CT and assuming a standard scenario for the formation of galaxies (Madau & Dickinson, 2014), finding high SFRs in galaxies at low redshifts seem more natural in spiral than elliptical galaxies (e.g., Bellstedt et al., 2020). This means that the consequences of the two solutions are significantly different: in the sfh2exp models, the SFRs would be characterized as normal or mildly starburst-like for late type spiral galaxies, while in the sfhdelayed models the high SFRs would be definitely starburst-like for the morphology (consistent with blue Elliptical galaxies or merging galaxies).

Continuing our analysis of the results, another important feature that is remarkable in Figure 5 is the similarity of the re-emission dust components in both models. This suggests that the WISE colours produced by the two SFH functions must also be similar. To check how well X-CIGALE reproduces the variations of these colours with the redshift, we trace in Figure 8 the medians of the models with $\chi^{2}<1$ in each bin, comparing them with the medians measured by WISE. For both models, the fits are relatively good, except maybe in the last two bins where the number of QSOs are very low (only 13 and 12 QSOs respectively). In particular, both models reproduce relatively well the inflection to the blue of the W1-W2 colour at $z\sim 1.4$. Coincidentally, this happens in Figure 5 passing from bin #5 to bin #6 (medians $z=1.38$ and $z=1.62$, respectively), where the stellar contributions in the two models become comparable: the SED becoming gradually bluer (flatter). Consequently, the inflection of W1-W2 colours in both models could be interpreted as the signature of an increase of SFR in the host of QSOs at high redshifts. Based on the literature (Elvis et al., 1994; Richards et al., 2009; Labita, Treves, & Falomo, 2008), this might be recognized as a common trait of quasar evolution.

Although both SFH functions predict high SFRs increasing with the redshifts, considering the possible difference in morphology, determining which SFR evolution is more realistic is not obvious. Comparing the two SFR solutions in Figure 7, the sfh2exp SFRs are generally higher than the sfhdelayed SFRs. The differences being quite significant at low redshifts, varying from 100 to a few 10s below and up to $z=1$ (first four bins), becoming slightly lower by only a few after bin #4 (median $z\sim 1.11$), the differences decreasing in the ultimate 5 bins.

To establish the reliability of each solution we used the option in X-CIGALE to generate two mock catalogues based on the best fit models for the two SFH functions (see examples in Mountrichas et al., 2021). In Figure 8, we compare the mock SFRs with the fitted SFRs for the two SFH functions. The principle of this test is simple: using the fitted values as guesses, running X-CIGALE should converge to these values rapidly, showing that the solutions are robust. This is easily confirmed for the sfh2exp, the linear regression having a coefficient of correlation of R $=0.89$. Although the sfhdelayed also yields relatively good fits, with coefficient of correlation R $=0.88$, the residuals suggest the solutions at low redshifts where Log(SFR) $<0$, are somewhat problematic, since the fitted SFRs are systematically below the mock SFRs and many QSOs show very large discrepant values at any redshift. This behavior is easy to explain: as the CT increases compared to the e-folding of the sfhdelayed function, the SFR cannot be kept as high as what the sfh2exp and mock values suggest. This also explains why the e-folding of the sfhdelayed function tend to be systematically higher than the e-folding of the sfh2exp. The difference between the two SFH functions has an important consequence on our interpretation. Although at high redshifts the results are degenerate, both requiring small e-folding and high SFR, the starburst like nature of QSOs hosts at low redshift (compare with Figure 9 in Kennicutt & Evans, 2012), implies that, independent of the morphology of the galaxy host, their SEDs are best reproduced using an sfh2exp function (e.g., Coziol, Doyon, & Demers, 2001). Consequently, one would expect many of the nearby HQWISE QSO hosts to show a perturbed morphology, consistent with interactions or mergers (e.g., Bahcall et al., 1997).

Based on the fact that the QSO hosts are starburst-like, and considering the solutions are degenerate at high redshifts, we decided to adopt the sfh2exp SFR in order to examine what could have been their evolution at different redshifts.

4.1 Possible presence of ULIRGs in QSOs at low redshifts

Adopting the sfh2exp SFRs, there is still one difficulty remaining, which is understanding why in Table 4 so many models (17%) yield $\chi^{2}>3$. Since this happens more frequently at low redshifts, one possible explanation, consistent with spiral host galaxies, is that we see the QSOs edge-on, relative to the torus of dust (typical of Type 2 AGN). Alternatively, an early phase of evolution during which the BH is enshrouded by a veil of dust, as expected for ULIRGs, might also explain the observations. The fact that we find many HQWISE QSOs in Figure 1 with W2-W3 colours comparable to QSOs at high redshifts, where the SFRs were confirmed by our models to be high, would support such interpretation. According to the ULIRGs hypothesis, however, one would also expect these QSOs to have both higher than normal IR fluxes and SFRs (e.g., Veilleux et al., 2009).

To test the above two possibilities, we have isolated all the QSOs with $\chi^{2}>3$ in Table 4 and run X-CIGALE, changing the $\psi$ angle using the smallest values possible. The medians in each redshift bin for these new models are reported in Table 5. Comparing col. 5 with col. 4, we do see a significant improvement in the fits but only until bin #6, after which the improvement becomes marginal, the new $\chi^{2}$, in fact, being worst in the last three bins.

In general, lowering $\psi$ produces redder WISE colours and slightly higher SFRs but no systematic changes of the other parameters (most specifically, those related to extinction). Although these changes could be seen as favoring the ULIRG hypothesis, the fact that we see a significant improvement of the fits only at low redshifts could rather be interpreted as consistent with the geometry hypothesis, since AGN in spiral host galaxies are expected to be more numerous at low than high redshifts. In quantitative terms, the fraction of spiral-like hosts or edge-on (EO) QSOs up to $z\sim 1.6$ would amount to $\sim 19$%, which is roughly consistent with what is usually observed in morphological studies of QSO hosts in nearby samples (e.g., Letawe, Letawe, & Magain, 2010).

Adopting the sfh2exp function and correcting for EO cases we trace the variation of SFRs as a function of the redshift in Figure 9. The models predict a significant increase of SFR by a factor 100 from $z=0.5$ to $z\sim 2.5$, with no convincing evidence for a plateau or turn-over.

5 Discussion: connecting the formation of the galaxy host to the formation of its SMBH

The two most important results using X-CIGALE for our sample of QSOs are: 1- independently of the SFH function (morphology), the e-folding is small, of the order 1 Gyrs or less, and 2- the SFRs are relatively high and increasing with the redshift. These results imply that the galaxy hosts of HQWISE QSOs formed their star rapidly in a phase similar to starburst. What remains to be determined, however, is how the formation of the stellar populations are related to the formation of the SMBHs.

Considering the rapid formation of both the SMBH and galaxy host stellar populations, and assuming the M${}_{BH}-\sigma$ relation is established early during the formation of the galaxy, it seems natural to assume the relation established for AGN at low redshift extend at high redshifts. What should we expect then for the mass of the galaxy host, or most specifically the mass of its bulge, M${}_{BH}-$Mbulge? According to Table 1 the SMBH in our sample of QSOs are quite massive, and consequently their galaxy host should also be massive. As a first approximation, we used the relation MBH-Mbulge recently determined by Ramsden et al. (2022) for a sample of SMBH with masses ranging from $10^{6}$ to $10^{10}$ M⊙:

[TABLE]

This yields masses ranging between $10^{10}$ and $\sim 10^{12.5}$ M⊙ (very few above this value). Using the masses of stars (young + old populations) produced by X-CIGALE we can compare them with the masses obtained using Equation 2. This is done in Figure 10. Based on the one-to-one relation, we conclude that X-CIGALE masses are comparable to those produced by the relation of Ramsden et al. (2022), which suggests that our SED solutions using X-CIGALE are consistent with the idea of a universal M${}_{BH}-\sigma$ (at least up to $z\sim 4$). In general, the X-CIGALE masses seem slightly over-estimated above 11.5 dex. For this reason, we choose to use the masses determined by Equation 2 for the rest of our discussion.

Tracing in Figure 11 the variation of SFR with the mass of the bulge, there is a clear trend for the SFR to increase with the mass, a factor 100 increase in mass corresponding to an increase by 1,000 in SFR. Because the BH mass in our sample of QSOs increases with the redshift (Figure 2), this suggests that more massive BHs form in more massive galaxies at high redshifts, a phenomenon known as downsizing in the literature and which is accepted as a trait of hierarchical structure formation (cf. Neistein, van den Bosch, & Dekel, 2006).

However, the dispersion in SFR is large, especially after $\sim 11.7$ dex in bulge mass. Note also that this dispersion does not depend on the emission line used, since we already observe it for Mg ii. What explains this dispersion is not clear. This could either imply a variation in SFR related to different phases of evolution, or different SFHs due to different host morphologies (due to differences in environment densities). At 11.7 dex in bulge mass the BH mass is $10^{9.25}\ {\rm M}_{\odot}$, which in Figure 2 is a typical value in bin #4 ($z\sim 1.1$), that is, close to where we observe the W1-W2 inflection and a possible change in SFH at high redshift, from late to early; due to the small e-folding and high SFR, the two solutions become indistinguishable at high redshifts.

In Figure 12 we examine how the SFR correlates with two parameters related to the AGN activity: the luminosity (upper panel) and Eddington ratio (lower panel). In general, and clearer than for the bulge mass, the SFR increases with the luminosity. There also seems to be a slight inflection in SFR above L${}_{AGN}=10^{46}$ erg s*-1*, the SFR increasing more rapidly at high luminosity. In Table 2, L${}_{AGN}=10^{46}$ erg s*-1* is once again typical in bin #4. At high redshifts, consequently, the higher the AGN activity the higher the stellar formation activity. Since higher AGN luminosity implies higher radiation pressure, this result is a strong argument against the AGN quenching hypothesis: contrary to expectation, AGN feedback related to radiation pressure does not quench SF in the AGN host galaxies.

Examining the lower panel, there is also a weak trend for the SFR to increase with the Eddington ratio. This is consistent with the relation between SFR and LAGN, since higher accretion rates produce higher AGN luminosities. The reason why the relation is weak is also easy to comprehend, since, according to the definition of ${nEdd}$, any increase in BHAR is balanced by an increase in BH mass. Considering this factor, we can conclude that more massive BH form in more massive galaxies and that this implies higher BHAR and SFR. This is another argument against AGN quenching SF in their hosts. Another important trend we can deduced from Figure 12 is that QSOs at high redshifts, those in yellow (identified as C iv), have clearly higher ${nEdd}$, implying that BHAR increases with the redshift.

From all the above relations, we conclude that the more massive the galaxies the more massive its BH, which translate into higher BHAR, and thus luminosity, and also higher SFR. This suggests that the mass of the galaxy host in QSOs is the primary parameter determining the AGN characteristics.

Following this line of thought, we trace the Eddington ratio as a function of the mass of the bulges in Figure 13. Only a few QSOs have super-Eddington ratios ($nEdd\geq 0$ dex). Those are Narrow line QSOs candidates candidates (NLQSOs; Krongold, Dultzin-Hacyan, & Marziani, 2001; Rakshit et al., 2021), where lower than normal BH masses (not necessarily supper-Eddington accretion rates) explain their high $nEdd$. In general, therefore, the majority of the HQWISE QSOs have $nEdd$ below the Eddington limit, suggesting that the BH growth in these QSOs is strictly limited by radiation pressure. How do we reconcile this result with the masses of the galaxy hosts? Obviously, there is no associated limit in mass: the upper envelop of the $nEdd$-Mbulge distribution peaks at a bulge mass of $10^{12}$ M⊙, then decreases afterward as the mass continues to rise. Coincidentally, in galaxy formation simulations, $10^{12}$ M⊙ is recognized as a limiting mass above which gas falling into an halo of dark matter would not have enough time to cool down and form stars (Wechsler & Tinker, 2018; Das, Pandey, & Sarkar, 2021).

Although in Figure 11 we do not see a clear cut in SFR above a mass of $10^{12}$ M⊙, except for a few NLQSOs candidates, the SFR does not seem to grow much higher either. Assuming the galaxy host of QSOs form their stellar populations rapidly, over 1 Gyrs or less according to their e-folding, and considering the increase of SFR with the mass of the galaxies, then the higher the galaxy mass the faster must have been its formation (Zhang, Yang, & Guo, 2021; Varma et al., 2022). Such rapid formation would deplete the reservoir of gas rapidly and this would limit the mass of the bulge. A depletion rate of gas growing with the mass might also explain the decrease of BHAR above $10^{12}$ M⊙ seen in Figure 11. Adopting Ramsden et al. (2022) relation, for example, the fraction of the mass of the BH to the mass of the bulge increases with the mass from 0.19% to 0.5%, suggesting faster BH growth, and thus more rapid depletion of gas to accrete as the mass of the bulge increases. This suggests that physical constraints on the formation process of galaxies could have limited the growth of the bulges and possibly also the growth in mass of their SMBHs.

According to the AGN quenching hypothesis, the regulation of BH mass by AGN winds and outflows driven by radiation is assumed to have limited the mass of their host galaxies by quenching their SF. This process was proposed as the main explanation for the MBH-$\sigma$ relation. However, considering the rapid formation of the host galaxies the relation could be inverted: physical constraints on the formation process of galaxies limit the growth of the bulges and of their SMBHs. One way to distinguish which path the QSO hosts in our sample followed consists in comparing the variation of specific BH accretion rate (sBHAR) with specific star formation rate (sSFR) in different redshift bins. This is done in Figure 14, starting with the lower bins (upper left panel) and ending with those in the highest bins (lower right panel). In this figure, the diagonal line represents the one-to-one relation, implying equal growth rates for the BH and its galaxy host. In our sample, very few QSOs are found on the diagonal, the majority being on the left side, suggesting BHs in QSOs grow more rapidly in mass than their host galaxies. Also indicated in Figure 14 are two different limits in sSFR identifying the Green valley and quenched region according to Bait, Barway, & Wadadekar (2017). Only a few QSOs in our sample cross the Green valley limit (only two cross the quenching limit), and they are all at low redshifts (bin #6 or lower). The general trends as the redshift increases is for the sSFR to get farther away from the quenched region and for sBHAR to increase. There is consequently no evidence of AGN quenching of SF in the HQWISE QSOs. Comparing Figure 14 with Figure 11 in McDonald et al. (2021), although the sBHAR of the HQWISE QSOs agree with those estimated by Dong & Wu (2016), 207 quasars detected in all three band of Herschel-SPIRE, our sSFR are much lower. Since the SFRs estimated by X-CIGALE are comparable to those estimated by these authors (see their Figure 6), the difference could only be in the masses of the bulges. This is confirmed in Figure 10, where the masses of our host galaxies are higher by a factor ten (compared with Figure 8 in Dong & Wu, 2016). This is due to a difference of sample, not of model. Comparing now with the RGs (radio loud AGNs) observed by Drouart et al. (2014), 70 galaxies at redshifts $1<z<5.2$, detected using the PACS and SPIRE instruments of Herschel Space Observatory, the main difference is in their sBHAR consistent with Super-Eddington QSOs, which are rare in the HQWISE QSOs. Since their BH masses are only slightly smaller than what we used, it is their higher BHARs that explain the difference. In a way, our solution makes more sense for radio quiet AGNs, since how a SMBH can accrete at higher rates than the Eddington limit is not physically explained. However, conditions could be more extreme in radio loud AGNs, which is what Drouart et al. (2014) considered.

Also comparing with Figure 11 of McDonald et al. (2021), the HQWISE QSOs at high redshifts have comparable sSFR as the Sy2s studied by Xu et al. (2015) at low redshift. However, the HQWISE QSOs have also higher sSFR than the mixed sample of local Seyfert studied by Diamond-Stanic & Rieke (2012). Once again, the main difference with these two samples seem to be in the masses of the galaxy hosts. Despite these differences, however, what is generally observed for AGNs is a sequence in sBHAR-sSFR parallel and to the left of the one-to-one relation suggesting that the AGN phenomenon is characteristics of a phase in galaxy evolution where the BHs grows in mass more rapidly than their galaxy hosts.

Considering the whole family of AGNs, we believe that the sequences in evolution left to the diagonal is physically significant. Depending how we estimate the stellar mass, the difference with a one-to-one relation growth could be smaller or larger, however, not enough to justify the AGN quenching hypothesis. Moreover, in the HQWISE QSOs the evolution of the sBHAR-sSFR sequence with redshift in Figure 14 is clear: after reaching the Eddington limit in bin #6 the trend is for QSOs to move to the left (lower sSFR) and down (lower sBHAR), which is consistent with a decrease of sSFRs and sBHAR as the galaxy and BH mass decrease and CT increases, which is consistent with downsizing. Considering the fast formation of the bulge, the present state of QSOs to the left of the one-to-one relation could imply that the sSFR must have decreased more rapidly than the sBHAR. Since we see no evidence of AGN quenching in our sample, a rapid formation of the bulge, implying higher SFRs a few Gyrs before we observe the QSOs (as indicated by the e-folding), would decrease the reservoir of gas rapidly, reducing the sSFR below the one-to-one relation with sBHAR. A difference in evolution time scales might also be easy to justify physically, since the BH growth in a galaxy not only depends on the reservoir of gas a galaxy attracts through its mass but also on funneling this gas towards the center of the galaxies and how fast this gas is accreted by the BH, two processes which could easily be limited/regulated by radiation pressure through AGN feedback.

In general, therefore, our results suggest that quenching of star formation through high astration rates (high star formation efficiencies) is the main process limiting the mass of galaxies during their formation. Consistent with downsizing, massive galaxies forming before less massive ones, we would thus naturally expect to see the efficiency in star formation to increase at high redshifts (Inayoshi et al., 2022; Yang et al., 2022), also possibly explaining massive galaxies forming at earlier epoch than previously expected according to the present cosmological paradigm (Robertson et al., 2022; Curtis-Lake et al., 2022; Finkelstein et al., 2022).

6 Conclusions

The modelling program X-CIGALE (Boquien et al., 2019; Yang et al., 2022) integrates different physical modules, representing the quintessence of our understanding about galaxy formation and evolution processes, to reconstruct from various data available in different passbands, from the FUV to the FIR, the SED of galaxies at different redshifts, allowing us to better understand their specific assembly histories. By applying X-CIGALE on a carefully selected SDSS-WISE sample of 1,359 QSOs, within a range $0<z\leq 4$, we were able to show that their host galaxies, 1- have high SFRs, consistent with starburst galaxies, that rise with the redshift and the mass of the galaxy hosts, and 2- based on their e-folding these hosts must have formed the bulk (69%) of their stars rapidly, within 50 to 1000 Myrs of their formation.

Some particularities of the SEDs also worth mentioning: 1- polar dust emission is present in all the QSOs at any redshift, 2- in 19% of the QSOs at $z<1.6$, small angle of view (20-30 degrees) are favored, consistent with a dust torus seen edge on, suggesting that 19% of the QSO galaxy hosts at low redshifts have a spiral morphology, or are consistent with hidden AGNs in low redshift ULIRGs.

Comparing the SFRs in the QSO hosts with the characteristics of their SMBHs, we were able to deduce that, 1- the SFR increases with the mass of the bulge and redshift, in a way that is consistent with downsizing, more massive galaxies, hosting more massive BHs, forming at higher redshifts than less massive galaxies, hosting less massive BHs, and 2- the SFR increases with the BHAR and AGN luminosity, suggesting there is no evidence of AGNs quenching SF in their hosts.

Finally, comparing the sSFRs with sBHARs at different reshifts, we found that, 1- the sSFR is lower than the sBHAR, which implies that QSOs are a special phase in the evolution of galaxies during which the BH grow more rapidly in mass than its host galaxy, and 2- the sSFR and sBHAR increase in the past, the latter being limited by luminosity pressure at the Eddington limit.

A simple interpretation in terms of variation of SFR consistent with the hierarchical model of galaxy formation and downsizing suggests that galaxy hosts of QSOs form their stars rapidly, exhausting their reservoir of gas over shorter time scales than the growth in mass of their SMBH, explaining why their sSFR are observed to be lower than their sBHAR. A higher efficiency of star formation in the past, rising with the mass of the galaxies, makes quenching by astration the dominant process by which the mass of galaxies is limited.

Acknowledgements

K. A. Cutiva-Alvarez acknowledges CONACyT for its support through grant CVU 940597. She would also like to thank Dr. Karla A. Alamo-Martínez for helping her with the installation of X-CIGALE, her friends and colleagues, at Guanajuato, Dr. Aitor C. Robleto-Orús for the open discussions about her work and helps in learning Pyhton, and at Guadalajara, Juan Pablo Gutierrez for his support and guidance in programming with python. Finally, she would like to acknowledge her compatriot and friend Dr. Andres Felipe Ramos, for clarifying some of her doubts in using X-CIGALE, as well as Dr. Guang Yang and Dr. Jianwei Lyu for their openness in discussing some aspects of her study with them, most specifically, emphasizing the importance of polar dust components in AGNs. For his part, J. P. Torres-Papaqui acknowledges DAIP-UGTO (Mexico) for grant support 0077/2021.

This research has made use of the VizieR catalogue access tool, CDS, Strasbourg, France (DOI : 10.26093/cds/vizier). The original description of the VizieR service was published in 2000, A&AS 143, 23. Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the U.S. Department of Energy Office of Science. The SDSS-III web site is http://www.sdss3.org/.

Data Availability

The data used for this research are all available through the VizieR service (DOI: 10.26093/cds/vizier).