H0LiCOW VII: cosmic evolution of the correlation between black hole mass and host galaxy luminosity

Abstract

Strongly lensed active galactic nuclei (AGN) provide a unique opportunity to make progress in the study of the evolution of the correlation between the mass of supermassive black holes (⁠|$\mathcal {M}_{\rm BH}$|⁠) and their host galaxy luminosity (L_host). We demonstrate the power of lensing by analysing two systems for which state-of-the-art lens modelling techniques have been applied to deep Hubble Space Telescope imaging data. We use (i) the reconstructed images to infer the total and bulge luminosity of the host and (ii) published broad-line spectroscopy to estimate |$\mathcal {M}_{\rm BH}$| using the so-called virial method. We then enlarge our sample with new calibration of previously published measurements to study the evolution of the correlation out to z ∼ 4.5. Consistent with previous work, we find that without taking into account passive luminosity evolution, the data points lie on the local relation. Once the passive luminosity evolution is taken into account, we find that black holes in the more distant Universe reside in less luminous galaxies than today. Fitting this offset as |$\mathcal {M}_{\rm BH}$|/L_host ∝ (1 + z)^γ, and taking into account selection effects, we obtain γ = 0.6 ± 0.1 and 0.8 ± 0.1 for the case of |$\mathcal {M}_{\rm BH}$|–L_bulge and |$\mathcal {M}_{\rm BH}$|–L_total, respectively. To test for systematic uncertainties and selection effects we also consider a reduced sample that is homogeneous in data quality. We find consistent results but with considerably larger uncertainty due to the more limited sample size and redshift coverage (γ = 0.7 ± 0.4 and 0.2 ± 0.5 for |$\mathcal {M}_{\rm BH}$|–L_bulge and |$\mathcal {M}_{\rm BH}$|–L_total, respectively), highlighting the need to gather more high-quality data for high-redshift lensed quasar hosts. Our result is consistent with a scenario where the growth of the black hole predates that of the host galaxy.

black hole physics, galaxies: active, galaxies: evolution

1 INTRODUCTION

It is commonly accepted that almost all the galaxies have a supermassive black hole (BH) in their centre, whose mass (⁠|$\mathcal {M}_{\rm BH}$|⁠) is known to be correlated with the host properties. The tight correlations are usually, but not uniquely, explained as the results of their co-evolution (e.g. Magorrian et al. 1998; Ferrarese & Merritt 2000; Gebhardt et al. 2001; Marconi & Hunt 2003; Häring & Rix 2004; Gültekin et al. 2009; Graham et al. 2011; Beifiori et al. 2012; Kormendy & Ho 2013; Park et al. 2015, hereafter P15) (see, however, Peng 2007; Jahnke & Maccio 2011, for a different view). A powerful way to explore the origin of this physical coupling and understand the role of active galactic nuclei (AGN) feedback in galaxy formation is to measure the correlations directly at a high redshift and determine how and when they emerged and evolved over cosmic time [e.g. Treu, Malkan & Blandford 2004; Salviander et al. 2006; Woo et al. 2006; Jahnke et al. 2009; Schramm & Silverman 2013, hereafter SS13; DeGraf et al. 2015).

The most common technique used to estimate |$\mathcal {M}_{\rm BH}$| beyond the local Universe (z > 0.1) is the so-called virial method, based on the properties of broad emission lines in type 1 AGN (Shen 2013; Peterson 2014). However, the bright source associated with the AGN makes the study of its host galaxy very difficult. Strong gravitational lensing (see e.g. Courbin, Saha & Schechter 2002; Schneider, Kochanek & Wambsganss 2006; Treu 2010; Treu & Ellis 2015, for reviews) stretches the host galaxy out from the wings of the bright point source as point spread function (PSF), providing a unique opportunity to infer its magnitude robustly (Peng et al. 2006, hereafter P06). However, in order to measure host luminosity (L_host) and construct the |$\mathcal {M}_{\rm BH}$|–L_host correlation from strongly lensed AGN, it is necessary to ensure that any systematic uncertainties associated with the gravitational lens model should be controlled to the desired level of accuracy.

Recently, Ding et al. (2017) studied the fidelity of the measurement of lensed AGN host brightness through a set of extensive and realistic simulations of Hubble Space Telescope (HST) observation and lens modelling. First, the mock images of the lensed AGNs in our sample (see Ding et al. 2017, table 1) were generated as realistically as possible. Secondly, the simulated AGN host galaxy images were reconstructed with the state-of-the-art lens modelling tool (glee ¹). Thirdly, by fitting the host magnitude with the software galfit (Peng et al. 2002) and comparing the inference to the input value, Ding et al. (2017) found that L_host can be recovered with better accuracy and precision than the uncertainty on single-epoch |$\mathcal {M}_{\rm BH}$| estimates (∼0.5 dex) for hosts as faint as 2–4 mag dimmer than the AGN itself.

In this paper, we apply our advanced techniques to two strongly lensed systems (i.e. HE0435−1223 and RXJ1131−1231), with excellent imaging data. The host galaxy luminosity is inferred from the lens detailed model developed as part of the H0LiCOW collaboration² with the goal of measuring cosmological parameters from gravitational time delays (Bonvin et al. 2017; Suyu et al. 2017). |$\mathcal {M}_{\rm BH}$| is inferred by applying a set of self-consistent calibrations of the virial method to the broad emission-line properties measured by Sluse et al. (2012). In addition, we combine our new measurements with a large sample of AGNs taken from the literature and consistently recalibrated, and study the evolution of the |$\mathcal {M}_{\rm BH}$|–L_host relation for 146 objects in the redshift range 0 < z < 4.5. It is still unclear whether the bulge or the total luminosity provides the tightest correlation with |$\mathcal {M}_{\rm BH}$| (Jahnke et al. 2009; Bennert et al. 2011b, hereafter B11; P15). Thus, we consider both of them in this study.

This paper is organized as follows. We briefly describe the sample selection in Section 2. The host galaxy surface photometry and the |$\mathcal {M}_{\rm BH}$| are inferred in Sections 3 and 4, respectively. In Section 5, we present our main result. Discussion and conclusion are presented in Sections 6 and 7. Throughout this paper, we adopt a standard concordance cosmology H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.30 and |$\Omega {_\Lambda } = 0.70$|⁠. Magnitudes are given in the AB system.

2 SAMPLE SELECTION

First, we analyse the two quadruply imaged AGN HE0435−1223 and RXJ1131−1231 (hereafter HE0435 and RXJ1131) with source redshifts at 1.693 and 0.654, respectively. Detailed information for these two systems is given by Suyu et al. (2017). Accurate lens models have been derived in an effort to measure cosmological parameters from gravitational time delays as described by Wong et al. (2017) and Suyu et al. (2013). These models provide the reconstructions of AGN hosts, from which in turn we estimate L_host.

Secondly, we combine and compare our new measurements with those by P06. P06 explored the |$\mathcal {M}_{\rm BH}$|–L_host relation based on 20 non-lensed AGNs and 31 gravitationally lensed AGNs (including HE0435 and RXJ1131). P06 is so far the only paper in which the |$\mathcal {M}_{\rm BH}$|–L_host relation has been comprehensively investigated using lensed AGNs observed with HST. We note that for the two systems in common, the HST images used in our work are much deeper than those used by P06, and the lens models are much more detailed. Also, P06 was based on NIC2 images, as opposed to the much more powerful and more modern cameras used in our work. Therefore, our measurements supersede those by P06 for these two systems. Furthermore, we exclude MG 2016+112 because it is a type II AGN (Koopmans et al. 2002) and the black hole mass using the virial method cannot be considered reliable. We also exclude the lens system B2045+265 used by P06 because of the incorrect redshift identification of the AGN spectrum by Fassnacht et al. (1999) (Nierenberg et al., in preparation).

