Mon. Not. R. Astron. Soc. 000, 1–13 (2013)

Printed 31 January 2014

(MN LATEX style file v2.2)

arXiv:1312.5610v2 [astro-ph.CO] 30 Jan 2014

Accounting for selection effects in the BH-bulge relations: noevidence for cosmological evolution
Andreas Schulze1,2⋆ and Lutz Wisotzki3
1 Kavli Institute for Astronomy and Astrophysics, Peking University, 100871 Beijing, China
2 Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU, WPI), Todai Institutes for Advanced Study,

The University of Tokyo, Kashiwa 277-8583, Japan
3 Leibniz-Institut für Astrophysik Potsdam (AIP), An der Sternwarte 16, 14482 Potsdam, Germany

Accepted 2013 December 17. Received 2013 December 17; in original form 2013 November 20

ABSTRACTThe redshift evolution of the black hole - bulge relations is an essential observational
constraint for models of black hole - galaxy coevolution. In addition to the observational challenges for these studies, conclusions are complicated by the influence of
selection effects. We demonstrate that there is presently no statistical significant evidence for cosmological evolution in the M• -bulge relations, once these selection effects
are taken into account and corrected for. We present a fitting method, based on thebivariate distribution of black hole mass and galaxy property, that accounts for theselection function in the fitting and is therefore able to recover the intrinsic black hole bulge relation unbiased. While prior knowledge is restricted to a minimum, we at least
require knowledge of either the sample selection function and the mass dependence of
the active fraction, or the spheroid distribution function and the intrinsic scatter inthe black hole - bulge relation. We employed our fitting routine to existing studies ofthe M• −bulge relation at z ∼ 1.5 and z ∼ 6, using our current best knowledge of the
distribution functions. There is no statistical significant evidence for positive evolution
in the M• − M∗ ratio out to z ∼ 2. At z ∼ 6 the current constraints are less strong,
but we demonstrate that the large observed apparent offset from the local M• −bulge
relation at z ∼ 6 is fully consistent with no intrinsic offset. The method outlined here
provides a tool to obtain more reliable constraints on black hole - galaxy co-evolution
in the future.
Key words: Galaxies: active - Galaxies: nuclei - quasars: general

1

INTRODUCTION

The tight correlation between the mass of a supermassive
black hole and the properties of its host galaxies’ spheroid
component is now well established in the local Universe (e.g. Magorrian et al. 1998; Ferrarese & Merritt
2000; Gebhardt et al. 2000; Tremaine et al. 2002;
Marconi & Hunt 2003; Häring & Rix 2004), indicatinga coeval growth history of galaxies and their central black
holes. However, the detailed role of this correlation for
our understanding of galaxy evolution and black hole
growth is still an open question. The standard theoreticalscenario explains this relation by joint triggering of starformation and black hole activity via major mergers, while
the black hole growth is self-regulated via active galactic
nuclei (AGN) feedback that shuts off star formation andquenches the accretion on to the black hole (e.g. Silk & Rees⋆

E-mail: andreas.schulze@ipmu.jp

1998; Kauffmann & Haehnelt 2000; Di Matteo et al. 2005;
Sijacki et al. 2007; Somerville et al. 2008). Alternatively,
it may be possible to reproduce the observed black hole
- galaxy correlations without the need for self-regulation
via AGN feedback, either as a natural consequence of acommon merger history (Peng 2007; Jahnke & Macciò
2011), or in a black hole growth scenario via gravitational
torques (Anglés-Alcázar et al. 2013).
An observational constraint to distinguish between
these alternative scenarios and to constrain individual theoretical models (e.g. Croton 2006; Robertson et al. 2006;
Lamastra et al. 2010; Booth & Schaye 2011; Dubois et al.
2011) is provided by the redshift evolution of the M• −bulge
relations. However, direct dynamical black hole mass measurements are not feasible at higher redshifts. While integrated constraints can be obtained from the study ofthe black hole mass function (BHMF) and the AGN luminosity function (Merloni et al. 2004; Hopkins et al. 2006;

2

A. Schulze and L. Wisotzki

1.5

Indeed, sample selection is a particular important issue
for these studies, as selection effects are almost inevitable(Lauer et al. 2007; Shen & Kelly 2010; Schulze & Wisotzki
2011; Volonteri & Stark 2011; Portinari et al. 2012). In particular, Lauer et al. (2007) argued that intrinsic scatter inthe M• -bulge relation and a steep exponential cutoff in thegalaxy distribution function will lead to biased apparent relations. This bias is exacerbated by the effects of a bright flux
limit in the employed AGN samples. In Schulze & Wisotzki(2011, hereafter SW11) we presented a common framework,
based on the bivariate distribution function of black holemass and spheroid property, to investigate and model theluminosity bias and other selection biases on the M• -bulge
relations. We additionally discussed an active fraction bias,
which is introduced if only active black holes are selected andif the probability to be in an active stage depends on blackhole mass. It is possible to account for these effects if (1) the
selection function of the respective sample is well known, and(2) the underlying distribution functions, such as the galaxy
distribution function, the active BHMF and the Eddington
ratio distribution function (ERDF), are known. However, inparticular at high redshift this is currently not the case.
Therefore, we are facing the problem that even if weare able to obtain reasonably large AGN samples with individual black hole mass estimates and host galaxy property
measurements, the interpretation of such results in terms ofevolution or non-evolution is challenging due to the sample
selection effects. In this paper we present a practical solution to this dilemma. We build on the framework presented
in SW11 and give a quantitative assessment of, and a correction procedure for sample selection effects on observations of
the M• -bulge relations. We outline a maximum likelihood
fitting approach that directly incorporates information on
the selection function into the fitting. We also restrict priorknowledge of the underlying distribution functions to a minimum. The method is able to reconstruct the true unbiased

M
offset in log

The underlying assumption for this approach is thatbroad line AGN obey the same relationship as non-active
galaxies (e.g. Gebhardt et al. 2000b; Woo et al. 2013).
For broad line AGN the black hole mass can be estimated employing the ’virial method’ (e.g. McLure & Jarvis
2002; Vestergaard & Peterson 2006). The main challenge
for these studies is to determine the properties of the
AGN host galaxy, targeting either the velocity dispersion(Shields et al. 2003; Salviander et al. 2007; Woo et al. 2008;
Canalizo et al. 2012; Hiner et al. 2012; Salviander & Shields
2013), the galaxy luminosity (Peng et al. 2006; Decarli et al.
2010; Bennert et al. 2010; Targett et al. 2012) or thestellar mass (McLure et al. 2006; Schramm et al. 2008;
Jahnke et al. 2009a; Merloni et al. 2010; Nesvadba et al.
2011; Cisternas et al. 2011). While early studies suggested
a clear trend of positive evolution in the M• -bulge relation,
they preferentially focused on luminous quasars and theirhosts, sampling from the bright end of the AGN luminosity function. More recent studies, using fainter AGN, tend tofind mild or no evolution (Jahnke et al. 2009a; Merloni et al.
2010; Cisternas et al. 2011; Salviander & Shields 2013;
Schramm & Silverman 2013). A compilation of current literature results is shown in Fig. 1.

•

Zhang et al. 2012), the only route for observational studies
on individual objects is to use broad line AGN samples.

1.0
0.5
0.0

Peng+06
Jahnke+09
McLure+06
Bennert+10
Decarli+10
Cisternas+11
Schramm+08

-0.5
0

1

2

3

z 4

Schramm+13
Merloni+10
Nesvadba+11
Targett+12
Riechers+08
Maiolino+07
Wang+105

6

Figure 1. Compilation of results obtained from previous studies
on the apparent offset in black hole mass at a given galaxy massfrom the local relation from Häring & Rix (2004). Overall, an apparent trend of an increasing M• /M∗ ratio is suggested. However,sample selection effects need to be considered to reconstruct the
intrinsic offset.

M• -bulge relation for an observed sample and thus largely
overcomes the limitations imposed by the sample selection.
The paper is organized as follows: In section 2 we
present the general maximum likelihood fitting method andits application to the M• -bulge relation for AGN. In section 3 we verify the robustness and uncertainties of themethod via Monte Carlo simulations. Section 4 illustrates
the method by its application to representative previous
studies and discusses its ramifications for the evolution of
the M• -bulge relation. We conclude in section 5.

2
2.1

