Beyond the Background: Gravitational Wave Anisotropy and Continuous Waves from Supermassive Black Hole Binaries

Pulsar timing arrays have found evidence for a low-frequency gravitational wave background (GWB). Assuming the GWB is produced by supermassive black hole binaries (SMBHBs), the next gravitational wave (GW) signals astronomers anticipate are Continuous Waves (CWs) from single SMBHBs and their associated GWB anisotropy. The prospects for detecting CWs and anisotropy are highly dependent on the astrophysics of SMBHB populations. Thus, information from single sources can break degeneracies in astrophysical models and place much more stringent constraints than the GWB alone. We simulate and evolve SMBHB populations, model their GWs, and calculate their anisotropy and detectability. We investigate how varying components of our semi-analytic model, including the galaxy stellar mass function, the SMBH--host galaxy relation ($M_\mathrm{BH}$--$M_\mathrm{bulge}$), and the binary evolution prescription impact the expected detections. The CW occurrence rate is greatest for few total binaries, high SMBHB masses, large scatter in $M_\mathrm{BH}$--$M_\mathrm{bulge}$, and long hardening times. The occurrence rate depends most on the binary evolution parameters, implying that CWs offer a novel avenue to probe binary evolution. The most detectable CW sources are in the lowest frequency bin for a 16.03-year PTA, have masses from $\sim\!\!10^9-10^{10}\mathrm{M}_\odot$, and are $\sim\!\!1$ Gpc away. The level of anisotropy increases with frequency, with the angular power spectrum over multipole modes $\ell$ varying in low-frequency $C_{\ell>0}/C_0$ from $\sim\!\!5\times 10^{-3}$ to $\sim\!\!2\times10^{-1}$, depending on the model; typical values are near current upper limits. Observing this anisotropy would support SMBHB models for the GWB over cosmological models, which tend to be isotropic.


INTRODUCTION
Supermassive black hole binaries (SMBHBs) are predicted to result from galaxy mergers.Two galaxies, each hosting a central supermassive black hole (SMBH) (Richstone et al. 1998), merge as predicted by hierarchical structure formation (Lacey & Cole 1993).Then, their SMBHs sink to the center of the merged galaxies via dynamical friction, become gravitationally bound, and form a binary with ∼pc separation.Stellar scattering and circumbinary disk torques harden the binary to small separations (∼10 −2 pc) (Begelman et al. 1980;Quinlan 1996;Cuadra et al. 2009;Kelley et al. 2017) beyond which they evolve primarily by emitting gravitational waves.
The superposition of these continuous waves (CWs) from many SMBHBs across the universe creates an incoherent stochastic gravitational wave background (GWB) (Rajagopal & Romani 1995;Jaffe & Backer 2003;Burke-Spolaor et al. 2019), like that for which pulsar timing arrays (PTAs) have recently found strong evidence (Agazie et al. 2023a;EPTA Collaboration et al. 2023;Reardon et al. 2023;Xu et al. 2023).In the likely scenario that the PTA-observed GWB is produced by SMBHBs (Agazie et al. 2023b;Antoniadis et al. 2023a), CWs from individual, loud SMBHBs are the next highly-anticipated GW signal PTAs could detect.PTA searches have yet to find a CW source (Agazie et al. 2023c;Antoniadis et al. 2023b), but simulation-based predictions suggest single source CWs could be detected within a few years of the GWB (Kelley et al. 2018).These single sources will likely brighten certain regions of the gravitational wave sky, inducing anisotropy in the background (Pol et al. 2022) before they can be individually resolved.
Cosmological models (cosmic inflation, phase transitions, cosmic strings, domain walls, etc.) for the GWB have also been suggested (Afzal et al. 2023;Antoniadis et al. 2023a).These are more likely to be isotropic.Thus, measuring anisotropy in the GWB would serve as compelling evidence for SMBHBs being the source.This anisotropy has been predicted using analytic (Mingarelli et al. 2013;Hotinli et al. 2019;Sato-Polito & Kamionkowski 2023), semi-analytic (Mingarelli et al. 2017), and simulation-based (Taylor & Gair 2013;Taylor et al. 2020;Bécsy et al. 2022;Agazie et al. 2023d) methods.We conduct the first study into what infor-mation content GW anisotropy contains about astrophysical models.Further, this paper offers the first look at how singlesource detection statistics and anisotropy are related.
Past works have predicted the amplitude and shape of the GWB using host galaxy populations generated from galaxy formation simulations (Kelley et al. 2017;Bécsy et al. 2022;Sykes et al. 2022;Bécsy et al. 2023), dark matter (DM) merger trees (Izquierdo-Villalba et al. 2022), galaxy catalogs (Mingarelli et al. 2017), or semi-analytic models (Sesana et al. 2008;Agazie et al. 2023b), and others have predicted single-source CWs using galaxy simulations (Kelley et al. 2018), DM merger trees (Sesana et al. 2009), and semianalytic SMBHB assembly models (Rosado et al. 2015).Such studies have historically focused on specific hardening processes (Kelley et al. 2018;Siwek et al. 2020) or accretion scenarios and SMBH-host galaxy relations (Sesana et al. 2005).To advance this field, we predict the parametric dependence of the likelihood and nature of low-frequency CW signals on the most complete SMBHB assembly and evolution models to date.This is the first systematic investigation of model parameter space and the information content of single CW sources.
We generate SMBHB populations using holodeck (Kelley et al. in prep.) as explained in §2.1, extract the loudest single sources, and calculate gravitational waves and binary properties of both the background and single sources as described in §2.2.We present the resulting characteristic strain spectra, total masses, and final comoving distances, for variations on several model components, including the galaxy stellar mass function (GSMF) ( §3.1.1),the SMBH-host relations ( §3.1.2),and the binary evolution ( §3.1.3).Then we calculate single source detection statistics for simulated PTAs using the methods described in §2.3.The resulting single source occurrence rates and predicted properties (mass, distance, and frequency) are given in §3.2.1 and §3.2.2, respectively.Finally, we calculate the GWB anisotropy from these SMBHB populations as described in §2.4, with the resulting angular power spectrum presented in §3.3.We discuss caveats to our model and future steps in §4 and summarize our key findings in §5.

