Constraints from dwarf galaxies on black hole seeding and growth models
with current and future surveys

By construction, SAMs are set up to explore a range of seeding and growth modes for a population of black holes as a function of redshift. Time slices from these SAMs can then be compared to measured black holes. In this work, we focus on the SAM developed by Ricarte & Natarajan (2018b) and explore new predictions for the low mass end of the local black hole mass function. In particular, we include the properties of the host galaxies explicitly in the scheme and do so flexibly by exploring multiple scaling relations between them and MBH properties. It is this new feature that allows us to delve more deeply into further predictions in combination with the Bayesian inference framework in this work.

Intriguingly, as first noted and discussed in detail in Ricarte & Natarajan (2018a), signatures of seeding survive to late times via the observed occupation fraction though there is dependence on the details of their mass assembly history. There are currently other SAMs from Somerville et al. (2008) & Bower et al. (2010) and updated versions thereof, that do not focus on seeding and include recipes for star and galaxy formation as well as MBH growth, and these models inevitably have many more tuneable parameters. Here, we stick with a pared-down, phenomenologically driven model that is primarily focused on seeding signatures from nearly local observations ($z\sim 0$). The Ricarte & Natarajan (2018a) SAM considers two seed models: either a massive and rare “heavy” seed typically produced from direct collapse following the seeding prescription of Lodato & Natarajan (2006); Volonteri et al. (2008), or abundant “light” seeds drawn from a power law initial mass function between 30 and 100 solar masses motivated by cosmological simulations Hirano et al. (2015). For each seed, the SAM explored in Ricarte & Natarajan (2018b) follows three different MBH growth prescriptions:

• •

  Power Law (PL) Model: In this model, major halo mergers trigger a burst of growth at the Eddington rate until MBHs reach the mass determined by the observed local $M_{\bullet}-\sigma_{*}$ relation. When MBHs are not in the burst mode, their accretion rates are drawn from a power law tuned to help reproduce low-redshift luminosity functions (Ricarte & Natarajan, 2018b).

• •

  AGN Main Sequence (AGNMS) Model: Major halo mergers trigger burst mode growth at the Eddington rate until MBHs reach the $M_{\bullet}-\sigma_{*}$ relation. When MBHs are growing secularly, their accretion rates are set to be $10^{-3}\ \times$ the star formation rate of their host. This mode is motivated by empirical relation determined by Mullaney et al. (2012). Since star formation is not actually tracked in this SAM, scaling relations are instead used to estimate stellar masses of galaxies as a function of halo mass and redshift (Ricarte & Natarajan, 2018b).

