© ESO, 2011

1. Introduction

Supermassive black holes (SMBHs) appear to be ubiquitous in the centers of massive galaxies (Kormendy & Richstone 1995). Furthermore, they show tight correlations with the properties of the spheroidal components of their host galaxies, e.g. with the stellar velocity dispersions (Ferrarese & Merritt 2000; Gebhardt et al. 2000; Tremaine et al. 2002; Gültekin et al. 2009), bulge luminosities and masses (Magorrian et al. 1998; Marconi & Hunt 2003; Häring & Rix 2004; Sani et al. 2011), or concentration indices (Graham et al. 2001). These relations have been established by direct dynamical measurement of the black hole mass in a few dozens of nearby galaxies, mainly using stellar dynamics (e.g. van der Marel et al. 1998; Emsellem et al. 1999; Gebhardt et al. 2003) or gas dynamics (e.g. Ferrarese et al. 1996; Marconi et al. 2001; Dalla Bontà et al. 2009).

Nevertheless, significant uncertainties in these relations persist. The high mass and the low mass regimes are still poorly constrained. There may also be systematic effects on the determined masses; for example, not accounting for triaxiality of the galaxy (van den Bosch & de Zeeuw 2010) or not including a dark matter halo in the dynamical models (Gebhardt & Thomas 2009; Schulze & Gebhardt 2011) can bias the dynamical measurement. Also the intrinsic scatter in these relations is still not well established (Gültekin et al. 2009). Furthermore, there is evidence that late-type galaxies and in particular pseudobulges do not follow the same relations as early-type galaxies (Hu 2008; Graham 2008; Greene et al. 2010; Kormendy et al. 2011; Jiang et al. 2011).

The black hole-bulge relations contain important information, as they imply a connection between the growth of a black hole and the evolution of its host galaxy. A common picture invokes AGN feedback to shut down star formation and self-regulate black hole accretion. Models building on this scenario are able to reproduce the local black hole-bulge relations in numerical simulations (Di Matteo et al. 2005; Sijacki et al. 2007; Booth & Schaye 2009) and semi-analytically (Kauffmann & Haehnelt 2000; Cattaneo et al. 2005; Croton et al. 2006; Bower et al. 2006; Somerville et al. 2008; Marulli et al. 2008). A correlation between black hole mass and host galaxy mass will also be tightened or may even be generated naturally within a merger driven galaxy evolution framework (Peng 2007; Hirschmann et al. 2010; Jahnke & Macciò 2011).

Essential constraints on the origin of the black hole-bulge relation can be inferred from their redshift evolution. There are several theoretical predictions, based on numerical simulations (Robertson et al. 2006; Hopkins et al. 2007; Di Matteo et al. 2008; Johansson et al. 2009; Booth & Schaye 2011; Dubois et al. 2011) as well as on semi-analytic models (Wyithe & Loeb 2003; Croton 2006; Malbon et al. 2007; Hopkins et al. 2009; Lamastra et al. 2010). Although the details are still far from being settled, they tend to predict an increase in the M_•/M_Bulge ratio with redshift, while only a weak or even negative evolution in the M_•-σ _∗ relation is expected.

These models need to be confronted with observations. Several different approaches have been followed in the last years to observationally constrain the evolution in the black hole-bulge relations. These range from more indirect arguments, to direct estimates of M_• and the respective bulge property. Constraints on integrated quantities can be gained from the black hole mass function (Merloni et al. 2004; Hopkins et al. 2006; Shankar et al. 2009; Somerville 2009; Kisaka & Kojima 2010; Willott et al. 2010). Bluck et al. (2011) studied X-ray selected AGN, employing Eddington ratio arguments, to constrain black hole-bulge coevolution.

Direct dynamical determinations of M_• are not feasible outside of the very local universe. Therefore, the most commonly adopted approach resorts to broad line AGN, for which black hole masses are easily accessible through the so-called “virial estimator” (e.g., McLure & Jarvis 2002; Vestergaard & Peterson 2006). The main challenge for the determination of the M_•-bulge relationship from AGN samples is the determination of the bulge properties, hampered by the bright nuclear point source of the AGN, which may well outshine the entire galaxy.

Evolution in the M_•-σ _∗ relation has been studied either by measuring stellar velocity dispersions directly (Woo et al. 2006, 2008; Shen et al. 2008) or by using the widths of narrow emission lines as surrogates for σ _∗  (Shields et al. 2003, 2006; Salviander et al. 2007). The M_•-M_Bulgerelation has been investigated from QSO host galaxy luminosities (Peng et al. 2006a,b; Treu et al. 2007; McLeod & Bechtold 2009; Decarli et al. 2010; Bennert et al. 2010; Targett et al. 2011), by estimating the stellar masses utilising colour information (Schramm et al. 2008; Jahnke et al. 2009; Merloni et al. 2010), or by dynamical mass measurements (Inskip et al. 2011). Also obscured AGN with detectable broad lines have been used (Sarria et al. 2010; Nesvadba et al. 2011), for which the determination of stellar masses is less problematic. McLure et al. (2006) inferred the M_•-M_Bulgeratio up to z = 2 by matching distributions of radio-loud QSOs and radio galaxies. At the highest redshifts, CO rotation curves have been used to determine dynamical masses for a few individual objects (Walter et al. 2004; Riechers et al. 2008, 2009; Wang et al. 2010). All these methods have their own advantages and drawbacks. Nevertheless, there is increasing evidence for apparent evolution, in a sense of high-z SMBHs being more massive at given bulge mass than in the local universe.

However, it has been suspected that these results may be biased due to sample selection effects. One such bias could arise from the fact that luminous AGN on average tend to have massive black holes. The intrinsic scatter in the BH-bulge relation then generates a bias towards a higher M_•/M_Bulge ratio. Such a bias has been pointed out by several authors (Adelberger & Steidel 2005; Fine et al. 2006; Salviander et al. 2007), but received major attention by the work of Lauer et al. (2007) who discussed it in some detail (but see also Peng 2010). A further effect recently discussed by Shen & Kelly (2010) is caused by the uncertainty in virial black hole mass estimates in connection with the steep decrease of the active black hole mass function.

In this paper we argue that besides the above-mentioned effects there are even more potential biases that need to be considered. While previous authors mostly focused on one particular source of bias, we aim at providing a comprehensive study of selection effects on observations of the M_•-M_Bulgerelation and its evolution. We also investigate how well these biases can be quantified and corrected for. The paper is organised as follows: in Sect. 2 we present our general framework to investigate the consequences of selection effects and show a first application to the quiescent black hole sample. Section 3 discusses several selection effects for AGN samples and their ramifications for the black hole-bulge relations. In Sect. 4 we take into account redshift evolution effects. We discuss the implications of our results on observational studies in Sect. 5. We finally conclude in Sect. 6. For the cosmological parameters we assume H₀ = 70 km s^-1 Mpc^-1, Ω_m = 0.3 and Ω_Λ = 0.7.

2. The local M_•-M_Bulge relation

2.1. General framework

We start with some basic thoughts on the BH-bulge relation and its sensitivity to selection effects, postponing specific aspects of AGN samples to the next section. We adopt the following convention: Black hole masses are given by μ = log M_•, while any of the various relevant bulge properties is denoted as s, with s = log σ _∗  or s = log M_Bulge, respectively. In some examples we focus on the M_•-M_Bulgerelation, but the arguments equally hold for the M_•-σ _∗ relation. We explicitly do not imply one relation to be more fundamental than the other.

The distribution of objects in the μ-s diagram is given by the bivariate distribution function of bulge property and black hole masses Ψ(s,μ). Thus Ψ(s,μ)dμds gives the number of objects per Mpc³ with galaxy property between s and s + ds and black hole mass between μ and μ + dμ. We assume the following parameterisation for most of the paper, (1)where Φ_s(s) is the spheroid distribution function, e.g. the spheroid mass function, and g(μ | s) gives the probability of finding the black hole mass μ, given s. If s and μ are correlated, as suggested by the observations, then g(μ | s) corresponds to this correlation. In the following we assume a linear relation μ = a + bs with log-normal intrinsic scatter σ, i.e. (2)This parameterisation for the bivariate distribution function is not the only one possible. We comment on the observational justification of this form in the Appendix. While our quantitative results depend on the adopted parameterisation, most of our qualitative results are independent of the specific choice of Ψ(s,μ).

The corresponding distribution functions for s and μ are given by marginalisation over the other variable, (3)(4)The integration over μ simply returns the spheroid mass function. The integration over s provides the quiescent black hole mass function (BHMF), equivalent to the common approach for its determination (e.g. Yu & Lu 2004; Marconi et al. 2004; Merloni 2004). The first equality in Eqs. (3) and (4) is independent of the special choice of Ψ(s,μ). In particular, as the galaxy distribution function is an observable, Eq. (3) sets a constraint on the bivariate distribution Ψ(s,μ). Note that the specific shape of Ψ(s,μ) has direct consequences for the quiescent black hole mass function.

