Do you see what I see? Exploring the consequences of luminosity limits in black hole-galaxy evolution studies

TL;DR

This study uses simulations with selection effects to understand how luminosity limits influence the observed properties of AGN populations and their host galaxies, revealing biases introduced by different selection criteria.

Contribution

It introduces a simple galaxy evolution simulation incorporating a new AGN activity prescription and models selection effects to better interpret observed AGN populations.

Findings

Luminosity threshold selection reproduces the observed AGN luminosity function.

Eddington ratio threshold captures a broader range of black hole growth.

Different selection methods yield distinct host galaxy property distributions.

Abstract

In studies of the connection between active galactic nuclei (AGN) and their host galaxies there is widespread disagreement on some key aspects stemming largely from a lack of understanding of the nature of the full underlying AGN population. Recent attempts to probe this connection utilize both observations and simulations to correct for a missed population, but presently are limited by intrinsic biases and complicated models. We take a simple simulation for galaxy evolution and add a new prescription for AGN activity to connect galaxy growth to dark matter halo properties and AGN activity to star formation. We explicitly model selection effects to produce an "observed" AGN population for comparison with observations and empirically motivated models of the local universe. This allows us to bypass the difficulties inherent in many models which attempt to infer the AGN population by inverting selection effects. We investigate the impact of selecting AGN based on thresholds in luminosity or Eddington ratio on the "observed" AGN population. By limiting our model AGN sample in luminosity, we are able to recreate the observed local AGN luminosity function and specific star formation-stellar mass distribution, and show that using an Eddington ratio threshold introduces less bias into the sample by selecting the full range of growing black holes, despite the challenge of selecting low mass black holes. We find that selecting AGN using these various thresholds yield samples with different AGN host galaxy properties.

Equations

$$M_{\ast }/M={\mathcal{E}}_{M_{vir}}\exp \left(-{\left(\frac{\Delta M_{vir}}{{\sigma }_{M_{vir}}}\right)}^2\right),$$

$$\frac{dt}{d\log {\text{L}}_{bol}}={\left(\frac{{\text{L}}_{bol}}{{\text{L}}_{Edd}}\right)}^{-\alpha }\exp \left(-{\text{L}}_{bol}/{\text{L}}_{\mathrm{Edd}}\right).$$

$${\text{L}}_{x,LMXB}(z)$$

$${\text{L}}_{x,HMXB}(z)$$

$$\log \text{L}/{\text{L}}_{[2-10keV]}=(m\pm dm){\text{L}}_{Edd}+(q\pm dq),$$

Full text

1. Introduction

In the past decade there has been significant progress in building a generalized model of the formation and evolution of galaxies and their host halos across cosmic time (See Silk & Mamon 2012 for a review). These galaxy formation models have illustrated the impact of stellar and AGN feedback on galaxy growth (e.g., Kauffmann et al. 1993; Springel et al. 2005a; Bower et al. 2006; Croton et al. 2006; Genel et al. 2014; Schaye et al. 2015; Volonteri et al. 2015; Khandai et al. 2015; Feng et al. 2016). However, the connection between galaxy growth and black hole growth is not well understood. The volume-averaged galaxy-black hole growth rate is consistent with the black hole-spheroid mass relationship (e.g., Heckman et al. 2004) and the masses of supermassive black holes (SMBH) are found to be correlated with host stellar bulge properties (Kormendy & Ho, 2013), however the physical processes connecting black holes to their galaxies are still uncertain (Alexander & Hickox, 2012).