THE METHOD
Maximum likelihood fit in consideration of a
selection function

While the main goal of the method presented here is to
study the M• −bulge relation, it is also applicable to moregeneral cases. Identifying and characterizing correlations between two observables is a fundamental aspect of astrophysical research. Given a joint distribution of two parameterswith non-zero intrinsic scatter, this correlation will often beaffected by the specific sample selection. We therefore firstoutline the general concept of the approach, before applying
it to the M• −bulge relation.
Let us assume that there exists a correlation between
two observational parameters x and y for a class of objects.
The space density of these objects is given by a bivariate
distribution function Ψ(x, y). Marginalizing Ψ(x, y) over x
gives the distribution function of the variable y, and viceversa. In the case of the M• −bulge relation x is for example
the stellar velocity dispersion log σ∗ and y is the black holemass log M• . Thus marginalizing their bivariate distributionfunction over log σ∗ gives the BHMF, while marginalizing
over log M• will provide the galaxy stellar velocity dispersion
function (SVDF).
We here assume the following parameterization for theInfering the intrinsic BH-bulge relation
bivariate distribution
Ψ(x, y) = g(y | x) Φx (x) ,(1)

where Φx (x) is the distribution function in x, e.g. the galaxy
SVDF, and g(y | x) gives the probability of finding the value
y given x. This is the intrinsic correlation between x andy. We here specifically assume a linear relation y = α + βx
with log-normal intrinsic scatter σ, i.e.
1(y − α − βx)2g(y | x) = √
.
(2)
exp −
2σ 2
2πσ
The bivariate distribution function Ψ(x, y) describes the intrinsic joint distribution function of x and y in the universe.
In practice, this will not be equal to the observed bivariate distribution for a sample, since this sample is almost
always affected by implicit or explicit selection effects. We
incorporate these selection criteria into a selection functionΩ(x, y, θ), where θ refers to the set of additional parameterspresent as selection criteria, such as redshift or luminosity.
The multivariate distribution Ψ(x, y, θ) for an observed sample is then given byΨo (x, y, θ) = Ω(x, y, θ) Ψ(x, y) pθ (θ) ,(3)

where pθ (θ) is a set of normalized distribution functionsof the parameters θ. The respective bivariate distribution
Ψo (x, y) is given by integration over θ. In general Ψo (x, y) 6=
Ψ(x, y). This means the correlation between x and y, e.g. theM• − σ∗ relation, fitted directly to the observations is notequal to the intrinsic correlation. Ignoring this inequality
will lead to a bias (see SW11).
In practice, we are mainly interested in the intrinsic correlation between x and y, i.e. g(y | x) and less in the full joint
distribution function. The information of this intrinsic correlation is contained in the conditional probability function,
but in general is not equal to it. There are two options todefine the observed correlation between x and y. First, therelation of y at a given x, and second the relation of x at agiven y. These are given by the conditional probabilitiesp(y | x, θ) = R
p(x | y, θ) = R

Ψo (x, y, θ)Ψo (x, y, θ) dy(4)Ψo (x, y, θ),Ψo (x, y, θ) dx(5)where p(y | x, θ) gives the probability of finding the propertyy (e.g. black hole mass) in a galaxy with properties x (e.g.
σ∗ ) and θ, while p(x | y, θ) is the probability of finding agalaxy with x for given y and θ. We will refer to these as
forward and inverse probabilities, respectively.
Let us assume the special case there are no selection
effects on y. When we use the bivariate distribution given
by Equation (1) the forward conditional probability is given
by pint (y | x) = g(y | x), i.e. the intrinsic underlying relation
is directly recovered in the conditional probability. For the
inverse conditional probability it follows in the absence ofany selection effects
pint (x | y) = R

g(y | x) Φx (x).
g(y | x) Φx (x)dx(6)

This is already a deviation from the intrinsic relation g(y | x),even without applying any selection criteria to the sample. In the context of the M• −bulge relation the stellar

3mass function exponentially decreases at the high mass end,
which together with the intrinsic scatter leads to a deviationof pint (M∗ | M• ) at the high mass end (Lauer et al. 2007,
SW11). This can already qualitatively explain the steeper
slope found in the M• −σ∗ relation for galaxies with dynamical black hole mass measurements when fitting the inverseM• −σ∗ relation p(σ∗ | M• ) instead of the relation p(M• | σ∗ )(Graham et al. 2011; Park et al. 2012).
In practice selection effects are present and can produce a bias, i.e. p(y | x) 6= g(y | x). However, if we know theselection function of our sample, we can incorporate thisknowledge directly into the fitting of the conditional probability to recover the intrinsic relation g(y | x). The maximum
likelihood technique provides a conceptually simple way toachieve this. This method aims atQminimizing the likelihood
function S = −2 ln L, with L = i li being the product of
the likelihoods for the individual measurements.
Both conditional probabilities, p(y | x) and p(x | y), canserve as likelihoods, i.e. we either minimize

Sx = −2NX
i=1

Sy = −2ln p(yi | xi , θi ) , or

N
X
i=1

ln p(xi | yi , θi ) ,

(7)

(8)

to obtain the best-fitting parameters of g(y | x). We here refer to the minimization of Sx as forward ML regression andof Sy as inverse ML regression. We emphasize that thesemaximum likelihood fits go beyond a simple linear regression, since we are not only fitting a model to the data, but
incorporate additional information about the data set andthe underlying distributions. Thereby, the used conditional
probability function can have an arbitrarily complex form.
We will refer to the more standard approach of fitting theconditional probability p(µ|s) without including the effectsof sample selection as ’simple regression’.
Which of the two maximum likelihood fit options (forward/inverse) is the preferred choice depends on the specificselection effects, and on the underlying distributions. For example, when there are no selection effects then forward regression is preferable, since the fitted conditional probability
is equal to the intrinsic relation g(y | x). In this case the maximum likelihood fit is equal to a simple regression approach,used in other studies to fit the M• − σ∗ relation for quiescent galaxies (e.g. Gültekin et al. 2009; Schulze & Gebhardt
2011; Park et al. 2012). The inverse ML regression (using
Equation 6) requires additional information about the distribution function Φx (x) and is therefore more complicatedto use. However, in the presence of selection effects on y
the inverse ML regression might become preferable, since it
can be unaffected by the selection effects, in contrast to theforward ML regression. We will illustrate this below for theM• −bulge relation for active galaxies.
While other applications of this method are possible,
in the following we focus on the specific case of studying
the M• −bulge relation at high redshift via broad line AGN
samples.

4

A. Schulze and L. Wisotzki

2.2

Application to the M• −bulge relation forbroad line AGN samples

regression then follows as

In SW11 we discussed the qualitative influence of sample
selection effects on observations of the M• −bulge relation
for broad line AGN samples in detail within our frameworkof the bivariate distribution function.
The M• −bulge relation is given by µ = α + βs, with
µ = log M• being the logarithmic black hole mass and s
is the logarithm of the host galaxy property, e.g. log σ∗ orlog Mbulge , and α and β are the normalization and the slopeof the relation. Thus in the above terminology x = s and y =µ. With the parametrization of equation (1), the bivariate
distribution function isΨ(s, µ) = g(µ | s) Φs (s) ,(9)where Φs (s, z) is the spheroid distribution function andg(µ | s, z) is the intrinsic M• −bulge relation. Integrating thisbivariate distribution function over s gives the total BHMF.
For an AGN sample the main selection criteria will be
the AGN flux. Therefore, the selection function will dependon AGN luminosity and redshift. The AGN luminosity isdetermined by the black hole mass M• and the accretionrate, the latter is represented by the Eddington ratio λ =Lbol /LEdd , i.e. log Lbol = µ + log λ + 38.1. This implies aselection function Ω(µ, λ , z), depending on black hole mass,
Eddington ratio and redshift.
The corresponding distribution function is the bivariate
distribution function of black hole mass and Eddington ratio
for active black holes Ψ̃(µ, λ), with its projections the active
BHMF and the ERDF. We here make the assumption thatthe bivariate distribution function of black hole mass andEddington ratio is separable into an independent BHMF
and ERDF, i.e.
Ψ̃(µ, λ, z) = pλ (λ ) Φ•,a (µ, z) ,(10)

