Toward determining the number of observable supermassive black hole shadows

We present estimates for the number of shadow-resolved supermassive black hole (SMBH) systems that can be detected using radio interferometers, as a function of angular resolution, flux density sensitivity, and observing frequency. Accounting for the distribution of SMBHs across mass, redshift, and accretion rate, we use a new semi-analytic spectral energy distribution model to derive the number of SMBHs with detectable and optically thin horizon-scale emission. We demonstrate that (sub)millimeter interferometric observations with ${\sim}0.1$ $\mu$as resolution and ${\sim}1$ $\mu$Jy sensitivity could access ${>}10^6$ SMBH shadows. We then further decompose the shadow source counts into the number of black holes for which we could expect to observe the first- and second-order lensed photon rings. Accessing the bulk population of first-order photon rings requires ${\lesssim}2$ $\mu$as resolution and ${\lesssim}0.5$ mJy sensitivity, while doing the same for second-order photon rings requires ${\lesssim}0.1$ $\mu$as resolution and ${\lesssim}5$ $\mu$Jy sensitivity. Our model predicts that with modest improvements to sensitivity, as many as $\sim$5 additional horizon-resolved sources should become accessible to the current Event Horizon Telescope (EHT), while a next-generation EHT observing at 345 GHz should have access to ${\sim}$3 times as many sources. More generally, our results can help guide enhancements of current arrays and specifications for future interferometric experiments that aim to spatially resolve a large population of SMBH shadows or higher-order photon rings.


INTRODUCTION
The observations and resulting images of the supermassive black hole (SMBH) in the M87 galaxy by the Event Horizon Telescope (EHT) collaboration (Event Horizon Telescope Collaboration et al. 2019a,b,c,d,e,f) represent the first steps in a new field of spatially resolved horizon-scale studies of black holes. The emission from around the SMBH in M87 takes the form of a bright ring surrounding a darker central "shadow," as expected from simple models of spherical accretion (Falcke et al. 2000;Narayan et al. 2019). A wide variety of simulated images of black hole accretion flows show that this ring generically has a diameter that is comparable to the theoretical curve bounding the photon capture cross-section of the time-reversed black hole (Event Horizon Telescope Collaboration et al. 2019e,f). General relativity predicts that the boundary of this cross-section should take on a nearly circular shape with a diameter of approximately five times the Schwarzschild radius (Bardeen 1973), and that this diameter should depend only weakly (to within ±∼4%) on the black hole's spin and inclination (Takahashi 2004;Johannsen & Psaltis 2010). These properties permit spatially resolved observations to constrain the black hole mass using measurements of the shadow size; EHT observations of M87 yielded a ∼10% mass mea-surement via this approach (Event Horizon Telescope Collaboration et al. 2019f).
Though the EHT has focused its attention thus far on only those black holes with the largest angular sizes as seen from Earth, almost all massive galaxies are expected to host SMBHs (Magorrian et al. 1998;Kormendy & Ho 2013). As the EHT and future facilities improve upon the angular resolution and flux density sensitivity of the first M87 observations, more SMBH shadows -and their corresponding constraints on the black hole masses -will become observationally accessible. Though new black hole mass measurements are valuable for individual galaxy studies, questions about SMBH formation and growth mechanisms and the degree to which they co-evolve with their host galaxies are most effectively addressed using large statistical samples of precisely-measured SMBH masses (Volonteri 2010;Heckman & Best 2014). To this end, it is natural to ask what observational requirements would be necessary to access large numbers of SMBHs with spatially resolved shadows.
In addition to mass measurements, sufficiently highresolution observations of SMBHs can also provide unique access to the black hole spin and potentially other spacetime properties. Hidden within the ring of emission seen by the EHT is an unresolved series of approximately concentric "photon rings," formed by rays that execute increasingly many orbits about the black hole prior to escaping (Darwin 1959;Luminet 1979;Gralla et al. 2019;Johnson et al. 2020). Each higher order photon ring -enumerated by the number n of half-orbits that the constituent photon trajectories make around the black hole -is expected to have an exponentially narrower angular width on the sky than the previous order. The lowest-order (n = 0, corresponding to direct emission) photon ring is impacted by specific details of the accretion flow (e.g., stochastic turbulent structure) that complicate precise spacetime constraints, while the geometric properties of higher-order rings contain the same spacetime information while being exponentially less impacted by such "astrophysical" contamination. Furthermore, interferometric observations naturally decompose the emission by spatial scale, meaning that with fine enough angular resolution the signal from n > 0 will dominate the interferometric response in a time-averaged image (Johnson et al. 2020;Gelles et al. 2021).
The goal of this paper is to determine the number of SMBH shadows and low-order photon rings that could be observed as a function of angular resolution, flux density sensitivity, and observing frequency. We assume that such observations will be carried out using (sub)millimeter-wavelength interferometry, and we take 230 GHz to be a characteristic observing frequency when not otherwise specified. In Section 2 we describe our formalism and input assumptions, which we use to compute the number and distribution of SMBH shadows in the universe as seen from Earth. In Section 3 we modify these shadow counts to reflect the flux density response expected when observing with interferometers, and we further decompose the total source counts into contributions from systems for which we could observe the n ≥ 0, n ≥ 1, and n ≥ 2 photon rings. In Section 4 we discuss the implications of the source count distributions for current and future telescope specifications. We summarize and conclude in Section 5. Throughout this paper we assume a flat cosmology with Ω m = 0.3, Ω Λ = 0.7, and H 0 = 70 km s −1 Mpc −1 unless otherwise specified.

POPULATION SOURCE COUNTS
Our goal is to estimate the number of black hole shadows that we could hope to observe. Concretely, we would like to determine the number N (θ r , σ ν ) of SMBHs that satisfy the following three conditions: 1. The shadow of the black hole has an angular size larger than some resolution threshold θ r .
2. The flux density of the horizon-scale emission exceeds some sensitivity threshold σ ν .
3. The emitting plasma is optically thin.
The first of these criteria is set primarily by the mass of and distance to the black hole, while the second two also depend on the mass accretion rate and the physical conditions in the accretion flow. The third criterion exists to ensure that we could identify a black hole shadow as such; i.e., an optically thick emission region could obscure the shadow even if the angular resolution and sensitivity would otherwise make it accessible.

Overview of strategy
Our strategy for determining N (θ r , σ ν ) starts by considering the global distribution of SMBHs as a function of mass M and redshift z, to which we then sequentially apply the above three criteria to narrow down the number of potentially detectable sources. The distribution Φ(M, z) is described by the black hole mass function (BHMF), which we discuss in Section 2.2.
For a given SMBH mass M and redshift z, applying our first criterion -i.e., that the angular shadow size ϑ is larger than some resolution θ r -amounts to requiring that the black hole mass exceed some minimum mass m 0 (z). A black hole of mass M situated at an angular diameter distance D A has an angular shadow size that is given by where R S is the Schwarzschild radius and the numerical prefactor √ 27 is determined by the shadow diameter for a Schwarzschild black hole (Hilbert 1917;Bardeen 1973). At a particular redshift z, the condition ϑ ≥ θ r corresponds to where m 0 (θ r ) is the critical mass for which a SMBH at redshift z has a shadow with angular size θ r , and where we have cast the expression in terms of the comoving distance, D(z) = (1 + z)D A (see also Bisnovatyi-Kogan & Tsupko 2018). Applying our second condition -i.e., that the flux density S ν be greater than some threshold σ ν -requires knowing the distribution P (S ν |M, z) of flux densities for a SMBH of mass M at redshift z. The flux density S ν (ν 0 ) observed at a frequency ν 0 is related to the emitted luminosity density L ν by (Peacock 1999), Here, L ν ([1 + z]ν 0 ) denotes the luminosity density evaluated at the redshifted frequency (1 + z)ν 0 , and we have assumed that the emission is isotropic. 1 L ν is determined by the spectral energy distribution (SED) of the source, which we model as described in Section 2.3 (with more comprehensive details provided in Appendix A). Within our SED model, L ν depends not only on the mass of the SMBH but also on its mass accretion ratė M , which we cast in terms of the Eddington ratio λ, Here,Ṁ Edd ≡ L Edd /ηc 2 is the Eddington mass accretion rate, and η is a nominal radiative efficiency that relates Edd to the Eddington luminosity L Edd ; for this paper, we take the radiative efficiency to be η = 0.1 (e.g., Yuan & Narayan 2014). Determining P (S ν |M, z) thus further requires knowledge of the Eddington ratio distribution function (ERDF), which we describe in Section 2.4.
Applying our third condition -i.e., that the horizonscale emission be optically thin at the observing frequency ν 0 -can also be achieved using our SED model, which provides an optical depth prediction for a SMBH with any given M , λ, and ν 0 . Practically, we can absorb this condition into the definition of the flux density distribution by considering only those systems that are optically thin, i.e., by determining P (S ν |M, z, τ ≤ 1). The fraction f (σ ν ) of SMBHs for which we could expect to detect the horizon-scale emission is then given by where σ ν is some specified sensitivity threshold. Combining all three criteria, we can compute the source counts expected for any choice of θ r and σ ν by integrating the global distribution over mass and redshift, Many of the results presented in this paper are derived from evaluating Equation 7. When computing this integral, we must keep in mind that m 0 (θ r ) is a function of z and that f (σ ν ) is a function of both z and M . Figure 1 illustrates the procedure we follow to determine N (θ r , σ ν ) for an example set of angular resolution and flux density thresholds (in this case, θ r = 1 µas and σ ν = 10 −5 Jy).

