Abstract

The lack of a strong correlation between AGN X-ray luminosity (LX; a proxy for AGN power) and the star formation rate (SFR) of their host galaxies has recently been attributed to stochastic AGN variability. Studies using population synthesis models have incorporated this by assuming a broad, universal (i.e. does not depend on the host galaxy properties) probability distribution for AGN specific X-ray luminosities (i.e. the ratio of LX to host stellar mass; a common proxy for Eddington ratio). However, recent studies have demonstrated that this universal Eddington ratio distribution fails to reproduce the observed X-ray luminosity functions beyond z∼1.2. Furthermore, empirical studies have recently shown that the Eddington ratio distribution may instead depend upon host galaxy properties, such as SFR and/or stellar mass. To investigate this further we develop a population synthesis model in which the Eddington ratio distribution is different for star-forming and quiescent host galaxies. We show that, although this model is able to reproduce the observed X-ray luminosity functions out to z∼2, it fails to simultaneously reproduce the observed flat relationship between SFR and X-ray luminosity. We can solve this, however, by incorporating a mass dependency in the AGN Eddington ratio distribution for star-forming host galaxies. Overall, our models indicate that a relative suppression of low Eddington ratios (λEdd≲0.1) in lower mass galaxies (M∗≲1010\-−11M⊙) is required to reproduce both the observed X-ray luminosity functions and the observed flat SFR/X-ray relationship.

One means of connecting SMBH and galaxy growth that has gained popular support over the past two decades is for accreting SMBHs (observed as active galactic nuclei, AGNs) to directly influence the star formation rates of their hosts (hereafter SFRs) via a variety of feedback mechanisms (see the review of Fabian, 2012, for details on the feedback mechanisms). However, a major difficulty in identifying the precise role that SMBH accretion has on influencing SFR is the stochastic nature of AGNs (e.g. Aird
et al., 2013; Hickox et al., 2014, see also § 4.3 of Stanley
et al. 2015). It is argued that this randomness is the reason why most studies that have explored the relationship between SFR and X-ray luminosity (a proxy for SMBH accretion power) have reported no evidence of a strong correlation between these parameters, at least for moderate luminosity AGNs which form the majority of the population (i.e. 1042<LX<1045ergs−1; e.g. Lutz
et al., 2010; Harrison
et al., 2012; Mullaney
et al., 2012; Rosario
et al., 2012; Rovilos
et al., 2012; Santini
et al., 2012; Azadi
et al., 2015; Stanley
et al., 2015). Recently, using optically selected AGNs, Stanley
et al. (2017) reported an enhancement of SFR among the highest luminosity AGNs (i.e. LX>1045 erg s−1). However, they demonstrated that this is a direct consequence of the most luminous AGNs residing in more massive host galaxies, meaning the enhanced SFR is a consequence of the relationship between stellar mass and SFR (hereafter referred to as the main sequence, MS; e.g. Daddi
et al., 2007; Rodighiero
et al., 2011; Schreiber
et al., 2015).

A major problem in investigating the relationship between AGN and SFR is that observations are often subjected to biases (e.g. flux limited samples, rarity of high accretion rate AGNs). Therefore, many studies investigating the connection between SMBH and galaxy growth have adopted a modelling approach (e.g. Aird
et al., 2013; Conroy &
White, 2013; Veale
et al., 2014; Caplar
et al., 2015; Jones et al., 2017; Weigel
et al., 2017). In these models, the stochastic nature of AGNs is often implemented via a broad specific AGN X-ray luminosity (i.e. the X-ray luminosity relative to the host stellar mass; LX/M∗) distribution that, once combined with a mass function, reproduces the observed X-ray luminosity function of AGNs. Overall, these studies find that below z∼1 the specific X-ray luminosity distribution is well-represented by a broad, universal (i.e. the same distribution, irrespective of SFR or stellar mass) distribution, described by a broken power-law (or similar; e.g. a Schechter 1976 function), the normalisation of which increases with increasing redshifts (e.g. Aird
et al., 2013; Veale
et al., 2014; Jones et al., 2017). These models are also often successful at reproducing the observed flat relationship between SFR and X-ray luminosity (Hickox et al., 2014; Veale
et al., 2014; Stanley
et al., 2015). However, Jones et al. (2017) recently found that using a universal broad distribution for the specific X-ray luminosity distribution of AGNs cannot reproduce the X-ray luminosity functions of Aird
et al. (2010) beyond z∼1.2. They argue for the need of a more complicated specific X-ray luminosity distribution of AGNs beyond that redshift. Furthermore, recent studies have reported that the specific X-ray luminosity distribution for AGNs in star-forming hosts significantly differs from that of AGNs in quiescent hosts (e.g. Georgakakis
et al., 2014; Aird
et al., 2017a; Wang
et al., 2017). Overall, they report a higher contribution from star-forming AGN hosts to the total specific X-ray luminosity distribution out to z∼2.

Motivated by the above findings, we explore a model for the specific X-ray luminosity distribution split between star-forming and quiescent galaxies. After demonstrating that this simple model cannot simultaneously reproduce both the X-ray luminosity functions and the flat relationship between SFR and X-ray luminosity out to z∼2, we introduce a mass dependency for the specific X-ray luminosity distribution of star-forming AGN hosts, as recently observed (e.g. Bongiorno
et al., 2016; Aird
et al., 2017a; Georgakakis
et al., 2017). Hereafter, we refer to our model specific X-ray luminosity distribution split between star-forming and quiescent galaxies as “mass-independent”, and to our model that includes the mass dependency for the specific X-ray luminosity distribution of AGNs in star-forming galaxies as “mass-dependent”.

In § 2 we describe our method to infer the specific X-ray luminosity distribution for both the mass-independent and the mass-dependent model. Subsequently, in the same section, we present how we incorporate host galaxy properties. We then report the results for the mass-independent model in § 3.1, and for the mass-dependent model in § 3.2. We discuss the general implications of our results in § 4, and conclude in § 5. Throughout, we assume the specific X-ray luminosity (i.e. LX/M∗) to be proportional to Eddington ratio (i.e. the ratio between bolometric AGN luminosity and Eddington luminosity, LAGN/LEdd), converting the X-ray 2-10 keV luminosity into bolometric AGN luminosity using a universal conversion factor of 22.4 (median value found in Vasudevan &
Fabian, 2007, and based on a local AGN sample with LX=1041\-−46ergs−1), and assuming a ratio of 0.002 between SMBH masses and stellar masses (Marconi &
Hunt, 2003). Therefore, we have,

λEdd=22.4×LX1038.1%
erg s−1×0.002M∗M⊙,

(1)

where λEdd is the Eddington ratio, and LX is the intrinsic (absorption-corrected) 2-10 keV X-ray luminosity. We stress that the important parameter we investigate is the specific X-ray luminosity (LX/M∗), and that the Eddington ratio is solely used for its greater familiarity.