One major challenge lies in understanding selection effects (e.g., color selection, obscuration, flux limits, and luminosity limits) in observed AGN samples. Selection effects become apparent as different methods across wavelengths intended to probe the same physical phenomena yield populations with different AGN and galaxy properties. It is difficult to correct for selection effects as it is not immediately apparent which effects, or even how many, are influencing what is observed. For example, the distribution of the Eddington ratio (the ratio of bolometric luminosity to the Eddington limit) in the X-rays is observed to be approximately power-law in shape and spans several orders of magnitude, while in the optical, the Eddington ratio distribution is observed to vary by galaxy age (Kauffmann & Heckman, 2009). This discrepancy between observed accretion rate distributions, where accretion rate refers to the specific accretion rate, i.e. Eddington ratio, can be resolved when selection effects in spectroscopic classification are taken into account (Jones et al., 2016). Additionally, AGN are preferentially found in the most massive, bulge dominated galaxies across all wavelengths. X-ray selected AGN, in particular, are observed to reside in more massive host dark matter halos compared to optical AGN (e.g., Hickox et al. 2009; Koutoulidis et al. 2013; Richardson et al. 2013). Furthermore, the probability of finding an AGN is tied to the host galaxy star formation rate: the likelihood of finding an AGN increases for star forming galaxies compared to quiescent galaxies for a given accretion rate (e.g., Azadi et al. 2015). Yet the observation that AGN preferentially reside in massive galaxies is likely an effect due to AGN selection above a certain luminosity threshold. For a given luminosity threshold, it is possible to probe lower in Eddington ratio for high mass AGN compared to low mass AGN, for which only the highest accreting systems would be selected (Aird et al., 2012). The strength of AGN clustering measurements are also shown to vary across wavelengths (e.g., Hickox et al. 2009; Mendez et al. 2016).

In recent years, there has been increased insight into the physical origin of the black hole accretion rate distribution. This distribution is expected to cover a wide dynamic range and has converged on a broad universal shape that may be dependent on star formation (e.g., Aird et al. 2012; Bongiorno et al. 2012; Chen et al. 2013; Hickox et al. 2014; Azadi et al. 2015). Further evidence for this broad distribution is presented by Veale et al. (2014) in which a variety of accretion distributions are used to recreate the observed quasar luminosity function (QLF). An individual AGN may vary on short timescales within the dynamic range of this AGN accretion rate distribution (e.g., Hickox et al. 2014; Schawinski et al. 2015), while galaxies vary on longer timescales, on order of $\gtrsim 100$ Myr (e.g., Alexander & Hickox 2012). This variability may be due to accretion disk instabilities, or feedback on the accreting material (e.g., Siemiginowska & Elvis 1997; Hopkins et al. 2005; Janiuk & Czerny 2011; Novak et al. 2011), however, most observational evidence of this variability is indirect (e.g., Schawinski et al. 2010; Keel et al. 2012; Schirmer et al. 2013; Sartori et al. 2016). Due to this variability (or “flickering”) each observation is a snapshot of that source’s current AGN accretion: an instantaneous value, rather than a description of the average black hole activity.

Further attempts have been made to understand the interplay between AGN and their host galaxies through theoretical modeling. Hydrodynamical models, in particular, investigate AGN accretion and explicitly spatially resolve feedback on the interior halo of the host galaxy. Recent models, such as eagle (Schaye et al., 2015), illustris (Vogelsberger et al. 2014; Genel et al. 2014), massiveblack-II (Khandai et al., 2015), and BlueTides (Feng et al., 2016) are built to simultaneously describe complex physical processes and dark matter halo growth providing spatial and phase information, although this usually involves sub-grid modeling and many free parameters. These hydrodynamical simulations have had success in modeling a variety of AGN properties (e.g., winds, gas properties, luminosity functions), but number statistics are sacrificed for better resolution due to the high computational cost. Alternatively, semi-analytic models separate baryonic physics from dark matter halo growth and sacrifice spatial information by averaging over the spatial scale which makes them less computationally expensive than hydrodynamic models. This also usually involves utilizing many free parameters that are often degenerate, such that model predictions are not unique and can be difficult to interpret (e.g., Henriques et al. 2009; Lu et al. 2011; Mutch et al. 2013a).