Model for SMBHB Populations
Using holodeck (Kelley et al. in prep.)we assemble a population of galaxy mergers with comoving volumetric number density η gal-gal ≡ dN gal-gal /dV c (Chen et al. 2019), 10) We direct the reader to Agazie et al. (2023b) for a full description of the semi-analytic model components, including the galaxy pair fraction P and galaxy merger time T gal-gal ,  B1.
a The fiducial values are calculated as the mean across the varying ranges, which correspond to the uniform priors in Agazie et al. (2023b).
both of which are power-law functions of galaxy stellar mass m ⋆1 , galaxy mass ratio q ⋆ and initial redshift z ′ .
The components of the model that we investigate in this paper are: (1) the normalization ψ 0 and characteristic mass m ψ,0 of the galaxy stellar mass function (GSMF) Ψ, (2) the dimensionless mass normalization µ and intrinsic scatter ϵ µ of the SMBH mass-bulge mass (M BH -M bulge ) relation, and (3) the binary lifetime τ f and 'inner regime' power-law index ν inner of the phenomenological hardening model, each of which are summarized below.We study the effects of each of these six parameters in isolation, by independently varying one parameter across the range listed in Table 1 while fixing the five other parameters to the fiducial values listed there and all other model components to the fiducial values in Agazie et al. (2023b) Table B1.We examine a wide range of parameter space corresponding to the same range explored in Agazie et al. (2023b).While this range extends beyond currently predicted uncertainties, several model components remain poorly constrained by observations, and exploring above and below likely values allows us to clearly identify the impact of each model component and where GWs provide constraining power.
GSMF -The GSMF is the number density of galaxies per decade of stellar mass that determines the initial distribution of galaxies.We represent the GSMF as a single Schechter function (Schechter 1976), where m ⋆1 is the primary galaxy stellar mass; Ψ 0 , M ψ , and α ψ are phenomenological functions parameterized over redshift as in Chen et al. (2019) such that The population of merged galaxies is derived from this distribution through the galaxy pair fraction and galaxy merger time specified in Agazie et al. (2023b) using a fixed powerlaw dependence on mass, redshift, and mass ratio.While the GSMF is well constrained by current astrophysics, varying our GSMF normalization captures changes in the coefficients normalizing the pair fraction and merger time as well, thus representing the overall number of merged galaxies.
SMBH-Host Relation -The SMBH masses are related to their host galaxies' bulge masses M bulge by assuming an M BH -M bulge relation defined by dimensionless mass normalization log 10 µ and power-law index α µ , in addition to a random normal distribution of log 10 scatter N 0, ϵ µ with standard deviation ϵ µ : The bulge mass is calculated as a fraction of the total galaxy stellar mass = f ⋆,bulge •m ⋆1 , with f ⋆,bulge = 0.615 based on empirical observations from Lang et al. (2014) and Bluck et al. (2014).
Applying the M BH -M bulge relation in Eq. 4, the number density of merged galaxies in Eq. 1 translates to the number density of SMBHBs by Hardening -To model gravitational waves detectable by PTAs, the population of SMBHBs must be evolved in time and separation to rest frame orbital frequencies f p corresponding to GW frequencies of a few nHz.This binary hardening is described in terms of a rate of decreasing separation, da/dt = (da/dt) gw +(da/dt) phenom , i.e. the sum of a GW component and a phenomenological component, A double power law allows for distinct asymptotic behavior in the small-separation 'inner' regime and large-separation 'outer' regime, distinguished by a critical break separation a c .H a is a normalization factor, calibrated for every binary such that it has a total binary lifetime from initial separation a init to coalescence at the innermost stable circular orbit a isco of This serves as a self-consistent approach to modeling binary evolution, without depending on assumptions about the binary hardening processes or galactic environment.We also investigate the effects of varying our four GSMF and M BH -M bulge parameters for the GW-only model as in Agazie et al. (2023b), which is not self-consistent because GWs alone cannot bring the binaries to small enough separations to emit nHz GWs, but serves as a useful comparison.

Binary Properties and Gravitational Waves
The analytic model described in §2.1 determines a comoving volumetric number density of SMBHBs ∂ 3 η ∂M∂q∂z , from which we calculate a continuous number of SMBHBs per mass M, ratio q, redshift z (at the time of GW emission), and log rest-frame orbital frequency ln f p (Sesana et al. 2008): This continuous distribution sets a fractional expectation value for the number of binaries.In reality, gravitational waves are produced by a discrete population of binaries.We generate random universe realizations of this population by selecting a number of binaries N(M, q, z, f ) in each parameter bin of ∆M, ∆q, ∆z, and ∆(ln f p ) from a Poisson distribution (P) centered at the aforementioned expectation value, We assume circular orbits for all binaries and assign them the M, q, z, and f p values corresponding to their bin centers.
These define their chirp mass M ≡ (m 1 m 2 ) 3/5 /M 1/5 = Mq 3/5 / (1 + q) 6/5 comoving distance d c , observer-frame GW frequency f = (2 f p )/(1 + z), and (sky and polarization averaged) GW strain amplitude of (Finn & Thorne 2000) Because the loudest single source may not be the most detectable, depending on sky position, inclination, etc., the ten loudest single sources (SS, i.e. with the greatest h s ) in each frequency bin are then extracted from this population.Their individual characteristic strains are calculated as (Rosado et al. 2015) Here, ∆ f is the frequency bin width and arises when considering a finite number of sources in finite frequency bins N ∼ f * T ∼ f /∆ f , over an observing duration T .
The GWB is then calculated as the sum of gravitational waves from all background (BG) binaries Here, we define BG binaries to include all but the single loudest at each frequency because the most immediate observational application of this work will be the detection of one CW source, before PTAs can resolve multiple of them.When considering the GW-only model without phenomenological hardening, this in combination with the GW hardening rate (da/dt) gw , leads to the h c,BG ∝ f −2/3 power law often used as a comparison point for the characteristic strain spectra.In reality, we expect deviations from this power law not only due to the phenomenological hardening before GWs dominate the evolution (Kocsis & Sesana 2011), but also due to the discretization of sources where the power-law would otherwise predict fractional binaries (Sesana et al. 2008).We also calculate the characteristic mass, mass ratio, redshift, comoving distance, separation, and angular separation for the SMB-HBs contributing most to the background at each frequency.
To do so, we perform an h c,BG -weighted average over all BG binaries emitting at that frequency.

Detection Statistics
Given the h c,SS and h c,BG spectra, we calculate SS and BG detection statistics following the formalism in Rosado et al. (2015).This includes the background signal-to-noise ratio (SNR BG ), and detection probability (DP BG ), and each individual source's SNR (SNR SS,i ) and detection probability (DP SS,i ).The probability of detecting any single source is then (Rosado et al. 2015) and the expected number of detections for that realization is In this prescription, single source detection probabilities are given by integrating over the F e -statistic from some threshold Fe to infinity, where Fe is set to give a false alarm probability (FAP) of 10 −3 .Even with no signal present, the area under this curve will produce a nonzero total detection probability (DP SS ) equal to the FAP.Thus, 10 −3 is the lower limit on detection probabilities calculated in Eq. 14 and should be treated effectively as 0.
In light of the strong evidence for a GWB in current PTA data (Agazie et al. 2023a;EPTA Collaboration et al. 2023;Reardon et al. 2023;Xu et al. 2023), we study the single source detection probability under the same conditions that are likely to produce measurable GWB evidence by calibrating every realization to DP BG = 0.50.We use a white-noiseonly simulated PTA of 40 pulsars at randomly assigned sky positions, 16.03 yr duration (corresponding to Agazie et al. 2023a), and 0.20 yr cadence ∆t.Our fiducial method of calibration is to vary the level of white noise S WN , given by the error in pulsar times of arrival (ToAs) σ (Rosado et al. 2015): until achieving 0.49 < DP BG < 0.51.We calculate ⟨N SS ⟩ using the same pulsar positions and σ, with all characteristic strains except that of the source in question considered additional noise, Then we normalize for small variations around DP BG = 0.50 with For one realization of BG and SS characteristic strain, we calibrate the PTA to the background, then create 100 'sky realizations'-the random position, inclination, polarization, and phase assignment for single sources-and conduct SS detection statistics for each.By calculating detection statistics for 10 single sources for each frequency of each realization, we allow for the most detectable to depend on both strain amplitude and random location/orientation.This is repeated for 500 'strain realizations' of h c,BG and h c,SS -those created by Poisson sampling in Eq. 10, each with their own BG-calibrated PTA-to create 50,000 combined 'strain+sky realizations.'Next, we predict the most likely frequencies of detection by calculating the DP SS,i -weighted average frequency of all n loudest single sources across all realizations of a given model: with weighted standard deviation: The likely frequency of detection is sensitive to the shape of the PTA noise.Thus, we explore different noise models inspired by realistic sensitivity curves and intrinsic pulsar red noise in §A.