Figure 1 shows Ψ(s,μ), Φ_s(s) and Φ_•(μ) for the local universe, employing our parameterisation. We estimated the spheroid mass function from the early type and late type stellar mass functions of Bell et al. (2003). We computed the spheroid mass function as sum of the early type mass function and the mass function of the bulge components of late type galaxies. Since the bulge fraction in disk galaxies is still somewhat uncertain, we simply assumed an average value of B/T = 0.3 to convert the late type mass function into a bulge mass function, consistent with current observations (e.g. Graham & Worley 2008; Gadotti 2009). This value also produces a black hole mass function that is consistent with determinations from the M_•-σ_∗and M_•-L relations (Marconi et al. 2004). For the M_•-M_Bulgerelation we used the relation by Häring & Rix (2004) and an intrinsic scatter of σ = 0.3 dex. For most of this paper we keep this value for the intrinsic scatter fixed; we comment on the consequences of changing it when appropriate.

The distribution shown in Fig. 1 would be obtained for a strictly volume limited sample of galaxies with M_Bulge and M_• measurements. In practice, the distribution will be affected by the way that the sample was constructed. This is formally accounted for by a selection function Ω, defined as the probability of observing an object of a given bulge mass, black hole mass and further possible selection criteria, such as for instance redshift or AGN luminosity. The selection function can thus be written as Ω(s,μ,θ), where θ refers to the set of additional parameters present as selection criteria. The observed apparent bivariate distribution is then given by (5)where p_θ(θ) is a set of normalised distribution functions of the parameters θ. Depending on Ω, the observed M_•-M_Bulgerelation can be significantly biased if the selection effects are not taken into account. However, this requires a proper knowledge of the selection function.

The bivariate distribution Ψ_o(s,μ), when normalised to one, represents the full probability distribution of the expected M_•-M_Bulge relation. However, it may happen that we are not so much interested in the full distribution, but rather in the mean relation. This is given by the mean black hole mass at a given bulge property (6)If no selection effects are present then  ⟨ μ ⟩ (s) will be identical to the input M_•-bulge relation,  ⟨ μ ⟩  = a + bs. Alternatively we can consider the mean bulge property at a given black hole mass (7)Even without any selection effects,  ⟨ s ⟩ (μ) does not simply correspond to the inverse relation. It deviates from it quite strongly at the high mass end where the space density in the galaxy distribution function is decreasing. This effect is illustrated as the red dot-dashed line in Fig. 1 and was already discussed by Lauer et al. (2007). It is not an observational bias but just a direct consequence of the different projections chosen, combined with the steep decrease of the galaxy distribution function. In general, under the presence of selection effects  ⟨ s ⟩ (μ) may change as well.

In observational studies of the evolution of the M_•-bulge relation, the sample is often compared to the local relation in terms of a single value, the offset from the local relation. The sample bias for this offset is given by (8)Alternatively, the mean offset of the galaxy property from the local relation is given by (9)where we use the inverse relation s = α + βμ with α =  −a/b and β = 1/b. This directly leads to  ⟨ Δs ⟩  =  −1/b ⟨ Δμ ⟩ , without any specific assumption on Ψ_o(s,μ). Thus the biases in  ⟨ Δμ ⟩  and  ⟨ Δs ⟩  should always be directly proportional to each other.

Already the well-studied samples of a few dozen local galaxies with reasonably well-measured dynamically-based black hole masses will invariably be affected by selection effects. However, these samples are very inhomogeneous, and presumably they cannot be described by a single well-defined selection function, which hampers a proper investigation of potential selection effects. Yu & Tremaine (2002) and Bernardi et al. (2007) reported on a possible bias in the dynamical black hole mass sample identified through finding a discrepancy in the σ_∗ − L relation compared to the SDSS. More recently, Gültekin et al. (2009) discussed a bias induced by culling the sample based on resolving the spheres of influence of the black holes. Before discussing selection effects on AGN samples in detail, we investigate this bias to illustrate the general applicability of Eq. (5).

2.2. Bias by sphere of influence resolution

Dynamical black hole mass determinations need to spatially resolve the region of gravitational influence of the central black hole on the stellar velocity distribution. The size of this region is commonly estimated by the “sphere of influence” of radius . It has been argued that dynamical black hole masses are unreliable for R_inf/d_res < 1 and therefore should be excluded from the sample (e.g. Ferrarese & Ford 2005). Here d_res is the spatial resolution of the observations. Gültekin et al. (2009) demonstrated that rejecting galaxies below a line leads to a bias in the M_•-bulge relations. Batcheldor (2010) pointed out that it would even be possible to artificially generate a M_•-σ_∗correlation via this selection effect, while Gültekin et al. (2011) excluded this possibility at least for the most extreme case. Gültekin et al. (2009) argued against culling the sample based on the sphere of influence to avoid this bias. However, even without any active rejection the same effect may be effectively present, either due to an implicit target selection or through the mere ability to detect a black hole. The condition R_inf/d_res ≥ 1 is no strict limit for the reliability of a black hole detection, but it is clear that for R_inf/d_res ≪ 1 no black hole can be detected dynamically.

We now illustrate the results of this bias. We assume that we can detect all black holes above a certain threshold in resolution of the sphere of influence, R_inf/d_res > r_min. We keep the angular resolution of our hypothetical survey fixed at 0.1. The selection function is then given by (10)where d is the distance of the galaxy. The observed bivariate probability distribution follows as (11)For the purpose of illustration we again use the M_•-M_Bulgediagram. We convert bulge masses into velocity dispersions assuming a fixed mass-to-light ratio for the r-band and using the σ _∗  − L_r relation from Bernardi (2007). We cover the distance range from 1 to 30 Mpc. In Fig. 2 we show the M_•-M_Bulgediagram for two thresholds in the sphere of influence resolution, r_min = 1 and r_min = 0.1. We additionally plot the “observed” mean M_•-M_Bulgerelation  ⟨ μ ⟩ (s) for both cases as thick black solid line and the “observed” mean relation  ⟨ s ⟩ (μ) as red dashed-dotted line. The bias induced by whether or not resolving the sphere of influence leads to a flatter slope, a higher normalisation and a smaller intrinsic scatter, as already shown by Gültekin et al. (2009). See also their work for an extensive Monte Carlo investigation of this bias.

While a significant bias is induced only for a high threshold r_min, the probability distribution is affected even for lower values, most severely at the low mass end. Also the relation  ⟨ s ⟩ (μ) is biased in both cases. The selection along lines of is clearly visible in the probability contours. However, it is hard to quantify how much the actually observed samples of black holes are actually biased by this selection effect, due to their inhomogeneous character. But these calculations show that selection effects are a concern in general, also for dynamical black hole mass measurements.

3. Biases of broad-line AGN samples

Selection effects are a major concern for samples selected as broad line AGN. For these objects the black hole mass can be easily estimated using the virial method. This makes them so valuable for studies of the evolution in the M_•-bulge relations, as at high redshift black hole masses cannot be determined by means of direct dynamical observations. We first discuss generic issues that are also inherent in low-z AGN samples. The application to higher redshifts will be covered in Sect. 4.

Samples of non-active galaxies are typically drawn from the galaxy luminosity function. In contrast, AGN samples are effectively drawn from the luminosity function of broad-line AGN. This can induce non-trivial selection effects on the sample, as already emphasized by Lauer et al. (2007). Firstly, there will be a luminosity bias. Flux limited AGN surveys will on average contain more luminous objects then an “ideal” volume-limited sample. Thus, the AGN will on average host more massive black holes than in the volume-limited case. The intrinsic scatter in the black hole-bulge relations, together with the steep decrease of the spheroid distribution function, will then cause a Malmquist type bias towards more massive black holes at a given spheroid mass. The strength of the bias depends on the luminosity limit, the shape of the spheroid distribution function, the intrinsic scatter in the M_•-bulge relation, and the distribution of Eddington ratios. We will further discuss and illustrate these details below.

A second effect, so far not discussed in the literature, can be called an “active fraction bias”. For any sample of broad line AGN, their black holes are not randomly drawn from the entire black hole population (as described by the BH mass function), because they are by definition active. Only a minority of all black holes is in an active state, quantified by the active fraction (or equivalently by the AGN duty cycle). If this active fraction is a function of M_• then the intrinsic scatter in the M_•-M_Bulgerelation will induce an additional bias, as will be shown below. If, for example, the active fraction decreases with increasing black hole mass, then for a given spheroid mass it will be more probable to find a smaller mass black hole in an AGN sample, causing a bias towards a lower M_•/M_Bulge ratio. Conversely, for an increasing active fraction we expect a positive bias, while for a constant (M_•-independent) active fraction no such bias will occur. Thus, this bias can work in both directions, depending on the black hole mass dependence of the active fraction. It will add to the luminosity bias.