In order to investigate the impact of selection effects on AGN observations, we have built a simple semi-numerical model for galaxy evolution and AGN accretion based on the observed connections between AGN and their host galaxies. In this work we take a “forward” modeling approach in which we simulate an intrinsic full galaxy population with a complete knowledge of its physical properties and apply known limits due to selection to compare to observations, rather than modeling unknown selection effects in observations. We combine the best characteristics of a semi-analytic simulation, a fast model with a low number of parameters, and the expository power of a straightforward, intuitive prescription for galaxy and black hole growth. We discuss our method for forward modeling the local AGN population in Section 2. This is broken into two subsections that discuss the galaxy formation model (2.1) and AGN prescription (2.2). In this paper we focus on the most simple of selection effects, namely luminosity limits that are primarily driven by the sensitivity limits of X-ray surveys. We compare our simulated AGN population to observations by imposing limits on luminosity and Eddington ratio in Section 3. A discussion of our results and summary are given in Section 4.

We utilize the same assumed cosmology as in Mutch et al. (2013b), in which a 1-year Wilkinson Microwave Anisotropy Probe (WMAP1; Spergel et al. 2003) cold dark matter (CDM) cosmology with $\Omega_{\text{m}}=0.25$, $\Omega_{\Lambda}=0.75$, and $\Omega_{\text{b}}=0.045$ is used. Likewise, all results are shown with a Hubble constant of $h=0.7$, where $h\equiv H_{0}/100\;\text{km s}^{-1}\;\text{Mpc}^{-1}$.

2. Method of Simulation

2.1. The Semi-Numerical Galaxy Formation Model

The host galaxies of our model AGN sample are built from a simulation of galaxies and dark matter halos using the formation history model of Mutch et al. (2013b). This model begins with a sample of dark matter halos from the N-body Dark Matter Millennium Simulation (Springel et al., 2005b), where the baryonic content of these halos is determined by the halo growth rate and the cosmological fraction of baryons compared to dark matter. The Millennium Simulation uses $N=10^{10}$ particles in a volume of $714\times 714\times 714$ Mpc and follows galactic evolution from $z=127$ to $z=0$ in 64 time steps of $\sim 200-350$ Myr.

Mutch et al. (2013b) combines two simple functions to connect galactic growth to the formation history of the Millennium Simulation dark matter halos; a baryonic growth function and a physics function. The baryonic growth function regulates the availability of baryonic material to be used by stars, mapping the dark matter growth history to the star formation rate (SFR). The physics equation represents the combined effects of internal and external physical processes (e.g. shock heating, supernova feedback, AGN feedback, galaxy mergers, tidal stripping, etc.) on the star formation rate, acting as an efficiency of baryonic matter consumption.

In this model, the conversion efficiency of the baryons into stars is based on a simple function (Equation 1) that represents baryon cooling, star formation, and both stellar and AGN feedback physics as one net component, rather than treating them individually which can be computationally costly or can introduce significant uncertainty due to a large number of physical parameters:

[TABLE]

where $\mathcal{E}_{M_{vir}}$ represents the baryonic matter conversion efficiency and $M_{vir}$ is the halo virial mass with standard deviation $\sigma_{M_{vir}}$. Treating the baryonic matter in this way allows for each galaxy’s formation to follow the growth history of its individual halo, capturing the true diversity of galaxy formation histories, rather than adding an artificial scatter as is done with many other methods. This approach has proven to be accurate for $z<2$ in reproducing the derived star formation rates as a function of halo mass (Conroy & Wechsler 2009; Béthermin et al. 2012). The simplicity is also important and makes the model general and flexible such that the functional form can be used for either the virial mass of the halo (Mvir) or the instantaneous maximum circular velocity (Vmax), which is more directly tied to the gravitational potential. Both Mvir and Vmax provide good fits and can be used interchangeably, although Vmax may be more tightly coupled with the stellar mass growth (Reddick et al., 2013).

The Mutch et al. (2013b) model is able to accurately recreate the local stellar mass function and its evolution to $z\sim 3$. It also has the capability to trace the stellar mass function as a function of stellar age, which allows for passive and star forming galaxies to be treated separately. For more information on this model, please refer to Mutch et al. (2013b).