GWB Anisotropy
We measure the anisotropy corresponding to each model by decomposing a simulated GW sky into spherical harmonics, as in Agazie et al. (2023d).To generate this GW sky, we create a HEALpix map (Górski et al. 2005) of h 2 c /∆Ω (Ω being solid angle, or equivalently, pixel area) at each frequency of each realization by assigning the single sources to random pixels and distributing the remaining h c,BG evenly among all the pixels.Then, we can decompose this sky into an angular power spectrum of multipole modes ℓ and m each accompanied by a coefficient a ℓm such that the total GW power is the sum of each a ℓm times the real-valued spherical harmonic Y ℓm .The anafast code, via healpy (Zonca et al. 2019) calculates these coefficients with the estimator (Gorski et al. 1999) where N pix = 12N 2 side is the number of pixels, indexed by p at positions γ p , and f (γ p ) is h 2 c /∆Ω in each pixel.Using anafast, we calculate the corresponding angular power spectrum where C ℓ represents the measure of fluctuations (i.e.anisotropy) on the angular scale θ ≈ 180 deg /ℓ.C 0 represents the purely isotropic average component, thus we normalize our results by C ℓ /C 0 .This method is tested for 10, 100, and 1000 loudest single sources per frequency bin, by which point our results are insensitive to the addition of more single sources.By placing these sources randomly and treating the remaining signal as purely isotropic, we do not weigh in possible correlations with large-scale-structure, making this a conservative estimate for anisotropy.We also test the resolution and find C ℓ to be indistinguishable for N side = 8 up to N side = 32.Thus, we adopt N side = 8 (N pix = 768) to efficiently calculate the spherical harmonic decomposition for each realization and present the results in §3.3.

Characteristic Strain and Binary Properties
In this section we present the characteristic strain h c , total mass M, and final comoving distance d c of GWB and CW sources as a function of frequency.The first column of Fig. 1 includes three models with varying GSMF normalization: ψ 0 = −3.5, −2.5, and −1.5, while all other parameters remain fixed to their fiducial values listed in Table 1.Information about the CW sources is shown in green, for these three models, respectively.This includes the 68% confidence intervals (CIs, shaded regions) across 500 realizations of the single loudest source at each frequency.The 95th percentile of these sources' h c,SS and M SS , and the 5th percentile of these sources' d c,SS is also shown (points).For comparison, h c,BG and the h c,BG -weighted average properties (⟨M⟩ BG and ⟨d c ⟩ BG ) of the background (all but the loudest single sources at each frequency) are shown in corresponding shades of grey, with dashed lines representing their medians.The same is shown for models of varying m ψ,0 in the second column of Fig. 1, for µ and ϵ µ in Fig. 2, and for τ f and ν inner in Fig. 3.The following three sections describe the physical scenarios producing these results for each model component: the GSMF ( §3.1.1),the M BH -M bulge relation ( §3.1.2),and binary evolution ( §3.1.3).

GSMF
The GSMF parameters shape the masses of galaxies, and thus their residing SMBHBs for fixed M BH -M bulge .The masses in each SMBHB directly determine its strain amplitude as h s ∝ M 5/3 (Eq.11).When considering the characteristic strain from many background sources or from a single source sampled probabilistically from a distribution, h c,SS and h c,BG also have mass dependence in their GW strains, their hardening timescales, and their number density, such that h c ∝ M 5/6 η(M, q, z).By varying GSMF parameters, we examine the η(M, q, z) component of the mass distribution.Changes to the distribution of SMBHB masses also result secondarily in small d c variations due to the mass dependence of hardening rate versus frequency (described below).
Following the Schechter GSMF, the number of SMBHBs decreases with increasing mass.When this expectation value approaches zero, few random realizations contain a source in that bin.When ψ 0 increases, the number of sources in every bin increases.After sampling, this translates to the loudest randomly realized sources having higher masses.
The background sees a similar increase in the mass of its dominating sources.In addition to this, it also has a larger number of contributing sources in every mass and frequency bin.Thus, h c,BG increases near-uniformly across all frequencies.Its amplification matches that of h c,SS at high frequencies and exceeds that of h c,SS at low frequencies.There, SMBHB numbers are the largest because sources harden more slowly at low frequencies.This is particularly true of high-mass sources, whose numbers dwindle at high frequencies where they harden quickly by emitting more GWs.Thus, scaling up the number of sources in every bin leads to many more massive sources contributing to the lowfrequency GWB.
The GSMF characteristic mass m ψ,0 (Fig. 1, right-hand side), sets where the expected number of binaries drops off.Thus, varying m ψ,0 only significantly impacts mass bins corresponding to the lowest end of our varying m ψ,0 range (m ψ,0 = 10.5) and above.Thereby, Fig. 1 shows that increasing m ψ,0 from 10.5 to 12.5 raises the M SS 68% CI most dramatically (∼ 1.6 dex) at low frequencies, where massive sources are more common, and more moderately (∼ 0.5 dex) at high frequencies.
The background also sees large increases in ⟨M⟩ BG at low frequencies because many loud binaries remain there, even after the loudest has been extracted, and little change in typical mass at high frequencies where few high-mass sources remain after single-source extraction.This amounts to h c,SS having an overall more significant increase than h c,BG .
While high mass sources dwindle at high frequencies due to fast GW hardening, lower mass sources are driven more quickly through the low frequency regime (for fixed binary lifetime) by environmental processes.Thus, a lower m ψ,0 flattens both h c,SS and h c,BG at low PTA-band frequencies because the loss of the most massive sources is exacerbated by the environmentally driven hardening of lower mass sources.The effects of environmental hardening (represented by ν inner ) are explored further in §3.1.3,but here we find that when the characteristic strain is dominated by less massive binaries, environmental hardening flattens the low-frequency spectrum more.
Regardless of the model, low-frequency sources tend to be nearer.This is because more massive sources (which are numerous at low frequencies) take longer to evolve to the PTA band, as explained in greater detail in §3.1.3.Longer evolution times let these massive sources reach the PTA band at smaller redshifts and closer distances.

