Dwarf AGNs from Variability for the Origins of Seeds (DAVOS): Properties of Variability-Selected AGNs in the COSMOS Field and Expectations for Rubin Observatory

Abstract

We study the black hole mass $-$ host galaxy stellar mass relation, $M_{\rm{BH}}-M_{\ast}$, of a sample of $z<4$ optically-variable AGNs in the COSMOS field. The parent sample of 491 COSMOS AGNs were identified by optical variability from the Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP) program. Using publicly-available catalogs and spectra, we consolidate their spectroscopic redshifts and estimate virial black hole masses using broad line widths and luminosities. We show that variability searches with deep, high precision photometry like the HSC-SSP can identity AGNs in low mass galaxies up to $z\sim 1$. However, their black holes are more massive given their host galaxy stellar masses than predicted by the local relation for active galaxies. We report that $z\sim 0.5-4$ variability-selected AGNs are meanwhile more consistent with the $M_{\rm{BH}}-M_{\ast}$ relation for local inactive early-type galaxies. This result is consistent with most previous studies of the $M_{\rm{BH}}-M_{\ast}$ relation at similar redshifts and indicates that AGNs selected from variability are not intrinsically different from the broad-line Type 1 AGN population at similar luminosities. Our results demonstrate the need for robust black hole and stellar mass estimates for intermediate-mass black hole candidates in low-mass galaxies at similar redshifts. Assuming that these results do not reflect a selection bias, they appear to be consistent with self-regulated feedback models wherein the central black hole and stars in galaxies grow in tandem.

1 Introduction

The observed local scaling relations between supermassive black hole (SMBH) mass $M_{\rm{BH}}$ and host galaxy properties (total galaxy stellar mass, bulge stellar mass, and bulge stellar velocity dispersion: $M_{\rm{BH}}-M_{\ast}$, $M_{\rm{BH}}-M_{\ast,\rm{bulge}}$, $M_{\rm{BH}}-\sigma_{\ast}$) and bulge luminosity in both active and inactive galaxies anchor our understanding of SMBH-host galaxy co-evolution (e.g., Magorrian et al. 1998; Haehnelt et al. 1998; Kormendy & Ho 2013; Reines & Volonteri 2015). These relations are usually interpreted as evidence for some form of self-regulated feedback in active galactic nuclei (AGNs). However, given that there is no consensus model of AGN feedback, observational studies of these host scaling relations across redshift and extending to lower luminosities are critical for placing new constraints on such models (e.g., Ricarte & Natarajan 2018). In addition to understanding feedback, the low-mass end of these scaling relations could also be sensitive to currently not well constrained initial mass distribution of SMBH seeds at high redshift and accretion physics (Volonteri & Natarajan, 2009; Natarajan, 2011).

Actively accreting AGNs and inactive local black hole populations appear to differ in the slope and amplitude of the $M_{\rm{BH}}-M_{\ast}$ scaling relation on which they lie (Reines & Volonteri, 2015). Scaling relations have been studied for accreting SMBH populations selected across wavelengths and redshifts. Moving beyond $z\sim 0$, Merloni et al. (2010) studied the $M_{\rm{BH}}-M_{\ast}$ relation at $z\sim 0.5-2.5$ for broad-line AGNs from the zCOSMOS survey (Lilly et al., 2007), finding evidence that the black hole to host galaxy stellar mass ratio, $M_{\rm{BH}}/M_{\ast}$, increases with redshift. Suh et al. (2020) and Zhuang et al. (2023) studied the $z\sim 0.5-2.5$ relation for X-ray selected AGNs selected from the Chandra-COSMOS Legacy Survey (Civano et al., 2016). Similarly, Ding et al. (2020) studied the $M_{\rm{BH}}-M_{\ast}$ relation for broad-line X-ray selected AGNs in deep fields. Li et al. (2021) measured the $M_{\rm{BH}}-M_{\ast}$ relation for SDSS quasars with black hole masses estimated from reverberation mapping. Mezcua et al. (2023) show the $M_{\rm{BH}}-M_{\ast}$ for seven $z\sim 0.4-0.9$ AGNs in dwarf galaxies from the VIPERS survey. With the notable exception of Suh et al. (2020), these results generally suggest that these intermediate and high redshift AGNs have over-massive black holes compared to the local ($z<0.055$) AGN relation of $M_{\rm{BH}}/M_{\ast}\sim 0.025\%$ (Reines & Volonteri, 2015). Instead, they more closely follow the relations for local inactive early-type galaxies.

Recently, exploring a higher redshift population, Pacucci et al. (2023) report that $z\sim 4-7$ quasars discovered by the James Webb Space Telescope (JWST; Harikane et al. 2023; Maiolino et al. 2023; Kocevski et al. 2023; Übler et al. 2023) have black holes masses $\sim 10-100$ times more massive compared to local AGNs with comparable stellar mass hosts. Kokorev et al. (2023) report a $z=8.5$ AGN with a $M_{\rm{BH}}/M_{\ast}$ ratio of at least $\sim 30$ percent discovered with JWST. Leveraging JWST and the Chandra X-ray Observatory, Bogdán et al. (2023); Natarajan et al. (2023); Goulding et al. (2023) discovered a $z\approx 10.1$ quasar UHZ1 with a $M_{\rm{BH}}/M_{\ast}\sim 1$. The ratio is also found to be skewed for the $z=10.6$ source GN-z11 (Maiolino et al., 2023). In fact, departure from local scaling relations at these extremely high redshifts $9<z<12$ is predicted to be a signature of the formation of heavy black hole seeds in the early Universe (Natarajan et al., 2017).

These observations of currently small samples of individual sources are strongly affected by selection biases, whereby AGNs selected only by luminosity, can produce a false evolution in host galaxy scaling relations (Lauer et al., 2007). Meanwhile, uncertainties from single-epoch virial black hole mass measurements can lead to the systematic over-estimation of black hole masses, especially at the high mass end of these relations (Shen & Kelly, 2010). Therefore, it is imperative to select AGNs with lower luminosities across redshifts in order to mitigate these selection biases (e.g., Izumi et al. 2019, 2021; Suh et al. 2020).

Optical variability is becoming an established and potentially powerful approach to identify AGNs in low luminosity sources and/or with low mass black holes in dwarf galaxies (Baldassare et al., 2018, 2020; Halevi et al., 2019; Guo et al., 2020; Burke et al., 2020, 2022, 2023; Ward et al., 2022). Besides, as noted earlier, the low mass end of these scaling relations may potentially encode additional information on black hole formation.

In this paper, we obtain black hole masses, redshifts, and reasonably robust stellar masses for variability-selected low luminosity AGNs selected by Kimura et al. (2020) from the Hyper Suprime-Cam (HSC) UltraDeep survey within the Cosmic Evolution Survey (COSMOS) field (Scoville et al., 2007). About 90 percent of the sources are detected in the X-ray (Kimura et al., 2020). Using these data, we measure the $M_{\rm{BH}}-M_{\ast}$ relation at $z\sim 0.5-3$. Our sample has a bolometric luminosity range of $L_{\rm{bol}}\sim 10^{44-47}$ erg s${}^{-1}$, comparable to the COSMOS X-ray selected sample of Suh et al. (2020). We find that this sample of $z\sim 0.5-3$ variability-selected AGNs have over-massive black holes compared to the local AGN relation of Reines & Volonteri (2015), broadly consistent with most previous studies of $z\gtrsim 0.5$ AGNs selected with other techniques (Merloni et al., 2010; Mezcua et al., 2023; Li et al., 2021; Zhang et al., 2023; Zhuang et al., 2023; Stone et al., 2023; Tanaka et al., 2024).

Our paper is organized as follows. In §2, we describe our procedure for obtaining the archival spectra and photometry and construction of a spectroscopic database for our variable AGN sample. In §3, we describe our procedure for fitting the spectral energy distribution (SED) to broad-band photometry and the reliability of our stellar mass estimates in the presence of an AGN. In §4, we describe our spectral fitting approach and resulting broad-line black hole mass estimates. In §5, we place the sources on the $M_{\rm{BH}}-M_{\ast}$ relation and compare with previous work. We discuss our results §6 and conclude in §7.

1.1 Hyper Suprime-Cam Subaru Strategic Program

In this work, we use the variability-selected AGNs from the UltraDeep Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP) program (Aihara et al., 2022) within the $\sim 1.5$ deg${}^{2}$ COSMOS field (Kimura et al., 2020). The UltraDeep HSC-SSP program has a deep single-epoch limiting magnitude of $r\sim 26$, and has therefore been used to study high redshift galaxies, supernova, and active galactic nuclei (AGNs; Kimura et al. 2020). Also see Zhong et al. (2022), who studied the morphologies of the hosts. The single-epoch photometric precision is significantly better than previous surveys like the Dark Energy Survey (DES) supernova program ($r\sim 24.5$; Kessler et al. 2015; Burke et al. 2022), and even better than what is planned for LSST Rubin without co-adding multiple exposures ($r\sim 24.5$; Ivezić et al. 2019) albeit with much a smaller survey volume. The very deep single-epoch photometry enabled identification of low luminosity variability from Type 1 AGNs up to $z\sim 4$ (Kimura et al., 2020), providing a more complete sample for detailed follow-up studies of the black hole mass – host galaxy relationships.