In the Mutch et al. (2013b) semi-numerical model, galaxies are assigned a star formation rate (SFR) following the growth of the dark matter halos and availability of the baryonic material. If the galaxy is not flagged as the central galaxy of a dark matter halo, it is assigned SFR$\;=0$. Since these leftover “satellite” galaxies are passive with low SFR, they may be treated with this simplification when investigating the local and global stellar mass function. However, in this work we are directly comparing to the observed SFR distribution and this simple assumption of assigning “satellite” galaxies a SFR$\;=0$ is no longer applicable. There are also a fraction $\sim 59\%$ of central galaxies with SFR$\;=0$ for which we calculate a SFR.

We have devised a simple prescription for repopulating galaxies with SFR$\;=0$ based on the distribution of specific star formation rate (sSFR; SFR/M*) of our simulated passive central galaxies. For central galaxies with SFR$\;=0$, the SFR are smoothed over three redshift bins in order to limit the number of central galaxies with SFR$\;=0$. Any galaxy with SFR$\;=0$, including central galaxies with SFR$\;=0$ after the smooth, are assigned a sSFR representative of a passive galaxy. The distribution of sSFR for central galaxies is well fit by two gaussians, so we assign an sSFR consistent with the passive gaussian distribution. Using this simple method, we can repopulate the SFR for our “satellite” and central galaxies in order to more accurately compare to observations. The repopulated sSFR drawn from the passive central galaxy sSFR distribution are low enough that they do not greatly impact the average doubling time for our simulated sample (e.g. for a stellar mass of $\sim 10^{10}\;M\odot$ the typical sSFR assigned is $10^{-11}\;\text{yr}^{-1}$, corresponding to a doubling time of $\sim 100$ Gyr, and as such, the corresponding increase in the stellar mass at each time step would be negligible).

2.2. A Simple Prescription for AGN Accretion

Since the Mutch et al. (2013b) galaxy evolution model synthesizes complicated physics down to a one-dimensional function, we are able to add complexity in the form of an AGN component and limit the number of free parameters, while remaining computationally tractable. The foundation of our simple prescription for AGN accretion is motivated by the work of Jones et al. (2016), in which it is assumed the instantaneous observed AGN luminosity is due to short-term variability (e.g., Alexander & Hickox 2012; Hickox et al. 2014) and follows an Eddington ratio ($\text{L}_{bol}/\text{L}_{Edd}$) distribution described by a Schechter function, a power law with an exponential cutoff near the Eddington Limit (Hopkins et al., 2009):

[TABLE]

An Eddington ratio distribution that takes the functional form of a Schechter function consists of only two free parameters: the slope of the power law $\alpha$, and the lower cut off $\text{L}_{cut}$ which sets the amplitude of the function such that the integral of the curve is one.

Our black hole masses are derived from the black hole-bulge relationship of Häring & Rix (2004), with total galaxy stellar mass from the model galaxies used as a proxy for bulge mass (Kormendy & Ho, 2013). We note that there is significant uncertainty and scatter in the relationship between black hole mass and total stellar mass, with marked differences observed for local ellipticals and AGN hosts (e.g., Reines & Volonteri 2015). We adopt the Häring & Rix (2004) relationship as it approximately bisects that found for ellipticals and AGN hosts, and broadly represents the correlation observed for the full galaxy population. Furthermore, we do not include a redshift evolution in our black hole mass calculations. There is widespread disagreement about whether a redshift evolution is present or if any observed evolution is caused by a sample selection bias (e.g. Decarli et al. 2010; Merloni et al. 2010; Schulze & Wisotzki 2011; Bennert et al. 2011; Woo et al. 2013; DeGraf et al. 2015; Shankar et al. 2016). We tested the impact of a moderate evolution (Merloni et al., 2010) on our analysis and found that our simulated luminosity functions for both an evolving and non-evolving black hole-galaxy relationship were consistent within their random scatter for redshifts up to $z\sim 3.5$.

For simplicity, we assume a constant black hole-galaxy relationship. We can thus calculate the Eddington Luminosity, $\text{L}_{Edd}=1.38\times 10^{38}\;\text{M}_{BH}$. We then select an instantaneous AGN bolometric luminosity from the Eddington ratio distribution. Both the galaxy and AGN components are defined in terms of bolometric luminosities; using scaling relationships these can be converted into different broad-band luminosities. In this paper we focus on the hard ($2-10$ keV) X-rays, although we note that this approach could be equally valuable in the $>10\;\text{keV}$, soft X-rays, infrared, and/or optical regime.