3.1. Luminosity limited samples

We now show analytically how these selection effects influence the M_•-bulge diagram. Broad line AGN are drawn from the AGN luminosity function Φ_L(l), with l ≡ log L_bol. The luminosity of an AGN is produced by mass accretion onto a SMBH, and it is therefore determined by the black hole mass and the mass accretion rate. The latter is expressed in normalised units by the Eddington ratio λ = L_bol/L_Edd. The Eddington luminosity is proportional to M_• , thus the bolometric luminosity is given by l = log λ + μ + 38.1. The AGN luminosity function Φ_L(l) (in logarithmic units) is then given by (12)By Φ_•,a(μ) we define the active BHMF, where active black holes are defined as those that contribute to the corresponding AGN luminosity function. Thus, when restricting the AGN sample to type 1 AGN (showing broad lines), the active BHMF includes only these type 1 AGN. The distribution function p_λ(l − μ) gives the probability of finding a black hole with mass μ given an AGN luminosity l, i.e. it corresponds to the normalised distribution function of Eddington ratios λ. We implicitly assume that the Eddington ratio distribution function is independent of black hole mass, i.e. p_λ(μ | l) = p_λ(l − μ) = p_λ(log λ). Type 1 AGN are observed to have Eddington ratios in the range 0.01 < L_bol/L_Edd < 1. When applying a lower threshold to the Eddington ratio, not all black holes are currently in an active state. The active fraction is defined as the ratio between the active and the total black hole population, estimated as (13)with Φ_•,a(μ) and Φ_•,q(μ) as the active and quiescent BHMFs, respectively. Therefore, the AGN luminosity function can be expressed as (14)There is thus a connection between the bulge property and the AGN luminosity, which is however smeared out by convolution with a set of additional distribution functions.

Following Eq. (5), the bivariate black hole-bulge distribution function for a luminosity limited sample is given by (15)If we assume a fixed lower luminosity limit l_lim, the selection function is (16)The bivariate distribution function is then (17)Hence, Φ_L(l) and Ψ_o(s,μ) are just different projections of a multivariate distribution function Ψ_o(s,μ,l). The non-trivial difference between Ψ(s,μ) and Ψ_o(s,μ) is the source of the bias. It is controlled by the three probability distributions, p_ac(μ), p_λ(l − μ), and g(μ | s). The distribution of p_ac(μ) regulates the active fraction bias, p_λ(l − μ) rules the luminosity bias and g(μ | s) controls the overall strength of the bias.

In Eq. (17) we can combine the parts governing the AGN  selection by defining the selection function, integrated over AGN luminosity, (18)and thus Ψ_o(s,μ) = Ω(μ)Ψ(s,μ).

The mean relations obtained from such a sample, affected by an AGN luminosity limit, are (19)(20)The relation  ⟨ μ ⟩ (s) is independent of the galaxy distribution function, but it is affected by the selection function. Thus,  ⟨ μ ⟩ (s)a + bs, as would be the case without selection effects. On the other hand  ⟨ s ⟩ (μ) is already intrinsically affected by the galaxy distribution function, but it is independent of the selection function. Therefore,  ⟨ s ⟩ (μ) contains information on the intrinsic relation, unaffected by the AGN luminosity bias and active fraction bias. This is an interesting fact, as it provides a potential route to study the evolution in the M_•-bulge relations without having to account for AGN-specific selection effects. However, even if  ⟨ s ⟩ (μ) is not biased by the selection effects, the mean offset  ⟨ Δs ⟩  is still biased, as shown by Eq. (9), because we have to integrate over the entire observed black hole mass distribution.

We now first illustrate the AGN bias in the M_•-bulge plane, and then investigate it in terms of a simple offset from the input relation. For this purpose, we need to know all the underlying distribution functions, for which we here adopt the z ≈ 0 values. For Φ_s(s) we again use the above estimated spheroid mass function. For g(μ | s) we assume a log-normal distribution with intrinsic scatter σ = 0.3 dex around the relation from Häring & Rix (2004). These then define the quiescent BHMF via Eq. (4). If we further know the active BHMF, the active fraction is given by Eq. (13). The active BHMF and the Eddington ratio distribution function (ERDF) p_λ(l − μ) for type 1 AGN in the local universe were determined by Schulze & Wisotzki (2010, hereafter SW10), and we use those values, adopting a modified Schechter function for the BHMF and a Schechter function for the ERDF. The SW10 active BHMF implies a significant decrease of the active fraction with increasing black hole mass, thus the active fraction bias will work towards a lower M_•/M_Bulge ratio, opposite to the luminosity bias. Note that this parameterisation by construction ensures the consistency with the broad line AGN luminosity function via Eq. (14).

In Fig. 3 we show the M_•-M_Bulge diagram for a luminosity limited local AGN sample, for two lower luminosity limits, l_min = 45 and l_min = 47. At z ~ 0 the latter case is rather unrealistic, as QSOs of these luminosities are extremely rare. However, as further discussed below, it is of importance for higher z observations, where such luminosities are quite common in QSO samples. A luminosity limit of l_min = 47 roughly corresponds to the SDSS magnitude limit of i ≃ 19 at z ≃ 3 and to the SDSS high-z limit of i = 20.2 at z ≃ 6.

For both luminosity limits the M_•-M_Bulge relation is biased towards larger black hole masses at a given spheroid mass. This is already true for l_min = 45, but it is greatly enhanced for the high luminosity limit. The solid line again shows the mean black hole mass for a given spheroid mass (computed via Eq. (6)), i.e. the apparent M_•-M_Bulge relation obtained for the given sample. It is clearly visible that in general the bias is stronger for lower masses. Specifically, it is strongest close to the luminosity limit of the survey, which on average corresponds to the lowest mass black holes. The red dashed-dotted line again shows the mean spheroid mass for a given black hole mass (computed via Eq. (7)). As noted, this relation is unaffected by the AGN selection and is identical to the one shown in Fig. 1. However, it is not identical to the intrinsic relation, as it turns upwards at the high mass end due to the decrease in the spheroid mass function. Therefore, an integration over  ⟨ s ⟩ (μ) with a stronger weight on the high mass end will lead to a bias for the mean offset from the local relation, i.e. for the total sample bias.

The sample bias is greatly enhanced for a brighter luminosity limit. The steep decrease of the spheroid mass function and the intrinsic scatter in the relation itself induce an increase in the bias compared to the flat section of the spheroid mass function. This is shown in Fig. 4, where the mean offset from the input relation for the whole AGN sample,  ⟨ Δμ ⟩  (computed via Eq. (8)), is plotted for a range of bolometric AGN luminosity limits, assuming the local distribution functions. Here we assumed our reference value for the intrinsic scatter, σ = 0.3 dex, but we also indicate the dependence of the bias on the intrinsic scatter. A smaller intrinsic scatter also leads to a smaller sample bias, while a larger scatter causes a larger bias. By the dotted line we additionally show the bias predicted by Lauer et al. (2007), using their Eq. (25) and the local bolometric type 1 AGN luminosity function from Schulze et al. (2009). The main difference between their result and ours is that we not only use the AGN luminosity function, but evaluate the whole set of underlying distribution functions. In particular, we also take into account the active fraction bias that reduces the total bias, due to the decrease of the active fraction towards higher black hole masses. As a caveat we note that for very high luminosities, our knowledge of the underlying distributions (active and passive BHMFs, AGN luminosity function and M_•-M_Bulgerelation) are observationally poorly determined.

We provide an additional illustration of these selection effects by Monte Carlo Simulations. We generated a large galaxy sample drawn from the spheroid mass function, and attributed black hole masses drawn from a log-normal distribution with mean from the Häring & Rix (2004) relation. We decided if the black hole is in an active state based on the probability p_ac(μ). If it is active, an Eddington ratio was drawn from the SW10 Eddington ratio distribution function, which also sets the bolometric luminosity of the AGN. Now a number of specific selection criteria were applied to this sample, with the results presented in Fig. 5. It is notable that for a bright luminosity limit most, if not all, AGN lie above the input M_•-M_Bulge relation, while they are intrinsically drawn from it. The applied luminosity limit effectively corresponds to an unsharp lower limit on black hole mass. But at the same time also many massive black holes are excluded that accrete only at low rates. For low spheroid masses only positive outliers from the relation can be detected, while for increasing spheroid masses also black holes on and below the relation can still be detected. This causes the bias at low M_Bulge.