Thirdly, we combine our new measurements with samples of non-lensed AGN that have been measured by members of our team using the same techniques as those applied here. The samples consist of 52 intermediate-redshift AGNs (0.36 < z < 0.57) summarized by P15, 27 distant AGNs (0.5 < z < 1.9) measured by B11 and SS13, and 19 local AGN measurements (Peterson et al. 2004; Bennert et al. 2010). It is worth noticing that they are so far the largest HST imaging samples which are carefully selected as moderate-luminosity AGN, for which the contrast between nucleus and host galaxies is much more favourable for the inference of L_host than for high-luminosity lensed quasars. Thus, their host luminosities are measured with high accuracy even without lensing.

Overall, our sample consists of two new lensed systems and active galaxies from the literature, including elliptical and spiral hosts with redshift up to 4.5. This total sample of 146 objects is the largest compilation of AGNs from HST which are cross-calibrated to study the |$\mathcal {M}_{\rm BH}$|–L_host relation. The objects and their basic properties are listed in Tables 1 and 2.

Table 1.

Properties of AGNs in the distant sample.

Object	Line(s) used	Redshift	log |$\mathcal {M}_{\rm BH}$|	log Lhost,R
			(M⊙)	(L⊙)
			±0.4 dex	±0.2 dex
HE 0435	MgII	1.693	8.61	10.96
RXJ 1131bulge	MgII/Hβ	0.654	8.26/8.41	10.58
RXJ 1131disc	MgII/Hβ	0.654	8.26/8.41	11.12
Lensed AGNs from P06
RXJ 1131	Hβ	0.66	7.90	11.02
SDSS 1226−0006	MgII	1.12	8.41	10.74
FBQ 0951+2635	MgII	1.24	8.57	10.25
CTQ 414	MgII/CIV	1.29	7.78/8.19	10.87
B 0712+472	MgII	1.34	7.44	10.90
SBS 0909+532	MgII	1.38	9.13	10.54
Q 0957+561	MgII/CIV	1.41	9.06/8.97	11.79
FBQ 1633+3134	MgII/CIV	1.52	8.84/8.91	11.08
SDSS 0924+0219	MgII	1.54	7.61	11.09
B 1030+071	MgII	1.54	8.13	11.06
SDSS 1335+0118	MgII	1.57	8.77	10.97
B 1600+434	MgII	1.59	7.56	10.98
HE 0047−1756	CIV	1.66	8.83	11.24
HE 0435	CIV	1.69	8.36	11.12
PG 1115+080	MgII/CIV	1.72	8.67/8.63	11.08
SBS 1520+530	CIV	1.86	8.60	10.82
HE 2149−2745	CIV	2.03	9.48	11.47
HE 1104−1805	CIV	2.32	9.03	11.48
Q 1017−207	CIV	2.55	8.88	11.69
H 1413+117	CIV	2.55	8.08	11.48
MG 0414+0534	Hβ	2.64	9.07	11.41
J 1004+1229	CIV	2.65	8.97	11.70
Q 0142−110	CIV	2.72	9.01	11.37
LBQS 1009−0252	MgII/CIV	2.74	8.51/8.70	11.48
RXJ 0911+0551	CIV	2.80	8.57	10.93
PMNJ 1632−0033	CIV	3.42	8.25	11.35
B 1422+231	MgII/CIV	3.62	8.93/9.34	11.60
BRI 0952−0115	CIV	4.5	8.80	11.95
Non-lensed AGNs from P06
PKS 0440−00	MgII	0.844	8.09	11.29
MGC 2214+3550A	MgII	0.88	8.76	11.13
MGC 2214+3550B	MgII	0.88	8.26	10.59
3C 422	MgII	0.942	9.04	11.64
PKS 0938+18	MgII	0.943	8.53	11.16
SGP5:46	MgII	0.955	8.03	10.92
LBQS 1009−0252c	CIV	1.63	8.88	11.34
RXJ 0921+4528A	MgII/CIV	1.65	9.09/8.71	11.12
RXJ 0921+4528B	MgII/CIV	1.65	8.82/8.4	11.38
SGP4:39	CIV	1.716	8.07	10.49
MZZ 11408	CIV	1.735	8.00	10.82
SGP2:36	CIV	1.756	8.87	11.26
MZZ 1558	CIV	1.829	8.80	10.94
SGP2:25	CIV	1.868	8.45	11.27
MZZ 4935	CIV	1.876	8.02	10.43
SGP3:39	CIV	1.964	8.69	11.38
SGP2:11	CIV	1.976	8.69	11.04
4C 45.51	MgII/CIV	1.992	8.60/8.48	12.18
MZZ 9592	CIV	2.71	8.47	11.49
MZZ 9744	CIV	2.735	8.52	11.09
Non-lensed AGNs from P15				(Lbulge/ Ltotal)
S09a	Hβ	0.3545	7.99	11.03/11.05
S10	Hβ	0.3513	8.44	10.60/11.03
S12	Hβ	0.3583	8.78	10.11/11.02
S21	Hβ	0.3546	8.93	10.38/11.30
S16	Hβ	0.3702	8.02	9.89/10.53
S23	Hβ	0.3513	8.82	10.27/11.04
S24	Hβ	0.3619	8.20	10.95/11.14
S26	Hβ	0.3692	8.00	10.47/10.76
S27	Hβ	0.3669	7.85	11.00/11.07
S01	Hβ	0.3594	8.15	10.41/10.87
S02	Hβ	0.3544	8.02	10.43/10.68
S03	Hβ	0.3584	8.17	9.64/11.00
S04	Hβ	0.3579	8.11	10.51/10.98
S05	Hβ	0.3535	8.73	9.75/10.94
S06	Hβ	0.3688	7.72	9.53/10.95
S07	Hβ	0.352	8.53	10.28/11.02
S08	Hβ	0.3586	7.74	10.02/10.89
S11	Hβ	0.3559	7.76	10.50/10.90
SS1	Hβ	0.3566	7.73	10.13/10.94
SS2	Hβ	0.3671	7.56	10.76/10.76
SS5	Hβ	0.3735	7.98	10.02/10.66
Non-lensed AGNs from P15				log Lbulge/log Ltotal
S31	Hβ	0.3506	8.19	10.76/11.00
SS6	Hβ	0.3588	7.39	9.67/10.47
SS7	Hβ	0.3613	7.67	10.08/10.74
SS8	Hβ	0.3655	7.85	9.96/10.92
SS9	Hβ	0.3702	7.95	10.82/10.82
SS10	Hβ	0.3658	8.10	10.61/10.81
SS11	Hβ	0.3731	7.83	10.04/10.83
SS12	Hβ	0.3629	8.15	10.54/10.66
SS13	Hβ	0.3743	7.69	10.60/10.60
S28	Hβ	0.3678	8.12	10.70/10.91
SS14	Hβ	0.3706	7.45	10.48/10.48
S29	Hβ	0.3574	7.95	10.03/10.69
SS18	Hβ	0.3585	7.51	9.81/10.62
W11	Hβ	0.565	7.95	10.73/10.73
W22	Hβ	0.5652	8.68	11.27/11.27
W12	Hβ	0.5623	8.94	10.49/10.98
W20	Hβ	0.5761	8.60	11.03/11.03
W16	Hβ	0.578	7.86	10.81/10.81
W8	Hβ	0.5712	8.74	11.01/11.01
W3	Hβ	0.576	8.76	10.28/10.88
SS15	Hβ	0.3593	7.44	10.36/10.36
W1	Hβ	0.5736	8.84	10.71/10.95
W4	Hβ	0.5766	8.28	11.01/11.01
W5	Hβ	0.5767	8.29	10.96/10.96
SS3	Hβ	0.3566	7.51	10.10/10.78
SS4	Hβ	0.3629	7.85	10.81/10.81
W17	Hβ	0.5617	8.31	10.12/10.80
W2	Hβ	0.572	9.07	10.91/10.91
W10	Hβ	0.5711	7.94	10.51/10.77
W14	Hβ	0.5617	8.68	10.97/10.97
W9	Hβ	0.5654	8.70	10.89/10.96
J033252−275119	MgII	1.227	8.87	10.53/11.28
J033243−274914	MgII	1.900	9.17	11.73/11.73
J033239−274601	MgII	1.220	8.24	11.22/11.22
J033226−274035	MgII	1.031	7.85	10.37/11.46
J033225−274218	MgII	1.617	8.08	11.48/11.48
J033210−274414	MgII	1.615	8.30	11.52/11.52
J033200−274319	MgII	1.037	7.75	10.47/10.47
J033229−274529	MgII	1.218	8.37	11.39/11.39
J123553+621037	MgII	1.371	8.27	10.73/11.66
J123618+621115	MgII	1.021	8.35	10.11/11.37
J123618+621115	MgII	1.450	8.77	11.45/11.45
158b	MgII	0.717	7.28	10.71/10.89
170	MgII	1.065	7.07	10.16/10.46
271	MgII	0.960	7.43	10.46/11.36
273	MgII	0.970	8.23	10.19/10.45
305	MgII	0.544	8.61	11.02/11.14
333	MgII	1.044	7.90	10.45/10.91
339	MgII	0.675	7.95	10.83/11.00
348	MgII	0.569	8.11	10.64/11.17
379	MgII	0.737	9.14	11.03/11.70
413	MgII	0.664	7.05	10.19/10.74
417	MgII	0.837	8.37	10.18/10.87
465	MgII	0.740	8.02	10.67/11.47
516	MgII	0.733	7.93	11.02/11.38
540	MgII	0.622	7.61	10.84/11.30
597	MgII	1.034	8.12	10.80/10.91
712	MgII	0.841	8.65	11.22/11.59