• •

  Broad-line Quasars (BLQ) Model: Major halo mergers (defined to be those with a halo mass ratio of at least 1:10) trigger burst mode growth. Additionally, at each time step the accretion rate is drawn from a distribution of Eddington rates observationally inferred for broad-line quasars in SDSS (Kelly & Shen, 2013). Specifically, we adopt a redshift-evolving log-normal distribution fit by Tucci & Volonteri (2017)¹¹1Here, the probability distribution of Eddington ratios $\lambda$ is given by $P(\lambda)=\exp(-(\ln\lambda-\ln\lambda_{c}(z))^{2}/(2\sigma(z)^{2})/(2\pi\sigma(z)\lambda)$, where $\log\lambda_{c}=\mathrm{max}(\{-1.9+0.45z,\log 0.03\})$ and $\sigma(z)=\mathrm{max}(\{1.03-0.15z,0.6\})$. No additional steady accretion mode for mass growth is added since the AGN luminosity function is adequately predicted with this prescription (Ricarte et al., 2019). However, as discussed in this work, this assumption fails when considering deep surveys of dwarf galaxies that are not selected a priori to be AGN.

In this model, only 10$\%$ of the halo mergers lead to black hole mergers. This is done to prevent the overgrowth of BHs in the most massive halos, whose BHs would otherwise grow mostly from BH-BH mergers at low-redshift. We do not expect this assumption to have an impact on the dwarf galaxies considered in this work, which have much lower merger rates. In our study, we do not include satellite halos, since the SAM does not include a prescription for environmental processes such as ram-pressure/tidal stripping, which in turn is expected to affect the stellar mass as well as the growth of MBH. To increase our sample size, we compile the output from three snapshots at: $z=0$, $z=0.1$, and $z=0.2$. This yields a total of 6,579 halos per seed $\times$ growth model, so the uncertainties from Poisson (small number) statistics on the predictions are small, even in the lower stellar mass bins. The scaling relations between galaxy and MBH properties evolve very marginally in this redshift range, and therefore no further modification is applied to the $z=0.1,0.2$ halos.

The SAMs provide the mass and mass accretion rates for MBHs at every snapshot. The simplest way to compute the luminosity is to convert the accretion rate into a bolometric luminosity with some efficiency $\epsilon_{f}$, and then use a bolometric correction factor to infer the luminosity in some specific energy/wavelength range. Like much of the literature, we assume the radiative efficiency of accretion $\epsilon_{f}$ = 0.1. The SAM is known to contain too little scatter (Ricarte & Natarajan, 2018a) since the AGN are essentially always on. In practice, a duty cycle would increase the scatter in the observed instantaneous luminosity of the black holes. Meanwhile, we know that relatively complete surveys like AMUSE (Gallo et al., 2008; Miller et al., 2012), with a limiting luminosity of $2\times 10^{38}\ \rm{erg\ s^{-1}}$ at 0.5-7.0 keV, find AGN in $\sim 40\%$ of galaxies with $10^{8}<M_{*}/M_{\odot}<10^{12}$. AMUSE is a remarkably uniform and complete survey of AGN in nearby galaxies performed with the Chandra X-ray telescope, observing 200 galaxies – some in the field, some in the Virgo galaxy cluster – to a limiting flux of $10^{38.3}\ \rm{erg\ s^{-1}}$, shown in Figure 1 with a dotted line. While more recent catalogs exist, they are compiled from the Chandra archive and are thus necessarily non-uniform in exposure (e.g. Gallo & Sesana, 2019). We measure the detected fraction of AGN in each SAM as:

	$\displaystyle w_{i}$	$\displaystyle=\frac{\Phi(M_{*,SAM})}{\Phi(M_{*,obs})}$		(1)
	$\displaystyle f_{d}$	$\displaystyle=\frac{\sum_{i=0}^{N_{a}}w_{i}}{\sum_{j=0}^{N_{g}}w_{j}}$		(2)

where $w_{i}$ is the relative abundance of the host galaxy stellar mass in the SAM compared to the observations. The sum in the numerator is over the $N_{a}$ galaxies with $L_{X}>0$ and in the denominator is over all $N_{g}$ galaxies. The weighting corrects for the different stellar mass functions of the SAMs and the observed sample. Both numerator and denominator include galaxies at $10^{8}<M_{*}/M_{\odot}<10^{12}$. The power-law and AGN-MS growth models predict the existence of many AGN below the AMUSE flux limit, shown as the dotted horizontal line in each panel.

If the SAMs are taken at face value, 75-90$\%$ of the galaxies in the AGN-MS and PL models would have had AGN detections in an AMUSE-like survey. We, therefore, add scatter to the $L_{X}$ in these SAMs until the detected fraction, as defined above, is $40\%$ for an AMUSE-like survey. These values are shown in Table 1, and the $L_{X}-M_{*}$ relations for the SAMs with the added scatter are shown in Figure 2.

Shen et al. (2020) measured the bolometric correction for the $0.5-2$ keV and $2-10$ keV X-ray bands and for the infrared B-band, as a function of bolometric luminosity. This corresponds to $1/k_{band}=L_{bol}/L_{band}=10-30$ in each of the X-ray bands and $4-7$ in the B-band, each with a scatter of about 50$\%$. For each black hole, we draw from the bolometric correction from a normal distribution with the mean and scatter based on its bolometric luminosity.

We focus first on how observational selection can impact the derived BHOF and explore the optimal future survey strategies that will enable potentially discriminating between seeding models. We assume that black hole masses are available from the full range of possible observational probes and documented empirical correlations - either from dynamical modeling of galaxy cores or other techniques. AGN, for instance, have been identified via their optical variability (Baldassare et al., 2018, 2020) or line ratio diagnostics (Polimera et al., 2022; Cid Fernandes et al., 2010; Kewley et al., 2013); masses are then measured via reverberation mapping (e.g. Peterson et al., 2004) or direct dynamical modeling (e.g. Walsh et al., 2013; Seth et al., 2014; Nguyen et al., 2017), and then used to construct scaling relations based on easily observable quantities like emission line strengths. Therefore, we proceed with the premise that we have in hand $M_{BH}$ measurements/upper limits for every single galaxy in our SAM survey sample. We first investigate the impact of the completeness of mass measurement on our inference of the BHOF.

Figure 3 shows the occupation fraction for different values of black hole mass completeness. If we assume in the best-case scenario that all black holes above $10^{3}\ \rm{M_{\odot}}$ can be detected, then the Pop III seeds produce a higher occupation fraction than DCBH models. This is particularly noticeable for galaxies with $M_{*}<10^{9}\ \rm{M_{\odot}}$. We note that this clearly reflects the efficiency and abundance of initial seeds, light seeds are more ubiquitous and are expected to form more efficiently compared to heavy seeds in the SAMs. If, however, we can only detect BHs more massive than $3\times 10^{4}\ \rm{M_{\odot}}$, then light seeds growing through the AGN-MS channel predicts $f_{\rm occ}=1$ for all stellar masses; for the other two growth channels, seed models become indistinguishable for $M_{*}<4\times 10^{8}\ \rm{M_{\odot}}$. If we can only find MBHs above $10^{6}\ \rm{M_{\odot}}$, information about seeding cannot be cleanly disentangled from accretion. These findings clearly show that the detectable present-day mass needs to be close to the initial assumed ”heavy” seed mass from high redshifts, namely, $10^{4-5}\ \rm{M_{\odot}}$, in order to be disentangled from the confounding effects of the accretion physics.

In summary, if the initial seeds have not grown much by any of the accretion modes or mergers, then their memory of seeding is retained in the BHOF as measured even at low redshift.

In practice, MBHs are usually detected when they are accreting and luminous. The scaling relations used to measure black hole masses regardless of wavelength deployed or method used, all suffer from large-intrinsic scatter (Woo & Urry, 2002) and are most reliable at X-ray wavelengths where the obscuration is the lowest. Therefore, we stick to X-ray observations for comparison with SAMs for the rest of this study. Figure 4 shows the fraction of galaxies with AGN detected at $0.5-2.0$ keV as a function of the survey detection limit. The BLQ prescription is ruled out because it predicts almost no AGN in dwarf galaxies. A relatively shallow survey of $10^{40}\ \rm{erg\ s^{-1}}$ would only be able to marginally distinguish between the AGN-MS and PL growth channels modes, with the former predicting 10-20$\%$ more detections in dwarf galaxies; if all AGN above $10^{38}\ \rm{erg\ s^{-1}}$ are found, the detected fraction is identical for both these growth channels and both growth channels. We need a survey that either directly detects, or can model the presence of, all AGN with luminosities down to $10^{36}\ \rm{erg\ s^{-1}}$ if the heavy and light seed models are to be clearly distinguished.

Figure 3 and 4 reveal that most of the differences between the models are apparent at the lowest halo masses. Firstly, this highlights the importance of using a parent sample of galaxies that is stellar mass complete down to $10^{8}\ \rm{M_{\odot}}$, pushing down an order of magnitude compared to current observational limits. Figure 5 therefore shows the occupation fraction of galaxies with stellar masses in the range $10^{8}-10^{10}\ \rm{M_{\odot}}$ as a function of the minimum detectable black hole mass. If all the information we have is the black hole and stellar masses, we can only tell different seed models apart if we find every MBH at least as massive as $10^{4}M_{\odot}$. If we can separately confirm that the growth channel is like the AGN-MS model, then we can distinguish seeds using all MBH above $10^{5}M_{\odot}$.

In practice, again, we almost always detect black holes when they are luminous. Figure 6 shows the fraction of galaxies with stellar mass $10^{8}<M_{*}/M_{\odot}<10^{10}$ with detectable AGN, as a function of survey limiting luminosity. Heavy and light seeds can be distinguished if the survey has a limiting luminosity better than $10^{37}\ \rm{erg\ s^{-1}}$.

Figure 7 shows the AGN bolometric luminosity function for the six models, with observations from Ueda et al. (2011) in black for comparison. The uncertainties incorporate 1000 draws from the uncertainty in the bolometric corrections; the BLQ model uncertainties are not shown because, due to its low fraction of active MBHs, has error bars down to 0 at all luminosities. The AGN-MS and PL models have diverging luminosity functions above $10^{42}\ \rm{erg\ s^{-1}}$, with the former model requiring more growth at late times to match the $M_{\bullet}-\sigma$ relation. The PL model is thoroughly consistent with the observations, although the AGN-MS model cannot be entirely ruled out due to the large scatter required to bring it into an agreement with the observed active fraction in the AMUSE survey.

As noted above, the BHOF depends on the minimum AGN luminosity that we can measure or model. Studies like Miller et al. (2015) compute the BHOF from the observed AGN luminosities and upper limits, using a model where the $L_{X}-M_{*}$ relation contains log-normal scatter. Therefore, this BHOF accounts for AGN even below the nominal luminosity limit of the survey, which is $10^{38}\ \rm{erg\ s^{-1}}$. Since they did not publish the scaling relations inferred by their Bayesian model, we re-measure it here. To recap, the original study assumed that the BHOF had the functional form of a sigmoid function:

$$f_{\rm occ}=\frac{1}{2}\left[1+\tanh\left({2.5^{|8.9-\log_{10}M_{*,0}|}\log_{10}\frac{M_{*}}{M_{*,0}}}\right)\right]$$

(3)

The tuning parameter $M_{*,0}$ dictates how rapidly the occupation fraction jumps from 0 to 1; in other words, a higher value for this parameter means that galaxies with lower stellar masses are less likely to host SMBH than their more massive counterparts. Put another way, a higher $M_{*,0}$ means a lower BHOF.

A Bayesian inference model then constrains $M_{*,0}$ simultaneously with the normalisation, slope and scatter of the $L_{X}-M_{*}$ relation:

$$\log(L_{X}/\rm erg\cdot s^{-1})=\alpha+\beta\times\log(M_{*}/M_{\odot})$$

(4)

This approach thus assumes 1) that the AGN duty cycle can be described in terms of the time spent below the detection limit if the $L_{X}-M_{*}$ scaling relation is a power law with log-normal scatter and 2) that we have no prior knowledge of this scaling relation. Again, Figure 2 shows the $L_{X}-M_{*}$ relation for each SAM, with points below the AMUSE detection threshold shown as translucent.