Many recent studies have used Population Synthesis Models (hereafter PSMs) to investigate the connection between AGNs and host galaxies. PSMs randomly draw properties from probability distribution functions (e.g. Aird
et al., 2013; Hickox et al., 2014; Stanley
et al., 2015; Jones et al., 2017). Using PSMs we are able to generate populations of N massive mock galaxies, with stellar masses M∗, drawn from a mass distribution, for which we randomly allocate Eddington ratios, λEdd, following an Eddington ratio distribution, and, using Eq. 1, derive X-ray luminosities, LX, the histogram of which is proportional to the model X-ray luminosity function. The key point is to optimise the Eddington ratio distribution until the model X-ray luminosity function matches the observed one (e.g. Aird
et al., 2013; Jones et al., 2017). This can be repeated at different redshifts to investigate any redshift dependence of the Eddington ratio distribution. Because it implies generating a new population of galaxies while optimising the Eddington ratio distribution, it can be computationally expensive. Therefore, we will use an analytical approach (similar to those of e.g. Veale
et al., 2014; Caplar
et al., 2015) by directly considering the various distributions as probability density functions. In this study, we fit to the X-ray luminosity functions of Aird et al. (2015) and use the stellar mass functions of Davidzon
et al. (2017). The latter are shown in Fig. 1 for various redshifts, out to z=2.5.

In our analytical approach, we first assume a model λEdd distribution which is described by a set of parameters θ (for example, for a power law θ={N,α}, where N is the normalisation and α the slope). Our aim is to iterate over θ to identify the solution that best fits the observed X-ray luminosity function. The model X-ray luminosity function is generated by analytically combining the observed mass function with the λEdd distribution given by the set of parameters θ. To achieve this, we treat the X-ray luminosity function as a probability distribution, p(LX). The probability of observing an AGN of luminosity LX=X in a galaxy of mass M∗ is given by1,

p(LX=X)≡p(λEdd=XM∗),

(2)

where p(λLedd=X/M∗) is the probability of λEdd=X/M∗ calculated with the current set of model parameters, θ. We derive the total probability of observing LX=X in all galaxies as the sum of these individual probabilities weighted by the galaxy mass function (which we also treat as a probability distribution, p(M∗)),

p(LX=X)=MMax∑Y=MMinp(λEdd=XY|M∗=Y)×p(M∗=Y),

(3)

where p(M∗=Y) is the probability of M∗=Y (or the mass function evaluated at Y). The total model luminosity function derived from a given set of parameters, θ, is then simply the combination of these probabilities (i.e. for LX=1041\-−1046 erg s−1 in bins of 0.1 dex).

By taking this approach we avoid having to model N galaxies (where N∼O(109) to generate sufficient numbers of luminous galaxies), and can simply split our mass function into O(103) bins between MMin=108M⊙ to MMax=1014M⊙ (where MMinandMMax are chosen to cover the full range of stellar masses). Since we adopt an iterative process to identify the best fitting parameters (see later in this subsection), this factor of O(106) reduction in the number of required calculations per iteration dramatically reduces the time taken to converge.

To optimise the parameters that define our Eddington ratio distribution, we performed maximum likelihood estimation. Assuming Gaussian distributions for the uncertainties, the log-likelihood of a set of parameters θ is defined by -0.5χ2, where χ2 is the chi-square between the model and the observed X-ray luminosity functions. We used flat proper (i.e. with finite boundaries) prior distributions for each parameter that defines our model, and checked the posterior distributions to ensure that they are not constrained in any way by the limits of our prior distributions. The parameter space (defined by the flat prior distributions) was explored using a Monte Carlo Markov Chain (MCMC) fully implemented in a python application programming interface, emcee2(Foreman-Mackey et al., 2013), that uses the affine invariant MCMC ensemble sampler of Goodman &
Weare (2010). We ran our MCMC code using 100 walkers for 10000 steps as a burn-in phase to find the global solution (i.e. these 10000 steps are then discarded), and for 5000 steps as a production phase (i.e. these 5000 steps are kept for analysis). The convergence of each chain was assessed using the Gelman &
Rubin (1992) test. We then extracted the best parameters by fitting a Gaussian function to each sampled posterior distribution, taking the mean and the standard deviation as the best estimate and the 1σ uncertainties, respectively, of each parameter. When the posterior distribution of a parameter was flat below or above a given value, we assumed that it is an upper or a lower limit, respectively (again, we ensured that values of these limits are not affected by our prior distributions). This optimisation method was applied to infer the Eddington ratio distribution for both our mass-independent and mass-dependent models.

Figure 1: Mass functions employed in our models. Top panel: mass functions for star-forming galaxies as reported in Davidzon
et al. (2017). Each colour corresponds to a different redshift bin, out to z=2.5 (see key). Bottom panel: same as top panel but for the quiescent galaxies.

2.1 A mass-independent λEdd distribution

As highlighted in the introduction, a number of studies have recently reported a difference in the λEdd distribution of quiescent and star-forming host galaxies (e.g. Georgakakis
et al., 2014; Wang
et al., 2017). We therefore take this “two-component” distribution as our starting point (i.e. we do not attempt to model a single universal λEdd distribution for all galaxies). Most studies use a broken power-law (or related functions) to represent the Eddington ratio distribution of AGNs (e.g. Aird
et al., 2013; Veale
et al., 2014; Jones et al., 2017; Weigel
et al., 2017). As such, we take a conservative approach by assuming a broken power-law for both the star-forming and the quiescent component of our λEdd distribution. Each broken power law is defined as,

p(λEdd)=A(λEddλbreak)γ1+(λEddλbreak)γ2,

(4)

where p(λEdd) is the probability of Eddington ratio, A is the normalisation, λbreak is the position of the break, and γ1 and γ2 are the slopes at low λEdd (i.e. below the break) and high λEdd (i.e. above the break), respectively. As such, each of our broken power law components are represented by four parameters, giving a total of eight free parameters to optimise for this mass-independent model (see § 2 for the optimisation method).

Following empirical results reported in e.g. Georgakakis
et al. (2014); Aird
et al. (2017a); Wang
et al. (2017), prior to optimisation, we assume that both the normalisation, ASF, and the position of the break, λSFbreak, of the star-forming component are each always higher than their quiescent analogues, AQui and λQuibreak, respectively. This also has the benefit of reducing the degeneracy of the model. However, it is important to stress that we do not assume any specific parameter values for the ratio between the two normalisations and the two positions of the breaks, we simply incorporate this information in the prior distributions by excluding some parts of the parameter space.

Incorporating the aforementioned assumptions, we optimised the eight free parameters of our mass-independent model. To investigate any redshift evolution of the Eddington ratio distribution, we performed the optimisation at z=0.3, z=0.5, z=0.7, z=1.0, z=1.2, z=1.5, z=1.7, z=2.0, z=2.3, z=2.5, z=2.5, z=2.7, z=3.0, and z=3.5, interpolating the mass functions of Davidzon
et al. (2017) split between star-forming and quiescent galaxies, and the total X-ray luminosity functions of Aird et al. (2015) at these redshifts. We stress that the Bayesian methodology adopted in Aird et al. (2015) to derive the observed X-ray luminosity functions (see Aird et al. (2015) for details on their method) means that they consider all X-ray detected sources when deriving their X-ray luminosity functions (which our analysis aims to fit). Therefore, there is no implied mass cut beyond the fact that their sources are detected in X-rays.

2.2 A mass-dependent λEdd distribution