Note: – Column 1: object ID. Column 2: Emission line used to estimate |$\mathcal {M}_{\rm BH}$|⁠. Column 3: redshift as listed in the literature. Column 4: |$\mathcal {M}_{\rm BH}$| calibrated from equation (2) using the corresponding lines. Column 5: Inferred rest-frame R-band luminosity not corrected for evolution. Note that all the host galaxies in P06 are assumed to be pure ellipticals.

^aID taken from P15.

^bID taken from SS13.

Table 2.

Properties of local AGNs.

Object	redshift	log |$\mathcal {M}_{\rm BH}$|	log Lbulge, R/log Ltotal,R
		(M⊙)	(L⊙)
			±0.2 dex
3C120	0.03301	7.71 ± 0.21	10.51/10.51
3C390.3	0.0561	8.43 ± 0.10	10.55/10.55
Ark120	0.03271	8.15 ± 0.06	10.23/10.69
Mrk79	0.02219	7.69 ± 0.12	9.62/10.01
Mrk110	0.03529	7.37 ± 0.11	9.42/9.93
Mrk279	0.03045	7.51 ± 0.11	9.87/10.32
Mrk335	0.02579	7.12 ± 0.11	9.61/10.00
Mrk590	0.02639	7.65 ± 0.07	9.98/10.45
Mrk817	0.03146	7.66 ± 0.07	9.34/10.32
PG0052+251	0.155	8.54 ± 0.09	11.34/11.34
PG0804+761	0.1	8.81 ± 0.05	10.81/10.81
PG0844+349	0.064	7.94 ± 0.18	10.54/10.54
PG1211+143	0.0809	8.13 ± 0.13	10.43/10.43
PG1226+023	0.15834	8.92 ± 0.09	11.64/11.64
PG1229+204	0.06301	7.83 ± 0.21	10.30/10.70
PG1411+442	0.0896	8.62 ± 0.14	10.58/10.58
PG1613+658	0.129	8.42 ± 0.20	11.47/11.47
PG1700+518	0.292	8.86 ± 0.10	11.51/11.51
PG2130+099	0.06298	7.55 ± 0.17	9.82/10.42

Note. Local AGN measurements, taken from Bennert et al. (2010). Following P15, we adopted virial factor as log f = 0.71. Note that Bennert et al. (2010) adopted log f = 0.74.

3 SURFACE PHOTOMETRY

In this section, we describe the measurement of host luminosity. For HE0435 and RXJ1131, we first derived their host magnitude from the reconstructed surface brightness maps in the source plane. Then, we inferred the rest-frame R-band luminosities based on their spectral energy distribution (SED). For the other AGNs, we collected and homogenized their luminosities from the literature.

3.1 Surface photometry of HE0435 and RXJ1131

We used the software galfit to model the reconstructions from Wong et al. (2017) and Suyu et al. (2013). The reconstruction of HE0435 was fitted as the Sérsic profile with n limited between 1 and 4. It has been tested that this prior on n does not bias the inference of magnitude (Ding et al. 2017). In the case of RXJ1131 a clearly visible residual image was present and the resulting parameters were physically acceptable when fitted with an additional profile; we concluded that the host galaxy is composed of a disc and a bulge. In this case, we fixed the reconstruction as two-component Sérsic profiles with n equal to 1 and 4, corresponding to exponential disc profile and de Vaucouleurs (1948) profile, respectively. Although the luminosities of lens systems are corrected from lensing magnification using a lens model, small differences exist between models of different groups. The derived magnification rarely differs by more than 20 per cent. According to detailed simulations presented by Ding et al. (2017), the inferred values of L_host can be recovered with sufficient accuracy and precision to study the |$\mathcal {M}_{\rm BH}$|–L_host relation using our approach. Finally, we derived the rest-frame R-band luminosity using a standard K-correction. These steps are described below in more detail for each system.

3.1.1 HE0435

HE0435 was imaged with HST/WFC3-IR through filter F160W from program HST-GO-12889 (PI: S. H. Suyu). Wong et al. (2017) produced a set of 12 reconstructions for this system, based on different assumptions, in order to estimate the amplitude of systematic errors associated with these choices. In nine out of 12 cases, the source plane resolution was set to 40 × 40 pixels. For the other three cases, a higher resolution of 50 × 50 pixels was adopted. The reconstructions were based on an image plane size of ∼1.9 arcsec².

By fitting each of the 12 reconstructions with single Sérsic profile, we summarized the inference and found the mean value and the scatter of the host magnitudes are m_host = 21.75 ± 0.13; the inferred effective radius and Sérsic index are R_eff = 0.82 ± 0.14 arcsec; n = 3.94 ± 0.14, as shown in Table 3. Furthermore, to test the type of the host galaxy, we fitted the reconstructions as two-component Sérsic profile. However, we obtained unphysical results and no improvements in the fit indicating that the host galaxy of HE0435 is consistent with being a pure elliptical. One example of the reconstruction and its corresponding galfit best fit is shown in Fig. 1 (panel a). We also note that there is a small structure at the lower left of the host. However, its brightness is negligible compared to the host which does not affect the inference of the L_host. Interestingly, this could correspond tidal features in the host galaxies. If true, the mergers could be related to triggered AGN activity. It is beyond the scope of this work to pursue this further, but it would be intriguing to simulate the hosts with merger signature and to see if they can be recovered in the source reconstruction.

Figure 1.

Illustration of the surface photometry study of HE0435, presented with the same stretch for each panel, based on HST/WFC3-IR images through filter F160W.

Table 3.

The inference of HE0435 and RXJ1131.

Object	Magnitude	Reff	Sérsic index (n)
		(arcsec)
HE0435	21.75 ± 0.13	0.82 ± 0.14	3.94 ± 0.14
RXJ1131disc	20.07 ± 0.06	0.84 ± 0.09	Fixed 1
RXJ1131bulge	21.81 ± 0.28	0.20 ± 0.08	Fixed 4

We can verify the accuracy of our result by carrying out simulations as described in our previous paper (Ding et al. 2017), using our inferred parameters as input. The observed and simulated HE0435 images are shown in Fig. 1 (panel b). By repeating the analysis on the simulated image, we recover the input value (input: m_host = 21.75 mag; output 21.88 mag) showing an accuracy much better than our target 1.25 mag (0.5 dex). We note that while in the simulations the PSF is assumed to be perfectly known, for the real data the PSF is inferred from the data using an iterative correction procedure (see Chen et al. 2016; Wong et al. 2017; Suyu et al., in preparation).

