Exploring Proxies for the Supermassive Black Hole Mass Function: Implications for Pulsar Timing Arrays

Supermassive black holes (SMBHs) reside at the center of every massive galaxy in the local Universe with masses that closely correlate with observations of their host galaxy implying a connected evolutionary history. The population of binary SMBHs, which form following galaxy mergers, is expected to produce a gravitational wave background (GWB) detectable by pulsar timing arrays (PTAs). PTAs are starting to see hints of what may be a GWB, and the amplitude of the emerging signal is towards the higher end of model predictions. Simulated populations of binary SMBHs can be constructed from observations of galaxies and are used to make predictions about the nature of the GWB. The greatest source of uncertainty in these observation-based models comes from the inference of the SMBH mass function, which is derived from observed host galaxy properties. In this paper, I undertake a new approach for inferring the SMBH mass function starting from a velocity dispersion function rather than a galaxy stellar mass function. I argue that this method allows for a more direct inference by relying on a larger suite of individual galaxy observations as well as relying on a more"fundamental"SMBH mass relation. I find that the resulting binary SMBH population contains more massive systems at higher redshifts than previous models. Additionally, I explore the implications for the detection of individually resolvable sources in PTA data.


INTRODUCTION
Supermassive black holes (SMBHs), with masses of 10 6 -10 10 M , reside in the nuclei of most, if not all, nearby galaxies. Shared evolution likely gives rise to the close correlation between the mass of the SMBH and several global properties of the host galaxy (Kormendy & Richstone 1995). Following a galaxy merger, the SMBHs from each galaxy sink to the center of the common merger remnant through interactions with the galactic gas, stars and dark matter (Begelman et al. 1980). Once the SMBHs become gravitationally bound, they will emit strong gravitational waves (GWs) in the nanohertz frequency band.
Utilizing a hierarchical framework of galaxy formation, many theoretical models of SMBH evolution have been used to infer the binary SMBH population and the resulting gravitational wave background (GWB) (e.g., Rajagopal & Romani 1995;Jaffe & Backer 2003;Wyithe & Loeb 2003;Sesana et al. 2008;Chen et al. 2019). Full cosmological simulations have also been used to great effect from the Millennium simulation (Sesana et al. 2009;Ravi et al. 2012) to the Illustris simulation (Kel- * NSF Astronomy and Astrophysics Postdoctoral Fellow ley et al. 2017; Siwek et al. 2020). In recent years, focus has turned toward modeling the GWB based on the plethora of well-constrained galaxy observations (Sesana 2013;McWilliams et al. 2014;Ravi et al. 2015;Mingarelli et al. 2017). Simon & Burke-Spolaor (2016, hereafter SB16) undertook an in-depth analysis of how the uncertainty in galaxy observations propagated into the resulting GWB prediction and found that the inference of the SMBH mass function from host galaxy parameters had the largest impact.
Over the past few decades, precision timing of millisecond pulsars has allowed the creation of a galactic-scale GW observatory, a pulsar timing array (PTA), which is sensitive to GWs at nanohertz frequencies (Alam et al. 2021;Perera et al. 2019). The brightest signature detectable by PTAs is the GWB, which results from the incoherent superposition of GWs from SMBH binaries that form in post-merger galaxies. While the first stages of galaxy mergers have been repeatedly identified through electromagnetic observations (e.g., Duncan et al. 2019;Stemo et al. 2020), bound SMBH binaries, at sub-parsec separations, remain elusive. In fact, PTAs offer one of the most direct avenues to observing bound SMBH binaries (Burke-Spolaor et al. 2019).
Recently, new PTA datasets are uncovering a signal which may be the first hints of the GWB (Arzoumanian et al. 2020;Goncharov et al. 2021;Chen et al. 2021;Antoniadis et al. 2022). If the signal is of an astrophysical origin, then confirmation of the signal's GW nature is expected in the next year or two (Pol et al. 2021). The emerging signal is towards the higher end of predicted amplitudes (Rosado et al. 2015;Kelley et al. 2017;Chen et al. 2019), which may indicate that binary SMBHs are more massive than previously thought (Middleton et al. 2021) and/or that the local number density of SMBHs is larger than expected (Casey-Clyde et al. 2021). To help understand these implications, it is worth revisiting the inference of the SMBH mass function from galaxy observations.
Inference of the SMBH mass function is predicated on measurements of host galaxy relationships in the local Universe (e.g., M • -M bulge and M • -σ). Previous work has used observations of the galaxy stellar mass function (GSMF) as the basis for the host galaxy relation used to infer the SMBH mass function (Sesana 2013, SB16), however this method requires some assumption to be made about the fraction of the total stellar mass of the galaxy contained in the bulge, f bulge . While the majority of contributions to the GWB is expected to be from bulge-dominated, early-type galaxies where f bulge ≈ 1 (Sesana 2013), it is increasingly unclear how galaxy morphology tracks SMBH mass at higher redshift (Chen et al. 2020), where galaxies tend to be smaller and clumpier (Conselice 2014).
A more direct inference could be made using a galaxy's stellar velocity dispersion (σ), which is measured from a galaxy's spectra. In addition to removing the need for an f bulge assumption, recent work suggests that M • -σ is a more "fundamental" relationship than M • -M bulge (e.g., Kormendy & Ho 2013;van den Bosch 2016;de Nicola et al. 2019;Marsden et al. 2020) and therefore more accurately predicts the SMBH mass. In particular the discovery of relic galaxies, such as NGC1271 and NGC1277 (Walsh et al. 2015(Walsh et al. , 2016, which formed at z > 2 and have evolved passively ever since (Ferré-Mateu et al. 2015), points to the limitations of the M • -M bulge relation at higher redshifts. These relic galaxies are much smaller than most present day early-type galaxies, which is consistent with the strong redshift evolution of galaxy size (van der Wel et al. 2014), however they host extremely massive SMBHs. While these galaxies are significant outliers on the M • -M bulge relation (Ferré-Mateu et al. 2015), they are consistent within the intrinsic scatter for the M • -σ relation (van den Bosch 2016).
Until recently there were no spectroscopic surveys dedicated to measuring the galaxy velocity dispersion function (VDF) at higher redshifts, however, the Large Early Galaxy Astrophysics Census (LEGA-C) survey now provides the first direct spectroscopic measurement of the VDF from 0.6 < z < 1.0 (Taylor et al. 2022).
The LEGA-C result is consistent with earlier attempts to indirectly measure the VDF using an inferred velocity dispersion, calculated from photometric data (Bezanson et al. 2011(Bezanson et al. , 2012, especially at large σ. Previous attempts to use the M • -σ relation in GWB calculations have started directly from the GSMF and used a galaxy stellar mass to σ conversion derived from observations in the local Universe (Sesana 2013;Sesana et al. 2016), however, these results neglected to include the strong redshift evolution in the galaxy size-mass relation (van der Wel et al. 2014). Thus, by incorporating more observational data, the VDFs used in this work provide fundamentally different starting points for calculating the SMBH mass function than anything that has been attempted before.
In this paper, I will use the VDFs from both Bezanson et al. (2012) and Taylor et al. (2022) as the basis for inferring the SMBH mass function and compare the resulting binary SMBH population and GWB signal to that from the standard GSMF inference method. In the next section I describe my methods, including how to infer the SMBH mass function with these different approaches. I show my results in the following section, and close with a summary.

METHODS
Following SB16, this paper calculates the characteristic strain spectrum h 2 c (f ) from a cosmic population of binary SMBHs by integrating over all the sources emitting in an observed gravitational wave frequency bin df multiplied by the square strain of each source in that bin (e.g., Sesana et al. 2008), where d 4 N is explicitly the number of binaries in a given redshift range dz, primary black hole mass range dM • , and black hole mass ratio range dq • , which are emitting in a given Earth-observed GW frequency range df , and h s is the polarization-and sky-averaged strain contribution from each binary (e.g., Thorne 1987), where M c = M • (q 3 • /(1 + q • )) 1/5 is the chirp mass of the binary, D c is the proper (co-moving) distance to the binary, and f r is the frequency of the GWs emitted in the rest frame of the binary. The Earth-observed GW frequency f = f r /(1 + z).
The number of binaries d 4 N is directly related to the comoving number density per unit redshift, primary black hole mass and binary mass ratio, The conversion above takes the number of binaries per co-moving volume shell, dV c , and converts it to the number of binaries per Earth-observed GW frequency bin, df , by first converting to redshift, and then to Earthobserved time. Once the binary hardens and decouples from the surrounding galactic environment, the orbital changes of the system become dominated by the emission of gravitational radiation at a rate of (Peters & Mathews 1963), For binaries in a circular orbit, the frequency of GWs emitted in the rest frame of the binary f r = 2f orb . The frequency dependence of h c is captured in both h s and dt/df , therefore by combining Eqns. 1-4, h c (f ) can be written as a simple power-law with dimensionless amplitude A yr referenced to the frequency of an inverse year f yr = 1 yr −1 , (Jenet et al. 2006).
However, this only applies for circular binaries in GWradiation driven inspirals. If the inspirals were instead driven by some environmental coupling in the lowest frequencies of the PTA band, this would cause a turn-over in the power-law spectral shape (Sampson et al. 2015). Environmental coupling can also impact the eccentricity of the binary and therefore the spectral shape of the characteristic strain spectrum in the PTA band (Ravi et al. 2014;Huerta et al. 2015). For the purposes of this project, where I am solely focused on the SMBH mass function, I ignore the specifics of environmental coupling as well as eccentricity and assume that all binaries remain in circular orbits while emitting GWs in the PTA band.

Determining the Number Density of Binary SMBHs
The number density of binary SMBHs is a combination of two functions: the SMBH mass function Φ • (z, M • ) and the binary SMBH merger rate There are currently no direct observational constraints on the demographics of binary SMBHs, so these functions are not explicitly known. Instead, they need to be inferred from other data. As in SB16, I use the demographics of galaxy mergers as a proxy. The binary SMBH merger rate is directly related to the galaxy merger rate, R • ∼ R Gal , but offset in redshift due to the binary evolution timescale. The standard way to calculate the galaxy merger rate, R Gal , directly from observations is by combining a galaxy close pair fraction, f pair , with a merger timescale, τ . The galaxy close pair fraction is an astronomical observable, while the merger timescale, which approximates the amount of time galaxy close pairs are observable, must be determined through simulations (SB16).
The SMBH mass function, Φ • , is more difficult to derive, and needs to be indirectly inferred by populating each galaxy with a central SMBH assuming some relationship between SMBHs and their host galaxies (e.g., M • -σ relation, M • -M bulge relation, etc.). The standard approach to inferring the SMBH mass function from galaxy surveys is to start from an observed galaxy stellar mass function (GSMF), then make some assumption about the fraction of the total stellar mass of the galaxy contained in the bulge, f bulge , and finally calculate the SMBH mass from the M • -M bulge relation (Sesana 2013;Ravi et al. 2015, SB16). Using this approach, the number density of binary SMBHs is directly inferred from the number density of galaxy mergers, where (Gal → •) represents the step where the SMBH is populated into the galaxy assuming some galaxy number density parametrization Φ Gal . When using a GSMF, Φ Gal = Φ GSMF (z, M ), where M is the total stellar mass of the primary galaxy. The companion galaxy then has a mass of q M M , where 0.25 < q M < 1 is set to capture major merger events, and the (Gal → •) step in Equation 7 utilizes the M • -M bulge relation and includes the determination of f bulge for each galaxy. While the specific prescription for f bulge can change the predicted amplitude of the GWB by a factor of two or more, the largest source of uncertainty in this inference of the SMBH mass function inference comes from uncertainty in the GSMF itself (SB16). This is to be expected given discrepancies over the number density of the most massive galaxies in the local Universe (Bernardi et al. 2013). In addition to the issues with properly identifying massive galaxies, there is increasing evidence that SMBH evolution and bulge evolution track differently at higher redshifts (z > 1) (Chen et al. 2020), however, there is no significant evidence that these different evolutionary tracks are caused by a change in the relationship between SMBH mass and other host galaxy properties (Shen et al. 2015;Suh et al. 2020). Ideally, one would like a more robust method for inferring the SMBH mass function, particularly something that doesn't require an f bulge determination and doesn't have as much dependency on the underlying mass-to-light (M/L) model used to determine stellar mass (Bernardi et al. 2013).
The stellar velocity dispersion (σ) appears to allow for a more direct inference of SMBH mass. Recently, it has been suggested that σ is a more fundamental property of a galaxy than it's stellar mass (van den Bosch 2016; de Nicola et al. 2019; Marsden et al. 2020), and that the M • -M bulge relation is just an extension of the more fundamental M • -σ relation (Wake et al. 2012). However, it has been difficult to obtain VDFs, specifically at higher redshifts (z > 0.3), due to a lack of complete spectroscopic surveys to use as a separate starting point for inferring the SMBH mass function. Instead, when σ has been explored in the context of simulating the binary SMBH population, it has been inferred from galaxy total stellar mass, and therefore, the results have unsurprisingly been consistent (Sesana 2013;Sesana et al. 2016). Instead, in an effort to reduce uncertainty in the inferred SMBH mass function, we are pursuing a novel approach to modeling M • using the velocity dispersion functions from both Bezanson et al. (2012), which are derived from an inferred velocity dispersion (σ inf ), and Taylor et al. (2022), which are derived from spectroscopic observations at 0.6 < z < 1.0.

The Velocity Dispersion Function
When using a VDF in Equation 7, Φ Gal = Φ VDF (z, σ), where σ is the velocity dispersion of the primary galaxy. The companion galaxy then has a velocity dispersion of q σ σ, where q σ is set to capture major mergers the same way that q M was in the case of GSMFs. Finally, the (Gal → •) step in Equation 7 now utilizes the M • -σ relation and no longer requires any additional steps as it did in the case of GSMFs. Taylor et al. (2022) measures the VDF directly through spectroscopic observations. Bezanson et al. (2011) uses a suite of photometric data to calculate an inferred VDF, which is shown to be comparable with a spectroscopic measurements in the local Universe.
In Bezanson et al. (2011), the inferred velocity dispersion, σ inf , is measured from a combination of a galaxy's photometric properties. Starting from the virial theorem and building from the observational evidence that for local galaxies stellar mass is proportional to dynam-ical mass (Taylor et al. 2010), the central velocity dispersion, σ, of a galaxy can be effectively predicted based on its size (effective radius, r e ), shape (Sèrsic indec, n), and total stellar mass (M * ).
where (Bertin et al. 2002), and since Taylor et al. (2010) has shown that K * (n) only weakly depends on mass, the average ratio of stellar to total mass is adopted and calibrated such that the median σ inf equals the median measured σ in SDSS data. Bezanson et al. (2011) finds a scatter in the σ inf -σ relationship of 0.06 dex, which has the effect of increasing the high-σ end of the σ inf distribution. This scatter must be properly taken into account when inferring the "true" velocity dispersion function (VDF) in order to not over predict the number density of massive galaxies. Numerical fits for the "true" VDF were not published in either Bezanson et al. (2011) or Bezanson et al. (2012), and instead were obtained through private communication with the author (Bezanson 2016). The fits are done using a modified Schechter function and, given the strong correlation between σ * , α, and β (Sheth et al. 2003), uncertainties on the fit were obtained by setting α and β to their maximum likelihood values and allowing σ * to absorb the combined uncertainties in those three parameters.

Observational Constraints
Unlike previous work (e.g., SB16), this paper is focused on keeping as many things constant in the calculation of the GWB as possible, in order to highlight the differences in the SMBH mass function and the implications for PTA experiments. As such, this work only uses the galaxy close pair fraction from Keenan et al. (2014) in it's calculation of R Gal , as opposed to Robotham et al. (2014) since SB16 showed them to be consistent. For τ , the same formulation as SB16 is used, which combines the mass and redshift dependence from Kitzbichler & White (2008) with the results from Lotz et al. (2011). Bezanson et al. (2012) calculates the VDF at 0.3 < z < 1.5 using photometric measurements of galaxy stellar mass, effective radii, and Sèrsic index for the galaxies in the NMBS COSMOS (Whitaker et al. 2011) and the UKIDSS UDS fields (Williams et al. 2009). In an attempt to make sure that the galaxy samples are as consistent as possible, this paper compares the VDF derived from the above σ inf data to the GSMF from Tomczak et al. (2014) since it is partially based on NMBS COSMOS. It is important to also note that the stellar mass determination in Bezanson et al. (2012) and Tomczak et al. (2014) use the same software pipeline, which should reduce the amount of bias inserted into these calculations from different M/L determinations. When using the spectroscopic VDF from Taylor et al. (2022) for 0.6 < z < 1.0, the inferred VDF is used at z > 1 since the inferred VDF matches the spectroscopic VDF.
For inference of the SMBH mass function at z < 0.3, this paper uses the GSMF in Moustakas et al. (2013) calculated from SDSS data since that is what was used in SB16, and it was also used as the local Universe comparison in Tomczak et al. (2014). To minimize introducing additional bias, the VDF from Bernardi et al. (2010) is used at z < 0.3, since it was directly derived from similar SDSS data as Moustakas et al. (2013).
To infer the SMBH in each galaxy, the M • -M bulge relation is used with the GSMF while the M • -σ relation is used with the VDF. The one place this paper deviates from holding many things constant is when determining which M

Calculating the GWB
As in SB16, Eqs. 1, 3, and 7 are combined to calculate h c . When using a GSMF, When using the VDF, as described in §2.2, (10) To calculate A 2 yr , the above integrals are used setting f = f yr .
The limits of integration for the above equations are chosen in order to capture the major galaxy mergers in the local Universe, which make up the majority of contributing systems to the GWB in the PTA band. For the GSMF (Eqn 9) in SB16, these are shown to be 0 < z < 1.5, 10 < log 10 M Gal < 12.5, and 0.25 < q < 1. Each histogram includes 3000 predictions, which 1000 coming from each of the three SMBH -host galaxy relations used. The 90% interval of the VDF predictions cover 0.24 dex with a mean at log10 Ayr = −14.74, while the GSMF predictions cover 0.4 dex over the 90% interval with a mean at log10 Ayr = −14.9. Both methods are able to reproduce the emerging signal, but the VDF predictions tend towards higher amplitudes than the GSMF predictions.
When using the VDF, it is imperative to ensure that the integration happens over a simiilar range of galaxy mergers. To do so, this work uses the size -mass relation from Williams et al. (2010) because it is based on the same UKIDSS sample that the VDF in Bezanson et al. (2012) is derived from. The corresponding limits of integration for the VDF are then 0 < z < 1.5, 1.85 < log 10 σ < 2.6, 0.67 < q σ < 1. While these ranges are smaller, it is important to remember that the corresponding range of σ values is much narrower for a given range of M . Integration limits derived from the sizemass relation in Shen et al. (2003), which is based on SDSS data, were also tested, but produced almost no change in the VDF predictions for h c .

RESULTS
Combining all of the observational constraints from §2.3 into the models in Equation 9 and Equation 10, I calculate predictions for the GWB. For A yr predictions, I generate 1000 realizations each for the VDF and GSMF under each SMBH -host galaxy relation. In each realization, all parameters are allowed to vary within their observed uncertainties.  Figure 1 for all three relations is shown in gray. The VDF predictions appear to not be sensitive to the choice of host galaxy relation in the same way that the GSMF predictions are, and as seen in Figure 1 the VDF predictions are always at a higher amplitude. However, the magnitude of the difference between the VDF and GSMF predictions is dependent on the specific host galaxy relation chosen. The largest difference is seen when using MM13 and the smallest difference is seen when using KH13.  3. 2D differential contribution to A 2 yr from two different aspects of the simulated binary SMBH population: binary chirp mass (Mc) shown on the x-axis and redshift shown on the y-axis. The columns show the VDF inference on the left (in blue) and the GSMF inference on the right (in orange). The rows represent two different SMBH -host galaxy relations with MM13 on the top and KH13 on the bottom. The limits of the colorbars are the same in all cases. The biggest difference between the VDF and GSMF is the contribution to the signal at higher redshift (z > 0.5). In the case of the GSMF inference, the population of binaries that contribute significantly to the GWB is contained to the local Universe, while the VDF inference shows significant contributions all the way out to z ∼ 1 and beyond. One can infer from this that the underlying SMBH mass function implied by the VDF includes more massive systems at higher redshifts then that implied by the GSMF. Figure 1 shows the relative distributions predicted for A yr from the two different methods (GSMF and VDF) employed in this paper. The predictions based on the VDF are both at a higher amplitude than those from the GSMF and cover a smaller spread implying that the underlying binary SMBH populations inferred by these two methods are quite different. As discussed in §2.3, I have endeavored to minimize any differences in observational inputs for these predictions, in order to highlight that the subsequent A yr differences are due primarily to differences in the method used to infer the SMBH mass function (discussed in §2.1), and not simply due to using different observational inputs to the model. The larger values of A yr predicted by the VDF are more in line with the amplitude of the common process starting to appear in PTA data, however, both methods are able to reproduce the signal emerging in PTA datasets. To highlight this, Figure 1 also includes the A yr constraints from both the NANOGrav 12.5yr dataset (Arzoumanian et al. 2020), as well as IPTA DR 2 (Antoniadis et al. 2022). Yet, without the requisite spatial correlations the origin of the signal appearing in PTA data remains unclear.

Predicted Amplitude of the GWB
Interestingly, the VDF predictions for A yr appear to be less sensitive to the choice of SMBH -host galaxy relation than the GSMF method. In Figure 2, the three panels show the A yr predictions for each SMBH -host galaxy relation used in this work. In all three cases, the VDF predictions are nearly identical and almost perfectly match the combined distribution seen in Figure 1. In contrast, the mean value of the GSMF prediction changes significantly depending on which SMBH -host galaxy relation is utilized. This behavior of the GSMF predictions has been previously noted in SB16.
One of my aims for pursuing this new method of inferring the SMBH mass function was to reduce the uncertainty on A yr that is present under the GSMF method. Some of the reduction in uncertainty likely comes from removing the intermediate step of calculating f bulge for each galaxy. f bulge covered a broad range of values due to limited observational evidence, however, most systems that strongly contribute to the GWB likely have f bulge ∼ 1 (Sesana 2013;Kelley et al. 2017). Some of the reduction in uncertainty is due to the narrower spread of measurements for the M • -σ relation in comparison to the M • -M bulge relation. Either way, it appears that by including observations of both a galaxy's size and Sérsic index along with stellar mass in the inference of the central SMBH, the VDF provides a more direct avenue to the SMBH mass function.
To help understand where the differences in the A yr predictions are coming from, Figure 3 shows the 2D differential contribution to A 2 yr from M c and z. The most striking difference is the contribution to A 2 yr that the VDF receives at higher redshifts (z > 0.5). Figure 3 also shows that the additional contribution at higher redshift is present regardless of which SMBH -host galaxy relation is used with the VDF method. Additionally, one sees that the different SMBH -host galaxy relations impact the predicted SMBH masses from each method in opposite ways. For instance, the relations in MM13 predict higher SMBH masses when using M • -σ than when using M bulge , and that is reversed when the relations in KH13 are used, with the GSMF method producing more massive SMBH systems, especially in the local Universe (z < 0.2). It is obvious just from looking at the different scaling in the colormaps between the two GSMF models that the predicted amplitude using MM13 is lower than KH13. This is in contrast to the similarities in the amplitude predictions from the VDF models, which can be seen in Figure 2. Overall, one can infer that the SMBH mass function implied by the VDF includes more massive SMBHs at higher redshifts than the mass function implied by the GSMF.
The results from using the VDF model are similar to those found in Casey-Clyde et al. (2021), which uses a quasar-based model to predict the binary SMBH population. In that case, the SMBH mass function has more support at larger mass values than the GSMF models from Sesana (2013) that they compare to, the amplitude of the GWB is more in line with the NANOGrav 12.5yr dataset, and the redshift contribution from the quasar-based model is more evenly distributed like the VDF model shown here.

Spectroscopic vs. Inferred VDF
As discussed in §2.2, this paper primarily uses the inferred VDF from Bezanson et al. (2012), however, given the recent publication of the spectroscopic VDF observed with the LEGA-C survey (Taylor et al. 2022), I calculate A yr using both functions. Figure 4 shows the relative distributions predicted for A yr under each model, however unlike earlier comparisons, this one only uses the dN19 SMBH -host galaxy relation. Because spectroscopic observations are only available for z < 1, I use the same inferred VDF functions for 1 < z < 1.5, however, Figure 3 shows that this redshift range is not a significant contributor to the overall amplitude of the GWB.
The consistency between the spectroscopic and the inferred VDF predictions is promising, implying that the inferred VDF is a good proxy for the spectroscopic Ayr predictions calculated using the inferred VDF (σ inf , black line) and the spectroscopic VDF (σspec, red line) to infer the SMBH mass function. In the case of the σspec predictions, spectroscopic data from SDSS (Bernardi et al. 2010) and the LEGA-C survey (Taylor et al. 2022) were used at z < 1, while the inferred VDF was used to fill in from 1 < z < 1.5. Given the consistency that the LEGA-C survey finds between the spectroscopic and inferred VDFs (Taylor et al. 2022), it is not surprising that these results are so consistent. While it would be preferred to use only spectroscopic VDFs for this work, the similarity in these results show that the inferred VDFs are an acceptable alternative.
VDF. As shown in the LEGA-C survey results (Taylor et al. 2022), the spectroscopic VDF is on the higher end of the error region covered by the inferred VDF. Thus, it is not surprising that the spectroscopic VDF predictions are consistent with the higher end of the inferred VDF predictions, and the smaller spread is directly attributable to the smaller error bars on the spectroscopic VDF.

Prevalence of Individual Sources
Beyond the amplitude of the GWB, the presence of more massive systems at higher redshifts predicted by the VDF method also has implications for the prevalence of individually resolvable, or continuous wave (CW), sources in PTA datasets. These systems are unique from the GWB, which is compromised of the GW contributions from the entire population of SMBH binary systems in the Universe. CWs are individual systems and once one is detected, it will provide a unique multimessenger opportunity to study a binary SMBH system. The CWs detectable by PTAs are different than those detectable by other GW experiments, because these systems are typically a long way away from coalescence, and so their frequency evolution should be negligible, allow-  Figure 5. Redshift distribution of individual binary SMBH systems at fGW = 8 nHz averaged over 100 different population realizations. On average, the VDF method (blue curve) predicts more than twice as many potential individual systems as the GSMF method (orange curve). These systems are expected to be detectable for the IPTA by 2025 (Xin et al. 2021). As seen in previous work, the GSMF method predicts that single sources will be found in the local Universe (z < 0.5) (Rosado et al. 2015). In contrast, the VDF method shows only a slight preference for systems at lower redshift. As seen in Figure 3, the larger mass systems at higher redshifts that are present in the VDF predictions directly translate to an increase in individually resolvable sources in PTA data for z > 0.5.
ing PTAs to probe a part of the SMBH binary evolution that is complimentary to other GW experiments. Figure 5 shows the redshift distribution of individual binary SMBH systems from one hundred different population realizations using the KH13 SMBH -host galaxy relations. I use this host galaxy relation since it predicts the most similar amplitudes of the GWB between the two methods, and CW prevalence is directly linked to GWB amplitude (Rosado et al. 2015). In each case, I have only looked at individual systems that are emitting GWs at a GW frequency close to 8 nHz in order to compare to the expected sky-averaged sensitivity of the IPTA by 2025 (Xin et al. 2021). On average, the VDF method predicts more than twice as many potential CWs as the GSMF, and the redshift distribution of sources is much broader with the VDF method than with the GSMF. This is to be expected, given the differences in the underlying SMBH populations seen in Figure 3.
The GSMF results shown here are very similar to those in Rosado et al. (2015), which finds that a GWB is more likely to detected before an individually resolvable source and that the first detected individual binary system is likely to be very close by (z 0.5). The VDF results are much more optimistic due to the increase in massive SMBH systems at higher redshift. Notably, the VDF method predicts that potentially resolvable sources are fairly evenly distributed across redshift, which is in contrast to other CW predictions, which find that the first detected system is likely to be very near by (Rosado et al. 2015;Kelley et al. 2018). The larger systems at higher redshifts have intriguing multi-messenger potential since if they were to have any variable AGN activity they would be promising multi-messenger candidates given near-future time-domain surveys (Charisi et al. 2021).

SUMMARY
In this paper, I undertake a new approach to inferring the SMBH mass function using observations of an inferred VDF and compare that to the standard inference method from a GSMF. I also compare the inferred VDF to the spectroscopic VDFs that are available within limited redshift ranges. I am careful to use as many of the same observations as possible between the two methods to ensure that the results are solely due to the different method of inference rather than different observational inputs. However, in order to not be biased by the choice of SMBH -host galaxy relation, which is arguably one of the most impactful parameters (Sesana 2013, SB16), I use three separate measurements that span the range of observed values. The VDF method predicts a binary SMBH population that is significantly more massive at higher redshifts than the GSMF method. This population is more inline with observed relic galaxies, like NGC1271 and NGC1277 (Walsh et al. 2015(Walsh et al. , 2016.
Interestingly, the VDF method produces A yr predictions that are at larger amplitudes and cover a narrower spread than the GSMF method, although both methods are able to reproduce the amplitude of the emerging signal in PTA data. The narrower spread in the VDF method's predictions is partially due to the smaller range of M • -σ measurements, compared to the M • -M bulge relation. This provides further evidence that the M • -σ relation is a more "fundamental" relation between host galaxy's and their SMBHs.
The largest difference between the two methods is from massive binaries at higher redshifts, where the VDF predicts significant contributions, while the GSMF does not. Additionally, the more massive binary SMBHs at higher redshifts predicted by the VDF method increase the number of potential individually resolvable systems in PTA data as well as broaden the redshift range where those systems are expected to be found. This is inline with recent studies of z ∼ 1 galaxy clusters, which appear to host potential CW systems (Wendt & Romani 2023).
It is worth noting that the results presented here are a preliminary investigation into using a VDF to infer the SMBH mass function. The fits to the VDF in Bezanson et al. (2012) are estimates and a full error analysis has yet to be conducted on the inferred velocity dispersion functions. However, the consistency between the LEGA-C spectroscopic VDF and inferred VDFs is a promising development. Additionally, the VDFs used in this paper only went out to z = 1.5 and as Figure 3 shows, there is reason to think that significant contributions to the GWB from the VDF method may lie at z > 1.5. Ideally, the VDF would reach out to z ∼ 2 − 3, and work is currently underway to provide it (Matt et al. in prep).
Finally, I want to reiterate some of the other initial caveats to this work. The model presented here is overly simplistic, looking only at binary SMBHs in circular orbits, ignoring environmental coupling, and assuming a quick, uniform binary evolution timescale. However, these results show that a more direct method of SMBH mass inference is available to binary SMBH populations based on galaxy observations and that the inferred SMBH mass function from this new method is fundamentally different than what has been used in the past. While more work is needed, these results hint towards there being a population of more massive SMBHs at higher redshift than previously thought.
I am grateful to the referee, Alberto Sesana, for comments and insights that greatly improved this work. This project was conceived and advanced through fruitful discussions with Sarah Burke-Spolaor, Julie Comerford, Joseph Lazio, Xavier Siemens, and Sarah Vigeland. The manuscript was improved by comments from Luke Kelley and Kayhan Gültekin. I am grateful to Rachel Bezanson for providing key insights to the velocity dispersion function. This work has been supported by the National Science Foundation (NSF) under award 1847938. Additionally, I am supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-2202388.