M BH -M bulge Relation
Fig. 2 shows h c , M, and d c for varying M BH -M bulge parameters.Given a population of host galaxies, the M BH -M bulge relation in Eq. 4 sets the masses of these galaxies' central SMBHBs.Increasing the relation's mass normalization (µ) shifts all SMBHBs to higher masses.This has negligible impact on low mass bins, where the number density is large and changes gradually with mass.However, at high masses, the number density drops off quickly with increasing mass.This means that shifting the expected number of sources of one bin to the next bin of increasing mass significantly raises the chances of realizing a source in that higher mass bin.Overall, this increases the odds of randomly sampling a source in any high-mass bin.
The left column of Fig. 2 shows that a ∼1 dex increase in M SS follows the increase in µ from 7.6 to 9.0 at low frequencies, where massive sources are numerous.The ∼1 dex increase in mass is not as great as the ∼1.6 dex increase in mass normalization because these higher mass binaries also evolve faster.At high frequencies, the 68% CIs see a more modest increase, which follows from the fact that these represent lower mass sources.The changes in h c,SS are approximately proportional to M 5/3 SS (see Eqs. 11 and 12) with deviations below this relation owing to unequal mass ratios.Meanwhile, at high frequencies, the background is minimally affected by µ due to the lack of massive sources, especially after the loud-est have been removed.Ultimately, across all frequencies, increasing µ raises h c,SS more than h c,BG .
Increasing scatter ϵ µ in the M BH -M bulge relation preferentially scatters sources to higher mass bins through Eddington Bias.Like µ, this increases the odds of sampling sources from the highest mass bins.The right column of Fig. 2 shows that an increase in ϵ µ from 0.0 to 0.9 increases all low-frequency h c , h c,SS slightly more than h c,BG .This preferential scattering effect becomes negligible at lower masses where the mass function flattens and number densities become large; thus, we see the M SS 68% CI and ⟨M⟩ BG medians both converge at high frequencies where lower mass sources dominate.Similarly to the m ψ,0 scenario, anytime low-mass sources are a greater contributor to the overall characteristic strain spectrum, the low-frequency end of the spectrum is also flatter, due to low-mass sources' hardening being sped up by environmental processes.
In addition to the preferential scattering systematically increasing masses, introducing scatter to the M BH -M bulge relation adds a second element of randomness beyond the Poisson bin sampling, widening the variance of random realizations.This is apparent as a slight widening of the 68% CI M SS and h c,SS regions and is more obvious in the SS 95th percentiles, which generally increase by ≳ 0.5 dex, even at high frequencies where the 68% CIs converge.
The M BH -M bulge relation slightly impacts the distribution of d c,SS at low frequencies, with the 68% CI including more distant sources when ϵ µ is higher.When scatter is low and arbitrarily high mass sources are less likely to be sampled, nearby sources tend to be the loudest.When the scatter is increased, the loudest source could instead be more massive, but further away.This effect occurs only at low frequencies because that is where ϵ µ causes a measurable change in M SS and h c,SS .We also see a frequency dependence of d c,SS matching that described in the last section: low-frequency single sources tend to be nearer because these sources are more massive, and thus take longer to reach the PTA band.

Hardening
In our phenomenological hardening model, all sources take the same total time from galaxy merger to SMBH coalescence, set by τ f .The way this hardening rate is distributed across the binary's lifetime depends on its mass, such that higher mass sources spend less time at high frequencies.This is because 1) more massive sources evolve more rapidly by GWs when they reach small separations and 2) more massive sources merge at larger separations.These combined effects require high-mass sources to evolve slower than low-mass sources at low frequencies, in order to meet the same fixed τ f .Thus, in our model, massive sources take longer to reach the PTA frequency band.Lower-mass sources will reach the PTA band more quickly, and then dwell there longer as they harden slower than their high-mass counterparts until they eventually coalesce.The slower evolution of massive binaries to PTA frequencies means they tend to be nearer when they emit in the PTA band.However, those that start at too low of redshifts will not reach separations small enough for PTA-band emissions by redshift zero.
When τ f is increased, this extends the evolution time of binaries at all masses, moving the entire populations to smaller final comoving distances.The left column of Fig. 3 shows this to be similarly true for both the background and single sources, with ⟨d c ⟩ BG decreasing by ∼ 1.0 dex and d c,SS by ∼1.1 dex when τ f is raised from 0.1 Gyr (light blue) to 11.0 Gyr (dark blue).Secondly, when the hardening time is extended, the most massive sources are unlikely to reach small enough separations to emit at PTA-detectable frequencies, producing a decrease in mass.This effect is largest at high frequencies, to which binaries take the longest to evolve.Thirdly, when the overall binary lifetime is short, even high mass binaries can be driven by non-GW hardening represented by ν inner , causing a low-frequency flattening or turnover in the h c spectrum as described in the ν inner analysis below.The characteristic strain has a greater proportional dependence on M SS than d c,SS , so the changes in h c,SS with varying τ f follow the changes in M SS .When considering the background, the filtration of massive sources applies to a larger number of sources at low frequencies, where massive sources are numerous.Thus, long τ f decreases h c,BG slightly more than it decreases h c,SS at low frequencies.
Per Eq. 7, ν inner sets the hardening rate in the smallseparation regime, with the asymptotic behavior of A flatter (less negative) ν inner increases the hardening rate at the lowest end of the PTA frequency regime before (da/dt)| gw dominates (see Fig. 3 in Agazie et al. 2023b).This represents faster hardening by processes like stellar scattering and circumbinary disk torques and produces the low-frequency turnover in both h c,BG and h c,SS apparent in the top panels of Fig. 3. Single sources are only impacted at low frequencies because there, when ν inner = 0, even the most massive sources will have phenomenological hardening dominate GW hardening.The same is true of massive binaries when the total binary lifetime is short.Lower mass binaries can be dominated by phenomenological hardening up to higher frequencies as their GW emission is weaker.Thus, the backgroundincluding contributions from lower mass binaries-sees a lower h c,BG across all frequencies for flat ν inner .
The bottom right panel of Fig. 3 shows the 68% CI of lowfrequency single sources to be more distant for flat ν inner than either of the other two cases.We attribute this to a selection bias: when ν inner is flat, massive sources evolve through the low end of the PTA band quickly.Thus, there are fewer of them.While any individual source is just as able and likely to reach small distances, having fewer of them decreases the odds of having an especially close source.These trends only appear in the 68% CI because the lower bounds of d c,SS for ν inner represent the less-common cases where one of the few loud sources just so happens to be nearby.This source can be as near under flat ν inner as it can for steep ν inner ; it is just less likely to exist in the first place.