2 Data Analysis

2.1 Spectra Database

Table 1: Sources of optical/NIR spectroscopic redshifts for the HSC-SSP variable AGNs.

Count	Reference	Name
177	Trump et al. (2009)	Magellan XMM AGN
134	Lilly et al. (2009)	zCOSMOS 10k-bright ${}^{a}$
78	Hasinger et al. (2018)	DEIMOS 10K ${}^{a}$
64	Cool et al. (2013)	PRIMUS ${}^{a}$
45	Hasinger et al. (2018)	DEIMOS 10K
45	Ahumada et al. (2020)	SDSS DR16 ${}^{a}$
35	Silverman et al. (2015)	FMOS COSMOS DR2${}^{a}$
31	Lilly et al. (2007)	zCOSMOS DR3
17	Marchesi et al. (2016)	Chandra COSMOS legacy
14	Momcheva et al. (2016)	3D-HST v4.1.5 ${}^{a}$
11	van der Wel et al. (2021)	LEGA-C DR3
6	Silverman et al. (2015)	FMOS-COSMOS
6	Straatman et al. (2018)	LEGA-C DR2 ${}^{a}$
5	Schulze et al. (2018)	FMOS-COSMOS AGN NIR
4	Damjanov et al. (2018)	hCOSMOS
4	Alam et al. (2015)	SDSS DR12
4	Pâris et al. (2014)	SDSS DR10 quasar cat
3	DESI Collaboration et al. (2023)	DESI EDR
3	Balogh et al. (2014)	GEEC2
3	Kartaltepe et al. (2015)	FMOS-COSMOS NIR
3	Monzon et al. (2020)	CLAMATO
2	Pâris et al. (2018)	SDSS DR14 quasar cat
2	Boutsia et al. (2018)	IMACS faint AGN
2	Le Fèvre et al. (2013)	VVDS DRFinal ${}^{a}$
1	Ahn et al. (2012)	SDSS DR9
1	Knobel et al. (2012)	zCOSMOS 20k Group
1	Brusa et al. (2009)	z$>$3 X-ray QSOs
1	Allevato et al. (2012)	z$<$1 X-ray AGN
1	Jamal et al. (2018)	VVDS reprocessed
1	Onodera et al. (2016)	Star-forming Galaxies NIR
1	Harrison et al. (2016)	KASHz
1	Koprowski et al. (2016)	SCUBA-2
1	Harrison et al. (2017)	KROSS
1	Straatman et al. (2018)	LEGA-C DR2
1	Ono et al. (2018)	GOLDRUSH
1	Masters et al. (2019)	C3R2 DR2
1	Rosani et al. (2020)	SMUVS Lya emitters
1	Mukae et al. (2020)	HETDEX LAEs/eBOSS QSOs
1	Lyke et al. (2020)	SDSS DR16 QSOs
1	Stanford et al. (2021)	C3R2 DR3
1	Masters et al. (2017)	C3R2 DR1 ${}^{a}$

Note. — ($a$) Data taken from the HSC PDR3 catalog of spectroscopic redshifts.

Figure 1: Redshift from Kimura et al. (2020) versus our updated redshifts from the HSC PDR3 redshift catalog or SIMBAD (blue circle symbols and orange triangle symbols). Photometric redshifts are from the COSMOS2020 catalog (gray square symbols; Weaver et al. 2022).

Figure 2: Bolometric AGN luminosity (calculated from 10 $\times$ the 2$-$10 keV X-ray luminosity; e.g. Duras et al. 2020) versus our updated redshifts. The corresponding redshift and bolometric luminosity distributions are shown.

Figure 3: Comparison of COSMOS2020 (Weaver et al., 2022) photometric redshifts with LePHARE (left) and EAZY (right) against the spectroscopic redshifts in this work. The LePHARE photometric redshifts includes AGN SED templates, and are therefore more appropriate for the HSC-SSP variable AGNs. The root mean square error (RMSE) values are shown in the upper left hand corner of both figure panels. The gray $y=x$ line is shown for comparison.

Kimura et al. (2020) have compiled spectroscopic redshifts for their HSC-SSP variable AGNs from the HSC ancillary data products. Specifically, they matched their variable AGNs to the HSC public data release 2 (PDR2) catalog of spectroscopic redshifts, and spectroscopic redshifts from the Chandra-COSMOS Legacy Survey and DEIMOS catalogs (Hasinger et al., 2018). This database includes spectroscopic redshifts overlapping with the HSC-SSP COSMOS field from zCOSMOS DR3 (Lilly et al., 2009), 3D-HST (Skelton et al., 2014; Momcheva et al., 2016), FMOS-COSMOS (Silverman et al., 2015; Kashino et al., 2019), VVDS (Le Fèvre et al., 2013), SDSS DR12 (Alam et al., 2015), PRIMUS DR1 (Coil et al., 2011; Cool et al., 2013). Kimura et al. (2020) also matched to the DEIMOS 10k Spectroscopic Survey Catalog (Hasinger et al., 2018), which was not included in the HSC PDR2 catalog of spectroscopic redshifts.

We have instead matched to the HSC PDR3 version of the same catalog, which includes updated SDSS redshifts from SDSS DR15 and additional spectroscopic redshifts overlapping with the HSC-SSP COSMOS field from the DEIMOS 10k sample (Hasinger et al., 2018) and LEGA-C DR2 (Straatman et al., 2018). Using an astroquery (Ginsburg et al., 2019) SIMBAD search, we identified additional spectroscopic redshifts from the COSMOS XMM-Newton AGN spectroscopic survey (Trump et al., 2009) and a several additional surveys in the literature that have been indexed by SIMBAD (Lilly et al., 2007; Brusa et al., 2009; Trump et al., 2009; Knobel et al., 2012; Allevato et al., 2012; Ahn et al., 2012; Pâris et al., 2014; Balogh et al., 2014; Kartaltepe et al., 2015; Alam et al., 2015; Silverman et al., 2015; Onodera et al., 2016; Harrison et al., 2016; Koprowski et al., 2016; Harrison et al., 2017; Jamal et al., 2018; Pâris et al., 2018; Hasinger et al., 2018; Boutsia et al., 2018; Damjanov et al., 2018; Schulze et al., 2018; Straatman et al., 2018; Ono et al., 2018; Masters et al., 2019; Rosani et al., 2020; Monzon et al., 2020; Mukae et al., 2020; Lyke et al., 2020; Stanford et al., 2021; van der Wel et al., 2021). We restrict the SIMBAD search to spectroscopic redshifts derived from reliable optical or near-infrared (NIR) spectra using the flags ${\rm RVZ\_WAVELENGTH}=={{}^{\prime}O^{\prime}}$ or ${\rm RVZ\_WAVELENGTH}=={{}^{\prime}N^{\prime}}$ and ${\rm rvz\_qual}!={{}^{\prime}E^{\prime}}$. The first constraint restricts the search to optical or NIR spectra. The second constraint excludes photometric redshifts. Finally, we download the publicly-available spectra from these sources (see Appendix A). We also obtained spectroscopic redshifts from the Dark Energy Spectroscopic Instrument (DESI) Early Data Release (DESI Collaboration et al., 2023), which have not yet been indexed by SIMBAD at the time of writing. We choose our fiducial spectroscopic redshifts according to the following priority: SIMBAD (316 sources), HSC PDR3 catalog of spectroscopic redshifts (93 sources), DESI (3 sources), and Chandra-COSMOS Legacy Survey spectroscopic redshifts (17 sources). But significantly inconsistent redshifts are resolved by hand after careful visual inspection (described in §2.2). If no spectroscopic redshift was found, we use the photometric redshifts from the COSMOS2020 catalog (50 sources; §2.3). There are also 12 sources with no COSMOS2020 match, and therefore with no redshift. Our updated redshifts are shown in Figure 1. The bolometric luminosity and redshift distributions are shown in Figure 2. A table showing the sources of the public spectroscopic redshifts for the HSC-SSP AGNs is shown in Table 1, and our updated redshifts are given in Table 3. To include cases where a single source has more than one available spectrum from different programs, we always repeat the matching between the HSC-SSP AGNs and the spectroscopic sample when downloading the spectra from publicly-available sources.

2.2 Inconsistent spectroscopic redshifts