In a realistic AGN sample there may be additional selection effects influencing the distribution of M_Bulge. For example, usually a minimum host-to-nucleus ratio is required for the AGN host galaxy to be detected. This will correspond to some sort of cutoff in M_Bulge, excluding low mass hosts. This additional selection effect will decrease the average bias for the sample, as the excluded objects are the most heavily biased. However, this will only be of importance if objects are actively eliminated during the sample construction based on these considerations. Therefore, it is important to understand and quantify the selection effects as precisely as possible for a realistic prediction of the average bias of the given sample.

3.2. Flux limited samples

So far, we have only discussed purely luminosity limited samples. This approximation is valid for samples spanning small ranges in redshift. However, more generally we have to consider samples defined by limitations in apparent flux. (Note that luminosity limits may be present in addition, because the AGN need to be identified as such and separated from their host galaxies.) We first define the multivariate distribution function Ψ_o(s,μ,l,z), as (21)The bivariate distribution function for a flux limited sample is then (22)where ΔV_c is the comoving volume within the redshift range.

At this place we ignore possible redshift dependences of the individual distribution functions, postponing their discussion to Sect. 4.2. We also assume that the selection function is purely defined by an AGN luminosity limit and an AGN flux limit, i.e. we neglect any potential dependence of the selection function on the bulge property s. The selection function Ω is then given by (23)where f is the bolometric flux in logarithmic units.

At high redshifts and for a narrow z range, this is almost identical to the luminosity limited case. However, it describes the more general case for realistic observations and remains valid when the sample covers a wider z range, including the low z regime. The application of this form to low redshift AGN samples is illustrated in the next subsection. Consequences for higher z are discussed in Sect. 4.

3.3. Application to the reverberation mapping AGN sample

We have demonstrated above that already the low z AGN population will be affected by the details of the sample selection. The main difference compared to high z is that for AGN in the local universe the luminosity limit is generally low and therefore the bias is less severe. Nevertheless, it is present, and as a test case we briefly estimate its influence on the AGN sample with black hole masses determined via reverberation mapping.

The technique of reverberation mapping (e.g. Blandford & McKee 1982; Peterson 1993) provides the most precise black hole mass estimates for type 1 AGN. It builds the foundation for the virial method through establishing a scaling relation between the size of the broad line region and the continuum luminosity (Kaspi et al. 2000, 2005; Bentz et al. 2009a). However, this method currently only provides measurements of M_• up to a scale factor depending on the (largely unknown) geometry and dynamics of the broad line region. The usual approach to determine this scale factor and thus fix the virial mass scale is to scale the reverberation mapped black hole masses of galaxies with stellar velocity dispersion measurements to the local M_•-σ_∗relation of quiescent galaxies (Onken et al. 2004; Woo et al. 2010). However, this implicitly assumes that the reverberation mapping (RM) sample follows the same relation as quiescent galaxies. If its M_•-σ_∗relation is biased by selection effects then this bias will propagate into the absolute calibration of the virial mass estimates.

Properly defining the selection function for the RM sample is difficult, due to the heterogeneous selection of the objects. The RM sample as such is certainly not a well defined sample. Nevertheless, in order to estimate potential systematics that may be inherent in the sample, we tried to estimate the selection function, as a kind of first order approximation. The full sample covers a wide range in luminosity, from bright quasars to moderate-luminosity Seyfert galaxies. Especially for the Seyfert galaxies there is no clear AGN luminosity dependence for the selection. However, there must be an implicit limit on the AGN luminosity and also on the AGN flux. First of all, the AGN has to be luminous enough to classify the galaxy as harbouring an active nucleus. Furthermore, the continuum and broad lines have to be sufficiently bright to enable spectroscopy and the measurement of a reliable time lag. For simplicity, we relate the low luminosity limit to the faintest AGN in the sample of Bentz et al. (2009a), NGC 4051, with log λL_λ,AGN(5100Å) = 41.9. This approximately corresponds to log L_bol = 42.8. Equally, we use the object with the lowest flux, PG 1229 + 204, to define the flux limit. This yields a bolometric flux limit of f_min =  −10.4 in logarithmic units, which is also in good agreement with the magnitude limit of the Palomar-Green (PG) quasars in the sample (Schmidt & Green 1983; Kaspi et al. 2000). We adopt a redshift range of 0.003 ≤ z ≤ 0.15, which contains the vast majority of the RM AGN. We also tested the effect of extending the range to z ≤ 0.3, including the full RM sample.

If we assume only a luminosity limit, the total sample bias  ⟨ Δμ ⟩  predicted by our model would be negligible (see Fig. 4). Indeed, when omitting the flux limit we find  ⟨ Δμ ⟩  = 0.003, as a lower limit to the bias. However, the luminosity limit is only important for z0.008. At higher z the flux limit dominates the selection function. Incorporating the flux limit we get  ⟨ Δμ ⟩  = 0.09 for z ≤ 0.15 ( ⟨ Δμ ⟩  = 0.12 for z ≤ 0.3). If we relax our conservative flux limit, the expected sample bias would decrease slightly. Decreasing the assumed luminosity limit has almost no effect, as we are dominated by the flux limit. Increasing the luminosity limit will slightly increase the bias.

In Fig. 6 we show the predicted M_•-M_Bulgediagram. At the high mass end no bias is expected, whereas a mild bias at low masses is present. Thus, the M_•-M_Bulgerelation is slightly affected.

We further tested this result with Monte Carlo simulations. We restricted our Monte Carlo sample to the same luminosity limit and flux limit and fitted the culled sample with a maximum likelihood method (see e.g. Gültekin et al. 2009), with slope b, normalisation a and intrinsic scatter σ as free parameters. From the restricted sample, we fitted 1000 random subsamples of 100 objects each. The distribution of the free parameters are shown in Fig. 7. We found mean values of a = 8.31, b = 1.12 and σ = 0.29, compared to input values of a = 8.2, b = 1.12 and σ = 0.3, i.e. we recover the predicted mean offset and confirm that the slope and scatter are not strongly affected.

Therefore, based on our model assumptions, we estimate that the reverberation mapping sample may be biased towards a high M_•/M_Bulge ratio by  ~ 0.1 dex. This would correspond to an underestimation of the virial method by the same amount, when using the scale factor normalised to the quiescent M_•-σ _∗ relation (Onken et al. 2004; Woo et al. 2010). However, due to the poorly defined selection function of the RM sample this can only be taken as a qualitative evaluation of possible selection effects.

3.4. On a black hole mass uncertainty bias

Shen & Kelly (2010) discussed an additional bias that will affect AGN samples with virial mass estimates. This is a Malmquist type bias that is expected to arise from the steep turn-down of the active BHMF and the intrinsic uncertainty of virial BH mass estimates. They argue that this is independent of the luminosity bias discussed by Lauer et al. (2007) and thus adds an additional bias to the observations. In our framework we are directly sampling from the active BHMF, thus we are already taking this effect into account. If we introduce an uncertainty in the mass estimates this will only increase the scatter, but it will not affect the mean relation  ⟨ μ ⟩ (s) or  ⟨ Δμ ⟩ . This can again be illustrated by the Monte Carlo simulations. Adding a scatter in the BH masses does not affect our sample selection, as it is by construction based on luminosities, which can be measured with relatively high precision. Thus the uncertainty just moves the black hole masses symmetrically towards lower or higher values, without affecting the mean.

It is also straightforward to see this from the distribution function Ψ_o(s,μ). We assume that the virial mass estimate (the observed mass μ_o) is given by a log-normal probability distribution with mean μ and dispersion σ_vm, the uncertainty in the virial mass estimate, (24)The bivariate distribution function for bulge property and virial black hole mass is then (25)The mean relation is (26)when integrating over μ_o. This is identical to Eq. (6), so no additional bias is introduced by the virial mass uncertainty.

The conceptual difference between our results and the work by Shen & Kelly (2010) is, firstly, that we do not treat the black hole mass bias independently, but within the total bias budget. Secondly, we account for the fact that in observations we do not sample directly from the true BHMF, but only from the apparent BHMF affected by the relevant luminosity limits. As discussed by SW10 in detail, this directly observable distribution suffers from incompleteness at the low mass end and therefore turns over towards low M_•.

In Fig. 8 we show the result for a luminosity limit of l_min = 47, assuming σ_vm = 0.3 dex. While the black hole mass distribution is broadened, the mean relation is unchanged.

3.5. On a spheroid uncertainty bias