Motivated by recent studies reporting that the λEdd distribution for star-forming galaxies is mass-dependent (e.g. Bongiorno
et al., 2016; Aird
et al., 2017a; Georgakakis
et al., 2017), we develop a second model that can accommodate three different λEdd distributions in three different mass bins for star-forming galaxies. As with our mass-independent model, we still split in terms of quiescent and star-forming galaxies, but for the star-forming component, we assume instead three different broken power-laws (see Eq. 4), one for each of our mass bins, which we define as: low mass (8<log(M∗/M⊙)<10), medium mass (10<log(M∗/M⊙)<11), and high mass (11<log(M∗/M⊙)<12). The choice of the boundaries for our mass bins is arbitrary and independent of any empirical studies, yet aims to cover a large range of stellar masses (i.e. 108<M∗/M⊙<1012).

Since assuming three broken power-laws generates 12 free parameters, we again rely on some assumptions based on the results of Aird
et al. (2017a) to help prevent degeneracies. First, we require that each broken power-law breaks at different λEdd values, in such a way that a higher mass bin breaks at a lower value of λEdd than its lower mass neighbour. Secondly, we assume that the three Eddington ratio distributions share the same slope at high Eddington ratios. Finally, we assume the normalisation of each of the three broken power-laws is such that the λEdd distribution above the break is always coincident (see Fig. 8). In making these assumptions, we reduce the parameter space to eight free parameters:

one normalisation, ASF, (from which the others are derived via our third assumption);

three break positions, λSFbreak, (one for each mass bin and ordered according to our first assumption);

three power-law slopes below the break, γSF1, (which are unconstrained);

a single shared power-law slope above the break, γSF2.

We stress that all these assumptions are inspired by empirical findings reported by Aird
et al. (2017a). We will demonstrate in § 3.1.1 that our previous model (i.e. the mass-independent model, see § 2.1) is able to reproduce the observed X-ray luminosity functions out to z∼2, and recreates both their star-forming and quiescent galaxy components in good agreement with observations, at least out to z∼1. As such, for this mass-dependent model, we assume that the Eddington ratio distribution for quiescent galaxies is unchanged from our previous model, and only update the star-forming component of the model Eddington ratio distribution by implementing the aforementioned mass-dependency. To do this, we optimise the eight free parameters that describe our mass-dependent Eddington ratio distribution for star-forming galaxies. Again, to investigate any redshift evolution of the mass-dependent Eddington ratio distribution, we performed the optimisation at z=0.1, z=0.3, z=0.5, z=0.7, z=1.0, z=1.3, z=1.5, and z=1.7, interpolating the mass functions for star-forming galaxies of Davidzon
et al. (2017) and the star-forming component of the X-ray luminosity functions as extracted from our mass-independent model. This second model relies on the ability of our first, mass-independent, model to split the X-ray luminosity functions in terms of star-forming and quiescent galaxies. As we show in § 3.1.1, our mass-independent model is unable to reproduce the total X-ray luminosity function beyond z∼2, and thus cannot be used to reliably split the luminosity function between star-forming and quiescent galaxies. As such, we cannot extend our mass-dependent model beyond z∼2.

2.3 Generating galaxies with AGNs

As we aim to investigate the connection between AGN accretion and host star-formation activity, we must model the star-forming properties of the host galaxies as well as the Eddington ratios of the AGNs. Since the distributions of specific SFRs (i.e SFR relative to the host stellar mass; sSFR) are now well-defined, we can simply adopt published models to generate our galaxy populations. As such, for this part of the model we use a PSM to attribute SFRs to our host galaxies. We use the PSM outlined in Bernhard
et al. (2014), which generates a population of mock star-forming galaxies using the mass functions of Ilbert
et al. (2013) and allocates SFRs using the relationship between stellar masses and SFRs, or Main Sequence (MS; e.g. Daddi
et al., 2007; Salim
et al., 2007; Elbaz
et al., 2011; Rodighiero
et al., 2011; Sargent et al., 2012; Schreiber
et al., 2015, for the MS) reported in Rodighiero
et al. (2011). The SFRs are then randomly scattered around the sSFR distribution of Sargent et al. (2012). This model successfully reproduces the observed infra-red and ultraviolet luminosity functions up to z∼6 (see Bernhard
et al., 2014, and references therein). However, for this work, we slightly update the empirical relationships used in Bernhard
et al. (2014). Briefly, we change the mass functions to the more recent ones reported in Davidzon
et al. (2017), we update the MS sSFR distribution to that reported in Schreiber
et al. (2015), and lastly, add a population of quiescent galaxies (using the quiescent galaxy mass functions of Davidzon
et al. 2017) with sSFRs at least a factor of ten below that of the MS (Ilbert
et al., 2013). We confirm that this updated model still reproduces the infra-red and ultraviolet luminosity functions. X-ray luminosities are then allocated using our Eddington ratio distributions derived in § 2.

Using this PSM, we generate a population of 33,925,192 galaxies (i.e. 32,939,834 – or 97 per cent – star-forming and 985,358 – or 3 per cent – quiescent galaxies), with stellar masses in the range 108<M∗/M⊙<1014. This corresponds to a 50 square degrees blank-field survey out to z=3. We use a very high (arguably un-physical) upper limit on our range of stellar masses to rule-out the possibility of being affected by any boundary effects. The large relative number of star-forming compared to quiescent galaxies in our model is related to the differences of their respective mass functions, particularly at masses M∗≲1010M⊙ and z≳1. As shown in Fig. 1, at these lower masses and higher redshifts, the mass function of star-forming galaxies steeply rises while that of quiescent galaxies sharply decreases. This leads to a large difference in the absolute number of star-forming galaxies compared to that of quiescent galaxies at lower masses (we have a cut at M∗=108M⊙ in our model). Above a stellar mass of M∗>1010M⊙, we have a total of 3,309,201 galaxies, among which 2,646,102 (80 per cent) are star-forming and 663,099 (20 per cent) are quiescent.

In this section we first present the results of our mass-independent model, followed by the results of our mass-dependent model. We demonstrate that while our mass-independent model can successfully reproduce the observed X-ray luminosity functions out to z≈2, it cannot at the same time reproduce the observed flat relationship between SFR and AGN luminosity (e.g. Rosario
et al., 2012; Stanley
et al., 2015). Instead, we show that this can be resolved by introducing a mass dependency in the λEdd distribution of AGNs hosted in star-forming galaxies.

3.1 The mass-independent model

X-ray luminosity functions

We show in Fig. 2 the fit to the X-ray luminosity functions for our mass-independent model. Up to z=1.75 our model X-ray luminosity functions are in very good agreement with those of Aird et al. (2015). However, beyond that redshift, our mass-independent model is unable to reproduce the observed X-ray luminosity functions. Fig. 2 also illustrates the model X-ray luminosity functions split between star-forming and quiescent galaxies. Overall, we find that the star-forming galaxies always dominate the X-ray luminosity functions at LX≳1042.5 erg s−1, which corresponds to the bulk of the AGN population. The contribution from quiescent galaxies in our mass-independent model, however, is only significant at lower X-ray luminosities (i.e. LX≲1042 erg s−1).