Using $\Delta z>0.1$ as our criterion, we found 36 sources with inconsistent spectroscopic redshifts between SIMBAD, HSC PDR3, DESI, and Chandra-COSMOS Legacy Survey catalogs and 20 inconsistent optical spectra between Magellan, zCOSMOS, DEIMOS, DESI, and LEGA-C spectral files, from which we will use to estimate BH masses. For each source with inconsistent spectroscopic redshifts, we plotted their publicly available spectra against a Type 1 quasar template spectrum (e.g., Vanden Berk et al. 2001; Richards et al. 2006) and attempted to identify which spectroscopic redshift is correct by visually matching the spectrum with the strong emission lines from the template. At least two of the authors have visually inspected and agreed to these corrected redshifts. In most cases the correct spectroscopic redshift can be obviously identified. For 7 of these sources, we were unable to determine the correct spectroscopic redshift due to poor spectral quality or lack of obvious emission lines or features in the spectrum. We adopt the COSMOS2020 photometric redshifts for these 7 sources. 28 of the clearly incorrect spectroscopic redshifts originate from the PRIMUS catalog. The PRIMUS spectra are not publicly available, preventing us from visually inspecting them. We suspect these incorrect PRIMUS redshifts are due to a misidentified emission line given the relatively low resolution of the prism spectra of $R=\lambda/\Delta\lambda\sim 40$ (Coil et al., 2011).

2.3 Photometric redshifts

For sources without a reliable spectroscopic redshift, we use the LePHARE (Arnouts et al., 1999; Ilbert et al., 2006) photometric redshifts from the COSMOS2020 catalog (Weaver et al., 2022). The LePHARE redshifts are derived by fitting galaxy and AGN templates to the photometry. When compared to EAZY (Brammer et al., 2008) photometric redshifts without AGN templates, we find that the LePHARE redshifts are generally much more reliable. Using the spectroscopic redshifts as our ground truth, we estimate the scatter in the photometric redshifts as $\sigma=1.4826$ MAD, where MAD is the median absolute deviation, which is robust to outliers (e.g., Burke et al. 2022). We find $\sigma=0.43$ for the LePHARE photometric redshifts and $\sigma=0.53$ for the EAZY photometric redshifts. A comparison between the two photometric redshifts for the HSC-SSP variable AGNs are shown in Figure 3.

3 SED fitting

Table 2: X-CIGALE SED fitting parameters.

Module	Parameter	Values
Star formation history: delayed model, $\mathrm{SFR}\propto t\exp(-t/\tau)$	$e$-folding time, $\tau$ (Gyr)	0.1, 0.5, 1, 5
	Stellar age, $t$ (Gyr)	0.5, 1, 3, 5, 7
Simple stellar population: Bruzual & Charlot (2003)	Initial mass function	Chabrier (2003)
	Metallicity, $Z$	0.02
Galactic dust attenuation: Calzetti et al. (2000) & Leitherer et al. (2002)	$E(B-V)$ of starlight for the young population	0.0, 0.5, 1.0, 1.5
	$E(B-V)$ ratio between the old and young populations	0.44
Galactic dust emission: Dale et al. (2014)	$\alpha$ slope in $dM_{\rm dust}\propto U^{-\alpha}dU$	1.5, 2.0, 2.5
AGN (UV–IR): SKIRTOR	Torus optical depth at 9.7 microns $\tau_{9.7}$	7.0
	Torus density radial parameter $p$ ($\rho\propto r^{-p}\mathrm{e}^{-q|\cos(\theta)|}$)	1.0
	Torus density angular parameter $q$ ($\rho\propto r^{-p}\mathrm{e}^{-q|\cos(\theta)|}$)	1.0
	Angle between the equatorial plan and edge of the torus	40${}^{\circ}$
	Ratio of the maximum to minimum radii of the torus	20
	Viewing angle $\theta$ (face on: $\theta=0^{\circ}$, edge on: $\theta=90^{\circ}$)	30${}^{\circ}$ (type 1)
	Power-law index of the disk UV/optical slope $\delta$	$-1.0$, $-0.5$, $-0.36$, 0.0, 0.5, 1.0
	UV AGN fraction $0.1-0.3\ \mu$m $f_{\rm{AGN}}$	$0-0.9$ (step 0.1), 0.95, 0.9999
	Extinction law of polar dust	SMC
	$E(B-V)$ of polar dust	0.0, 0.1, 0.3, 0.5
	Temperature of polar dust (K)	100
	Emissivity of polar dust	1.6
X-ray:	AGN photon index $\Gamma$	1.8
	Maximum deviation from the $\alpha_{\rm{ox}}$-$L_{\rm{2500}\AA}$ relation	0.2

Note. — See Yang et al. (2020, 2022) for details of x-cigale SED parameters and models.

We obtained deep, broad-band photometry by matching to version 2.2 of the COSMOS2020 catalog (Weaver et al., 2022). We identified 455 matches with the COSMOS2020 catalog. We use the photometry derived from The Farmer pipeline (Weaver et al., 2023), which is based on source modeling using The Tractor software (Lang et al., 2016). We use included deblended photometry extracted from HSC-SSP PDR2 optical imaging, UltraVISTA near-infrared imaging (McCracken et al., 2012; Moneti et al., 2023), Spitzer/IRAC mid-infrared images from the Cosmic Dawn Survey (Moneti et al., 2022), and near and far ultraviolet imaging from the COSMOS GALEX catalog (Zamojski et al., 2007). The X-ray 0.5$-$10 keV photometry is from the Chandra COSMOS Legacy Survey (Civano et al., 2016; Marchesi et al., 2016). Inclusion of the optical to mid-infrared photometry is essential to constrain the star-formation and reprocessed dust emission. Additionally, including the X-ray photometry or even upper limits can help constrain the AGN emission (Yang et al., 2020), helping to eliminate some degeneracies between star-formation and AGN emission which can lead to spurious stellar mass estimates.

We use version 2022.1 of the cigale code (Burgarella et al., 2005; Noll et al., 2009; Boquien et al., 2019; Yang et al., 2020, 2022) to perform SED fitting. This version of the code includes X-ray fitting modules from x-cigale, with detailed AGN emission models that have been extensively tested with both galaxy and AGNs (Yang et al., 2020, 2022). These codes work by imposing a self-consistent energy balance constraint between different emission and absorption mechanisms across the EM spectrum. A large grid of models is computed and fitted to the data, allowing for an estimation of the star formation rate (SFR), stellar mass, and AGN contribution via a Bayesian-like analysis of the likelihood distribution.

We use a delayed exponential star formation history and vary the $e$-folding time and age of the stellar population assuming solar metallicity. Nevertheless, different choices of the initial mass function, stellar population models, and star formation histories can introduce systematic uncertainties of $\sim 0.3$ dex (Conroy, 2013). Zou et al. (2022) found that different parametric star formation histories results in systematic differences in stellar mass of $\sim 0.1$ dex for a sample of $z=0-6$ AGN. We adopt the commonly-used Chabrier (2003) initial stellar mass function with the stellar population models of Bruzual & Charlot (2003) and adopt the nebular emission template of Inoue (2011). We use the Leitherer et al. (2002) extension of the Calzetti et al. (2000) model for reddening due to dust extinction, and the Draine et al. (2014) updates to the Draine & Li (2007) model for dust emission. Finally, we adopt the SKIRTOR clumpy two-phase torus AGN emission model (Stalevski et al., 2012, 2016) allowing for additional polar extinction. We do not force the UV part of the SED to be completely dominated by the AGN component, because this can result in over-estimation of the stellar masses. We assume Type-1-like inclination angle of $30$ deg for the HSC-SSP variable AGNs. Our choice of considering only a single viewing angle close to the average values for Type 1 AGNs is justified by previous studies, which find that different viewing angles were largely degenerate with the average values of 30 and 70 degrees for Type 1 and 2 AGNs, respectively (e.g., Mountrichas et al., 2021; Ramos Padilla et al., 2022).

Table 2 shows our cigale input parameters, with example SED fitting results shown in Figure 4 (Yang et al., 2020). Our stellar masses range from $M_{\ast}\sim 10^{9-11.5}\ M_{\odot}$. Inclusion of far-infrared photometry could further improve the results by constraining the dust emission component and reducing degeneracies in the modeling, given cigale’s energy conservation principle. Unfortunately, there are difficulties with spatially associating sources given the large PSF sizes in the far-infrared. Our stellar masses could also be refined using a future COSMOS2020 super-deblended catalog with far-infrared photometry (Jin et al. in preparation).

3.1 Reliability of stellar mass estimates

Figure 4: Left panel: Example SED fitting result for an AGN-dominated source (Kimura et al. (2020) ID = 2). Including a stellar emission component does not significantly improve the fit in this source. Therefore, the stellar emission and associated parameters (e.g., stellar mass, star formation history) cannot be reliably constrained. Right panel: Example SED fitting result for a non AGN-dominated source (Kimura et al. (2020) ID = 3). The inclusion of a stellar emission component improves the best-fit $\chi_{\nu}^{2}$, and the resulting best-fit stellar emission is a significant excess of the AGN continuum emission near $\sim 1.2$ $\mu$m, enabling a reliable stellar mass estimate. The uncertainties on the stellar mass here are estimated by cigale and do not include uncertainties beyond the parameter choices in Table 2. The vertical red line is the region near $1.2$ $\mu$m, where the AGN emission is at a minimum, and the red number is the ratio of SF emission (blue) from the SF$+$AGN model to total emission (black) of an AGN-dominated model (gray). Both redshifts are spectroscopic.