Following P06, we made no corrections for dust extinction of the host galaxy because they are likely to be small for a pure elliptical. The observed magnitudes were then transformed to rest-frame R band by applying K-correction with Sbc template spectrum using Coleman, Wu & Weedman (1980) templates. We used the Sbc template because the stellar populations cannot be older than a few Gyr at this redshift and the local elliptical template would be too red. Nevertheless, since the HE0435 is observed through the F160W filter, which roughly corresponds to the rest-frame R band at z ∼ 1.5, the K-correction is only weakly dependent on the assumed SED (see fig. 7 in P06), and does not contribute significantly to the error budget. Finally, the best inferred value of L_host of HE0435 in rest-frame R band is log L_host = 10.96, which is very close to the one inferred in P06 (i.e. log L_host = 11.12).

Lens models based on archival HST Advanced Camera for Surveys (ACS) in the filters F555W and F814W are also available from Wong et al. (2017). Unfortunately, due to the short exposure time, the signal-to-noise ratio of the reconstructed host images in these bands is insufficient to infer the luminosity robustly in these bands and study the colours of the host. Thus they are not considered in this study.

3.1.2 RXJ1131

RXJ1131 is imaged with HST/ACS through filter F814W. A set of seven source resolutions including 50 × 50, 52 × 52, 54 × 54, 56 × 56, 58 × 58, 60 × 60 and 64 × 64 pixels was selected when modelling the host image into source plane (Suyu et al. 2013), with a frame size of ∼2.9 arcsec².³

As noted by Suyu et al. (2013), all the reconstruction of the host shows a compact peak near the centre (see Fig. 2, panel (a), left panel), exhibiting the boundary line between the dominated area of bulge and disc which indicates the host galaxy is a spiral galaxy. Similarly, Claeskens et al. (2006) reconstructed the host of RXJ1131 and found it to be a spiral, disc-dominated galaxy. Thus, we fitted the reconstructions as two-component Sérsic profiles, and the inferred properties of the disc are m_disc = 20.07 ± 0.06 mag, R_eff|$\__{\rm {disc}}=0.84\pm 0.09$| arcsec and the properties of the bulge are m_bulge = 21.81 ± 0.28 mag and R_eff|$\__{\rm bulge}=0.20\pm 0.08$| arcsec, as summarized in Table 3. An example of the reconstruction and the best-fitting image is shown in Fig. 2 (panel a).

Figure 2.

Illustration of the surface photometry study of RXJ1131, presented with the same stretch for each panel, based on HST/ACS images through filter F814W.

In the simulations of Ding et al. (2017), the host of RXJ1131 was assumed to be a single Sérsic profile with the magnitude between 19.0 and 20.5. In this work, we simulate a more realistic two-component profile, with key parameters (i.e. m_host and R_eff) equal to the inferred values. The real and mock RXJ1131 images are shown in Fig. 2 (panel b). We first use a single Sérsic profile to fit the reconstruction, but applying this model is a poor representation with an obvious residual in the central image (i.e. Fig. 2, panel c, left). This result suggests the lens model of RXJ1131 reconstructs the host with a sufficiently high resolution to distinguish a bulge+disc model from a single component. Fitting with two-component Sérsic profile, we find that the residual map is much improved and both components can be reconstructed accurately with our data and analysis techniques: input m_disc = 20.07 mag and m_bulge = 21.80 mag; inferred values are m_disc = 20.37 mag and m_bulge = 22.07 mag.

As for HE0435, we derived the rest-frame R-band magnitude using a standard K-correction. At the redshift of RXJ1131, the conversion to R-band magnitude depends significantly on the adopted SED. Therefore, we determined the K-correction directly from the colour of lensed host arc, based on the multi-band SED fitting available in the archive (GO-9744; PI: C. S. Kochanek). The final estimations are Δmag_disc(R−F814W) ≈ −0.3 and Δmag_bulge(R−F814W) ≈ −0.7. For detail, see Appendix A.

3.2 Surface photometry for the literature samples

In this section we describe our inference of the rest-frame R-band luminosity for the P06 and P15 samples.

P06 used the galfit (for non-lensed source) and lensfit⁴ (for lensed source) softwares to infer the brightness of the AGN hosts, describing the host galaxy as single Sérsic profile. In P06, they reported a single value of luminosity for each object, suggesting that the host galaxies are ellipticals. However, in our analysis, we find that RXJ1131 is a spiral galaxy which suggests the approach in P06 may not be accurate for all the host galaxies. We will return to this issue in Section 6. Their measurements of absolute magnitude are presented by P06 (tables 3 and 4 therein) in rest-frame R-band, Vega system. Thus, we transfer to AB system using m_{AB, R} − m_{Vega, R} = 0.21 (Blanton & Roweis 2007).

Similarly, for the P15 sample, which includes the samples from B11 and SS13, the host galaxy was fitted as an n = 4 profile to model the bulge component; an exponential disc profile was added if deemed necessary. The rest-frame V-band luminosity is derived (see P15, table 4, column 3) by applying the K-correction with an early-type galaxy template spectrum. The same template is taken; we converted their results to rest-frame R band. As the scatter in V − R colours is small, the associated uncertainty is estimated to be 0.16 mag (i.e. 0.06 dex in luminosity). Likewise, the luminosities for 19 local active galaxies are converted to rest-frame R band.

Having obtained the R-band mag, the luminosity is derived by log L_R/L_{R, ⊙} = 0.4(M_{R, ⊙} − M_R), where M_{R, ⊙} = 4.61 (Blanton & Roweis 2007). We summarized the homogenized R-band luminosities in Tables 1 and 2.

4 BLACK HOLE MASS

Assuming that the dynamics of the broad-line region (BLR) is dominated by the gravity of the central supermassive black hole, |$\mathcal {M}_{\rm BH}$| can be derived by applying the so-called virial method, based on the size of the BLR (R_BLR) and the line-of-sight velocity width (ΔV) which can be inferred in turn from continuum luminosity and emission-line width, respectively. Usually, the C_IV(λ1549), Mg_II(λ2798) and Hβ(λ4861) emission-line width and their local continuum luminosities λL_λ(1300 Å), λL_λ(3000 Å) and λL_λ(5100 Å) are used, respectively.

Sluse et al. (2012), P06 and P15 used different lines and different calibrations of the virial method. Thus, we need to cross-calibrate them in order to avoid any systematic bias between the samples.

We first choose the recipe of P15 as the baseline:

\begin{eqnarray} \log \left(\frac{\mathcal {M}_{\rm BH}} {{\rm M}_{\odot }}\right)& = & 7.536 + 0.519 \log \left(\frac{\lambda L_{5100}}{10^{44}\, \rm erg\, s^{-1}}\right) \nonumber \\ &&+\, 2 \log \left(\frac{\sigma _{\rm{H\beta }}}{1000\, \rm km\, s^{-1}}\right). \end{eqnarray}

(1)

Then, we align the self-consistent recipes (including emission lines using Hβ and Mg_II) from McGill et al. (2008) with this baseline by adding a small constant to the intercept (i.e. −0.144). In order to cross-calibrate the C_IV-based estimator, we exploit the nine AGNs in our sample for which both Mg_II and C_IV are available. We take the C_IV recipe from P06 and add a small constant intercept (i.e. −0.331) to match on average the value inferred from Mg_II. Overall, we adopt the following virial formalism:

\begin{eqnarray} \log \left(\frac{\mathcal {M}_{\rm BH}}{\mathrm{M}_{{\odot }}}\right)& = & a+b \log \left(\frac{ \lambda L _{\lambda _{\rm line}}}{10^{44}{\rm erg\, s^{-1}}}\right) \nonumber \\ &&+\, 2 \log \left(\frac{\rm FWHM(line)}{1000\, {\rm km\, s^{-1}}}\right), \end{eqnarray}

(2)

with a {C_IV, Mg_II, Hβ} = {6.322, 6.623, 6.882} , b{C_IV, Mg_II, Hβ} = {0.53, 0.47, 0.518}, and λ_line{C_IV, Mg_II, Hβ} = {1350, 3000, 5100}. Having achieved a consistent cross-calibration, |$\mathcal {M}_{\rm BH}$| is estimated by adopting the emission-line properties measured by Sluse et al. (2012), P06 and P15.