Beyond z∼2, our model cannot reproduce the observed X-ray luminosity functions. At the position of the break, our model X-ray luminosity functions at these redshifts are at least a factor of three below the observed ones. We also find that the contribution from quiescent galaxies is very low (by at least a factor of ten) compared to that of the star-forming galaxies. To explain the reasons for this, we examine the stellar mass functions of Davidzon
et al. (2017; see Fig. 1). They report that the contribution from quiescent galaxies to the mass function decreases as the redshift increases (the knee of the mass function for quiescent galaxies is at least a factor of ten in normalisation below that of the star-forming galaxies at z≳2; see Fig. 1). In our model, this lower contribution from the quiescent galaxies can be compensated by an increase in the normalisation of the model Eddington ratio distribution for quiescent galaxies (as our model relates the mass function to the Eddington ratio distribution). This would increase the contribution from quiescent galaxies to the total model X-ray luminosity functions in order to match the observed ones. However, in our model, the normalisation of the quiescent galaxy component cannot exceed that of the star-forming component (see assumptions in § 2.1). While looking at the normalisations of both the star-forming (ASF) and the quiescent (AQui) component in our model, we find that beyond z∼2, they are similar (i.e. ASF∼AQui). As a consequence, the contribution from quiescent galaxies to the total X-ray luminosity functions cannot increase enough to match the observed X-ray luminosity functions. If we relax this assumption, our model finds that the Eddington ratio distribution beyond z∼2 is fully dominated by quiescent galaxies, with the normalisation of the Eddington ratio distribution for star-forming galaxies to that of quiescent galaxies as low as ∼0.01 at z=2.25 and ∼0.001 at z=2.5. This is at odds with empirical measurements of the Eddington ratio distributions of e.g. Georgakakis
et al. (2014) and Wang
et al. (2017), that suggest a more equal splitting between quiescent and star-forming galaxies, as well as infrared studies that find large numbers of high-redshift AGNs residing in star-forming galaxies (e.g. Mullaney
et al., 2015; Stanley
et al., 2015; Netzer et al., 2016).

Finally, as a further check, we include in Fig. 2 the measured X-ray luminosity functions separated into star-forming and quiescent galaxies for z<1 as reported in Georgakakis
et al. (2014). Despite not including this information during our optimisation, we find good agreement between our model and these observed X-ray luminosity functions of star-forming and quiescent galaxies, increasing our confidence in the model at these lower redshifts.

Figure 2: Fit to the X-ray 2-10 keV luminosity functions of Aird et al. (2015) when optimising our mass-independent λEdd distribution model. Each panel represents a different redshift bin out to z=3, and as indicated at their bottom right-hand side. The downward and upward empty black triangles are the X-ray luminosity functions of Aird et al. (2015) for the soft and the hard band, respectively. The thick orange line in each panel is for the model total X-ray luminosity function given by our mass-independent model, the blue dashed line is for that of the star-forming component and the red triple dashed line shows the quiescent component. The blue empty stars and the red empty circles at z<0.90 are for the star-forming and the quiescent component of the X-ray luminosity functions, respectively, as reported in Georgakakis
et al. (2014) out to z∼1. We hatched the two panels at z=2.25 and z=3.0 to indicate that our mass-independent λEdd distribution model is unable to reproduce the observed X-ray luminosity functions at these redshifts (see § 3.1.1).

The Eddington ratio distribution

We show in § 3.1.1 that our mass-independent model is able to reproduce the observed X-ray luminosity function out to z∼2. We now present in Fig. 3 the Eddington ratio distributions at various redshifts after optimisation of our mass-independent model (normalised to coincide at λEdd= 10−3). The shape of our total Eddington ratio distribution (see left-hand panel in Fig. 3) is consistent with a broken power-law, as requested by our model and in agreement with previous studies (e.g. Aird
et al., 2013; Jones et al., 2017). However, we find a significant difference in the slope below the break (γ1 in Eq. 4) of the Eddington ratio distribution for star-forming hosts compared to that of quiescent hosts (see central and right-hand panels in Fig. 3 for the Eddington ratio distribution of star-forming and quiescent host galaxies, respectively). We report that below z∼2, the slope at low-λEdd for the star-forming component is rising (although we are only able to place upper limits). As a consequence, our model Eddington ratio distribution of AGNs in star-forming galaxies is consistent with a peaky3 distribution, similar to the light-bulb shaped distribution sometimes explored in past studies (e.g. Veale
et al., 2014; Stanley
et al., 2015). This suggests that AGNs hosted in star-forming galaxies have typical λEdd higher than AGNs in quiescent hosts, yet the λEdd of AGNs in quiescent galaxies span a broader range of values. Finally, we notice an overall shift of the knee of the Eddington ratio distribution to higher λEdd values as redshift increases. This implies a typical AGN activity that increases as redshift increases. At z≳2, we find that the shape of the Eddington ratio distribution of AGNs in star-forming galaxies becomes consistent with a broken power-law, in contrast to the peaky Eddington ratio distribution reported for lower redshifts. However, as demonstrated in § 3.1.1 our mass-independent model is unable to reproduce the X-ray luminosity functions at these higher redshifts. We report in Table 1 the redshift evolution of the parameters that define these Eddington ratio distributions.

Figure 3: The total mass-independent Eddington ratio distribution (left-hand panel), its star-forming component (central panel), and its quiescent component (right-hand panel) evaluated at z=0.2, z=0.5, z=1.0, z=1.5, z=2.0, and z=2.5 (colour coded; see key). The downward arrows indicate the overall effects of the upper and lower limits found in the parameters that describes these distributions. For the star-forming and quiescent panels, we also indicate with faint grey lines the corresponding total Eddington ratio distribution. The normalisation is such that the total Eddington ratio distribution at z=0.2 integrates to unity, applying an arbitrary cut at λEdd=10−7. For ease of comparison, we then re-normalise all the distributions at higher redshifts such that their values at log(λEdd)=-3.0 coincide with that at z=0.2, hence the arbitrary unit of the ordinate axis.

Parameters

Intercepts

Slopes

log(ASF)

−2.07±0.19

0.13±0.07

for z<1.7

−1.15±0.03

for z>1.7

log(AQui)

−2.80±0.23

0.37±0.10

for z<1.7

−1.20±0.11

for z>1.7

log(λSFbreak)

−2.36±0.06

0.65±0.02

log(λQuibreak)

−2.62±0.16

0.66±0.06

γSF1

≲−3.0

0.0

for z<1.7

−0.96±0.51

0.25±0.12

for z>1.7

γQui1

0.28±0.05

0.07±0.02

γSF2

2.29±0.04

0.01±0.02

for z>1.7

0.15±0.03

for z>1.7

γQui2

≳2.7

0.0

for z<1.7

3.75±2.33

−0.53±0.63

for z>1.7

Notes: Slopes and intercepts are given for an evolution as a function of (1+z). The intercept for z>1.7, when assuming a break in the z evolution of the parameter, is given by the continuity at z=1.7 (i.e. [intercept for z>1.7]=(1+1.7)×([slope for z<1.7]−[%slopefor$z$>1.7])+[intercept for z<1.7]).