Figure 5: Results of applying our stellar mass estimation technique to photometry of the the mostly host-dominated local AGN sample of Reines & Volonteri (2015). Our stellar masses ($y$-axis) are in very good agreement with those of Reines & Volonteri (2015), estimated from mass-to-light ratios after subtracting the AGN contribution from the photometry ($x$-axis). The dotted gray lines of $\pm 0.3$ dex are shown to guide the eye about the solid gray $y=x$ line. All stellar masses are reliably estimated following our procedure This figure demonstrates that our reliable stellar masses are reasonable when the AGN contribution does not dominate the SED.

Figure 6: Validation of the stellar masses from COSMOS2020 broad-band SED fitting ($x$-axis; this work) to stellar masses from COSMOS-Webb and HST archival imaging ($y$-axis; Zhuang et al. 2023). The blue points are sources with reliable COSMOS2020-based SED stellar masses, and the gray points have unreliable stellar masses. The stellar masses from COSMOS-Webb and HST imaging are based on SED fitting after spatially decomposing the AGN+host emission. Our stellar masses are consistent to $\sim 0.3$ dex (1$\sigma$ scatter, shown as dashed gray lines about the solid gray $y=x$ line), which is comparable to the expected systematic error of $\sim 0.3$ dex on SED-based stellar masses (Conroy, 2013). This figure demonstrates that our reliable stellar masses are reasonable, even when compared to those measured from AGN+host decomposition with space-based imaging. The good agreement for our sample could partially be due to the larger number of sources with lower AGN luminosities (compared to the host) at a given stellar mass. Our sample of HSC-SSP AGNs includes many more sources outside the COSMOS-Webb and HST archival imaging footprint.

Table 3: Properties of COSMOS variable AGNs from SED and spectral fitting analysis.

ID	RA	DEC	i-mag	$z_{{\rm best}}$	$z_{{\rm ph}}$	$\log L$	$\log M_{\rm{BH}}$	$\log M_{\rm{BH,err}}$	$\log M_{\ast}$	$\log M_{\ast,{\rm{err}}}$	${\rm{SF}}_{\rm{ex}}$	$\chi^{2}_{\nu}$	class
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)
	deg	deg	AB mag			$\log$ erg s${}^{-1}$	$\log M_{\odot}$	$\log M_{\odot}$	$\log M_{\odot}$	$\log M_{\odot}$
1	150.74386	2.20245	22.71	1.58	1.58				10.64	0.12	1.1	3.0
2	150.73557	2.19957	20.36	3.5026	-1.0	46.68	9.04	0.11	11.53	0.13	0.8	1.9	BL
3	150.73353	2.15646	20.88	0.979	1.03				11.4	0.11	2.1	0.7	NeV
4	150.79702	2.13888	21.01	0.5727	0.51				10.78	0.2	3.5	1.3
⋮
491	150.03524	2.72781	21.04	0.5093	0.52				10.49	0.06	5.4	0.7

Note. Column (1): Identifier from Table 4 of Kimura et al. (2020). Column (2): RA. Column (3): Dec. Column (4): $i$ band AB magnitude from Kimura et al. (2020). Column (5): Best redshift. Column (6): Photometric redshift from COSMOS2020 LePHARE. Column (7): Bolometric luminosity from the spectral continuum. Column (8): Virial black hole mass. Column (9): Inferred cigale stellar mass. Column (10): Excess SF over an AGN dominated model (if $>1.2$, stellar masses are considered reliable). Column (11): Best-fit reduced cigale model $\chi^{2}$ (recommend $<5$). Column (12): Visual classification of spectrum (see §4.3). All uncertainties are $1\sigma$ statistical errors from fitting. This table is published in its entirety in the published version. Only a portion is shown here.

For galaxies with non-negligible AGN emission, the AGN emission must be constrained or subtracted in order to properly model the star-formation emission in the galaxy and obtain a reliable stellar mass. High resolution optical/NIR imaging combined with source profile fitting can be used to subtract the AGN point source emission from the underlying host galaxy. When high resolution imaging is not available or feasible to analyze in large quantities, SED fitting to broad-band catalog photometry can be used with some caveats.

The UV/optical emission from an unobscured AGN accretion disk can be strongly degenerate with stellar emission. When the SED is dominated by an unobscured AGN, the stellar emission component is swamped by the AGN emission and the resulting stellar emission parameters (e.g., stellar mass, star formation history) cannot be reliably constrained (e.g., Merloni et al. 2010; Ciesla et al. 2015). We employ a model comparison technique to determine whether the stellar emission component, and by extension stellar mass, can be constrained by SED fitting. The model comparison test works as follows. First, we fit the SED using an AGN-dominated model by setting $f_{\rm{AGN}}=0.9999$, where $f_{\rm{AGN}}$ is the AGN fraction computed between observed-frame $0.5-1$ $\mu$m. This wavelength is where the AGN emission is near a minimum, assuring the overall SED is totally dominated by the AGN emission. For technical reasons, the AGN fraction cannot be set to exactly unity in the cigale code (Yang et al., 2020). We have rounded the AGN fraction up to 1 for clarity in the presentation of this paper. Then, we fit the SED using a mixed AGN+stellar emission model by allowing the AGN emission to vary, only requiring the UV/blue optical emission to be AGN-dominated, i.e., $f_{\rm{AGN}}>0.9$ between observed-frame 0.1$-$0.3 $\mu$m (Suh et al., 2020). This is justified by evidence that most of our spectra show significant AGN continuum emission at these wavelengths, and the fact that the sources have been identified from optical band light curves. Optical variability, presumably arising primarily from variations in the AGN continuum luminosity, tends to be stronger in the bluer bands as well. Finally, we consider the stellar masses reliable for only those SEDs with significant stellar emission at rest-frame 1.2 $\mu$m that cannot be explained by a totally AGN-dominated model.

To determine whether the stellar emission is constrained, we compute the ratio of the total model fit for the AGN+stellar emission model and the AGN-dominated model. We consider the stellar masses reliable if this ratio is greater than 1.2, as justified in Appendix C. This method is a combination of similar approaches used in the literature (Merloni et al., 2010; Suh et al., 2020; Burke et al., 2022). Our approach eliminates SEDs that are totally degenerate with AGN-dominated emission at rest-frame 1.2 $\mu$m, where the AGN emission is close to a minimum. We found that simply comparing the reduced $\chi^{2}$ values was less reliable, because $\chi^{2}$ values are sensitive to the photometric uncertainties on the data, model parameter choices (i.e., over-fitting stellar emission), and does not always indicate an improved fit in the optical-NIR region where the stellar emission is strongest and from where the stellar masses are derived. Example SEDs from our sample with reliable and unreliable stellar mass estimates are shown in Figure 4. The scatter in the recovered stellar masses is typically $\sim 0.2$ dex. This uncertainty is fully not taken into account in the cigale code, which tends to under-estimate the uncertainties in stellar mass. We add 0.2 dex in quadrature to our cigale uncertainties throughout the figures in this paper.

3.2 Validation of stellar mass estimation method against the local AGN sample

Reines & Volonteri (2015) measured the $M_{\rm{BH}}-M_{\ast}$ relation for a sample of $z<0.055$ AGNs selected from SDSS optical spectroscopy. Their stellar masses are estimated using mass-to-light ratios as function of host galaxy color (Zibetti et al., 2009). They estimated the AGN contribution to the integrated photometry by re-scaling a mock AGN spectrum to match the continuum luminosity, then convolving the scaled mock spectrum with the SDSS filter throughput curves. They then subtracted the contribution from AGN emission from the photometry before applying the mass-to-light ratios. Five sources were excluded from their analysis due to being dominated by AGN emission at the $\gtrsim 50\%$ level.

As a consistency check, we ran our SED fitting procedure on the Reines & Volonteri (2015) sample of local AGNs using Petrosian flux photometry from the NASA Sloan Atlas catalog of GALEX UV and SDSS optical photometry (Blanton et al., 2011), shown in Figure 5. All of the stellar masses are reliably estimated with our cigale procedure and we find a very good agreement with their results across the entire range of stellar masses. Of course, the photometry of the local AGN sample are comparatively host-dominated compared to our HSC-SSP AGNs. Nevertheless, this demonstrates that our stellar masses are highly reliable when the AGN contribution does not completely dominate the SED. This also ensures that discrepancies between the local $M_{\rm{BH}}-M_{\ast}$ relation and those at higher redshifts are not due to errors in the stellar mass measurements of local AGNs.

3.3 Validation of stellar masses using JWST and HST spatially-decomposed and PSF-subtracted stellar masses