Galactic X-ray emission is produced by a combination of high mass X-ray binaries (HMXB), low mass X-ray binaries (LMXB), and hot gas (e.g., Hornschemeier et al. 2005; Lehmer et al. 2010; Fragos et al. 2013). In this work we ignore hot gas since X-ray binaries dominate the X-ray emission above $1.5\;\text{keV}$ (Lehmer et al., 2016). We use the Lehmer et al. (2016) scaling relationships given by

[TABLE]

where $\log\alpha_{0}=29.30\pm 0.28$, $\gamma=2.19\pm 0.99$, $\log\beta_{0}=39.40\pm 0.08$, $\delta=1.02\pm 0.22$, and the contribution from HMXB scales with star formation rate (SFR) while LMXB scales with stellar mass.

The AGN X-ray luminosity is derived from the bolometric AGN luminosity using the Lusso et al. (2012) observed relationship for Type 1 AGN,

[TABLE]

where $m\pm dm=0.752\pm 0.035$ and $q\pm dq=2.134\pm 0.039$. We use the function for Type 1 AGN to model the implicit relation for our full sample. Rather than utilizing the Lusso et al. (2012) Type 2 AGN relationship for an “obscured” percentage of our full sample, we add in obscuration using a different method to lessen the influence of potential observational biases incurred in observing Type 2 AGN.

We randomly assign some of our model AGN to be obscured based on the relationship between LX and the obscured fraction from Merloni et al. 2014 (See their figure 9). For AGN that are “unobscured”, we assign a column density of $\log N_{H}=20\;\text{ cm}^{-2}$, while those selected to be “obscured” have column densities drawn from the NuSTAR-informed $N_{H}$ distribution (Lansbury et al. 2015; Figure 13b). Based on these column densities and X-ray spectral models (e.g., Lansbury et al. 2014) we then add X-ray absorption to our obscured sample. We run our simulation both with and without this prescription for obscuration and find that adding obscuration does not significantly alter the results of our analysis, other than extending the lower cutoff of the AGN X-ray luminosities from $\text{L}_{X}=10^{37.5}\;\text{erg s}^{-1}$ to $\text{L}_{X}=10^{35.5}\;\text{erg s}^{-1}$ (with obscuration) and altering the amplitude of our AGN X-ray luminosity function (see Section 3.1). There is no change for the highest luminosities ($\text{L}_{X}>10^{44}\;\text{erg s}^{-1}$). We thus have a simulated galaxy and AGN X-ray luminosity for every object in our sample and can investigate the properties of this X-ray emission to compare to observations.

3. Comparison of the Simulated AGN Population with X-ray Observations

In our simulation we built bolometric AGN and galaxy luminosities which are then scaled into X-ray luminosities. At this point our X-ray sample represents an intrinsic snapshot of the full AGN population. In order to compare our simulation to observations we must forward model the selection effects we would expect to see in these observed results.

In the literature, there are a range of different parameters used to define an AGN. In the X-ray band, AGN are most often selected based on an observed luminosity threshold. This threshold is often set at the transition where AGN emission begins to dominate host galaxy emission ($\text{L}_{X}\sim 10^{41.5}\;\text{erg s}^{-1}$). Despite the ease of this method for selecting the brightest AGN, by definition it does not select less luminous AGN, especially when these AGN have luminosities comparable to their host galaxy (e.g., highly star forming galaxies).

Alternatively, an AGN may be defined by its accretion rate using colors, broad lines, or by the Eddington ratio. Typically, this threshold is set at an Eddington ratio of $\lambda\gtrsim 0.01$ (e.g., Hopkins et al. 2009). Since the Eddington ratio is closely linked to the black hole mass, selecting an AGN in this way often relies on accurate measurements of the black hole, which may be difficult to acquire for large samples. Additionally, it becomes challenging to separate AGN accretion from other accretion processes at low luminosities (e.g., Lehmer et al. 2016). An Eddington ratio limit, however, is better for selecting low-luminosity AGN compared to a luminosity limited sample, as demonstrated by our following analysis. In this work we are able to fully explore this threshold since each of our simulated AGN has a known Eddington ratio drawn from our Schechter function distribution.