At least at high redshifts, not only the estimated black hole masses are significantly uncertain, but also the observed spheroid properties. For example, the emission line width of [O iii] has been used several times as surrogate of σ _∗ . However, this is at best on average a reliable estimator, with a dispersion of at least  ~0.2 dex (Nelson 2000). Furthermore, spheroid masses estimated for AGN host galaxies suffer from several uncertainties, including the difficulties of performing a proper AGN-host decomposition, of deblending bulge and disk components, and of converting host galaxy luminosities into stellar masses, using either a fixed mass-to-light ratio or colour information. These uncertainties may well accumulate to a total error of similar magnitude as for the black hole masses, or more. In that case, again a Malmquist type bias is possible.

We model the observed spheroid property s_o by a log-normal probability distribution with mean s and dispersion σ_s, analog to Eq. (24). The bivariate distribution function is then (27)and the mean relation is (28)The total sample bias  ⟨ Δμ ⟩  is not affected by the uncertainty, as can easily be verified. However, the predicted M_•-bulge relation is affected, as shown in Fig. 9. At the low mass end the mean M_• increases compared to the prediction without scatter, while at the high mass end it decreases. To understand this trend we have to recall that we are sampling from the luminosity-limited AGN host galaxy spheroid mass function. Due to the luminosity limit the density of sampled objects decreases at the low mass end, similarly as in the active BHMF (cf. SW10). Because of this decrease at the low mass end more objects are scattered from slightly larger masses to lower masses than the other way around. This produces an excess of higher true s at a given observed s_o, which in turn leads to an increase in the average μ at the given s_o. At the high mass end the reverse happens. Due to the steep decrease towards higher masses, more black holes are scattered to higher s_o, having on average a lower μ, as they intrinsically have a lower s than observed. In total, the observed M_•-bulge relation flattens, without changing its normalisation.

For the mean spheroid mass at a given black hole mass the opposite effect happens. Uncertainties in the spheroid masses do not change the mean relation  ⟨ s_o ⟩ (μ), but uncertainties in the black hole masses can strongly affect the mean relation  ⟨ s ⟩ (μ_o), causing a steepening of the observed relation. This is indicated by the red lines in Figs. 8 and 9.

Finally, we note that in general the convolution with the uncertainty should really be done with Ψ(s,μ) rather than with Ψ_o(s,μ), i.e. before applying the selection criteria. In our case of a purely luminosity limited sample, the selection function does not depend on s_o or μ_o. Therefore, our approach is justified and no additional bias is introduced for the sample in total. In case of a more complicated selection function that depends on s_o or μ_o this might no longer be the case. An additional bias by the measurement uncertainty then becomes possible which will depend on the details of the selection function.

3.6. The slope in low-redshift AGN samples

The black hole-bulge relations determined for low-redshift AGN samples usually plays an important role in the search for possible evolution in the relations. But they are also interesting in themselves, as it is not at all clear from a physical perspective that AGN host galaxies should follow exactly the same relation as largely inactive galaxies. Notice that the current calibration of the virial method ensures by design that AGN and non-AGN relations have no global offset from each other. But it has been noted that the slope of the M_•-bulge relation for some local AGN samples appeared to be flatter than in the relation of black holes with dynamical black hole mass measurements (Greene & Ho 2006; Kim et al. 2008; Bentz et al. 2009b; Woo et al. 2010). We were interested to see whether this flattening can be explained within our framework by selection effects, or whether an intrinsic difference in the relations between active and inactive galaxies has to be invoked.

In Fig. 10 we show the AGN samples for which such evidence for a flatter slope of the M_•-bulge relations has been diagnosed. We compare them with predictions of the bivariate distribution function Ψ_o(s,μ), adopting plausible numbers for the selection function. For the M_•-σ_∗relation we used the stellar velocity dispersion function from Sheth et al. (2003) and assume our approximation of the selection function for the reverberation mapping sample from Sect. 3.3. For the M_•-L_V relation we adopted the galaxy luminosity function from Kochanek et al. (2001), converted to the V-band and into a spheroid luminosity function following Marconi et al. (2004). To roughly model the selection function we assume an AGN luminosity limit of log L_bol = 10⁴⁵ erg s^-1.

Adopting these selection functions, we show in Fig. 10 the predicted mean relations  ⟨ μ ⟩ (s) as thick solid lines. The bias caused by the AGN selection depends strongly on s and is largest at low s. Thus the relation does not only change in offset but also flattens compared to the true relation. Furthermore, measurement uncertainties in the bulge property lead to a further flattening of the slope, as discussed in 3.5. Indeed, for active galaxies the bulge property is usually measured with larger uncertainty than for inactive galaxies. We computed the expected mean relations  ⟨ μ ⟩ (s) adding reasonable values for the measurement errors in the respective spheroid properties to the bivariate distributions. These are shown in Fig. 10 as dotted black lines.

Considering our rough approximations to the respective sample selection functions we find a very good agreement between our prediction of a flatter M_•-bulge relation and the observations. Thus we argue that these selection effects are able to reproduce the observed flattening in the slope. There is no reason to assume that AGN have generally a different black hole-bulge relation than inactive galaxies.

We note in passing that the M_•-L_V relation seems to disagree with our bivariate probability distribution at the high mass end. This may reveal a problem with our assumption, but it may also just indicate that the M_•-L_V relation is not the same for inactive galaxies and AGN host galaxies. In fact that is even expected to be so, since it is known that host galaxies of high luminosity AGN have on average younger stellar populations (e.g. Kauffmann et al. 2003; Jahnke et al. 2004; Vanden Berk et al. 2006) compared to inactive galaxies. Thus, if they obey the same M_•-M_bulge relation, they actually should deviate in the M_•-L_Bulge relation.

4. Evolution in the M_•-M_Bulge relation

4.1. Evolution in a flux limited sample

How will selection effects influence studies that test for redshift evolution in the black hole-bulge relations? If exactly the same sample selection criteria are applied at high z as for a z ≈ 0 comparison sample, no bias will be present, at least to first order (see Sect. 4.2). Any change in the M_•-bulge relations could then be taken as evidence for their evolution. However, this simple situation is usually not given. The local comparison is commonly provided either directly by the quiescent relation, or by a local type 1 AGN comparison sample. While the former clearly possesses different selection criteria, also the latter is not automatically selected in the same manner. This has to be taken into account when comparing low and high redshift results.

To illustrate this point we derive the sample bias for a flux limited sample at a given redshift, (29)with the multivariate distribution function Ψ_o(s,μ,l,z) given by Eq. (21). We thus assume these distributions, for the moment, to be non-evolving with redshift (we discuss the more general case below). Equation (29) is equivalent to the bias for a luminosity limited sample, but with a redshift dependent luminosity limit.

Instead of using a bolometric flux limit as before, we now adopt an optical flux limit, for example in the B band. For simplicity we assume a simple power law K-correction with spectral index α =  −0.44 (Vanden Berk et al. 2001). The B band luminosity is converted to bolometric luminosity using the bolometric corrections of Marconi et al. (2004). In Fig. 11 we provide the expected sample bias at a given z for different apparent magnitude limits. The deeper the respective survey, the lower is the expected bias. For a fixed flux limit, the expected bias increases with redshift, which can mimic an evolutionary trend.

Furthermore, to have a comparable bias between a local sample with, say, B_min ≈ 19 mag and a sample at z ≈ 2, the high redshift sample has to be complete down to B_min ≈ 24 mag. The use of a local AGN sample as comparison sample for higher z studies is already an improvement over the use of the quiescent black hole sample, but it does not ensure the absence of a selection bias. Ideally, the local comparison sample should be matched in AGN luminosity to the high z sample.

On the other hand, if the selection function of the respective sample at high z is known, the expected bias can be computed, and any clear offset from this prediction could be interpreted as evidence for evolution in the M_•-bulge relations. However, we show in the following subsection that unfortunately even this approach is not fully unbiased.

4.2. AGN evolution biases

So far we ignored any possible explicit redshift dependence of the underlying distribution functions in Eq. (21). However, at least some of them must evolve with redshift. The stellar mass function, and thus also the spheroid mass function, are certainly changing with z (e.g. Bundy et al. 2005; Franceschini et al. 2006; Pozzetti et al. 2007; Ilbert et al. 2010). Furthermore, it is well known that the AGN population itself strongly evolves, apparent in the evolution of the AGN luminosity function (e.g. Ueda et al. 2003; Hasinger et al. 2005; Richards et al. 2006; Bongiorno et al. 2007; Croom et al. 2009). Not only the normalisation and typical luminosity change with z but also the shape of the luminosity function. At low z the faint end of the QSO luminosity function steepens, known as “AGN cosmic downsizing”. While in the local universe the QSO luminosity function shows only mild curvature (Schulze et al. 2009), a prominent break is present in the high z luminosity function. This directly implies evolution in the active black hole mass function, the Eddington ratio distribution function, or most probably in both, leading to a change of the predicted bias with redshift, even for a fixed luminosity limit. Practically speaking it is generally not straightforward to judge if an observed evolution in the M_•-bulge relations is caused by evolution in the intrinsic relations, or whether the trend is spurious and caused by evolution in the underlying distribution functions.