To validate our stellar masses in cases where the AGN contribution to the SED is more substantial, we compare our stellar masses to those estimated by Zhuang et al. (2023). These authors measured stellar masses for broad-line X-ray selected AGNs with BH masses from Suh et al. (2020) using archival HST and JWST COSMOS-Webb imaging. Their stellar masses from COSMOS-Webb and HST imaging are based on SED fitting with cigale after spatially decomposing the AGN+host emission. Our stellar masses are consistent to $\sim 0.4$ dex (1$\sigma$ scatter, shown as dashed gray lines about the solid gray $y=x$ line) from Figure 6, which is comparable to the expected systematic error of $\sim 0.3$ dex on SED-based stellar masses (Conroy, 2013). The single outlier in the lower left hand corner of the plot has an excess SF model ratio of 1.33 from our SED fitting reliability approach, hence, only a very slight bump near 1.2 $\mu$m from stellar emission. Both sources have inconsistent photometric and spectroscopic redshifts with no obvious emission line features in their spectra. No significant systematic offset is found between our stellar masses and those of Zhuang et al. (2023). This is encouraging, as it indicates that our COSMOS2020-based stellar masses are not strongly biased due to AGN contamination. We attribute this to our model-comparison technique for selection of reliable stellar masses based on the strength of the stellar emission relative to the AGN emission. Of course, our stellar masses are not as robust as those measured from high resolution HST and JWST imaging.

3.4 Detection limits

Figure 7: Host galaxy stellar mass versus redshift for HSC-SSP variable AGNs with reliable stellar mass estimates from broad-band SED fitting using cigale. Each AGN is shaded by its HSC $i$-band apparent magnitude. Symbols with black borders have spectroscopic redshifts from our redshift database. Symbols without black borders have photometric redshifts. The gray curve is the theoretically predicted detection limit.

In Figure 7, we show our derived stellar mass estimates versus redshift for the HSC-SSP AGNs. The gray curves show the theoretically-predicted BH mass detection limits following Burke et al. (2023). These BH mass horizon curves are computed assuming a limiting detectable variability amplitude of 0.1 magnitudes and the modified correlations between optical variability amplitude and BH mass (e.g. MacLeod et al. 2010) as described in Burke et al. (2023). The predicted BH mass detection limits are derived assuming a typical amplitude of variability of 0.1 mag and a photometric precision given by Ivezić et al. (2019) with limiting magnitude of $m_{5}=25.5$.

4 Spectral Measurements

4.1 Flux calibration and data cleaning

Figure 8: Demonstration of absolute flux calibration of our spectra. Each data point is a single HSC SSP AGN with multiple spectra—one from zCOSMOS and at least one from another source. The bolometric luminosity from zCOSMOS spectral continua is plotted against the bolometric luminosity from other spectra (see Table 1) before (left) and after (right) performing absolute flux calibration to the COSMOS2020 photometry. Our resulting flux-calibrated spectra have no systematic offset with zCOSMOS and the fluxes are more consistent as evident by the significantly reduced scatter. The right panel only includes those with successful flux calibration.

We downloaded publicly-available spectra for our parent sample of HSC-SSP AGNs from Table 1. In total, we have nearly 1000 spectra of varying quality and wavelength coverage. Below, we describe our procedure for calibrating and cleaning the data. A reasonable estimate of the uncertainties on the data is essential for the least squares minimization and $\chi^{2}$ estimation. Unfortunately, the flux uncertainties (error spectrum) are not provided for the DEIMOS and zCOSMOS data. In these cases where the error spectrum is not provided, we estimate the uncertainties using the median absolute deviation of the flux spectrum. Data quality (e.g. artifact, spectral gap) masks are used wherever possible. The fitting results of each spectrum is visually inspected. When multiple spectra are available for the same source, both with a valid BH mass estimate, we pick the spectrum with the highest signal-to-noise ratio.

Accurate flux calibration is important for obtaining virial BH masses, which depends on the continuum or broad-line luminosity. In order to minimize systemic differences between the spectra from different surveys and instruments, we first perform absolute flux calibration for each spectrum. We integrate the spectrum over the available COSMOS2020 bands to generate synthetic photometry. Then, we scale by the error-weighted mean ratio between synthetic photometry and HSC photometry. This procedure is commonly adopted in the literature (e.g. Mallery et al. 2012). Figure 8 demonstrates a significantly improved consistency between the spectral calibration for sources with multiple spectra from different surveys/instruments after performing absolute flux calibration with the COSMOS2020 photometry. Other issues in the data reduction, such as residual instrumental sensitivity, are difficult to correct. In addition, we have not corrected for the effect of variability in the spectra. We estimate a 1$\sigma$ uncertainty of $\sim 10$ percent in the final absolute flux calibration from the scatter in Figure 8.

4.2 Spectral modeling

We fit the continuum and emission lines in each 1D spectrum using a modified version of the publicly-available PyQSOFit code (Guo et al., 2018; Shen et al., 2019), used to measure SDSS quasar properties (Shen et al., 2011; Wu & Shen, 2022). However, unlike the SDSS quasar sample, we found that many of our spectra had a significant underlying host galaxy component. This is not surprising given our much fainter quasar sample, which tend to be more host dominated, compared to the magnitude limit of SDSS quasar spectra of $i<20$. Therefore, we first perform a quasar/host galaxy decomposition using principle component analysis (PCA) with host galaxy templates (Bruzual & Charlot, 2003). After subtracting any significant host galaxy component from the spectrum, the quasar continuum is modeled as a blue power-law plus a 3rd-order polynomial for reddening. Fe II emission templates (Vestergaard & Wilkes, 2001) are fitted if including them improves the reduced $\chi^{2}$ of the continuum fitting by 20 percent. The total model is a linear combination of the continuum and single or multiple Gaussians for the emission lines. Since uncertainties in the continuum model may induce subtle effects on measurements for weak emission lines, we first perform a global fit to the emission-line free region to better quantify the continuum. We then fit multiple Gaussian models to the continuum-subtracted spectra around the H$\beta$ and Mg II emission line complex regions locally.

Variability is efficient at identifying low luminosity AGNs. Generally speaking, the host galaxy continuum component is significant in the spectra HSC-SSP variable AGNs. It is difficult to correctly decompose the underlying AGN continuum from the host galaxy given the lower AGN luminosities and lower continuum signal-to-noise ratios of the spectra owing to their higher redshifts and faintness. For this reason, we chose to use the broad-line luminosities instead of AGN continuum luminosities to estimate BH masses.

4.3 Visual Spectral Classification

In addition to our automated broad-line detection and fitting approach, we have visually inspected each spectrum to identify AGN signatures and present initial classifications. We classify sources with at least one broad line feature (“broad line”), sources without strong broad-line features but with Ne V emission line indicating a probable AGN (“Ne V”), and sources without either AGN feature (“host-dominated”), or low $S/N$ spectrum (“noisy”). Due to varying quality of the spectrophotometric calibrations, we do not attempt to identify AGN continuum features for non broad-line/Ne V emitting sources. We additionally identify J095835.9+015156 (ID=280) as a broad absorption line (“BAL”) quasar. 252 of our sources have AGN features (broad line or Ne V line), leaving 128 without obvious AGN-like features in their spectra (host-dominated or noisy). The classifications are presented in Table 3. We caution that the absence of AGN spectral features that we have identified does not necessarily imply the absence of an AGN. For example, the spectrum could be in a host-dominated state due to variability, or the spectrum simply does not cover any of the broad emission lines or the Ne V line given the redshift.

4.4 Black Hole Masses

Following Shen et al. (2011), we estimate the BH masses from broad-emission lines (e.g., Greene & Ho 2005) using the single-epoch virial method. This method assumes that the broad-line region (BLR) is virialized and uses the continuum luminosity and broad-line FWHM as a proxy for the BLR radius and virial velocity respectively. Under these assumptions, the BH mass can be estimated by:

$$\begin{split}\log{\left(\frac{M_{\rm{BH}}}{M_{\odot}}\right)}=a+b&\log{\left(\frac{L_{\rm{br}}}{10^{44}\ \rm{erg\ s}^{-1}}\right)}\\ &+2\log{\left(\frac{\rm{FWHM}_{\rm{br}}}{\rm{km\ s}^{-1}}\right)}\end{split}$$ (1)

where $L_{\rm{br}}$ and FWHM${}_{\rm{br}}$ are the broad-line luminosity and full-width-at-half-maximum (FWHM) with an intrinsic scatter of $\sim 0.4$ dex in BH mass. The coefficients $a$ and $b$ are empirically calibrated against local AGNs with BH masses measured from reverberation mapping. We adopt the calibrations (Vestergaard & Peterson, 2006) from in Shen et al. (2011) derived by Shaw et al. (2012):

$$(a,b)=(1.63,0.49),\quad\rm H\beta$$ (2)

$$(a,b)=(1.70,0.63),\quad\rm Mg\ II$$ (3)

$$(a,b)=(1.52,0.46),\quad\rm C\ IV.$$ (4)