3.1. The AGN Luminosity Function

We first compare our model AGN population to the observed AGN X-ray luminosity function (XLF). This distribution is particularly useful for investigating AGN evolution but observations are uncertain at faint luminosities and high redshifts due to poor number statistics and contamination (e.g., Nandra et al. 2005; Aird et al. 2010). To mitigate these problems, luminosity limits are commonly placed on samples, thereby adding a known selection effect.

3.1.1 The AGN luminosity function at $z\;=0$

We compare our model initially to the luminosity-dependent density evolution (LDDE) model of the hard X-ray luminosity function from Aird et al. (2010). We are able to “tune” our model to match the Aird et al. (2010) XLF at $z=0$ by varying our model AGN Eddington ratio distribution $\alpha$ and lower cutoff parameters. We approximately match the distribution with $\alpha=0.8$ and $\text{L}_{cut}=-6.0$ $\text{L}_{bol}/\text{L}_{Edd}$. This Schechter function distribution deviates from the parameters presented in Jones et al. (2016) ($\alpha=0.4$ and $\text{L}_{cut}=-3.75$ $\text{L}_{bol}/\text{L}_{Edd}$). We expect that alpha will be steeper since our distribution is describing AGN activity in both the active and quiescent galaxy population (Gabor & Bournaud, 2013), whereas in Jones et al. (2016), we focus on the Eddington ratio distribution of star forming galaxies alone. It is worthwhile to note that when comparing these theoretical distributions at $z=0$ with observations that are higher redshifts, an exact comparison is difficult due to the evolution of the XLF (e.g., Ueda et al. 2014).

The luminosity function that we calculate with our simulated AGN population appears as a double power law with its break at $\text{L}_{X}\sim 10^{44}\;\text{erg s}^{-1}$ (solid grey line, Figure 1). This characteristic shape is influenced by the black hole mass function and is consistent with what is modeled by Aird et al. (2010, orange, Figure 1). At low luminosities we observe a turnover due to the galaxy dominating the X-ray emission and the minimum contribution from star formation as imposed by our black hole mass limit. Since our simulation relies on the random selection of the accretion rate from the Eddington ratio distribution, as well as randomly assigning which galaxies are obscured, we ran a bootstrap resampling to demonstrate the variance of our simulation due to these random selections. We show these random assignments and limited volume do not introduce a significant uncertainty to our model XLF (grey area, Figure 1).

We can further examine the effects of selecting AGN based on two different thresholds on the “observed” XLF. The Eddington ratio limited sample (solid purple line, Figure 1) yields a similar wide distribution in luminosity ($39.5\lesssim\log\text{L}_{X}\lesssim 45.5$) compared to the full model galaxies (solid grey line, Figure 1). We note that in practice, it is difficult to select AGN below $\text{L}_{X}\sim 10^{41}\;\text{erg s}^{-1}$ due to contamination from other accretion processes (e.g., Lehmer et al. 2016). Imposing a luminosity limit (solid green line, Figure 1), as we would expect, causes the XLF to truncate below $\text{L}_{X}=10^{41.5}\;\text{erg s}^{-1}$ as it is defined, which is consistent with the LDDE model of Aird et al. (2010). At high luminosities ($\text{L}_{X}\gtrsim 10^{44}\;\text{erg s}^{-1}$), where all AGN are accreting at high Eddington rates, we find that our models, including those that have imposed limits, are consistent with the Aird et al. (2010) XLF at $z=0$. This suggests that our simple prescription for AGN accretion is valid to first order for describing the full AGN population at $z=0$.

3.1.2 An evolving AGN luminosity function