For the 19 local AGNs, rather than using continuum luminosity, R_BLR was derived from time lags between continuum and emission-line variations (Peterson et al. 2004). Thus, same as P15, we adopt the reverberation-mapping |$\mathcal {M}_{\rm BH}$| measurements with virial factor (log f = 0.71; Bentz et al. 2009; Park et al. 2012), noting that they are the anchor for the virial method and thus are inherently self-consistent.

|$\mathcal {M}_{\rm BH}$| estimates are listed in Tables 1 and 2, together with details on the emission line used. For RXJ1131, since the estimated |$\mathcal {M}_{\rm BH}$| using Mg_II and Hβ are very similar, we adopt their average. Moreover, we note that the values of |$\mathcal {M}_{\rm BH}$| for HE0435 and RXJ1131 inferred in this paper are larger than the estimates by P06 (Δlog |$\mathcal {M}_{\rm BH}$| =0.25 and 0.44 for HE0435 and RXJ1131, respectively) due to the fact that the properties of the emission lines of these two systems have been revised upwards by Sluse et al. (2012) based on data of superior quality.

5 RESULTS

Following P15 and Ding et al. (2017), for the distant objects, we adopt total uncertainty for L_host and |$\mathcal {M}_{\rm BH}$| of 0.2 dex (∼0.5 mag) and 0.4 dex, respectively.

5.1 The observed |$\mathcal {M}_{\rm BH}$|–L_host relation

The |$\mathcal {M}_{\rm BH}$|–L_bulge and |$\mathcal {M}_{\rm BH}$|–L_total relations defined by our samples are shown in Fig. 3 (panels a and b). There is a clearly positive correlation between |$\mathcal {M}_{\rm BH}$| and L_host as in local samples. For a direct comparison to local samples, we fit the local |$\mathcal {M}_{\rm BH}$|–L_host relation as:

\begin{eqnarray} \log \left(\frac{\mathcal {M}_{\rm BH}}{10^{7}\,\mathrm{M}_{{\odot }}}\right)= \alpha + \beta \log \left(\frac{L_R}{10^{10}\,\mathrm{L}_{{\odot }}}\right). \end{eqnarray}

(3)

Using a Markov Chain Monte Carlo (MCMC) process we derive α = 0.68 ± 0.18; β = 0.74 ± 0.09 for the |$\mathcal {M}_{\rm BH}$|–L_bulge and α = 0.33 ± 0.22; β = 0.95 ± 0.15 for the |$\mathcal {M}_{\rm BH}$|–L_total, with intrinsic scatter σ_int ∼ 0.25 for both of them. Consistent with previous work (e.g. P06 and P15), the observed correlation at a high redshift is nearly identical to the local. This is perhaps surprising, considering that both the black hole mass and host galaxy luminosity are expected to evolve over cosmic time. For example, in a minimal evolution toy model, the elliptical galaxies and their black hole are formed at a high redshift and evolve passively thereafter. Thus, we expect L_host to fade over time, owing to ageing stellar populations. To allow a direct comparison to the local samples, we considered this scenario in the next section.

Figure 3.

Illustration of observed (top) and evolution-corrected (bottom) correlations of |$\mathcal {M}_{\rm BH}$|–L_bulge  (left) and |$\mathcal {M}_{\rm BH}$|–L_total  (right). For distant AGNs, the redshifts are colour coded. The local data and their linear fitting (using an MCMC process) are coloured in grey (1σ region) with the best-fitted coefficients in blue colour. We use the star symbol to highlight our new lense-based measurements of HE0435 and RXJ1131. The total uncertainty for L_host and |$\mathcal {M}_{\rm BH}$| of distant AGNs is adopted to be 0.2 dex (∼0.5 mag) and 0.4 dex, respectively.

5.2 The passive evolution-corrected |$\mathcal {M}_{\rm BH}$|–L_host relation

In order to test the passive evolution scenario, we correct the observed L_host to its expected value at z = 0 by accounting for the ageing of the stellar populations. It has been shown that the evolution of the mass-to-light ratio of early-type galaxies can be effectively described as that of a single burst stellar population formed at appropriate redshifts (e.g. Treu et al. 2005). In order to represent the uncertainty in the star formation history we consider a range of single burst models formed at z_f equal to 2, 3 and 5.⁵ We choose to parametrize the evolution with the functional form dmag_R = δ_mdlog (1 + z), i.e.

\begin{eqnarray} d\log L_{\rm R}/d\log (1+z)=\delta , \end{eqnarray}

(4)

with δ = −δ_m/2.5, so that

\begin{eqnarray} \log (L_{R,0})=\log (L_{R}) - \delta \log (1+z). \end{eqnarray}

(5)

For this parametrization, we derive that δ_m ≃ −3.7 ± 0.2 (i.e. δ = 1.48 ± 0.08) provides a good representation of typical star formation histories.

This formalism is more accurate when considering a broad range in redshift with respect to adopting a single slope as a function of dmag/dz as done by P06 and P15. For a direct comparison, we also plot the passive evolution correction as a function of redshift in Fig. 4. Note that our chosen functional form describes well the P15 form at z < 1 and the P06 form at z ∼ 3 redshift, while avoiding the extreme corrections at a very high z implied by previous parametrizations. Furthermore, our chosen functional form facilitates the analysis of the |$\mathcal {M}_{\rm BH}$|–L_host evolution in the following way. Combining equation (3) with the passive evolving correction, i.e. equation (5), and adding γ^΄ term which describes the evolution of the correlation between |$\mathcal {M}_{\rm BH}$| and observed L_host, leads to the following formalism:

\begin{eqnarray} \log \left( \frac{\mathcal {M}_{\rm BH}}{10^{7}\,\mathrm{M}_{{\odot }}}\right)& = & \alpha + \beta \log \left(\frac{L_{R,0}}{10^{10}\,\mathrm{L}_{{\odot }}}\right) \nonumber \\ &&+\,\beta \delta \log (1+z) +\gamma ^{\prime } \log (1+z). \end{eqnarray}

(6)

In this equation, βδ represents the effects of passive evolution. The evolution at a fixed present-day luminosity is given by γ = γ^΄ + βδ. In this way the effects of the passive evolution correction can be easily separated and a different passive evolution model can be applied to the data, if desired. In our specific case, since we derived β = 0.74 ± 0.09 for the |$\mathcal {M}_{\rm BH}$|–L_bulge relation, the passive evolution term corresponds to approximately βδ = 0.74 × 1.48 ∼ 1.0, neglecting the effects of scatter and errors. Likewise, the passive evolution term is βδ = 0.95 × 1.48 ∼ 1.4 for the |$\mathcal {M}_{\rm BH}$|–L_total relation.

Figure 4.

Illustration of the comparison of the passive evolution correction adopted by P15, P06 and in this work. Note that all the samples in P15 are at low redshift (z ≲ 1). Thus, dmag_V/dz ≃ −1.55 is derived by assuming z_f = 2, which is appropriate at these redshifts. P06 adopted dmag_R/dz ≃ −0.8 by assuming z_f = 5.

The resulting |$\mathcal {M}_{\rm BH}$|–L_host relation after applying the passive evolution correction is shown in Fig. 3 (panels c and d). Clearly, after the correction, the high-redshift samples are offset with respect to the local samples, indicating a tendency of BH in the more distant Universe to reside in less luminous hosts at a fixed |$\mathcal {M}_{\rm BH}$|⁠. This tendency is consistent with previous work, and also consistent with the studies of the |$\mathcal {M}_{\rm BH}$|-σ_* (stellar velocity dispersion) and |$\mathcal {M}_{\rm BH}$|–|$\mathcal {M}_*$| (stellar mass) correlations, which do not require correction for passive evolution (Treu et al. 2004; Woo et al. 2006, 2008; Bennert et al. 2011a).

We fit the offset in black hole mass at a fixed passively evolved luminosity as a function of redshift in the form:

\begin{eqnarray} \Delta \log \mathcal {M}_{\rm BH}= \gamma \log (1 + z), \end{eqnarray}