Right away, we can rule out the BLQ prescription as the primary growth channel, because this predicts that only $3-4\%$ of the galaxies at $z<0.2$ contain AGN, far below the $40\%$ detection fraction in AMUSE. This prescription was motivated from observations of the brightest quasars, the high luminosity end of the SDSS quasar LF. Since the abundance of bright quasars drops precipitously at low redshifts, the BLQ SAM is self-consistent in predicting low active fractions at low-z. The fact that this differs with observations of low-mass galaxies shows that this growth channel, which is bursty and merger-driven, cannot be important in dwarf galaxies. Next, we aim to distinguish between the AGN-MS and PL growth channels.

Both the power-law and AGN-MS, models predict very high (75-99$\%$) occupation fractions. From Figure 4, we know that we need to apply a luminosity cut to the SAMs before comparing them to observations. For the Miller et al. (2015) study, this is not the same as the limiting luminosity of the observations. Rather, the BHOF in that study includes all the AGN whose luminosity falls within log-normal scatter from the mean value of the power-law $L_{x}-M_{*}$ relation. Therefore, we start by measuring this $\sigma$.

The probability of a parameter set $\alpha,\beta,\sigma,M_{*,0}$ for detections is defined as:

$$\begin{split}p(\log(L_{X})|\log(M_{*}))=f_{occ}(M_{*},M_{*,0})\times\\ N(\log(L_{X})|\mu=(\alpha+\beta\times\log(M_{*}),\sigma=\sigma)\end{split}$$

(5)

where $N(x|\mu,\sigma)$ is the normal probability distribution function centered at $\mu$ with standard deviation $\sigma$. For a $n$-th non-detection, the probability that there is a black hole is given by the latent variable $I_{n}$, which in Miller et al. (2015) depended only on the scaling relation and occupation fraction.

$$\begin{split}L_{X,\rm Pred}&=\alpha+\beta\times\log(M_{*})\\ \Phi&=C(L_{\rm lim}-L_{X,\rm pred}/\sigma)\\ p_{BH}(L_{\rm lim},M_{*})&=\Phi/((1-f_{occ}(M_{*},M_{*,0})+f_{occ}\times\Phi)\end{split}$$

(6)

where $C$ is the cumulative probability of the normal distribution. This says that a non-detection means either that there is no black hole, or there is a black hole and it is on but the luminosity is below the detection limit $L_{\rm lim}$. The latent variable $I_{n}$ is then set to either 1 or 0 by drawing from a binomial distribution with $p_{\rm true}=p_{BH}$. For further details of the inference pipeline, we refer the reader to §2 of Miller et al. (2015).

As a test, we run the AMUSE sample through an MCMC pipeline four times, each time using a different scaling relation between $L_{X}-M_{*}$ in Table 1, inferred from Figure 1. The resulting posterior on $M_{*,0}$, the stellar mass at which occupation fraction is $50\%$, is shown in Figure 8. Surprisingly, these posteriors are almost indistinguishable from each other.

The best-fit value of $\sigma$ for AMUSE is 2; as seen in Table 1, this is very comparable to the total scatter required in each SAM to match the AMUSE detection fraction. Thus, the BHOF inferred assuming an $L_{X}-M_{*}$ relation with such a large scatter would capture all the AGN in the SAMs with both the PL and AGN-MS growth channels, and no luminosity cut needs to be applied on the SAM MBHs.

Finally, Figure 9 compares the AMUSE constraints applied to each of the four SAMs. The observations are clearly consistent with the heavy seed models growing via either the PL or AGN-MS channels. The light seed models, even with large scatter in their $L_{X}-M_{*}$ relation, cannot be reconciled with an X-ray observed sample like AMUSE.

A key result of our study is that distinguishing between MBH seeding and growth models relies crucially on detecting the least massive MBH (or least luminous AGN) in galaxy samples that are stellar mass complete down to at least $10^{8}\ \rm{M_{\odot}}$. We show here how such a survey could be constructed by following up on existing surveys.

The RESOLVE survey (Eckert et al., 2015) was designed to be complete in baryonic mass down to $10^{9.3}\ \rm{M_{\odot}}$, and thus to even lower stellar masses, within $z<0.025$ ($d<110$ Mpc). Figure 10 shows the spatial (left) and stellar mass (right) distribution of the RESOLVE galaxies, with the latter demonstrating stellar mass completeness above $5\times 10^{8}\rm{M_{\odot}}$. The orange histogram shows only the galaxies above $DEC=30^{\circ}$, where emission from the Milky Way is negligible, and even faint extra-galactic AGN could be detected. In principle, a survey could include galaxies at $|DEC|<30^{\circ}$ as long as they are in the direction opposite to the Galactic Center; we are being conservative in our recommendation here to account for telescope scheduling difficulties, which may make it difficult to point at extragalactic fields at all times. The RESOLVE survey contains $>10,000$ such galaxies.

Using this subset of galaxies as a parent sample, we create three different mock catalogs. Each assumes a different value of $M_{*,0}$ and thus a different BHOF. Where a MBH exists, we assign an X-ray luminosity drawn from the AMUSE observed $L_{X}-M_{*}$ relation, including the scatter. We pass it through the Bayesian inference pipeline described above, assuming survey luminosity limits of $10^{36}\ \rm{erg\ s^{-1}}$ and $10^{38}\ \rm{erg\ s^{-1}}$. Figure 11 shows the posterior distributions recovered for each of the mocks. Higher occupation fractions are harder to rule out than lower ones, due to the large scatter inherent in the $L_{X}-M_{*}$ relation. While broad, the mode of each distribution matches the input value for $L_{lim}=10^{36}\ \rm{erg\ s^{-1}}$. For shallower surveys, if the intrinsic BHOF is high, the survey would return a flat posterior on $M_{*,0}$ even for this high level of completeness. On the one hand, this suggests that the uncertainties obtained from the AMUSE survey may be biased low. On the other hand, it confirms that higher values of $M_{*,0}$ would have been ruled out with high confidence.

On the other hand, we have checked that the eFEDS survey, despite its uniform coverage overlapping with $\sim$500 times more galaxies than AMUSE and detection of several new AGN candidates in dwarf galaxies (Latimer et al., 2021a), yields flat posteriors on $M_{*,0}$, and that this situation is not expected to improve even with the deepest, polar regions of eRASS8. The key takeaway from this analysis is that survey depth, and the resulting completeness in stellar mass and MBH mass/luminosity, are much more powerful than survey volume when measuring the underlying BHOF.

A key challenge at low MBH masses, and therefore luminosities, is confusion with other sources such as Utra- or HyperLuminous X-ray sources (ULXs/HLXs) (Swartz et al., 2004, 2011), which are likely the high Eddington-ratio end of the accreting stellar black hole population (Stobbart et al., 2006; Sutton et al., 2013), although some studies argue that the brightest of these may indeed be intermediate-mass black holes (Miller et al., 2003; Bellovary et al., 2010; Sutton et al., 2012). ULXs may be differentiated from IMBHs by a lack of radio emission or a position offset from the galaxy center (Thygesen et al., 2023), although Bellovary et al. (2021) showed that in the shallow potentials of dwarf galaxies, MBHs could live off-center. Multi-band X-ray observations can enable spectral studies to further differentiate between MBHs and high-accreting stellar black holes, and multi-epoch observations may find different variability scales. Further contamination at the faintest ends may arise from populations of X-ray binaries (XRBs), although this can potentially be modeled via low-scatter scaling relations with the stellar mass and star-formation rate (Lehmer et al., 2010, 2019).

Whether the black hole selection is by mass (Figure 3) or luminosity (Figure 4), seeding models cannot be distinguished unless we find every MBH $\gtrsim 10^{4}\ \rm{M_{\odot}}$. This finding was also independently arrived at in recent work by Burke et al. (2022) using future variability surveys. These smallest MBHs have been recently discovered using optical variability (Baldassare et al., 2018) and line diagnostics Baldassare et al. (2015); Mezcua & Domínguez Sánchez (2020); Polimera et al. (2022). Emission line widths, such as $H-\alpha$, can yield black hole masses, although these can easily be confused with ongoing star formation.

MBHs are, however, notorious for the wide variety and variability of their Spectral Energy Distributions (SEDs). This means that we do not yet have a complete understanding of how the emission from the reprocessed material feeding a black hole is distributed at different wavelengths. Below, we summarize some of the opportunities and challenges at various wavelengths for future low redshift surveys.

In this idealised study, we have examined in detail the effects of seed and growth channels in isolation. It is likely that black holes may grow form via multiple channels. Spinoso et al. (2022) used a semi-analytic model to populate halos from the Millennium-II with both heavy and light seeds, and evolve their growth to z=0 using a two-phase prescription, corresponding to the observed quasar and radio modes of AGN. They find that with such a multiple-channel-seeding model, the slope of the black hole-stellar mass relation is flatter for dwarf galaxies than their more massive counterparts. We see this in Figure 1 in the models that grow with a power law ERDF, respective of the seeding model, although the difference is starker for heavy seeds. Interestingly, this heavy-seed, power-law ERDF model is also the one that is most consistent with the AMUSE observations. This is a promising sign for future X-ray surveys, suggesting that we may find more dwarf AGN by simply pushing to slightly lower luminosities.