We further explore the evolution with redshift of the XLF in our model. AGN X-ray luminosities for galaxies in earlier (higher-z) snapshots from the Mutch et al. (2013b) simulation are calculated following the $z=0$ prescription outlined in Section 2.2 (also using an Eddington ratio distribution given by a Schechter function). We calculate the minimum chi-squared fit between our model XLF and the Aird et al. (2010) AGN XLF for a selection of nine redshift bins in order to determine our best fit parameters for the Schechter function $\alpha$, and lower cutoff (Figure 2). Over this redshift range, $\alpha$ varies between $0.3$ and $0.8$, while the lower cutoff varies between $-5.0$ and $-7.0$. We find that our XLF are qualitatively consistent with the Aird et al. (2010) XLF within these redshift bins to $z\sim 1.2$ for our best fit parameters. This confirms the applicability of the model at $z\lesssim 1$ for further comparisons to observed sSFR, mass, and halo occupation distribution results.

3.2. Host Galaxy SFR and Mass

The Mutch et al. (2013b) simulation tracks the star formation rate of our model galaxies, thus we can calculate the specific star formation rate (sSFR; ratio of the star formation rate to the stellar mass). The full simulated distribution in sSFR-stellar mass (grey contours, Figure 3) is consistent with the shape of the eagle hydrodynamical theoretical distribution for $z=0$ (Guo et al., 2016) as well as the observed sSFR-stellar mass distributions from Schiminovich et al. (2007) ($0.01<z<0.25$) and the GAMA survey (Bauer et al., 2013, $0.05<z<0.32$). Thus, with our simple prescription, we are able to match more complex models to first order.

By design, our Eddington ratios are tied directly to the stellar mass through our black hole mass calculations. Thus, applying our Eddington ratio cut to the full sample does not change the overall shape of the distribution in sSFR-stellar mass but simply limits the total number of AGN selected in the sample.

We now apply a luminosity limit to “define” our AGN and compare the resulting distributions to the X-ray observations from Mendez et al. (2016). Mendez et al. (2016) utilizes an X-ray luminosity threshold of $\text{L}_{X}=10^{41.5}\;\text{erg s}^{-1}$ to select AGN from $0.2<z<1.2$. In order to more accurately compare these results to our model output at $z=0$, we correct the sSFR of this observed comparison sample using the Lee et al. (2015) galaxy main sequence evolution ($\text{SFR}_{0}\propto(1+z)^{4.12\pm 0.08}$) at $z=0.2$. We find comparable “observed” distributions that are truncated at higher stellar mass compared to our intrinsic AGN population; we are able to match observations that are different from our intrinsic distribution by using appropriate limits. We see a slight shift of our luminosity limited model to higher stellar masses compared to the Mendez sample as shown by the histograms.

Our overestimation of the stellar masses may in part be due to the connection between AGN accretion and star formation, such that the AGN fraction decreases for passive galaxies (e.g., Chen et al. 2013; Azadi et al. 2015). AGN suppression in high mass systems is not yet included in this simple model.

With the addition of the galaxy main sequence evolution to the Mendez et al. (2016) sample, our limited model sSFR appears consistent with observations. Using a straightforward and simple semi-numerical model, we are able to build comparable intrinsic sSFR-stellar mass distributions, as well as broadly match what is observed with appropriate limits.

3.3. The Halo Occupation Distribution of AGN

We further investigate the host dark matter halos of our simulated AGN and compare our simple model to results from the observational clustering analysis of Richardson et al. (2013) as well as the simulation of Chatterjee et al. (2012). We compute the halo occupation distribution (HOD) for luminosity and Eddington ratio limited samples with respect to our intrinsic full sample. As shown in the color-mass analysis, our Eddington ratio limited sample (magenta, Figure 4) follows a distribution consistent with the intrinsic sample and thus the distribution of the ratio appears flat on the HOD. This suggests that by selecting AGN based on Eddington ratio we are capturing the black hole growth across all mass limits.