Table 1: Redshift evolution of the parameters that describe the Eddington ratio distribution for star-forming and quiescent galaxies, given by Eq. 4. SF and Qui labels stand for star-forming and quiescent, respectively. The slopes and intercepts are given for an evolution in (1+z). Uncertainties on the parameters are the standard deviation of 1000 Monte Carlo realisations (i.e. 1σ).

SFR in bins of X-ray luminosities

As previously highlighted, a key test of any model that describes the Eddington ratio distribution is whether it can reproduce other observed features of the AGN population besides just the (X-ray) luminosity functions. Here, we test our PSM (which incorporates pre-defined sSFR distributions and our own optimised mass-independent Eddington ratio distribution; see § 2.3) by assessing whether it can reproduce the observed flat relationship between SFR and X-ray luminosity as reported in many observational studies (e.g. Rosario
et al., 2012, 2013; Stanley
et al., 2015; Bernhard et al., 2016). We do this by calculating the mean-average SFRs of the AGNs in our PSM in bins of 0.5 dex in X-ray luminosity, taking into account both star-forming and quiescent galaxies (see § 2.3). As shown in Fig. 4, our mass-independent model predicts a strong correlation between these SFRs and X-ray luminosities (with an average slope of ∼0.6), at least up to z∼2. This is in good agreement with empirical results for our lowest redshift bin (i.e. z=0.35), but strongly contradicts the flat observed relationship between SFR and X-ray luminosity at higher redshifts (e.g. Rosario
et al., 2012; Stanley
et al., 2015). It is important to note that other studies also find a similar strongly increasing relationship between SFR and X-ray luminosity when employing a peaky (i.e. similar in shape than our peaky distribution; see central panel in Fig. 3) Eddington ratio distribution for star-forming galaxies (e.g. Veale
et al., 2014; Stanley
et al., 2015). We will show in § 4 that this strong correlation is due to the AGN host mass distribution and the SFR/stellar-mass relationship (i.e. MS) induced by the peaky Eddington ratio distribution found for our mass-independent model.

As SFR is correlated to stellar mass via the MS, it is possible that the empirical flat relationship between SFR and X-ray luminosity is a consequence of an observational bias toward higher mass galaxies. To test this, we recalculate from our model the average SFR in bins of X-ray luminosity, but now excluding galaxies with M∗<109.5M⊙. This mass limit roughly corresponds to that reached in the highest redshift bins of Stanley
et al. (2015). We show in Fig. 4 that although this mass cut causes the model SFR/X-ray relationship to rise slightly at z=2.0 and LX≲1043.5 erg s−1, it is not sufficient to reproduce the observed flat relationship between SFR and X-ray luminosity.

Figure 4: The mean-average SFR in bins of X-ray 2-10 keV luminosity at z=0.35, z=0.65, z=1.15, and z=2. The black thick lines shows the correlation predicted by our PSM that incorporates AGNs using our mass-independent Eddington ratio distribution split between star-forming and quiescent galaxies. The error bars show the 1σ uncertainties on the mean for average SFR. The thick black dashed line is identical than the thick black line but for galaxies with M∗≳109.5M⊙. Various symbols are for a compilation of published SFR/X-ray relationships (see keys; i.e. Rosario
et al., 2012, 2013; Stanley
et al., 2015; Bernhard et al., 2016). The orange dashed line at z=1.15 shows the relationship reported in Silverman
et al. (2009). Finally, the black dashed line represents the correlation displayed by AGN-dominated systems in Netzer (2009). Our mass-independent model predicts a significant correlation between SFR and X-ray luminosity which is at odds with empirical results for z≳0.5.

3.2 The mass-dependent model

Since our mass-independent model fails to reproduce the X-ray luminosity functions at z≳2 and is unable to reproduce the flat relationship between SFR and X-ray luminosity (at z>0.5), we consider whether relaxing the mass-independence requirement (as suggested by recent studies, e.g. Bongiorno
et al. 2016; Aird
et al. 2017a; Georgakakis
et al. 2017) can help to resolve these issues. As explained in § 2.2, for this model we assume that the quiescent component remains unchanged from the mass-independent model meaning that we only fit the star-forming X-ray luminosity functions at z<2. These are extracted from the results of our mass-independent model, and is justified by the ability of our mass-independent model to reproduce the X-ray luminosity functions of star-forming and quiescent galaxies to z∼1 from Georgakakis
et al. (2014).

X-ray luminosity functions

We show in Fig. 5 the fit to the star-forming component of the X-ray luminosity functions using our mass-dependent Eddington ratio distribution model to z=1.75. Again, we do not extend our mass-dependent model beyond z∼2 since it relies on the ability of our previous – mass-independent – model to split the X-ray luminosity functions in terms of star-forming and quiescent components (see § 2.2). Out to z∼2, our mass-dependent model does a good job at fitting the star-forming X-ray luminosity functions in all our redshift bins. Since our mass-dependent model assumes an λEdd distribution split into three different mass bins (see § 2.2), we are also able to derive the X-ray luminosity function of each mass bin (see Fig. 5). In doing so, our model predicts that the highest mass galaxies (i.e. M∗>1011M⊙) dominate the star-forming X-ray luminosity functions across almost the full range of X-ray luminosities (i.e. 1041<LX<1046 erg s−1). The exception being in our lowest redshift bin (i.e. z=0.11) where 1010<M∗/M⊙<1011 galaxies contribute marginally more than our highest mass galaxies at LX≳1043.5 erg s−1. In general, however, these medium mass galaxies only contribute significantly (i.e. >10 per cent) at luminosities above the knee of the X-ray luminosity functions (i.e. LX>1044 erg s−1). By contrast, the lowest mass galaxies (i.e. M∗<1010M⊙) provide almost no contribution to the X-ray luminosity functions of star-forming galaxies at LX≳1042.5erg s−1 (i.e. <10 per cent), although we are only able to place upper limits on this contribution from lower mass galaxies.

Our mass-dependent model suggests that AGNs hosted by star-forming galaxies with M∗>1011M⊙ dominate the X-ray luminosity functions. This is consistent with Georgakakis
et al. (2017), which demonstrates that AGNs with LX>1041\-−42 erg s−1 significantly contribute to the AGN host stellar mass function for masses M∗>1011M⊙. This is also consistent with results from Bongiorno
et al. (2016) which show that the host stellar mass function of AGNs with LX>1043 erg s−1 peaks at M∗∼1011M⊙, out to z=2. However, the vast majority of AGN host galaxies with LX<1043 erg s−1 display typical stellar masses in the range M∗∼1010\-−11M⊙(e.g. Bongiorno
et al., 2016; Aird
et al., 2017b), which contrasts with our findings of a model X-ray luminosity function dominated by host galaxies with M∗>1011M⊙. However, our model is limited by the choice for the boundaries of our mass bins. Therefore, should we change our highest mass bin to incorporate galaxies with M∗∼1010.5M⊙, we would expect to find that galaxies with M∗>1010.5M⊙ dominate the X-ray luminosity functions, in agreement with recent observational studies (e.g. Bongiorno
et al., 2016; Aird
et al., 2017b; Wang
et al., 2017).