(7)

where |$\Delta \log \mathcal {M}_{\rm BH} = \log ( \frac{\mathcal {M}_{\rm BH}}{10^{7}\,\mathrm{M}_{{\odot }}}) -\alpha -\beta \log (\frac{L_{R,0}}{10^{10}\,\mathrm{L}_{{\odot }}}),$| and obtain γ = 0.75 ± 0.11 for the |$\mathcal {M}_{\rm BH}$|–L_bulge and γ = 0.95 ± 0.11 for the |$\mathcal {M}_{\rm BH}$|–L_total relations, as shown in Fig. 5 (panels a and b). We also obtain γ^΄ = −0.14 ± 0.11 for the |$\mathcal {M}_{\rm BH}$|–L_bulge and γ^΄ = −0.26 ± 0.12 for the |$\mathcal {M}_{\rm BH}$|–L_total correlation, when not taking into account the passive evolving correction. As expected, the difference γ − γ^΄ is consistent with the effects of the passive evolving correction, i.e. βδ ∼ 1.0 for the |$\mathcal {M}_{\rm BH}$|–L_bulge and βδ ∼ 1.4 for the |$\mathcal {M}_{\rm BH}$|–L_total relation.

Figure 5.

Illustration of the offset in log|$\mathcal {M}_{\rm BH}$| for a given L_bulge (left) and L_total (right) as a function of redshift, after passive evolution correction. Top panels correspond to the fitting using the whole sample. We also highlight the subsamples from SS13. Bottom panels correspond to the fitting excluding the samples from P06 and SS13. The red solid line represents the best-fitting trend for all distant objects as a function of Δlog |$\mathcal {M}_{\rm BH}$| = γlog (1 + z), with the 1σ region confidence range shaded in grey. The orange band is the intrinsic scatter of local linear relation.

We conclude by noting that this fit does not take into account selection effects, which are discussed in the next section.

5.3 Selection effects

From Fig. 3, we can see that at a high redshift we preferentially study systems with the larger |$\mathcal {M}_{\rm BH}$| and L_total. This is expected as observational samples tend to be flux limited and thus favour the high-luminosity tail (and hence typically high |$\mathcal {M}_{\rm BH}$|⁠) of the distribution. Like many other instances in astronomy, it is essential to take into account the selection function when estimating the evolution of the black hole mass host galaxy correlations (Lauer et al. 2007; Treu et al. 2007; Bennert et al. 2011a; Schulze & Wisotzki 2014; P15).

Following P15, we take selection effects into account by using a Monte Carlo simulation method based on the methodology introduced by Treu et al. (2007) and Bennert et al. (2010). The simulated samples are generated from a combination of the local active BH mass function from Schulze & Wisotzki (2010) and the local |$\mathcal {M}_{\rm BH}$|–L_host relation from Bennert et al. (2010) with Gaussian random noise added as a function of the two free parameters γ and intrinsic scatter of the correlation σ_int. Note that the scatter is assumed to be independent of redshift in our description. For each object, the likelihood of the observed |$\mathcal {M}_{\rm BH}$| with a given L_host is calculated from the simulated sample at the given γ and σ_int, and taking into account whether the object would be selected or not based on our sensitivity. Finally, by adopting uninformative uniform (flat) prior or lognormal prior from Bennert et al. (2010, σ_int = 0.21 ± 0.08), the posterior distribution function of γ and σ_int is evaluated. Selection effects are modelled in the same way for the lensed-quasar sample, neglecting any second-order effects related to lensing magnification. We note, however, that these effects are small (Collett & Cunnington 2016) and magnification-related biases should affect the quasar and host galaxy in a similar manner, thus moving objects mostly along the |$\mathcal {M}_{\rm BH}$|–L_host correlation and not away from it.

Taking into account selection effects, the results of the inference are shown in Fig. 6. The fitted values of γ are 0.6 ± 0.1 (⁠|$\mathcal {M}_{\rm BH}$|–L_bulge) and 0.8 ± 0.2 (⁠|$\mathcal {M}_{\rm BH}$|-L_total), almost independent of the choice of prior. These values are consistent with the previous inference in Sections 5.1 and 5.2.

Figure 6.

Posterior distribution function given the entire data set for a model with evolution in the form Δlog |$\mathcal {M}_{\rm BH}$| = γlog (1 + z) with intrinsic scatter σ_int, taking into account selection effects. The |$\mathcal {M}_{\rm BH}$|–L_bulge  (top) and |$\mathcal {M}_{\rm BH}$|-L_total  (bottom) correlations with flat (left) and lognormal prior (right) are shown.

Interestingly, the intrinsic scatter of the correlations is found to be consistent with typical values inferred for local samples (0.3–0.4 dex). This result is consistent with the hypothesis that well-defined correlations exist at the redshifts probed by our sample, and indicates that we have not significantly underestimated our errors at a high z. It would be beneficial to study how the selection bias changes as a function of some key factors such as the values of L_host and |$\mathcal {M}_{\rm BH}$|⁠, the level of the uncertainties and the redshift distribution of the samples. However, this topic is trivial in this study as we obtained consistent inference by either or not talking selection effects into account. Moreover, to study this relation quantitatively requires considerate tests and simulations. Thus, we leave it for the future study.

6 DISCUSSION

In this section, we first estimate the importance of potential systematic errors in Section 6.1. Then, we carry out a detailed comparison with previous observational work in Section 6.2. Finally, we discuss how our measurements fit into our understanding of galaxies and BHs co-evolution in Section 6.3.

6.1 Systematic errors

We have combined our new measurements with ones taken from the literature in order to increase the sample size and reduce statistical uncertainties. Even though we have restricted our analysis to the samples that have been analysed in the most similar manner to our new data and we have cross-calibrated the black hole mass estimators, there are still some residual differences.

First, P06 obtained the luminosity of one galaxy by combining the fluxes together, even though some of them may include a disc component (e.g. RXJ1131). According to morphological studies of AGN host galaxies, the fraction of spiral/elliptical hosts of AGN is approximately one third (Kocevski et al. 2012), with the exact value depending on |$\mathcal {M}_{\rm BH}$| and luminosity. Thus, it is possible that P06 overestimates the bulge component of some of the host galaxies. The total luminosity should be less affected by this bias, even though not completely immune.

Furthermore, the subsample by SS13 included in the compilation by P15 was X-ray selected as opposed to optically selected like the rest of the non-lens sample (some of the lenses are radio-selected.). This difference in selection could potentially lead to a systematic difference between the two samples.

In order to estimate these systematic uncertainties, we repeat the analysis by excluding the P06 and SS13 samples. This reduced sample will have significantly less statistical power, owing to the reduced size and redshift coverage, but should be more robust with respect to the systematic uncertainties discussed above. Given this reduced sample, we obtain γ = 0.88 ± 0.28 for the |$\mathcal {M}_{\rm BH}$|–L_bulge and γ = 0.51 ± 0.28 for the |$\mathcal {M}_{\rm BH}$|–L_total relations, as shown in Fig. 5 (panels c and d). Moreover, we use the same approach to study the selection effects and obtain the consistent inference, as illustrated in Fig. 7. Even though as expected the uncertainties are larger than for the full sample, the results are statistically mutually consistent at 1σ level. To facilitate the comparison between different γ, we summarize our inference in Table 4. We conclude that our inferred trends are not dominated by systematic differences between the samples, and systematic uncertainties of this kind are smaller than the random ones.

Figure 7.

Same as Fig. 6, using the reduced sample.

Table 4.

The summary for the different inference of γ.

Sample	Account for	|$\mathcal {M}_{\rm BH}$|-Lbulge	|$\mathcal {M}_{\rm BH}$|-Ltotal
	selection effects?
Entire	No	0.75 ± 0.11	0.95 ± 0.11
Entire	Yesa	0.6 ± 0.1	0.8 ± 0.1
Exclude P06, SS13	No	0.88 ± 0.28	0.51 ± 0.28
Exclude P06, SS13	Yes	0.7 ± 0.4	0.2 ± 0.5