A prediction of the sample bias at given redshift requires, apart from a well defined selection function, also knowledge about these underlying distribution functions, the spheroid mass function, active BHMF and ERDF. In the local universe these distributions are at least reasonably well established (SW10), but for the high redshift universe the situation is currently less comfortable. The distribution function of velocity dispersions is reasonably well established for the local universe (Sheth et al. 2003), but rather uncertain beyond that (Chae 2010; Bezanson et al. 2011). The total stellar mass function is observed up to z ≈ 4 (Fontana et al. 2006), but galaxy mass functions for different morphological types have been determined only up to z ≈ 1.4 (Bundy et al. 2005; Franceschini et al. 2006; Ilbert et al. 2010). The spheroid mass function itself is again basically unknown for higher z. To enable an analysis of how evolution influences the biases, we now derive very rough estimates of these distribution functions for higher redshifts. These evolving distributions are meant as illustrative only; a more detailed investigation would be beyond the scope of the current paper.

An upper limit to the spheroid mass function is given by the total stellar mass function. We use the parametric fit to the stellar mass function by Fontana et al. (2006) for this purpose. A lower limit is given by the mass function of elliptical galaxies. To derive this distribution for arbitrary z we assume the same elliptical-to-total ratio as for the local mass function from Bell et al. (2003), and apply this correction to the total mass function at higher redshift. This is clearly an oversimplification, as the relative contribution of elliptical galaxies seems to decrease with increasing z, at least to z ~ 1 (e.g. Bundy et al. 2005). However, at higher redshifts the low mass end that is most affected by this correction is not well determined and may be underestimated if also the distribution of mass-to-light ratios changes. For the present purpose we hold that our simple approximation serves as a reasonable lower limit of the spheroid mass function. Additionally, we define a case in between these limits, using the same conversion of total mass function to spheroid mass function as assumed above for redshift zero (Sect. 2.1).

There has been significant progress in the determination of the active BHMF and Eddington ratio distribution function at high redshifts in the last years (Vestergaard et al. 2008; Vestergaard & Osmer 2009; Kelly et al. 2010; Shen & Kelly 2011). But these results mainly cover the bright end of the luminosity function and thus the high mass end of the BHMF, while the low mass end is still poorly determined; also the systematics are not fully understood. Here we restrict ourselves to the use of a mass function that is consistent with current observations.

We now used the local active BHMF and ERDF as well as the redshift evolution of the type 1 AGN luminosity function (LF) as constraints to predict the distribution functions at higher redshift. To achieve a smooth redshift evolution of the BHMF and ERDF we explored several arbitrary, but reasonable evolution models for these distribution functions. We fixed the zero point of the two distribution functions to their local values (as given by SW10), and fitted their redshift evolution to the optical LF (employing Eq. (12)). For the LF we used the luminosity dependent density evolution model from Bongiorno et al. (2007), which is based on the SDSS at the bright end and on the VVDS at the faint end. Optical B band magnitudes were converted to bolometric luminosity using the bolometric corrections from Marconi et al. (2004). Alternatively we also explored modelling the BHMF and ERDF directly from the observed LF in individual redshift bins, which lead to the same qualitative results.

As long as only the LF is used as constraint there is clearly a degeneracy between evolution in the BHMF and in the ERDF, respectively, while presumably both are evolving. We tried to bracket this degeneracy by exploring two extreme cases. Firstly, we assumed a constant ERDF and let the BHMF change with z. Secondly, we fixed the active BHMF and let the ERDF evolve with z.

The first case, a non-evolving ERDF, can be regarded as an upper limit to the expected AGN evolution bias. The downsizing of the AGN LF directly corresponds to a downsizing in the active BHMF. In Fig. 12 we show the best fit BHMFs, ERDFs, and the corresponding AGN LFs, for a few representative redshifts. Here we adopted a mass-dependent density evolution model to model the redshift dependence of the active BHMF, inspired by the LDDE model used by Bongiorno et al. (2007) to parameterise the evolution in the AGN LF. In order to reproduce the observed AGN LF, the space density at the high-mass end has to increase substantially. Thus, the mass dependence of the active fraction is reduced, which in turn changes the active fraction bias. Because this effect tends to counteract the AGN luminosity bias, a reduced active fraction bias leads to an increase of the expected overall sample bias compared to the redshift zero case. This is illustrated in Fig. 13, where we show the expected sample bias  ⟨ Δμ ⟩  as a function of the applied lower luminosity limit for various redshifts.

The second case, a fixed BHMF and an evolving ERDF, provides a lower limit to the expected AGN evolution bias. In this scenario the active fraction is only mildly evolving with z, through the evolution of the spheroid mass function. On the other hand, the ERDF needs to evolve strongly, in the sense of an increasing average log λ, to satisfy the constraints from the AGN LF. For a higher average log λ the same luminosity limit corresponds to a lower average black hole mass limit, thus reducing the expected bias. Therefore, this assumption will predict the lowest sample bias. However, the best fit model enforces an unreasonably high space density of objects accreting at super Eddington rates and also provides a poorer fit to the AGN LF. Furthermore, it disagrees with the picture of anti-hierarchical black hole growth (Marconi et al. 2004; Merloni 2004; Merloni & Heinz 2008). Thus, this case is not physically plausible. We also experimented with intermediate evolution scenarios, where a mild evolution of the ERDF balances the increasing bias through the evolution of the BHMF to some degree. For all physically plausible scenarios, we found at least some increase of the bias compared to the local case. We conclude that the estimated sample bias for the local universe is a lower limit for high z samples.

In Fig. 14 we show the expected bias at a fixed magnitude limit as a function of redshift, equivalent to Fig. 11, but now allowing for an evolving active BHMF and a constant ERDF (our case 1 above). As discussed, the resulting bias can be taken as an upper limit. Within the assumptions made in our model, only observations that show a clear excess on top of this predicted bias would constitute firm evidence for real evolution in the M_•-bulge relations.

While evolution is usually expressed in terms of the normalisation of the relation, it is also possible that there is evolution in its intrinsic scatter (Merloni et al. 2010). The strength of the sample bias is strongly affected by the amount of intrinsic scatter in the M_•-bulge relations. Therefore, an offset of the M_•/M_bulge ratio on top of the expected selection bias can equally be interpreted as being caused by an increased bias due to an increased intrinsic scatter, rather than as a true offset in the zero point of the respective M_•-bulge relation. To illustrate this effect, we arbitrarily assume a redshift evolution for the intrinsic scatter as σ(z) = σ_z = 0(1 + z)^0.5 with σ_z = 0 = 0.3, giving approximately σ(z = 2) ≈ 0.5. In Fig. 15 we show the resulting redshift evolution of the bias for fixed magnitude limits. The predicted scatter is strongly enhanced in this case, especially at high redshifts.

An increase in the intrinsic scatter is also some kind of evolution in the relations. However, its interpretation is of course different from a change in the zero point. Whereas the latter would suggest a phase where black holes grow stronger than galaxies, or the other way around, the former would be consistent with a true coeval growth on average, in which the correlation is tightened from a rather loose one to the tight correlation we observe today. Indeed, such an evolution in scatter towards a tighter relation is expected in various models (e.g. Peng 2007; Volonteri & Natarajan 2009; Lamastra et al. 2010). Given sufficiently high quality data and large enough samples, it should be possible in the future to estimate the amount of evolution of the intrinsic scatter in the relation.

5. Discussion

Several observational studies of type 1 AGN samples found tentative evidence for an increase in the M_•/M_Bulge ratio with redshift. This increase is often parametrised as an evolutionary behaviour M_•/M_Bulge ∝ (1 + z)^γ, with values for γ up to 2.1 reported (McLure et al. 2006). As discussed in this paper, these observational studies are invariably affected by selection effects. We now apply our formalism to a few published studies. We estimate approximate distribution functions and predict the resulting sample biases. Note that we do not aim at deriving detailed corrections, or at presenting an exhaustive discussion of the relevant literature. Rather we wish to highlight the magnitudes of possible systematic biases.

Specifically, we explored two scenarios for the distribution functions that serve as upper and lower limits for the magnitude of our predicted selection bias.

• Model 1: We used the local distributionfunctions presented in Sect. 3.1 throughout theentire redshift range, i.e. we ignore any effect of redshiftevolution. As discussed above, this will serve as a lower limit to theexpected bias even if the scenario is unrealistic.