Black hole mass function
Any evaluation of Equation 7 requires a choice of BHMF, which commonly takes the form where dN is the number of SMBHs in the mass range (M, M + dM ) and the comoving volume range (V, V + dV ). 2 For our purposes it is more useful to work with Φ(M, z), the number of black holes in the redshift range (z, z + dz) (see Equation 1), which is related to Φ by BHMF provides Φ( , ) Apply flux density and optical depth criteria Apply angular resolution criterion and integrate Repeat for many values of Figure 1. Flowchart illustrating the strategy for determining N (θr, σν ) (see Section 2.1) using an example case of θr = 1 µas and σν = 10 −5 Jy. Panels (a), (b), and (c) show the three primary inputs: the black hole mass function, the SED model, and the Eddington ratio distribution function, respectively. The BHMF provides the global distribution of SMBHs as a function of M and z (see Section 2.2), the SED model predicts the emitted flux density and optical depth for every M , λ, and z (see Section 2.3), and the ERDF provides the distribution of Eddington ratios λ (see Section 2.4). In panel (d) the ERDF and the SED model are used to determine the fraction f (σν ) of sources that simultaneously have flux densities exceeding σν (left plot in the panel) and are optically thin (i.e., τ < 1; right plot in the panel); in both plots, darker colors indicate a larger fraction. The combined fraction, as a function of M and z, is then used in panel (e) to modify the global SMBH distribution from the BHMF. In panel (f) we further apply the requirement that the angular shadow size exceed θr, which can be cast as a minimum mass m0(θr) at every z; N (θr, σν ) is then determined by integrating over the region outside of the gray shaded area. Finally, panel (g) illustrates that this procedure can be repeated for many other values of both θr and σν (see Section 2.5).
Estimating the BHMF from observations is difficult because astronomical surveys are inevitably incomplete in ways that impose poorly-known selection functions on the SMBH count in any mass bin, and because there are currently no SMBH mass measurement techniques that are both precise and broadly-applicable (Kelly & Merloni 2012). Many variants of the BHMF thus exist in the literature (e.g., Salucci et al. 1999;Aller & Richstone 2002;Marconi et al. 2004;Greene & Ho 2007;Lauer et al. 2007;Natarajan & Treister 2009;Kelly & Shen 2013). Recognizing that no single one of these BHMFs is likely to be uniquely correct, in this paper we consider two different BHMF prescriptions -which we will refer to as our "lower" and "upper" BHMFs -that aim to capture a reasonable range of possibilities.
We take as our lower BHMF the phenomenological model developed by Shankar et al. (2009) and shown in the left panel of Figure 2. This BHMF is evolved selfconsistently forward in time within a continuity equation formalism (Cavaliere et al. 1971;Small & Blandford 1992) tuned to match an estimate of the bolometric AGN luminosity function based primarily on the Xray observations compiled by Ueda et al. (2003). The Shankar et al. (2009) BHMF is a function of both M and z, covering SMBH masses in the range 10 5 -10 9.5 M and redshifts in the range 0-6. To account for the known existence of SMBHs with masses exceeding 10 9.5 M (e.g., Event Horizon Telescope Collaboration et al. 2019f), we extrapolate the BHMF using a power law with an exponential cutoff, The index and normalization of the power law are determined for every z by fitting the BHMF values between 10 9 M and 10 9.5 M . Natarajan & Treister (2009) argued on empirical and theoretical grounds for the existence of an upper mass limit for SMBHs at every cosmic epoch. First, using a physical argument based on selfregulation, they showed that when the accretion energy of a growing SMBH back-reacts with the gas flow and exceeds the binding energy of the feeding disk, it leads to the BH stunting its own growth and results in an upper limit for its mass (see also King 2016). Empirically, such a limit is expected from the observed SMBH mass -bulge luminosity relation when the relation is extrapolated to the bulge luminosities of bright central galaxies in clusters (Magorrian et al. 1998). Natarajan & Treister (2009) showed that consistency between the optical and X-ray BHMFs requires an upper mass limit for local SMBHs that is on the order of ∼10 10 M . Calibrating their estimates using the more recent observational measurements of the ∼1.7 × 10 10 M SMBH in NGC 1600 (Thomas et al. 2016), we determine an exponential cutoff mass of 3.5 × 10 10 M . The extrapolated portion of the BHMF is plotted using dashed lines in the left panel of Figure 2.
As a counterpart to the model-based lower BHMF, we also consider an upper BHMF derived empirically using the UniverseMachine stellar mass function (SMF) from Behroozi et al. (2019). The UniverseMachine SMF is constructed as part of a comprehensive model for galaxy growth spanning redshifts 0 ≤ z ≤ 10 and accommodating many observational constraints, including among them a number of observational SMFs determined in various bands (Baldry et al. 2012;Ilbert et al. 2013;Moustakas et al. 2013;Muzzin et al. 2013;Tomczak et al. 2014;Song et al. 2016). From the Uni-verseMachine SMF, we convert from stellar mass M * to SMBH mass M using the scaling law from Kormendy & Ho (2013) as done in Ricarte & Natarajan (2018), After converting from stellar to SMBH mass, we convolve the SMBH mass distributions with a Gaussian kernel with a 0.3-dex FWHM to account for the intrinsic scatter in the scaling relations. The resulting upper BHMF is shown in the center panel of Figure 2. Relative to the lower BHMF, the upper BHMF predicts systematically more SMBHs at low to intermediate redshifts (i.e., z 3) and at all masses, though at the highest redshifts the lower BHMF predicts more SMBHs with M 10 9.5 M (see right panel of Figure 2). The low-redshift behavior of the lower BHMF agrees well with a BHMF derived from the UniverseMachine SMF using the McConnell & Ma (2013) scaling relation (see also Saglia et al. 2016) 3 . To remain conservative in our estimates, throughout this paper we treat the lower BHMF as our fiducial case and use it for all computations and figures unless otherwise specified; we use the upper BHMF primarily to determine plausible uncertainty ranges for computed values. For this paper, we treat both BHMFs as being nonzero only in the range 0 ≤ z ≤ 6 and 10 5 ≤ M ≤ 10 11 M .
For the analyses carried out in this paper, the highmass end of the BHMF is most important. To assess the fidelity of the high-mass end of the lower and upper BHMFs, we compare their predictions against the number of known massive SMBHs in the local universe. In this regard, the MASSIVE galaxy survey provides a convenient comparison point because it is a volume-limited survey targeting massive early-type galaxies with stellar masses above 10 11.5 M within a distance of 108 Mpc, or z ≈ 0.02 (Ma et al. 2014). To date, four 4 SMBHs in this volume have dynamically measured masses at or above M87's M ≥ 6.5 × 10 9 M : M87 (Event Horizon Telescope Collaboration et al. 2019f), NGC 1600 (Thomas et al. 2016), NGC 3842, andNGC 4889 (Mc-Connell et al. 2011). Our lower and upper BHMFs predict that the number of SMBHs within z ≤ 0.02 and M > 6.5 × 10 9 M should be ∼5 and ∼29, respectively, which are consistent with the MASSIVE survey results. The specific behavior of the BHMF at low masses is less important because these black holes do not contribute significantly at the angular resolutions and flux densities of most interest for this paper.