Figure 5: Fit to the star-forming component of the X-ray 2-10 keV luminosity functions using our mass-dependent Eddington ratio distribution model up to z=1.75. Each panel shows a different redshift bin as indicated on their bottom right-hand side. The black line in each panel shows the star-forming component of the X-ray luminosity function (derived from our mass-independent model), and the orange line shows the best fit given by our mass-dependent model. The dashed, single-dot dashed, and triple-dot dashed lines show the X-ray luminosity functions for our highest (11<log(M∗/M⊙)<12), medium (10<log(M∗/M⊙)<11), and lowest (8<log(M∗/M⊙)<10) mass bin, respectively. The black arrows indicate the overall effects of upper and lower limits on the parameters that define our mass-dependent Eddington ratio distribution. The blue stars and the red circles show the X-ray luminosity functions for the star-forming and quiescent galaxies, respectively, as reported in Georgakakis
et al. (2014). The light grey downward and upward triangles show the total X-ray luminosity functions of Aird et al. (2015) for the soft and the hard band, respectively.

SFR in bins of X-ray luminosities

Having confirmed that our mass-dependent model is able to reproduce the X-ray luminosity functions for star-forming galaxies up to z∼2, we now consider whether the corresponding PSM reproduces the observed flat relationship between SFR and X-ray luminosity (e.g. Rosario
et al., 2012; Stanley
et al., 2015). We show in Fig. 6 the mean-average SFR of AGNs split into 0.5 dex-wide bins of X-ray luminosity at similar redshifts as those of Stanley
et al. (2015). Again, we include quiescent galaxies in our model when calculating these averages. Contrary to our mass-independent model, we find that our mass-dependent model predicts a flat relationship between SFR and X-ray luminosity in very good agreement with observations (see Fig. 6). It is important to stress that the model was not optimised to recreate the flat SFR/X-ray luminosity relationship. Therefore, we conclude that the mass-dependent model is able to reproduce the X-ray luminosity functions for star-forming galaxies (with a good agreement with observations at least up to z∼1) while also independently reproducing the observed flat relationship between SFR and X-ray luminosity out to z∼2.

Figure 6: The mean-average SFR in bins of X-ray 2-10 keV luminosities predicted by our mass-dependent model (black thick line) at z=0.35, z=0.65, z=1.15, and z=2.0. The error bars represent the 1σ uncertainties on the mean for average SFR. The various symbols are from a compilation of published relationships (see keys; i.e. Rosario
et al., 2012, 2013; Stanley
et al., 2015; Bernhard et al., 2016). The orange dashed line shows the relationship between average SFR and X-ray luminosity found in Silverman
et al. (2009), and the black dashed line represents the relationship for AGN-dominated systems, as reported in Netzer (2009). Although it is not implemented in the model, our mass-independent model successfully reproduces the observed flat SFR/X-ray relationship out to z=2.

The evolution of the λEdd distribution

Having demonstrated that our mass-dependent model is able to reproduce both the X-ray luminosity functions for star-forming galaxies and the flat SFR/X-ray luminosity relationship to z∼2 we now explore how the eight parameters (see § 2.2) that define our mass-dependent model change with redshift. We show in Fig. 7 the redshift evolution of these parameters. We report that the overall normalisation of the star-forming component (ASF) increases very slightly with increasing redshift. The position of the break (λSFbreak) in each mass bin shifts toward higher λEdd values. This is consistent with the overall increase of the position of the break observed in our mass-independent model (i.e. all mass bins collapsed together; see § 3.1.2), suggesting a higher typical accretion rate for AGNs at higher redshifts. Regarding the slopes at lower λEdd, (γSF1), we find that it does not evolve with redshift for any of our mass bins. However, although γSF1 for our highest mass bin is well constrained, we have large uncertainties associated with γSF1 for our medium mass bin, and upper limits for our lowest mass bin (see Fig. 7). These weak constraints on γSF1 for our medium and lowest mass bins are a direct consequence of their low contributions to the total Eddington ratio distribution for star-forming galaxies (see Fig. 8), as was hinted by their low contribution to the total X-ray luminosity functions (see § 3.2). We also find that the shared slope at high λEdd (γSF2) does not change significantly with redshift. Overall, our findings suggest that there is a suppression of low-λEdd (i.e. λEdd≲0.1) for AGNs hosted in lower mass galaxies (i.e. M∗≲1010\-−11M⊙). Finally, we performed a linear fit of the redshift evolution of each parameter, using the standard deviation of 1000 Monte Carlo realisations for the uncertainties on the fitting parameters. We report the results of this fit in Table 2.

Each of the aforementioned trends can be seen in the evolution of the overall Eddington ratio distribution, which we plot out to z=2 in Fig. 8. In this figure, we also plot the Eddington ratio distributions of quiescent host galaxies found in our mass-independent model at similar redshifts. In Fig. 8, for the left-hand panel, each distribution is normalise such that the total Eddington ratio distribution (i.e. the combination of the star-forming and quiescent component) integrates to unity after applying an arbitrary cut at λEdd= 10−7. For the right-hand panels, each total distribution is also normalised such that it integrates to unity after applying an arbitrary cut at λEdd= 10−7.

Parameters

Intercepts

Slopes

log(ASF)

-1.57±0.08

0.03±0.04

log(λlowmassbreak)

>-2.53

1.27

log(λmediummassbreak)

-2.05±0.15

0.68±0.08

log(λhighmassbreak)

-2.60±0.08

0.70±0.04

γlowmass1

<0.5

0.0

γmediummass1

-1.0±2.67

-0.12±1.30

γhighmass1

0.45±0.05

-0.04±0.02

γ2

2.21±0.06

-0.03±0.03

Notes: Slopes and intercepts are given for an evolution as a function of (1+z). The intercept for z>1.7 when assuming a break in the z evolution of the parameters is given by the continuity at z=1.7 (i.e. [intercept for z>1.7]=(1+1.7)×([slope for z<1.7]−[%slopefor$z$>1.7])+[intercept for z<1.7]).

Table 2: Redshift evolution of the parameters that describe the Eddington ratio distribution of star-forming galaxies for our mass-dependent model. The slopes and intercepts are given for an evolution as a function of (1+z). Uncertainties on the fitting parameters are extracted via 1000 Monte Carlo realisations (i.e. 1σ).Figure 7: Redshift evolution of the parameters that define our mass-dependent Eddington ratio distribution for star-forming galaxies at z=0.1, z=0.3, z=0.5, z=0.7, z=1.0, z=1.3, z=1.5, and z=1.7. Error bars show the 3σ uncertainties. The upward and downward arrows indicate the lower and upper limits, respectively. The two top panels present the shared parameters among our three mass bins, i.e. the normalisation, logASF, and the slope at high Eddington ratio, γSF2. The left-hand panels are the three different break positions λSFbreak, one for each mass bin (as indicated in the bottom left-hand side in each panel), and the right-hand panels are the same but for that of the lower Eddington ratio slope, γSF1. The black lines show the best fit with redshift of each parameter.Figure 8: The Eddington ratio probability distributions of star-forming and quiescent galaxies in our mass-dependent model at z=0.2, z=0.5, z=1.0, z=1.5, and z=2.0. The left-hand panel shows both the star-forming (plain lines) and the quiescent components (single-dashed lines) of the Eddington ratio probability distribution. The normalisation is such that the total Eddington ratio distribution at each redshift integrates to unity after applying an arbitrary cut at λEdd= 10−7. The right-hand panels are, for each redshift (indicated on their top right-hand corner), the contribution of the different mass bins to the Eddington ratio distributions of AGNs hosted in star-forming galaxies. In each of these panels, the dashed line is for the highest mass bin (11<log(M∗/M⊙)<12), the single-dotted dashed line is for the medium mass bin (10<log(M∗/M⊙)<11), and the triple-dotted dashed line is for the lowest mass bin (8<log(M∗/M⊙)<10). The downward arrows for the slope at low-λEdd in our lower mass bin indicate that we have upper limits. The grey area illustrates the large uncertainties for the slope at low-λEdd in our medium mass bin (see Fig. 7). We normalised each of the total Eddington ratio distribution such that it integrates to unity, applying an arbitrary cut at λEdd=10−7.