Our luminosity limited sample, however, deviates from our model intrinsic sample so the HOD takes on a positive slope with a turnover at $M\sim 10^{12}M_{\odot}$. We compare this shape to the Richardson et al. (2013) distribution at $z\sim 1.2$ and find a similar turnover, albeit at higher halo masses. The Chatterjee et al. (2012) HOD at $z=1.0$ similarly exhibits a positive slope. While these two HODs represent populations at varying halo mass and differing redshifts, they provide a general shape and range with which we may check the validity of our local universe luminosity limited sample distribution. We find that our model HOD with luminosity limits built from a simple semi-numerical model attains the general shape of a theoretically motivated HOD as well as the observationally informed HOD. This reinforces the result that a simple prescription can describe the observed AGN population to first order and by varying the criteria that define our AGN population, we are selecting AGN with different host galaxy and halo properties.

4. Summary

Our goal has been to build a simulation of the full AGN population and their host galaxy properties using a straightforward semi-numerical model, in order to investigate the impact of AGN selection criteria on the “observed” population. We begin with a sample of dark matter halos from the Millennium N-body simulation (Springel et al., 2005b) and connect galactic growth to the dark matter halo formation history using the Mutch et al. (2013b) formation history semi-numerical model. Our AGN accretion model is motivated by the results of Jones et al. (2016) and is straightforward and computationally inexpensive. It treats AGN accretion as an instantaneous rate selected from an Eddington ratio distribution that takes the functional form of a Schechter function. Using a variety of galaxy and AGN scaling relationships, we convert our simulated bolometric luminosities into X-ray luminosities to better compare to observations. We assign AGN to be obscured or unobscured following the prescription outlined in Section 2.2.

With our simulated full AGN population, we explicitly model selection effects by defining our “observed” AGN samples based on typical observational limits; a luminosity threshold of $\text{L}_{X}>10^{41.5}\;\text{erg s}^{-1}$ and an Eddington ratio limit of $\lambda\gtrsim 0.01$. We then compare our intrinsic AGN population, as well as the two limited samples, to a variety of published relationships. We are able to recreate the Aird et al. (2010) XLF using a Schechter function Eddington ratio distribution with a power-law slope of $\alpha=0.8$ and lower cutoff of $\text{L}_{cut}=-6.0$. The sSFR-M* distribution of our intrinsic AGN population is broadly consistent with that of the eagle hydrodynamical simulation at $z=0$ (Guo et al., 2016) as well as the observational distributions of Schiminovich et al. (2007) ($0.01<x<0.25$) and Bauer et al. (2013) ($0.05<z<0.32$), while our luminosity limited AGN population matches the observed X-ray AGN of Mendez et al. (2016). We further examine the consequences of these selection criteria by comparing our model HODs to the observational analysis from Richardson et al. (2013) and theoretical model of Chatterjee et al. (2012) and find that our luminosity limited sample yields the same general shape of the HOD, while our Eddington limited sample maintains the distribution of our intrinsic sample. We find that our straightforward prescription for AGN accretion is able to recreate both theoretical and observational relationships for the “observed” AGN population, once AGN selection criteria are taken into account.

AGN are often defined using a luminosity or Eddington threshold. We have shown that each criterion selects different “observed” AGN properties and host galaxy parameters from our intrinsic full AGN population. Selecting AGN based on an Eddington ratio limit yields a wider range of AGN and host galaxy properties compared to a luminosity limit and may better represent the underlying AGN population. A luminosity limited sample is often driven by the sensitivity of surveys; our luminosity limited model deviates most from the intrinsic full population. This selection effect can be partially mitigated by selecting AGN properties that are directly related to the AGN accretion rate. This work highlights the difficulty in correcting for even the most simple of selection effects in AGN X-ray observations.

In the future, this analysis can be extended to higher redshifts and tied to star formation using the full history of the Mutch et al. (2013b) model and our simple prescription for AGN accretion. Modeling the AGN with additional selection techniques such as luminosity thresholds at various wavelengths and photometric color selection will be particularly useful for comparison with survey data.

We thank our collaborators as well as the referee for constructive comments that improved the paper. This work was supported in part by the National Aeronautics and Space Administration under Grant Number NNX15AU32H issued through the NASA Education Minority University Research Education Project (MUREP) through the NASA Harriett G. Jenkins Graduate Fellowship activity. This research has made use of NASA’s Astrophysics Data System.