Spectral energy distribution model
Given the global distribution of SMBHs across mass and redshift, Equation 7 selects only the fraction f (σ ν ) that have optically thin emission with a flux density that exceeds the sensitivity threshold. This fraction is defined in Equation 6, and it results from integrating over the distribution of flux densities P (S ν |M, z) at a given M and z. The first piece of information we need to compute this integral is an SED model, which will permit us to determine the flux density S ν (λ|M, z) corresponding to a particular choice of Eddington ratio, black hole mass, and redshift, and also to assess when the observed emission will be optically thin.
Observational constraints on SMBH growth indicate that SMBHs spend the majority of their time accreting at well below the Eddington rate (Hopkins et al. 2006). At these low accretion rates, the material in the vicinity of the black hole is thought to follow the advectiondominated accretion flow (ADAF) solution to the hydrodynamic equations describing viscous and differentiallyrotating flows around black holes (Narayan & Yi 1995a;Narayan et al. 1998;Yuan & Narayan 2014). An ADAF accretion disk has a two-temperature structure in which the ion temperature is greater than the electron temperature. The electrons are able to cool via a combination of synchrotron, bremsstrahlung, and inverse Compton radiation, which together define the SED for the observed emission.
For SMBHs observed in the radio to submillimeter wavelength range, as relevant for this work, the SED is dominated by synchrotron and Compton emission. Ma-hadevan (1997, hereafter M97) provides a convenient formalism for computing the gross spectral properties of an ADAF system given a black hole mass M and accretion rateṀ (see also Narayan & Yi 1994, 1995b. We use a modified version of the M97 formalism for the SED models in this paper, and Appendix A provides a detailed description of our updated model. We note that this SED model only considers emission from the accretion flow, and it does not incorporate a jet component. Our SED model provides an estimate of the emitted luminosity density L ν as a function of frequency for any input values of M and λ. Given a particular redshift z, we convert L ν to S ν using Equation 4. We determine whether the system is optically thin by comparing the rest-frame observing frequency, (1 + z)ν 0 , to the peak synchrotron frequency in the source, ν p (see Equation A.16). So long as (1 + z)ν 0 ≥ ν p , we consider the system to be optically thin.

Eddington ratio distribution function
The last piece of information we need to compute the integral in Equation 6 is an ERDF, which provides a probabilistic description of what fraction of SMBHs should be accreting at any particular Eddington rate λ. In this paper, we consider every SMBH to be active at some level, rather than considering the accretion to have only binary "on" and "off" states. We thus dispense with the notion of a "duty cycle" often adopted for AGN (or equivalently, we take the duty cycle to be unity), and we instead work exclusively in terms of an Eddington ratio distribution function (e.g., Merloni & Heinz 2008) to account for the differences in accretion rates.
There is emerging evidence that luminous ("Type 1"; unobscured; λ 10 −2 ) and low-luminosity ("Type 2"; obscured; λ 10 −2 ) AGN follow different distributions (Kauffmann & Heckman 2009;Trump et al. 2011;Weigel et al. 2017). Though the ERDF for luminous AGN appears to be consistent with a log-normal distribution (Lusso et al. 2012), there is no clear consensus in the literature on a specific form for the ERDF of low-lumminosity AGN (LLAGN). Different authors have used variants that include a power-law (Aird et al. 2012;Bongiorno et al. 2012), a Schechter function (Hopkins & Hernquist 2009;Cao 2010;Hickox et al. 2014), and a log-normal (Kauffmann & Heckman 2009;Conroy & White 2013). Additionally, while there seems to be broad agreement on a power-law behavior towards low Eddingtion ratios in the local Universe (i.e., z 1), few observational constraints currently exist for the ERDF of LLAGN at z 1. We proceed with a form for the ERDF adapted from the analytic prescription used by Tucci & Volonteri (2017) and updated using the more recent measurements from Aird et al. (2018). For their ERDF, Tucci & Volonteri (2017) used a Schechter function with an exponential cutoff value of λ = 1.5, but for our purposes (i.e., LLAGN with λ 1) only the power-law component of the ERDF is relevant. Furthermore, the LLAGN portion of the ERDF from Tucci & Volonteri (2017) was constructed to match the low-redshift behavior from Hopkins & Hernquist (2009), Kauffmann &Heckman (2009), andAird et al. (2012). None of these previous papers included observational constraints for AGN accreting below λ ≈ 10 −5 . To avoid the strong dependence on the low-end cutoff that comes from continuing the power law to arbitrarily small values, we posit instead that the distribution breaks (as in, e.g., Weigel et al. 2017). Specifically, we modify the power-law ERDF from Tucci & Volonteri (2017) such that it flattens out for Eddington ratios smaller than some value λ 0 . That is, we have where P (λ) is the probability density per unit logarithmic interval in λ, λ min and λ max are the lowest and highest permitted values, and the coefficient A is constructed such that the distribution integrates to unity: In this paper, we use values of λ min = 10 −10 , λ 0 = 10 −5 , and λ max = 10 −2 (see Section A.4). In addition to permitting the power-law index α to evolve with redshift, we also allow for additional evolution with SMBH mass, .
(14) Here, a(M ) encodes the mass dependence of the powerlaw index. Though there is some prior observational evidence indicating that the ERDF is approximately independent of SMBH mass (Kauffmann & Heckman 2009;Kelly & Shen 2013;Weigel et al. 2017), recent measurements by Aird et al. (2018) found that more massive SMBHs tend to be accreting at higher rates. We thus treat a(M ) as being essentially bimodal, with low-mass SMBHs having one power-law index value and high-mass SMBHs having another, and we use a logistic function to smoothly vary a(M ) between these two extremes, Here, a lo describes the power-law index at small masses, a hi describes the power-law index at large masses, M 0 denotes the midpoint mass, and ∆ is the logistic width in log(M ) that controls how quickly the transition from the low-mass regime to the high-mass regime occurs. We determine the values of these four parameters by fitting Equation 12 to the Aird et al. (2018) measurements; our fitting procedure is described in Appendix B. We find best-fit values of a lo = 0.55, a hi = 0.20, log(M 0 ) = 7.5, and ∆ = 0.3, and the resulting ERDF is shown in Figure 3. Equation 12 defines the probability P (λ) per unit log(λ) for any particular SMBH to be accreting at the rate λ. Given some specified M and z, we determine the probability P (S ν |M, z) by numerically sampling from P (λ) and using our SED model (see Section 2.3) to associate each sample with a particular S ν . Efficient sampling of P (λ) can be achieved by transforming a random variable x that is distributed according to a unit uniform distribution through the inverse cumulative distribution function (CDF) of Equation 12. This inverse CDF is given by which we can use to generate random samples distributed according to Equation 12. The associated distribution of S ν provides an estimate of P (S ν |M, z), which we then integrate per Equation 6 for the purposes of evaluating Equation 7.

The number of black hole shadows
Putting it all together, Figure 4 shows the result of evaluating Equation 7 over a range of values for both the angular resolution threshold θ r and the flux density sensitivity σ ν at an observing frequency of ν 0 = 230 GHz. The top panel shows the source counts predicted without imposing the optical depth condition, while the bottom panel restricts the sources to those that satisfy τ ≤ 1 (see Equation 6). Each point in both panels of Figure 4 is computed from an integral over the remaining (M, z) space. These plots thus represent an observation-independent prediction about the character of the SMBH population; namely, how many SMBHs are expected to have angular shadow sizes in excess of θ r , horizon-scale flux densities at 230 GHz greater than σ ν , and (in the case of the bottom panel) an optically thin accretion flow. An approximate analytic description of the resulting N (θ r , σ ν ) is provided in Appendix C.
The two panels of Figure 5 show the behavior of N (θ r , σ ν ) in the limit as σ ν = 0 (left panel) and θ r = 0 (right panel); these limits correspond approximately to one-dimensional slices through the top panel of Figure 4 along the horizontal and vertical axes, respectively. The black curve in the left panel shows N (ϑ > θ r , σ ν = 0), while the colored curves show the contribution from SMBHs in different mass ranges. At large θ r we see that the source counts follow the N ∝ θ −3 r behavior expected from simple volume scaling. The upturn around θ r ≈ 1 µas occurs because this is the resolution threshold below which the most massive SMBHs can be seen at any redshift (because of the turnover in angular diameter distance at z ≈ 1.6), and the re-flattening at smaller θ r is caused by the finite redshift coverage of the BHMF. The black curve in the right panel shows N (θ r = 0, S ν > σ ν ), while the colored curves again split out the contribution by SMBH mass. Throughout most of the space we see the source counts climbing volumetrically as the flux density decreases, following 10 −2 10 −1 10 0 10 1 10 2 θ r (µas) The solid contours start with the thick contour indicating a count of N = 1 and then increase by factors of 10 towards the lower left, while the dashed contours each decrease by a factor of ten towards the upper right. Bottom: Same as the top panel, but with the additional restriction that the sources must be optically thin.
Cosmological effects become noticeable at the lowest σ ν values, where the curve starts flattening out owing to a combination of the luminosity distance increasing more rapidly as well as the finite redshift coverage of the BHMF.

INTERFEROMETRIC SOURCE COUNTS
The analysis performed in the previous section predicts the source counts corresponding to the population of SMBHs that adhere to the three criteria specified at the beginning of Section 2. We now aim to estimate a subtly different quantity: the number of shadowresolved sources that could be observed by a telescope with angular resolution θ r and flux density sensitivity σ ν . This conceptual distinction is relevant because the telescopes that we expect to be carrying out spatially resolved studies of black hole shadows in the foreseeable future are radio interferometers. While the source counting analysis performed in Section 2 uses the SED model detailed in Appendix A to determine the flux density expected from any particular SMBH, this SED model only provides an estimate for the total (i.e., spatially integrated) horizon-scale flux density. However, an interferometric baseline is only sensitive to flux on specific spatial scales, determined by the length of the baseline and the wavelength of light being observed. Though in this paper we do not explore specific methods for estimating shadow diameters, sparse interferometric observations have previously been used to constrain the shadow diameter for M87 under the assumption that the source is ring-like (Doeleman et al. 2012;Wielgus et al. 2020). In this section, we thus investigate the prospects for detecting SMBH shadows on an individual interferometric baseline.

Flux density seen by a single baseline
We base our expectations for the horizon-scale emission structure from a SMBH on the observational and theoretical understanding of the M87 system. Johnson et al. (2020) provide an approximate analytic expression for the expected flux density of the photon ring emission as a function of baseline length for optically thin emission, which we adapt to take the following form: Here, S 0 is the total flux density (i.e., the value provided by the SED model, given the redshift of the SMBH), ϑ is the angular diameter of the photon ring (which for Left: The black curve shows the total number of SMBHs with shadows larger than some threshold angular resolution θr as a function of that threshold; this curve approximately corresponds to a horizontal cut through the bottom part of the top panel of Figure 4. The upper axis indicates the minimum mass of a black hole for which the corresponding angular resolution would permit that black hole to be spatially resolved at any redshift. Right: The black curve shows the total number of SMBHs with horizon-scale flux densities larger than some threshold value σν as a function of that threshold; this curve approximately corresponds to a vertical cut through the left part of the top panel of Figure 4. In both panels the source counts for different choices of black hole mass binning are shown as colored curves. our purposes is given by Equation 2), u is the length of the baseline in units of wavelengths, W 0 is the FWHM angular thickness of the lowest-order (i.e., n = 0) photon ring, and η = 1−e −π is a normalizing prefactor. We assume W 0 = ϑ/5 (Event Horizon Telescope Collaboration et al. 2019e,f). Equation 17 is shown as the gray curve in Figure 6.
On long baselines (i.e., u 1/ϑ), the bandwidthaveraged flux density will be given bȳ which is smaller by a factor 2/π than the envelope of Equation 17 Figure 7. Unlike in Figure 4, the source counts no longer monotonically increase as angular resolution improves (i.e., as θ r decreases), because Equation 19 ensures that longer baselines see lower flux den-sities from any given SMBH. An analytic approximation for the resulting N (θ r , σ ν ) is provided in Appendix C. Figure 8 shows the integrand from Equation 7 plotted over the domain of integration for several values of θ r and σ ν , providing the distribution of observable SMBHs as a function of M and z. The number of objects generally increases with increasing redshift (at fixed mass) and with decreasing mass (at fixed redshift), though the density peaks at z ≈ 2 for the smallest values of σ ν . For certain configurations, such as θ r = 0.1 µas and σ ν = 10 −5 Jy, the impact of Equation 19 is visually apparent as a lack of monotonicity in the source counts with increasing redshift (at fixed mass). This behavior reflects the fact that a fixed baseline becomes sensitive to emission from larger spatial scales around a particular SMBH as that SMBH is moved to larger distances; i.e.,S increases as ϑ decreases. On certain intervals in z, this flux increase associated with smaller ϑ is more than sufficient to compensate for the flux decrease associated with the increased distance to the SMBH. shadow n = 1 n = 2 Figure 6. Fraction of the total source flux density that can be detected on long baselines for the photon ring model described in Section 3.1, shown here for a ϑ = 40 µas diameter. The gray curve shows Equation 17, and the dashed black line tracks the envelope of this function. The red curve shows a running average of the gray curve across a 2% fractional observing bandwidth, and the dashed blue curve shows Equation 19. The vertical cyan lines show the resolution criteria used for the shadow (θr = ϑ) and for the first two orders of photon ring (θr = 2wn−1).

The expression in Equation 19
for the horizon-scale flux density contains contributions from all orders of photon rings, and in Figure 6 we can see that rings of different order are expected to dominate the observed flux density on different baseline length intervals. Depending on the value of u relative to 1/ϑ, a telescope may thus be primarily sensitive to emission from photon rings with n > 0. To determine the number of sources from which we expect to be able to detect higher-order photon rings, we can decompose the total source counts into bins corresponding to which order of photon ring dominates the emission.
We take as our resolution requirement to "see" the nth sub-ring that θ r ≤ 2w n−1 , where w n = W n / 8 ln(2) is the Gaussian width corresponding to the FWHM W n (Equation 18). This angular resolution requirement can be re-cast as a mass threshold m n for a given redshift, analogous to Equation 3; for n > 0, we have 10 −2 10 −1 10 0 10 1 10 2 θ r (µas) where m 0 is defined in Equation 3. Figure 6 marks the n > 0 and n > 1 resolution thresholds using vertical dashed cyan lines. To ensure that the emission is optically thin enough to see down to the nth sub-ring, we further impose a more stringent condition on the optical depth of τ ≤ τ n = 1 n + 1 .
By replacing the lower mass limit m 0 in Equation 7 with m n , and by replacing the τ ≤ 1 condition in Equation 6 with τ ≤ τ n , we can compute the source counts associated with objects for which a photon ring of order n or greater is detectable. Figure 9 shows these source counts for the first three orders of photon ring at observing frequencies of 86 GHz, 230 GHz, 345 GHz, and 690 GHz, corresponding to standard atmospheric transmission windows (Thompson et al. 2017). At each observing frequency we see qualitatively similar behavior: the source counts corresponding to the higher-order photon rings look approximately like scaled-down versions of the n ≥ 0 counts. For each additional order, the same source count value is achieved at an angular resolution threshold that is approximately ∼20 times finer and a sensitivity thresh- , showing the distribution of the number of shadow-resolved and optically thin SMBHs that can be seen by a single baseline as a function of redshift and black hole mass. Each panel shows a different choice of θr and σν , and all panels assume an observing frequency of 230 GHz. The total number of black holes, integrated over M and z, is given in the lower left-hand corner of each panel. The colorscale maps to the logarithm of the source number density (i.e., the number of sources per unit logarithmic interval in M and z), and the black contours enclose 50%, 90%, 99%, and 99.9% of the total source count. All panels share the same horizontal and vertical axis ranges, which are explicitly labeled in the bottom-left panel.
old that is approximately ∼100 times fainter than was necessary at the previous order. The angular resolution increment is associated with the factor e −π ≈ 1/23 in Equation 18 that sets the angular size ratio between consecutive photon rings. The flux density increment comes from a combination of a similar exponential suppression factor (a factor e −π from the summand of Equation 19) as well as the fact that the flux density profile is being observed on baselines that are typically a factor of e π longer, thereby incurring an additional flux density factor of e −π/2 ≈ 1/5 from the u −1/2 proportionality in Equation 19.
The evolution of the source counts with frequency primarily affects the required sensitivity, with higherfrequency observations achieving the same source counts at a higher value of σ ν than lower-frequency observations. The flux density threshold required to detect a particular number of objects is about one order of magnitude smaller at 86 GHz than at 690 GHz; i.e., about an order of magnitude better sensitivity -in terms of Jy -is required at 86 GHz than at 690 GHz. The angular resolution requirement does not show substantial evolution with frequency across this range.

The impact of baseline projection
The analysis presented in this section thus far has assumed that an interferometric baseline can observe the entire sky with the same angular resolution. However, in reality any physical baseline between two stations will have a different projected length as seen from different locations in the sky. The resolving power of the baseline will thus be a function of source location on the sky, which means that the number of black hole shadows a baseline can detect per unit solid angle will also vary across the sky.
For a particular baseline, we can define a spherical coordinate system (θ, φ) such that θ is a polar angle measured from the axis defined by the baseline orientation and φ is measured azimuthally around this axis. θ r (θ) is then the effective angular resolution of the baseline when projected toward a source at a sky position with polar angle θ, where b is the baseline length in units of the observing wavelength and θ r,0 = 1/b is the angular resolution achieved when θ = π/2 (i.e., the finest resolution achievable by the baseline). Denoting the number density of sources per unit solid angle as d 2 N dΩ , we can express the total number of sources observable by this baseline as where we have explicitly indicated that the number density is a function of the angular resolution, θ r (θ, b), and we have assumed that sources are distributed isotropically on the sky such that there is no φ dependence.
The integral is carried out over the solid angle Ω vis on the sky that is visible to the baseline.
To illustrate the impact of this geometric effect on source counts, we consider the concrete example of an interferometric baseline formed between two spacebased antennas, each of which can see the entire sky. In this case, the function d 2 N/dΩ is given simply by N (θ r , σ ν )/4π, and the domain of integration for Equation 23 will be all (θ, φ); Figure 10 shows the result of this evaluation. Relative to the source counts in Figure 7, at large values of θ r,0 (e.g., ∼20 µas) the source counts in Figure 10 are reduced because some fraction of the sky is not observed with sufficient angular resolution to see SMBHs with shadow sizes that are close to θ r,0 . The magnitude of this reduction is modest, amounting to a factor of 3π/16 ≈ 0.59 for uniformly distributed sources in flat space (see Equation D.3 with α = 4). However, a much more pronounced impact can be seen in the region of fine angular resolution and poor sensitivity (e.g., the region around θ r,0 ≈ 10 −1 µas and σ ν ≈ 10 −4 Jy), where the source counts in Figure 10 are significantly increased relative to Figure 7. This difference arises because N (θ r , σ ν ) climbs rapidly toward larger θ r in this region, and so the coarser angular resolutions arising from baseline projection provides access to many SMBHs that a baseline with a fixed angular resolution of θ r,0 across the entire sky would be unable to see. In this region of the (θ r,0 , σ ν ) space, the impact of baseline projection is to increase the accessible number of SMBH shadows by several orders of magnitude.
While Equation 23 provides the source counts appropriate for a fixed baseline, in real-world arrays the baseline will typically be changing orientation with time. For instance, a spaceborne antenna forming a baseline with another antenna situated on the Earth would execute a complete revolution once every orbital period, as observed by a distant source. One effect of this rotation is to make a larger fraction of the sky observable with the finest resolution than would otherwise be possible with just the instantaneous configuration, up to a unit fraction if both stations are spaceborne and thus can view the entire sky. The net impact of rotating the baseline 10 −2 10 −1 10 0 10 1 10 2 θ r (µas) . Each panel shows a plot analogous to that in Figure 7, but decomposed into the number of sources for which we could expect to detect the first three orders of photon ring. Each row shows this decomposition for one of four observing frequencies, with the frequency labeled in the upper left-hand corner of each panel. For each row of panels, the left panel shows the number of sources for which we could detect any order of photon ring at all, while the center and right panels show the number of sources for which we could detect the first-and second-order photon rings, respectively. In each panel, the drawn contours mark the same source count values as those in Figure 4. All panels share the same horizontal and vertical axis ranges, which are explicitly labeled in the bottom-left panel. is to bring more SMBH shadows into view than would be accessible by a static baseline. Appendix D provides a more detailed exposition of the sampling behavior of such a baseline as it rotates.

DISCUSSION
Our general strategy for carrying out the various source counting analyses presented in this paper is laid out in Section 2.1 and illustrated in Figure 1. To recap: • We start with the BHMF, which describes the global distribution Φ(M, z) of SMBHs across mass and redshift.
• Using our SED model and a prescription for the distribution of SMBH accretion rates (i.e., the ERDF), we determine the fraction f (σ ν ) of objects for which the horizon-scale emission is both optically thin and has either a total flux density (in Section 2) or a resolved flux density (in Section 3) exceeding some threshold σ ν .
• We then integrate the product f (σ ν )Φ(M, z) over M and z, excluding objects with shadow sizes smaller than some angular resolution threshold θ r (see Equation 7). . Power-law fits to the ridge-line in N (θr, σν )defined as the location of maximum σν for every fixed Nfor four different observing frequencies; this ridge-line can be seen as the turnover in the contours in Figure 7. Configurations of (θr, σν ) that fall below and to the right of the ridge-line can most effectively increase N by improving angular resolution (i.e., by decreasing θr), while configurations that fall above and to the right of the ridge-line can increase N by improving sensitivity (decreasing σν ) or by increasing θr. For reference, we mark the approximate specifications of the EHT (i.e., θr = 20 µas, σν between 1 mJy and 0.5 Jy) by a shaded gray region.
ing relations expected if the number of sources grows with the accessed volume (see Appendix C); for example, ∼hundreds of sources are predicted to be resolvable with an angular resolution of ∼1 µas and a flux density sensitivity of ∼1 mJy. For interferometric observations, we find that the number of detectable SMBH shadows generally increases as the angular resolution θ r and sensitivity σ ν improve, but that the gradient of N (θ r , σ ν ) changes orientation throughout the parameter space (see Figure 7). At large θ r and small σ ν , the source counts increase exclusively toward smaller θ r ; at large σ ν and small θ r , the source counts increase both toward smaller σ ν and toward larger θ r . The gradient changes orientation from pointing primarily toward smaller θ r to pointing primarily toward smaller σ ν around a ridge-line in the (θ r , σ ν ) space that approximately follows a power law σ ν ∝ θ 2.2 r ; Figure 11 shows power-law fits to this ridge-line for four different observing frequencies. At an observing frequency of 230 GHz, we find a best-fit power law of This expression can be used to estimate the angular resolution and sensitivity corresponding to an effective "Pareto front" 6 in source counts, whereby (θ r , σ ν ) pairs living on this curve are in some sense maximally economical. That is, to access the same number of sources using a different set of (θ r , σ ν ) would require improving either the sensitivity or the angular resolution. Table 1 provides estimates for the number of SBHMs with shadow sizes and optically thin horizon-scale flux densities that live on the ridge-line approximated by Equation 24; Table 2 lists the same for the number of sources we could expect to detect using telescopes with different resolution and sensitivity thresholds.

The case of M87 and the EHT
As of the writing of this paper, the SMBH in M87 is the only one whose shadow size (∼40 µas) and horizonscale flux density (∼0.5 Jy at 230 GHz) have been di-6 A "Pareto front" is the set of locations within a space of interest that satisfy the property that no one condition can be relaxed without making another more stringent. In our case, the "Pareto front" constitutes the set of locations in (θr, σν ) space where neither the angular resolution threshold nor the flux density threshold can be increased (i.e., made less demanding) without requiring a decrease in the other, while still being sensitive to the same number of objects. rectly imaged 7 (Event Horizon Telescope Collaboration et al. 2019a,b,c,d,e,f). M87 thus presents a natural test case against which to compare our source counts predictions from Section 2. Our model predicts that the number of SMBHs having ϑ > 40 µas and S ν > 0.5 Jy should be between ∼0.03 and ∼0.23 for the lower and upper BHMF prescriptions, respectively. Compared against the 1 object known to adhere to the chosen criteria, our model is systematically underpredicting the prevalence of M87. This underprediction may be explained at least in part if the local density of galaxies exceeds the cosmic mean, as suggested by, e.g., Dálya et al. (2018) 8 , which violates our model assumption of a homogeneous distribution of SMBHs. However, any such overdensity likely does not explain a discrepancy larger than a factor of ∼2, indicating that we may simply be finding ourselves on the high end of sampling variance. We thus expect that using the existence and properties of M87 to extrapolate the number of SMBHs with smaller shadows or weaker flux densities will result in systematically over-optimistic predictions; i.e., more sources will be predicted than our modeling suggests the real Universe likely contains. Similarly, the EHT is currently the only telescope to have successfully carried out shadow-resolved observations of a SMBH. The number of sources that the EHT is able to resolve and detect the shadows for thus presents a test case against which to compare our source counts predictions from Section 3. The EHT currently relies on observing with ALMA as part of the array, and during the 2017 observing campaign that led to the published M87 black hole images, ALMA itself required in-beam sources with flux densities of ≥0.5 Jy to perform the array phasing necessary for it to participate in VLBI observations (Matthews et al. 2018). For the purposes of estimating source counts, this phasing threshold effectively 7 The second shadow-resolved black hole that the EHT has targeted -the Milky Way SMBH Sgr A* -does not present a relevant comparison for this work because it is located in our own Galaxy, and it therefore does not fit within our modeling framework. In addition, Sgr A* has an additional observing constraint beyond those given in Section 2: it is heavily scattered by the ionized interstellar medium along its line of sight, so high-resolution observations must be conducted at correspondingly high frequencies of ν > ∼ 1 THz/ θr/1 µas (e.g., Lo et al. 1998;Bower et al. 2006;Johnson et al. 2018). The scattering is significantly weaker for sources off the Galactic plane (such as M87), requiring only ν > ∼ 30 GHz/ θr/1 µas (e.g., Cordes & Lazio 2002;Johnson & Gwinn 2015). Thus, interstellar scattering is unlikely to significantly affect our estimates for observable source counts. 8 We note that the overdensity in Dálya et al. (2018) is driven almost entirely by the existence of the Virgo cluster, and there are other indications (e.g., Tully et al. 2019;Böhringer et al. 2020) that when considering a somewhat larger volume (out to ∼100 Mpc) the local Universe may actually be underdense.
sets the sensitivity limit of the EHT. In this case, our model predicts that for θ r = 20 µas and σ ν = 0.5 Jy we should expect to resolve and detect up to ∼0.4 sources, similar to the projected number based on the above extrapolation using M87 as a benchmark. However, the 0.5 Jy phasing threshold has since been relaxed by permitting the transfer of phase corrections to faint targets from nearby but bright out-of-beam calibrators, and even the on-source phasing threshold can potentially be lowered through refinement of the phasing algorithm. Moving forward, the EHT may thus be able to observe much fainter targets. In a best-case scenario in which the phasing threshold is reduced to mJy levels, these improvements could permit the nominal sensitivity of the EHT to be used for source count estimates. Observing at 230 GHz, the EHT achieves θ r ≈ 20 µas and σ ν ≈ 10 −3 Jy, for which our model predicts the number of accessible SMBHs to be between ∼0.6 and ∼5.7 for the lower and upper BHMF prescriptions, respectively. We thus predict that the EHT could potentially gain access to approximately an order of magnitude more shadow-resolved sources by improving its effective sensitivity to mJy levels in this way.

Implications for array design
More generally, the behavior of N (θ r , σ ν ) -in particular, the behavior of its gradient -has implications for how an existing array can be most efficiently augmented to increase the number of accessible black hole shadows. As mentioned in Section 4.1, the EHT is currently operating with an angular resolution of θ r ≈ 20 µas and an effective flux density sensitivity between σ ν ≈ 10 −3 Jy and σ ν ≈ 0.5 Jy. This sensitivity range straddles the Pareto front for θ r ≈ 20 µas (see Figure 11), such that with σ ν ≈ 0.5 Jy the EHT array could most significantly increase the number of horizon-resolved black hole targets through improvements in sensitivity. However, once the sensitivity improves beyond the Pareto threshold of ∼70 mJy then the EHT will require enhanced angular resolution to increase the source counts further. For instance, at a fixed sensitivity of σ ν = 10 −3 Jy, an order-of-magnitude improvement in the angular resolution would yield an increase of roughly two orders of magnitude in the number of detectable black hole shadows; in contrast, while keeping the angular resolution fixed at 20 µas, arbitrary improvements in sensitivity beyond 10 −3 Jy would not yield many additional sources.
In practice, an Earth-based array like the EHT is limited to a maximum physical baseline length of one Earth diameter, meaning that any significant angular resolution improvements must come from increasing the observing frequency. A near-future aspiration for the EHT (Event Horizon Telescope Collaboration et al. 2019b), and a defining capability for the next-generation EHT (ngEHT; Doeleman et al. 2019;Raymond et al. 2021) will be to observe at a frequency of 345 GHz. At a fixed long-baseline sensitivity of σ ν = 10 −3 Jy, we expect that the effective 50% improvement in angular resolution over the current EHT should correspond to a factor of ∼3 increase in the number of detectable black hole shadows. In contrast, at a fixed θ r = 20 µas, the doubling of the baseline sensitivity that the ngEHT is expected to provide will only increase the source counts by ∼10%.
While angular resolution may ultimately limit the number of observable black hole shadows for groundbased interferometers like the EHT and ngEHT, sensitivity is expected to be the limiting factor for many prospective interferometers that network with spacebased stations. For instance, a baseline connecting a station on Earth to one located at the Earth-Sun L2 Lagrange point -such as may be possible using the proposed Millimetron (Kardashev et al. 2014) or Origins (Wiedner et al. 2021) space telescopes -would have a finest 230 GHz angular resolution of θ r ≈ 0.2 µas. At this resolution, we expect that a sensitivity of σ ν 10 −4 Jy would be required to detect even a single object. To achieve this sensitivity level, a 10-meter dish observing at 230 GHz as part of a baseline with the phased ALMA array would require a time-bandwidth product of ∼3 × 10 12 (e.g., 3 minutes of on-source integration time using 16 GHz of bandwidth), which is already larger than achieved by the EHT. Improving the sensitivity to 10 −5 Jy would require a time-bandwidth product that is two orders of magnitude larger still (e.g., 2 hours of on-source integration time using 32 GHz of bandwidth), and pushing to 10 −6 Jy would require an additional two orders beyond that (e.g., 5 days of onsource integration time using 64 GHz of bandwidth). Achieving the ∼10 −6 Jy Pareto front flux density corresponding to a ∼0.2 µas angular resolution thus imposes demanding sensitivity and stability requirements, and we expect that the number of sources accessible using long ( 1 Earth diameter) Earth-space baselines will be sensitivity-limited rather than resolution-limited.

SUMMARY AND CONCLUSIONS
Motivated by the success of the EHT and the promise of next-generation radio interferometric facilities, we have presented a framework for estimating the number of black hole shadows that are expected to be observationally accessible to different telescopes. Given assumptions about the distribution of SMBHs across mass, accretion rate, and redshift, we use a semi-analytic ADAF-based SED model to derive estimates for the number of SMBHs with detectable and optically thin horizon-scale emission as a function of angular resolution, flux density sensitivity, and observing frequency. Using a simple analytic prescription for the interferometric flux density distribution expected from black hole photon rings, we further decompose the SMBH source count estimates into the number of objects for which we could expect to observe first-and second-order photon rings.
Our main findings can be organized into two categories. First, we provide the following characterizations of the SMBH population: • Figure 4 shows the distribution of observationally accessible SMBH shadows, predicting that large numbers (>10 6 with ∼0.1 µas resolution and ∼1 µJy sensitivity) of objects should have resolvable horizon-scale emission at (sub)millimeter wavelengths.
• Figure 7 shows the angular resolution and sensitivity that an interferometer would require to observe the black hole shadows for this same population of SMBHs.
• For any particular choice of angular resolution and sensitivity, the population density of SMBHs with observable shadows generally increases toward higher redshifts and toward smaller black hole masses (see Figure 8). As a consequence, a majority of observable shadows are expected to have angular sizes that fall close to the resolution limit.
• The bulk population of SMBHs with observable n = 1 photon rings starts to become accessible at angular resolutions of 2 µas and flux density sensitivities of 0.5 mJy (see Figure 9 and Table 2). Similarly, the n = 2 population is accessible for angular resolutions of 0.1 µas and flux density sensitivities of 5 µJy.
We also consider the implications of these findings for current and future interferometric facilities: • The current effective sensitivity of the EHT is insufficient to maximally utilize its angular resolution. We predict that as many as ∼5 additional horizon-resolved sources could become accessible by improving the effective sensitivity of the EHT from ∼0.5 Jy to <70 mJy. ALMA should be sufficiently sensitive to achieve phased observations on sources with flux densities at this level, so an important next step will be to identify the specific sources that then become accessible. • Once the effective sensitivity of the EHT improves beyond the ∼tens of mJy level, a large (i.e., order-of-magnitude) additional increase in the number of observable black hole shadows can only be achieved by improving the angular resolution. We predict that an ngEHT observing at 345 GHz should have access to ∼3 times as many sources as the EHT observing at 230 GHz.
• Future telescopes that observe with 1 µas angular resolution, such as could be achieved using Earth-space interferometry, will require flux density sensitivities of 1 mJy to detect large numbers of black hole shadows.
In carrying out our analyses we have produced a library of synthetic SEDs and several tables of source counts 9 , as well as the code used to generate each SED 10 . The source count tables provide the predicted number of black hole shadows, n = 1 photon rings, and n = 2 photon rings accessible using different combinations of angular resolution, flux density sensitivity, and frequency. These resources may be useful for determining the specifications of future telescopes that aim to observe a large population of SMBH shadows or higher-order photon rings. Once such observations have been carried out, the predictive framework developed in this paper could be inverted so that the source counts become inputs rather than outputs, in turn providing constraints on the distribution of SMBH masses and accretion rates across cosmic history.

ACKNOWLEDGMENTS
We thank Gary Melnick for motivating conversations that sparked initial interest in pursuing this project. We also thank Avery Broderick, Tim Davis, Jason Dexter, and the anonymous referee for constructive comments that improved the quality of the paper. Support for this work was provided by the NSF through grants AST-1952099, AST-1935980, AST-1828513, AST-1440254, AST-1816420, and OISE 1743747, and by the Gordon and Betty Moore Foundation through grant GBMF-5278. This work has been supported in part by the Black Hole Initiative at Harvard University, which is funded by grants from the John Templeton Foundation and the Gordon and Betty Moore Foundation to Harvard University. NN acknowledges funding from Nucleo Milenio TITANs (NCN19−058). Yuan, F., & Narayan, R. 2014, ARA&A, 52, 529, doi: 10.1146 APPENDIX A. SED MODEL We use a spectral energy distribution (SED) model for ADAFs that largely follows the formalism presented in M97, though we introduce a number of modifications that update the SED to align it better with more recent work. In this section we detail these modifications, propagate them into the relevant expressions from M97 and Narayan & Yi (1995a, hereafter NY95), and describe the resulting SED model.
In determining the form of the SED, the primary equation we aim to solve is one of energy balance between the heating and cooling of the electrons in the flow. Following M97 Eq. (8), we have where Q + is the total viscous heating rate, δ is the fraction of this heating rate that goes directly to the electrons, Q ie is the rate of energy transfer from the ions to the electrons, Q − is the total radiative cooling rate of the electrons, and we have introduced an additional term Q adv,e that accounts for the electron energy that gets advected into the black hole. We note that in the extremely low accretion regime considered here, energy loss from neutrino cooling is negligible. The radiative cooling term is given by: where P synch , P compt , and P brems correspond to the power emitted in synchtrotron, inverse Compton, and bremsstrahlung radiation, respectively. It is the combined contributions from these three emission processes that ultimately constitute our model SED.
The emission processes of interest for this paper depend on the electron temperature, which is determined selfconsistently such that the total heating and cooling satisfy Equation A.1. The left panel of Figure 13 shows the various model contributions to the electron heating and cooling as a function of electron temperature for an example M87-like system, and the right panel shows the corresponding predicted SED as a function of frequency. The right panel of Figure 14 shows the derived temperatures as a function of m andṁ. Table 3 provides a list of the various parameters used in the SED model, and Figure 12 shows example SEDs.

A.1. Flow equations
We take the underlying accretion flow properties to be described by the self-similar models developed by NY95, in which the relevant parameters are the black hole mass M , the accretion rateṀ , the radius R, the viscosity parameter α, the ratio of gas to magnetic pressure β, 11 and the fraction f of viscously dissipated energy that gets advected into the black hole. Following M97, we use scaled quantities, where R S = 2GM/c 2 is the Schwarzschild radius,Ṁ Edd = L Edd /ηc 2 is the Eddington accretion rate, and we have taken the radiative efficiency η to be 0.1. Here, the difference between Equation A.3c and M97 Eq. (4) comes from our adoption of the radius-dependent accretion rate from Blandford & Begelman (1999), which accounts for outflowing material via a radial dependence of the mass accretion rate with power-law index s.
The self-similar equations describing the accretion flow, M97 Eq. (5), become ρ =Ṁ 4πHαc 1 √ GM R = 6.022 × 10 −5 g cm −3 α −1 c −1 1 m −1ṁ 0 r −(3/2)+s , (A.4a) where ρ is the mass density, B is the magnetic field strength, n e is the number density of electrons, α is the disk viscosity parameter (Shakura & Sunyaev 1973), H is the disk scale height (we have followed M97 in setting H = R), µ e = 1.14 is the mean molecular weight (NY95), and c 1 ≈ 0.5 and c 3 ≈ 0.3 are constants defined in NY95 and specified in Table 3. We adopt a power-law radial profile for the electron temperature T e of the form T e = T e,0 r 1−t , (A.5) with t ≤ 1. From NY95 Eq. (2.16), the two-temperature accretion flow must satisfy T i + 1.08T e = 6.66 × 10 12 K (1 + β) −1 β c 3 r −1 , where T i is the ion temperature. Setting T i = T e at some maximum radius r = r max yields an expression for t, t = 1 ln (r max ) ln 6.66 × 10 12 K βc 3 2.08T e,0 (1 + β) , such that T i > T e for all r < r max .

A.2. Heating
The plasma in an ADAF is heated by viscous forces, with the total heating rate per unit volume denoted as q + . Some fraction δ of this energy gets deposited into the electrons, while the remaining fraction (1 − δ) heats the ions. The ions can transfer thermal energy to the electrons via Coulomb collisions, with the rate of this transfer denoted by q ie , and the electrons can radiate energy away at a rate q − . Taken altogether, energy balance yields advected energy rates of for the ions (q adv,i ) and electrons (q adv,e ). The ion heating is driven by viscous dissipation, while the dominant electron heating source depends on the accretion rate; at high accretion rates the ion-electron heating is dominant, whereas at low accretion rates the viscous heating is more important. NY95 give an expression for the viscous heating rate per unit volume, where f is the fraction of viscously dissipated energy that gets advected into the black hole. The total (i.e., volumeintegrated) viscous heating rate is then given by where r min and r max are the minimum and maximum radius, respectively. 12 The heating rate per unit volume of the electrons from Coulomb interactions with protons is given by Stepney & Guilbert (1983), where θ e = kT e /m e c 2 is the dimensionless electron temperature, θ i = kT i /m p c 2 is the dimensionless ion temperature, ln(Λ) ≈ 20 is a Coulomb logarithm, K n represents a modified Bessel function of the nth order, and we have assumed n e = n i . In the second line we have adopted the approximation from M97. 13 We note that in evaluating the prefactor in Equation A.11 we have followed NY95 and multiplied by an additional factor of 1.25 to account for the ions containing a mixture of roughly 75% hydrogen and 25% helium. The volume-integrated ion-electron heating rate is given by Q ie = 3.236 × 10 17 cm 3 m 3 rmax rmin q ie r 2 dr, (A.12) which does not have an analytic form and so must be integrated numerically.

A.3. Cooling
The observed emission in radio and (sub)millimeter bands is dominated by synchrotron radiation, but the primary electron cooling mechanisms also include bremsstrahlung and inverse Compton radiation. Each of these emission mechanisms contributes to q − , and each depends on the electron temperature T e .

A.3.1. Synchrotron emission
We use a form for the synchrotron spectrum from (Mahadevan et al. 1996, see also NY95; M97), which assumes an isotropic distribution of relativistic electrons. The synchrotron spectral emissivity is given by synch,ν = 4.43 × 10 −30 erg s −1 Hz −2 4πn e ν K 2 (1/θ e ) M (x M ), (A.13) 12 We note and correct an error in the original expression for Q + from M97 Eq. (9), for which the exponent of the c 3 term should be 1/2 rather than 1. 13 We note and correct an error in the original expression for q ie from M97 Eq. (10), for which the exponent of the r term should be −3 rather than −1. and x M is a dimensionless frequency, .13 assumes optically thin emission, but below some critical frequency ν c (which is a function of radius) we expect the synchrotron to be optically thick and thus described by a blackbody spectrum. We follow M97 and determine ν c (r) by equating emission within a volume of radius r to the Rayleigh-Jeans blackbody emission from a spherical surface at that radius, We numerically solve Equation A.16 for ν c at R min and R max , yielding a peak frequency ν p at R min (with luminosity L p ) and a minimum frequency ν m at R max (with luminosity L m ); the left panel of Figure 14 shows how ν p changes with m andṁ. We take the synchrotron spectrum to be blackbody (i.e., optically thick) at frequencies below ν m , optically thin with an emissivity described by Equation A.13 at frequencies above ν p , and a power law at intermediate frequencies.
That is, The total emitted synchrotron power is then the integral of L ν,synch over frequency, which we evaluate numerically.

A.3.2. Bremsstrahlung emission
We use an expression for the bremsstrahlung emission that follows M97 Eq. We integrate both of the above expressions numerically.

A.3.3. Inverse Compton emission
We follow M97 in considering Comptonization only of synchrotron photons emitted predominantly at the peak frequency ν p , for which the spectrum in the temperature range of interest is expected to be a power law, The power-law index α c is determined by both how frequently photons are scattered (which is determined by the optical depth of the scattering process) and by how much a photon gets amplified during a typical scattering event.
We use an expression for the optical depth to electron scattering τ es adapted from NY95 Eq. (2.15), We take the mean amplification factor A from M97 Eq. (32) (originally inspired from Rybicki & Lightman 1979), (A.26) which together with τ es determines the power-law slope for the Compton emission, The total Compton power will then be given by the integral of L ν,compt up to the maximum final frequency of a Comptonized photon (ν f = 3kT e,0 /h),

. Electron advection
The M97 model assumed that T e T i , and it therefore ignored electron energy advection. This assumption was reasonable for the parameters considered in that paper, particularly the choice of δ = m e /m p . However, the modern view is that δ is much larger (≈ 0.3; see Yuan & Narayan 2014). Such large values of δ make electrons significantly hotter, especially at very lowṁ, and so energy advection in electrons can no longer be ignored.
When electron advection is included, M97 Eq. (8) gains an additonal term Q adv,e and becomes Equation A.1. This advective cooling term is given by where, following the approach described in Narayan & Yi (1994), we express the specific heat at constant volume C V in terms of an effective γ CV . Substituting in Equation A.4b and Equation A.5 and differentiating with respect to R yields Sadowski et al. (2017) provide an accurate fitting function for γ CV , which we write as γ CV = 20 2 + 8θ e + 5θ 2 e 3 (8 + 40θ e + 25θ 2 e ) . (A.32) Noting further that n e = ρ/µ e m p andṀ = −4πR 2 vρ, we finally obtain Q adv,e = 1.013 × 10 26 erg s −1 K −1 mṁ 0 T e,0 rmax rmin 1 − t γ CV − 1 − 3 2 + s r s+t−2 dr, (A.33) which we integrate numerically. We note that there are conditions under which Equation A.33 can yield a negative value for Q adv,e ; in these cases, we impose Q adv,e = 0.

A.4. Maximum mass accretion rate
The ADAF solution ceases to exist above some critical mass accretion rate,ṁ crit , where the accretion flow is no longer advection-dominated (NY95; M97). Within the context of our SED model, this condition manifests as a maximum accretion rate above which there is no equilibrium temperature (i.e., the heating and cooling curves never cross). We numerically determine a value ofṁ crit ≈ 10 −1.7 , and so in this paper we only work with values ofṁ 0 ≤ 10 −2 .
In addition to the criticalṁ above which no ADAF solutions exist, there is also a softer threshold accretion rate above which solutions do exist but our assumed input radiative efficiency of η = 0.1 is no longer consistent with the output of the SED model. Figure 15 shows the predicted radiative efficiency from the model as a function ofṁ 0 ; we take the model radiative efficiency to be the ratio of the bolometric luminosity, L bol = ∞ 0 (L ν,synch + L ν,compt + L ν,brems ) dν, (A.34) to the accretion rate equivalent luminosity,Ṁ c 2 . Regardless of the input mass, the output radiative efficiency exceeds the assumed input value forṁ 0 10 −2.5 . Though this inconsistency reflects a physical limitation of the model, we note that given the ERDF prescription used in this paper (see Section 2.4) it impacts only a small fraction of SMBHs (<1% for most M and z, reaching a peak of ∼5% for M > 10 9 M and z > 5). 1.4 1.42 × 10 9 c 3 αc 1 (1+β) 1.5 × 10 9 s2 1.19 × 10 −13 xM . . . Note-A list of the parameters used for the SED model. Certain parameters in the model take on the default values listed here, while others must either be specified as inputs (e.g., m,ṁ) or else are internally computed as part of the model (e.g., Te, t). ) for example SMBH masses of 10 6 and 10 9 M , respectively, as labeled in the legend at the lower left. All panels share the same horizontal and vertical axis ranges, which are explicitly labeled in the left panel.