CW Detection Occurrence Rate
Fig. 4 shows the single source detection probability (DP SS ) and the background detection probability (DP BG ) as a function of each varying parameter, for a fixed-PTA configuration.This 'fixed-PTA' method involves calibrating the PTA's noise level so that the median h c,BG across all realizations of the mean-parameter model (i.e. with the fiducial parameter values listed in Table 1) has a 50% DP BG .For example, the top left panel represents varying ψ 0 , so the PTA is calibrated to the median h c,BG across 500 realizations of the ψ 0 = −2.5 model.This PTA is used throughout the rest of the varying ψ 0 with phenomenological hardening analysis.The resulting DP SS medians are represented by a solid green line, as well as 50% and 95% DP SS CI as green-shaded regions.The resulting DP BG medians are represented by a dashed darker green  1).The PTA is calibrated once for each panel's mean phenomenological model, and once for each panel's mean GW-only model.The single source detection probability is shown in color, for varying GSMF parameters in green, varying M BH -M bulge parameters in orange, and varying hardening parameters in blue.The DP SS medians are represented by solid lines and the 68% and 95% CIs are shaded.The background detection probability (DP BG ) is given in darker shades of the same colors, with medians as dashed lines and 68% and 95% CI shaded.
line, with 50% and 95% CI as darker green shaded regions.The rest of the panels show the same, but for varying m ψ,0 , µ, ϵ µ , τ f , and ν inner as labeled.
In all cases, the DP SS medians remain below DP BG , consistent with the expectation for GWB detection to occur before CW detection (e.g.Rosado et al. 2015).DP BG is remarkably well constrained (95% CI spanning ≲ 0.5 dex), while DP SS 95% CI often range all the way from ∼10 −3 to ∼10 0 .These 95% CIs of DP SS can exceed DP BG in a few corners of parameter space, most notably for ψ 0 ≲ −2.3, τ f ≳ 5 Gyr, and ν inner ≳ −0.75.Thus, low m ψ,0 , long τ f , and flat ν inner are disfavored.Recall from §2.3 that calculating DP SS by integrating over the F e -statistic with zero signal present still produces a nonzero detection probability equal to the FAP, hence the floor of DP SS ≥ 10 −3 .
Although the variance between DP SS realizations is large, there are clear trends in how both DP SS and DP BG depend on the model parameters.For m ψ,0 , µ, and ϵ µ , the DP SS medians behave similarly to the DP BG medians, just at lower values.The greatest difference in DP SS and DP BG behavior occurs for the hardening parameters, both of which decrease DP BG by 3 dex, but only decrease DP SS by ≲0.7 dex.ψ 0 also shows significantly less impact on h c,SS than h c,BG , the h c,SS medians increasing only by ∼1.3 dex compared to DP BG increasing by ∼3 dex.
In Fig. 4, there are a wide range of DP BG s.However, following recent PTA results, there is considerable evidence for a GWB signal.Thus, in Fig. 5 we calibrate a PTA independently for each realization of each set of parameters.This shows how single source detection depends on each parameter, for a fixed confidence in the GWB.Because this 'realization-calibrated' method is informed by current evidence for the GWB and allows for a more nuanced exploration of parameter space, we present this as our fiducial method for identifying where in parameter space single sources are most/least likely to be detected, with the key results being those of Fig. 5.Meanwhile, Fig. 4 is useful to distinguish effects due to background calibration from direct effects on single sources.
Fig. 5 also includes the GW-only model, which use similar 'realization-calibrated' PTAs, with the resulting ⟨N SS ⟩ medians in dash-dotted light grey and 68% CI shaded in light grey.Note that the τ f and ν inner panels include constant GWonly data because the GW-only model has no τ f or ν inner to vary.The rest of the GW-only results follow the same trends as the phenomenological cases but are lower by up to 1 dex.This effect is similar to having a very steep ν inner : both involve GW hardening dominating the entire PTA band.Without phenomenological hardening speeding up the evolution, sources dwell in the PTA band longer, increasing the number of binaries contributing to the total h c,BG , while the individual loudest remain unaffected (see §3.1.3).Since DP BG is calibrated to 50%, we see the changes in ⟨N SS ⟩.
Given that PTA's have not yet detected a CW, very long τ f and flat ν inner , both of which predict ⟨N SS ⟩ ∼ 1, are unlikely.This is independent, but consistent with the short τ f constrained by GWB data in Agazie et al. (2023b) and Antoniadis et al. (2023a).Agazie et al. (2023b) also favors flatter ν inner , with −0.4 as their maximum-likelihood posterior.There are clear trends in the medians and 68% CI of ⟨N SS ⟩ for all other model parameters, as well.Thus, CW detections can inform and constrain our astrophysical models for SMBHB populations and evolution, beyond the constraints placed by measuring the GWB amplitude.
GSMF -The top left panel of Fig. 5 shows that ⟨N SS ⟩ decreases smoothly from ψ 0 = −3.5 to ψ 0 = −1.5, as a result of the more significant increases in h c,BG than h c,SS for increasing GSMF normalization.Therefore, a larger total number of galaxies in the universe increases the likelihood of any GW detection, but disfavors CW detection for a fixed DP BG .The top right panel of Fig. 5 shows ⟨N SS ⟩ increases with m ψ,0 , indicating that the single source detectability increases more with GSMF characteristic mass than background detectability does.
M BH -M bulge Relation - §3.1.2and Fig. 2 show that increasing µ raises h c,SS slightly more than it raises h c,BG .Thus, the middle left panel of Fig. 5 shows a subtle positive trend in ⟨N SS ⟩ versus µ.In Fig. 4 it is clear that ϵ µ affects DP BG the least of the six parameters because the impact of scatter on h c,BG is minor.Rather, ϵ µ primarily affects the edge cases of the loudest/most massive sampled sources, i.e., those extracted as single sources.Thus, there is dramatic growth in DP SS and moderate growth in ⟨N SS ⟩ when raising ϵ µ from 0.0 to 0.9.
Hardening -Fig. 4 shows DP BG decreases with τ f , except at short τ f where the loss of massive binaries due to failure to reach the PTA band becomes less significant (until a slight uptick occurs at 0.1 Gyr where decreasing distance  1).The different colors correspond to the same models as in Figs. 1, 2, and 3, where green, orange, and blue represent the single sources for the GSMF, M BH -M bulge relation, and hardening parameters respectively, and shades of grey represent the h c,BG -weighted average values.The 10 loudest single sources at each frequency are used for the DP SS -weighted number densities, and all but these 10 loudest are used for the DP BG -weighted number densities outweighs the loss of massive sources).Single sources become nearer, but are less affected by the filtration of binaries, leading ⟨N SS ⟩ in the bottom left panel of Fig. 5 to increase to ≳1 at long τ f .In the bottom right panel of Fig. 5, ⟨N SS ⟩ increases most dramatically between −0.75 ≲ ν inner ≲ 0. Thus, CW detections (assuming fixed GWB confidence) become increasingly likely when phenomenological hardening processes speed up small separation (a ≲ a c = 10 2 pc) binary evolution.Any single source detection would make it highly unlikely that ν inner is steep.