where pλ (λ , z) is the ERDF, normalized to one andΦ•,a (µ, z) is the active BHMF. This assumption is motivatedby results on the local bivariate distribution function of µand λ (Schulze & Wisotzki 2010) and by the observation of
an Eddington ratio distribution of AGN that is independentof stellar mass (Aird et al. 2012; Bongiorno et al. 2012), suggesting also an independence of µ [but see Kelly & Shen
(2013) for evidence of a µ dependence of the ERDF]. We
comment on the consequences of relaxing this assumption
below.
The active BHMF is related to the total BHMF via theactive fraction, which gives the probability for a black hole to
be in an active (type 1 AGN) state as a function of black holemass and redshift. It is defined by the ratio of active BHMF
to total BHMF, pac (µ, z) = Φ•,a (µ, z)/Φ•,t (µ, z). The latteris directly linked to the bivariate distribution of µ and s(Equation (9)).
Following Equation (3), the multivariate probability
distribution Ψo (s, µ, λ , z) of galaxy property, black holemass, Eddington ratio and redshift is then given byp(µ | s, λ , z)==Ψo (s, µ, λ , z)Ψo (s, µ, λ , z) dµΩ pac (µ, z) g(µ | s, z)R.
Ω pac (µ, z) g(µ | s, z) dµR(12)The ERDF pλ (λ ) and the galaxy distribution function Φs (s)cancel out and are not required for the fitting. We require
only knowledge of the selection function and of the massdependence of the active fraction to obtain the intrinsic
M• −bulge relation from the fit to the conditional probability. The free parameters α, β and σ in the intrinsic relation g(µ | s, z) are obtained by minimizing Equation 7. If the
ERDF is mass dependent it will not cancel out and knowledge of the ERDF is additionally required.
Second, the conditional probability for the inverse ML
regression isp(s | µ, λ , z)

=
=

R
R

Ψo (s, µ, λ , z)Ψo (s, µ, λ , z) ds(13)g(µ | s, z) Φs (s).
g(µ | s, z) Φs (s) ds

This is identical to the case without AGN selection effects(Equation 6), i.e. the conditional probability p(s | µ) is notaffected by the AGN selection. We already discussed thisfact in SW11, but we will use it here to recover the intrinsic
relation by the maximum likelihood fit. The only quantitythat has to be known to recover the intrinsic relation is thegalaxy distribution function.
Both fitting options have their pros and cons. Whilethe latter is largely unaffected by AGN selection effects, it
is susceptible to the shape of the galaxy distribution function
and also more susceptible to the intrinsic scatter in the relation itself. The former does not require any knowledge aboutΦs (s), but demands a proper understanding of the sample
selection and also of the underlying distribution function of
AGN or at least the active fraction for the AGN sample.
Both regression methods constitute an improvement forthe determination of the intrinsic M• −bulge relation, compared to a standard regression that does not include knowledge of the selection function or galaxy distribution function,
as we will illustrate below.

2.3

Incorporating measurement uncertainties

In practice, the situation becomes more complicated, since
black hole mass and galaxy property measurements havemeasurement uncertainties. Virial black hole masses are associated with a non-negligible uncertainty of ∼ 0.3 − 0.4 dex(Vestergaard & Peterson 2006; Park et al. 2012), while also
AGN host galaxy properties at high z have significant measurement errors. These lead to second-order dependenciesof the conditional probabilities. In SW11 we already studied the consequences of measurement uncertainties on thebivariate distribution and the M• −bulge relations. Here, we
discuss the ramifications of measurement uncertainties ondV the presented fitting approach.
Ψo (s, µ, λ , z) = Ω(µ, λ , z) pac (µ, z) pλ (λ , z) g(µ | s, z) Φs (s, z),
We assume that the black hole mass estimated by thedz(11)virial method µo is drawn from a log-normal probability
where Ω(µ, λ , z) is the AGN selection function of the respecdistribution around the true black hole mass µ, with distive sample. The conditional probability for the forward ML
persion σvm , representing the uncertainty in the virial mass

Infering the intrinsic BH-bulge relation
estimate,g(µo | µ) = √1(µo − µ)
exp −
2
2σvm
2πσvm

2



.

(14)

For the galaxy property s we also assume a log-normal error
distribution g(so | s) with dispersion σs .
The bivariate distribution function for bulge property
and virial black hole mass is thenZΨo (so , µo ) = g(so | s) g(µo | µ) Ψo (s, µ) dsdµ .
(15)
The bivariate distribution function is smoothed out by themeasurement uncertainties. The conditional probability forthe forward ML regression isp(µo | so , λ , z) =RΨo (s, µ, λ , z) g(so | s) g(µo | µ) dsdµR=Ψo (s, µ, λ , z) g(so | s) g(µo | µ) dsdµdµoRΩ pac (µ, z) g(µ | s, z) Φs (s) g(so | s) g(µo | µ) dsdµR.
Ω pac (µ, z) g(µ | s, z) Φs (s) g(so | s) dsdµ(16)Assuming we know the measurement uncertainties σvm andσs , these can be directly incorporated into the fit. However,
we now introduce a dependence on the galaxy distributionΦs (s) into the forward ML regression which was canceledout before. To avoid this dependence, we can make the simplifying assumption Φs (s) ≈ Φs (so ). This approximation isjustified for a small measurement error in s and in particular in the flat part of the galaxy distribution function. It
will break down for high stellar velocity dispersion or highstellar mass, i.e. at the exponential decline of the galaxy distribution function. It can be applied to samples dominatedby lower mass galaxies. Equation 16 then simplifies toRΩ pac (µ, z) g(µ | s, z) g(so | s) g(µo | µ) dsdµRp(µo | so , λ , z) ≈.
Ω pac (µ, z) g(µ | s, z) g(so | s) dsdµ(17)
We will refer to both options as forward ML regression withmeasurement errors.
The conditional probability for the inverse ML regression is given byp(so | µo , λ , z) =RΨo (s, µ, λ , z) g(so | s) g(µo | µ) dsdµR=Ψo (s, µ, λ , z) g(so | s) g(µo | µ) dsdµdsoRΩ pac (µ, z) g(µ | s, z) Φs (s) g(so | s) g(µo | µ) dsdµR.
Ω pac (µ, z) g(µ | s, z) Φs (s) g(µo | µ) dsdµ(18)

For this case, the dependence of the conditional probabilityon the active fraction and on the selection function is notcanceled out as it was the case without measurement uncertainty in µ. This effectively prevents ignoring the AGN
selection function for the fit of the M• -bulge relations forobservational samples.
We can again derive an approximation to the probability distribution by making the simplifying assumptions
pac (µ) ≈ pac (µo ) and Ω(µ) ≈ Ω(µo ), for which Equation (18)
reduces toRg(µ | s, z) Φs (s) g(so | s) g(µo | µ) dsdµRp(so | µo , λ , z) ≈.
g(µ | s, z) Φs (s) g(µo | µ) dsdµ(19)
For this approximation no knowledge of the AGN selection

5function and the active fraction is required anymore. The
assumption on the selection function depends on the detailsof the sample selection, but in general will be true in mostcases, unless the object is close to the selection threshold(e.g. to the flux limit of the survey). The assumption on theactive fraction requires an active fraction with a weak or no
black hole mass dependence. This is, for example, the caseat z ∼ 1.5, as discussed in Section 4.1, but not in the local
universe (Schulze & Wisotzki 2010). Therefore, in practicesome idea about the active fraction is required in any caseto judge the validity of this approximation.
We will refer to these maximum likelihood fits as inverse
ML regression with measurement errors.
We here presented several options to fit for the underlying relation using a given sample with measurement uncertainties. Which of these is the best choice for the specificsample will depend on the details of the sample selection
and our knowledge of this selection function and the galaxy
distribution function. Formally, the most precise results are
obtained using the conditional probabilities for either theforward or inverse ML regression given by Equations 16 and18. This requires knowledge of the selection function, theactive fraction and of the galaxy distribution function. To
restrict the assumptions for the respective fit, we can also
use the approximations given by Equations 17 and 19. While
not precise, this approach can be justified if the selection
function and/or the underlying distributions are only poorly
known.

3

MONTE CARLO SIMULATIONS