^aUsing lognormal prior.

In this work, |$\mathcal {M}_{\rm BH}$| estimates are derived using the C_IV, Mg_II and Hβ emission lines. However, the C_IV and Mg_II lines are usually in outflow (Baskin & Laor 2005; Richards et al. 2011; Denney 2012) and therefore may not be dominated by the gravity of the central |$\mathcal {M}_{\rm BH}$| and result in biased |$\mathcal {M}_{\rm BH}$| estimates, especially for the C_IV line (Trakhtenbrot & Netzer 2012). Following McGill et al. (2008), the potential bias has been mitigated by cross-calibrating the |$\mathcal {M}_{\rm BH}$| estimates based on the different lines. As a further sanity check, we fitted the γ using only Hβ-based samples. We note that this Hβ sample is very similar to the subsample excluding those of P06 and SS13, and in fact the results are similar (γ = 1.10 ± 0.36 for the |$\mathcal {M}_{\rm BH}$|–L_bulge and γ = 0.7 ± 0.37 for the |$\mathcal {M}_{\rm BH}$|-L_totalrelation). We conclude that any potential residual bias related to the use of lines other than Hβ is smaller than statistical uncertainties or biases related to sample selection.

6.2 Comparison with previous work

P15, using a sample of 79 active galaxies, inferred the following evolutionary trends: γ = 0.9 ± 0.7 for the |$\mathcal {M}_{\rm BH}$|-L_bulge and γ = 0.4 ± 0.5 for the |$\mathcal {M}_{\rm BH}$|-L_totalrelation. These are consistent with our inference, although their uncertainties are much larger, owing to the smaller sample size and reduced high-redshift coverage. A similar result was obtained by P06, where they found that the ratio between |$\mathcal {M}_{\rm BH}$| and |$\mathcal {M}_*$| was approximately four times larger at z ∼ 2–4 than today (i.e. γ ∼ 0.8 − 1.2). The consistency between their measurements and ours is expected since the overall samples in this work are mostly composed of the samples by P15 and P06, even though there are some differences in the rest-frame bands chosen for photometry (we and P06 adopt rest-frame R, while P15 adopts rest-frame V.), in the passive evolution correction, and in the black hole mass calibration.

The cosmic evolution of the |$\mathcal {M}_{\rm BH}$|–L_bulge relation is a topic of intense debate in the literature. Many works have reported an evolutionary signal based on different relations including the |$\mathcal {M}_{\rm BH}$|–L_host (e.g. Treu et al. 2007; Bennert et al. 2010), the |$\mathcal {M}_{\rm BH}$|-|$\mathcal {M}_*$| (e.g. B11; McLure et al. 2006; Jahnke et al. 2009; Decarli et al. 2010; Cisternas et al. 2011; Trakhtenbrot et al. 2015) and the |$\mathcal {M}_{\rm BH}$|–σ_* (e.g. Woo et al. 2006, 2008) correlations. Nevertheless, other observational studies (e.g. Shields et al. 2003; Greene & Ho 2005; Komossa & Xu 2007; Shen et al. 2008) found no evidence for evolution. In Shankar et al. (2016), they find serious biases in the |$\mathcal {M}_{\rm BH}$|-|$\mathcal {M}_*$| relation and prove that σ_* is more fundamental than any other variable. However, in Shankar, Bernardi & Haiman (2009), they show that there is no evolution in the |$\mathcal {M}_{\rm BH}$|–σ_* relation once one accounts for the ages of local galaxies and the Sołtan argument. Moreover, Schulze & Wisotzki (2011, 2014) concluded that there is no statistically significant evidence for evolution once these selection effects are taken into account and corrected. Taking a different approach, DeGraf et al. (2015) used the results of the high-resolution numerical simulation MassiveBlackII to compare the observed and intrinsic evolution of the black hole mass host galaxy correlations and reproduced the evolutionary trend of the relation. Consistent with other considerations, they also found that the observed samples display steeper slopes than random ones, suggesting the selecting effects can exhibit faster evolution than a random sample. Similarly, by generating Monte Carlo realizations of the |$\mathcal {M}_{\rm BH}$|-σ_* relation at z = 6, Volonteri & Stark (2011) also found that due to selection bias the’observable’ subsample would suggest an average positive evolution even when the intrinsic correlation is characterized by no or negative evolution at a high redshift. These studies highlight once again the importance of taking selection effects into account.

Clearly, the absence of evidence does not imply evidence of absence, and one way to make progress is to improve the precision and accuracy of the measurement. In our work, we attain much higher precision than previous work owing to the enlarged sample, including lensed quasars. Thanks to the large sample size, even when selection effects are taken into account, the evolutionary trend is detected at high significance (γ ≠ 0 at more than 5σ). However, using a reduced sample by excluding the subsamples from P06 and SS13, we obtain a smaller evolutionary trend, with larger uncertainties. These results are consistent at 1σ level (see Table 4), and highlight the importance of studying a larger sample of high-redshift lenses with the state-of-the-art data.

6.3 Implication for the co-evolution of black holes and their host galaxies

Our results are consistent with a scenario in which BHs in the distant Universe typically reside in lower stellar mass galaxies than today, assuming that the passively evolved luminosity tracks approximately stellar mass (see SS13; Bennert et al. 2011a, for a consistent direct measurement based on stellar mass determination). In order to end up on the local final relation, the stellar mass of the host galaxy would have to grow faster than |$\mathcal {M}_{\rm BH}$|⁠.

An interesting clue to the physical mechanism driving the evolution could perhaps be found by comparing the inferred evolution for the correlation between |$\mathcal {M}_{\rm BH}$| and the total host galaxy luminosity, and that with the bulge luminosity. We found those two to be comparable within the uncertainties. Previous work found the |$\mathcal {M}_{\rm BH}$|–|$\mathcal {M}_*$|_bulge correlation to evolve somewhat faster than |$\mathcal {M}_{\rm BH}$|-|$\mathcal {M}_*$|_total, albeit at low statistical significance, suggesting that one of the mechanisms at work is the build up of the bulge component from stars in the disc (Croton 2006). It is difficult to perform a direct comparison because of the fact that the P06 sample did not attempt bulge–disc decomposition, and we have also assumed a single passive evolution trend for the entire galaxy. Both effects could potentially suppress the differences between the evolution of the bulge and total luminosity with respect to the black hole mass.

Also, our sample extends to much larger redshift than that of B11. One possible explanation of this possible tension is that the dominant evolutionary mechanisms change with redshift. At a low redshift (z ≲ 1), the growth of the bulge is dominated by the secular evolution with the redistribution of disc stars while at a high redshift (z ≳ 1), the growth is dominated by major mergers (see Bennert et al. 2010, for a similar conjecture).

To settle this issue, it is crucial to obtain high-quality data and model large samples of lens systems, so that robust bulge to total luminosity decompositions can be carried out. It would also be beneficial to obtain multi-colour data to estimate directly stellar mass, and ideally stellar kinematic information to distinguish pressure supported systems from rotationally supported ones.

Since the L_host of ellipticals would not change when considering the bulge and the total, we examine the offset using the sample limited to spiral galaxies. In our sample, there are nine local spirals and 41 distant spirals, excluding SS13. Fitting the offset with this subsample, we obtain γ = 2.15 ± 0.41 and 1.18 ± 0.41 for |$\mathcal {M}_{\rm BH}$|–L_bulge and |$\mathcal {M}_{\rm BH}$|-L_total relations, respectively, which are larger than the previous inference listed in Table 4 for the entire sample. This difference could suggest the spiral galaxies are undergoing a more rapid evolution than the ellipticals in order to end up on the local final relation. However, we caution that this result should be taken with a grain of salt, given the small sample size of the local disc comparison sample.

7 SUMMARY