4.1 Further checks to our mass-dependent model

As a further check of our mass-dependent model, we compared the model Eddington ratio distribution at z=1 to the empirical one reported in Wang
et al. (2017). As our various assumptions for this model are based on observations from Aird
et al. (2017a), comparing our mass-dependent Eddington ratio distribution to that empirical of Wang
et al. (2017) – instead of Aird
et al. (2017a) – constitutes a more independent test for our model. We show this comparison in Fig. 9. We find very good agreement between our model Eddington ratio distributions to that of Wang
et al. (2017) at z=1 for both star-forming and quiescent component (at least for log(λEdd)<-0.5). This strongly supports our mass-dependent model for the Eddington ratio distribution of star-forming galaxies (since the mass-independent model was predicting a peaky Eddington ratio distribution for star-forming galaxies which would contrasts with empirical results from Wang
et al. 2017). We predict that the low-λEdd end of the Eddington ratio distribution for star-forming galaxies is dominated by galaxies with M∗≳1011M⊙. However, Wang
et al. (2017) reported that their sample of AGNs have typical stellar masses between 1010.5<M∗/M⊙<1011 (once corrected for the differences in the IMFs between Salpeter 1955 for this work and Chabrier 2003 for Wang
et al. 2017), which is half a dex below our highest mass bin. As mentioned in § 3.2.1, this discrepancy in the host stellar masses between our model and observations is a consequence of the chosen boundaries for each of our mass bins.

Figure 9: Comparison of our mass-dependent model Eddington ratio distribution to empirical results of Wang
et al. (2017) at z=1. The blue thick line is the total model Eddington ratio distribution for star-forming galaxies, while the dashed line, the dotted-dashed line, and the triple-dotted dashed line indicate our highest, medium, and lowest mass bin contributions, respectively. The blue circles are empirical data of the star-forming component from Wang
et al. (2017) at z∼1, while the red circles are that of the quiescent component.

As a final check of our mass-dependent model, we show in Fig. 10 how the normalised average SFR (i.e. SFR relative to that of the MS) changes with λEdd in our mass-dependent model. Although it is not incorporated in the optimisation, we predict a slight enhancement of average normalise SFR at higher λEdd (i.e. λEdd≳1) compared to that of their lower λEdd counterpart (at least at z≳1.2). This slight enhancement of normalised average SFR at higher λEdd is also observed in Bernhard et al. (2016). In Fig. 10 we only consider star-forming galaxies for our model since we do not have a good prescription for the normalised SFRs of quiescent galaxies. This is a possible reason for the discrepancy at lower λEdd between our model and the empirical results (i.e. that does contain quiescent galaxies). We also stress that in Fig. 10 our highest redshift bin, i.e. 1.8<z<2.9 constitutes an extrapolation of our model. However, the results are still consistent with observations (i.e. a slight enhancement of the normalised average SFR at higher λEdd).

Figure 10: The normalised average SFR versus the Eddington ratio up to z∼3. The blue lines show the prediction for our mass dependent model. the orange circles show the results report in Bernhard et al. (2016). The grey area represent the scatter around the MS as reported in Schreiber
et al. (2015), the dashed black line being the MS.

One important result of our mass-independent model is that the Eddington ratio distribution for star-forming galaxies is better represented by a peaky distribution (i.e. a narrow distribution), rather than the more common broken power-law (i.e. a broader distribution). To explain this, we refer to Caplar
et al. (2015; see also Weigel
et al. 2017), where they demonstrate that the steepness of the faint-end slope of the X-ray luminosity function is determined by either the low-mass end of the mass function or the low-λEdd slope of the Eddington ratio distribution (since combined together to form the X-ray luminosity function) whichever is the steeper. Since that the faint-end of the observed X-ray luminosity function of Aird et al. (2015) flattens at higher redshifts whilst, in contrast, the low-mass end of the mass function for star-forming galaxies steepens with redshift (see Fig. 1), our mass-independent model attempts to reduce the contribution from star-forming galaxies at lower λEdd, leading the peaky Eddington ratio distribution for star-forming galaxies found by our mass-independent model. However, beyond z∼2, the low-mass end of the mass function for star-forming galaxies of Davidzon
et al. (2017) is already too steep compared to the flat faint-end of the X-ray luminosity function, such that our model is unable to reproduce the observed X-ray luminosity function. The steepening of the mass functions combined with the flattening of the X-ray luminosity functions suggests a mass-dependent Eddington ratio distribution for star-forming galaxies.

We further saw that our mass-independent model failed to reproduce the observed flat SFR/X-ray relationship (see Fig. 4). To explain why this happens, we show in Fig. 11 which galaxy masses populate different parts of the SFR/X-ray plane. This plot shows that the mass-independent model produces a strong correlation between X-ray luminosity and stellar mass. This is a direct result of the narrower peaky (in contrast to a broader broken power law) Eddington ratio distribution for star-forming galaxies that this model demands in order to reproduce the observed X-ray luminosity functions, at least out to z∼2. The narrowness of this Eddington ratio distribution means that a galaxy of a given mass can only produce an AGN within a limited range of luminosities. When we then include the correlation between SFR and stellar mass for MS galaxies (e.g. Schreiber
et al., 2015), the consequence is a correlation between SFR and X-ray luminosity. This correlation between SFR and X-ray luminosity was also noted by Veale
et al. (2014) when using their light-bulb Eddington ratio model, which is similar in shape to our peaky distribution for star-forming host galaxies.

Overall, there is a large discrepancy with our mass-independent model that demands a narrow distribution for the Eddington ratio distribution of star-forming galaxies in order to reproduce the X-ray luminosity functions, but also requires a broader Eddington ratio distribution to reproduce the flat SFR/X-ray relationship. This strongly suggests that the Eddington ratio distribution must be somehow mass-dependent, such that low Eddington ratios are suppressed in low mass galaxies to be able to reproduce the X-ray luminosity functions while still being broad to reproduce the flat relationship between SFR and X-ray luminosity. Our mass-dependent model provides further constrains on how the Eddington ratio distribution for star-forming galaxies should be mass-dependent.