B. MASS DEPENDENCE OF THE EDDINGTON RATIO DISTRIBUTION FUNCTION
As described in Section 2.4, in this paper we take the ERDF to have a broken power-law functional form (see Equation 12) with a power-law index α that evolves with both SMBH mass M and redshift z. Specifically, we adopt the redshift evolution prescription from Tucci & Volonteri (2017) (see Equation 14), and we add to it an evolution with SMBH mass (see Equation 15). For the mass evolution of the ERDF power-law index, we choose a logistic function in log(M ) such that low-mass SMBHs (i.e., those with masses below some value M 0 ) see a power-law index a lo while high-mass SMBHs (i.e., those with masses above M 0 ) see a power-law index a hi . The specific functional form of Equation 15 ensures that a(M ) transitions smoothly between the low-and high-mass regimes, with a logarithmic width that is set by the parameter ∆.
To determine the values of the ERDF parameters a lo , a hi , M 0 , and ∆, we rely on the observational constraints provided by Aird et al. (2018). Aird et al. (2018) determined the distribution of specific SMBH accretion rates λ si.e., the accretion rate relative to the stellar mass of the galaxy, rather than to the mass of the SMBH -by fitting a Bayesian mixture model to X-ray observations of ∼10 5 near-infrared-selected galaxies. This sample includes a mix of star-forming, quiescent, and AGN-dominated galaxies, and it spans a range ∼10 8.5 -10 11.5 M in stellar mass and ∼0.3-4 in redshift. We use two different prescriptions to convert from stellar mass to SMBH mass, corresponding to the two BHMF prescriptions described in Section 2.2. For our fiducial choice of the lower BHMF from Shankar et al. (2009), we adopt the stellar-to-SMBH conversion used by Aird et al. (2018) themselves, which is given simply by M * = 500M ; we use this fiducial prescription for all figures and values in this paper unless otherwise specified. For the instances in which we quote a range of values corresponding to the lower and upper BHMFs, for the upper BHMF we convert from λ s to λ using the same stellar-to-SMBH conversion as in Section 2.2 (i.e., Equation 11).
Given the observed P (λ) as a function of M and z, we determine the best-fit ERDF parameters by minimizing the squared logarithmic differences between the Aird et al. (2018) empirical model and Equation 12. We restrict our fitting to the region of parameter space between −4 ≤ log(λ) ≤ −1, with the low-λ cutoff determined by the observational limitations and the high-λ cutoff determined by our interest in LLAGNs. The resulting best-fit values for the fiducial case are a lo = 0.55, a hi = 0.20, log(M 0 ) = 7.5, and ∆ = 0.3, and Figure 16 shows a comparison between the best-fit ERDF and the Aird et al. (2018) model. For the upper BHMF prescription, the only parameter that changes is M 0 , for which we find a best-fit value of log(M 0 ) = 7.8.  Figure 4, but with the source counts predicted by the best-fitting analytic approximation (Equation C.2) overplotted in cyan contours. Right: Same as Figure 7, but with the source counts predicted by the best-fitting analytic approximation (Equation C.5) overplotted in cyan contours.
from Section C.1 while also accounting for this brightness temperature limit, we can modify Equation C.1 using an exponential cutoff to smoothly suppress the the high brightness temperature emission, N (θ r , σ ν ) ≈ e −T b /(10 10 K) θ r 40 µas 3 + σ ν 0.1 Jy 3/2 −1 . (C.4) Here, θ r should now be understood to represent a single-baseline angular resolution, and we have adjusted the σ ν normalization to match the 230 GHz flux density observed from M87 on long baselines (Event Horizon Telescope Collaboration et al. 2019c). We have set the brightness temperature cutoff to 10 10 K because we are selecting for SMBHs that are optically thin and which therefore should not typically saturate the brightness limit. We note that magnetohydrodynamic simulations of the M87 system also exhibit brightness temperatures that peak between 10 10 K and 10 11 K (Event Horizon Telescope Collaboration et al. 2019e). Following these expectations, we expand on the results of Section C.1 and fit a simple functional form to the source counts of Here, we have fixed the exponents of the θ r and σ ν terms to the values expected from initial considerations and further motivated by the fitting results of Section C.1. We also incorporate the understanding from Section 3.2 into the scale-setting parameters θ r,n and σ ν,n for the nth sub-ring, which are defined to be θ r,n = θ r,0 e −nπ σ ν,n = σ ν,0 e −3nπ/2 . (C.6) We do not have an a priori expectation for the scaling behavior of the brightness temperature T b,n in each sub-ring, so we simply include it as an additional parameter, The model has five free parameters -θ r,0 , σ ν,0 , T b,0 , µ, and C -which we fit in the same manner described in Section C.1. We find best-fit parameter values of θ r,0 = 23.2 µas, σ ν,0 = 0.17 Jy, T b,0 = 2.2 × 10 8 K, µ = 0.50, and C = 2.39; this best fit is shown in the right panel of Figure 17. The fit quality is similar to that in Section C.1, with the analytic source counts throughout most of the (θ r , σ ν ) space agreeing to better than an order of magnitude with the numerical results. The deviations become worse than an order of magnitude at small values of θ r 0.1 µas and σ ν 10 −7 Jy, as well as wherever the source counts contain substantial contributions from n > 0 photon rings.  Figure 18. Left: Schematic diagram of black holes sampled by a single baseline at one instant. The dashed black line represents the baseline, the solid gray line represents the line of sight to an example black hole, the dotted black line represents the projection of the baseline perpendicular to that line of sight, and the dotted blue curve represents the sampling function (with radial distance from the baseline center proportional to the number of black holes sampled). Middle: Same as left, but after the baseline evolves, rotating by an angle ∆χ, sweeping through the volume enclosed by the dashed red line. Dotted grey lines show the decomposition performed in Equation D.7. Both left and middle diagrams correspond to a two-dimensional cross-section of a three-dimensional surface; for a black hole distribution that follows d 3 N dΩdθr ∝ θ −α r with α = 2, the instantaneous sampling surface is exactly a torus. Right: Cumulative fraction of black holes sampled by sky density functions with increasingly steep dependence on resolution.