Frequency
The comoving distance d c versus mass M of the most detectable single sources are presented in Fig. 6.Contour lines show the 0.5, 1.0, 1.5, and 2.0 σ contours of the SNRweighted number of sources in each bin of M and d c , with the same panel and color convention for the six parameters as in Fig. 5.The light, medium, and dark shades represent each model parameter's lowest, mean, and highest values, respectively.Solid colored lines represent the single sources, while the h c,BG -weighted M and d c are shown in dashed grey contours for comparison.
Across all parameters, the 2σ region spans masses of 2 × 10 8 M ⊙ ≲ M SS ≲ 1.5 × 10 10 M ⊙ , while the meanparameters model (with the Table 1 fiducial values) peaks at ∼ 8 × 10 9 M ⊙ .These results are consistent with the chirp masses Rosado et al. (2015) predict a 20-year IPTA could most likely detect, all peaking in probability between 10 9.5 M ⊙ and 10 10 M ⊙ .Long hardening times significantly decrease the peak masses, down to ∼ 10 9 M ⊙ by ruling out nearby massive sources that would not have time to reach PTA-band separations.
The mean-parameters model predicts single-source d c peaking in SNR-weighted number at ∼800 Mpc, with a 1.5σ region spanning 70 Mpc ≲ d c,SS ≲ 2100 Mpc.Only a small region of this parameter space overlaps with IPTA redshift predictions by Rosado et al. (2015), which are mostly beyond ∼ 2100 Mpc (z ≳ 0.56).These distances are only impacted by the GSMF and M BH -M bulge relation parameters insofar as having more massive binaries can allow more distant sources to be the most detectable, as shown by the correlation between increasing mass and greater distance in Fig. 6.Yet, the background distances are not affected, so single sources break the degeneracy in distances for different GSMF and M BH -M bulge parameters.The same is true of masses for varying m ψ,0 , µ, and ϵ µ .This shows how properties of single sources can break degeneracies between our model components in a way that is not possible using GWB detection alone.
The hardening parameters, on the other hand, impact both single source and background distances by setting the time it takes binaries to reach the PTA band, as described in §3.1.3.This is most obvious for the hardening time, where longer τ f leaves the loudest sources at peak distances as close as ∼ 250 Mpc by the time they are emitting nHz GWs.Similarly, when ν inner is very flat, the hardening timescales at separations just before the PTA band are small, meaning single sources reach these frequencies more quickly and thus at larger distances, with the SNR-weighted number peaking at ∼1500 Mpc.
The DP-weighted average CW frequency across all single sources in all realizations ⟨ f SS ⟩ is presented as a function of each varying parameter in Fig. 7.The color and panel convention follows that of Fig. 5 with the GW-only model again represented by dash-dotted grey lines.Shaded regions represent one standard deviation above and below the median in log space, calculated according to Eq. 20.In all phenomenological cases regardless of parameters, the CW frequency most likely to be detected by our 16.03-yearPTA is For this 16.03-yearPTA, ⟨ f SS ⟩ is generally in the lowest frequency bin, in agreement with similar PTA duration predictions in Rosado et al. (2015), because white-noise-only PTA models give a monotonic decrease in DP SS versus frequency.If h c,SS continues to increase with decreasing frequency, the loudest sources will likely remain in the lowest frequency bin.However, h c,SS may instead plateau at low frequencies, moving the average detection frequency to a specific value where the SS strains are maximized relative to the combined noise of the PTA and GWB (Kelley et al. 2018).Including red noise decreases the detection probability of the lowest frequency sources, thus moving the ⟨ f SS ⟩ to higher frequencies.Given that pulsars typically have some intrinsic red noise (Agazie et al. 2023e), the white-noise-only ⟨ f SS ⟩ predictions should be treated as lower limits on the predicted frequency of first CW detection.We explore the effects of varying red noise models on these predictions in §A.
The GW-only model mostly predicts similar frequencies to the phenomenological models, but with some increase up to ∼0.1 yr −1 at high ψ 0 and low m ψ,0 , µ, and ϵ µ -everywhere ⟨N SS ⟩ is low.When ⟨N SS ⟩ is low, the background likely becomes a more significant source of red noise, pushing the most detectable sources to higher frequencies.This also allows for more variation in the highest-DP frequency between realizations, hence the larger weighted standard deviation in the low ⟨N SS ⟩ regions of parameter space.

Anisotropy in the Gravitational Wave Background
We calculate C ℓ up to ℓ max = 8 for each of the models presented in Fig. 6.The resulting C ℓ s are indistinguishable for each ℓ ≥ 1 of a given model, consistent with similar holodeck predictions marginalizing over many semianalytic models in Agazie et al. (2023d), predictions using cosmological dark matter simulations to assign SMBHBs (Taylor & Gair 2013), and the analytic shot noise approximation for GWB anisotropy in Sato-Polito & Kamionkowski (2023).Thus, we present C 1 /C 0 versus frequency as a proxy for any C ℓ /C 0 in Fig. 8, including results from the lowest, mean, and highest variation of each of the six model parameters.These are calculated using the 1000 loudest sources (solid line medians and shaded 68% CI) and 10 loudest sources (dotted line medians) in each frequency bin.Remarkably, using the 10 loudest and the 1000 loudest sources give C 1 /C 0 medians and standard deviations both within ≤ 0.16 dex of each other at any frequency, with average differences of just ∼0.03 dex and ∼0.02 dex, respectively.Thus, in our models, anisotropy in the GWB is determined by ≤ 10 loudest sources in any given frequency bin.
The medians span 6 × 10 −3 ≲ C ℓ /C 0 ≲ 2 × 10 −1 at low frequencies and 2 × 10 −1 ≲ C ℓ /C 0 ≲ 6 × 10 −1 at high frequencies.This increase in anisotropy with increasing frequency is expected because dwindling numbers of massive sources make the background h c,BG drop off more than individual sources' h c,SS , until there is hardly even a "background" at high frequencies.These results overlap with the C ℓ /C 0 of the Sato-Polito & Kamionkowski (2023) models favoring a small number of binaries with high mass and low redshift.However, the remaining three of their five models, which reproduce the same characteristic strain using a larger number of lower mass, higher redshift binaries, predict levels of anisotropy up to five orders of magnitude below ours, at low frequencies.Meanwhile, most of our confidence in-tervals overlap with the Taylor & Gair (2013) calculation of C ℓ>0 /C 0 ∼ 0.1 for a dark matter simulation-based SMBHB population, marginalized over many frequencies.Agazie et al. (2023d) place Bayesian upper limits of C 1 /C 0 ≲ 2 × 10 −1 (circles with dashed lines in Fig. 8).Most 68% CI overlap or nearly reach these upper limits, suggesting that if the GWB is produced by SMBHBs, anisotropy may be detected in the near future, and the lack thereof could place stringent constraints on our parameter space.In fact, very long hardening time and flat hardening index predict median anisotropy levels above the current upper limits.This disfavors these corners of parameter space and supports the idea that single sources are particularly useful for constraining binary evolution.
Fig. 9 shows C 1 /C 0 for the lowest five frequency bins, as a function of each varying parameter, using just ten loudest per frequency bin.In comparing Fig. 9 to Fig. 5, it is evident that the models with the greatest C ℓ /C 0 correspond to those with the highest ⟨N SS ⟩.This is because increasing ⟨N SS ⟩ and increasing anisotropy both stem from cases where the loudest single sources become more dominant.The greatest model-dependent changes in anisotropy are for long τ f and flat ν inner .Both these scenarios produce very high C ℓ /C 0 (≳ 0.2) at the lowest end of the PTA band, and then are nearly constant with frequency.The significant increases in C ℓ /C 0 at low frequencies correspond to the scenarios in Fig. 3 where h c,BG decreases significantly and h c,SS sees less change.

