The NANOGrav 15 yr data set shows evidence for the presence of a lowfrequency gravitationalwave background GWB. While many physical processes can source such lowfrequency gravitational waves here we analyze the signal as coming from a population of supermassive black hole SMBH binaries distributed throughout the Universe. We show that astrophysically motivated models of SMBH binary populations are able to reproduce both the amplitude and shape of the observed lowfrequency gravitationalwave spectrum. While multiple model variations are able to reproduce the GWB spectrum at our current measurement precision our results highlight the importance of accurately modeling binary evolution for producing realistic GWB spectra. Additionally while reasonable parameters are able to reproduce the 15 yr observations the implied GWB amplitude necessitates either a large number of parameters to be at the edges of expected values or a small number of parameters to be notably different from standard expectations. While we are not yet able to definitively establish the origin of the inferred GWB signal the consistency of the signal with astrophysical expectations offers a tantalizing prospect for confirming that SMBH binaries are able to form reach subparsec separations and eventually coalesce. As the significance grows over time higherorder features of the GWB spectrum will definitively determine the nature of the GWB and allow for novel constraints on SMBH populations.
['Gravitational waves', 'Supermassive black holes', 'Galaxy evolution', '678', '1663', '594', 'Astrophysics - High Energy Astrophysical Phenomena', 'Astrophysics - Cosmology and Nongalactic Astrophysics', 'General Relativity and Quantum Cosmology']
The NANOGrav 15 yr Data Set Constraints on Supermassive Black Hole Binaries from the Gravitationalwave Background
{'Constraints on Supermassive Black Hole Binaries from the Gravitational Wave Background': "G abriella A gazie , 1 A kash A numarlapudi , 1 A nne M. A rchibald , 2 P aul T. B aker , 3 B ence B' ecsy , 4 L aura B lecha , 5 A lexander B onilla , 6, 7 A dam B razier , 8, 9 P aul R. B rook , 10 S arah B urke -S polaor , 11, 12 R and B urnette , 4 R obin C ase , 4 J. A ndrew C asey -C lyde , 13 M aria C harisi , 14 S hami C hatterjee , 8 K aterina C hatziioannou , 15 B elinda D. C heeseboro , 11, 12 S iyuan C hen , 16 T yler C ohen , 17 J ames M. C ordes , 8 N eil J. C ornish , 18 F ronefield C rawford , 19 H. T hankful C romartie , 8, ∗ K athryn C rowter , 20 C urt J. C utler , 21, 15 D aniel J. D'O razio , 22 M egan E. D e C esar , 23 D allas D e G an , 4 P aul B. D emorest , 24 H eling D eng , 4 T imothy D olch , 25, 26 B rendan D rachler , 27, 28 E lizabeth C. F errara , 29, 30, 31 W illiam F iore , 11, 12 E mmanuel F onseca , 11, 12 G abriel E. F reedman , 1 E miko G ardiner , 32 N ate G arver -D aniels , 11, 12 P eter A. G entile , 11, 12 K yle A. G ersbach , 14 J oseph G laser , 11, 12 D eborah C. G ood , 33, 13 K ayhan Gultekin , 34 J effrey S. H azboun , 4 S ophie H ourihane , 15 K ristina I slo , 1 R oss J. J ennings , 11, 12, † A aron J ohnson , 1, 35 M egan L. J ones , 1 A ndrew R. K aiser , 11, 12 D avid L. K aplan , 1 L uke Z oltan K elley , 32 M atthew K err , 36 J oey S. K ey , 37 N ima L aal , 4 M ichael T. L am , 27, 28 W illiam G. L amb , 14 T. J oseph W. L azio , 21 N atalia L ewandowska , 38 T yson B. L ittenberg , 39 T ingting L iu , 11, 12 J ing L uo , 40, ‡ R yan S. L ynch , 41 C hung -P ei M a , 32, 42 D ustin R. M adison , 43 A lexander M c E wen , 1 J ames W. M c K ee , 44, 45 M aura A. M c L aughlin , 11, 12 N atasha M c M ann , 14 B radley W. M eyers , 20, 46 P atrick M. M eyers , 15 C hiara M. F. M ingarelli , 33, 13, 47 A ndrea M itridate , 48 P riyamvada N atarajan , 49, 50 C herry N g , 51 D avid J. N ice , 52 S tella K och O cker , 8 K en D. O lum , 53 T imothy T. P ennucci , 54 B enetge B. P. P erera , 55 P olina P etrov , 14 N ihan S. P ol , 14 H enri A. R adovan , 56 S cott M. R ansom , 57 P aul S. R ay , 36 J oseph D. R omano , 58 J essie C. R unnoe , 14 S hashwat C. S ardesai , 1 A nn S chmiedekamp , 59 C arl S chmiedekamp , 59 K ai S chmitz , 60 L evi S chult , 14 B rent J. S hapiro -A lbert , 11, 12 X avier S iemens , 4, 1 J oseph S imon , 61, § M agdalena S. S iwek , 62 I ngrid H. S tairs , 20 D aniel R. S tinebring , 63 K evin S tovall , 24 J erry P. S un , 4 A bhimanyu S usobhanan , 1 J oseph K. S wiggum , 52, † J acob T aylor , 4 S tephen R. T aylor , 14 J acob E. T urner , 11, 12 C aner U nal , 64, 65 M ichele V allisneri , 21, 35 S arah J. V igeland , 1 J eremy M. W achter , 66 H aley M. W ahl , 11, 12 Q iaohong W ang , 14 C aitlin A. W itt , 67, 68, 12 D avid W right , 69 O livia Y oung , 27, 28", 'T he NANOG rav C ollaboration': "Table B1. Astrophysical parameters of our semi-analytic population models. Units are denoted where relevant; all other parameters are defined to be dimensionless. For libraries, we denote uniform distributions with U (min , max) and normal distributions with N (mean , std. dev.).Table B2. Summary of semi-analytic SMBH binary population models used in this work. Model parameters are defined in § 3.2, and their fiducial values and assumed prior distributions are given in Table B1. The first model ( Phenom + Uniform ), indicated in boldface, is what we refer to as our fiducial model. Throughout the text, when the assumed priors can be omitted from the model name without loss of clarity, we simply refer to the models as Phenom , GWOnly , Phenom-Ext , and GWOnly-Ext . The latter two (also referred to as the 'extended models') are discussed in Appendix C.  \nFigure C1. Semi-Analytic Model posteriors comparing our fiducial six-dimensional Phenom model (solid) to our much larger Phenom-Ext parameter space (dashed). The standard version of the Phenom library with uniform priors (blue) is also compared to a version which is fit against the 15 yr GWB spectra using the same astrophysically motivated priors from the Phenom + Astro distributions (orange). The priors are also shown for the astrophysically motivated case (grey).  \n<!-- image -->  \nthe library grid points themselves. In this way we can weight the input parameters by the resulting likelihoods to obtain posteriors without using MCMC to dynamically explore the domain. This approach allows us to examine the e ff ects of freeing additional parameters.  \nFigure C1 compares the parameter posteriors for our fiducial library against the larger parameter spaces, including a twelvedimensional phenomenological version ( Phenom-Ext , orange dashed), and ten-dimensional GW-only version ( GWOnly-Ext , purple dashed). Blue solid lines show the standard Phenom fits with uniform priors, while orange solid lines show the same six-dimensional library, but fitting with priors taken from the Phenom-Ext distributions.  \nPosteriors from the six-dimensional phenomenological models are entirely consistent with the ten- & twelve-dimensional Phenom-Ext libraries when using the same priors. Setting the additional parameters to fixed values does not bias the resulting measurements, nor does it lead to under-estimating the width of posterior distributions. This is likely the case because the current 15 yr NANOGrav spectral measurements include large uncertainties. As the data improves, it will become more important to fully explore the parameter space.  \nAs previously discussed, for many parameters the shape of the priors does significantly impact the recovered posteriors. This is particularly noticeable in the GSMF parameters ( ψ 0, and m ψ, 0), and the M BHM bulge parameters ( µ , and ϵµ ) where there is significant degeneracy between parameters that broadly change the amplitude of the GWB spectrum.", 'ABSTRACT': 'The NANOGrav 15 yr data set shows evidence for the presence of a low-frequency gravitational-wave background (GWB). While many physical processes can source such low-frequency gravitational waves, here we analyze the signal as coming from a population of supermassive black hole (SMBH) binaries distributed throughout the Universe. We show that astrophysically motivated models of SMBH binary populations are able to reproduce both the amplitude and shape of the observed low-frequency gravitational-wave spectrum. While multiple model variations are able to reproduce the GWB spectrum at our current measurement precision, our  \n45  \n67  \nresults highlight the importance of accurately modeling binary evolution for producing realistic GWB spectra. Additionally, while reasonable parameters are able to reproduce the 15 yr observations, the implied GWB amplitude necessitates either a large number of parameters to be at the edges of expected values, or a small number of parameters to be notably di ff erent from standard expectations. While we are not yet able to definitively establish the origin of the inferred GWB signal, the consistency of the signal with astrophysical expectations o ff ers a tantalizing prospect for confirming that SMBH binaries are able to form, reach sub-parsec separations, and eventually coalesce. As the significance grows over time, higher-order features of the GWB spectrum will definitively determine the nature of the GWB and allow for novel constraints on SMBH populations.  \nKeywords: Gravitational Waves (678) - Supermassive black holes (1663) - Galaxy evolution (594)', '1. INTRODUCTION': "Strong observational evidence suggests that most, if not all, massive galaxies contain supermassive black holes (SMBHs) at their centers (Richstone et al. 1998). Additionally, hierarchical structure formation causes frequent galaxy mergers (Ostriker & Hausman 1977; White 1980; Lacey & Cole 1993), naturally leading to the formation of SMBH binaries, which may also merge (Begelman et al. 1980; Milosavljevi'c & Merritt 2001). At the last stages of their evolution, these binaries produce strong nanohertz gravitational wave (GW) emission that can be targeted by pulsar timing arrays (PTAs), which systematically monitor a large number of millisecond pulsars. By detecting coherent deviations in the times of arrival of pulsar signals, PTAs can observe a stochastic gravitational wave background (GWB) from the superposition of many unresolved binaries, as well as individually resolved sources on top of the background (Burke-Spolaor et al. 2019; Taylor 2021).  \nThe North American Nanohertz Observatory for Gravitational Waves (NANOGrav) 12.5-year data set showed evidence of a common-spectrum red noise process consistent with a GWB (Arzoumanian et al. 2020). This result was confirmed by the Parkes Pulsar Timing Array (PPTA; Goncharov et al. 2021), the European Pulsar Timing Array (EPTA; Chen et al. 2021), and the International Pulsar Timing Array (IPTA; Antoniadis et al. 2022). The NANOGrav 15 yr data set shows that the common uncorrelated red noise (CURN) signal discovered in Arzoumanian et al. (2020) persists with greater significance and is now detected in a larger number of pulsars (Agazie et al. 2023b, hereafter NG15gwb). Additionally, for the first time, there is evidence of inter-pulsar correlations following the characteristic Hellings-Downs (HD) pattern (Hellings & Downs 1983) expected for an isotropic  \nGWB. Careful analyses of the detection significance give false-alarm probabilities of ≈ 10 -4 -10 -3 ( ≈ 3 σ ).  \nIn this paper, we investigate whether the NANOGrav 15 yr results can be explained as a stochastic GWB produced by a cosmic population of SMBH binaries. While SMBH binaries have long been expected to produce such a background, a wide variety of alternative models exist, many of which invoke new physics that departs from the standard model and Λ Cold Dark Matter ( Λ CDM) cosmology. We refer the reader to Afzal et al. (2023, hereafter NG15newphys) for an analysis of the NANOGrav 15 yr results in the context of new-physics models, such as cosmic inflation, scalarinduced GWs, domain walls, cosmic strings, and first-order phase transitions.", '1.1. The galaxy-SMBH connection': 'Our understanding of galaxy formation and evolution has rapidly progressed in the last few decades. This includes the definitive and now direct observation of SMBHs in galaxy centers (Ghez et al. 1998; GRAVITY Collaboration et al. 2018; Event Horizon Telescope Collaboration et al. 2019, 2022). The mass of the central SMBH strongly correlates with global properties of the host galaxy (e.g., the stellar velocity dispersion of the galactic bulge, the bulge mass and luminosity), with tight correlations spanning several orders of magnitude in SMBH mass (Dressler 1989; Kormendy 1993; Magorrian et al. 1998; Gebhardt et al. 2000; Tremaine et al. 2002; Haring & Rix 2004; Gultekin et al. 2009a; Kormendy & Ho 2013; McConnell & Ma 2013; Saglia et al. 2016). These trends strongly imply coordinated evolution between SMBHs and their host galaxies, which may be driven by a variety of mechanisms such as galaxy mergers, secular dynamics, stellar feedback, and feedback from active galactic nuclei (AGN; Di Matteo et al. 2005; Hopkins et al. 2008; Somerville et al. 2008). SMBHs are believed to play particularly significant roles in shaping the structure of massive galaxies (Croton et al. 2006; Fabian 2012; Vogelsberger et al. 2014; Schaye et al. 2015; Weinberger et al. 2017), but many fundamental aspects, such as the formation channels  \nof SMBH seeds in the early Universe or how AGN feedback  \nshapes the host galaxies, are still poorly constrained via observations. The relevant physical processes are also very difficult to model theoretically, as they span size scales from galaxies ( ∼ 10 kpc) to SMBH event horizons ( ∼ 10 -5 pc). Similar challenges limit our ability to directly model the process of SMBH binary formation and evolution.', '1.2. SMBH binary evolution': "The formation of SMBH binaries begins with the merger of two galaxies, each hosting a central SMBH. At di ff erent stages of the evolution of the SMBH pair, di ff erent physical processes dominate energy and angular momentum extraction, which drives the binary to closer separations (Begelman et al. 1980; seeDe Rosa et al. 2019 for a recent review). Initially, the SMBHs are a gravitationally unbound pair (a dual SMBH) falling towards the center of the merging host (Barnes & Hernquist 1992) via dissipative 'hardening' processes, such as dynamical friction (Chandrasekhar 1943; Antonini & Merritt 2012). Once the mass enclosed within the orbit is comparable to the mass of the binary (typically at ∼ parsec-scale separations), the two black holes become a gravitationally bound pair (a SMBH binary; Merritt & Milosavljevi'c 2005). At these separations, the timescale for GW-driven inspiral is generally still longer than the Hubble time, and their GW frequencies are orders of magnitude below those that PTAs can probe.  \nThe astrophysical environment of the binary is therefore crucial for bringing these systems to the PTA bands and ultimately to their final coalescence. Scattering of individual stars that pass close to the SMBHs can extract energy and angular momentum from the system, hardening the binary orbit (Yu 2002). In some cases, the supply of stars on close orbits may be insu ffi cient, and the binary would fail to merge within a Hubble time (Begelman et al. 1980). However, this so-called 'final-parsec problem' has a number of potential theoretical solutions (e.g., Berczik et al. 2006; HolleyBockelmann & Sigurdsson 2006; Khan et al. 2011; HolleyBockelmann & Khan 2015). Similarly, in gas-rich systems, circumbinary gas disks can also catalyze the binary evolution (Escala et al. 2005; Dotti et al. 2007; Haiman et al. 2009), but the e ffi ciency of this process or whether the gas pushes the binary inward or outward is still unclear (Mu˜noz et al. 2019; Moody et al. 2019; Du ff ell et al. 2020; Siwek et al. 2023).  \nIf a binary stalls for longer than the time between successive galaxy mergers, a second galaxy could bring a third SMBH into the system. Triple SMBH interactions can greatly reduce the timescale for a SMBH binary merger and may also cause the ejection of the lightest SMBH from the system (Saslaw et al. 1974; Volonteri et al. 2003a; Ho ff man & Loeb 2007; Bonetti et al. 2016, 2018a). Once a SMBH binary reaches a su ffi ciently small separation the GWs will dominate its evolution, carrying away energy and angular  \nmomentum and leading the SMBHs to coalescence (Peters &Mathews 1963).", '1.3. Electromagnetic signatures of SMBH binaries and multi-messenger prospects': "Many studies have used electromagnetic observations of AGN to find candidate SMBH pairs and binaries (for reviews, see: Komossa 2006; Popovi'c 2012; De Rosa et al. 2019; Bogdanovi'c et al. 2022). Dual AGN, i.e., galaxies with two unbound, actively accreting SMBHs, have been identified at ≳ kiloparsec separations (e.g., Koss et al. 2012; Chen et al. 2022, and references therein). However, spatially resolving the two SMBHs becomes increasingly challenging as their separation decreases. Spectroscopic features, such as the kinematic o ff set of AGN narrow lines, can also be used to identify AGN in merging galaxies (e.g. Comerford et al. 2009; Comerford & Greene 2014). To date, only one parsecscale pair has been confirmed with very long baseline interferometry (VLBI; Rodriguez et al. 2006; Bansal et al. 2017) despite large-scale searches (Burke-Spolaor 2011; Breiding et al. 2021).  \nElectromagnetic searches for sub-parsec SMBH binaries typically focus on features that encode the binary's orbital motion on the temporal or spectral variability of AGN. Searches for o ff set broad emission lines have been used to identify several hundred candidates (Tsalmantza et al. 2011; Eracleous et al. 2012; Shen et al. 2013; Ju et al. 2013), but this method is subject to false positives and other limitations (Gezari et al. 2007; Runnoe et al. 2015, 2017; Pflueger et al. 2018; Kelley 2021). Periodically variable light curves (Farris et al. 2014; D'Orazio et al. 2015; Bowen et al. 2018; D'Orazio & Di Stefano 2018) have yielded a similar number of candidates (Graham et al. 2015; Charisi et al. 2016; Liu et al. 2019), though these samples are likely to also suffer significant contamination (Vaughan et al. 2016; Charisi et al. 2018; Sesana et al. 2018; Kelley et al. 2019b; Xin et al. 2020). Despite these challenges, the advent of large time domain surveys with the Vera Rubin observatory (Ivezi'c et al. 2019), combined with multi-wavelength observations and increasing PTA sensitivity to the GWB, o ff ers exciting opportunities for deriving multiple independent constraints on SMBH populations (e.g. Kelley et al. 2019a; Bogdanovi'c et al. 2022). The prospects for low-frequency multimessenger astrophysics are discussed further in § 5.", '1.4. The Astrophysical Imprint on the Gravitational Wave Background': 'All of the binary inspiral processes discussed above are imprinted on the GWB created by a population of SMBH binaries. Therefore, studying the GWB constitutes an important channel to obtain significant and novel insights on galaxy and binary mergers. For example, interactions with the binary environment and orbital eccentricities impact the shape of the  \nGWB spectrum (Sesana 2013a). Stellar- and gas-driven binary hardening will cause a flattening or turnover of the lowfrequency GWB spectrum, relative to the single power law predicted for GW-only evolution (Kocsis & Sesana 2011). The primary e ff ect of eccentricity is to boost GW emission to higher frequencies, owing to the emission of GWs at higher harmonics beyond twice the binary orbital frequency, which dominates for circular orbits (Enoki & Nagashima 2007). However, at extreme eccentricities ( ≳ 0 . 9 -0 . 95), close pericentric passages drive very rapid binary inspiral, leading to an overall attenuation of GWB amplitude at all frequencies (e.g., Kelley et al. 2017b).  \nGW observations will also probe the history of SMBH mass growth. The GWB depends strongly on the distribution of binary chirp masses, M , given by  \nM = ( m 1 m 2) 3 / 5 M 1 / 5 = M q 3 / 5 (1 + q ) 6 / 5 , (1)  \nwhere q ≡ m 2 / m 1 < 1 is the binary mass ratio, M = m 1 + m 2 is the total binary mass, and m 1 and m 2 are the masses of each SMBH. As a result, the GWB is intimately related to the SMBH mass function through its dependence on the chirp mass, and in turn to the scaling relations of SMBH mass with host galaxy properties. These relations are well studied in the local Universe but are unconstrained at higher redshifts, and thus detailed studies of the GWB will provide a novel path to probing these relations.  \nPrevious stochastic GWB constraints have been used to probe the SMBH binary population, by comparing with theoretical predictions for SMBH binary formation and evolution. All GWB results were strictly upper limits until the NANOGrav 12.5-year data set, but the limits were still potentially constraining. The constraints were especially informative when combined with electromagnetic observations of binary AGN candidates (Sesana et al. 2018; Holgado et al. 2018; Inayoshi et al. 2018; Nguyen et al. 2020). After the PPTA upper limit at 2.8 nHz (Shannon et al. 2013), it was first suggested that this ruled out a large range of SMBHbinary model space (Shannon et al. 2015). However, Middleton et al. (2016) showed that the upper limits were consistent with a wide variety of plausible astrophysical models and that, in general, upper limits alone would be relatively unconstraining until they were about an order of magnitude smaller. Subsequent work showed the importance of analyzing many pulsars and accounting for their red noise and systematic errors in solar-system ephemerides when establishing PTA upper limits (Arzoumanian et al. 2018; Vallisneri et al. 2020; Hazboun et al. 2020; Johnson et al. 2022). Since the 12.5-year NANOGrav data set showed evidence for a common red-noise process consistent with (but not unambiguously attributable to) GWs, the measurement was shown  \nto be consistent with a population of SMBH binaries with reasonable properties (Middleton et al. 2021).', '1.5. Astrophysical modeling of SMBH-binary GWB': 'Over the last few decades, a variety of di ff erent approaches have been used to model populations of SMBH binaries, 1 with a wide range of predictions for the resulting GWB amplitude. (See Appendix A for a summary of these model predictions and a comparison with the NANOGrav 15 yr results). Many of these studies start from either semi-analytic galaxy evolution models to obtain galaxy merger rates (Rajagopal & Romani 1995) or halo merger-trees with added galaxies (Menou et al. 2001; Sesana et al. 2004), onto which a SMBH binary population model can be imposed. In lieu of physically modeling environmentally driven SMBH binary evolution, galaxy mergers are often directly linked to the formation of a close SMBH binary emitting GWs at PTA frequencies, and a power-law form is assumed for the GWB (e.g., Phinney 2001; Ja ff e & Backer 2003; Wyithe & Loeb 2003; Enoki et al. 2004; Simon & Burke-Spolaor 2016). Some semi-analytic models also include prescriptions for physical processes that cause GWB spectra to deviate from a pure power law, such as interactions of the binary with the gaseous and stellar environment of its host galaxy, discreteness of the binary population, and orbital eccentricity (e.g., Sesana et al. 2008, 2009a; Sesana 2013b; Ravi et al. 2014; McWilliams et al. 2014; Ryu et al. 2018; Bonetti et al. 2018b; Chen et al. 2020). Versions of the semi-analytic model approach have also been applied to catalogs of specific galaxies or quasars from observations (Simon et al. 2014; Rosado & Sesana 2014; Mingarelli et al. 2017; Casey-Clyde et al. 2022).  \nAn alternative to the semi-analytic modeling approach is to trace galaxy and SMBH evolution directly in cosmological hydrodynamics simulations (e.g., Kulier et al. 2015; Salcido et al. 2016; Kelley et al. 2017a,b, 2018; Volonteri et al. 2020; Siwek et al. 2020; Curyło & Bulik 2022). This approach has the advantage of providing detailed information about the internal structures of galaxies and how they interact with SMBHs via AGN fueling and feedback. However, cosmological hydrodynamical simulations are very computationally expensive compared to semi-analytic models, and even the highest-resolution simulations must rely on sub-grid prescriptions to model unresolved processes, including SMBH accretion, mergers, and feedback. Each of these complementary approaches therefore o ff ers benefits and drawbacks, and importantly, each introduces certain systematics in their predictions for binary populations. In this work, we adopt  \nFigure 1. NANOGrav 15 yr GWB free-spectrum posteriors translated into the square root of timing residual power ( ρ , left panel) and characteristic strain ( hc , right panel). The HD-correlated free-spectrum measured while simultaneously fitting for monopole-correlated (MP), dipole-correlated (DP), and uncorrelated red noise noise (CURN) free spectra ( HD-w / MP + DP + CURN ; gray violins, left-side) is compared against the HD-DMGP model in which dispersion-measure variations are modeled using Gaussian Processes (green violins, right-side). The black dotted lines show idealized power-law spectra ( ρ 2 ∝ f -13 / 3 and hc ∝ f -2 / 3 ) fit to the median posterior value for the amplitude obtained from the HD-w / MP + DP + CURN model in NG15gwb. Over plotted are the best-fitting, simulated GWB spectra from models of SMBH binary populations produced in this analysis. Two models are shown, one which includes environmentally driven binary evolution (blue) and another that assumes GW-only evolution (purple). Both models are able to reproduce the data, while the environmentally driven model produces a slightly better fit. We conclude that the observed GWB spectrum is consistent with astrophysically motivated expectations from populations of SMBH binaries.  \n<!-- image -->  \na semi-analytic modeling approach to SMBH binary population synthesis and defer the use of cosmological hydrodynamics simulations for future work.', 'SUMMARY & OUTLINE': "Figure 1 shows the GWB spectrum recovered from the 15 yr NANOGrav data, along with the best fitting simulated GWB spectra produced in this work. In § 2 we summarize the NANOGrav 15 yr data set that forms the observational basis for this analysis, and the GWB spectra derived from it (grey and green 'violins'). In § 3, we describe our methods of modeling populations of SMBH binaries and calculating the GWB spectra that they would produce. There, we also detail the approach that we use to compare our simulations to the 15 yr data. Our best-fitting models (colored curves) are presented in § 4.  \nWe find that astrophysically motivated models of SMBH binary populations are able to accurate reproduce the observed GWB spectrum (§ 4.1 & 4.2). We focus our analysis on two population models. One includes a self consistent prescription for environmentally driven binary evolution (blue), and the other assumes GW-only evolution (purple) which is still commonly used in the literature. Both models are able to fit the data, while the environmentally driven case produces a slightly better match-particularly to the lowest frequency bin. We present the binary evolution parameters favored by 15 yr spectra fits for both models (§ 4.3). While the posterior distributions are broadly consistent with astrophysical  \nexpectations, parameters tend to be shifted towards values that produce larger GWB amplitudes than was previously most-favored. Generally higher binary masses or densities, or highly e ffi cient binary mergers are required to produce the observed amplitudes. The characteristics of the implied binary populations are presented in § 4.4.  \nOur results are discussed in the context of the field in § 5, along with highlights for the near future of low-frequency GWastronomy.  \nThroughout this paper we assume a WMAP9 cosmology with Ω m = 0 . 228, Ω b = 0 . 0472, and H 0 = 0 . 6933 km s -1 Mpc -1 .", '2. PULSAR TIMING ARRAY DATA': "This work is based on the NANOGrav 15 yr data set, which includes 68 pulsars, 67 of which have a baseline of at least 3 years and are included in the GWB analysis. The complete description of the data set can be found in Agazie et al. (2023a, hereafter NG15), while the detector characterization and noise modeling of individual pulsars is described in Agazie et al. (2023c, hereafter NG15detchar). The detailed description of the Bayesian search for the GWB is presented in NG15gwb. Here, we briefly summarize the measurement of the GWB spectrum from the NANOGrav data, focusing on the pieces which are necessary for the astrophysical interpretation presented in this paper.  \nPTA collaborations systematically monitor millisecond pulsars and record the times of arrival (TOAs) of their radio  \npulses. For each pulsar, a timing model is constructed, which estimates various factors a ff ecting the TOAs including its astrometry (sky position, proper motion, and parallax), its spin period and spin period derivative, and binary parameters for pulsars with companions. Additionally, variations in the ionized interstellar medium along the line of sight, also known as the dispersion measure (DM), are included in our model. The analysis of each pulsar provides a best-fit estimate for the timing residuals, r ( t ), which are the di ff erences between the TOAs and the timing model. For more on the construction of the timing residuals in NANOGrav's data set, see § 4 of NG15.  \nAll red noise processes, including the GWB itself, are modeled with a Fourier basis computed on the TOAs, as discussed in § 2 of NG15gwb. The frequencies are fi = i / T obs, where T obs = 16 . 03 yr is the time between the first and last TOAincluded in this data set 2 . The search for a GW signal is performed by constructing the cross-correlations of residuals between pairs of pulsars, a and b , i.e.,  \n⟨ ra ( t ) rb ( t ) ⟩ ∝ Z Sab ( f ) d f , (2)  \nwhere f is the observer-frame GW frequency and Sab is the timing-residual cross-correlated power spectral density,  \nSab ( f ) = Γ ( ξ ab ) Φ ( f ) . (3)  \nHere, Φ ( f ) is the power spectral density (PSD) of the timing residuals describing the spectrum of the process that is common among all pulsars, and Γ is the overlap reduction function, which describes the induced correlation between a pair of pulsars as a function of their angular separation, ξ ab . The timing residual PSD is related to the characteristic GW strain, hc ( f ), by  \nΦ ( f ) = hc ( f ) 2 12 π 2 f 3 . (4)  \nThe overlap reduction function is given by Γ ( ξ ab ) = δ ab for a CURN model, and by the characteristic HD pattern in the case of an isotropic GWB (Hellings & Downs 1983).  \nAspecific spectral shape is typically prescribed to the common red noise process (see § 3.1). Traditionally, a power law has been used, and the detailed spectral analysis presented in NG15gwb shows support for this idealized, simple model. However, deviations appear at a variety of frequencies, which may skew the determination of a spectral slope (see Figure 1 and Figure 6 in NG15gwb). It is therefore important to model the individual Fourier coe ffi cients independently rather than enforcing a specific spectral shape on the PSD. The resulting 'free spectrum' provides a minimally modeled Bayesian  \nspectral characterization of PTA data. The free spectrum recovers the posterior of the common red-noise power spectrum at all sampling frequencies, and is parameterized by the coe ffi cient ρ , where ρ 2 i = Φ ( fi ) / T obs is the power in the crosscorrelated timing residuals.  \nFigure 1 shows the free-spectrum posteriors both in terms of the square-root, timing residual power ( ρ ) and converted into GW characteristic strain ( hc ). At the current signal-tonoise ratio, the spectral characterization of the signal is uncertain, and the recovered HD-correlated GWB signal may be impacted by non-HD-correlated noise. For our astrophysical interpretation, we adopt the free-spectrum posteriors from the 15 yr HD-correlated free-spectrum modeled simultaneously with additional monopole-correlated (MP), dipole-correlated (DP), and uncorrelated red noise (CURN). This model, which we refer to as HD-w / MP + DP + CURN (grey), provides the most conservative constraints on the recovered GWB spectrum. As an additional comparison, we also analyze the HD-correlated free-spectrum posteriors utilizing an alternate model for DM variations, which we refer to as HD-DMGP (green)-described in detail in § 5.1 of NG15gwb.  \nFigure 1 also shows the median posterior amplitude value for the idealized power-law fit to the HD-w / MP + DP + CURN model. Even though the power-law model provides an illustrative example for quick model comparison, we do not include it in the astrophysical interpretation, because it fails to encapsulate the full range of information contained in the free-spectrum posteriors.  \nThe number of Fourier components used in an analysis is typically chosen based on the preference of the data for various red noise processes (e.g., in Arzoumanian et al. 2020 the CURN model preferred only 5 frequencies, while in NG15gwb that number has increased to 14). While the CURN model prefers 14 Fourier components in the 15 yr data set, the HD-correlated free spectra posteriors provide strong constraints only in the lowest 5 frequency bins, and thus only those bins are used in this analysis. However, we find no di ff erence in our results if we expand to using the full 14 frequencies.", '3. METHODS': 'Our goal is to constrain the properties of the underlying SMBH binary population that can produce a GWB consistent with the NANOGrav 15 yr data. Our approach consists of three main components described below, and depicted schematically in Figure 2.', 'SMBH Binary Population Synthesis Simulations (§ 3.1 - 3.3)': "-We generate 'libraries' of SMBH binary populations and their GW signals, exploring a large range of the binary formation / evolution parameter space. For this, NANOGrav has developed a flexible framework for SMBH binary popula-  \ntion synthesis called holodeck (Kelley et al. 2023, in prep), 3 which allows us to explore the binary population models and encompass systematic uncertainties. Within holodeck , we determine the number density of the cosmic population of SMBH binaries using semi-analytical models based on observationally constrained properties of galaxies and galaxy mergers. Using a SMBH-host relation, specifically the correlation between the mass of the SMBH and the mass of the stellar bulge, i.e. M BHM bulge, we assign SMBH masses to the mergers and calculate the binary evolution from large separations down to the GW regime. From each population, we compute the GWB signals they would produce.  \nInterpolation of the Population Synthesis Models (§ 3.4) -The simulated GWB spectra are sampled at discrete points of the multi-dimensional binary population parameter spaces that we explore. We refer to the collection of simulated spectra for a given parameter space as a 'library.' We then use Gaussian processes (GPs) to interpolate between the populationsynthesis simulations and predict the shape of the GWB spectrum for any point of the parameter domain. This is necessary as the population simulations are too computationally expensive to run live while fitting against the NANOGrav data.", 'Fitting Population Synthesis Models Against PTA Data (§ 3.5)': '-We use a Markov Chain Monte Carlo (MCMC) approach to fit the trained GPs against the input free-spectrum posteriors from NG15gwb, generating posterior distributions of the binary population model parameters. From these, we constrain the di ff erent SMBH binary populations and evolutionary scenarios that could produce the observed GWB.', '3.1. GWs from SMBH binary populations': "The GWB spectrum can be calculated as the integrated GWemission of individual binaries throughout the Universe. The characteristic strain of the GWB over a given logarithmic interval of frequency can be expressed as (Phinney 2001; Wyithe & Loeb 2003)  \nh 2 c ( f ) = Z dMdqdz ∂ 4 N ∂ M ∂ q ∂ z ∂ ln fp h 2 s ( fp ) . (5)  \nThe sky- and polarization-averaged GW spectral strain from a single, circular binary h s can be related to a binary's total GWluminosity, L GW as (Finn & Thorne 2000)  \nh 2 s ( f ) = G c 3 L GW GLYPH<16> 2 π fp GLYPH<17> 2 d 2 c = 32 5 c 8 ( G M ) 10 / 3 d 2 c GLYPH<16> 2 π fp GLYPH<17> 4 / 3 . (6)  \nHere, dc is the comoving distance to a source at redshift z . Because the GW frequency is twice the orbital frequency  \nFigure 2. Schematic of the analysis pipeline used in this study. Our population models are composed of galaxy mergers combined with SMBH masses and a prescription for binary evolution. Galaxy merger rates are derived from the combination of a galaxy stellarmass function (GSMF), galaxy pair fraction (GPF) and galaxy merger time (GMT). Binary evolution can either follow the selfconsistent, phenomenological approach, or it can assume GW-only evolution. Libraries contain a large number of binary populations and their resulting GWB signals, which are calculated for varying uncertain physical parameters in both evolution scenarios. Gaussian processes are used to interpolate across the library parameters when fitting against the NANOGrav 15 yr data. Fits can be performed using broad, uniform priors or with more tightly constrained, astrophysically motivated priors. Fits are performed against the 15 yr data using GWB spectra derived from the HD-w / MP + DP + CURN or HD-DMGP models (see text).  \n<!-- image -->  \nfor circular binaries, the observer-frame GW frequency f can be related to the rest-frame orbital frequency fp as f = 2 fp / (1 + z ). Throughout this paper, we take the chirp mass M (Equation 1), and by extension the total binary mass M , to be intrinsic rest-frame properties of the binary.  \nIn practice, it is much more convenient to calculate a comoving volumetric number density of binaries η ≡ dN / dVc , and use this quantity to infer the full population (Rajagopal &Romani 1995; Ja ff e & Backer 2003; Sesana et al. 2008)  \n∂ 4 N ∂ M ∂ q ∂ z ∂ ln fp = ∂ 3 η ∂ M ∂ q ∂ z ∂ t ∂ ln fp ∂ z ∂ t ∂ Vc ∂ z , = ∂ 3 η ∂ M ∂ q ∂ z · τ ( fp ) · 4 π c (1 + z ) d 2 c . (7)  \nHere τ ( fp ) ≡ fp / ( d fp / dt ) is the binary hardening timescale, the rest-frame duration that a binary spends in a given loga-  \nhmic interval of frequency. Equation (7) connects redshift evolution to the time-evolution of binary sources over frequencies. For a circular binary evolving purely due to GW emission, the rate of semi-major axis change and the hardening timescale are given by (Peters 1964)  \n∂ a ∂ t GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> gw = -64 G 3 5 c 5 m 1 m 2 M a 3 , τ gw ≡ ∂ t ∂ ln a = 5 96 G M c 3 ! -5 / 3 GLYPH<16> 2 π fp GLYPH<17> -8 / 3 . (8)  \nCombining the above equations with the comoving volume of a light-cone (e.g., Hogg 1999),  \n∂ Vc ∂ z ∂ z ∂ t = 4 π c (1 + z ) d 2 c , (9)  \ngives the idealized expression for a GWB produced by circular, GW-only driven SMBH binaries (Phinney 2001):  \nh 2 s,circ ( f ) = 4 π 3 c 2 (2 π f ) -4 / 3 (10) Z dMdqdz ∂ 3 η ∂ M ∂ q ∂ z ( G M ) 5 / 3 (1 + z ) 1 / 3 .  \nThis motivates the common expression for GWB spectra as a power law of the form,  \nh c( f ) = A yr · GLYPH<16> f / yr -1 GLYPH<17> -α , (11)  \nwhere A yr is the GWB amplitude referenced at a frequency of 1 yr -1 , and in the idealized case, α = 2 / 3. Because the timing-residual power spectral density of a GW signal is related to the characteristic GW strain by Equation (4), this ideal, power-law form of the GWB can be expressed relative to a reference frequency f ref equivalently as,  \nΦ ( f ) = A 2 12 π 2 f f ref ! -γ f -3 ref . (12)  \nNote that we have defined the power-law indices to be positive quantities such that hc ∝ f -α and Φ ∝ f -γ . The powerlaw indices are therefore related as γ = 3 + 2 α , such that the idealized, GW-only index is γ = 13 / 3.  \nRealistic GWB spectra can deviate substantially from a power law, primarily due to the following three e ff ects:  \nInteractions with the binary environment -Astrophysical processes that extract energy and angular momentum from the binary (e.g., via stellar and gaseous interactions) can accelerate its frequency evolution relative to the GW-only evolution. Therefore, any binary hardening via processes other than GW emission will necessarily result in an attenuation of the GWB compared to a purely GW-driven spectrum, as  \nbinaries spend less time emitting GWs in a given frequency interval. This e ff ect is particularly important at low frequencies ( f ≪ 1 yr -1 ) where binaries can more easily couple to their local galactic environments (Begelman et al. 1980; Kocsis & Sesana 2011), and where GW emission is weaker. In fact, coupling between SMBHs and their astrophysical environments is required for binaries to reach the PTA band within a Hubble time. The question is thus whether the resulting flattening (or turnover) in the GWB spectrum occurs within the PTA band or at frequencies too low to currently be accessible.  \nDiscreteness of the binary population -Equation (7) assumes a continuous distribution of SMBH binaries across the ( M , q , z , f ) parameter space. At low frequencies ( f ≲ 1 yr -1 ), the hardening timescale is very long, and a large number of binaries contribute to the GWB-making this approximation valid. At higher frequencies ( f ≳ 1 yr -1 ), however, the hardening timescale becomes shorter and the typical number of binaries producing the bulk of the GWB energy in a given frequency bin approaches unity (Sesana et al. 2008). In this regime, a continuous distribution overestimates the GWB signal. Properly accounting for the finite number of sources in each frequency bin therefore results in a steeper GWB spectrum at high frequencies (Ibid.). While a given overall amplitude of the GWB can be produced by either a larger number of lower-mass SMBH binaries or a smaller number of higher-mass binaries, these di ff erences change the frequency at which discreteness becomes important. As a result, they change the location and severity of the highfrequency spectral steepening.  \nOrbital eccentricity -Unlike circular binaries that emit GWs at exactly twice the orbital frequency, eccentric binaries emit GW energy at all integer harmonics. This leads to GW energy being moved from lower frequencies to higher frequencies (Enoki et al. 2004). Additionally, smaller pericenter distances tend to increase the rate of binary inspiral. These factors produce a variety of e ff ects, including a spectral turnover at low frequencies, a flatter spectrum at higher frequencies, and a 'bump' in between (Enoki & Nagashima 2007; Sesana 2013a; Huerta et al. 2015; Chen et al. 2017; Kelley et al. 2017b). However, for these e ff ects to be substantial, very large eccentricities ( e ≳ 0 . 9) are necessary at very small separations (well within the PTA band). 4 Since this is not ex-  \ncted to be the case, we restrict the current analysis to circular binaries.  \nThese e ff ects highlight the additional information encoded in the deviations of the GWB spectra from a pure power law and the importance of careful modeling of the binary population. The above considerations also demonstrate the need for explicit integration of the binary evolution that includes environmental interactions, the discreteness of binaries, and their expected cosmic variance.", '3.2. SMBH Binary Population Synthesis': 'Many previous works have constructed model populations of SMBH binaries and obtained predictions for the resulting GWB. These have produced a wide range of predictions for the GWB amplitudes, which are summarized in Table A1. These SMBH binary population models generally involve either semi-analytic models or cosmological hydrodynamic simulations. In this work, we focus only on semi-analytic models and defer the exploration of binary populations from cosmological simulations to a future study.  \nIn general, three key components are required for modeling the binary populations responsible for the GWB: (i) galaxy masses and merger rates (§ 3.2.1), (ii) SMBH masses based on a galaxy-host relationship (§ 3.2.2), and (iii) a binary evolution prescription (§ 3.2.3). We choose particular parameterizations for each of these components following Chen et al. (2019), described below, which are implemented in the holodeck code. A large number of free parameters are required for any such type of population synthesis calculation-more than can be meaningfully fit by the existing data. We therefore identify key parameters to vary in the models considered here and adopt standard literature values for the rest. In our analysis, we use these models to construct numerous di ff erent libraries of binary populations and we explore the impact of varying these key parameters on the resulting GWB spectra.  \nFor the binary evolution, we consider libraries using both a phenomenological binary inspiral model (dubbed Phenom ) and a naive GW-only inspiral scenario ( GWOnly ). In both cases, the libraries vary two parameters that determine galaxy number density (the normalization ψ 0 and turnover mass m ψ, 0) and two parameters describing the M BHM bulge relationship (the normalization µ and intrinsic scatter ϵµ ). The Phenom library includes two additional parameters describing the total binary lifetimes τ f , and the binary hardening rate at small separations ν inner. These models and parameters are described in detail in the following sections. Table B1 lists all of the parameters, giving fiducial values when they are fixed and the prior distributions for those which are varied. Table B2 summarizes our di ff erent libraries and which parameters are varied in each. In Appendix C, we also com-  \nainst larger, \'extended\' models in which additional parameters are varied ( Phenom-Ext & GWOnly-Ext ).  \n3.2.1. The Galaxy Merger Rate  \nThe number density of galaxy mergers ( η gal-gal) can be expressed (Chen et al. 2019) in terms of a galaxy stellar-mass function (GSMF; Ψ ), a galaxy pair fraction (GPF; P ), and a galaxy merger time (GMT; T gal-gal):  \n∂ 3 η gal-gal ∂ m ⋆ 1 ∂ q ⋆ ∂ z = Ψ ( m ⋆ 1 , z \' ) m ⋆ 1 ln(10) P ( m ⋆ 1 , q ⋆, z \' ) T gal-gal( m ⋆ 1 , q ⋆, z \' ) ∂ t ∂ z \' . (13)  \nThis distribution is calculated in terms of the stellar mass of the primary galaxy m ⋆ 1, the stellar mass ratio ( q ⋆ = m ⋆ 2 / m ⋆ 1 ≤ 1), and the redshift z . Because the galaxy merger spans a finite timescale ( T gal-gal) and corresponding redshift interval, we distinguish between the initial redshift at which a galaxy pair forms ( z \' = z \' [ t ] at some initial time t ) and the redshift at which the system becomes a post-merger galaxy remnant ( z = z [ t + T gal-gal]). An additional delay is required for binaries to reach the PTA frequency band, which is characterized in § 3.2.3.  \nThe GSMF is defined as  \nΨ ( m ⋆ 1 , z \' ) ≡ ∂η⋆ ( m ⋆ 1 , z \' ) ∂ log 10 m ⋆ 1 , (14)  \ni.e., the di ff erential number-density of galaxies per decade of stellar mass. The implementation used in this analysis described the GSMF in terms of a single Schechter function (Schechter 1976):  \nΨ ( m ⋆ 1 , z ) = ln (10) Ψ 0 · " m ⋆ 1 M ψ # αψ exp -m ⋆ 1 M ψ ! , (15)  \nwhere we have introduced Ψ 0, M ψ , and αψ as new variables. In order to allow the GSMF to vary with redshift, we parameterize these quantities as  \nlog 10 GLYPH<16> Ψ 0 / Mpc -3 GLYPH<17> = ψ 0 + ψ z · z , log 10 GLYPH<16> M ψ/ M ⊙ GLYPH<17> = m ψ, 0 + m ψ, z · z , αψ = 1 + αψ, 0 + αψ, z · z , (16)  \nsuch that each of these quantities has a simple linear scaling with redshift. This introduces six new dimensionless parameters into our models, corresponding to the normalization ( ψ 0, m ψ, 0, & αψ, 0) and slope ( ψ z , m ψ, z , & αψ, z ) of the redshift scaling. In all of the analysis presented here, the latter three are always kept fixed at the fiducial values specified in Table B1. The GSMF normalization and characteristic mass parameters ψ 0 and m ψ, 0 are allowed to vary in our fiducial Phenom library, while αψ, 0 is additionally varied in Phenom-Ext .  \nThe GPF and GMT are defined as  \nP ( m ⋆ 1 , q ⋆, z \' ) ≡ ∂ ∂ q ⋆ N ⋆, pairs( m ⋆ 1 , q ⋆, z \' ) N ⋆ ( m ⋆ 1 , z \' ) ! , (17)  \nT gal-gal( m ⋆ 1 , q ⋆, z \' ) ≡ Z a ⋆, f a ⋆, i GLYPH<2> ˙ a ⋆ ( m ⋆ 1 , q ⋆, z \' ) GLYPH<3> -1 da ⋆, (18)  \nwhere ˙ a ⋆ denotes the rate at which the merging galaxies\' separation decreases. The GPF describes the number of observable galaxy pairs relative to the number of all galaxies. The GMT is the duration over which two galaxies can be discernible as pairs, from an initial separation a ⋆, i at which they are associated with one another, until a final separation a ⋆, f after which they are no longer distinguishable as separate galaxies. These two distributions are typically determined empirically based on the detection of galaxy pairs in observational surveys and thus depend on observational definitions and selection criteria (e.g., Conselice et al. 2008; Mundy et al. 2017; Snyder et al. 2017; Duncan et al. 2019).  \nIn practice, we parameterize P ( m ⋆ 1 , q ⋆, z \' ) and T gal-gal( m ⋆ 1 , q ⋆, z \' ) as redshift-dependent power laws of m ⋆ 1, q ⋆ and z following Chen et al. (2019):  \nP ( m ⋆ 1 , q ⋆, z \' ) = P 0 m ⋆ 1 10 11 M ⊙ ! α p (1 + z ) β p q γ p α p = α p , 0 + α p , z · z γ p = γ p , 0 + γ p , z · z . (19)  \nT gal-gal( m ⋆ 1 , q ⋆, z \' ) = T 0 m ⋆ 1 10 11 M ⊙ / h ! α t (1 + z ) β t q γ t α t = α t , 0 + α t , z · z γ t = γ t , 0 + γ t , z · z . (20)  \nAs shown in Table B1, the values of the corresponding parameters are kept fixed to standard literature values in our fiducial model ( Phenom ), but parameters governing the scaling of the GPF and GMT with redshift ( β p , 0 & β t , 0), the scaling of the GPF with mass ratio ( γ p , 0), and the GMT normalization ( T 0) are allowed to vary in our extended models, as described in more detail below.', '3.2.2. The SMBH-Host Relation': 'The number density of galaxy mergers given in Equation (13) is a distribution that describes the number of galaxy pairs as a function of galaxy properties. We assume a oneto-one correspondence between galaxy pairs and SMBH binaries and adopt a SMBH-host relationship to translate from galaxies to SMBHs. In this analysis, we restrict ourselves to the M BHM bulge relationship, which relates the galaxy stellar bulge mass to the SMBH mass for each component of the binary as (Marconi & Hunt 2003)  \nlog 10 ( M BH / M ⊙ ) = µ + αµ log 10 M bulge 10 11 M ⊙ ! + N GLYPH<16> 0 , ϵµ GLYPH<17> . (21)  \nHere N (0 , ϵµ ) denotes normally distributed random scatter with a mean of zero and standard deviation of ϵµ (in dex). This relation depends on three model parameters that are allowed to vary in our analyses: the dimensionless BH mass normalization ( µ ), the intrinsic scatter ( ϵµ , in dex), and the  \npower-law index αµ (which is varied only in our extended models, and is dimensionless). A fraction of the galaxy stellar mass is in the stellar bulge component ( M bulge = f ⋆, bulge · m ⋆ ), which we take to be f ⋆, bulge = 0 . 615 based on empirical bulge fraction measurements of massive galaxies from Lang et al. (2014) and Bluck et al. (2014). Using the M BHM bulge relationship, we transform the number density of galaxy mergers to a number density of SMBH binaries via  \n∂ 3 η ∂ M ∂ q ∂ z = ∂ 3 η gal-gal ∂ m ⋆ 1 ∂ q ⋆ ∂ z ∂ m ⋆ 1 ∂ M ∂ q ⋆ ∂ q . (22)  \nEquation (22) provides an expectation value for the number of binaries in a point ( M , q , z ) in parameter space. To discretize the SMBH binary population, and also measure the e ff ects of cosmic variance, we assume that the true number of binaries in any given spatial volume is Poisson-distributed. We then integrate the di ff erential number of binaries over finite bins of parameter space to obtain the expected number of binaries in each bin. We generate multiple realizations by drawing many times from a Poisson distribution ( P ) centered at that value. Finally, we sum over parameter-space bins to calculate the resulting GWB spectrum. In practice, we implement Equation (5) as:  \nh 2 c ( f ) = X M , q , z , f P ∂ 4 N ∂ M ∂ q ∂ z ∂ ln fp ∆ M ∆ q ∆ z ∆ ln f ! h 2 s ( fp ) ∆ ln f . (23)', '3.2.3. Binary Evolution': "The final component for constructing the GWB is the most uncertain: the binary evolution from the initial galaxy merger until the eventual SMBH coalescence. Typically, interactions with the astrophysical environment (i.e., stars and gas in the host galaxy) are required to bring an SMBH binary into the PTA band within a Hubble time. For example, the high-mass binaries ( M ≳ 10 9 M ⊙ ) that dominate the GW signals in the PTA band must reach separations of ∼ 0 . 1 pc before GW emission becomes dominant and drives e ffi cient inspiral. Those binaries enter the NANOGrav band (currently: ∼ 1 / 15 yr -1 ≈ 2 nHz) when they reach separations of ≈ 0 . 05 pc - only a factor of two smaller. This immediately implies that the environmental processes may play a non-negligible role in binary evolution, even after the binaries reach the NANOGrav band.  \nEven if environmental hardening is e ff ective in bringing binaries to the PTA-detectable frequencies, binary lifetimes can still be many billions of years, and a large fraction of binaries may stall (Kelley et al. 2017a). This can lead to binaries reaching the PTA band at substantially lower redshifts than those at which their respective galaxy mergers occurred.  \nDetailed modeling of environmentally driven binary evolution can introduce dozens of free parameters, even when  \nSMBH and galaxy parameters are known a priori . Ultimately, many of these parameters become significantly degenerate in determining the resulting shape of the GWB spectrum and the properties of the SMBH binaries producing it. For this reason, we focus this analysis on a 'phenomenological' model that is designed to capture the overall e ff ects of more explicit binary evolution, while introducing only a small number of free parameters. In these models, the hardening timescale is parameterized in terms of the evolution of the binary semi-major axis a as  \nda dt GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> phenom = Ha · a ac ! 1 -ν inner · 1 + a ac ! ν inner -ν outer . (24)  \nThe hardening timescale is thus a double power law, with a break at the critical separation ac , and asymptotic behaviors of:  \ndt d ln a ( a ≪ ac ) ∼ a ν inner , (25)  \nin the 'inner' (small-separation) regime, and  \ndt d ln a ( a ≫ ac ) ∼ a ν outer , (26)  \nin the 'outer' (large-separation) regime. Hardening rates are added linearly, such that the total rate of evolution when also including GW emission (Equation 8) is given by da / dt = [ da / dt ] phenom + [ da / dt ] GW . We assume a fixed value of ν outer = + 2 . 5 in all of our analysis, motivated by detailed literature models of dynamical-friction-driven evolution of SMBH binaries (Kelley et al. 2017a). ν inner, which controls the hardening rate of binaries as they approach and enter the PTA band, is allowed to vary in our models.  \nIn addition to the two power-law indices ( ν inner , ν outer), and the characteristic break separation ( ac ), the normalization ( Ha ) is calculated such that the total lifetime of the binary matches a target τ f , i.e.,  \nτ f = Z a isco a init da dt ! -1 da , (27)  \nwhere a init is the initial binary separation and a isco ≡ 6 GM / c 2 is the innermost stable circular orbit, where we consider the two SMBHs to have merged. While this expression for a isco is based on the test-particle approximation ( q ≪ 1), the true value should di ff er by less than a factor of two (Flanagan & Hughes 1998) for low SMBH spins, and the contribution to the total lifetime is always negligible for a ∼ a isco. The total lifetime τ f is a key parameter that we vary in our models.  \nAt numerous po14ints in our analysis we compare the self-consistent, phenomenological model (in the Phenom and Phenom-Ext libraries) against a model where binaries decay only due to GW emission ( GWOnly and GWOnly-Ext libraries). In the GW-only model, we take the redshift (and  \n]  \nFigure 3. SMBH binary hardening timescales and GW spectra are shown for varying values of phenomenological binary evolution model parameters. The top panel shows the hardening timescale ( τ = dt / d ln a ), with the black, horizontal lines corresponding to total binary lifetimes. Solid lines correspond to binary lifetimes τ f = 0 . 1 Gyr, while dashed are 1 . 0 Gyr. The small separation hardening rate power-law index, ν inner, is also varied which changes the amount of time binaries spend at sub-parsec separations. The vertical dotted line shows the separation at which an M = 10 9 M ⊙ system reaches f = 1 / 16 . 03 yr. The bottom panel shows the resulting GWB characteristic strain. The dotted line shows GW-only driven evolution. For clarity, the top panel includes only binaries with 3 × 10 8 M ⊙ < M ≤ 3 × 10 9 M ⊙ and 0 . 1 < q ≤ 1 . 0. In both panels, the shaded regions denote the inter-quartile range. Other parameters of these population are fiducial values for the Phenom set of models (§ 3.3). Variation in the binary evolution parameters significantly impacts the shape and low-frequency amplitude of the GWB spectrum.  \n<!-- image -->  \nthus source distance) to be the post-galaxy-merger redshift, without an additional delay, and set the binary evolution time in Equation (7) to be that of GW-only evolution (i.e., Equation 8). This model is not self-consistent, as GW-only evolution is unable to bring binaries to the PTA band within a Hubble time. It is nonetheless a useful comparison, because the GW-only assumption is often still used in the literature and tends to produce the highest GWB amplitudes.  \nFigure 3 shows the binary evolution and GWB spectra resulting from the phenomenological evolution model. Total binary lifetimes of 0.1 and 1 Gyr are plotted with solid and dashed lines respectively, while varying small-separation power-law indices ( ν inner) are shown with di ff erent colors. In each panel, the median and 50% interquartile range of binaries are shown. Note that in the top panel, only binaries with 3 × 10 8 M ⊙ < M ≤ 3 × 10 9 M ⊙ and 0 . 1 < q ≤ 1 . 0 are shown. In the environmentally driven regime (larger separations), their hardening rate is determined such that their total life-  \ns the target value. The narrow interquartile regions in the environmental regime reflect the small variations in hardening rate required to produce the target total lifetime, for this range of masses.  \nThe GW hardening rate, which dominates at small separations, is determined entirely by the binary masses for our assumption of circular orbits. In the phenomenological model, the hardening rate at larger separations is determined such that the total inspiral time matches the input binary lifetime. This means that shorter lifetime populations are forced to transition into the GW-driven regime at smaller separations. The power-law indices also a ff ect the transition point by determining which separations dominate the binary's evolution time. More positive values of ν inner lead to flatter evolution trajectories with less and less time spent at sub-parsec binary separations. For reference, the dotted vertical line shows the separation at which a binary with M = 10 9 M ⊙ enters the 15 yr NANOGrav frequency band. The two models with ν inner = -0 . 5 and ν inner = -1 . 0 lead to environmentally driven evolution until su ffi ciently small separations so that the resulting GWB spectral turnover ( bottom panel) is clearly visible in simulated 15 yr spectra.  \nWe note that a value of ν inner = -1 . 0 is well-motivated by numerical stellar scattering experiments of closely bound SMBH binaries (Sesana 2010; Sesana & Khan 2015). However, the true rate of environmental hardening for close binaries will depend on the stellar distribution in a given host galaxy, as well as on the role of gas-driven binary evolution, motivating the choice to allow ν inner to vary in our models.  \nThe ν inner = -0 . 5 variation produces a substantial attenuation of the GWB: up to a 50% decrease in characteristic strain (75% reduction in GW power). Even though the 1 . 0 Gyr lifetime model (dashed lines) qualifies as e ffi cient and rapid binary evolution, the overall amplitude of the GWB at all frequencies is ∼ 10% lower than the GW-only model. This is because a fraction of the binaries, specifically those which formed within a look-back time of 1 . 0 Gyr, are unable to reach the PTA band before redshift zero, and thus do not contribute to the observable GWB. This figure highlights that only in a narrow region of parameter space do realistic GWB spectra match predictions from GW-only models, but the di ff erences between the self-consistent and GW-only models can be subtle.", '3.3. Libraries of SMBH Binary Populations and GWB Spectra': 'With the models described above, we use the holodeck code to calculate libraries of SMBH binary populations and their resulting GWB spectra. In each parameter space that we explore, we include a large number of sample points in the space in addition to many realizations of populations and spectra at each point. We use the same GW frequency bins as  \nthe 15 yr NANOGrav data ( fi = i / 16 . 03 yr -1 = i × 1 . 98 nHz, see § 2) to calculate spectra. Here we present the two primary libraries used in our analysis: Phenom and GWOnly , and outline some of their features. In Table B2 we summarize the parameters that are varied in each library, while the full list of model parameters (including fiducial values for fixed parameters and assumed prior distributions for varied parameters) is given in Table B1.  \nTens of free parameters are required in these models, many of which are poorly constrained either observationally or theoretically. In addition, many of them are formally degenerate in their e ff ects on the resulting GWB spectra (e.g., Chen et al. 2019). For this reason, we adopt as our fiducial library, Phenom , a model with six parameters { ψ 0 , m ψ, 0 , µ, ϵµ, τ f , ν inner } that produce SMBH binary populations and GWB spectra that e ff ectively span the broader model uncertainties. More specifically, this library varies the normalization and turnover mass of the GSMF ( ψ 0 and m ψ, 0,), along with the normalization and scatter of the M BHM bulge relationship ( µ and ϵµ ). The library also utilizes the phenomenological binary evolution model and varies the SMBH binary lifetime, τ f , and the hardening power-law index at small separations ν inner.  \nThe di ff erences in GWB spectra for systematic variations in Phenom model parameters are shown in Figure 4. The overall amplitude of the GWB spectrum varies most significantly in the left and top-middle panels, indicating that the GWB amplitude is most sensitive to parameters determining SMBH masses ( m ψ, 0, µ ) and the SMBH binary number density ( ψ 0). In the lower-middle panel, we see that increasing scatter in the M BHM bulge relationship ( ϵµ ) also increases the GWB amplitude. This owes to the fact that larger scatter increases the e ff ective SMBH masses through Eddington bias: because low-mass SMBHs are more numerous, their scatter towards higher masses outnumbers the scatter of the rarer, higher-mass SMBHs towards lower values. Notice that variations in SMBH mass, GSMF turnover mass, and M BHM bulge scatter (parameterized by µ , m ψ, 0, and ϵµ , respectively) all produce qualitatively similar changes in the GWB spectra. Higher masses preferentially increase the low-frequency amplitudes, thereby steepening the spectra at higher frequencies ( f ≳ 1yr -1 ). This occurs as rare high-mass binaries contribute less to the GWB at higher frequencies, due to their rapid evolution at smaller separations (see § 3.2). The ϵµ parameter shows this frequency-dependent e ff ect even more prominently, as it preferentially a ff ects the highest SMBH mass bins where the gradient in SMBH number density with respect to mass is steepest.  \nThe shape of the spectrum at low frequencies is determined by the binary hardening rate da / dt (as introduced in § 3.2.3), which includes the interaction of binaries with their nuclear galactic environments. Recall that in our models, the binary lifetime is an input parameter, one that is varied in the top-  \nFigure 4. Comparison of GWB spectra for systematic variations in the parameters of our fiducial range of models ( Phenom ). In each panel, the indicated parameter is varied while the other parameters are fixed at typical values. The purple dot-dashed line shows the spectrum from a fiducial GW-only model, while the black dotted line shows a pure f -2 / 3 power law as a reference. Variations in the GSMF parameters ( ψ 0-the dimensionless GSMF normalization, and m ψ, 0-the dimensionless GSMF turnover mass) are shown in green. Variations in the M BHM bulge parameters ( µ -the dimensionless M BHM bulge normalization), and ϵµ -the M BHM bulge scatter in dex) are in orange. Variations in the SMBH binary lifetime ( τ f in units of Gyr) and hardening power-law index ( ν inner) are shown in blue. The numeric label indicates the value of the parameter for that particular curve, and the fixed values for all other parameters are: ψ 0 = -2 . 5, m ψ, 0 = 11 . 5, µ = 8 . 25, ϵµ = 0 . 5 dex, τ f = 1 . 0 Gyr, and ν inner = -1 . 0. The solid lines show the median for each parameter value from 10,000 realizations. The shaded regions indicate the 16 th and 84 th percentiles of the distribution. For clarity, we only plot the shaded regions for every other parameter value. Apart from the mostly degenerate mass parameters µ and m ψ, 0, which nonetheless have a significant influence on the GWB spectrum, each of these parameters impacts the GWB spectra in distinct ways. This indicates the promise of GW observations for constraining the SMBH binary population.  \n<!-- image -->  \night panel of Figure 4 and kept fixed at τ f = 1 . 0 Gyr in all other panels. Consequently, for a given τ f , binaries of different masses enter the GW regime at di ff erent frequencies. Some of the variations in low-frequency spectral shape seen when the mass-determining parameters ( µ , ϵµ , & m ψ, 0) are varied and hardening parameters are kept fixed can therefore be attributed to τ f being the same for all binary masses. For our models with fixed total binary lifetimes, populations with lower masses tend to have stronger low-frequency turnovers as lower-mass binaries enter the GW regime at higher frequencies. Equivalently, lower mass systems spend more time at higher frequencies, meaning that their environmentally driven evolution must have proceeded even faster at lower frequencies.  \nWhile there is some degeneracy across all parameters, only the two parameters that directly a ff ect the average binary mass ( m ψ, 0 and µ ) produce mostly degenerate spectral  \nchanges. As mentioned above, even the M BHM bulge scatter parameter ( ϵµ ), which also changes the average binary mass, is noticeably distinct. This speaks to the possibility of independently constraining multiple parameters with a su ffi -ciently high signal-to-noise ratio, even without appealing to additional information content such as sky anisotropy, individual continuous-wave sources, or electromagnetic counterparts and other multi-messenger constraints.  \nAs introduced in § 3.2.3, in addition to Phenom , we also use a library with the same variations in GSMF and M BHM bulge parameters, but using GW-only evolution instead of the phenomenological model. We refer to this fourdimensional parameter space { ψ 0 , m ψ, 0 , µ, ϵµ } as GWOnly . Figure 5 shows a comparison of GWB spectra with variations in the Phenom model parameters versus spectra from the GWOnly models. For longer binary lifetimes ( ≫ Gyr), fewer systems are able to coalesce, and the GWB amplitude  \nFigure 5. Comparison of GWB spectra for systematic variations in two of the six model parameters ( τ f and ν inner) normalized to a model in which binaries evolve only because of GW emission. Colors and line styles are as in Fig. 4. The phenomenological hardening models show substantial deviations from the GW-only hardening models, especially for parameters far from our fiducial values of τ f = 1 Gyr and ν inner = -1. In particular, the spectra are suppressed at the lowest frequencies for larger values of ν inner and produce a more pronounced turnover.  \n<!-- image -->  \nis noticeably diminished at all frequencies. For shorter lifetimes ( ≲ Gyr), non-GW hardening is still important at low frequencies within the 15 yr NANOGrav band. This leads to binaries evolving faster than the GW-only prediction, fewer binaries existing at these frequencies, and thus attenuated GW emission producing a low-frequency turnover (Kocsis & Sesana 2011; Ravi et al. 2014). Moderate inspiral times ( ∼ Gyr) produce the closest match between the phenomenological and GW-only models, but still show a slight turnover in addition to an amplitude ∼ 10% lower at all frequencies.', '3.4. Interpolation of Population Synthesis Models with Gaussian Processes': "In order to infer the properties of SMBH binary populations that are consistent with the GWB, we need to compare the theoretically expected GWB spectra from holodeck with the observed NANOGrav data. Previous work used analytic expressions for this, e.g., by fitting the GWB spectra with a single power law (for a population of circular binaries purely driven by GWs) or a broken power law (to capture the turnover produced by environmental interactions; Sampson et al. 2015). However, the properties of the SMBH binary population are only indirectly extracted from these fits, and disentangling potential covariances between population parameters is challenging. To overcome this limitation, Taylor et al. (2017) developed a modeling framework that directly links the properties of the GWB spectrum to the binary population parameters by training Gaussian processes (GPs) on simulated GWB spectra from population-synthesis models. Here we adopt this approach to interpolate the strain of the GWBacross simulated holodeck libraries, generated in dis-  \ncrete points of the binary parameter space, to accurately predict the GWB spectrum at any point in the space.  \nGPs provide a powerful interpolation method that parameterizes noisy data in terms of a multivariate Gaussian distribution with a mean vector and covariance function (see Aigrain & Foreman-Mackey 2022, for a review). The covariance functions can be custom built from a suite of versatile kernel functions allowing for quick adaptability to a variety of complex parameter spaces. While GPs are not sparse and lose e ffi ciency in high-dimensional spaces (e.g., greater than a few dozen), one key advantage of GP regression is that it provides an estimate of the uncertainty in the interpolation process (i.e., the prediction is probabilistic). Importantly, one can use this in an iterative process to adapt and improve the fitting. Additionally, the GP uncertainty can be propagated forward to our final statistics, allowing for a full marginalization over the interpolation uncertainties.  \nThe GPs are trained on holodeck GWB spectra using the George GP regression library (Ambikasaran et al. 2015a), as in Taylor et al. (2017) and Arzoumanian et al. (2018). To capture fluctuations that arise directly from the discrete nature of the binary population, we train the GPs at each sampling frequency of the GW spectrum, fi . Since GP regression assumes that the interpolated quantity (here the strain of the GWB) is smooth with respect to the interpolation variables (here the model parameters). The use of an independent GP at each frequency thus enforces smoothness in the GW spectrum across model parameter space at a given frequency, but not across frequencies. Because the binary population is independent at each frequency, smoothness across frequencies is not expected in general. Two separate GPs are trained per frequency, one on the median values of log 10 GLYPH<16> h 2 c ( fi ) GLYPH<17> , and one on its standard deviation. This allows us to predict both the typical value and the typical spread of the strain, and to account for the uncertainty in each value's interpolation separately.  \nWeselect the training set (i.e., the library generation points that make up our model grid) from our multi-dimensional parameter space using Latin hypercube (LHC) sampling (e.g., see Taylor & Gerosa 2018, and references therein). 5 This o ff ers an e ffi cient method to generate a near-random set of parameter values, representative of the entire parameter space, with relatively few points. Since we aim to explore high-dimensional spaces, this type of sampling is necessary to keep the total number of simulations computationally tractable.  \nFigure 6. Accuracy of trained Gaussian process (GP) interpolants for our fiducial, six-dimensional parameter space Phenom , aggregated over the first 5 frequency bins used in our primary analyses. The top panel shows the error ( ≡ predicted -true) on the median log 10 h 2 c | predicted (solid lines) for the training set (gray) and validation set (green), while the middle panel shows the error on the standard deviation (dashed lines). Contours contain 20%, 50%, and 90% of the population. The bottom panel shows the distributions of errors for both validation sets. The vertical lines bound the 50% and 90% region of the errors. The GPs predict 99.4% (98.5%) of the validation set for medians (standard deviations) within 10% of the actual value. Both the training and validation sets contain 2,000 sample points, with medians and standard deviations calculated over 2,000 realizations, which we have found to be more than su ffi cient in accurately training the GPs.  \n<!-- image -->  \nThe training of the GPs proceeds as follows. Using the LHC method, we draw s samples in the binary parameter space. For each sample, we produce r realizations of the GWB spectrum using holodeck and calculate the median and standard deviation of log 10 ( h 2 c ) at each frequency fi . These means and standard deviations constitute the inputs for the training of the two GPs. For each point in the training set, GPs require the value of the quantity on which they are trained (here the median or standard deviation of log 10 [ h 2 c ( fi )]) and optionally its uncertainty. Including uncertainties on the input values helps to avoid overfitting. We adopt the standard uncertainties for sample mean and sample standard deviation. When training on the median, we estimate the uncertainty as the standard deviation divided by r 1 / 2 ; and for training on the standard deviation, the uncertainty is given by the standard deviation divided by [2( r -1)] 1 / 2 .  \nTo test the performance of the GPs, we create a validation set with points in the SMBH binary parameter space that were not included in the training set. For each validation point, we calculate the median and standard deviation of log 10 ( h 2 c ), both with GPs and with holodeck simulations. For comparison purposes, we label the value obtained from GPs as the 'predicted' value, while the holodeck values are considered to be the 'true' value. Based on this, an error (i.e., predicted minus true) can be calculated for the GP interpolation performance.  \nWe used this approach to test a variety of kernels (i.e., covariance functions) along di ff erent directions in parameter space to determine which combination most accurately captured each parameter's response to changing the GWB. We determined that two types of kernels were necessary: rational quadratic kernels for the phenomenological timescale τ f and the hardening power-law index ν inner, and squared exponential kernels for the remaining parameters. An iterative process of checking the performance of the GPs was used to determine the necessary number of LHC sample points and holodeck realizations to converge on a su ffi cient accuracy level. The performance of the GP trained on the median values is more sensitive to the choice of the number of sample points, while the performance of the GP trained on the standard deviations is more sensitive to the choice of the number of realizations. We found that training on s = 2,000 LHC samples with r = 2,000 holodeck realizations at each sample point was more than enough to acquire the desired accuracy level. Figure 6 shows the response of GPs trained on the Phenom library for the 5 frequencies used in our analysis. The reconstruction is quite accurate, with 99.4% (98.5%) of the test set cases for medians (standard deviations) falling within 10% of the actual value-significantly smaller than the standard deviation across spectra realizations.", '3.5. Fitting Simulated GWB Spectra to PTA Observations': "In Arzoumanian et al. (2018), once GPs were trained, they were inserted into the full PTA likelihood calculation in order to obtain posteriors on the SMBH binary population parameters. As PTA data sets have grown in size and with the new discovery of HD correlations, the likelihood computation time has increased. As such, inserting two GPs into the 15 yr data set's likelihood calculation in order to obtain posteriors for Phenom was not an e ffi cient analysis approach.  \nInstead, we use the ceffyl package (Lamb et al. 2023) to fit the interpolated GWB spectra to the previously computed free-spectrum posteriors of the cross-correlated timing-residual PSD. Fitting on intermediate PTA analysis products, such as the free-spectrum posteriors, o ff ers a substantial speed-up by factors of 10 2 -10 4 compared to directly fitting the full likelihood of timing residuals. Importantly, the resulting posterior distributions of GW spectral model pa-  \nrameters achieved by ceffyl have been found to be nearly identical to those obtained from the full likelihood approach.  \nIn detail, we expand the likelihood, L GLYPH<16> ⃗ d | ⃗ Θ GLYPH<17> , where ⃗ d is the PTA data (e.g., the TOAs) and ⃗ Θ are the SMBH binary population parameters (e.g., the parameters from Phenom ), by inserting an intermediate data product such as the freespectrum posteriors (log 10 ( ρ i )). Then, instead of directly calculating the fit of a GWB spectrum (generated by the trained GPs for a given draw of SMBH binary population parameters) to the TOAs, we compute the probability that a given GWBspectrum is supported by the free-spectrum posteriors. The expanded likelihood function is now given by  \nL GLYPH<16> ⃗ d | ⃗ Θ GLYPH<17> ∝ Nf Y i = 1 Z d GLYPH<0> log 10 ρ i GLYPH<1> p GLYPH<16> log 10 ρ i | ⃗ d GLYPH<17> p GLYPH<16> log 10 ρ i | ⃗ Θ GLYPH<17> , (28)  \nwhere Nf is the number of Fourier components used in the GWB analysis (5 or 14, but see also § 2 and NG15gwb), p GLYPH<16> log 10 ρ i | ⃗ d GLYPH<17> is the posterior probability density of log 10 ( ρ i ) (i.e., the free-spectral posteriors) which are represented by highly optimized kernel density estimators, while p GLYPH<16> log 10 ρ i | ⃗ Θ GLYPH<17> is the probability of log 10 ρ i given a GWB spectrum from the trained GPs. Since the GPs are trained on the median and standard deviation of the characteristic strain log 10 ( h 2 c ), these provide the mean and variance of a Gaussian when calculating p GLYPH<16> log 10 ρ i | ⃗ Θ GLYPH<17> . The above likelihood is sampled through MCMC techniques to obtain the resultant posteriors on ⃗ Θ .  \nWhile all of the libraries generated for GP training draw uniformly from the SMBH binary population parameter space, when we perform the MCMC analysis, we have the opportunity to place di ff erent priors onto each parameter. For the analysis in this paper, we utilize two distinct prior setups: a uniform prior and a set of astrophysical priors based on galaxy observations (e.g., see Table B1). When relevant, we denote the prior distribution shape in combination with the library designation as e.g., Phenom + Uniform or Phenom + Astro (see Table B2).", '4. RESULTS': 'We simulate populations of SMBH binaries using a phenomenological ( Phenom ) and GW-only ( GWOnly ) model. We create holodeck libraries of GWB spectra at fixed points of the SMBH binary parameter space and interpolate them with GPs. We fit the models to the 15 yr freespectrum posteriors considering the HD-w / MP + DP + CURN as the fiducial 15 yr NANOGrav results for this analysis (but we also fit the HD-DMGP posteriors for comparison) using both uniform and astrophysically motivated priors (see Table B1). As shown in Table B2, the Phenom library is fit against the data using both uniform priors and as-  \nFigure 7. Power-law amplitude ( A ) and spectral index ( γ ) from purely power-law fits to HD free-spectrum model posteriors from the 15 yr data set, compared to simulated GWB spectra from holodeck libraries. Data set fits include both the 15 yr HDw / MP + DP + CURN and HD-DMGP models for comparison. Fits to the lowest 5-frequency bins of spectra from two holodeck libraries are shown: the self-consistent phenomenological binary evolution model ( Phenom ) and the purely GW-driven evolution model ( GWOnly ). We show 1, 2, & 3 σ contours for each. The analytic, GWB PSD power-law index of γ = 13 / 3 is shown as reference (black, dashed). The spectral shape of the HD signal present in the 15 yr data set is broadly consistent with expectations for a GWB from binary SMBH population. The amplitude is towards the higher end of predictions and the recovered spectral index deviates from the idealized power law in similar ways as the phenomenological binary evolution model.  \n<!-- image -->  \ntrophysically informed priors ( Phenom + Uniform and Phenom + Astro ), while the GWOnly library is fit only with uniform priors ( GWOnly + Uniform ). Our results are summarized as follows.  \nIn all of our analysis, we find that the NANOGrav 15 yr data set is consistent with a GWB produced by a population of SMBH binaries. In the first, most simplified approach, power-law fits 6 to both the observed GWB spectrum and those from simulations produce amplitudes and spectral indices that overlap in the 2- and 3σ regions depending on model (§ 4.1). The remainder of this section presents the results of our systematic approach of fitting simulated SMBH binary populations to the data, which yield more realistic GWBspectra that match the 15 yr results (§ 4.2). From these fits we obtain posterior distributions on uncertain astrophys-  \nl parameters of the SMBH population synthesis models (§ 4.3) and make predictions from our models for the properties of the population of SMBH binaries that produce the GWobservations (§ 4.4).', '4.1. Comparison of idealized power-law fits to GWB spectra': 'The approach of fitting simple power-law models to the GWB is a common one in the literature. While idealized power-law fits to GWB spectra neglect most of the information imprinted by astrophysical processes on the background, they are e ff ective in broadly examining the consistency between simulated binary populations and PTA data sets. Therefore, we carry out this straightforward analysis as a first check of the Phenom and GWOnly libraries before implementing the full methodology described in § 3.4 - 3.5. In practice, we constrain the amplitude, A , and slope, γ , of an idealized power-law GWB spectrum (in timing-residual PSD; Equation 12) with a non-linear least-squares fits to the GWB spectra from each realization of the binary population from the Phenom and GWOnly libraries using the five lowest-frequency bins. We then compare these to the results of power-law fits of the 15 yr data, illustrating their overlap in the A -γ parameter space.  \nIn Figure 7, we show the range of GWB amplitudes and spectral indices (in timing-residual PSD) based on these fits. We see that the amplitudes and power-law indices vary significantly across the simulated GWB spectra. Even for the GWOnly models, which match the premise of the analytic γ = 13 / 3 ( α = 2 / 3) models, our simulations yield indices typically varying from 4-5 . 5 in the 95% credible region. Recall that even GW-only binary evolution with circular orbits does not produce a pure power-law spectrum, owing to the steepening of the spectrum at higher frequencies where the finite number of binaries in each frequency bin becomes important. The slight o ff set of the GWOnly models towards steeper values of γ > 13 / 3 reflects this higherfrequency spectral steepening caused by finite-number effects. The Phenom libraries, which self-consistently model the e ff ects of environmental interactions on binary evolution, produce much wider ranges of spectral indices as disparate as 2-7, with lower values corresponding to shallower characteristic strain spectra (i.e., increasing across the lowest five frequency bins). Note that we exclude from Figure 7 a small number ( ∼ 1%) of Phenom samples in which all binaries stall and thus produce zero GWB.  \nWe also show the power-law parameter posteriors for the fiducial HD-w / MP + DP + CURN free-spectrum posteriors, and the HD-DMGP , which we use for comparison. While the 15 yr free-spectrum posteriors are not perfectly fit by powerlaw models (see Figure 1), the di ff erences between these two models highlight that the measurement of a spectral index is  \nparticularly sensitive to choices of fit in the 15 yr data, and to features in particular frequency bins (NG15gwb).  \nWe use the A -γ fits as a general measure of parameter space coverage. Figure 7 demonstrates that the range of simulated populations is able to reproduce the measured GWB within the two-sigma curve of the Phenom library, and between the two- and three- sigma curves of the GWOnly library.', '4.2. The GWB is Consistent with Expectations from Populations of SMBH Binaries': 'The consistency between the 15 yr NANOGrav data set and GWBs produced by SMBH binaries is best supported by an analysis of the full range of astrophysical information contained in the free-spectrum posteriors. Fitting the GPs trained on the GWOnly and Phenom libraries to the 15 yr data (with uniform and astrophysically motivated priors) facilitates a comparison of observations to the GWB spectra from SMBH binary populations that is agnostic to any particular spectral model (including a power law).  \nFigure 8 shows GWB spectra produced by our simulated SMBH binary populations that accurately fit the 15 yr HD-w / MP + DP + CURN free spectrum. As mentioned above, the Phenom library is fit both with uniform and astrophysically informed priors ( Phenom + Uniform and Phenom + Astro ), while the GWOnly library is fit only with uniform priors ( GWOnly + Uniform ). Thin curves show 200 random draws of the binary parameter posterior distributions for each of the above models, with thick lines denoting the maximum likelihood spectra for each model.  \nBoth libraries are able to fit the GWB within the 15 yr posteriors. However, the GWOnly spectra have more di ffi culty matching the data, as indicated by their preference for the edges of the 15 yr free-spectrum posteriors in the highest and lowest frequency bins, and the best fit spectrum missing the highest probability regions of the 15 yr GWB data. As a comparison, power-law fits are shown for the idealized γ = 13 / 3 ( α = 2 / 3) spectral indices obtained from analytic calculations of SMBH binaries (Phinney 2001). The GWOnly models, which more closely resemble these analytic estimates, tend to be steeper than the bulk of the 15 yr distributions. In contrast, the maximum-likelihood spectra and likelihood draws from the Phenom model exhibit noticeable spectral turnovers to match the 15 yr data. While these results are suggestive of a low-frequency turnover or flattened spectrum, they are still consistent with an α = 2 / 3 power law and the associated GW-driven evolution.', '4.3. Parametric Constraints on SMBH Binary Models': "The MCMC exploration of the likelihood in Equation (28) returns constraints (posterior distributions) on the parameters of the population synthesis models based on the observed  \nFigure 8. GWBspectra from simulated SMBH binary populations that best fit the 15 yr free-spectrum data. The left panel shows the square-root of cross-correlated timing-residual power ( ρ ) and the right panel shows characteristic strain ( hc ). GP interpolated spectra are shown, with thick lines showing maximum likelihoods and thin lines showing 200 random draws from the posteriors. The Phenom library uses self-consistent binary evolution models, while GWOnly assumes purely GW-driven evolution. For the former, fits using uniform priors ( Phenom + Uniform ) are compared against more informed, astrophysically motivated ones ( Phenom + Astro ). Power-law fits to the 15 yr spectra with A yr = 2 . 1 × 10 -15 are also shown as dotted black lines.  \n<!-- image -->  \nGWB spectrum. The peaks of the marginalized posteriors indicate the most likely values of the parameter space that the binary population must occupy in order to produce the 15 yr free-spectrum posteriors. Figure 9 shows the posteriors of these binary population parameters for the Phenom binary evolution model. Results are compared for di ff erent prior choices, i.e. the Phenom + Uniform and the Phenom + Astro fits. Owing to the substantial uncertainty in the GWB spectrum at NANOGrav's current sensitivity, the posteriors are sensitive to the assumed priors, and only weak parameter constraints can be made.  \nHowever, we can identify some general trends among the preferred parameter values. The measured amplitude of the GWBstrongly prefers a combination of e ffi cient mergers occurring in high-mass systems. The data favor short binary lifetimes ( τ f ), high GSMF number densities ( ψ 0), and high characteristic masses ( m ψ, 0, µ ). It is worth noting that the range of priors in the Phenom + Uniform fits is quite wide compared to typical values adopted in the astronomical literature (see references with Table B1). While the parameters in the Phenom + Astro model are more constrained, the posteriors are still fairly broad. Because our models utilize simplified analytic prescriptions for each physical component, we use broader parameter distributions in Phenom + Uniform in part to introduce some added flexibility. None the less, Figure 9 demonstrates that the posteriors almost uniformly favor parameters that produce larger GWB amplitudes (e.g., see also Figure 4). This suggests that the amplitude of the GWB inferred from the 15 yr data set is di ffi cult to reach with standard values of some astrophysical parameters.  \nVery long binary lifetimes are disfavored for both Phenom + Uniform and Phenom + Astro . A large fraction of binaries with such long lifetimes would fail to reach the NANOGrav frequency band, resulting in lower GWB amplitudes, inconsistent with the GW data. Flatter values of the hardening rate power-law index, ν inner ≳ -1 . 0, are also preferred, as they produce spectral turnovers in the lower frequency bins (see Figure 3) resembling what is seen in the 15 yr data. Steeper values of ν inner correspond to binaries that spend more time at ∼ 10 -2 - 1 pc separations and transition into the GW-dominated regime earlier, at frequencies below the PTA band. Flatter values of ν inner correspond to very e ffi cient inspiral through this range of separations, leading to environmentally driven evolution even in the lower PTA band which produces noticeable GWB attenuation and a sharp low-frequency spectral turnover (see Figure 5). We note that this is dependent on the parameterization of our binary evolution model and the parameters varied in the Phenom library. Steeper evolution profiles at very large separations ( ∼ kpc), i.e., larger values of ν outer could similarly produce low-frequency turnovers (but this is not explored in this paper, since ν outer is kept fixed throughout). In either case, e ffi cient binary inspiral in the environmentally driven regime produces the most noticeable spectral turnovers.  \nThe posteriors for GSMF and M BHM bulge parameters differ noticeably for uniform versus astrophysical priors, unlike the binary inspiral parameters. The posteriors for the normalization and the characteristic mass of the GSMF ( ψ 0 and m ψ, 0) favor values at the higher end of the prior range, especially in the Phenom + Uniform fits, in which the galaxy number densities are pushed against the edges of the prior. We  \nFigure 9. Binary evolution parameter posteriors from fitting against the 15 yr HD-w / MP + DP + CURN free spectrum. The parameters correspond to the phenomenological library Phenom : binary lifetime ( τ f ), the hardening power-law index ( ν inner), the GSMF normalization ( ψ 0) and characteristic mass ( m ψ, 0), and the M BHM bulge normalization ( µ ) and scatter ( ϵµ ). The fiducial holodeck spectra library is fit to the lowest five frequencies utilizing two separate priors. Phenom + Uniform uses flat priors across all parameter spaces (blue, solid-line contours), while Phenom + Astro (orange, dashed-line contours) uses constrained priors on the GSMF and M BHM bulge parameters that encapsulate observational measurements for these parameters. The contours in two dimensions correspond to the 1, 2, and 3 σ regions, and the shaded regions in the one-dimensional plots are 1 σ regions. Although individual parameters are only weakly constrained, the data strongly prefer e ffi cient mergers in high-mass systems.  \n<!-- image -->  \nalso see higher values of these posteriors when binary lifetimes are longer, such that larger fractions of binaries stall before reaching the PTA band.  \nUnsurprisingly, the GSMF characteristic mass ( m ψ, 0) is almost entirely degenerate with the M BHM bulge mass normalization ( µ ), as indicated by the diagonal band in the respective 2D posterior. The scatter in the M BHM bulge relationship ( ϵµ ),  \nhowever, shows di ff erent trends. This is likely due to two factors. First, increasing ϵµ primarily increases the GWB amplitude in the lower frequency bins (Figure 4), as larger scatter preferentially increases SMBH masses and higher masses are more prevalent at lower frequencies (discussed more below). Secondly, larger values of ϵµ also produce significant variance across multiple population realizations which may decrease the aggregated likelihoods when calculating fits.  \nIn Figure 10, we compare the one-dimensional distributions of parameter priors versus posteriors for the Phenom + Uniform , Phenom + Astro , and GWOnly + Uniform models. Particularly in the case of the M BHM bulge parameters, we see that the posteriors closely follow the priors. While still consistent with the priors, the GSMF parameters are pushed noticeably towards higher values, even for the Phenom + Astro fits.  \nFigure 10 also shows fits using the GWOnly library. Note that this library does not include the τ f or ν inner parameters by definition. The posterior distributions for the GSMF and M BHM bulge parameters are generally consistent in both the Phenom and GWOnly models, with only weak constraints. While GWOnly shows the same preference for high values of ψ 0 as the phenomenological model, the preference is less pronounced. This is likely due to the decrease in GWB power at the lowest frequencies when a spectral turnover is induced, which is consistent with the covariance seen in Figure 9 between binary lifetime and GSMF normalization.  \nIt is important to note that these parametric constraints must be interpreted in the context of the semi-analytic binary evolution models used to generate the binary populations and corresponding GWB spectra. For example, the usage of a fixed-time phenomenological binary evolution model is forcing a particular relationship between typical binary masses and the degree of low-frequency spectral turnover. Another model, in which the degree of environmental coupling scales di ff erently with binary mass (or similarly, host-galaxy properties; e.g., Kelley et al. 2017a), may produce di ff erent dependencies and thus di ff erent posteriors. We are also assuming a fixed M BHM bulge relationship for all redshifts, while the canonical M BHM bulge relationship in the literature is specifically calibrated to the local Universe. Our values of µ and ϵµ for both the Phenom and GWOnly models should thus be interpreted as 'redshift-averaged' quantities.", '4.4. Inferred Properties of the SMBH Binary Populations': "While a large amount of information is encoded in the GWBspectra, there are numerous degeneracies - particularly in the current low signal-to-noise regime. For example, given a particular GWB spectral shape, a certain GWB amplitude can be produced by a large number of lower-mass SMBH binaries, or a small number of higher-mass SMBH binaries. To determine the characteristic properties of SMBH bina-  \nries contributing to the GWB, we calculate the distribution of GW-weighted average binary parameters. We use 1000 draws from the posteriors of the Phenom + Uniform fits. For each draw, the h 2 c -weighted parameters are calculated over 100 realizations. This gives a distribution of average parameters for each draw and each realization, which are plotted in Figure 11. As in all of our analysis, we fit binary population models to only the lowest five frequency bins in the 15 yr data set. However, in order to better visualize the trends in binary population parameters with GW frequency, Figure 11 shows the Phenom + Uniform library priors and posteriors for ten frequency bins.  \nThe GWB is characterized by the most massive SMBH binaries in the Universe with M ≳ 10 8 . 5 M ⊙ , and extending to just above 10 10 M ⊙ at the lowest frequency bins. At higher frequencies, as binaries evolve more quickly and fewer binaries occupy each frequency bin, these most massive systems become rarer and the typical masses decrease. Because of the trend in mass, the typical separations of binaries decreases more rapidly than f -2 / 3 as would be expected for a fixed mass. The binary total masses are the most strongly constrained parameters when comparing between the library priors and the posteriors. This is unsurprising given (a) the strong dependence of GW strain of binary mass, and (b) the numerous varying model parameters that a ff ect the masses. The mass-ratio distributions, on the other hand, are nearly constant across the band, and narrowly localized for both the priors and posteriors. Typical binary mass ratios are almost entirely above q ∼ 0 . 5. Note that this is determined primarily by our fiducial parameters for mass-ratio dependence in the GPF and GMT 7 . For the latter in particular, the GMT scales as q -1 , which strongly disfavors extreme mass-ratio pairs.  \nAcross the PTA band, the binaries producing the GWB are typically at many hundreds to a few thousands of comoving Mpc ( z ≈ 0 . 15 - 0 . 9). Average redshift posteriors are higher than the priors due to fits preferring shorter binary lifetimes. Binary separations are tightly constrained by the strong constraints on the binary masses. Typical separations are just below 10 -1 pc at the lowest frequency bin ( ≈ 2 nHz), down to just below 10 -2 pc at the tenth frequency bin ( ≈ 20 nHz). Projecting these separations at the cosmological distances of the binaries leads to angular separations of tens of µ as.  \nHaving explored the distributions of GWB-weighted binary properties from our fiducial model (Figure 11), we now examine which binary parameter ranges contribute most to the GWB signal. In Figure 12, we show the fraction of the GWB contributed by di ff erent portions of the binary pop-  \nFigure 10. One-dimensional binary evolution parameter posteriors from fitting against the 15 yr HD-w / MP + DP + CURN free spectrum. Priors (dotted) and posteriors (solid) are shown for the Phenom model with both uniform and astrophysical priors. Note that the astrophysical and uniform priors are identical for τ f and ν inner (see Table B1). The M BHM bulge parameters are the least constrained (i.e., the posteriors largely resemble the priors), while the GSMF parameters are the most constrained. This is consistent with our expectations given each parameter's impact on the simulated GWB (see Figure 4). However, these 1D projections can miss significant 2D constraints due to various projection e ff ects (e.g., the combined M BHM bulge mass normalization and m ψ, 0 2D space is much more constrained than each parameter individually, as seen in Figure 9). Overall, there is broad consistency between the posteriors for the Phenom models, regardless of which prior is used, showing that these constraints are coming directly from the GWB spectrum itself and are not based solely on our binary population model construction. Distributions are also shown for the GWOnly model, which gives posteriors broadly consistent with the phenomenological models. However, they are peaked at slightly lower parameter values, indicating an overall lower amplitude for their simulated GWB. This is most likely due to this model's rigidity, which does not allow for significant deviations from a power-law-like GWB and leads to an inability to capture the behavior of the lowest frequency bin in the 15 yr free spectrum (see Figure 8).  \n<!-- image -->  \nulation. Using the same posteriors as were sampled previously, we plot the fraction of h 2 c that is contributed by binaries at each binary parameter value (solid lines). The results are also separated by GW frequency (colors). We compare these GWB fractions to the total number of binaries in our simulated populations that are emitting at a given frequency (dashed lines). There is a stark di ff erence between the number of binaries emitting at each frequency as a function of M tot and their relative contribution to the GWB signal. Lower-mass binaries are far more common, but the signal is dominated by the rare, massive black holes near M tot = 10 9 . 5 -10 10 M ⊙ . For example, binaries between M tot = 10 9 . 2 -10 10 . 4 M ⊙ make up 2 . 6 × 10 -4 of all binaries emitting in the lowest frequency band, but they make up 73% of the signal at that frequency. Similarly, we are preferentially sensitive to the largest mass ratios.  \nAlthough the bulk of the GWB signal is made up of binaries at z > 0 . 4, it is most sensitive to the nearest binaries, relative to the underlying SMBH binary population. Since, all else being equal, the binaries that come from later galaxy mergers (lower z init) will enter the PTA band at later times (lower z final), the GWB source population is biased towards these lower-redshift sources. These biases a ff ect all frequencies relatively uniformly, with the mass bias slightly more  \npronounced at lower frequencies. This reflects the increased total number of binaries at low frequencies, which allows for the rarer, higher-mass SMBH binaries to dominate, as also seen in Figure 11.  \nNote the distinction between Figure 11, which shows h 2 c -weighted binary parameter distributions, and Figure 12, which shows the fractional contribution to the GWB of binaries with given parameters. These figures are very closely related, but Figure 11 shows the representative properties of GWB binaries, while Figure 12 shows the fractional GWB contribution of actual binaries.  \nOur simulated populations also contain individually loud, high-mass binaries that can contribute substantially to the GWB. These sources, apparent in Figure 12 as spikes in the GWB fraction at high total masses, are likely the types of sources that will be detectable as continuous-wave signals by PTAs (Arzoumanian et al. 2023a; Agazie et al. 2023a). While they can occur at a range of frequencies, they are typically expected mid-band ( f ∼ 3 -30 nHz) where the overall GWB amplitude has dropped somewhat but a sizeable population of binaries remains (Kelley et al. 2018; B'ecsy et al. 2022).  \nThe dashed lines in the leftmost panel of Figure 12 indicate the general shape of the mass function of binaries contributing to the GWB. In Figure 13, we explicitly calculate  \nseparation  \nFigure 11. Distributions of GWB-weighted binary parameters versus GW frequency for populations drawn from the Phenom + Uniform library priors (left-hand violins in gray), as well as from the posteriors after fitting to the first five frequency bins of the 15 yr NANOGrav data (right-hand violins in blue). For each population sample, the h 2 c -weighted averages of each parameter are calculated, and the distributions of those averages are plotted here. The GWB favors high total binary masses, especially in the lowest frequency bins, and typical binary separations range from ∼ 10 -1 -10 -2 pc between 2 and 20 nHz.  \n<!-- image -->  \nthe implied mass function of SMBH binaries at frequencies within the PTA band for several redshifts ( z final). Distributions are shown for populations drawn from a sample of the Phenom + Uniform parameter space priors and from our 15 yr spectral fit posteriors. For comparison, we also show the mass functions for non-binary SMBHs at each redshift, and in the rightmost panel, we show the fraction of SMBHs at each redshift that are in binaries. At z ≥ 1 . 0, the implied mass function is consistent with a wide range of values, indicating weak PTA constraints on high-redshift SMBH binary populations. This partly reflects the steep drop in binary fraction with increasing redshift; binaries that emit in PTA bands are ∼ 10 - 100 times rarer at z = 1 . 5 than at z = 0 . 25.  \nBecause of this, the binary number density increases over time (note that the non -binary number density also increases with time as required for monotonic SMBH growth, but the  \nevolution of the non-binary mass function is slight compared to the binary mass function evolution). The shape of the binary mass function also significantly evolves, such that at lower redshifts we see a much clearer turnover at M ∼ 10 9 M ⊙ that more closely traces the shape of the non-binary mass function. The binary mass function is also much more tightly constrained at low redshifts than at high redshifts, especially at the high-mass end. These ≳ 10 9 M ⊙ binaries at z < 1 are precisely the objects that comprise the bulk of the GWB (Figure 12).  \nFigures 11 - 13 examine the properties of binaries that emit in PTA bands. In Figure 14, we show the fraction of all binaries that reach the lowest frequency bin of the 15 yr data before redshift zero. Although binaries in the PTA bands will not reach coalescence on human timescales, the binary lifetime from PTA frequencies to merger is significantly shorter than a Hubble time. 8 We therefore use the fraction of binaries reaching PTA frequencies as a proxy for the fraction of systems that coalesce entirely before redshift zero.  \nWhether a given binary coalesces or stalls is determined by its formation redshift combined with the binary evolution time to reach the PTA band. Fitting to the 15 yr data strongly favors short binary lifetimes, which drives the di ff erence between priors and posteriors. The redshift ( bottom panel) at which the coalescing fraction reaches zero marks the redshift at which the look-back time of the Universe matches the binary lifetime of the model. The median posterior value of this lifetime is ≈ 2 . 8 Gyr, corresponding to a redshift of z ≈ 0 . 25-where the median coalescing fraction reaches zero. The gradual increase of coalescing fraction with redshift after this point is due to the additional delays from the GMTs. Coalescing fractions reach unity once the combined binary lifetimes and GMTs are longer than the look-back time.  \nThe gradual increase in coalescing fraction with mass ratio is due primarily to the GMT's strong dependence on mass ratio ( T gal-gal ∼ q -1 ⋆ ) in the fiducial model. The more gentle decline in coalescing fraction at the highest total masses is due to more extreme mass-ratio systems, as seen in the comparison between dashed lines ( q > 0 . 2) and solid lines (all q ). This is caused by a combination of increased GMTs, and also the increased binary coalescing times within the PTA band-produced by more extreme mass-ratio systems that have longer GW-driven inspiral times.  \nWe compare the Phenom model to a higher-dimensional library ( Phenom-Ext ) which includes variations in two of the GMT parameters, as well as others, in Appendix C and Figure C1. Generally, all of the recovered posteriors are consis-  \nFigure 12. Comparison of binary parameters by their fractional contribution to the GWB ( h 2 c ; solid), and their total number of binaries (dashed). Each line is colored by the GW frequency at which the binaries emit. The third panel ( z init), corresponds to the 'initial' galaxy merger redshift, while the fourth panel ( z final) is the redshift at which the binary is emitting GWs. This figure demonstrates that the GW signal is produced by a relatively small and highly biased sub-sample of a much larger population.  \n<!-- image -->  \nFigure 13. Left four panels show comoving number density per logarithmic interval of mass. Densities are shown for binary black holes in the NANOGrav frequency band at the redshift indicated in the upper-right panel for the priors (grey dotted line) and the posteriors (blue solid line) derived from our Phenom population models. Densities are also shown for the total SMBH population (teal dash-dotted line). Shaded regions show the 68% distributions. The rightmost panel shows the binary fraction as a function of mass for the four selected redshifts. The implied mass functions of our posteriors prefer relatively high density of black holes larger than M = 10 9 M ⊙ , and our posteriors are overall more confined than our priors at lower redshifts.  \n<!-- image -->  \ntent between the di ff erent libraries for parameters that they have in common, suggesting that our choices of fiducial parameters are su ffi ciently representative of the binary evolution parameter space. The posteriors on the additional parameters themselves are generally broad. The exception is the GMT parameters which, like the phenomenological evolution parameters, strongly favor shorter lifetimes as a way of producing higher GWB amplitudes.", '5. DISCUSSION': "The NANOGrav PTA has detected a common-spectrum correlated stochastic process that is consistent with an astrophysical GWB. In our 15 yr data set (NG15; NG15detchar), we find evidence of the Hellings-Downs correlations that would definitively mark this signal as GW in origin (NG15gwb). In this paper, we have presented analysis of the NANOGrav 15 yr data set under the assumption that  \nthese data represent a GWB produced by SMBH binaries. With reasonable choices of astrophysical parameters governing galaxy masses, galaxy mergers, SMBH masses, and SMBH binary inspiral timescales, we are able to reproduce the inferred GWB amplitude and spectral shape. We find that the data are suggestive of a GWB spectral turnover at low frequencies, as expected for binary inspiral driven by astrophysical environments. However, the broad free-spectrum posteriors from the 15 yr data are still consistent with the canonical α = 2 / 3 ( γ = 13 / 3) power law expected for GWdriven inspiral.  \nFigure A1 compares the posteriors for the GWB amplitude and spectral index inferred from the 15 yr data with a wide variety of GWB model predictions in the literature (see also Table A1). Although the inferred GWB amplitudes are within the range of some of these model predictions, they  \nFigure 14. Fraction of binary systems that reach the PTA band before redshift zero in our simulated populations. The top panel shows the total binary mass M , the middle panel shows the binary mass ratio q , and the bottom panel shows the initial binary redshift (defined as the galaxy merger redshift) z init. Solid lines show median values over 1000 samples from the priors (grey) and posteriors (blue) after fitting to the 15 yr data. The shaded regions correspond to the 50% inter-quartile range. The dashed lines correspond to binary subsets, with q > 0 . 2 in the top panel, with M > 10 8 M ⊙ in the middle panel, and both cuts in the bottom panel. Fits to the 15 yr data strongly favor shorter binary lifetimes, allowing larger fractions of binaries to coalesce. The nearly step-function behavior in the bottom panel marks the redshift at which the look-back time matches the sample's binary-lifetime parameter.  \n<!-- image -->  \nlie at the high-amplitude end of this range. The implied GWB amplitude from the NANOGrav data therefore indicates that SMBH binary model parameters di ff er from standard expected values, although still remaining within reasonable bounds.  \nIn this analysis, we have generated simulated populations of SMBH binaries and GWB spectra and fit them to the observed 15 yr signal. Our fiducial models explore a sixdimensional space of binary evolution parameters. Relative to typically assumed values for these parameters, our results indicate that the inferred GWB amplitude could be achieved with short binary hardening timescales, higher galaxy number densities (translating to higher galaxy merger rates), or higher normalization of the M BHM bulge relation. This may be accomplished if multiple parameters di ff er somewhat  \nfrom standard expectations, or if a small number of parameters di ff er more significantly.  \nOur models also demonstrate that the GWB signal is strongly dominated by the most massive, high-mass-ratio SMBH binaries, even among the subset of SMBH binaries emitting in PTA bands. The binaries contributing to the GWB form at typical redshifts of z ≈ 0 . 15 - 0 . 9. Their typical separations (which are tightly constrained via the SMBH masses) range from ∼ 0 . 1 - 0 . 01 parsec; this corresponds to binary angular separations of tens of µ as. Owing to the short binary lifetimes preferred by the 15 yr data set, most of these binaries will merge by z = 0-the coalescing fraction is near unity for binaries that form by z ∼ 0 . 25. In addition, we note that our simulated binary populations contain loud, high-mass continuous-wave sources that could be detected above the GWB.  \nBecause we are currently in the low signal-to-noise regime of GWB observations, we are still limited in our ability to make stringent parametric constraints. In this analysis, our constraints on the binary population inferred from the GWB spectral shape and amplitude are dependent on both our choices of priors and on which 15 yr GWB measurements are used. NANOGrav continues to collect data from an ever increasing number of pulsars. 43 pulsars were included in the 12.5-year 9 analysis (Arzoumanian et al. 2020), versus 67 pulsars in the current 15 yr analysis. As of summer 2023, we are timing roughly 75 pulsars with a total baseline of over 17 years. Also, NANOGrav data are currently being combined with those from other PTAs to create a new IPTA data set that will contain over 100 pulsars (Antoniadis et al. 2022). These e ff orts will improve our GWB measurement accuracy, along with our ability to constrain SMBH binary physics. The theoretical forecasts from Pol et al. (2021), for example, suggest that analyses using the future NANOGrav 20-year data set would be far more constraining than the 15 yr data. In that case, the authors found that subtle di ff erences in the degree of environmental coupling could be distinguished based on spectra with nearly identical reference amplitudes, but di ff ering low-frequency spectral shapes.  \nWhile we are unable to definitively attribute the inferred GWBsignal to SMBH binaries at the current signal-to-noise, we show that all of the signal's features are consistent with binaries. Nonetheless, many other possible origins of the GWB have been proposed, as detailed in NG15newphys. It is worth emphasizing that SMBH binaries must necessarily form throughout the Universe as a natural product of galaxy mergers. If the inferred GWB is not dominantly produced by SMBH binaries, the lack of GW signal from inspiraling SMBH binaries must somehow be accounted for. One  \npossibility is that SMBH binaries usually stall outside of the GW-driven regime of inspiral, which could occur if gas- and stellar-driven processes are insu ffi cient to bring binaries to the GW-dominated regime. If indeed the so-called 'final parsec problem' (e.g., Begelman et al. 1980) lengthens most inspiral timescales to a Hubble time or longer, the resulting GWB from binaries could be attenuated to amplitudes well below the inferred 15 yr signal. Even in this pessimistic case, multiple studies have suggested that triple-SMBH interactions would still produce a detectable GWB signal (Volonteri et al. 2003b; Ho ff man & Loeb 2007; Ryu et al. 2018; Bonetti et al. 2018b). In either case, this would also imply the existence of a large population of stalled SMBH binaries in the local Universe.  \nAdditional data is required to resolve the origin of the GWB. One of the strongest distinguishing features between di ff erent source models is the significantly higher degree of anisotropy for binaries as opposed to new physical processes (NG15newphys). From the 15 yr data set, the first limits on anisotropy have now been placed (Agazie et al. 2023b). While the limits are still consistent with astrophysical expectations for binary populations (e.g. Sato-Polito & Kamionkowski 2023), they will become significantly more constraining over time (Ali-Ha¨ımoud et al. 2021; Pol et al. 2022).  \nEventually, individual 'continuous-wave' GW sources will also become distinguishable above the GWB, if it is indeed produced by binaries. Di ff erent models have produced a variety of expectations for the plausibility of continuous-wave source detection in the near future (Sesana et al. 2009b; Rosado et al. 2015; Mingarelli et al. 2017; Kelley et al. 2018; B'ecsy et al. 2022). A search for continuous-wave sources has yielded improved upper limits on their occurrence rates in both the 12.5-year (Arzoumanian et al. 2023b) and 15 yr data (Agazie et al. 2023a). A continuous-wave detection would present the exciting possibility of multi-messenger detections: GWs from a single SMBH binary for which an electromagnetic counterpart could be identified. Such a multimessenger source would provide a wealth of information about the origin of low-frequency GWs, the astrophysical environment of SMBH binaries, and SMBH accretion processes (Kelley et al. 2019a).  \nIn the next decade, the Laser Interferometer Space Antenna (LISA) will begin operation in the ∼ millihertz band, sensitive to the merger of SMBH binaries with masses in a range between ∼ 10 4 -10 8 M ⊙ out to z ∼ 10 (Amaro-Seoane et al. 2023). While LISA promises to reveal the elusive formation mechanism of massive black-hole seeds in the early Universe, this range of binary masses is more poorly constrained and far more challenging to model than the higher-mass PTA binaries. The approaches and analyses developed for studying PTA GW sources will be crucial for paving the way for  \nLISA science. The identification of LISA electromagnetic counterparts will also be much more di ffi cult, further motivating the development of techniques for PTA sources and generalizing them to signals at lower masses.  \nIf the GWB signal is indeed produced by astrophysical binaries, it will be the first proof that SMBH binaries do indeed form, evolve to sub-parsec separations, and eventually coalesce. These systems will join the new landscape of multimessenger GW astrophysics, o ff ering the opportunity to study the most extreme and energetic environments in the Universe and to probe the closely coupled co-evolution of galaxies and their nuclear engines. If instead the GWB has a di ff erent cosmological origin, it may provide answers to the most outstanding questions in fundamental physics that challenge the standard model and Λ CDM. In either case, PTAs have cracked open the era of low-frequency GW astronomy.", 'ACKNOWLEDGMENTS': "Author contributions. An alphabetical-order author list was used for this paper in recognition of the fact that a large, decade timescale project such as NANOGrav is necessarily the result of the work of many people. All authors contributed to the activities of the NANOGrav collaboration leading to the work presented here, and reviewed the manuscript, text, and figures prior to the paper's submission. Additional specific contributions to this paper are as follows. L.Z.K. proposed and coordinated the project. Development of the holodeck population modeling framework was led by L.Z.K., with contributions from A.C-C., D.W., E.C.G., J.M.W., K.G., M.S.S., & S.C. Libraries were generated by J.S., K.G., L.B., L.Z.K., & N.G-D. Development of the Gaussian process modules was led by J.S., with contributions from D.W. & J.M.W. Development of the fitting infrastructure was led by W.G.L., with contributions from J.S. Gaussian process training was performed by D.W. & J.S., with diagnostics performed by J.M.W. The 15 yr GWB data products were provided by the NANOGrav Detection Working Group. The analysis was performed by J.S., K.G., & L.Z.K., who produced most of the figures, and by L.B. Additional analysis code and interpretation were provided by D.J.D. and E.C.G. The collection of literature predictions, and production of the literature review figure, was performed by M.C. The paper was written by J.C.R., J.S., K.G., L.B., L.Z.K., & M.C. A.Mi, B.B., D.J.D., J.S.H., M.T.L., M.V., S.J.V., S.R.T., T.J.W.L., & T.D. provided feedback on the analysis and manuscript.  \nAcknowledgments. The NANOGrav collaboration receives support from National Science Foundation (NSF) Physics Frontiers Center award numbers 1430284 and 2020265, the Gordon and Betty Moore Foundation, NSF AccelNet award number 2114721, an NSERC Discovery Grant, and CIFAR. NANOGrav is part of the International Pulsar Timing Array (IPTA).  \nThe Arecibo Observatory is a facility of the NSF operated under cooperative agreement (AST-1744119) by the University of Central Florida (UCF) in alliance with Universidad Ana G. M'endez (UAGM) and Yang Enterprises (YEI), Inc. The Green Bank Observatory is a facility of the NSF operated under cooperative agreement by Associated Universities, Inc. The National Radio Astronomy Observatory is a facility of the NSF operated under cooperative agreement by Associated Universities, Inc.  \nThis work was conducted in part using the resources of the Advanced Computing Center for Research and Education (ACCRE) at Vanderbilt University, Nashville, TN. This work received computational support from UCF's Advanced Research Computing cluster, operated by the University of Central Florida. This research used the Savio computational cluster resource provided by the Berkeley Research Computing program at the University of California, Berkeley (supported by the UC Berkeley Chancellor, Vice Chancellor for Research, and Chief Information O ffi cer). This research was supported in part through computational resources and services provided by Advanced Research Computing at the University of Michigan, Ann Arbor. This work was also conducted in part using the Thorny Flat HPC Cluster at West Virginia University (WVU), which is funded in part by National Science Foundation (NSF) Major Research Instrumentation Program (MRI) Award number 1726534, and West Virginia University. This work utilized the Alpine high performance computing resource at the University of Colorado Boulder. Alpine is jointly funded by the University of Colorado Boulder, the University of Colorado Anschutz, Colorado State University, and the National Science Foundation (award 2201538). This work utilized the Blanca condo computing resource at the University of Colorado Boulder. Blanca is jointly funded by computing users and the University of Colorado Boulder.  \nL.B. acknowledges support from the National Science Foundation under award AST-1909933 and from the Research Corporation for Science Advancement under Cottrell Scholar Award No. 27553. P.R.B. is supported by the Science and Technology Facilities Council, grant number ST / W000946 / 1. S.B. gratefully acknowledges the support of a Sloan Fellowship, and the support of NSF under award #1815664. The work of R.B., R.C., D.D., N.La., X.S., J.P.S., and J.T. is partly supported by the George and Hannah Bolinger Memorial Fund in the College of Science at Oregon State University. M.C., P.P., J.C.R. and S.R.T. acknowledge support from NSF AST-2007993. M.C. and N.S.P. were supported by the Vanderbilt Initiative in Data Intensive Astrophysics (VIDA) Fellowship. K.Ch., A.D.J., and M.V. acknowledge support from the Caltech and Jet Propulsion Laboratory President's and Director's Research and Development Fund. K.Ch. and A.D.J. acknowledge support from  \nthe Sloan Foundation. Support for this work was provided by the NSF through the Grote Reber Fellowship Program administered by Associated Universities, Inc. / National Radio Astronomy Observatory. Support for H.T.C. is provided by NASA through the NASA Hubble Fellowship Program grant #HST-HF2-51453.001 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. K.Cr. is supported by a UBC Four Year Fellowship (6456). M.E.D. acknowledges support from the Naval Research Laboratory by NASA under contract S-15633Y. T.D. and M.T.L. are supported by an NSF Astronomy and Astrophysics Grant (AAG) award number 2009468. E.C.F. is supported by NASA under award number 80GSFC21M0002. G.E.F., S.C.S., and S.J.V. are supported by NSF award PHY-2011772. K.A.G. and S.R.T. acknowledge support from an NSF CAREER award #2146016. The Flatiron Institute is supported by the Simons Foundation. S.H. is supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE-1745301. N.La. acknowledges the support from Larry W. Martin and Joyce B. O'Neill Endowed Fellowship in the College of Science at Oregon State University. Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). D.R.L. and M.A.Mc. are supported by NSF #1458952. M.A.Mc. is supported by NSF #2009425. C.M.F.M. was supported in part by the National Science Foundation under Grants No. NSF PHY-1748958 and AST-2106552. A.Mi. is supported by the Deutsche Forschungsgemeinschaft under Germany's Excellence Strategy - EXC 2121 Quantum Universe - 390833306. P.N. acknowledges support from the BHI, funded by grants from the John Templeton Foundation and the Gordon and Betty Moore Foundation. The Dunlap Institute is funded by an endowment established by the David Dunlap family and the University of Toronto. K.D.O. was supported in part by NSF Grant No. 2207267. T.T.P. acknowledges support from the Extragalactic Astrophysics Research Group at Eotvos Lor'and University, funded by the Eotvos Lor'and Research Network (ELKH), which was used during the development of this research. S.M.R. and I.H.S. are CIFAR Fellows. Portions of this work performed at NRL were supported by ONR 6.1 basic research funding. J.D.R. also acknowledges support from start-up funds from Texas Tech University. J.S. is supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST2202388, and acknowledges previous support by the NSF under award 1847938. C.U. acknowledges support from BGU (Kreitman fellowship), and the Council for Higher Education and Israel Academy of Sciences and Humanities (Excellence fellowship). C.A.W. acknowledges support from CIERA, the", 'Facilities: Arecibo, GBT, VLA': 'Software: acor, astropy (Astropy Collaboration et al. 2018), ce ff yl (Lamb et al. 2023), chainconsumer (Hinton 2016), cython (Behnel et al. 2011), emcee (ForemanMackey et al. 2013), enterprise (Ellis et al. 2020), enterprise extensions (Taylor et al. 2021), George (Ambikasaran et al. 2015b), holodeck (Kelley et al. 2023, in prep), jupyter (Kluyver et al. 2016), kalepy (Kelley 2021), matplotlib (Hunter 2007), numba (Lam et al. 2015), numpy (van der Walt et al. 2011), ptmcmc (Ellis & van Haasteren 2017), scipy (Virtanen et al. 2020),', 'A. GWB PREDICTIONS IN THE LITERATURE': 'In Figure A1 and Table A1, we summarize model predictions of GWB amplitudes from the literature. For each reference, amplitude, predictions are cited at frequencies of GLYPH<0> 1 yr GLYPH<1> -1 and GLYPH<0> 10 yr GLYPH<1> -1 . While it is most common to reference GWB amplitudes', 'NANOG rav 15 yr A strophysical I nterpretation': 'Table A1. Literature predictions for the amplitudes A yr and A 10yr of the GWB at frequencies of GLYPH<0> 1 yr GLYPH<1> -1 and GLYPH<0> 10 yr GLYPH<1> -1 , respectively. Column (1) gives the literature reference, columns (2) & (3) give the 16 th & 84 th percentiles for the uncertainty region of the corresponding prediction of A yr, and columns (4) & (5) give the 16 th &84 th percentiles for A 10yr predictions. Figure A1 provides a visual comparison of these predictions to the NANOGrav 15 yr results.  \nat f = 1 yr -1 , NANOGrav constraints are derived primarily at much lower frequencies. Small deviations in power-law indices at low frequencies can lead to large changes in the corresponding A yr amplitudes. The amplitudes at f = GLYPH<0> 10 yr GLYPH<1> -1 are much more representative of current PTA constraints. Figure A1 also shows the posterior distributions for the inferred GWB amplitudes obtained from the 15 yr data NANOGrav when assuming a power-law model; results are shown for both the HD-w / MP + DP + CURN and the HD-DMGP models. Numerous literature model predictions overlap with the inferred GWB amplitudes from the 15 yr data, especially when comparing the amplitude at f = GLYPH<0> 10 yr GLYPH<1> -1 . Even so, the 15 yr results lie at the higher-amplitude end of the predicted ranges.', 'B. SEMI-ANALYTIC MODEL PARAMETERIZATIONS': 'In § 3.2, we presented the semi-analytic models and the underlying equations used in this paper. In Table B1, we detail all of the model parameters, showing the fiducial values for fixed parameters as well as the prior distributions for parameters that are varied when fitting the models to the 15 yr data. These models are summarized Table B2.', 'C. HIGHER-DIMENSIONAL PARAMETER SPACES': 'We have shown that in our fiducial, six-dimensional Phenom parameter space, a wide range of semi-analytic model parameters are consistent with current measurements of the GWB. In this library, a large number of additional parameters are held fixed to astrophysically motivated values. When fitting to the 15 yr GWB data, uniform priors on the included parameters are typically used. We have compared these results to fits of the same parameter space, but adopting more informed priors based on the astronomical literature.  \nCurrently it is not feasible to run MCMC fits using Gaussian processes that have been trained to significantly larger parameter spaces. However, we have generated higher dimensional libraries and directly evaluated them against the 15 yr GWB spectra at', 'REFERENCES': 'Afzal, A., et al. 2023, ApJL, doi: 10.3847 / 2041-8213 / acdc91  \nAgazie, G., Alam, M. F., Archibald, A. M., et al. 2023a, Astrophys.  \nJ. Lett., doi: 10.3847 / 2041-8213 / acda9a  \nAgazie, G., Anumarlapudi, A., Archibald, A. M., et al. 2023b,  \nAstrophys. J. Lett., doi: 10.3847 / 2041-8213 / acdac6  \n-. 2023c, Astrophys. J. Lett., doi: 10.3847 / 2041-8213 / acda88  \nAgazie, G., et al. 2023a, ApJL  \n-. 2023b, ApJL  \n- Aigrain, S., & Foreman-Mackey, D. 2022, arXiv e-prints, arXiv:2209.08940, doi: 10.48550 / arXiv.2209.08940\n- Alam, M. F., Arzoumanian, Z., Baker, P. T., et al. 2021, ApJS, 252, 4, doi: 10.3847 / 1538-4365 / abc6a0\n- Ali-Ha¨ımoud, Y., Smith, T. L., & Mingarelli, C. M. F. 2021, PhRvD, 103, 042009, doi: 10.1103 / PhysRevD.103.042009\n- Amaro-Seoane, P., et al. 2023, Living Rev. Rel., 26, 2, doi: 10.1007 / s41114-022-00041-y\n- Ambikasaran, S., Foreman-Mackey, D., Greengard, L., Hogg, D. W., & O\'Neil, M. 2015a, IEEE Transactions on Pattern Analysis and Machine Intelligence, 38, 252, doi: 10.1109 / TPAMI.2015.2448083\n- -. 2015b, IEEE Transactions on Pattern Analysis and Machine Intelligence, 38, 252, doi: 10.1109 / TPAMI.2015.2448083\n- Antoniadis, J., Arzoumanian, Z., Babak, S., et al. 2022, MNRAS, 510, 4873, doi: 10.1093 / mnras / stab3418  \nAntonini, F., & Merritt, D. 2012, ApJ, 745, 83,  \ndoi: 10.1088  \n/  \n0004-637X  \n/  \n745  \n/  \n1  \n/  \n83  \n- Arzoumanian, Z., Baker, P. T., Brazier, A., et al. 2018, ApJ, 859, 47, doi: 10.3847 / 1538-4357 / aabd3b\n- Arzoumanian, Z., Baker, P. T., Blumer, H., et al. 2020, ApJL, 905, L34, doi: 10.3847 / 2041-8213 / abd401\n- Arzoumanian, Z., Baker, P. T., Blecha, L., et al. 2023a, arXiv e-prints, arXiv:2301.03608, doi: 10.48550 / arXiv.2301.03608 -. 2023b, arXiv e-prints, arXiv:2301.03608,  \ndoi: 10.48550  \n/  \narXiv.2301.03608  \n- Astropy Collaboration, Price-Whelan, A. M., Sip"ocz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847 / 1538-3881 / aabc4f\n- Bansal, K., Taylor, G. B., Peck, A. B., Zavala, R. T., & Romani, R. W. 2017, ApJ, 843, 14, doi: 10.3847 / 1538-4357 / aa74e1  \nBarnes, J. E., & Hernquist, L. 1992, ARA&A, 30, 705,  \ndoi: 10.1146  \n/  \nannurev.aa.30.090192.003421  \n- B\'ecsy, B., Cornish, N. J., & Kelley, L. Z. 2022, ApJ, 941, 119, doi: 10.3847 / 1538-4357 / aca1b2\n- Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Nature, 287, 307, doi: 10.1038 / 287307a0\n- Behnel, S., Bradshaw, R., Citro, C., et al. 2011, Computing in Science and Engineering, 13, 31, doi: 10.1109 / MCSE.2010.118\n- Berczik, P., Merritt, D., Spurzem, R., & Bischof, H.-P. 2006, ApJL, 642, L21, doi: 10.1086 / 504426\n- Bluck, A. F. L., Conselice, C. J., Buitrago, F., et al. 2012, ApJ, 747, 34, doi: 10.1088 / 0004-637X / 747 / 1 / 34\n- Bluck, A. F. L., Mendel, J. T., Ellison, S. L., et al. 2014, MNRAS, 441, 599, doi: 10.1093 / mnras / stu594\n- Bogdanovi\'c, T., Miller, M. C., & Blecha, L. 2022, Living Reviews in Relativity, 25, 3, doi: 10.1007 / s41114-022-00037-8\n- Bonetti, M., Haardt, F., Sesana, A., & Barausse, E. 2016, MNRAS, 461, 4419, doi: 10.1093 / mnras / stw1590  \nZhu, X.-J., Cui, W., & Thrane, E. 2019, MNRAS, 482, 2588, doi: 10.1093 / mnras / sty2849  \nZrake, J., Tiede, C., MacFadyen, A., & Haiman, Z. 2021, ApJL, 909, L13, doi: 10.3847 / 2041-8213 / abdd1c'}