Following Shen et al. (2011), depending on the wavelength coverage, redshift, and line $S/N$, we adopt a fiducial or preferred “best” BH mass following the ordering above. We only consider BH masses when the broad line component is detected over the residual with a $S/N>2$ as defined in Burke et al. (2024). When more than one spectrum exists with a valid BH mass for a single source, we adopt the spectrum with the highest median per-pixel $S/N$.

4.5 X-ray Properties

Figure 9: 2$-$10 keV X-ray luminosities against our measured broad-line black hole masses. The points are colored by their stellar mass when deemed reliable. We assume an uncertainty of $0.4$ dex on the black hole masses.

Kimura et al. (2020) show that $\sim 90$ percent of their sources are detected in the X-ray. For non X-ray detected sources, their X-ray stacking analysis favors an AGN emission origin for the X-ray stacked sample. We have re-done the matching and plot the 2$-$10 keV X-ray luminosities from the Chandra-COSMOS Legacy Survey catalog (Civano et al., 2016) against our broad-line black hole masses in Figure 9. Kimura et al. (2020) demonstrated that the X-ray luminosities are too high to be explained by X-ray binary populations. Assuming a typical bolometric correction of $L_{\rm{bol}}/L_{2-10\ \rm{{keV}}}=10$ (Duras et al., 2020), the sources have a typical (median) Eddington ratio of $\sim 0.05$, slightly lower than the median Eddington ratio of $\sim 0.1$ of the sample from Suh et al. (2020). Our estimated Eddington ratios and black hole masses are broadly reasonable given the correlation between these parameters and optical variability amplitude (e.g. MacLeod et al. 2010) and the typical variability amplitude of the HSC-SSP AGNs of $\sim 0.1$ magitudes.

5 $M_{\rm{BH}}-M_{\ast}$ relation

Figure 10: $M_{\rm{BH}}-M_{\ast}$ relation for the HSC-SSP AGNs with reliable BH mass estimates from this work, and reliable stellar masses from this work (left) or from AGN+host decomposition from HST and COSMOS-Webb archival imaging from Zhuang et al. (2023) (right). Points are shaded by redshift. We assume uncertainties of $\sim 0.4$ dex on the BH masses and $\sim 0.4$ dex on the COSMOS-Webb stellar masses. The local $M_{\rm{BH}}-M_{\ast}$ relation for AGNs is plotted as the solid blue line (Reines & Volonteri, 2015) and for inactive elliptical galaxies as the dash-dotted orange line (Greene et al., 2020). Lines of constant $M_{\ast}=10M_{\rm{BH}}$ and $M_{\ast}=100M_{\rm{BH}}$ are plotted in dashed and dotted gray lines, respectively.

Figure 11: Black hole mass versus redshift for HSC-SSP variable AGNs with reliable broad-line BH mass estimates from broad-band SED fitting using cigale. Each AGN is shaded by its HSC $i$-band apparent magnitude. The gray curve is the theoretically predicted detection limit, which is deeper than the BH mass detection limit.

We show the $M_{\rm{BH}}-M_{\ast}$ relation for our sample of variability-selected AGNs with reliable stellar masses from COSMOS2020 photometry and virial black hole masses measured from our spectroscopic database in Figure 10. We also show the $M_{\rm{BH}}-M_{\ast}$ relation using the HST and JWST stellar masses from Zhuang et al. (2023) and virial black hole masses measured from our spectroscopic database. In either case, the $M_{\rm{BH}}-M_{\ast}$ relation for our sample of $z\sim 0.5-4$ (median redshift of $\sim 1.5$) AGNs is more consistent with the relation for local inactive elliptical galaxies (Greene et al., 2020) than local ($z<0.055$) AGNs (Reines & Volonteri, 2015). However, the variability-selected AGNs show a possible trend toward the local AGN relation at lower redshifts. In contrast, the Zhuang et al. (2023) sources are highly consistent with the relation for local inactive elliptical galaxies. This difference could be due to having more $z\lesssim 1$ and lower luminosity AGNs (with relatively more host-dominated SEDs) in our figure panel.

At face value, this result could indicate an evolution in the $M_{\rm{BH}}-M_{\ast}$ relation for AGNs with redshift. Such an evolution could be connected to the host star-formation or black hole accretion activity (the densities rates of both peak near $z\sim 2$). On the other hand, Bongiorno et al. (2012) studied the star-formation rates and AGN activity in a sample of AGNs selected from X-ray and optical spectroscopy up to $z\sim 3$ in the COSMOS field, finding no strong evidence for a connection between the AGN activity and star-formation processes in their host galaxies. However, Hickox et al. (2014) point out that the weak correlations between observed AGN properties and host SFR or stellar mass could be explained by the much shorter $\sim$ Myr timescales of AGN activity compared to the $\sim 100$ Myr timescales of star formation in galaxies. These weak trends (see also Lutz et al. 2010; Bonfield et al. 2011) could also partly be due to large scatter in the stellar mass and SFR estimates (see §6.2).

5.1 Selection Bias

Measurements of the $M_{\rm{BH}}-M_{\ast}$ relation are strongly influenced by selection biases at low and high redshift (Lauer et al., 2007; Shen & Kelly, 2010; Shankar et al., 2019) which we discuss below. The probability of a source being included in our $M_{\rm{BH}}-M_{\ast}$ relation sample is given by the probability that (a) the source is variable and detected, (b) the source has a spectrum and the broad-line is detected in the spectrum, (c) the probability that the stellar mass measurement is reliable. Our spectroscopic sample contains spectra from a variety of survey and quasar targeting programs with different targeting criteria, flux limits, and spectral sensitivities. This makes it difficult quantify the effect of the selection biases on our result. Lauer et al. (2007) and Shen & Kelly (2010) point out several selection biases that can result in a false evolution of the $M_{\rm{BH}}-M_{\ast}$ relation. Inevitably, we have selected a sample of sources based on AGN activity, which biases our sample toward more luminous and larger black hole masses as redshift increases. Given the spectral sensitivity limit of $\sim 10^{8}M_{\odot}$ and lack of obvious dependence on the stellar mass reliability with stellar mass (Figure 9), the resulting $M_{\rm{BH}}-M_{\ast}$ relation seems unlikely to be dominated by selection biases (b) and (c). Given other authors have found a similar relation for X-ray and spectroscopically-selected AGNs at similar redshifts for higher AGN luminosities (Merloni et al., 2010; Zhuang et al., 2023), it seems unlikely that the variable population is probing a significantly different parameter space of the Type 1 AGN population. A model of the detectable AGN population at these redshifts may be able to quantify effects of the selection biases on our particular sample (e.g., Pacucci et al. 2023), but requires knowing the true intrinsic scatter of the $M_{\rm{BH}}-M_{\ast}$ relation. Modeling these selection biases is well beyond the scope of this work.

6 Discussion

Figure 12: Spectral fitting results for Kimura et al. (2020) ID 290 ($M_{\ast}\sim 10^{10}M_{\odot}$). The DESI spectral data (black) is plotted with the model fitting components: polynomial+power-law continuum (orange), total line model flux (blue), broad line model flux (red), narrow line model flux (green) fitted host+QSO template (pink) and fitted host template (magenta). In this case, the PCA-derived host flux is small. The resulting broad-line black hole mass is $M_{\rm{BH}}\sim 10^{8}M_{\odot}$ at $z=0.7300$.

6.1 IMBHs in low stellar mass galaxies at similar redshifts

Intriguingly, our results in combination with estimates at higher redshift suggest redshift evolution in the $M_{\rm{BH}}-M_{\ast}$ relation. For example, the recently proposed model of Pacucci & Loeb (2024), motivated by higher redshift over-massive black holes seen in JWST data (e.g., Harikane et al. 2023; Maiolino et al. 2023; Kocevski et al. 2023; Übler et al. 2023; Kokorev et al. 2023; Bogdán et al. 2023; Natarajan et al. 2023; Goulding et al. 2023), predict that at lower redshifts $z\sim 0.5-3$ black holes are expected to be $3-10$ times more massive than the local AGN relation, similar to what is seen in our sample. In the (Pacucci & Loeb, 2024) model, it is argued the $M_{\rm{BH}}/M_{\ast}$ ratio sets the average star formation efficiency in the galaxy, until the galaxy is no longer able to efficiently form stars, bringing it into agreement with the local $M_{\rm{BH}}-M_{\ast}$ relation. Earlier models of BH-galaxy coevolution proposed by Wyithe & Loeb (2003); Caplar et al. (2018) predict similar levels of redshift evolution in the $M_{\rm{BH}}-M_{\ast}$ relation attributed to self-regulated feedback before quenching at a critical value of $M_{\rm{BH}}/M_{\ast}$. While all these models attempt to modulate the SFRs to account for stellar assembly in BH host galaxies, they do not address the key issue of the relationship between the BH accretion rate and the SFR, which likely holds the key to co-evolution. The detailed dependence of this ratio $\dot{M}_{\rm{acc}}/{\dot{M}_{\ast}}$ on feedback; the environment and gas content of galaxies remains to be understood.