DISCUSSION
We present the dependence of single-source detection statistics and anisotropy on astrophysical model parameters.These models include several assumptions to keep in mind.First, we assume circular orbits for all binaries.Allowing for eccentricity may move some GW energy from lower to higher frequencies, and could have different impacts on the loudest single sources versus the background -a subject worth further investigation (Siwek et al. in prep.).Another caveat is that our hardening model prescribes a fixed hardening time for all binaries, regardless of mass or redshift.This is a useful approximation to self-consistently examine how changing overall binary evolution impacts GW models without adding too many degrees of freedom, but there is no reason that these hardening times should not be mass-dependent.Thus, we suggest allowing binary lifetimes to depend on mass as a potential way to expand upon this hardening model.A third caveat to this model is that the SMBH-host galaxy relations use empirical measurements of local galaxies.These relations can be improved as more EM data is collected, particularly about more distant galaxies.The relations are based entirely on bulge mass and could be expanded by including velocity dispersion (Matt et al. 2023;Simon 2023).
Degeneracies and covariances between model components serve as key challenges in making GW-based constraints.This work shows the potential for using single sources to distinguish between parameters that are degenerate in shaping the GWB.By considering multiple parameters, one can raise single source detection probability while maintaining the GWB amplitude by increasing the masses of binarieswhether through the GSMF m ψ,0 or M BH -M bulge µ-and decreasing the total number of binaries.The previous examples translate similarly to anisotropy.Fig. 6 demonstrates that the mass and distance of loud single sources can further break degeneracies, as discussed in §3.2.2.While we defer an exploration of specific degeneracies to future work incorporating PTA data, the discrepancy in how ⟨N SS ⟩ and anisotropy versus the GWB amplitude vary implicitly shows that single sources can carry different information from the GWB.Thus, we can make the strongest constraints by combining them.
Another challenge in making any conclusions based on single source detection statistics is the 3 dex spread of ⟨N SS ⟩ 95% CIs.This is a result of the fact that CW detections depend upon the random chance of a particularly massive binary happening to be nearby.This randomness limits the precision of any single source predictions or parametric constraints by semi-analytic models.Incorporation of galaxy catalogs may allow for more narrowly constrained predictions as to what single sources could be realized in our universe.
EM data may also inform the level of GWB anisotropy in our universe.By placing the sources randomly and treating the background as purely isotropic, we make conservative estimates for anisotropy.However, one might predict that more SMBHBs will be emitting PTA-band GWs in regions of higher cosmic density.A future step would be to study possible correlations between GWB anisotropy and galaxy clustering or large-scale structure.On the other hand, if just a few loudest sources at distances of ∼ 1000 Mpc entirely determine anisotropy (as we have found), then these correlations seem less likely because the placement of individual sources is random on scales large enough to treat the universe as purely isotropic.Regardless, this would be an interesting hypothesis to test.
By varying PTA noise to make each model produce 50% DP BG , we comprehensively explore the detection statistics of a wide parameter space, including models that produce low GWB amplitudes.The next step to build on this parameter space exploration is to condition our models on current measurements of the GWB amplitude.With these GWB-conditioned models, we can fully explore the multidimensional parameter space to explicitly include covariances and use realistic PTAs to calculate our detection statistics, as opposed to calibrating to a fixed DP BG .The resulting background detection statistics will serve as a check on the GWB constraints set by (Agazie et al. 2023b).Then, we will test whether the current lack of CW detections (Agazie et al. 2023c) and upper limits on anisotropy (Agazie et al. 2023d) can further constrain these model parameters.We expect that long τ f and flat ν inner will be most easily constrained by single source detection statistics and anisotropy, given that these are the regions of parameter space where both ⟨N SS ⟩ and C ℓ /C 0 are highest and ⟨N SS ⟩ has the lowest variance.

CONCLUSIONS
In this work, we develop an approach for modeling CWs distinguishable from a background of SMBHBs, their sources' properties, and their corresponding GWB anisotropy.We develop a detection statistics pipeline that calibrates a simulated PTA to a 50% probability of detecting the background and calculates the expected number of single-source detections under those settings.Our primary conclusions are the following: 1. GW anisotropy and CW detections (or lack thereof) convey specific information about the astrophysics governing galaxy population, their SMBHBs, and binary evolution.This anisotropy and CW information can break model degeneracies and allow for much more stringent constraints than possible with the GWB alone.
2. CWs are increasingly likely to be observed for low GSMF normalization ψ 0 , high GSMF characteristic mass m ψ,0 , high M BH -M bulge mass normalization µ and intrinsic scatter ϵ µ , long hardening time τ f , and flat small-separation hardening index ν inner .
3. Anisotropy in the GWB, represented by the angular power spectrum C ℓ over multipole modes ℓ, is determined in our models by ≤ 10 loudest sources at each frequency and is the same for any 1 ≤ ℓ ≤ ℓ max .Models vary in their low-frequency predictions for C ℓ>0 /C 0 from ∼5 × 10 −3 to ∼2 × 10 −1 and converge to ∼3 × 10 −1 at high frequencies.
4. Models with greater single-source detection probability tend to have higher anisotropy, often exceeding current upper limits of C 1 /C 0 ≲ 0.2 at low frequencies.Thus, not detecting anisotropy could strongly constrain our parameter space.
5. Of all our model components, binary evolution shows the greatest promise for constraints using single sources and anisotropy.Long hardening time and flat ν inner give the greatest probability of CW detection for a fixed GWB confidence and high anisotropy even at low frequencies.
6.The most detectable single sources are found in the lowest frequency bin for a 16.03-year PTA with masses ranging from ∼10 9 M ⊙ to ∼3 × 10 10 M ⊙ and final comoving distances ranging from ∼250 Mpc to ∼2500 Mpc.The most detectable frequency has little dependence on the model but increases with greater pulsar red noise.
Only the hardening parameters have a demonstrable impact on these sources' final comoving distances, with longer τ f and steep ν inner resulting in the closest sources.
6. ACKNOWLEDGMENTS The authors thank Jeff Hazboun, Jeremy Baier, Bence Bécsy, and Neil Cornish for helpful discussions throughout the development of this paper, particularly regarding CW detection statistics.We also thank Nihan Pol for insight into the anisotropy methods and interpretation.Finally, we thank the anonymous referee for thoughtful and constructive feedback.The primary work in this paper uses a white-noise-only simulated PTA.We also consider the impact of red noise models, parameterized by an amplitude A RN at reference frequency f ref = 1 yr −1 and power-law index γ RN , on our detection probabilities and DP SS -weighted average frequencies.The red-noise PTAs are calibrated by fixing γ red and the ratio Q of red noise to white noise while allowing the total noise amplitude to vary.Fig. 10 shows the resulting ⟨N SS ⟩ as a function of the six varying model parameters for the fiducial white noise only model (black), red-noise with spectral index γ RN = −1.5 (purple), and red-noise with spectral index γ RN = −3.0(red).For the red-noise models we include ratios of Q = 0.01 (dash-dotted), Q = 1.0 (solid), and Q = 100 (dashed).The addition of red noise generally raises the relative single source detection probability, because it makes the GWB less distinguishable from the noise, such that the total noise calibration must be lower for a 50% DP BG .In fact, when steep (γ RN = −3.0)red noise dominates (Q = 100), the median ⟨N SS ⟩ is ≳ 0.1 for all model variations except the highest ψ 0 s.The increase in ⟨N SS ⟩ from previously low regions results in an overall flattening in parametric dependence, but maintains the sign of the derivative, i.e. whether ⟨N SS ⟩ is increasing or decreasing with each parameter.
Fig. 11 shows that adding red noise also raises the most detectable CW frequency because lower frequency sources  are drowned out.Moderate red noise (γ RN = −1.5, purple) maintains the lack of dependence on model parameters seen in the white noise cases.However, adding steep red noise (γ RN = −3.0,red) with a ratio of Q ≳ 1 creates a dependence of ⟨ f SS ⟩ on each varying parameter.For all but ν inner , this ⟨ f SS ⟩ dependence tends to follow the opposite trend of ⟨N SS ⟩.When ⟨N SS ⟩ is low, there is greater noise from the PTA calibration, especially at low frequencies, pushing ⟨ f SS ⟩ higher.