To illustrate the method and test its reliability we carriedout Monte Carlo simulations of the M• − M∗ relation. For
our simulations we chose as representative redshift z = 1.5.
For the galaxy spheroid mass function we assume the estimate given in SW11, which is based on the stellar massfunction from Fontana et al. (2006), applying a rough correction for the bulge-to-total ratio to translate it into aspheroid mass function. We set up the simulations by firstdrawing a large number of objects from the spheroid massfunction at z = 1.5. Black hole masses are attributed according to the M• − Mbulge relation from McConnell & Ma(2013) with log-normal intrinsic scatter σ = 0.34. We here
assumed a constant active fraction, consistent with recent
results on the active BHMF at z ∼ 1.5 (Nobuta et al. 2012;
Kelly & Shen 2013, Schulze et al. in prep.) and empiricalmodels (Merloni & Heinz 2008; Shankar et al. 2009). For theactive black holes an Eddington ratio was drawn from thez ∼ 1.5 ERDF from Nobuta et al. (2012). Black hole massand Eddington ratio define the bolometric AGN luminosity.
We convert Lbol into i-band magnitude using the AGN spectral energy distribution (SED) from Richards et al. (2006a)
and the K-correction from Richards et al. (2006b). We applied three flux limits to the simulated samples, imin =20, 22
and 24 mag, respectively, and generated samples of various
sizes. For each flux limit and sample size we generated 1000
Monte Carlo realizations.
For simplicity, we first discuss the results without measurement errors. We will discuss the effect of adding measurement errors further below. Each simulation realization
was fitted by a forward and an inverse ML regression. Fur-A. Schulze and L. Wisotzki

6

µ

N

s

8.2

8.4

α

8.6

8.8

0.5

1.0

β

1.5

0.1

0.2

0.3

σ

0.4

0.5

8.2

8.4

α

8.6

8.8

0.5

1.5

1.0

β

0.1

0.2

0.3

σ

0.4

0.5

Figure 2. Demonstration of biases due to selection effects and their removal by our ML fitting approach. For a set of Monte Carlosimulations we assumed the ’true’ M• -bulge relation at z = 1.5 to be described by a normalization α, a slope β, and an intrinsic scatter
σ; these input values are marked by the vertical red lines in each panel. For 1000 simulated random samples of 100 objects each, with aflux limit of i < 22 mag, the panels show the best-fitting parameters from simple linear regression (green dashed histograms) and from
our ML regression method (black histograms – left subpanels: forward regression, right subpanels: inverse regression). The horizontal
offsets between the red lines and the centroids of the dashed histograms indicate clearly the presence of biases, while the estimates using
our ML regression method are largely unbiased.

s

1.6

8.8

0.4

σ

10

50 100

N

Obj

500

)fixed

8.2

8.4

8.2
0.6

i

0.2

1.0

0.8

α

min =20
min =22
imin =24
i

8.4

β,σ

( )fixed
β

σ

0.3

8.6

1.2

(

( )fixed
α

8.6

σ

8.8

1.4

0.4

8.0
10

50 100

N

500

Obj

8.0
10

50 100

N

Obj

500

10

50 100

N

500

Obj

Figure 3. Recovery of intrinsic relation parameters by using our forward ML regression method, accounting for selection effects. The
Monte Carlo simulations are the same as for Fig. 2, but for a range of sample sizes and flux limits. Each panel shows median (symbols)
and 68% confidence intervals (dashed lines) for the best-fitting results of normalization α, slope β, and intrinsic scatter σ, respectively.
The horizontal black line indicates the input value. Both α and β can be recovered without any systematic bias even from small samples,
whereas σ tends to be somewhat underestimated unless the sample is large enough.

thermore, for each regression we fitted the realizations witha different set of free parameters. We first kept all threeparameters (α, β, σ) free. Secondly we fixed σ and finally
we also fixed β to the input value, leaving only the normalization as free parameter. Additionally, we fitted eachrealization with a simple linear regression, i.e. not including
the selection effects, both with µ and with s as dependent
variable.
In Fig. 2, we illustrate the strength of our proposed
fitting approach (shown as black histograms, left for the forward ML regression and right for the inverse ML regression),compared to a simple linear regression (shown as green histograms). We use 1000 realizations of a sample of N = 100
objects at z = 1.5 with i < 22 mag. While the dispersion on
the fitting parameters is increased, the bias present in thesimple regression approach is removed. In this example we
leave all three parameters free in the fit.
In Fig. 3 we show the median and 68% confidence limits of the 1000 realizations for a forward ML regression fordifferent flux limits and sample sizes. The magnitude limitsof i = 20, 22, 24 mag correspond to bolometric luminosity

limits of log Lbol = 46, 45, 44.3, respectively. The left panel
shows that we systematically underestimate σ for a small
sample size when keeping it as a free parameter. This is dueto the fact that the fitting procedure determines the sample
standard deviation, which is known to be an underestimateof the population standard deviation. This underestimate forsmall sample sizes is already present in the conventional fit
without including the selection function. Since the intrinsic
scatter is an important parameter, when the selection effects
are included this underestimate is enhanced. Also the bestfitting solutions for the normalization α and the slope β will
be biased. Therefore, we suggest to fix the intrinsic scatterfor small sample sizes, NObj . 50, to avoid this bias. In thefollowing we fix the intrinsic scatter σ in the maximum likelihood fit. In this case, the slope β and normalization α arerecovered unbiased irrespective of sample size, as shown in
the other three panels of Fig. 3. Increasing the sample size
mainly reduces the statistical errors, shown by the dashedlines. Fixing also the slope and keeping the normalization
as only free parameter further reduces the statistical error,
while the input is recovered in the mean. We further see that

Infering the intrinsic BH-bulge relation
µ

1.6

8.8

0.4

σ

50 100

N

8.4

8.2
0.6

i

10

)fixed

8.2

i

0.2

1.0

0.8

α

min =20
min =22
imin =24

8.4β,σ( )fixed
β

σ

0.3

8.6

1.2

(

( )fixed
α

σ

8.8

1.4
8.6

7

0.4

8.0

500

10

50 100

Obj

N

500

8.0
10

Obj

50 100

N

500

10

Obj

50 100

N

500

Obj

Figure 4. Same as Fig. 3, but for an inverse ML regression. Again this verifies the reliability of the method to reconstruct the intrinsic
relation.

R

s

s

µ

R

)fixed

8.8

8.4

α

(

β,σ

8.6

min =20
min =22
imin =24

8.2

i
i

8.0
10

50 100

N

500

10

Obj

50 100

N

Obj

500

10

50 100

N

500

Obj

Figure 5. Reliability of our fitting method when measurement errors on the black hole mass and spheroid mass are present. We again
show the best-fitting normalization α (for fixed slope and intrinsic scatter) from Monte Carlo simulations at z = 1.5 for different flux
limits, but now with measurement errors added. The different panels show different variants for the fitting, as explained in the text. Leftpanel: using a forward ML regression with measurement errors. Middle panel: using a forward ML regression, but without accounting
for measurement errors in the fitting of the relation. Right panel: using an inverse ML regression without accounting for measurement
errors. The filled symbols show the median value for different applied magnitude limits, the solid lines show the 1σ confidence region forthe simulations with measurement errors, while the dashed line shows the 1σ confidence region of the fit when no measurement errors
were added to the simulated data.

the error is basically independent of the flux limit. However,this assumes perfect knowledge of the selection function andunderlying distribution. Deeper samples are less susceptible
to uncertainties in these and therefore come closer to thestatistical errors shown here.
In Fig. 4 we show the same for the inverse ML regression. Keeping the intrinsic scatter as a free parameter also
leads to an underestimate for small sample sizes and weagain suggest to fix it for small samples, NObj . 50, to avoid
a bias in the slope and normalization. But even when fixing
σ we find in our simulations a small bias in α and β for verysmall samples, NObj . 10. However, for such small samples
the errors on the best fit are also large. A joint determination of slope and normalization of the relation requires areasonably large sample. When fixing the slope β and only
determining the normalization α, an unbiased result withsmall error can be achieved even for relatively small sample
sizes, as demonstrated in the right panel of Fig. 4.
We now investigate the consequences of adding mea-

surement errors in the virial black hole masses and in thehost galaxy properties to the simulations. We here test ourprescription for including these uncertainties, as presented
in section 2.3 and compare them to a regression without
measurement errors. We added a log-normal random errorof dispersion σvm = 0.3 dex to the black hole mass and alog-normal error of dispersion σs = 0.3 dex to the spheroid
mass of our Monte Carlo realizations at z = 1.5. For eachrealization we fixed β and σ and only fit for the normalization α. We fitted these with the forward ML regression
with measurement errors, i.e. using Equation 16. The difference to the inverse ML regression with measurement errors
is marginal, in particular when the errors in µ and s are of
the same magnitude. To test the importance of the full correction for measurement errors, we also fitted the same realisations with the forward and inverse ML regression without
measurement errors.
The results are shown in Fig. 5. We show the median(filled symbols), the 68% confidence limits (solid lines) andA. Schulze and L. Wisotzki

4

APPLICATION TO OBSERVATIONALSTUDIESWe now apply our method to a few observational studies
from the literature that investigated the evolution in theM• −bulge relations. This is meant for illustration of themethodology and constraining power of the approach. An
exhaustive investigation of the literature on evolution in theM• −bulge relations in consideration of selection effects isbeyond the scope of this paper. We note that such an investigation would also be hampered by the fact that many
previous literature studies did not use well-defined samples,or at least samples for which a selection function (howevercomplicated) could be defined.
4.1

The M• − M∗ relation at z ∼ 1.5

First, we reinvestigated the QSO sample from Merloni et al.
(2010). They studied the M• − M∗ relation for 89 broad line
AGN from zCOSMOS in the redshift range 1.06 < z < 2.19
with estimated stellar masses for the host galaxies (including
upper limits). The selection function is to first order definedby an IAB -band flux limit of 22.5 mag.
To fit their M• − M∗ relation we need additional information about the underlying distribution functions. We
use the same assumptions as above for the Monte Carlo realizations, i.e. the stellar mass function from Fontana et al.
(2006), and a constant active fraction. To investigate the degree of evolution, we assume a redshift evolution in the zeropoint a(z) = a0 + γ log(1 + z). For the local relation we use

⊙]