D. SAMPLING FUNCTION FOR A ROTATING BASELINE
Following Section 3.3, we would like to determine how the total number of sources that a particular baseline can resolve changes as that baseline rotates. Equation 22 provides the angular resolution θ r that a baseline of length b (in units of the observing wavelength) has when viewed from a source at a sky position with polar angle θ as measured from the tip of the baseline (see the left panel of Figure 18); there is no φ dependence because sources are assumed to be distributed isotropically. For any particular baseline, the number density of resolvable sources per unit solid angle can be written as which can then be used as in Equation 23 to determine the total number of resolvable sources for any fixed baseline b.
For the example calculations presented in this appendix, we will assume that the integrand of Equation D.1 follows a power law with index α and coefficient A, 2) The number of black holes instantaneously sampled is then given by an integral over solid angle, expressed in Equation 23 but now evaluated explicitly for the power law: where Γ is the gamma function. In order to compare the cumulative sampling of different power laws, we must normalize so that after the baseline rotates by 180 • , the number of black holes sampled is equal; this rotation corresponds to sweeping the largest projected spacing across the entire sky. Because the largest projected spacing is the only relevant quantity, we can normalize simply by requiring for all A, b, and α. We express the number of black holes sampled as a fraction of the total, removing the dependence on A and b: It is then straightforward to observe that the fraction of black holes in the sky instantaneously observed by a single baseline decreases with α. However, it is less straightforward to compute the cumulative number of unique sources a baseline resolves over the course of some rotation through space, because the sources sampled at each baseline orientation are partially redundant with the source sampled at prior baseline orientations. Put another way, the instantaneous sampling of black holes provided by a single baseline is given by a sweeping of the projected baseline through φ, whereas the sampling over a change in orientation is a sweep through a new angle, which we call χ. The angle χ is related to θ by χ = tan −1 (sin φ tan θ) .