Figure 11: Stellar mass distributions (colour coded) in the SFR/X-ray luminosity parameter space, using our PSM and assuming our mass-independent model for the Eddington ratio distribution. Each panel shows a different redshift, as indicated on their bottom-right hand side. Each discrete colour indicates a stellar mass density in bins of log(M∗)=0.5 (see colour-bar).

4.3 Caveats to our models

While deriving the λEdd distribution, our model is limited by the functional form chosen for each component of the λEdd distribution and the associated assumptions used to avoid degeneracies during fitting (see § 2). An obvious caveat to our model is that we cannot explore the full range of possible Eddington ratio distributions. However, we find that, should we relax the assumptions made prior to the fit of the X-ray luminosity functions, our solutions are too degenerate to form a set of meaningful parameters for the Eddington ratio distributions, meaning that the X-ray luminosity functions lack the information needed to constrain the Eddington ratio distribution and its mass-dependency. Using these assumptions is therefore a requirement in our modelling approach. As stressed in § 2, all our assumptions are motivated by the findings of recent studies (e.g. Jones et al., 2017; Aird
et al., 2017a), and the results predicted by the mass-dependent model are in good agreement with many recent observational studies, i.e. the observed X-ray luminosity functions out to z∼2 – including when split in terms of star-forming and quiescent galaxies (e.g. Georgakakis
et al., 2014; Aird et al., 2015), the flat relationship between SFR and X-ray luminosity (e.g. Rosario
et al., 2012; Stanley
et al., 2015), the empirical Eddington ratio distribution split between star-forming and quiescent component at z=1 (Wang
et al., 2017), and the enhancement of the star-forming properties at higher λEdd(Bernhard et al., 2016). As such while we do not claim that our Eddington ratio distributions are universal, they provide a simple means by which to explore the AGN-galaxy connection.

One of the main caveats for our mass-dependent model is that we are limited by the choice of the stellar mass bins. As a consequence, we find that our dominant component for the Eddington ratio distribution and the X-ray luminosity function is our highest mass bin with M∗>1011M⊙. This is qualitatively consistent with observations (e.g. Georgakakis
et al., 2014; Aird et al., 2015; Wang
et al., 2017), yet roughly half a dex above the observed typical host mass for AGNs. Investigating this in more detail would require additional mass bins. However, this would also generates a larger number of free parameters, and therefore degeneracies. We have chosen these mass bins to be broad enough to cover a wide range of masses, yet probing any mass dependency of the Eddington ratio distribution for AGNs in star-forming hosts.

4.4 Extending our mass-dependent model to higher redshifts

We recall that we fit the mass-dependent model to the X-ray luminosity functions of star-forming galaxies only, which was isolated using our mass-independent model. The results of our mass-dependent model therefore depends on the reliability of our mass-independent model to separate the X-ray luminosity functions into star-forming and quiescent components (see § 2). Thankfully, we can verify this up to z∼1 using the results reported in Georgakakis
et al. (2014). However, beyond z∼1 we do not have the empirical results to perform this check, while at z≳2 our mass-independent model is actually inconsistent with observations (i.e. too many AGNs in quiescent galaxies). Faced with this situation where we cannot reliably exploit the results of our mass-independent model, we now investigate whether we can fit the observed X-ray luminosity functions at z=2.25 using both a mass-dependent star-forming component and a mass-independent quiescent component λEdd distributions. This results in a total of 12 free parameters to optimise (i.e. eight parameters for the star-forming component and four parameters for the quiescent one). We use the same assumptions for the λEdd distribution of star-forming galaxies as described in § 2, and demand that the normalisation of the Eddington ratio distribution for quiescent galaxies lies below that of star-forming galaxies (i.e. consistent with the lower redshift bins).

We show in Fig. 12 the results of this new model. Contrary to our mass-independent model, we are now able to reproduce the total observed X-ray luminosity functions at z=2.25 by placing the majority of AGNs in star-forming galaxies. In particular, we find that the X-ray luminosity functions at z=2.25 is dominated by the highest mass star-forming galaxies, with a smaller contribution from medium mass star-forming galaxies. The contribution from our lowest mass bin is consistent with zero. Similarly, we can only place upper limits on the contribution from quiescent galaxies to the total X-ray luminosity functions at z=2.25. We conclude, therefore, that while we are able to reproduce the total X-ray luminosity functions using a mass-dependent λEdd distribution for star-forming galaxies (and ensuring that it is dominated by star-forming galaxies), this solution is degenerate with the level of contribution from quiescent galaxies. To break this degeneracy will require the separation of the high-redshift X-ray luminosity functions into quiescent and star-forming components (i.e. as performed by Georgakakis
et al. (2014) at z<1).

Figure 12: Fit of the X-ray luminosity functions at z=2.25 of Aird et al. (2015) for our model assuming a mass-dependent Eddington ratio distribution for star-forming galaxies and a mass-independent one for that of quiescent galaxies. The black downward and upward triangles show the X-ray luminosity functions of Aird et al. (2015) for the soft and the hard band, respectively. The orange line shows the total X-ray luminosity functions derived using our model, the blue lines are for that of the star-forming component in each different mass bin (see keys), and the quiescent component is shown with a red line. Downward black arrows indicate upper limits.

Motivated by recent results reporting a different Eddington ratio distribution for star-forming and quiescent galaxies (Wang
et al., 2017; Aird
et al., 2017a), we attempt to constrain these distributions by using an analytical model to fit the observed X-ray luminosity functions of Aird et al. (2015), assuming the galaxy mass functions of Davidzon
et al. (2017). Our first model assumes two mass-independent λEdd distributions: one for star-forming galaxies and another for quiescent galaxies. After optimisation, we find that this model is able to reproduce the X-ray luminosity functions out to z∼2 (see § 3.1.1), but demands a peaky distribution for the Eddington ratio distribution of star-forming galaxies. Despite this, our mass-independent model fails to reproduce the observed flat SFR/X-ray relationship when incorporated into a PSM for galaxies (see § 3.1.3). We argue that this failure arises because a mass-independent model places too many low-luminosity AGNs in low mass galaxies (see § 4). To resolve this problem, we develop a second model in which the Eddington ratio distribution for star-forming galaxies is allowed to be different in three mass bins (motivated by observation from e.g. Aird
et al. 2017a; Georgakakis
et al. 2017). By suppressing low Eddington ratio AGNs in low mass galaxies, this mass-dependent model is able to reproduce the X-ray luminosity functions for star-forming galaxies while simultaneously reproducing the observed flat relationship between SFR and X-ray luminosity (see § 3.2.2). Finally, we also find that the mass-dependent model is consistent with empirical Eddington ratio distribution for star-forming and quiescent galaxies at z=1, and is able to reproduce the slight enhancement of normalised average SFR at higher λEdd, as reported in Bernhard et al. (2016).

Overall, we conclude that, in our model, a suppression of lower AGN accretion activity in low mass galaxies is required to reproduce both the observed X-ray luminosity functions and the observed flat SFR/X-ray relationship.

We thank the referee for the useful comments that help to significantly improve the quality of the paper. EB thanks C. Tadhunter and D. Alexander for the helpful discussions. EB thanks the University of Sheffield for the support grant.