We presented a new measurement of the co-evolution of supermassive black holes and their host galaxies. First, we carried out a new analysis of two strongly lensed quasars, HE0435−1223 and RXJ1131−1231. By using the state-of-the-art lens models by Wong et al. (2017) and Suyu et al. (2013), we found that the host galaxies of HE0435 and RXJ1131 are well described by an elliptical and spiral surface brightness density profile, respectively. Then, we measured the host galaxy magnitude and tested for potential biases by carrying out realistic simulations following the procedure outlined by Ding et al. (2017). We found that the bias of our inference of L_host is small (0.1−0.2 mag) and that we can recover the host image precisely even if the host has multiple components (see Fig. 2, panel c). We estimated |$\mathcal {M}_{\rm BH}$| by using a set of self-consistent single-epoch estimators based on the quasar emission-line properties as measured by Sluse et al. (2012).

Secondly, we combined our measurements with the published ones from the literature (P15; P06), thus expanding our sample to 146 active galaxies up to z = 4.5. We have taken care of using self-consistent recipes to re-derive the black hole mass estimates and convert all the luminosities self-consistently to the rest-frame R band.

Our main findings can be summarized as follows.

• The observed correlations – without correction for evolution – are consistent with those observed in the local Universe.

• The data are inconsistent with a passive evolution scenario. By correcting the host galaxy rest-frame luminosity to z = 0, we find that galaxies are underluminous for a given |$\mathcal {M}_{\rm BH}$|⁠, even neglecting growth by accretion.

• The passively evolved correlations are well described by a relationship of the form Δlog |$\mathcal {M}_{\rm BH}$| = γlog (1 + z) with γ = 0.6 ± 0.1 and γ = 0.8 ± 0.2, respectively, at fixed bulge and total host luminosity, taking into account selection effects.

Considering that stellar populations must fade as they get older, and considering that similar results have been found when studying the correlations between |$\mathcal {M}_{\rm BH}$| and host galaxy velocity dispersion (Treu et al. 2004; Woo et al. 2006, 2008) and stellar mass (Jahnke et al. 2009; SS13; B11), we are forced to conclude that the co-evolution of galaxies and black holes is non-trivial in the sense that systems do not stay on the correlation as they evolve. At least for active galaxies in the range of black hole and stellar masses that can be analysed with current technology, it appears that the growth of the black hole predates that of the bulge (Croton 2006). However, given the complexity and variety of processes involved, direct comparisons with detailed numerical simulations are needed to further our understanding of the co-evolution of black holes and their hosts. Recent cosmological simulations including some prescriptions for black hole growth and feedback have been shown to reproduce the observations at least at z < 1 (DeGraf et al. 2015). It will be interesting to carry out similar detailed comparisons, taking into account errors and observational selection functions, for a variety of models (e.g. Sijacki et al. 2015; Taylor & Kobayashi 2016; Volonteri et al. 2016) and extending to higher redshifts. These comparisons will provide a powerful test of the various recipes that have been adopted to describe accretion and star formation physics at sub-grid level in numerical simulations.

Looking at the future, the sample of lensed quasars that can be analysed with high fidelity is going to grow. Currently, ultra deep HST imaging data have been obtained for six additional strongly lensed systems⁶ and their analysis will be described in a forthcoming paper. The sample of lensed quasars and their hosts that can be studied at high fidelity is likely to continue to grow as more such systems are discovered in wide field imaging and spectroscopic surveys (e.g. Agnello et al. 2015; More et al. 2016; Schechter et al. 2017; Ostrovski et al. 2017).

Acknowledgements

This work is based in part on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with programs 9744, 12889, 14254. Financial support was provided by NASA through grants from the Space Telescope Science Institute.

We are grateful to Vivien Bonvin, Geoff C.-F. Chen, Frederic Courbin, Matthew A. Malkan, Cristian E. Rusu, Jong-Hak Woo and Andreas Schulze for useful comments and suggestions that improved this manuscript. We thank Chien Peng for his help with the estimates of black hole mass. XD is supported by the China Scholarship Council. TT acknowledges support by the Packard Foundations through a Packard Research Fellowship and by the NSF through grants AST-1450141 and AST-1412315. SHS gratefully acknowledges support from the Max Planck Society through the Max Planck Research Group. CER acknowledges support from the NSF grant AST-1312329. DS acknowledges funding support from a Back to Belgium grant from the Belgian Federal Science Policy (BELSPO). KCW and DP are supported by an EACOA Fellowship awarded by the East Asia Core Observatories Association, which consists of the Academia Sinica Institute of Astronomy and Astrophysics, the National Astronomical Observatory of Japan, the National Astronomical Observatories of the Chinese Academy of Sciences, and the Korea Astronomy and Space Science Institute. VNB gratefully acknowledges assistance from a National Science Foundation (NSF) Research at Undergraduate Institutions (RUI) grant AST-1312296. Note that findings and conclusions do not necessarily represent the views of the NSF. VB acknowledges the support of the Swiss National Science Foundation (SNSF).

1

Developed by Suyu & Halkola (2010) based on Suyu et al. (2006) and Halkola et al. (2008).

2

H₀ lenses in COSMOGRAIL's Wellspring, http://www.h0licow.org/.

3

Suyu et al. (2014) updated the model of RXJ1131. Given the similarity in the composite and power-law model by Suyu et al. (2014), a similar time delay distance is obtained (within ∼2 per cent, and hence spatial scaling of the source due to mass–sheet degeneracy). This means the inference of total flux of the host should be unchanged to within ∼4 per cent.

4

lensfit is a version of galfit that has been extended to fit lensed host galaxies while optimizing the mass model for the lens galaxy. For details, see P06.

5

Stellar evolution is calculated with Galaxev (Bruzual & Charlot 2003), based on Padova-1994 stellar evolutionary tracks, assuming Salpeter IMF, solar metallicity and no dust attenuation.

6

WFI2033−4723 and HE1104−1805 from H0LiCOW program HST-GO-12889 (PI: Suyu); SDSS1206+4332, HE0047−1756, SDSS0246−0825 and HS2209+1914 as part of Program HST-GO-14254 (PI: Treu).

APPENDIX A: THE K-CORRECTION FOR THE RXJ1131 HOST

We apply the K-correction to the observed magnitudes to obtain the rest-frame R-band magnitude. At the redshift of RXJ1131 (z_s = 0.654), the conversion from F814W to R-band magnitude depends on the adopted SED templates, as shown in Fig. A1 (panel a). Therefore, we directly determine the K-correction of the disc component through SED fitting with the multi-band images (F555W, F814W, F160W), available in the archive (GO-9744; PI: C. S. Kochanek).

Figure A1.

Illustration of the K-correction for RXJ1131. (a) K-correction from observed F814w band to the rest-frame R-band magnitude, as a function of stellar population age. The colours are calculated based on Galaxev stellar evolution track with metallicities of Z = 0.4Z_⊙, 1.0Z_⊙, 2.5Z_⊙ for E-type and Sb-type galaxies. (b) Colour map (the rest-frame R magnitude – observed F814W magnitude) of RXJ1131, calculated by SED fitting with three broad-band imaging (F555W, F814W and F160W). For bulge region, a direct measurement of SED, is affected by the residual AGN contamination, and hence half blue and half red. For the disc region, i.e. around the area where lensing-distorted spiral arm patterns and star-forming regions are clearly visible, Δmag is approximately −0.3.

Briefly, we fit the SED for each pixel by using Fast (Kriek et al. 2009) based on Galaxev stellar evolution track (Bruzual & Charlot 2003), assuming the solar metallicity, exponentially declining star formation history and Calzetti et al. (2000) dust extinction law, while the redshift is fixed to the spectroscopic one. The error for each pixel is calculated based on empty regions of images. We then derive the rest-frame R-band magnitude by using the best-fitting template for each pixel, and see the offset from the observed F814W magnitude. As shown in Fig. A1, panel (b), Δmag_disc(R −F814W) ≈ −0.3 is an appropriate estimate for the disc region. For the bulge, a direct measurement of SED is affected by the residual AGN contamination, and hence the arc corresponding to the bulge is half blue and half red. Thus, we assume that the age of stellar populations is >3 Gyr, where the correction hardly changes (see Fig. A1, panel a), and adopts Δmag_bulge(R −F814W) ≈ −0.7.