(D.6)
χ is defined so that it aligns with θ when φ = π/2. An illustration of χ and the redundant sampling are shown (for the φ = π/2 cross-section) in the middle panel of Figure 18. The rotation of the baseline in space can then be described by a change in the axial angle ∆χ. Computing Equation 23 for the total sources sampled after this rotation (shown by the dashed red surface in the middle panel of Figure 18) can be simplified by breaking the integral into two regions: first, the longest projected baseline sweeps out a partial spheroid, while the rest of the baselines form arcs that intersect at the cusps of the dashed red line in Figure 18. We refer to the sources sampled by the partial spheroid as N s and those sampled by the cusped curves as N 4 . As ∆χ increases, N s increases and N 4 decreases. The cusped curves have a four-fold symmetry, so we integrate over a convenient curve (that between θ = π/2 and θ = π), and Equation 23  is the value of θ at the first leading cusp. The geometry of this computation is shown between the dotted gray lines in the middle panel of Figure 18. The right panel of Figure 18 shows the cumulative source sampling for several values of the power-law index α from Equation D.5. We use α = 2 for the schematic diagrams in Figure 18 because the sampled surface around an instantaneous baseline in this case reduces to a torus; α = 4 corresponds to a sky density of black holes that scales volumetrically (i.e., proportional to θ −3 r ), which is closer to the actual behavior. For plausible values of α, we find that the difference between a fully swept sampling of the sky and the instantaneous baseline sampling is not more than a factor of 2.