• Model 2: Here we incorporated redshift evolution in the underlying distribution functions, in particular in the AGN population. We used the model discussed in Sect. 4.2, assuming a mass dependent density evolution for the BHMF and a non-evolving ERDF. This model provides an approximate upper limit to the selection bias.

5.1. Merloni et al. (2010)

A relatively large and well-defined sample was investigated by Merloni et al. (2010). They used an I band limited set of objects drawn from zCOSMOS (Lilly et al. 2007) in the redshift range z =  [1.06,2.19] , with I_AB < 22.5, for which they estimated stellar masses by fitting single stellar population models to the observed multiband SEDs. The results are shown in Fig. 16. While the individual data points show considerable scatter, they are on average offset from the z = 0 relation by  ⟨ Δμ ⟩ _obs = 0.34. Merloni et al. (2010) considered also selection biases, largely following the argument given by Lauer et al. (2007), as one possible origin of this offset. They concluded that their observations indicated either evolution in M_• /M_Bulge  or a significant increase of intrinsic scatter on the M_•-M_Bulgerelation.

We modelled their sample selection and applied the two above mentioned models for the distribution functions to predict the expected sample bias. The correspondingly modified bivariate probability distributions are overplotted in the left and middle panels of Fig. 16. Neither reproduces the distribution of data points very well, but it seems that an intermediate model might actually do a reasonable job. In the right-hand panel of the figure we compare the expected average sample offset due to selection biases, as a function of redshift, with the data points (individual and also binned into three redshift ranges). The binned mean offsets agree quite well with the range of values covered by our bias model. The prediction for a M_• /M_Bulge  ratio evolving with redshift as γ = 1, on the other hand, lies slightly above the binned data.

We pointed out in Sect. 3.1 that the relation of mean stellar masses at given black hole mass is to first order unaffected by the AGN selection bias and depends mainly on the spheroid mass function. We binned the data points in M_•  accordingly and compare this binned relation in Fig. 17 with the prediction for the null hypothesis (i.e., consistency with the local relation). We also computed the same predicted relations for an evolving M_• /M_Bulge  ratio (again with γ = 1), and also for a higher intrinsic scatter in M_• /M_Bulge  of 0.5 dex (but no evolution). Both curves are offset from the null hypothesis towards higher BH masses, albeit at slightly different tilts. We performed a simple χ² test to compare the binned data with the predicted distributions. The lowest value of χ², and in fact the only one close to a reduced χ² of unity, was obtained for the null hypothesis. We thus conclude that this sample provides no statistically significant evidence for evolution in M_• /M_Bulge  or its scatter, and that the apparent offsets found by Merloni et al. (2010) are consistent with the expected magnitude of AGN sample selection effects.

5.2. Jahnke et al. (2009)

Jahnke et al. (2009) studied the M_•-M _∗  relation for a sample of 10 AGN from COSMOS (Hasinger et al. 2007), using optical and NIR HST images to estimate stellar masses of the host galaxies. Their sample is X-ray selected, but the spectroscopic follow-up requirement leads to a similar optical flux limit as for the work by Merloni et al. (2010; see Trump et al. 2009). As the covered redshift range is very similar as well, we simply added the results by Jahnke et al. (2009) to Fig. 16. Jahnke et al. (2009) found their sample offset Δμ to be close to zero and concluded that their data are fully consistent with the local relation. Figure 16 suggests that when accounting for selection effects, these data actually might even imply a negative evolution. It is interesting to note that also the small X-ray selected sample by Sarria et al. (2010), covering a similar redshift range finds no, or even a negative evolution. However, notice that both studies (and also Merloni et al. 2010) estimated only total galaxy masses instead of bulge masses. Jahnke et al. (2009) discussed possible consequences of this important difference for the true evolution of the M_•-M_Bulgerelation.

5.3. Salviander et al. (2007)

We now consider the study by Salviander et al. (2007), who studied the M_•-σ _∗  relation in the redshift range 0 < z < 1.2, using the [O iii] line width as surrogate of σ _∗  for z < 0.8, and the [O ii] line width for 0.4 < z < 1.2. Their sample is drawn from the Sloan Digital Sky Survey Data Release 3 (SDSS DR3; Abazajian et al. 2005), which does not constitute a well-defined sample (see e.g. Richards et al. 2006). Furthermore, they excluded a large fraction of objects based on quality cuts. To model their sample we simply adopted the flux limit of the SDSS main quasar sample, i < 19.1 (Richards et al. 2002). In Fig. 18 we show our predicted redshift evolution for the sample bias and compare them to the observations (which we augmented by similar low-z measurements by Shen et al. 2008). For z < 0.6, their data points are slightly below, for z > 0.6 slightly above our prediction for Δμ(z) under the null hypothesis of no evolution in the M_•-σ _∗  relation. We also computed Δμ(z) again for the case of evolution with γ = 1.

In their discussion of possible biases in their results, Salviander et al. (2007) attributed  ~ 0.25 dex of their offset to selection effects, which roughly agrees with our estimate at z ~ 0.7. They also discussed another bias caused by their signal-to-noise requirement for the narrow lines, which may accumulate to an additional  ~ 0.15 dex at high z. Taking this additional effect into account, the data appear again to be largely consistent with the null hypothesis of no evolution.

5.4. High redshift luminous QSOs

The largest apparent offsets from the local M_•-bulge relations were so far observed for very luminous QSOs at high redshifts, z3 up to z ~ 6 (e.g. Walter et al. 2004; Wang et al. 2010; Targett et al. 2011). Because of the limitation to very luminous objects, significant selection effects are expected. At these high redshifts our knowledge of the required distribution functions becomes extremely vague, and we can do no more than guess the magnitude of the sample biases. It is nevertheless instructive to see how large they can become. We illustrate this by using the small sample of three hyperluminous QSOs around z ~ 3 from Schramm et al. (2008). These quasars originate in the Hamburg/ESO Survey (Wisotzki et al. 2000) which has an approximate magnitude limit of B17.5. Stellar masses were estimated from NIR imaging that provided rest-frame B − V colours.

Schramm et al. (2008) obtained M_• /M_∗ values that were an order of magnitude larger than found in the local universe. In Fig. 19 we superimpose this small sample over the predicted bivariate probability distribution assuming no evolution. The data are well captured by this distribution; in fact, it is expected that luminous QSOs would nearly always be found well above the local relation. Thus, there is again no evidence even at z ~ 3 that the M_•-M_Bulgerelation is any different from today.

However, the error bars of these and similar measurements are considerable, and some level of evolution would clearly be equally consistent with the data. To see what the expected signature for very strong evolution might be we show in Fig. 19 also the probability distribution for evolution with γ = 2, a value that has been suggested by some authors in the recent past. In that case, offsets of almost 2 order of magnitude would be expected; the data by Schramm et al. (2008) are clearly inconsistent with such a scenario.

The general features discussed with this example are applicable to other studies of bright QSOs at high z. In particular, the z ~ 6 SDSS main quasar sample (Fan et al. 2001) has a similar luminosity limit, and sample biases at least as strong can be expected. However, we make no attempt to guess the relevant AGN distribution functions at z ~ 6 and thus have to refrain from making a speculative statement about possible biases for these highest redshift objects.

6. Conclusions

We investigated the ramifications of sample selection effects for the observed black hole-bulge relations. Our starting point is the bivariate probability distribution of galaxy bulge properties and black hole masses, which ideally describes the true underlying relation. However, this relation is recoverable only from a strictly volume-limited sample with completely measured galaxy properties and black hole masses. In real observations the probability distribution is inevitably modified. We incorporate the modification of the bivariate distribution by the use of an appropriate selection function which represents the observational process of sample definition and measurement. In the presence of selection effects and scatter, the true M_•-bulge relations are generally not recovered by the observations.

Others have recognised the importance of selection effects before, and sometimes tried to include them in their budget (Salviander et al. 2007; Merloni et al. 2010; Bennert et al. 2010; Lamastra et al. 2010). However, different papers focussed on different effects, and sometimes with quite simplifying assumptions. In this paper we go one step further and provide a common framework in which all kinds of selection effects on the M_•-bulge relations can be investigated. We also identified some important sources of bias that were not discussed in the literature so far.

Already the widely used M_•-bulge relations derived from observations of nearby galaxies with dynamical black hole mass measurements do not strictly reflect the underlying intrinsic relations. We showed how the bivariate probability distribution is changed due to explicit or implicit selection against galaxies where the central black hole sphere of influence is not well resolved. The resulting bias can be quite strong, almost of the order of the intrinsic scatter, when a high threshold is adopted. It becomes weaker when the threshold is lowered and more galaxies are included, but this comes at the penalty of increasing the observational errors. Measuring BH masses in galaxies with poorly resolved spheres of influence also puts higher demands on the dynamical modelling, such as the inclusion of Dark Matter halos.