Figure 1 .
Figure1.Characteristic strain (top row), total mass (middle row), and final comoving distance (bottom row), for varying ψ 0 (left column) and m ψ,0 (right column).Shaded green regions represent 68% CI of the single loudest source at each frequency, with markers to indicate the 95th percentiles for h c,SS and M SS and 5th percentiles for d c,SS .Dashed lines represent the median background (all but the single loudest source per frequency) characteristic strain h c,BG and the h c -weighted background properties (⟨M⟩ BG and ⟨d c ⟩ BG ).m ψ,0 increases from -3.5 to -2.5 to -1.5 and ψ 0 increases from 10.5 to 11.5 to 12.5 for darkening shades of green/grey.

Figure 2 .
Figure2.Same as Fig.1but for the M BH -M bulge parameters: increasing µ (left column) from 7.60 to 8.30 to 9.00 and increasing ϵ µ (right column) from 0.00 to 0.45 to 0.90 for darkening shades of orange (single sources) and grey (background).

Figure 3 .
Figure 3. Same as Fig. 1 but for the phenomenological hardening parameters: increasing τ f (left column) from 0.10 Gyr to 5.55 Gyr to 11.00 Gyr and flattening ν inner (right column) from -1.50 to -0.75 to 0.00 for darkening shades of blue (single sources) and grey (background).

Figure 4 .
Figure 4. Detection probability for a PTA calibrated to the median DP BG of the mean parameter model (i.e. using the fiducial values in Table1).The PTA is calibrated once for each panel's mean phenomenological model, and once for each panel's mean GW-only model.The single source detection probability is shown in color, for varying GSMF parameters in green, varying M BH -M bulge parameters in orange, and varying hardening parameters in blue.The DP SS medians are represented by solid lines and the 68% and 95% CIs are shaded.The background detection probability (DP BG ) is given in darker shades of the same colors, with medians as dashed lines and 68% and 95% CI shaded.

Figure 5 .
Figure 5. Expected number of single-source detections ⟨N SS ⟩ for a white-noise-only PTA calibrated independently to a 50% background detection probability for each parameter and realization.⟨N SS ⟩ is given as a function of varying GSMF parameters (top row) in green, M BH -M bulge parameters (middle row) in orange, and phenomenological hardening parameters (bottom row in blue) for the phenomenological hardening model.The medians are solid lines and the 68% and 95% CI are shaded.⟨N SS ⟩ for the GW-only hardening model has medians represented by dash-dotted lines and 68% CI shaded in grey.These are replaced in the bottom row by constant values corresponding to the fiducial GSMF and M BH -M bulge model because with GW-only hardening there are no τ f and ν inner to vary.

Figure 6 .
Figure6.SNR-weighted number density of single sources' final comoving distance versus total mass.The contours represent 0.5, 1.0, 1.5, and 2.0 σ contours for 3 variations of a single parameter, while all other parameters are fixed at their mean values.The middle shade in each plot refers to the mean-parameters model (i.e. with the fiducial values in Table1).The different colors correspond to the same models as in Figs.1, 2, and 3, where green, orange, and blue represent the single sources for the GSMF, M BH -M bulge relation, and hardening parameters respectively, and shades of grey represent the h c,BG -weighted average values.The 10 loudest single sources at each frequency are used for the DP SS -weighted number densities, and all but these 10 loudest are used for the DP BG -weighted number densities

Figure 7 .
Figure 7. DP-weighted frequency of the loudest single sources, as a function of each varying parameter while the rest of the parameters are fixed.Colored regions and solid lines represent the 1σ regions and means for the phenomenological hardening model, while grey regions and grey dash-dotted lines represent the same for the GWonly hardening model.

Figure 8 .
Figure 8. Anisotropy in terms of C ℓ /C 0 for the first spherical harmonic mode as a function of frequency, for varying astrophysical models.Medians (solid lines) and 68% CI (shaded) correspond to the model of the same panel and color in Fig. 6, as labeled.By these methods, any C ℓ>0 /C 0 is indistinguishable up to ℓ max = 8, so the C 1 /C 0 data plotted represents any C ℓ>0 /C 0 distribution.Bayesian upper limits on C 1 /C 0 from Agazie et al. (2023d) are plotted for comparison as purple circles.

Figure 9 .
Figure 9. C ℓ>0 /C 0 medians and 68% CI for the first five frequency bins (1.98, 3.95, 5.93, 7.91, and 9.88 nHz) as a function of each varying model parameter.The data plotted uses ℓ = 1, but is indistinguishable from any other 1 ≤ ℓ ≤ ℓ max .The panels correspond to the same parameters as in Fig. 4, Fig. 5, and Fig. 7.These C ℓ /C 0 s are calculated using just 10 loudest sources in each frequency bin, which Fig. 8 shows sufficiently reproduces the anisotropy in our model calculated using 1000 loudest sources.The Agazie et al. (2023d) upper limit on C 1 /C 0 in the lowest frequency bin is denoted by a horizontal dashed blue line.

Figure 10 .
Figure 10.Expected number of single-source detections when the background detection probability is calibrated to 50% by varying the pulsar noise levels.The noise is set by a fixed ratio Q = S RN ( f ref )/S WN ( f ref ) between white noise and red noise at reference frequency f ref = 1 yr −1 , and fixed red noise index of γ RN = −1.5 (purple), or −3.0 (red).Noise ratios of 0.01, 1.0, and 100.0 are represented by dash-dotted, solid, and dashed lines, respectively, and the 68% CI of the white noise-only results are shaded.

Figure 11 .
Figure 11.DP SS,i -weighted frequency of the loudest single sources as a function of varying model parameters.The color and line style conventions follow those of Fig. 10, with shaded regions showing the weighted standard deviation of the white noise-only results.

Table 1 .
Astrophysical parameters of the model components investigated in this paper, while the rest remain fixed to the fiducial values in Agazie et al. (2023b) Table