9.5

9.0

• [

M M

also the 68% confidence limits for the realizations with zeroerror in µ and s (dashed lines). In the left panel the result for the forward ML regression with measurement errors
is presented. The input value is recovered in the median,
while the 1σ confidence region (solid lines) is increased by
an approximately constant factor of ∼ 2 for our simulation,compared to the case without measurement errors (dashedlines). The fitting method provides unbiased results if selection function, underlying distributions and the value of themeasurement error are all known. In the middle and rightpanels of Fig. 5 we show the results of applying the forward
ML regression and inverse ML regression without measurement errors to the simulations that have errors added. Whilein the forward ML regression the zero-point α is slightly biased towards a lower value, the inverse ML regression results
are fully consistent with the input value, only showing an increase in the 1σ confidence range. Thus, not fully including
the measurement errors has only a small, second order effecton the best-fitting solution.
The Monte Carlo simulations verify the reliability ofthe presented maximum likelihood fitting approach to recover the intrinsic M• − M∗ relation in the median from asample affected by sample selection effects. A simple linear
regression would provide a biased result in these cases. The
method is still reliable if there are measurement uncertainties in the black hole mass and galaxy property. While themost precise results are obtained when the uncertainty iswell known and taken into account, over- or underestimating or even ignoring the measurement uncertainty will only
have a small effect on the recovery of the intrinsic relation.

8.5

log

8

8.0

7.5

10.0

M M

10.5

log

11.0

∗ [

11.5

⊙]

Figure 6. Redshift evolution fit to the AGN sample from
Merloni et al. (2010). We fitted their data (blue circles) to theconditional probability distributions p(yo | xo , λ , z), determining
the evolution term γ. The red solid line shows the local relationfrom McConnell & Ma (2013), and the black dashed line gives
the best-fitting relation at the mean redshift of the sample. The
grey regions are the 1σ, 2σ and 3σ confidence regions.

first the relation by Häring & Rix (2004), with an intrinsic
scatter of σ = 0.3, consistent with the work by Merloni et al.
(2010), and second the recent results by McConnell & Ma(2013), with σ = 0.34. Upper limits are incorporated as described in Appendix A.
For the forward ML regression (using Equation 16), assuming measurement errors of 0.2 dex in the stellar massand 0.3 dex in the black hole mass estimates, we get γ =0.38+0.20−0.21 , using the Häring & Rix (2004) local relation. Using the relation from McConnell & Ma (2013) as zero-point,
we find γ = −0.35+0.22−0.23 (i.e. formally a negative evolution).
For the inverse ML regression and when using approximations for the probability distribution (Equations 17 and 19)
we get consistent results. The null hypothesis of no evolution in the M• − MBulge relation for both cases lies within
2 − 3σ of our results, thus it is not rejected with statisticalsignificance. The best fit relation, using the Häring & Rix(2004) local relation, is shown in Fig. 6.
Note that the given errors are statistical errors from the∆S contour, while uncertainties in the assumptions for thezero-point, the flux limit and the distribution functions are
not included. This implies that the uncertainty estimates
above are rather optimistic, since they only include statistical uncertainties and not additional uncertainties in theunderlying assumptions. The likelihood for the null hypothesis of no evolution is then even higher than the 2 − 3σmentioned above when those uncertainties are also includedin the confidence determination.
To get an estimate for the effect of uncertainties in theunderlying distribution function, we varied the assumed stellar mass function from Fontana et al. (2006) within theirInfering the intrinsic BH-bulge relation

4.2

High redshift QSOs

Studies of the most luminous QSOs, observed at highredshift found the largest apparent offsets from theM• −bulge relations, by a factor of 10 or more (Walter et al.
2004; Maiolino et al. 2007; Ho 2007; Riechers et al. 2008;
Schramm et al. 2008; Wang et al. 2010; Targett et al. 2012).
At first glance this seems to provide the most convincing evidence for strong evolution in the M• −bulge relations. However, these luminous QSO studies are only sampling fromthe bright end of the luminosity function. For this kind of
objects the bias caused by selection effects is most severe,
as discussed in SW11. Furthermore, current sample sizes atz > 3 are small, inhibiting firm statistical conclusions.
We now endeavor to test quantitatively if the strong
apparent offset found in these studies provides any statistically significant evidence for an intrinsic offset, once weaccount for sample selection. We focus on the highest redshift regime currently studied z ∼ 6 and the constraints thatcan be drawn on evolution in the M• −bulge relations.
We used the sample from Wang et al. (2010), and restricted it to the six objects selected from the Sloan Digital Sky Survey (SDSS) z > 5.7 quasar survey (Fan et al.
2000, 2001, 2003, 2006), which approximately might be considered a well-defined sample. For this sample we assume
a selection function given by the survey magnitude limitof z ∗ < 20.2 mag. Black hole masses are for some of theobjects estimated from the broad lines and for the rest
from the bolometric luminosity, assuming Eddington limited
accretion. Bulge masses are derived from dynamical massmeasurements (we assume an average inclination angle of40◦ ), and stellar velocity dispersions were estimated from

10.0

9.5

log M • [M ⊙ ]quoted uncertainties. Furthermore, we varied the bulge-tototal ratio between our default correction and no correction
at all, to investigate the importance of the transformationof the galaxy mass function to a spheroid mass function.
The uncertainty in the stellar mass function adds an uncertainty of ±0.05 to the best-fitting γ, while not applying anbulge-to-total correction increases the best-fitting solution
by 0.04.
Even more important is a precise calibration of the local
zero-point, as it is illustrated by the opposite trends of evolution found for the local relations by Häring & Rix (2004)
and McConnell & Ma (2013), which is mainly caused by theincrease in the zero-point of the M• − MBulge relation of
McConnell & Ma (2013).
While the local zero-point of the M• − MBulge relation
increased in the most recent revised determinations of the relation (McConnell & Ma 2013; Kormendy & Ho 2013), thef factor and thus the calibration to the AGN samples
is not strongly affected by this revision (Park et al. 2012;
Woo et al. 2013). Using this revised M• −MBulge relation significantly reduces the apparent trend of positive evolution,a fact recently noted by Kormendy & Ho (2013). Combinedwith the correction for the bias caused by sample selection
this even suggests or is at least consistent with a negative
evolution of the M• − MBulge relation where galaxy bulges
grow first and the black holes have to catch up.
In summary, we see no statistically significant evidence
for evolution in the M• − MBulge relation based on this dataset, consistent with our more qualitative analysis in SW11.

9

9.0

8.5

8.0

10.0

log Mdyn [M ⊙ ]

10.5

11.0

11.5

12.0

1.8

log σ ∗ [km s−1 ]
2.0

2.2

2.4

2.6

Figure 7. Redshift evolution fit of the M• − MBulge relation (leftpanel) and M• − σ∗ relation (right panel) to the main SDSS QSO
sample given in Wang et al. (2010). We used a forward regression,accounting for selection effects and including measurement errors.
The blue circles show the data, the red solid line is the local
relation, the dashed line is the best-fitting relation at the mean
redshift of the sample. The grey areas indicate confidence regionsof 1σ, 2σ and 3σ, respectively.

the width of molecular CO line emission. There are still significant systematic uncertainties associated with these measurements and how they relate to galaxy masses and velocity dispersion measured for local galaxies. However, we here
only discuss the statistical aspects of the sample.
The dynamical range in black hole mass and AGN luminosity covered by this sample is small and there is basicallyno apparent correlation between M• and either σ∗ or Mdyn .
We therefore fixed the slope and intrinsic scatter and only
fit for an evolution in the normalization, again parametrized
by a parameter γ.
For the stellar mass function we assume the mass function for star forming galaxies at z ∼ 5 from Lee et al.
(2012). We again assume a constant active fraction. Using the forward regression with measurement errors we findγ M∗ < 1.10 at 1σ confidence, with a best fit of γ M∗ = −0.02
for the M• − MBulge relation, using the local relation from
Häring & Rix (2004). We also varied the mass function according to their uncertainties, as given in Lee et al. (2012).
The best-fitting γ M∗ ranges from −0.74 to 0.59 with an
1σ confidence interval [−3.24, +1.63]. For the local relationfrom McConnell & Ma (2013) as zero-point we find a best
fit of γ M∗ = −0.58 with γ M∗ < 0.71 at 1σ confidence. At
any rate the inferred value of γ is always consistent with
zero.
Assuming a SVDF is more problematic. The SVDF
is well determined for z = 0 (Sheth et al. 2003; Chae
2010) and there are estimates out to z ≈ 1.5 (Chae
2010; Bezanson et al. 2011), but it is observationally unconstrained at higher redshifts. On the other hand, it is notclear at all to what degree the CO line width traces the stellar velocity distribution at these redshift. We here simply
try to bracket our ignorance about the SVDF at z ∼ 6 andinvestigate their consequences on the results for the M• − σ∗relation. Note that for the correction only the shape of the
distribution function is important, in particular the locationof the Schechter function cutoff, not its normalization.
We started from the evolution of the stellar mass function between z = 0 and z ∼ 5, comparing the mass functions
from Bell et al. (2003) and Lee et al. (2012). They show a

10

A. Schulze and L. Wisotzki

decrease of the break of the mass function of ∼ 0.6 dex. Assuming virialized systems, we can convert this into a shift for
the SVDF. However, this also requires a parameterization of
the size evolution of galaxies. At z ∼ 2 massive galaxies are
more compact by a factor of 3−5, compared to local galaxies
of similar mass (e.g. Daddi et al. 2005; Trujillo et al. 2007;
van Dokkum et al. 2008). Assuming similar compact galaxies at z ∼ 5−6 this approximately balances out the evolution
in the mass function, i.e. no evolution in the break of the
SVDF would be required, consistent with observations out
to z = 1.5 (Bezanson et al. 2011). We use the local SVDF
from Chae (2010) and converted it to a rough guess for the
shape of the SVDF at z ≈ 6, testing two extreme evolution
scenarios: (1) no effective evolution in the SVDF, due to an
approximate balance between mass evolution and size evolution of galaxies, (2) no size evolution, leading to a decreasein the break of the SVDF, consistent with the mass function
evolution.
For the first case we find γ σ∗ < 0.17 at 1σ confidence,with a best fit of γ σ∗ = −1.44, compared to the local M• −σ∗relation of McConnell & Ma (2013), while the second casegives γ σ∗ < 1.76 at 1σ confidence, with a best fit of γ σ∗ =0.25. In Fig. 7, we show the data and the best-fitting results
for the M• − Mdyn and the M• − σ∗ relation.
While the uncertainties in the parameters are large, due
to both the small number statistics and the uncertainty inthe underlying distributions, the results are fully consistent
with no evolution in the M• −bulge relations. Thus, evenwhile these z ∼ 6 QSOs show large apparent offsets fromthe local relation, this cannot be taken as evidence for intrinsic evolution in the M• −bulge relations up to high z.
The apparent offset is consistent with being fully caused by
AGN sample selection effects. Even negative evolution in theM• −bulge relations at z ∼ 6 is fully consistent with these
data.
Besides current methodical challenges in determining
black hole masses and galaxy properties at z ∼ 6, the significance on the M• −bulge relation at high z is currently
limited by small sample sizes and their bright flux limits.
Ongoing and future surveys, like UKIDSS, Pan-STARRS,VISTA, SuMIRe/Hyper-Suprime Cam, LSST will significantly increase the number of known z ∼ 6 QSOs andfuture sub-mm observations of their host galaxies, in particular using the Atacama Large Millimeter/submillimetre
Array (ALMA), will open a new window to the study of
their host galaxies and are able to obtain dynamical masses
from molecular line emission for a reasonably large sampleof high-z QSOs (Wang et al. 2013; Willott et al. 2013).
4.3

Current constraints on evolution of theM• −bulge relation

While we here do not aim at a full investigation of theliterature on the evolution of the M• −bulge relation, wewant to summarize the current constraints we have on evolution or non-evolution. We base this on samples that are
representative for the best data at their respective redshift, in terms of depth and well defined sample selection. We include the samples from Merloni et al. (2010) and
Wang et al. (2010) discussed above. Additionally, we augmented this by the samples from Cisternas et al. (2011) atz ∼ 0.5 and Targett et al. (2012) at z ∼ 4.

Cisternas et al. (2011) studied a sample of 32 AGN at
0.3 < z < 0.9, selected from XMM -COSMOS and report basically no offset from the local relation from Häring & Rix(2004) for their sample. Using the XMM -COSMOS soft Xray flux limit, the stellar mass function from Fontana et al.
(2006), and assuming a constant active fraction, we find an
insignificant bias for this sample. Thus, our inferred intrinsic relation is also fully consistent with the local relationfrom Häring & Rix (2004). However, compared to the local
relation from McConnell & Ma (2013) the sample is offset
by −0.25 dex. This might indicate a systematic offset in
the normalization of the virial black hole mass estimates,for example due to selection effects (SW11). In the following we use the M• −bulge relation from Häring & Rix (2004)
as local zero-point as a conservative constraint on possible
positive evolution.
Targett et al. (2012) could resolve the host galaxies for
two luminous QSOs at z = 4.16, selected from the SDSS,
with an additional selection criteria Mi < −28 mag. To account for the selection effects, we applied this absolute magnitude limit and assume the mass function for star-forming
galaxies at z ∼ 4 from Lee et al. (2012) and a constant active
fraction. The best-fitting result for the forward regression isγ = 0.21+1.36−2.50 .
In Fig. 8 we summarize our study of the intrinsic evolution in the M• −bulge relation for these four samples. The
apparent offset from the local relation is shown by the open
symbols, while our inferred intrinsic offset is given by thefilled symbols. While there is an apparent evolution trend,
accounting for selection effects we find an intrinsic offset
that is fully consistent with a non-evolving M• −bulge relation out to z ∼ 6. However, we note that at z > 2 theuncertainties on the M• −bulge relation are still substantial
and future work in particular in this redshift regime is necessary to obtain robust results. For reference we also indicate
the case of γ = 0.5, predicted by some numerical simulations (Di Matteo et al. 2008; Booth & Schaye 2011), by thedashed dotted line in Fig. 8. Current data are not able torobustly distinguish between a γ = 0.5 scenario and a nonevolving M• −bulge relation.

5

CONCLUSIONS

While the redshift evolution of the M• −bulge relation
provides important constraints on black hole-galaxy coevolution, its observational study is hampered by severe selection effects.
We presented a novel fitting method to the M• −bulge
relations and their evolution. This method corrects for thebias on the M• −bulge relation caused by sample selection
and is able to recover the intrinsic relation in the fitting. It
does so by incorporating the selection function of the observational sample and the underlying distribution function
into the conditional probability distribution of black holemass at given bulge property or vice versa. This probability
distribution is derived from the bivariate distribution function of M• and bulge property. The intrinsic relation is reconstructed via a maximum likelihood fit of this conditional
probability distribution to the data.
Since we are not modeling or fitting the full bivariate
distribution function but only the conditional probability

Infering the intrinsic BH-bulge relation

offset in log

M

•

1.5

1.0

0.5

0.0

-0.5

0

1

2

3

z

4

5

6

Figure 8. Redshift evolution of the inferred intrinsic M• −bulge
relation (filled symbols), compared to the apparent relation (opensymbols). We here studied the samples by Cisternas et al. (2011)at < z >= 0.56 (blue diamond), Merloni et al. (2010) at < z >=1.5 (red circle), Targett et al. (2012) at < z >= 4.2 (green triangle) and Wang et al. (2010) at z ∼ 6 (magenta square). The error
bars indicate the statistical uncertainty of the best fit, including
uncertainties in the underlying distributions. The dashed dottedline indicates a evolution model with γ = 0.5.

distribution, we can reduce the set of prior knowledge orassumptions to correct for the selection bias, compared toour discussion in SW11. We presented two routes to fit theconditional probability: (1) using the probability of findingM• at a given bulge property (forward regression), and (2)using the probability of finding a specific value for the bulge
property at a given M• (inverse regression). For the firstcase it is essential to properly understand the selection of
the sample and have a fair estimate of the black hole massdependence of the active fraction. For the second case thebulge distribution function and the intrinsic scatter of theM• −bulge relation need to be known. When measurement
errors on M• and the bulge property cannot be ignored,
all of these prior assumptions have to be known for a fullcorrection.
We extensively tested the fitting method on Monte
Carlo simulations and verified the reliability of the approachto recover evolution in the normalization unbiased. We also
illustrated the method by applying it to literature data fromz ∼ 0.5 to z ∼ 6. The sample of Merloni et al. (2010) atz ∼ 1.5 already puts significant constraints on the redshift evolution of the M• − M∗ relation. Assuming the local relation from Häring & Rix (2004), we can constrain theevolution in the M• − M∗ relation to γ < 0.9 out to redshift z ∼ 2 at 3σ confidence. This is in general agreementwith our conclusions in SW11 and with recent results atlower redshifts (e.g. Cisternas et al. 2011; Zhang et al. 2012;
Salviander & Shields 2013; Schramm & Silverman 2013).
Taking these empirical clues together, it becomes apparent
that there is currently no evidence for significant positive
evolution in the M• −bulge relation out to z ∼ 2. Either
black holes and galaxies grow almost co-eval since z ∼ 2, orgalaxies grow faster and the black holes would have to catchup to end up at the local relation at z = 0. The latter is11also suggested by observations of z ∼ 2 SMGs (Borys et al.
2005; Alexander et al. 2008; Carrera et al. 2011). Within
0 < z < 2 most of the black hole growth is governed by
secular processes and not by major mergers (Cisternas et al.
2011b; Schawinski et al. 2011; Kocevski et al. 2012). These
processes probably grow the black hole at the same rate
as the spheroid component. This changes at z > 2, beyond the peak of AGN activity, where major mergers
likely become the dominating black hole growth mechanism(Draper & Ballantyne 2012). Therefore also a change in theevolution of the M• − M∗ relation could be expected atz > 2. Currently, this redshift regime is still poorly covered by observations, and usually restricted to the most luminous, and hence most biased, QSOs. We illustrated thepoor constraints that can be drawn from current observations, employing a z ∼ 6 QSO sample based on Wang et al.
(2010). While the sample shows a large apparent positive
offset from the local relation, we find here that there is no
intrinsic offset, i.e. it is fully consistent with no or even negative evolution when sample selection is taken into account.
We therefore conclude that even at z > 2 and in particular
at z ∼ 6 we currently have no statistically significant evidence for evolution in the M• − M∗ relation. We confirmed
these results by additionally studying the low redshift sample from Cisternas et al. (2011) and the small z ∼ 4 sample
from Targett et al. (2012).
Future studies will reduce the observational uncertainties and increase the sample size and will allow a better
assessment of the M• − M∗ relation at these early cosmic
epochs. We will discuss the requirements on the design ofthese future studies of the M• − M∗ relation in this redshift
regime in a companion paper.
While selection effects severely complicate the observational study of the M• −bulge relation at high redshift,these can be included in the analysis as we have demonstrated here. We anticipate that the method presented here
will contribute to an improved understanding of black hole galaxy co-evolution. However, essential prerequisites are (1)the use of a well defined sample, i.e. the selection function forthe sample has to be under control, and (2) knowledge of theunderlying demographics of galaxies and active black holesover the black hole mass and galaxy property range coveredby the observations. Given the expected future increase in
the size of well defined samples and a better understanding of the underlying black hole and galaxy demographics,
we anticipate that we will be able to trace back the cosmic
history of the M• −bulge relations in the future.

ACKNOWLEDGEMENTS
We thank Knud Jahnke for fruitful discussions on thiswork. We also thank the anonymous referee for an efficientand constructive report that strengthened our conclusions.
We acknowledge support by the Deutsche Forschungsgemeinschaft under its priority programme SPP1177, grant
Wi 1369/23. AS acknowledges support by the China Postdoctoral Scientific Foundation, grant 2012M510256.

12

A. Schulze and L. Wisotzki

REFERENCES
Aird, J., Coil, A. L., Moustakas, J., et al. 2012, ApJ, 746,
90
Alexander, D. M., Brandt, W. N., Smail, I., et al. 2008, AJ,
135, 1968
Anglés-Alcázar, D., Özel, F., & Davé, R. 2013, ApJ, 770, 5
Bell, E. F., McIntosh, D. H., Katz, N., & Weinberg, M. D.
2003, ApJS, 149, 289
Bennert, V. N., Treu, T., Woo, J., et al. 2010, ApJ, 708,
1507
Bezanson, R., van Dokkum, P. G., Franx, M., et al. 2011,
ApJ, 737, L31
Bongiorno, A., Merloni, A., Brusa, M., et al. 2012, MNRAS, 427, 3103
Booth, C. M. & Schaye, J. 2011, MNRAS, 413, 1158
Borys, C., Smail, I., Chapman, S. C., et al. 2005, ApJ, 635,
853
Canalizo, G., Wold, M., Hiner, K. D., et al. 2012, ApJ, 760,
38
Carrera, F. J., Page, M. J., Stevens, J. A., et al. 2011,
MNRAS, 413, 2791
Chae, K.-H. 2010, MNRAS, 402, 2031
Cisternas, M., Jahnke, K., Bongiorno, A., et al. 2011, ApJ,
741, L11
Cisternas, M., Jahnke, K., Inskip, K. J., et al. 2011, ApJ,
726, 57
Croton, D. J. 2006, MNRAS, 369, 1808
Daddi, E., Renzini, A., Pirzkal, N., et al. 2005, ApJ, 626,
680
Decarli, R., Falomo, R., Treves, A., et al. 2010, MNRAS,
402, 2453
Di Matteo, T., Colberg, J., Springel, V., Hernquist, L., &Sijacki,D. 2008, ApJ, 676, 33
Di Matteo, T., Springel, V., & Hernquist, L. 2005, Nature,433, 604
Draper, A. R., & Ballantyne, D. R. 2012, ApJ, 751, 72
Dubois, Y., Devriendt, J., Slyz, A., & Teyssier, R. 2012,
MNRAS, 420, 2662
Fan, X., Narayanan, V. K., Lupton, R. H., et al. 2001, AJ,
122, 2833
Fan, X., Strauss, M. A., Richards, G. T., et al. 2006, AJ,
131, 1203
Fan, X., Strauss, M. A., Schneider, D. P., et al. 2003, AJ,
125, 1649
Fan, X., White, R. L., Davis, M., et al. 2000, AJ, 120, 1167
Ferrarese, L. & Merritt, D. 2000, ApJ, 539, L9
Fontana, A., Salimbeni, S., Grazian, A., et al. 2006, A&A,
459, 745
Gebhardt, K., Bender, R., Bower, G., et al. 2000, ApJ, 539,
L13
Gebhardt, K., Kormendy, J., Ho, L. C., et al. 2000, ApJ,
543, L5
Graham, A. W., Onken, C. A., Athanassoula, E., &Combes, F. 2011, MNRAS, 412, 2211
Gültekin, K., Richstone, D. O., Gebhardt, K., et al. 2009,
ApJ, 698, 198
Häring, N. & Rix, H.-W. 2004, ApJ, 604, L89
Hiner, K. D., Canalizo, G., Wold, M., Brotherton, M. S.,& Cales, S. L. 2012, ApJ, 756, 162
Ho, L. C. 2007, ApJ, 669, 821
Hopkins, P. F., Robertson, B., Krause, E., Hernquist, L.,& Cox, T. J. 2006, ApJ, 652, 107
Jahnke, K., Bongiorno, A., Brusa, M., et al. 2009a, ApJ,
706, L215
Jahnke, K., Elbaz, D., Pantin, E., et al. 2009b, ApJ, 700,
1820
Jahnke, K. & Macciò, A. V. 2011, ApJ, 734, 92
Kauffmann, G. & Haehnelt, M. 2000, MNRAS, 311, 576
Kelly, B. C., & Shen, Y. 2013, ApJ, 764, 45
Kocevski, D. D., Faber, S. M., Mozena, M., et al. 2012,
ApJ, 744, 148
Kormendy, J., & Ho, L. C. 2013, arXiv:1304.7762
Lamastra, A., Menci, N., Maiolino, R., Fiore, F., & Merloni, A. 2010, MNRAS, 405, 29
Lauer, T. R., Tremaine, S., Richstone, D., & Faber, S. M.
2007, ApJ, 670, 249
Lee, K.-S., Ferguson, H. C., Wiklind, T., et al. 2012, ApJ,
752, 66
Magain, P., Letawe, G., Courbin, F., et al. 2005, Nature,
437, 381
Magorrian, J., Tremaine, S., Richstone, D., et al. 1998, AJ,
115, 2285
Maiolino, R., Neri, R., Beelen, A., et al. 2007, A&A, 472,
L33
Marconi, A. & Hunt, L. K. 2003, ApJ, 589, L21
Marconi, A., Risaliti, G., Gilli, R., et al. 2004, MNRAS,
351, 169
McConnell, N. J., & Ma, C.-P. 2013, ApJ, 764, 184
McLure, R. J. & Jarvis, M. J. 2002, MNRAS, 337, 109
McLure, R. J., Jarvis, M. J., Targett, T. A., Dunlop, J. S.,& Best, P. N. 2006, MNRAS, 368, 1395
Merloni, A., Rudnick, G., & Di Matteo, T. 2004, MNRAS,
354, L37
Merloni, A., Bongiorno, A., Bolzonella, M., et al. 2010,
ApJ, 708, 137
Merloni, A. & Heinz, S. 2008, MNRAS, 388, 1011
Nesvadba, N. P. H., De Breuck, C., Lehnert, M. D., et al.
2011, A&A, 525, A43+
Nobuta, K., Akiyama, M., Ueda, Y., et al. 2012, ApJ, 761,
143
Park, D., Kelly, B. C., Woo, J.-H., & Treu, T. 2012, ApJS,
203, 6
Peng, C. Y. 2007, ApJ, 671, 1098
Peng, C. Y., Impey, C. D., Rix, H., et al. 2006, ApJ, 649,
616
Portinari, L., Kotilainen, J., Falomo, R., & Decarli, R.
2012, MNRAS, 420, 732
Richards, G. T., Lacy, M., Storrie-Lombardi, L. J., et al.
2006, ApJS, 166, 470
Richards, G. T., Strauss, M. A., Fan, X., et al. 2006, AJ,
131, 2766
Riechers, D. A., Walter, F., Brewer, B. J., et al. 2008, ApJ,
686, 851
Robertson, B., Hernquist, L., Cox, T. J., et al. 2006, ApJ,
641, 90
Salviander, S., Shields, G. A., Gebhardt, K., & Bonning,
E. W. 2007, ApJ, 662, 131
Salviander, S., & Shields, G. A. 2013, ApJ, 764, 80
Schawinski, K., Treister, E., Urry, C. M., et al. 2011, ApJ,
727, L31
Schramm, M., Wisotzki, L., & Jahnke, K. 2008, A&A, 478,
311
Schramm, M., & Silverman, J. D. 2013, ApJ, 767, 13

Infering the intrinsic BH-bulge relation
Schulze, A., & Gebhardt, K. 2011, ApJ, 729, 21
Schulze, A. & Wisotzki, L. 2010, A&A, 516, A87+
Schulze, A. & Wisotzki, L. 2011, A&A, 535, A87
Shankar, F., Bernardi, M., & Haiman, Z. 2009, ApJ, 694,
867
Shen, Y. & Kelly, B. C. 2010, ApJ, 713, 41
Sheth, R. K., Bernardi, M., Schechter, P. L., et al. 2003,
ApJ, 594, 225
Shields, G. A., Gebhardt, K., Salviander, S., et al. 2003,
ApJ, 583, 124
Sijacki, D., Springel, V., Di Matteo, T., & Hernquist, L.
2007, MNRAS, 380, 877
Silk, J. & Rees, M. J. 1998, A&A, 331, L1Somerville, R. S., Hopkins, P. F., Cox, T. J., Robertson,
B. E., & Hernquist, L. 2008, MNRAS, 391, 481
Targett, T. A., Dunlop, J. S., & McLure, R. J. 2012, MNRAS, 420, 3621
Tremaine, S., Gebhardt, K., Bender, R., et al. 2002, ApJ,
574, 740
Trujillo, I., Conselice, C. J., Bundy, K., et al. 2007, MNRAS, 382, 109
van Dokkum, P. G., Franx, M., Kriek, M., et al. 2008, ApJ,
677, L5
Vestergaard, M. & Peterson, B. M. 2006, ApJ, 641, 689
Volonteri, M., & Stark, D. P. 2011, MNRAS, 417, 2085
Walter, F., Carilli, C., Bertoldi, F., et al. 2004, ApJ, 615,
L17
Wang, R., Carilli, C. L., Neri, R., et al. 2010, ApJ, 714, 699
Wang, R., Wagg, J., Carilli, C. L., et al. 2013, ApJ, 773,
44
Willott, C. J., Omont, A., & Bergeron, J. 2013, ApJ, 770,
13
Woo, J.-H., Treu, T., Malkan, M. A., & Blandford, R. D.
2008, ApJ, 681, 925
Woo, J.-H., Schulze, A., Park, D., et al. 2013, ApJ, 772, 49
Zhang, X., Lu, Y., & Yu, Q. 2012, ApJ, 761, 5

APPENDIX A: ACCOUNTING FOR UPPER
LIMITSFor AGN samples, per definition the black hole is always
detected and its mass is estimated via the virial method.
However, the host galaxy is not always detected, as it canbe outshined by the AGN. But, if at least upper limits on
the host galaxy property can be set, these objects can be
included in the fitting. Otherwise their exclusion needs tobe reflected in the selection function. We include upper limits on x in the following manner, similar to Gültekin et al.
(2009)1 .
First, we compute the probability U (su , µ) that a blackhole with mass µ is located in a galaxy with their respective
property greater than su . For the forward regression without
measurement errors we simply assume a uniform distribution in s,
R∞RΩ pac (µ, z) g(µ | s, z) dssp(µ | su , λ , z) = Ru.
(A1)Ω pac (µ, z) g(µ | s, z) dµ1

We assume that every AGN is harboured by a galaxy, i.e.
there are no “naked“ quasars (See e.g. Magain et al. (2005);
Jahnke et al. (2009b) on this issue)

13

For the inverse regression without measurement errors wegetR∞g(µ | s) Φs (s) dssp(su | µ, z) = Ru.
(A2)g(µ | s) Φs (s)ds

When we include measurement errors in s and µ, the probability distribution for the forward regression is
R∞RΩ pac (µ, z) g(µ | s, z) Φs (s) g(µo | µ) dsdµsp(µo | so,u , λ , z) = u RR.
Ω pac (µ, z) g(µ | s, z) Φs (s) dsdµ(A3)
The likelihood for an upper limit at 1σ confidence level isli = δ1 U (su , µ) + (1 − δ1 )(1 − U (su , µ)) ,(A4)where U (su , µ) is one of the above conditional probabilitiesand δ1 = 0.16 is the probability to find the galaxy having s
greater than su . For an upper limit at a different confidence
level this probability changes accordingly.
This paper has been typeset from a TEX/ LATEX file preparedby the author.