For the main part of this paper we investigated in detail a number of selection effects for type 1 AGN samples, as such datasets are essential as probes for redshift evolution in the M_•-bulge relations. We identified a variety of effects:

• Active fraction bias: The probability for a black hole to be in anactive (type 1 AGN) state depends on the blackhole mass. At least in the local universe, very high mass SMBHsare much more commonly quiescent than lower mass black holes,which have a higher probability to be actively accreting. Incomparison with an ideal volume-limited set of galaxies, an AGNsample will therefore contain a higher fraction of low mass blackholes. We have shown that this selection effect implies a biastowards a low M_•/M_Bulge ratio. This source of bias has not been discussed before.

• Luminosity bias: The selection of AGN based on their luminosities implies a preferential selection of black holes with higher than average masses, especially when the luminosity limit is a high one. This effect has already been discussed and analysed by others (e.g. Salviander et al. 2007; Lauer et al. 2007). Our analysis revealed that the luminosity bias depends sensitively on the details of the physically underlying distribution functions, in particular the active black hole mass function and the Eddington ratio distribution.

• Biases from measurement uncertainties: Both black hole mass and bulge property measurements suffer from measurement and calibration errors. In combination with selection limits they may introduce additional biases.

• AGN evolution bias: The AGN distribution functions that regulate the magnitude of the bias are evolving themselves as a function of redshift, as implied by the evolution of the AGN luminosity function. This evolution will change the expected sample bias with redshift.

These effects alter the bivariate probability distribution and can bias the conclusions drawn from observational studies. It is however important to clarify what is meant by saying that measurements of the M_•-bulge relations are biased. These relations can be represented in a number of ways that differ in complexity. (1) The most comprehensive approach would be to investigate the full bivariate distribution of galaxy properties and black hole masses. This distribution is modified by the AGN selection, and adopting the observed relation as the measurement thus constitutes a bias with respect to the true relation. (2) The most basic approach involves computing simply the mean offset of the entire sample from the local relation. We showed that this quantity is equally biased by the AGN selection. (3) An intermediate way is to look at the projected M_•-bulge relation, defined as the mean black hole mass for a given galaxy property. This relation also suffers from selection biases, caused by excluding black holes with certain properties from the sample during the observational process. (4) Alternatively, one can evaluate the inverse projected relation, thus the mean galaxy property at a given black hole mass. While for our parameterisation of the bivariate distribution this projection deviates from the true relation even without any selection effects, it is not modified by the AGN selection effects outlined above. Therefore this relation has the potential to be used as an unbiased estimator for the M_•-bulge relation, as long as additional selection on the host galaxy properties do not play a role. Among these four outlined approaches, options (2) and (3) are most commonly adopted in the literature.

It would clearly be desirable if we could properly model and correct the observations for the known sample selection effects. In order to do so, it is necessary to know the relevant underlying distribution functions, i.e. the distribution of spheroid properties such as velocity dispersions or masses; the active fraction as a function of BH mass, or alternatively the active black hole mass function; and the Eddington ratio distribution function. Thanks to recent observational progress, these distributions are now reasonably well known for the local universe; but they are still rather ill-constrained for higher redshifts.

We employed our formalism together with a set of best-guess distribution functions to investigate how much existing studies of the M_•-bulge relation in low-redshift AGN samples might suffer from selection biases. For the heterogeneous sample of AGN with reverberation mapping data we estimated a mean offset of  ~ 0.1 dex in M_• at given bulge property, chiefly due to the necessary flux limitations imposed in the process of selecting these objects for observations. This is a small but not entirely negligible effect, especially in view of the important role of this sample to serve as absolute calibration of the virial method for estimating back hole masses. We also confirm that slope and scatter are not significantly affected.

Considering the more general case of low-redshift AGN samples with measured bulge properties and BH masses estimated by the virial method, we find that the observed M_•-bulge relations suffer from a strongly mass-dependent offset from the true relation. As a consequence, the observed relation of M_• at given bulge property is flattened noticably. Such a flattening has been noted several times by observers (Greene & Ho 2006; Kim et al. 2008; Bentz et al. 2009b; Woo et al. 2010); it can be explained in our framework without having to invoke intrinsic differences in the BH-bulge relation between active and inactive galaxies.

Extending this type of analysis to higher redshifts is hampered by our poor current knowledge of the various relevant distribution functions. In order to constrain the evolution in the M_•-bulge relations it is crucial to not only obtain larger observed samples, but also to improve on measuring the underlying distributions. For example, not knowing the mass dependence of the active fraction at high redshifts will fundamentally limit our ability to correct for the biases associated with this effect. Furthermore, the selection function of the studied sample needs to be known; for an ill-defined heterogeneous sample the chances to properly account for the selection effects are heavily reduced.

We briefly reviewed a number of recent observational attempts to study the M_•-bulge evolution. The result is sobering: In no case do we find statistically significant evidence for an evolving M_•-bulge relation. While the observed apparent offsets in M_• /M_Bulge  from the local relation can be quite large, the sample selection bias estimated from our formalism is typically of the same magnitude. This does of course not exclude the possibility of real evolutionary effects that might even be already visible in some of the data; but our ability to distinguish between real and artificial trends caused by selection biases is presently not sufficient.

A possible route to circumvent several of the most problematic issues connected with AGN selection is the approach to focus on the mean galaxy property at a given black hole mass. However, while this projected distribution is formally unaffected by the AGN-specific selection biases, it also does not directly yield the intrinsic relation. In order to reconstruct the true relation in this case it is necessary to know the spheroid distribution function and the intrinsic scatter in the M_•-bulge relation. Neither is well constrained at high redshifts, but these quantities are observables, so that this approach may become very powerful in the future.

Acknowledgments

We thank Chien Peng for fruitful discussions on this manuscript. We acknowledge support by the Deutsche Forschungsgemeinschaft under its priority programme SPP1177, grant Wi 1369/23.

References

Appendix A: Validation of the bivariate probability distribution

Our chosen parameterisation of the bivariate probability distribution Ψ(s,μ), given as (A.1)(cf. Eq. (1)), follows common practice, but it is not the only possible solution. In particular, this parameterisation produces the visual impression of an upturn in the contours towards high values of s, corresponding to an upwards curved ridge line of the distribution. This upturn is simply a consequence of the fact that the space density of galaxies is steeply declining towards high masses or velocity dispersions. Since the quantitative conclusions in this paper depend on the appropriateness of the adopted form for Ψ(s,μ), it is worthwhile to test if our parameterisation is actually consistent with current observations. Here we perform such a test.

In Fig. A.1 we superimpose our intrinsic bivariate probability distribution Ψ(s,μ), which includes an estimate of the relevant bulge property distribution function, with data of galaxies for which dynamical black hole mass measurements are available. We consider the distribution of M_• versus the three most commonly measured galaxy properties: stellar velocity dispersions σ_∗, K band absolute magnitudes M_K, and stellar bulge masses M_bulge. For the M_•-σ_∗ relation we incorporated the stellar velocity dispersion distribution function from Sheth et al. (2003) and compare the result to the observational sample from Gültekin et al. (2009). For the M_•-M_K relation we employed the 2MASS K-band luminosity function from Kochanek et al. (2001) and show the sample from Hu (2009) for comparison, and for the M_•-M_bulge relation we use our estimate for the spheroid mass function, based on galaxy mass functions from Bell et al. (2003), and compare this with the sample from Häring & Rix (2004).

It is important to note that in this comparison we intentionally do not consider any sample selection effects. Thus, the distribution of data points in the three diagrams does not, and cannot, agree with the overall probability distribution. The test that we want to perform focuses on the upper right corner of each diagram, and the question is: do the data follow the mentioned upturn in the ridge line? While there is of course considerable observational and intrinsic scatter in the data, a visual inspection of Fig. A.1 clearly confirms that indeed there is such a trend.

We can perform a more quantitative test of the adopted ansatz by separately looking at the two univariate projections of the bivariate distribution,  ⟨ μ ⟩ (s) (Eq. (6) – represented by the thick solid lines in Fig. A.1), and  ⟨ s ⟩ (μ) (Eq. (7) – shown by the red dashed-dotted lines. To compare these curves with the measurements we computed the mean μ in bins of s and the mean s in bins of μ from the data, shown in Fig. A.1 as black open squares and red open triangles, respectively. Towards high values of s there is a clear trend for the triangles to be located above and to the left of the squares, and the predicted curves trace the points quite well. We performed a simple χ² test on all points above a certain minimum value of s or μ, respectively, and found that the binned data are consistently described by the predicted curves, with a reduced χ² always of order unity independently of the chosen cutoff. We conclude that our adopted form for the bivariate probability distribution Ψ(s,μ) is consistent with the current observations.

All Figures