In order to build volume and given that the density and accretion rates of black holes peak near $z\sim 2$, several studies have identified IMBH candidates in dwarf galaxies beyond $z=0$. The higher than expected black hole masses than the local AGN relation, whether due to selection effects or intrinsic evolution in the $M_{\rm{BH}}-M_{\ast}$ relation, unfortunately casts doubt on the IMBH nature of previously-identified $L_{2-10\ \rm{{keV}}}\sim 10^{43-44}$ erg s${}^{-1}$ X-ray AGNs at similar redshifts (Mezcua et al., 2018, 2019; Zou et al., 2023). Simply extrapolating the local inactive $M_{\rm{BH}}-M_{\ast}$ relation from Figure 10, we expect an AGN with stellar mass of $M_{\ast}\sim 10^{10}M_{\odot}$ to have a black hole mass of $M_{\rm{BH}}\approx 10^{7.2}M_{\odot}$ rather than $M_{\rm{BH}}\approx 10^{6.4}M_{\odot}$ at similar redshifts.

We have visually inspected the spectra for dwarf galaxies in our sample with $M_{\ast}<10^{10}M_{\odot}$ with reliable stellar mass estimates from SED fitting. One source, Kimura et al. (2020) ID 290 ($M_{\ast}=10^{9.6\pm 0.9}M_{\odot}$), has a sufficiently-detected broad Mg II line with a black hole mass estimate of $M_{\rm{BH}}\sim 10^{7.89\pm 0.05}M_{\odot}$ (statistical fitting error) at $z=0.73$ (Figure 12). Five do not have a spectrum, and the remainder are either too host-dominated or noisy to obtain a broad-line black hole mass estimate. This highlights the need for higher $S/N$ spectroscopy for detecting broad lines in these relatively more host-dominated and fainter variability-selected AGNs in lower mass galaxies.

6.2 Biases in AGN Stellar Mass Estimation and the $M_{\rm{BH}}-M_{\ast}$ Relation

Figure 13: Comparison of stellar masses estimated by Suh et al. (2020) (left) and Bongiorno et al. (2012) (right) and our stellar masses when our stellar mass estimates are considered reliable (blue symbols) or unreliable (black symbols). The solid gray line is the line of $y=x$, and the dotted line falls 1 dex below the line of $y=x$.

Suh et al. (2020) found a sample of $z\sim 0.5-2.5$ X-ray selected AGN to be consistent with the local AGN relation of Reines & Volonteri (2015), contrary to our findings and the findings of several other authors (Merloni et al., 2010; Mezcua et al., 2023; Li et al., 2021; Zhang et al., 2023; Zhuang et al., 2023; Stone et al., 2023; Tanaka et al., 2024). Zhuang et al. (2023) have used black hole masses taken from Suh et al. (2020), and their black hole masses are consistent with ours given the systematic uncertainties in virial mass estimates and the spectral calibration. In order to scrutinize the possible origin of this discrepancy, we compare our stellar masses to those from Suh et al. (2020) in the left panel of Figure 13. We find that the Suh et al. (2020) stellar masses are over-estimated on average compared to our stellar masses and those estimated from COSMOS-Webb AGN-host decomposition by Zhuang et al. (2023). In some cases, the Suh et al. (2020) stellar masses are too large by up to an order of magnitude. The sources with unreliable stellar masses due to being swamped by AGN emission are more likely to be over-estimated by Suh et al. (2017, 2019, 2020). We attribute this discrepancy to inadequate separation of star-formation and AGN-dominated SEDs. There is also a slight tendency for even the reliable stellar masses to be over-estimated. It is possible this originates from their requirement that the UV part of the SED be AGN dominated, which can under-count the contribution from young stars even when the UV AGN fraction is low. The on-average over-estimated stellar masses by Suh et al. (2020) would explain why their sources appear to fall on the local AGN $M_{\rm{BH}}-M_{\ast}$ relation, in contrast with our results. We conclude that the Suh et al. (2020) stellar masses are significantly over-estimated on average. This highlights the extreme caution that must be taken when estimating stellar masses for AGNs without decomposing the AGN+host emission.

Bongiorno et al. (2012) studied the host galaxy properties if both Type 1 and 2 AGNs in the COSMOS field selected from optical spectroscopic and X-ray selection. They derive SED fitting-based stellar masses and SFRs using a similar approach to differentiate between AGN and star-formation dominated SEDs in the NIR. We compare our stellar masses to those from Bongiorno et al. (2012) in the right panel of Figure 13. We only consider sources with constrained stellar mass measurements according to Bongiorno et al. (2012). First, it is evident that a large fraction of the Bongiorno et al. (2012) stellar mass measurements are considered unreliable according to our analysis, having a low ratio between the total model fit for our AGN$+$stellar emission SED model and the AGN-dominated model. Second, it appears that the Bongiorno et al. (2012) stellar masses are systematically lower than our stellar masses below our $M_{\ast}\sim 10^{10.8}M_{\odot}$ for those stellar masses that we consider to be reliable. To reassure that this is not an issue with our stellar masses, we confirm that a substantial fraction (17/96 with $\Delta M_{\ast}>0.5$ dex) of the Bongiorno et al. (2012) stellar masses are also under-estimated compared to those from COSMOS-Webb AGN-host decomposition (Zhuang et al., 2023). Bongiorno et al. (2012) have used the AGN template from Richards et al. (2006), while the cigale AGN models are more flexible while retaining energy balance. Bongiorno et al. (2012) note that their SED fitting results tend to over-estimate the AGN component, which would in-turn under-estimate the stellar emission component. This could likely be the source of under-estimated stellar masses.

6.3 Forecasts for Rubin Observatory

The expected single-epoch limiting magnitude of $r\sim 24.5$ for LSST is more shallow than the HSC-SSP limiting magnitude of $r\sim 25.5$. The single-epoch liming magnitude per exposure depends mostly the exposure time, sky brightness, seeing, and airmass. The single-epoch liming magnitude scales logarithmically with the effective “visit” exposure time as $\propto 1.25\ \log_{10}(t_{\rm{vis}})$ (Ivezić et al., 2019). For example, co-adding an additional 20 visits in the deep drilling fields (in e.g., groups of a few nights) would increase the single-epoch imaging depth by about 1.6 magnitudes, potentially matching or out-performing HSC-SSP in depth. We recommend taking this approach early-on in the survey, which will significantly expand the discovery space to identify faint variables in early Rubin data. After several years of operation, these light curves may be sufficiently long to allow for estimation of the BH mass from the variability timescale (Burke et al., 2021). This would potentially circumvent the need for spectroscopic follow-up to estimate broad-line BH masses.

7 Conclusions

Using the sample of variable AGNs selected from the HSC-SSP COSMOS field (Kimura et al., 2020), we have obtained improved photometric redshifts and multi-wavelength photometry from the COSMOS2020 catalog and used SED fitting to estimate their stellar masses. We have devised an approach for determining the reliability of the stellar mass estimates using the stellar emission strength at 1.2 $\mu$m, where the AGN emission is expected to be at a minimum. After collecting a database of publicly-available spectra from the literature, we measured their virial black hole masses from the detected broad emission lines when $S/N$ permitted. We compared our stellar masses to those from those estimated from AGN-host decomposition using HST and JWST imaging (Zhuang et al., 2023) and measured the $M_{\rm{BH}}-M_{\ast}$ relation for our variable AGNs at $0.5\lesssim z\lesssim 3$. Our results concur with previous findings using other AGN samples of more massive black holes at a given stellar mass than the local AGN relation would suggest. These results suggest that AGNs selected from optical variability are not vastly different from samples of AGNs selected from broad lines at similar redshifts at fixed luminosity.

Using these results as a proxy for LSST Rubin capability, we have demonstrated that black holes with $M_{\rm{BH}}\sim 10^{8}M_{\odot}$ are detectable out to at least $z\sim 4$ in $M_{*}\sim 10^{11}M_{\odot}$ host galaxies using optical variability. Future work combining Rubin-selected AGNs and Roman host galaxy imaging will vastly increase the sample size of sources with reliable stellar masses estimate from AGN$+$host decomposition. Follow-up analysis of this sample will apply Scarlet source deblending (Melchior et al., 2018) to high-resolution time-resolved HSC COSMOS imaging in order to extract the SEDs of the variable AGN and their host galaxies, as well as the host galaxy morphologies, providing further details on host galaxy properties for AGN-dominated cases (Ward et al., in prep). High resolution imaging from JWST could be incorporated into this analysis framework to further improve the AGN and host decomposition, especially at high redshifts.

We are currently performing a similar analysis with variable AGNs from the Dark Energy Survey deep fields (Burke et al., 2022) and an ongoing repeat imaging survey of DES and Rubin deep drilling fields with the Dark Energy Camera (Zhuang et al. in prep). This catalog is not as deep as HSC-SSP, with a single-epoch photometric precision of $g\sim 24.5$, but is a larger area of $\sim 4.6$ deg${}^{2}$, which will allow us to fill in the space at $z<1.5$. This regime is interesting because it corresponds to the epoch when the globally averaged SFR has declined revealing that the overall gas supply available in galaxies for both star formation and accretion have dwindled and the merger rate of galaxies - additional mechanism for injecting gas into galactic nuclei - has also declined. Besides, this epoch lies conveniently in between the local and intermediate-redshift regimes.

We have not investigated narrow emission line ratio diagnostics in this paper. Given the varying wavelength coverage and redshifts of the spectra, the majority of the spectra do not cover both the H$\beta$ and H$\alpha$ spectral complexes. The majority of the sources have AGN features in their spectra (either a broad line detection or Ne V emission line), but a significant fraction are either too noisy or host-dominated to detect any obvious AGN features. Some of the latter could be false positives (i.e., non-AGN galaxy or transient interlopers). This work highlights the challenges of obtaining sufficient spectroscopy to investigate low luminosity AGNs at the redshifts that LSST Rubin will unveil. However, we are optimistic on the capability of optical variability to identify a relatively lower luminosity samples of Type 1 AGNs at low and intermediate redshifts.

Appendix A DATA AVAILABILITY

The HSC PDR3 redshift catalog is available at:
https://hsc-release.mtk.nao.ac.jp/doc/index.php/catalog-of-spectroscopic-redshifts__pdr3/.

The COSMOS2020 catalogs are available at: https://cosmos2020.calet.org.

The publicly-available spectra were collected from the following web pages:

1. zCOSMOS, VUDS, DEIMOS, Magellan, 3D-HST: https://irsa.ipac.caltech.edu/data/COSMOS/spectra/

2. 3D-HST: https://archive.stsci.edu/prepds/3d-hst/

3. C3R2: https://koa.ipac.caltech.edu/Datasets/C3R2/

4. FMOS-COSMOS: http://member.ipmu.jp/fmos-cosmos/FC_spectra.html

5. LEGA-C: https://www.eso.org/sci/publications/announcements/sciann17120.html

6. DESI: https://data.desi.lbl.gov/public/edr/

7. FMOS: https://member.ipmu.jp/fmos-cosmos/FMOS-COSMOS.html

8. LEGA-C: https://users.ugent.be/~avdrwel/research.html

Appendix B STELLAR MASS RECOVERY TESTS

Figure 14: Recoverability of stellar mass in AGNs using mock SEDs with varying multiwavelength coverage generate by cigale. The difference between the log true stellar mass and log recovered stellar mass (in dex) is shown versus the ratio of the AGN continuum emission of an AGN-dominated SED and the total continuum emission of an mixed star-formation $+$ AGN SED calculated near $1.2$ $\mu$m. Points are colored by their true AGN fraction calculated between $0.5-1$ $\mu$m. Sources with a very high AGN fraction ($f_{\rm AGN}\gtrsim 0.8$) are generally well-fitted by an AGN-dominated model which makes inferring the stellar mass difficult. Sources with true $f_{\rm AGN}\sim 1$ are very prone to catastrophically-underestimated stellar masses. In this case, the AGN emission swamps the star-formation emission, and the stellar mass is not well constrained. The inclusion of additional wavelength coverage helps to distinguish between the two different models, constraining the AGN emission, and leading to a more robust stellar mass estimate.

Degeneracies between AGN and stellar emission can lead to highly unreliable stellar mass estimates from SED fitting, depending on the strength of the AGN emission compared to the underlying stellar emission (i.e., the AGN fraction $f_{\rm AGN}$). Although cigale can infer the AGN fraction, it has been shown that quantity cannot be reliably inferred except in AGN-dominated sources ($f_{\rm AGN}\gtrsim 0.8$). For example, a young stellar population can be fitted to a UV/optical AGN continuum emission, leading to catastrophically-underestimated stellar masses. In order to overcome this degeneracy, we separate our sample into sources with “reliable” and “unreliable” stellar mass estimates using a model comparison test. Our approach is summarized in §3.1, and we describe it in detail below:

1. Fit the observed SED with a totally AGN-dominated model by setting $f_{\rm{AGN}}=0.9999$, computed over the wavelength range $0.5-1$ $\mu$m.

2. Fit the observed SED with an AGN$+$SF model with $f_{\rm{AGN}}$ as a free parameter varying between $0-1$, computed at observed-frame $0.1-0.3$ $\mu$m. This model assures that the UV component of the SED is AGN-dominated, but allows the redder optical and near-infrared part of the SED to be AGN or SF dominated.

3. Compute the ratio of the the best-fit total continuum emission from the AGN$+$SF model over the best-fit AGN continuum emission from the AGN-dominated model at rest-frame 1.2 $\mu$m (where the AGN emission is at a minimum).

4. Recovered stellar masses from the AGN$+$SF model are considered reliable if its best-fit reduced $\chi^{2}<5$ and the excess stellar emission from (3) is greater than 1.2, as justified using mock tests.

In order to test this procedure, we use the “savefluxes” mode of cigale to generate mock photometry at varying AGN fractions. We allow the stellar and AGN parameters as in Table 2. We then run the mock photometry through the same cigale fitting procedure and stellar-mass reliability tests as our real data. The recovered stellar masses given the true mock photometry is shown in Figure 14. The value of 1.2 is chosen to eliminate $>95\%$ of the catastrophic failures and does not depend strongly on the wavelength coverage (or redshift) of the SED. Although increased wavelength coverage does help increase the number of sources that meet this reliable stellar mass criteria.

The obvious limitation of this test is that it does not include systematic uncertainties in the recovered stellar masses beyond the parameter choices in the mock SEDs we have generated. Nevertheless, our results are broadly consistent with detailed tests of stellar mass estimates from SED fitting of AGNs using independently-calculated photometry from detailed star formation histories (Ciesla et al., 2015). We have also used stellar masses of our AGNs from spatially-decomposed imaging in the literature as an independent check on our stellar masses in Figure 6.

Appendix C RELIABILITY OF SPECTRAL CONTINUUM MEASUREMENTS

Figure 15: Black hole masses estimated from the broad-line luminosity versus ($y$ axis) AGN continuum luminosity ($x$ axis). The black hole masses from the AGN continuum luminosity are susceptible to a few catastrophically under-estimated masses when the AGN continuum is not well constrained.

AGNs with bolometric luminosities $L_{\rm{bol}}\lesssim 10^{45}$ erg s${}^{-1}$ tend to have significant contribution from the host galaxy (Shen et al., 2011; Kimura et al., 2020). Reliably constraining the quasar continuum luminosity is essential to obtaining a virial black hole mass estimate using the prescriptions that use the continuum luminosity as a proxy for the BLR luminosity (Shen et al., 2011). PCA decomposition can constrain the quasar continuum for sources below $L_{\rm{bol}}\sim 10^{45}$ erg s${}^{-1}$ for high quality spectra. However, our spectral quality (calibration and $S/N$) vary considerably depending on the instrument and spectral reduction. In order to test whether the quasar continuum luminosities are well-constrained, we computed the virial black hole masses using both the broad-line luminosity (Equation 1) or continuum luminosity approaches using the equation:

$$\begin{split}\log{\left(\frac{M_{\rm{BH}}}{M_{\odot}}\right)}=a+b\log{\left(\frac{\lambda L_{\lambda}}{10^{44}\ \rm{erg\ s}^{-1}}\right)}+2\log{\left(\frac{\rm{FWHM}_{\rm{br}}}{\rm{km\ s}^{-1}}\right)}\end{split}$$ (C1)

where $\lambda L_{\lambda}$ and FWHM${}_{\rm{br}}$ are the continuum luminosity and broad-line full-width-at-half-maximum (FWHM) with an intrinsic scatter of $\sim 0.4$ dex in BH mass. We adopt the calibrations (Vestergaard & Peterson, 2006) used in Shen et al. (2011):

$$(a,b)=(0.910,0.50),\quad\rm H\beta$$ (C2)

$$(a,b)=(0.740,0.62),\quad\rm Mg\ II$$ (C3)

$$(a,b)=(0.660,0.53),\quad\rm C\ IV.$$ (C4)

If the quasar continuum luminosities are well constrained, we expect the two approaches to yield consistent results within systematic uncertainties between the two prescriptions. We found that the black hole masses from the AGN continuum luminosity are susceptible to a few catastrophically under-estimated masses when the AGN continuum is not well constrained, as shown in Figure 15. For consistency and simplicity, we adopt black hole masses estimated from the broad line luminosity from Equation 1.

Appendix D SN 1000+0216

One source, SN 1000+0216 (ID=451), is reported as a superluminous supernova (SLSN) by Cooke et al. (2012) based on variability from years 2005$-$2008. The $z=3.8993$ source does not have obvious SN or AGN features in the spectrum presented in Cooke et al. (2012). It is possible the spectrum was taken in a host-dominated state, i.e., when the AGN luminosity happened to be at a minimum. Since the source is still variable in 2014$-$2017 HSC-SSP light curves, the quasar interpretation may be more likely. A detailed study of this source and further scrutiny of its light curve is required to reach a firm conclusion.