## Nanohertz gravitational waves from a quasar-based supermassive black hole binary population model as dark sirens

Si-Ren Xiao, a Yue Shao, a Ling-Feng Wang, b Ji-Yu Song, a Lu Feng, c,a Jing-Fei Zhang a and Xin Zhang a,d,e, ∗

- a Liaoning Key Laboratory of Cosmology and Astrophysics, College of Sciences, Northeastern University, Shenyang 110819, China
- b School of Physics and Optoelectronic Engineering, Hainan University, Haikou 570228, China c College of Physical Science and Technology, Shenyang Normal University, Shenyang 110034,
- China
- d MOE Key Laboratory of Data Analytics and Optimization for Smart Industry, Northeastern University, Shenyang 110819, China
- e National Frontiers Science Center for Industrial Intelligence and Systems Optimization, Northeastern University, Shenyang 110819, China

E-mail: xiaosiren@stumail.neu.edu.cn, shaoyue@stumail.neu.edu.cn, wanglf@hainanu.edu.cn, songjiyu@stumail.neu.edu.cn, fengluu@foxmail.com, jfzhang@mail.neu.edu.cn, zhangxin@mail.neu.edu.cn

Abstract. Recently, several pulsar timing array (PTA) projects have detected evidence of the existence of a stochastic gravitational wave background (SGWB) in the nanohertz frequency band, providing confidence in detecting individual supermassive black hole binaries (SMBHBs) in the future. Nanohertz GWs emitted by inspiraling SMBHBs encode the luminosity distances of SMBHBs. They can serve as dark sirens to explore the cosmic expansion history via a statistical method to obtain the redshift information of GW sources' host galaxies using galaxy catalogs. The theoretical analysis of the dark siren method relies on the modeling of the population of SMBHBs. Using a population model consistent with the latest SGWB observations is essential, as the SGWB provides significant information about the distribution of SMBHBs. In this work, we employ a quasar-based model, which can self-consistently account for the SGWB amplitude, to estimate the population of SMBHBs. We constrain the Hubble constant using the mock GW data from different detection cases of PTAs in the future. Our results show that a PTA consisting of 100 pulsars with a white noise level of 20 ns could measure the Hubble constant with a precision close to 1% over a 10-year observation period, and a PTA with 200 pulsars may achieve this goal over a 5-year observation period. The results indicate that modeling the SMBHB population significantly influences the analysis of dark sirens, and SMBHB dark sirens have the potential to be developed as a valuable cosmological probe.

∗ Corresponding author.

## Contents

| 1   | Introduction            | Introduction                       |   1 |
|-----|-------------------------|------------------------------------|-----|
| 2   | Methodology             | Methodology                        |   3 |
|     | 2.1                     | Subcatalog of 2MASS galaxy catalog |   3 |
|     | 2.2                     | Population models of SMBHBs        |   4 |
|     |                         | 2.2.1 Major merger model           |   5 |
|     |                         | 2.2.2 Quasar-based model           |   5 |
|     | 2.3                     | GW detection and Fisher analysis   |   7 |
|     | 2.4                     | Inference of the Hubble constant   |  10 |
| 3   | Results and discussions | Results and discussions            |  13 |
| 4   | Conclusion              | Conclusion                         |  16 |

## 1 Introduction

The study of cosmology has entered the era of precision cosmology [1, 2]. Five of the six basic parameters of the ΛCDM model have been precisely constrained by the Planck cosmic microwave background (CMB) data with constraint precisions better than 1% [3]. However, as the precision of cosmological parameter measurements improves, some puzzling tensions have come to light. The Hubble constant ( H 0 ) is a crucial cosmological parameter describing the universe's current expansion rate. Its values inferred from the Planck CMB observation assuming the ΛCDM model [3] and the model-independent local-universe distance ladder measurement [4] are in more than 5 σ tension, which has become a critical crisis in cosmology [5, 6] and has been widely discussed (see, e.g. [5-35]). While the Hubble tension indicates the potential for new physics beyond the ΛCDM model, there remains a lack of consensus on a valid extended cosmological model to address this tension. Additionally, it is imperative to explore new cosmological probes capable of independently measuring H 0 .

Compact binary coalescence (CBC) events, a category of gravitational wave (GW) sources, have been discovered by the LIGO-Virgo-KAGRA collaboration in recent years, with approximately a hundred events detected. The luminosity distances of CBC events can be determined through the analysis of GW waveforms, and these events are commonly referred to as 'standard sirens' [36, 37]. Standard sirens can provide absolute distance measurements, avoiding poorly understood calibration processes inherent in the distance ladder method. If the redshifts of GW sources are also determined, the distance-redshift relation could be established to study the expansion history of the universe (see, e.g. [38-84]). The binary neutron star coalescence event, GW170817, gives H 0 = 70 +12 -8 km s -1 Mpc -1 with a precision of 14% [85]. Its redshift was determined by directly identifying the electromagnetic (EM) signals produced when the binary merged [86, 87]. Such standard sirens are known as bright sirens, with only one case detected to date.

For the standard sirens without EM counterparts (usually referred to as dark sirens), we can cross-match them with galaxy catalogs and statistically infer their redshifts (see, e.g. [88-104]). Recently, using 46 dark sirens with signal-to-noise ratios (SNRs) over 11 from GWTC-3 and the GLADE+ K -band galaxy catalog [105, 106], researchers achieved a ∼ 19%

constraint on H 0 . Combined with the bright siren GW170817, the constraint precision of H 0 improves to ∼ 10% [98].

The frequencies of GWs generated by the GW sources with varying masses are hugely different, which requires GW detectors sensitive to diverse frequency bands. The current LIGO-Virgo-KAGRA detector network [107-109] can detect only stellar-mass CBCs at around hundreds of hertz. In the future, precise measurements of GWs in various bands will become possible. In the frequency band from a few hertz to several thousand hertz, the third-generation ground-based GW detectors, e.g., the Cosmic Explorer [110] and the Einstein Telescope [111], will significantly enhance the detection capabilities for stellar-mass CBCs. In the decihertz band, the Decihertz Interferometer Gravitational-Wave Observatory [112-114] and lunar-based detectors [115, 116] could detect intermediate-mass black hole binaries and binary white dwarfs. In the millihertz band, GWs from massive black hole binaries and extreme mass ratio inspirals could be detected by space-borne GW detectors such as the Laser Interferometer Space Antenna [117-119], Taiji [120-123], and TianQin [124129]. In the nanohertz (nHz) band, the pulsar timing arrays (PTAs) could detect continuous GWs from the inspiraling supermassive black hole binaries (SMBHBs) resulting from galaxy mergers.

The pulsars constructing PTAs are millisecond pulsars (MSPs) and have highly stable rotation periods. Their radio signals can be periodically detected by radio telescopes on the Earth with extremely high precision. GWs passing between the MSPs and the Earth will alter the proper distances between them. This leads to a discrepancy between the actual and expected arriving times of radio signals, which is known as the timing residuals. By analyzing the time residuals of PTAs, we can extract information about nHz GWs. Recently, North American Nanohertz Observatory for Gravitational Waves (NANOGrav) [130, 131], European PTA (EPTA) + Indian PTA (InPTA) [132, 133], Parkes PTA (PPTA) [134, 135], and Chinese PTA (CPTA) [136] have reported substantial evidence for the existence of a stochastic GW background (SGWB) [137]. Such evidence boosts confidence in making more advances in detecting nHz GWs. For instance, detecting individual SMBHBs has become one of the next important missions of PTA observations [138-147]. Although extracting an individual GW signal from the superimposed GW background is challenging, this goal may be achievable with the addition of more advanced radio telescopes, especially in the era of the Square Kilometre Array (SKA).

Recently, Yan et al. [148] and Wang et al. [149] indicated that utilizing SMBHBs as standard sirens could help to improve the measurement precisions of cosmological parameters. Yan et al. [148] pointed out that the EM counterparts of SMBHBs can be observed during its inspiraling stage and forecast the constraint on the equation of state of dark energy with 154 SMBHB bright sirens observed by SKA-era PTAs. Wang et al. [149] investigated the dominant factors affecting the detection capability of PTAs for SMBHBs and found that the root mean square (RMS) of white noise, σ t , is more dominant than the number of MSPs in affecting the detections of SMBHBs. They also proposed that detecting nHz GWs from individual SMBHBs has significant implications for exploring the properties of dark energy and precisely determining the Hubble constant.

In the aforementioned papers, the SMBHB standard siren mock data are simulated based on the galactic major merger model proposed by Mingarelli et al. [150] (hereafter the M17 model). Recently, Casey-Clyde et al. [151] indicated that most major merger models predicted characteristic amplitudes of SGWB lower than that derived from the NANOGrav 12.5-yr data. In this context, they proposed a quasar-based SMBHB population model (here-

after the C21 model), considering that galaxy mergers could trigger quasars, and SMBHBs arising from galaxy mergers might be linked to quasar populations. This model can selfconsistently predict the local number density of SMBHBs given the amplitude of the SGWB.

Based on the C21 model and the common-process signal in the NANOGrav 12.5-yr dataset, the number of SMBHBs in the local universe is roughly five times higher than that predicted by the M17 model. It is evident that there is a significant difference in the number of SMBHBs predicted by these two models. Notably, the number of SMBHB dark sirens significantly influences the constraints on cosmological parameters, especially when the localization of dark sirens is poor [98]. Nevertheless, previous studies have not forecast the capability of SMBHB dark sirens of constraining cosmological parameters based on the C21 model. Consequently, a crucial question arises: Based on the C21 model that is more consistent with SGWB observations, what role can SMBHB dark sirens actually play in measuring cosmological parameters?

To answer this question, we simulate the mock SMBHB dark siren data based on the C21 quasar-based model and the SGWB signal inferred from the NANOGrav 15-yr dataset [130]. It needs to be emphasized that we employ the latest NANOGrav 15-yr dataset rather than the 12.5-yr dataset to ensure that our SMBHB population estimate is consistent with the most up-to-date PTA observations. We analyze the ability of SKA-era PTAs to detect the local-universe SMBHBs based on various PTA detection cases, and then combine the mock dark siren data with the 2 Micron All Sky Survey (2MASS) Extended Source Catalog [152] to perform cosmological parameter estimations. In this work, we focus only on the ability of SMBHB dark sirens to measure H 0 . This is because our aim is to investigate how SMBHB dark sirens can contribute to resolving the Hubble tension, and we primarily mock SMBHB dark sirens in the local universe, where the distance-redshift relationship is sensitive only to H 0 . Therefore, we use the SMBHBs dark sirens to infer H 0 and fix all other cosmological parameters to the Planck 2018 values.

This paper is organized as follows. Section 2 introduces the methodology used in this work. In Section 3, we report our constraint results and make relevant discussions. The conclusion is given in Section 4. Unless otherwise stated, we adopt the system of units in which G = c = 1 throughout this work.

## 2 Methodology

## 2.1 Subcatalog of 2MASS galaxy catalog

Massive galaxies are hypothesized to have a higher probability of hosting SMBHBs [153]. To select massive galaxies in the local universe, we first estimate their K -band absolute magnitudes by M K = m K - × 5 log 10 ( d L ) -25 -0 11 . × A v , where d L represents the luminosity distance to the galaxy. The NANOGrav Collaboration [154] considered the impact of peculiar velocities on redshift measurements, allowing for precisely calculating the distances of all galaxies in the 2MASS Redshift Survey (2MRS) catalog [155]. The apparent magnitude m K and the extragalactic extinction A v can be obtained from the 2MRS catalog. We implement a K -band cutoff with M K ≤ -25 for selecting galaxies with stellar mass M ⋆ ≥ 10 11 M ⊙ . Based on this criterion, we select 13172 galaxies from the 2MRS catalog as potential host galaxies for SMBHBs. These candidate host galaxies constitute a new catalog, referred to as the 2MASS catalog in the following text.

Since the 2MASS catalog is incomplete at higher redshifts, we further consider its completeness. The completeness fraction, P com ( S ), can be expressed as

$$P _ { \text{com} } ( \mathcal { S } ) \equiv \frac { N _ { \text{cat} } ( \mathcal { S } ) } { \bar { n } _ { \text{gal} } V _ { \text{c} } ( \mathcal { S } ) },$$

where S = S ( z, ˆ ; ∆ Ω z, ∆Ω) represents a volume that is a portion of a cone. This cone is defined by the central redshift z with a range from z -∆ z to z +∆ , and the central angle z ˆ Ω with a solid angle ∆Ω. ¯ n gal represents the density of galaxies per comoving volume. V c is the comoving volume of S , and N cat is the number of galaxies in the catalog within S . We assume that the completeness of the 2MASS catalog is isotropic and count the number of galaxies in each redshift bin. The redshift range from 0 to 0.1 is equally divided into 30 bins. Additionally, we assume that galaxies are uniformly distributed within the comoving volume and consider the 2MASS catalog complete at z &lt; 0 04. . Therefore, we can derive the completeness of the 2MASS catalog as a function of redshift using Eq. (2.1), which is shown in figure 1.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00609v2-figure-3f5f77fe-a459-48d7-bcea-9c12c59f4786.jpg)

z

Figure 1 . The red solid lines represent the subcatalog completeness defined in Eq. (2.1) as a function of redshift. The histogram shows the number of galaxies within each redshift bin in the subcatalog. The blue dashed line represents the number of galaxies expected within each redshift bin assuming a uniform distribution of galaxies in comoving volume.

## 2.2 Population models of SMBHBs

The differential comoving number density of SMBHBs can be described as d Φ 3 BHB / d M BH 1 , d d , z q where M BH 1 , is the mass of the larger SMBHB, and q is the mass ratio defined as q = M BHB /M BH 1 , -1 with M BHB being the total binary mass. q is assumed to follow a log-normal distribution, p q ( ), centered at log 10 q = 0 with a standard deviation of σ log 10 q = 0 5 dex, . within the range [0.25, 1.0], similar to the 'hiq' model described in Sesana et al. [156]. While the galactic major merger models can effectively estimate the distribution of SMBHBs, these models are based on the assumptions about the evolution processes of the galaxies and black holes, the galaxy merger rates, and the relationship between the galaxy mass and the black hole mass, etc. However, these assumptions vary greatly among different major merger models, leading to significantly different SMBHB number densities in simulations. It should be noted that the SGWB characteristic amplitude predicted by the M17 model was lower than that of NANOGrav 12.5-yr data. This issue must be addressed if the SGWB

signal is considered to be contributed by SMBHBs. We note that the C21 model [151] can self-consistently incorporate the SGWB characteristic amplitude as a parameter to infer the distribution of SMBHBs, and the number of nearby SMBHBs inferred from the NANOGrav 12.5-yr data is more than five times the number estimated by the M17 major merger model. Therefore, we choose the C21 model to estimate the distribution of SMBHBs and compare it with the M17 major merger model [148, 149].

## 2.2.1 Major merger model

We convert M K to the stellar masses for early-type galaxies using the empirical equation, log 10 M ⋆ = 10 58 . -0 44( . M K + 23) [157]. The masses of SMBHBs in these galaxies are determined via the M BHB -M bulge relationship from McConnell et al. [158], where the bulge mass, M bulge , is approximated to the stellar mass M ⋆ . The mass ratio q of the binary is drawn from a log-normal distribution, p q ( ), centered at log 10 q = 0 and with the standard deviation of σ log 10 q = 0.5 dex, within the range [0.25, 1.0]. Then, the probability of the j -th galaxy hosting an SMBHB in the PTA band, p j , is given by

$$p _ { j } ( M _ { * }, z ) = \frac { t _ { c, j } } { T _ { \text{life} } } \int _ { \mu _ { * } \geq 0. 2 5 } \, \mathrm d \mu _ { * } \frac { \mathrm d N } { \mathrm d t } \left ( M _ { * }, \mu _ { * }, z ^ { \prime } \right ) T _ { \text{life} },$$

where t c ,j = (5 256)( / πf low ) -8 3 / ( GM /c c 3 ) -5 3 / is the time to coalescence of the SMBHB in the j -th galaxy, with f low = 1 nHz being the lower limit of the PTA band. T life = T df + T sh is the lifetime of the SMBHB, where T df and T sh are the dynamical friction timescale [159] and the stellar hardening timescale [160], respectively. d N/ t M ,µ ,z d ( ⋆ ⋆ ′ ) represents the galaxy merger rate derived from the Illustris cosmological simulation project [161, 162], where M ⋆ denotes the stellar mass of the galaxies, µ ⋆ denotes the progenitor stellar mass ratio, and z ′ denotes the redshift at which the galaxies merge, calculated at T z with Planck cosmological parameters [163].

The number of SMBHBs in the galaxy catalog is N SMBHB = ∑ j p j for one simulation. Weperform numerous simulations to obtain a distribution of N SMBHB . From this distribution, we select the simulation that yields the N SMBHB value closest to the maximum likelihood. We find that there are 213 sources emitting GWs in the PTA band in the 2MASS catalog. Then, we randomly select 213 galaxies from the 2MASS catalog as the host galaxies of SMBHBs according to the probability p j .

## 2.2.2 Quasar-based model

The characteristic amplitude h c of the GW signal generated by the population of merging SMBHBs can be expressed as [138, 164, 165]

$$h _ { \text{c} } ^ { 2 } ( f ) = \frac { 4 } { \pi f ^ { 2 } } \iiint \text{d} M _ { \text{BH}, 1 } \text{d} z \text{d} q \frac { \text{d} ^ { 3 } \Phi _ { \text{BHB} } } { \text{d} M _ { \text{BH}, 1 } \text{d} z \text{d} q } \frac { 1 } { 1 + z } \frac { \text{d} E _ { \text{gw} } } { \text{d} \ln f _ { \text{r} } },$$

where d Φ 3 BHB / d M BH 1 , d z d q represents the differential comoving number density of SMBHBs, d E gw / dln f r is the GW energy emitted per unit logarithmic frequency change, and f r = f (1 + z ) is the GW frequency in the rest frame of the source. For a population of SMBHBs in circular orbits, the GW energy per unit logarithmic frequency change is [164, 166]

$$\frac { \mathrm d E _ { \text{gw} } } { \mathrm d \ln f _ { r } } = \frac { ( \pi f _ { r } ) ^ { 2 / 3 } } { 3 } M _ { c } ^ { \ 5 / 3 },$$

where the SMBHB chirp mass is given by M c 5 3 / = q (1+ ) q 2 M 5 3 / BHB . Therefore, the differential comoving number density of SMBHBs is the main parameter influencing the value of h c .

Since some quasars could be driven by the mergers of black holes, the populations of quasars and black holes might be correlated. Therefore, the C21 model makes a simplified assumption that the ratio of the differential comoving number density of quasars to that of black holes remains constant:

$$\frac { d ^ { 3 } \Phi _ { B H B } } { d M _ { B H, 1 } d z d q } \in \frac { d ^ { 2 } \Phi _ { Q } } { d M _ { B H } d z } p ( q ).$$

That is to say, the differential comoving number density of SMBHBs is proportional to the differential comoving number density of quasars, d Φ 2 Q / d M BH d , and an assumed mass z ratio distribution p q ( ). The number density of quasars is derived from the quasar luminosity function and Eddington ratio distribution of the lifetime of the quasar from the galaxy merger simulations in Hopkins et al. [167]. To estimate the number density of quasars, we convert the differential quasar triggering rate [168], d 2 Φ Q / d M BH d t r , into a redshift distribution using the ΛCDM model with WMAP9 data,

$$\frac { d ^ { 2 } \Phi _ { Q } } { d M _ { \text{BH} } d z } = \frac { d ^ { 2 } \Phi _ { Q } } { d M _ { \text{BH} } d t _ { r } } \frac { d t _ { r } } { d z }.$$

According to Eqs. (2.5) and (2.6), we can renormalize the SMBHB population with both the SGWB characteristic amplitude and the quasar number density distribution,

$$\frac { d ^ { 3 } \Phi _ { \text{BHB} } } { d M _ { \text{BH,1} } d z d q } = \frac { \Phi _ { \text{BHB,0} } } { \Phi _ { \text{Q,0} } } \frac { d ^ { 2 } \Phi _ { \text{Q} } } { d M _ { \text{BH} } d t _ { r } } \frac { d t _ { r } } { d z } \frac { p ( q ) } { \int _ { q \geq 0. 2 5 } d q p ( q ) },$$

with

$$\Phi _ { B H B, 0 } = \iint _ { \substack { M _ { B H, 1 } \geq M _ { B H, \det } \\ q \geq 0. 2 5 } } d M _ { B H } d q \frac { d ^ { 3 } \Phi _ { B H B } } { d M _ { B H, 1 } d z d q } \Big | _ { z = 0 },$$

$$\Phi _ { Q, 0 } = \int _ { M _ { B H } \geq M _ { B H, \det } } d M _ { B H } \frac { d ^ { 2 } \Phi _ { Q } } { d M _ { B H } d t _ { r } } \frac { d t _ { r } } { d z } \Big | _ { z = 0 } \,.$$

In this model, the characteristic amplitude of the SGWB, h c , is treated as an input parameter. The Bayesian posterior median of the SGWB characteristic amplitude from the NANOGrav 15-yr dataset is 2 4 . × 10 -15 [130], assuming a spectral index of 13 / 3, as expected for SMBHB inspirals. We substitute h c = 2 4 . × 10 -15 into Eq. (2.3) and use Eqs. (2.3)-(2.9) to calculate the number of SMBHBs. We estimate that there are approximately 688 SMBHBs in the redshift range (0, 0.1). It should be added that using different values of h c results in different numbers of SMBHBs but the C21 model is applicable to the observational signals from any PTA project.

The possible origins of the current signal observed by PTAs can be broadly classified into astrophysical and cosmological sources. In this study, we adopt an astrophysical interpretation, using h c as an input to the C21 model under the assumption that the NANOGrav 15-yr signal is attributed to an ensemble of inspiraling SMBHBs. However, the

true origin of the PTA signals remains uncertain. While an incoherent superposition of GWs from numerous SMBHBs provides a natural astrophysical explanation, several alternative cosmological scenarios have been proposed. These include primordial GWs generated by cosmic inflation [169], scalar-induced GWs [170, 171], GWs produced during cosmological first-order phase transitions [172, 173], GW radiation from cosmic strings [174, 175] and domain walls [176], among others. Each of these mechanisms could contribute to the GW signal in the nHz frequency range. The current observational data are insufficient to distinguish these potential sources, making future high-precision observations crucial for uncovering the true origin of this signal.

Considering the completeness of the 2MASS catalog, we use Eq. (2.1) to calculate the total completeness fraction of the 688 SMBHBs. This completeness fraction actually refers to the proportion of the host galaxies included in the 2MASS catalog relative to the total host galaxies. We find that it is approximately equal to 490. Therefore, we randomly select 490 galaxies based on their P com values from the 2MASS catalog as SMBHBs' host galaxies. We assume that the host galaxies of the other 198 SMBHBs are not included in the 2MASS catalog and ignore their contributions to the subsequent calculations for dark sirens.

Figure 2 . Distributions of h 0 , f 0 , M BHB , and z for simulated SMBHBs. The blue triangles and the orange circles correspond to the C21 and M17 SMBHB population models, respectively.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00609v2-figure-2169b4b7-020e-4111-a63c-8411aed139c3.jpg)

## 2.3 GW detection and Fisher analysis

The timing residuals induced by a single GW source located in the direction represented by the unit vector Ω and measured at time ˆ t on Earth can be expressed as

$$s ( t, \hat { \Omega } ) = F ^ { + } ( \hat { \Omega } ) \Delta A _ { + } ( t ) + F ^ { \times } ( \hat { \Omega } ) \Delta A _ { \times } ( t ).$$

Here, F + ( ˆ ) Ω and F × ( ˆ ) Ω are the antenna pattern functions [148], which can be expressed as

$$F ^ { + } ( \hat { \Omega } ) = & \frac { 1 } { 4 ( 1 - \cos \theta ) } \left \{ \left ( 1 + \sin ^ { 2 } \delta \right ) \cos ^ { 2 } \delta _ { p } \cos \left [ 2 \left ( \alpha - \alpha _ { p } \right ) \right ] - \sin 2 \delta \sin 2 \delta _ { p } \cos \left ( \alpha - \alpha _ { p } \right )$$

$$F ^ { \times } ( \hat { \Omega } ) = & \frac { 1 } { 2 ( 1 - \cos \theta ) } \left \{ \cos \delta \sin 2 \delta _ { p } \sin \left ( \alpha - \alpha _ { p } \right ) - \sin \delta \cos ^ { 2 } \delta _ { p } \sin \left [ 2 \left ( \alpha - \alpha _ { p } \right ) \right ] \right \},$$

$$F ^ { + } ( \hat { \Omega } ) = & \frac { 1 } { 4 ( 1 - \cos \theta ) } \left \{ \left ( 1 + \sin ^ { 2 } \delta \right ) \cos ^ { 2 } \delta _ { \mathrm p } \cos \left [ 2 \left ( \alpha - \alpha _ { \mathrm p } \right ) \right ] - \sin 2 \delta \sin 2 \delta _ { \mathrm p } \cos \left ( \alpha - \alpha _ { \mathrm p } \right ) \\ & + \cos ^ { 2 } \delta \left ( 2 - 3 \cos ^ { 2 } \delta _ { \mathrm p } \right ) \right \}, \\ F ^ { \times } ( \hat { \Omega } ) = & \frac { 1 } { \varrho / 1 \quad \infty \alpha \lambda } \left \{ \cos \delta \sin 2 \delta _ { \mathrm p } \sin \left ( \alpha - \alpha _ { \mathrm p } \right ) - \sin \delta \cos ^ { 2 } \delta _ { \mathrm p } \sin \left [ 2 \left ( \alpha - \alpha _ { \mathrm p } \right ) \right ] \right \},$$

where ( α δ , ) and ( α p , δ p ) represent the right ascension and declination of the GW source and MSP, respectively. θ represents the opening angle between the GW source and MSP,

$$\cos \theta = \cos \delta \cos \delta _ { p } \cos \left ( \alpha - \alpha _ { p } \right ) + \sin \delta \sin \delta _ { p }.$$

In Eq. (2.10), ∆ A { + , ×} ( ) = t A { + , ×} ( ) t -A { + , ×} ( t p ) [143] represents the difference between the Earth term A { + , ×} ( ) t and the pulsar term A { + , ×} ( t p ). t p is the time when the GW passes MSP and can be expressed as t p = t -d p (1 -cos θ /c ) , with d p representing the pulsar distance. For evolving SMBHBs in circular orbits, we have

$$A _ { + } ( t ) = \frac { h ( t ) } { 2 \pi f ( t ) } \{ ( 1 + \cos ^ { 2 } \iota ) \cos 2 \psi \sin [ \phi ( t ) ] + 2 \cos \iota \sin 2 \psi \cos [ \phi ( t ) ] \},$$

$$A _ { \times } ( t ) = \frac { h ( t ) } { 2 \pi f ( t ) } \{ ( 1 + \cos ^ { 2 } \iota ) \sin 2 \psi \sin [ \phi ( t ) ] - 2 \cos \iota \cos 2 \psi \cos [ \phi ( t ) ] \}.$$

Here, ψ is the GW polarization angle, ϕ is the phase of the GW, ι is the inclination angle of the GW, and the GW strain amplitude h t ( ) can be formulated as

$$h ( t ) = 2 \frac { ( G M _ { c } ) ^ { 5 / 3 } } { c ^ { 4 } } \frac { ( \pi f ( t ) ) ^ { 2 / 3 } } { d _ { \text{L} } },$$

where the frequency of GW signals is calculated by the following formula [177],

$$f ( t ) = \pi ^ { - 1 } \left ( \frac { G M _ { c } } { c ^ { 3 } } \right ) ^ { - 5 / 8 } \left [ \frac { 2 5 6 } { 5 } ( t _ { c } - t ) \right ] ^ { - 3 / 8 }.$$

Here, t c is uniformly sampled from [100 yr, 26 Myr]. This range represents the time to coalescence of a SMBHB with M BH 1 , = 10 9 M ⊙ and q = 1 from 1 nHz to 100 nHz. h 0 and f 0 represent the values of h t ( ) and f ( ) t when t = 0, respectively. Figure 2 shows the distribution of mock SMBHBs on the f 0 -h 0 , f 0 -M BHB , z M -BHB planes. The inclination angle, the polarization angle, and the initial phase are randomly sampled in the ranges cos ι ∈ [ -1 1], , ψ ∈ [0 , 2 π ], and ϕ 0 ∈ [0 , 2 π ], respectively. The GW phase, ϕ t ( ), can be expressed as

$$\phi ( t ) = \phi _ { 0 } + \frac { 1 } { 1 6 } \left ( \frac { G M _ { c } } { c ^ { 3 } } \right ) ^ { - 5 / 3 } \left \{ ( \pi f _ { 0 } ) ^ { - 5 / 3 } - [ \pi f ( t ) ] ^ { - 5 / 3 } \right \}.$$

The square of SNR of a GW signal detected by a PTA is defined as

$$\rho ^ { 2 } = \sum _ { j = 1 } ^ { N _ { p } } \sum _ { i = 1 } ^ { N } \left [ \frac { s _ { j } \left ( t _ { i } \right ) } { \sigma _ { t, j } } \right ] ^ { 2 },$$

where N is the total number of the data points of timing residuals for each MSP, N p is the number of MSPs, s j ( t i ) is the timing residual induced by the GW signal in the j -th MSP at time t i , and σ t,j is the RMS of the white noise of the j -th MSP. The SNR threshold is set to 10 in the simulation.

As this paper focuses on the impact of SMBHB population models on dark sirens, we consider only white noise as the noise source, without taking into account red noise or SGWBs. According to the analyses in other papers, taking into account red noise and SGWB could reduce SNRs of GW signals [178], and red noise could worsen the measurement precision of cosmological parameters by 3-5 times [149]. Therefore, we can roughly estimate that, red noise and SGWB would also worsen the constraint precision of H 0 presented in this

paper by a few times. A more detailed analysis of red noise and SGWB will be addressed in our future works.

We utilize the Fisher information matrix for parameter estimations. For a PTA comprising N p independent MSPs, the Fisher information matrix can be expressed as

$$\Gamma _ { a b } = \sum _ { j = 1 } ^ { N _ { p } } \sum _ { i = 1 } ^ { N } \frac { \partial s \left ( t _ { i } \right ) } { \sigma _ { t, j } \partial \theta _ { a } } \frac { \partial s \left ( t _ { i } \right ) } { \sigma _ { t, j } \partial \theta _ { b } },$$

where θ denotes the free parameters to be estimated.

For each GW source, the timing residuals depend on N p +8 system parameters, including eight parameters of the GW source (i.e., M ,α,δ,ι, ψ, ϕ , f c 0 0 , d L ) and the distances of MSPs d p , j ( j = 1 2 , , · · · , N p ). With additional EM information priors, the estimation precision of the parameter can be improved by adding to the appropriate diagonal element of the Fisher information matrix, i.e.,

$$\Gamma _ { a b } \rightarrow \Gamma _ { a b } - \left \langle \frac { \partial ^ { 2 } \ln P \left ( \theta _ { i } \right ) } { \partial \theta _ { a } \partial \theta _ { b } } \right \rangle,$$

where P θ ( i ) is the prior distribution of the parameter θ i [179]. Since the disk direction is randomly distributed in a 4 π solid angle, the prior distribution of the inclination angle ι is given by P ι ( ) ∝ sin( ). ι Consequently, we can modify the Fisher information matrix element Γ ii to Γ ii → Γ ii + 1 sin 2 ι with θ i = . ι

The 1 σ absolute error of the parameter θ a can be estimated by

$$\Delta \theta _ { a } = \left ( \Gamma ^ { - 1 } \right ) _ { a a } ^ { 1 / 2 }.$$

The angular resolution ∆Ω can be calculated according to the uncertainties of the parameters,

$$\Delta \Omega = 2 \pi \sqrt { ( \cos \delta \Delta \delta \Delta \alpha ) ^ { 2 } - \left ( \cos \delta \frac { ( \Gamma ^ { - 1 } ) _ { \delta \alpha } } { \Delta \delta \Delta \alpha } \right ) ^ { 2 } }.$$

For the parameter d L , we first estimate its instrumental error by Eq. (2.21). Besides, we also consider the error arising from weak lensing, which is expressed by the fitting formula [180, 181],

$$\Delta d _ { \text{L} } ^ { \text{lens} } ( z ) = d _ { \text{L} } ( z ) \times 0. 0 6 6 \left [ \frac { 1 - ( 1 + z ) ^ { - 0. 2 5 } } { 0. 2 5 } \right ] ^ { 1. 8 }.$$

The total error of d L is expressed as

$$\Delta d _ { \mathrm L } = \sqrt { \left ( \Delta d _ { \mathrm L } ^ { \text{inst} } \right ) ^ { 2 } + \left ( \Delta d _ { \mathrm L } ^ { \text{lens} } \right ) ^ { 2 } }.$$

The distances of MSPs are expected to be measured via the timing parallax method in the SKA era. We estimate the uncertainties of MSPs' distances by using Eq. (18) in Lee et al. [142],

$$\sigma _ { d _ { \mathrm p } } \simeq \frac { 2. 3 4 } { \cos ^ { 2 } \beta _ { \mathrm p } } \left ( \frac { N } { 1 0 0 } \right ) ^ { - 1 / 2 } \left ( \frac { d _ { \mathrm p } } { 1 \ k p c } \right ) ^ { 2 } \left ( \frac { \sigma _ { t } } { 1 0 \ n s } \right ) \, \mathrm p c,$$

where β p represents the ecliptic latitudes of MSPs. In the Fisher matrix, we take into account the priors of MSPs' distances using Eq. (2.20). In this context, P d ( p ) is a Gaussian distribution with mean d p and the standard deviation σ d p . Under this assumption, the derivatives of GW frequencies can be obtained, and the chirp masses of the SMBHBs can be further measured. This enables us to break the degeneracy between M c and d L , and therefore precisely infer d L from GW amplitudes. The distributions of the d L relative errors and the angular resolutions of simulated SMBHBs versus SNR based on 20, 100, and 200 MSPs with σ t = 20 ns and 10-year observation are displayed in figure 3. The sky locations of MSPs are randomly selected from the Australia Telescope National Facility pulsar catalog [182].

## 2.4 Inference of the Hubble constant

In this work, we adopt the flat ΛCDM model to make the cosmological analysis. The Hubble parameter, which quantifies the rate of expansion of the universe at redshift z , can be expressed as

$$H ( z ) = H _ { 0 } \sqrt { \Omega _ { m } ( 1 + z ) ^ { 3 } + 1 - \Omega _ { m } }, \quad \ \ ( 2. 2 6 ) \\ \pi _ { \mu } \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \end{matrix}$$

where H 0 is the Hubble constant and Ω m is the current matter density parameter. The luminosity distance at redshift z is given by

$$d _ { \text{L} } ( z ) = c ( 1 + z ) \int _ { 0 } ^ { z } \frac { 1 } { H ( z ^ { \prime } ) } \text{d} z ^ { \prime },$$

where c is the speed of light in vacuum. Throughout this work, we adopt the flat ΛCDM model to generate the dark siren data with the best-fit values of Ω m = 0 3111 and . H 0 = 67 66 km s . -1 Mpc -1 given by the Planck 2018 TT, TE, EE+lowE+lensing+BAO data [3]. Based on the flat ΛCDM model, we utilize the GW signals from SMBHBs as dark sirens to infer H 0 , fixing the other cosmological parameters to the Planck 2018 results.

Using the Bayesian method to infer the cosmological parameters, the posterior probability distribution of H 0 can be expressed as

$$p \left ( H _ { 0 } \, | \, \mathcal { D } _ { \text{GW} }, \mathcal { D } _ { \text{EM} } \right ) = \frac { \prod _ { i = 1 } ^ { N _ { \text{GW} } } p \left ( \mathcal { D } _ { \text{GW}, i }, \mathcal { D } _ { \text{EM}, i } \, | \, H _ { 0 } \right ) p ( H _ { 0 } ) } { \gamma ( H _ { 0 } ) },$$

where D GW and D EM represent the detected GW data and the EM data, respectively. p H ( 0 ) is the prior distribution of H 0 , and N GW is the number of GW events. p ( D GW ,i , D EM ,i | H 0 ) is the likelihood function of a single GW event, which can be written as

p

(

GW

,i

,

EM

,i

H

0

)

$$= \int _ { \Phi } p \left ( \mathcal { D } _ { \text{GW}, i }, \mathcal { D } _ { \text{EM}, i }, d _ { \text{L} }, \alpha, \delta, z \ | \ H _ { 0 } \right ) \text{d} d _ { \text{L} } \text{d} \delta \text{d} z$$

D

D

|

$$= \overset { \mu } { \int } p \left ( \mathcal { D } _ { \text{GW}, i } \, | \, d _ { L }, \alpha, \delta \right ) p \left ( \mathcal { D } _ { \text{EM}, i } \, | \, z, \alpha, \delta \right ) p \left ( d _ { \text{L} } \, | \, z, H _ { 0 } \right ) p _ { 0 } \left ( z, \alpha, \delta \, | \, H _ { 0 } \right ) d d _ { L } d \alpha d \delta d z$$

$$= \overset { \kappa } { \int } _ { \kappa } p \left ( \mathcal { D } _ { \text{GW}, i } \left | \, d _ { \text{L} }, \alpha, \delta \right ) p \left ( \mathcal { D } _ { \text{EM}, i } \left | \, z, \alpha, \delta \right ) \delta \left ( d _ { \text{L} } - \hat { d } _ { \text{L} } \left ( z, H \right ) \right ) p _ { 0 } \left ( z, \alpha, \delta \right | H _ { 0 } \right ) \text{d} d _ { \text{L} } \text{d} \delta \text{d} z$$

$$& = \overset { \nu } { \int } p \left ( \mathcal { D } _ { \text{GW}, i } \left | \, \hat { d } _ { \text{L} } \left ( z, H _ { 0 } \right ), \alpha, \delta \right ) p \left ( \mathcal { D } _ { \text{EM}, i } \left | \, z, \alpha, \delta \right ) p _ { 0 } \left ( z, \alpha, \delta \left | \, H _ { 0 } \right ) \text{dad} \delta d z. \\ & \quad. \int \infty \quad \ \ + \hat { d } \ \ + \ \pi \ \times \ \ - \ \times \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ - \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \ \ + \$$

p ( D GW ,i | ˆ d L ( z, H 0 ) , α, δ ) is the likelihood of the GW data, which can be expressed as

$$p \left ( \mathcal { D } _ { \text{GW}, i } \, | \, \hat { d } _ { \text{L} } \left ( z, H _ { 0 } \right ), \alpha, \delta \right ) \in e ^ { - x ^ { 2 } / 2 },$$

Figure 3 . Distributions of the d L relative errors and the angular resolutions of simulated SMBHBs versus SNR for different PTA detection cases. The SMBHBs are simulated based on the M17 model (upper panel) and the C21 model (lower panel), respectively. From left to right, the three columns correspond to different PTA detection cases with N p set to 20, 100, and 200, respectively. Observations are conducted over ten years with σ t = 20 ns. The red dotted vertical and horizontal lines denote the thresholds of ρ = 10 and ∆ d L /d L = 1, respectively. The purple dots represent the GW events with SNR &gt; 10 and ∆ d L /d L &lt; 1, which are detectable and suitable as dark sirens.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00609v2-figure-730d0f3d-36e3-4d33-90e6-3698f42a868a.jpg)

with χ 2 = ( x -x gw ) T C -1 ( x -x gw ), where C is the 3 × 3 covariance matrix only relevant to ( d L , α, δ ), obtained from the Fisher information matrix. x = ( ˆ d L ( z, H 0 ) , α, δ ) represents the

three-dimensional (3D) position in the sky. x gw = ( d L ,i , α , δ i i ) represents the 3D position of the GW source. Considering the bias between the true value and the posterior median from the actual observation, we use the Fisher matrix to determine the Gaussian distribution of x gw and randomly sample x gw from this distribution. We establish the boundary of the GW source's localization volume based on the 3D measurement errors of x gw . The galaxies whose positions satisfying χ 2 ≤ 11 34 (corresponding to 99% confidence) can be considered as the . potential host galaxies of the GW source. The number of potential host galaxies for a GW source is defined as N in , which is dependent on the prior values of H 0 . We display N in as a function of SNR in figure 4, fixing H 0 = 67 66 km s . -1 Mpc -1 . Due to the lack of EM signals in dark sirens, p ( D EM ,i | z, α, δ ) can be set to a constant value [89].

For the potential candidate host galaxies, p 0 ( z, α, δ | H 0 ) in Eq. (2.29) can be expressed as

$$\begin{array} {$$

where f com is the completeness fraction of the catalog, which can be expressed as

$$f _ { \text{com} } = \frac { N _ { \text{cat} }$$

where N cat is the number of galaxies in the galaxy catalog, and V max is the entire comoving volume bounded by the upper limit of redshift, i.e., z max = 0 1. . p cat ( z, α, δ ) is the observed galaxy distribution function provided by the galaxy catalog, which is expressed as

$$p _ { \text{cat} } ( z, \alpha, \delta ) = \frac$$

where N represents the Gaussian distribution with mean z j and the standard deviation σ z,j = (1 + z ) √ ⟨ v p 2 ⟩ /c [103, 183]. √ ⟨ v p 2 ⟩ is the standard deviation of the radial peculiar velocity of a galaxy, set to 500 km s -1 [184]. The expression for p miss is given by

$$p _ { \text{miss} } ( z, \alpha, \delta \, | \$$

The normalization term γ H ( 0 ) in Eq. (2.28) can be written as

$$\gamma \left ( H _ { 0 } \right ) = \int p _ { \text$$

where p GW det and p EM det represent the GW detection probability and the EM detection probability, respectively. The GW detection probability can be expressed as

$$p _ { \text{det} } ^ { \text{GW} } \left ( d _ { \text{L} } \left ( z, H _ { 0 } \right ), \alpha, \delta \right ) = \int _ { \mathcal { D } _ { \text{GW} } > \mathcal { D } _ { \text{GW} } ^ { \text{th} } } p \left ( \mathcal { D } _ { \text{GW} } \, \left | \, d _ { \text{L} } \left ( z, H _ { 0 } \right ), \alpha, \delta \right ) \mathrm d \mathcal { D } _ { \text{GW} },$$

where D th GW denotes the SNR threshold for GW detections and is set to 10. The EM detection probability can be expressed as p EM det ( z, α, δ ) ∝ H ( z max -z ), where H is the Heaviside step function and z max = 0 1. .

M17 model

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00609v2-figure-0dc7112d-652d-40b7-b94b-bbae503ee704.jpg)

C21 model

Figure 4 . Correlations between SNR and the number of potential host galaxies N in for the mock GW events based on the M17 model (upper panel) and the C21 model (lower panel) under different PTA detection cases. The dots with different colors and shapes represent different PTA detection cases with N p = 20 100 and 200, respectively. The observation period is 10 years and , , σ t is set to 20 ns.

| Detection case   | T span [yr]   | M17 model   | M17 model   | C21 model   | C21 model   |
|------------------|---------------|-------------|-------------|-------------|-------------|
| Detection case   | T span [yr]   | N s         | ε ( H 0 )   | N s         | ε ( H 0 )   |
|                  | 5.0           | 49          | 7.75%       | 174         | 3.56%       |
|                  | 10.0          | 63          | 4.49%       | 237         | 2.91%       |
|                  | 15.0          | 67          | 3.10%       | 274         | 2.27%       |
|                  | 20.0          | 80          | 2.58%       | 295         | 1.82%       |
| σ                | 5.0           | 105         | 2.21%       | 377         | 1.48%       |
| σ                | 10.0          | 130         | 1.78%       | 421         | 1.18%       |
| σ                | 15.0          | 138         | 1.54%       | 451         | 0.97%       |
| σ                | 20.0          | 150         | 1.44%       | 461         | 0.91%       |
| σ t              | 5.0           | 133         | 1.72%       | 439         | 1.11%       |
| σ t              | 10.0          | 152         | 1.36%       | 457         | 0.88%       |
| σ t              | 15.0          | 164         | 1.18%       | 481         | 0.77%       |
| σ t              | 20.0          | 173         | 1.11%       | 487         | 0.68%       |

Table 1 . Relative errors of H 0 using the mock dark sirens from different PTA detection cases based on the M17 model and the C21 model, respectively. Here, N s represents the number of detected SMBHBs with ρ &gt; 10 and ∆ d L /d L &lt; 1 for different T span and N p .

## 3 Results and discussions

In this section, we compare the two SMBHB population models, discuss the GW parameter estimations, and present the constraint results for H 0 from various PTA detection cases. As shown in figure 2, the h 0 , z , M BHB , and f 0 distributions of SMBHBs simulated from the two population models are essentially similar. Nevertheless, the numbers of simulated SMBHBs for the two population models are 490 and 213, respectively, showing that the C21 model has a noticeably larger number of SMBHBs than the M17 model. We consider the 10-year observations based on the C21 model as the representative example of PTA observations for the following discussions.

In figure 5, we present the distributions of ρ , ∆Ω, and ∆ d L /d L based on the 20, 100, and 200 MSPs with σ t = 20 ns. Note that the calculations are based on 213 (M17 model) and

M17 model

Figure 5 . Distributions of ρ , ∆Ω, and ∆ d L /d L for simulated SMBHBs in a 10-year observation with σ t =20 ns. The upper and lower panels correspond to the M17 and C21 models, respectively. The lines in different colors and line styles represent different PTA detection cases with N p = 20 100 and 200 , , respectively.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00609v2-figure-a55396bc-ee92-4bfd-80ef-cd938a2a7fda.jpg)

490 (C21 model) simulated SMBHBs at z &lt; 0 1. . As shown in the left panel of figure 5, fewer available MSPs result in lower SNRs, thus showing worse localization abilities. We note that the localization capability with 100 MSPs is much better than that of 20 MSPs, but 200 MSPs are slightly better than 100 MSPs. This indicates that 100 MSPs seem to be sufficient and more economical, as further increasing the number of MSPs does not significantly enhance the localization ability.

For the mock SMBHBs from the two population models, we consider PTA cases of monitoring 20, 100, and 200 MSPs with σ t = 20 ns to detect them. The observation time spans are T span = 5 10 15 20 years, respectively, with a cadence of 2 weeks that is the typical , , , cadence used in current PTAs. The detected SMBHBs are used as dark sirens to constrain H 0 .

Figure 6 illustrates the relative errors for H 0 as a function of T span . Table 1 shows the relative error of H 0 , defined as ε H ( 0 ) = σ H ( 0 ) /H 0 . We observe that the PTAs with

Figure 6 . The constraint precision of H 0 from the mock dark sirens as a function of the PTA observation time span. The dotted and solid lines represent the constraint results from the M17 and C21 models, respectively. Different colors represent N p = 20 100 and 200 respectively. , , Here, σ t is set to 20 ns.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00609v2-figure-50398b7f-da7d-4aa0-9aa5-2ab1e986a248.jpg)

N p = 20 give the poorest measurement precision for H 0 . Based on the C21 model, the constraint precision of H 0 is 3.56% for a 5-year observation, improving to 1.82% over a 20year observation. These values are worse than the constraint precision obtained from current distance ladder measurements. Increasing the number of MSPs to 100 significantly enhances the precision of H 0 compared to using 20 MSPs. With a 5-year observation, the constraint precision of H 0 is 1.48% and can be improved to 0.91% with a 20-year observation; both values are close to or better than 1%. Furthermore, monitoring 200 MSPs provides even better constraints on H 0 , achieving precisions ranging from 1.11% to 0.68% over 5 to 20 years of observations. Considering that SKA is expected to begin operations around 2030, our optimistic estimate ( σ t = 20 ns, N p = 200) suggests that PTA dark sirens can achieve a measurement precision of approximately 1% for H 0 over a 10-year observation, which would be around 2040. Nevertheless, under a more conservative noise assumption and considering the impact of the SGWB, the first results are more likely to be obtained after 2050. In figure 6, we plot the current H 0 values obtained from the distance ladder measurements and the Planck CMB data. It should be noted that when the precision of the PTA dark siren method in measuring H 0 reaches its expected level in the future, the measurement accuracy of the distance ladder method will be also significantly enhanced. Although it is difficult for the PTA dark siren method to surpass the distance ladder method in terms of precisely measuring H 0 , it can still provide assistance in resolving the Hubble tension as an independent late-universe probe.

It is worth noting that our results are influenced by the value of h c . Since h 2 c ∝ Φ BHB 0 , , a larger h c corresponds to a greater inferred number of SMBHBs. This increases the number of detectable SMBHBs serving as dark sirens, and further enhances the measurement precision of H 0 . Compared with the NANOGrav 12.5-yr data, the h c value inferred from the NANOGrav 15-yr dataset suggests a larger SMBHB population. However, the h c values reported by NANOGrav [130], EPTA+InPTA [132], PPTA [134] and CPTA [136] in 2023 are generally consistent. Therefore, using different h c values provided by various PTA projects

will not alter our main conclusions.

Figure 7 presents the constraints on H 0 obtained from individual GW events as well as the combined results. The colored lines represent the constraints from individual events and the black line represents the combined result of all events. There is a notable discrepancy between the two SMBHB population models in terms of the H 0 posterior distributions given by the individual events. In the M17 model, there are several events whose SNRs far exceed the others, but in the C21 model, there is no significant difference in the SNRs given by all the events. This discrepancy can also be seen in figure 2: Several SMBHBs (represented by the orange circles) in the M17 model have significantly higher masses and GW strains than the others. This is due to the fact that the mechanisms of the two models in predicting the SMBHB populations are substantially different. The M17 model predicts the masses of SMBHBs based on the apparent magnitude of the galaxies in the galaxy catalog. A few of these galaxies have significantly larger masses, and the SMBHBs inhabiting them could produce the GW signals with notably higher SNRs than the other SMBHBs. However, for the C21 model, its SMBHB distribution is derived from the SGWB signal, which is the overall superposition of many GW events and does not reflect the information of individual sources. As a result, the SMBHB distribution predicted by the C21 model is statistical rather than individualistic. Therefore, in the C21 model, no event exhibits a SNR that is significantly higher than other events.

PTA dark sirens primarily rely on low-redshift ( z &lt; 0 1) GW sources and galaxy cata-. logs. Within this range, the distance-redshift relation can be approximated by Hubble's law and has limited dependence on late-time cosmological extended models with new physics, such as time-varying dark energy [185, 186], interacting dark energy [26, 187], phantom dark energy [188, 189], decaying dark matter [190, 191], and spatial curvature [192, 193]. Consequently, these additional parameters cannot be effectively constrained by the low-redshift PTA dark sirens. If these extended models are considered, the additional introduced parameters would exhibit degeneracy with H 0 , thereby slightly diminishing the measurement precision of H 0 .

Compared with the dark siren analysis in Wang et al. [149], under the same detection conditions and population models, the constraint precision of H 0 is improved in this paper. This is mainly due to the fact that we simulate SMBHBs over a broader range of redshift (0 &lt; z &lt; 0 1), . resulting in a larger number of GW sources available as dark sirens. We further consider the effect of incompleteness for galaxies at z &gt; 0 04. . Wang et al. [149] also indicated some advantages of nHz GW standard sirens over the standard sirens in other frequency bands. Here, we want to emphasize another advantage: Since SMBHBs are all hosted in supermassive galaxies, we can directly exclude those galaxies with low masses when using the dark siren method, and only need to consider the supermassive galaxies within the localization volume of the GW source. Compared with the stellar-mass dark sirens, it can greatly reduce the computational costs and the uncertainties caused by spurious host galaxies in Bayesian inference.

## 4 Conclusion

Several PTA projects have reported substantial evidence for the existence of a SGWB in the nHz band. In addition to SGWB, the GW signals from individual SMBHBs may also be observed in the future. Such individual GW events can be used as dark sirens to constrain cosmological parameters. In the previous studies on dark sirens, the simulations of SMBHBs

were based on the M17 model, which cannot fit the observed SGWB signal well. In this paper, we use the C21 model, which can self-consistently predict the local number density of SMBHBs given the amplitude of the SGWB, to investigate the potential of SMBHB dark sirens in measuring the Hubble constant. According to our results, the following conclusions can be drawn:

- (i) In the C21 model, a PTA with 100 MSPs over a 10-year observation period could measure H 0 with a precision close to 1%. If the number of MSPs increases from 100 to 200, a 5-year observation could measure H 0 with a precision close to 1%.
- (ii) The C21 model predicts a larger number of SMBHBs than the M17 model. The numbers of SMBHBs derived from the M17 and C21 models are 490 and 213, respectively.
- (iii) The SMBHB distribution in the C21 model is derived from the SGWB signal, which is the overall superposition of the GW signals from individual events. Therefore, the mass distribution of SMBHBs predicted by the C21 model is more statistical than those predicted by the M17 model.

Our conclusions indicate that modeling the SMBHB population has a great impact on the analysis of SMBHB dark sirens. In the future, more accurate measurements of the SGWB signal will not only give us a better understanding of the distribution and evolution history of SMBHBs, but also have important implications for measuring cosmological parameters via the dark siren method.

## Acknowledgments

We thank Shang-Jie Jin for helpful discussions. This work was supported by the National SKA Program of China (Grants Nos. 2022SKA0110200 and 2022SKA0110203), the National Natural Science Foundation of China (Grants Nos. 12473001, 11975072, 11875102, 11835009, 12305069, and 12305058), the National 111 Project (Grant No. B16009), the Program of the Education Department of Liaoning Province (Grant No. JYTMS20231695), and the Natural Science Foundation of Hainan Province (Grant No. 4240N215).

## References

- [1] WMAP collaboration, First year Wilkinson Microwave Anisotropy Probe (WMAP) observations: Determination of cosmological parameters , Astrophys. J. Suppl. 148 (2003) 175 [ astro-ph/0302209 ].
- [2] WMAP collaboration, First year Wilkinson Microwave Anisotropy Probe (WMAP) observations: Preliminary maps and basic results , Astrophys. J. Suppl. 148 (2003) 1 [ astro-ph/0302207 ].
- [3] Planck collaboration, Planck 2018 results. VI. Cosmological parameters , Astron. Astrophys. 641 (2020) A6 [ 1807.06209 ].
- [4] A.G. Riess et al., A Comprehensive Measurement of the Local Value of the Hubble Constant with 1 km s -1 Mpc -1 Uncertainty from the Hubble Space Telescope and the SH0ES Team , Astrophys. J. Lett. 934 (2022) L7 [ 2112.04510 ].
- [5] L. Verde, T. Treu and A.G. Riess, Tensions between the Early and the Late Universe , Nature Astron. 3 (2019) 891 [ 1907.10625 ].
- [6] A.G. Riess, The Expansion of the Universe is Faster than Expected , Nature Rev. Phys. 2 (2019) 10 [ 2001.03624 ].

- [7] M. Li, X.-D. Li, Y.-Z. Ma, X. Zhang and Z. Zhang, Planck Constraints on Holographic Dark Energy , JCAP 09 (2013) 021 [ 1305.5302 ].
- [8] J.-F. Zhang, J.-J. Geng and X. Zhang, Neutrinos and dark energy after Planck and BICEP2: data consistency tests and cosmological parameter constraints , JCAP 10 (2014) 044 [ 1408.0481 ].
- [9] J.-F. Zhang, Y.-H. Li and X. Zhang, Cosmological constraints on neutrinos after BICEP2 , Eur. Phys. J. C 74 (2014) 2954 [ 1404.3598 ].
- [10] J.-F. Zhang, Y.-H. Li and X. Zhang, Sterile neutrinos help reconcile the observational results of primordial gravitational waves from Planck and BICEP2 , Phys. Lett. B 740 (2015) 359 [ 1403.7028 ].
- [11] M.-M. Zhao, D.-Z. He, J.-F. Zhang and X. Zhang, Search for sterile neutrinos in holographic dark energy cosmology: Reconciling Planck observation with the local measurement of the Hubble constant , Phys. Rev. D 96 (2017) 043520 [ 1703.08456 ].
- [12] R.-Y. Guo and X. Zhang, Constraints on inflation revisited: An analysis including the latest local measurement of the Hubble constant , Eur. Phys. J. C 77 (2017) 882 [ 1704.04784 ].
- [13] R.-Y. Guo, J.-F. Zhang and X. Zhang, Can the H 0 tension be resolved in extensions to Λ CDM cosmology? , JCAP 02 (2019) 054 [ 1809.02340 ].
- [14] W. Yang, S. Pan, E. Di Valentino, R.C. Nunes, S. Vagnozzi and D.F. Mota, Tale of stable interacting dark energy, observational signatures, and the H 0 tension , JCAP 09 (2018) 019 [ 1805.08252 ].
- [15] R.-Y. Guo, J.-F. Zhang and X. Zhang, Inflation model selection revisited after a 1.91% measurement of the Hubble constant , Sci. China Phys. Mech. Astron. 63 (2020) 290406 [ 1910.13944 ].
- [16] R.-G. Cai, Editorial , Sci. China Phys. Mech. Astron. 63 (2020) 290401.
- [17] E. Di Valentino, A. Melchiorri, O. Mena and S. Vagnozzi, Nonminimal dark sector physics and cosmological tensions , Phys. Rev. D 101 (2020) 063502 [ 1910.09853 ].
- [18] E. Di Valentino, A. Melchiorri, O. Mena and S. Vagnozzi, Interacting dark energy in the early 2020s: A promising solution to the H 0 and cosmic shear tensions , Phys. Dark Univ. 30 (2020) 100666 [ 1908.04281 ].
- [19] M. Liu, Z. Huang, X. Luo, H. Miao, N.K. Singh and L. Huang, Can Non-standard Recombination Resolve the Hubble Tension? , Sci. China Phys. Mech. Astron. 63 (2020) 290405 [ 1912.00190 ].
- [20] X. Zhang and Q.-G. Huang, Measuring H 0 from low-z datasets , Sci. China Phys. Mech. Astron. 63 (2020) 290402 [ 1911.09439 ].
- [21] Q. Ding, T. Nakama and Y. Wang, A gigaparsec-scale local void and the Hubble tension , Sci. China Phys. Mech. Astron. 63 (2020) 290403 [ 1912.12600 ].
- [22] R.-G. Cai, Z.-K. Guo, L. Li, S.-J. Wang and W.-W. Yu, Chameleon dark energy can resolve the Hubble tension , Phys. Rev. D 103 (2021) 121302 [ 2102.02020 ].
- [23] E. Di Valentino, O. Mena, S. Pan, L. Visinelli, W. Yang, A. Melchiorri et al., In the realm of the Hubble tension-a review of solutions , Class. Quant. Grav. 38 (2021) 153001 [ 2103.01183 ].
- [24] E. Abdalla, G.F. Abellan, A. Aboubrahim et al., Cosmology intertwined: A review of the particle physics, astrophysics, and cosmology associated with the cosmological tensions and anomalies , JHEAp 34 (2022) 49 [ 2203.06142 ].
- [25] H. Li and X. Zhang, A novel method of measuring cosmological distances using broad-line regions of quasars , Sci. Bull. 65 (2020) 1419 [ 2005.10458 ].

- [26] L.-F. Wang, J.-H. Zhang, D.-Z. He, J.-F. Zhang and X. Zhang, Constraints on interacting dark energy models from time-delay cosmography with seven lensed quasars , Mon. Not. Roy. Astron. Soc. 514 (2022) 1433 [ 2102.09331 ].
- [27] S. Vagnozzi, F. Pacucci and A. Loeb, Implications for the Hubble tension from the ages of the oldest astrophysical objects , JHEAp 36 (2022) 27 [ 2105.10421 ].
- [28] S. Vagnozzi, Consistency tests of Λ CDM from the early integrated Sachs-Wolfe effect: Implications for early-time new physics and the Hubble tension , Phys. Rev. D 104 (2021) 063524 [ 2105.10425 ].
- [29] S. Vagnozzi, New physics in light of the H 0 tension: An alternative view , Phys. Rev. D 102 (2020) 023518 [ 1907.07569 ].
- [30] L. Feng, D.-Z. He, H.-L. Li, J.-F. Zhang and X. Zhang, Constraints on active and sterile neutrinos in an interacting dark energy cosmology , Sci. China Phys. Mech. Astron. 63 (2020) 290404 [ 1910.03872 ].
- [31] M.-X. Lin, W. Hu and M. Raveri, Testing H 0 in Acoustic Dark Energy with Planck and ACT Polarization , Phys. Rev. D 102 (2020) 123523 [ 2009.08974 ].
- [32] A. Hryczuk and K. Jodglyph[suppress] lowski, Self-interacting dark matter from late decays and the H 0 tension , Phys. Rev. D 102 (2020) 043024 [ 2006.16139 ].
- [33] L.-Y. Gao, Z.-W. Zhao, S.-S. Xue and X. Zhang, Relieving the H 0 tension with a new interacting dark energy model , JCAP 07 (2021) 005 [ 2101.10714 ].
- [34] Z.-W. Zhao, J.-G. Zhang, Y. Li, J.-M. Zou, J.-F. Zhang and X. Zhang, First statistical measurement of the Hubble constant using unlocalized fast radio bursts , arXiv e-prints (2022) [ 2212.13433 ].
- [35] L.-Y. Gao, S.-S. Xue and X. Zhang, Dark energy and matter interacting scenario to relieve H 0 and S 8 tensions* , Chin. Phys. C 48 (2024) 051001 [ 2212.13146 ].
- [36] B.F. Schutz, Determining the Hubble Constant from Gravitational Wave Observations , Nature 323 (1986) 310.
- [37] D.E. Holz and S.A. Hughes, Using gravitational-wave standard sirens , Astrophys. J. 629 (2005) 15 [ astro-ph/0504616 ].
- [38] N. Dalal, D.E. Holz, S.A. Hughes and B. Jain, Short grb and binary black hole standard sirens as a probe of dark energy , Phys. Rev. D 74 (2006) 063006 [ astro-ph/0601275 ].
- [39] C. Cutler and D.E. Holz, Ultra-high precision cosmology from gravitational waves , Phys. Rev. D 80 (2009) 104009 [ 0906.3752 ].
- [40] W. Zhao, C. Van Den Broeck, D. Baskaran and T.G.F. Li, Determination of Dark Energy by the Einstein Telescope: Comparing with CMB, BAO and SNIa Observations , Phys. Rev. D 83 (2011) 023005 [ 1009.0206 ].
- [41] R.-G. Cai and T. Yang, Estimating cosmological parameters by the simulated data of gravitational waves from the Einstein Telescope , Phys. Rev. D 95 (2017) 044024 [ 1608.08008 ].
- [42] R.-G. Cai and T. Yang, Standard sirens and dark sector with Gaussian process , EPJ Web Conf. 168 (2018) 01008 [ 1709.00837 ].
- [43] R.-G. Cai, T.-B. Liu, X.-W. Liu, S.-J. Wang and T. Yang, Probing cosmic anisotropy with gravitational waves as standard sirens , Phys. Rev. D 97 (2018) 103005 [ 1712.00952 ].
- [44] W. Zhao, B.S. Wright and B. Li, Constraining the time variation of Newton's constant G with gravitational-wave standard sirens and supernovae , JCAP 10 (2018) 052 [ 1804.03066 ].
- [45] L.-F. Wang, X.-N. Zhang, J.-F. Zhang and X. Zhang, Impacts of gravitational-wave standard siren observation of the Einstein Telescope on weighing neutrinos in cosmology , Phys. Lett. B 782 (2018) 87 [ 1802.04720 ].

- [46] X.-N. Zhang, L.-F. Wang, J.-F. Zhang and X. Zhang, Improving cosmological parameter estimation with the future gravitational-wave standard siren observation from the Einstein Telescope , Phys. Rev. D 99 (2019) 063510 [ 1804.08379 ].
- [47] J.-F. Zhang, H.-Y. Dong, J.-Z. Qi and X. Zhang, Prospect for constraining holographic dark energy with gravitational wave standard sirens from the Einstein Telescope , Eur. Phys. J. C 80 (2020) 217 [ 1906.07504 ].
- [48] H.-L. Li, D.-Z. He, J.-F. Zhang and X. Zhang, Quantifying the impacts of future gravitational-wave data on constraining interacting dark energy , JCAP 06 (2020) 038 [ 1908.03098 ].
- [49] Y.-F. Cai, C. Li, E.N. Saridakis and L. Xue, f ( T ) gravity after GW170817 and GRB170817A , Phys. Rev. D 97 (2018) 103513 [ 1801.05827 ].
- [50] Z. Chang, Q.-G. Huang, S. Wang and Z.-C. Zhao, Low-redshift constraints on the Hubble constant from the baryon acoustic oscillation 'standard rulers' and the gravitational wave 'standard sirens' , Eur. Phys. J. C 79 (2019) 177.
- [51] M. Du, W. Yang, L. Xu, S. Pan and D.F. Mota, Future constraints on dynamical dark-energy using gravitational-wave standard sirens , Phys. Rev. D 100 (2019) 043535 [ 1812.01440 ].
- [52] J.-h. He, Accurate method to determine the systematics due to the peculiar velocities of galaxies in measuring the Hubble constant from gravitational-wave standard sirens , Phys. Rev. D 100 (2019) 023527 [ 1903.11254 ].
- [53] W. Yang, S. Vagnozzi, E. Di Valentino, R.C. Nunes, S. Pan and D.F. Mota, Listening to the sound of dark sector interactions with gravitational wave standard sirens , JCAP 07 (2019) 037 [ 1905.08286 ].
- [54] W. Yang, S. Pan, E. Di Valentino, B. Wang and A. Wang, Forecasting interacting vacuum-energy models using gravitational waves , JCAP 05 (2020) 050 [ 1904.11980 ].
- [55] X. Zhang, Gravitational wave standard sirens and cosmological parameter measurement , Sci. China Phys. Mech. Astron. 62 (2019) 110431 [ 1905.11122 ].
- [56] R.R.A. Bachega, A.A. Costa, E. Abdalla and K.S.F. Fornazier, Forecasting the Interaction in Dark Matter-Dark Energy Models with Standard Sirens From the Einstein Telescope , JCAP 05 (2020) 021 [ 1906.08909 ].
- [57] J.-F. Zhang, M. Zhang, S.-J. Jin, J.-Z. Qi and X. Zhang, Cosmological parameter estimation with future gravitational wave standard siren observation from the Einstein Telescope , JCAP 09 (2019) 068 [ 1907.03238 ].
- [58] H.-Y. Chen, Systematic Uncertainty of Standard Sirens from the Viewing Angle of Binary Neutron Star Inspirals , Phys. Rev. Lett. 125 (2020) 201301 [ 2006.02779 ].
- [59] S.-J. Jin, D.-Z. He, Y. Xu, J.-F. Zhang and X. Zhang, Forecast for cosmological parameter estimation with gravitational-wave standard siren observation from the Cosmic Explorer , JCAP 03 (2020) 051 [ 2001.05393 ].
- [60] H.-Y. Chen, P.S. Cowperthwaite, B.D. Metzger and E. Berger, A Program for Multimessenger Standard Siren Cosmology in the Era of LIGO A+, Rubin Observatory, and Beyond , Astrophys. J. Lett. 908 (2021) L4 [ 2011.01211 ].
- [61] A. Mitra, J. Mifsud, D.F. Mota and D. Parkinson, Cosmology with the Einstein Telescope: No Slip Gravity Model and Redshift Specifications , Mon. Not. Roy. Astron. Soc. 502 (2021) 5563 [ 2010.00189 ].
- [62] Z.-W. Zhao, Z.-X. Li, J.-Z. Qi, H. Gao, J.-F. Zhang and X. Zhang, Cosmological parameter estimation for dynamical dark energy models with future fast radio burst observations , Astrophys. J. 903 (2020) 83 [ 2006.01450 ].

- [63] S.-J. Jin, L.-F. Wang, P.-J. Wu, J.-F. Zhang and X. Zhang, How can gravitational-wave standard sirens and 21-cm intensity mapping jointly provide a precise late-universe cosmological probe? , Phys. Rev. D 104 (2021) 103507 [ 2106.01859 ].
- [64] S.-J. Jin, R.-Q. Zhu, L.-F. Wang, H.-L. Li, J.-F. Zhang and X. Zhang, Impacts of gravitational-wave standard siren observations from Einstein Telescope and Cosmic Explorer on weighing neutrinos in interacting dark energy models , Commun. Theor. Phys. 74 (2022) 105404 [ 2204.04689 ].
- [65] S.-J. Jin, T.-N. Li, J.-F. Zhang and X. Zhang, Prospects for measuring the Hubble constant and dark energy using gravitational-wave dark sirens with neutron star tidal deformation , JCAP 08 (2023) 070 [ 2202.11882 ].
- [66] S.-J. Jin, S.-S. Xing, Y. Shao, J.-F. Zhang and X. Zhang, Joint constraints on cosmological parameters using future multi-band gravitational wave standard siren observations* , Chin. Phys. C 47 (2023) 065104 [ 2301.06722 ].
- [67] S.-J. Jin, R.-Q. Zhu, J.-Y. Song, T. Han, J.-F. Zhang and X. Zhang, Standard siren cosmology in the era of the 2.5-generation ground-based gravitational wave detectors: bright and dark sirens of LIGO Voyager and NEMO , arXiv e-prints (2023) [ 2309.11900 ].
- [68] S.-J. Jin, Y.-Z. Zhang, J.-Y. Song, J.-F. Zhang and X. Zhang, Taiji-TianQin-LISA network: Precisely measuring the Hubble constant using both bright and dark sirens , Sci. China Phys. Mech. Astron. 67 (2024) 220412 [ 2305.19714 ].
- [69] J.-Z. Qi, S.-J. Jin, X.-L. Fan, J.-F. Zhang and X. Zhang, Using a multi-messenger and multi-wavelength observational strategy to probe the nature of dark energy through direct measurements of cosmic expansion history , JCAP 12 (2021) 042 [ 2102.01292 ].
- [70] C. Ye and M. Fishbach, Cosmology with standard sirens at cosmic noon , Phys. Rev. D 104 (2021) 043507 [ 2103.14038 ].
- [71] M.-D. Cao, J. Zheng, J.-Z. Qi, X. Zhang and Z.-H. Zhu, A New Way to Explore Cosmological Tensions Using Gravitational Waves and Strong Gravitational Lensing , Astrophys. J. 934 (2022) 108 [ 2112.14564 ].
- [72] J.M.S. de Souza, R. Sturani and J. Alcaniz, Cosmography with standard sirens from the Einstein Telescope , JCAP 03 (2022) 025 [ 2110.13316 ].
- [73] L.-F. Wang, S.-J. Jin, J.-F. Zhang and X. Zhang, Forecast for cosmological parameter estimation with gravitational-wave standard sirens from the LISA-Taiji network , Sci. China Phys. Mech. Astron. 65 (2022) 210411 [ 2101.11882 ].
- [74] E.O. Colg´in, a Probing the Anisotropic Universe with Gravitational Waves , in 17th Italian-Korean Symposium on Relativistic Astrophysics , 3, 2022 [ 2203.03956 ].
- [75] A. Dhani, S. Borhanian, A. Gupta and B. Sathyaprakash, Cosmography with bright and Love sirens , arXiv e-prints (2022) [ 2212.13183 ].
- [76] Y.-J. Wang, J.-Z. Qi, B. Wang, J.-F. Zhang, J.-L. Cui and X. Zhang, Cosmological model-independent measurement of cosmic curvature using distance sum rule with the help of gravitational waves , Mon. Not. Roy. Astron. Soc. 516 (2022) 5187 [ 2201.12553 ].
- [77] L.-G. Zhu, L.-H. Xie, Y.-M. Hu, S. Liu, E.-K. Li, N.R. Napolitano et al., Constraining the Hubble constant to a precision of about 1% using multi-band dark standard siren detections , Sci. China Phys. Mech. Astron. 65 (2022) 259811 [ 2110.05224 ].
- [78] M. Califano, I. de Martino, D. Vernieri and S. Capozziello, Exploiting the Einstein Telescope to solve the Hubble tension , Phys. Rev. D 107 (2023) 123519 [ 2208.13999 ].
- [79] W.-T. Hou, J.-Z. Qi, T. Han, J.-F. Zhang, S. Cao and X. Zhang, Prospects for constraining interacting dark energy models from gravitational wave and gamma ray burst joint observation , JCAP 05 (2023) 017 [ 2211.10087 ].

- [80] P.-J. Wu, Y. Shao, S.-J. Jin and X. Zhang, A path to precision cosmology: synergy between four promising late-universe cosmological probes , JCAP 06 (2023) 052 [ 2202.09726 ].
- [81] T. Han, S.-J. Jin, J.-F. Zhang and X. Zhang, A comprehensive forecast for cosmological parameter estimation using joint observations of gravitational waves and short γ -ray bursts , Eur. Phys. J. C 84 (2024) 663 [ 2309.14965 ].
- [82] J.-G. Zhang, Z.-W. Zhao, Y. Li, J.-F. Zhang, D. Li and X. Zhang, Cosmology with fast radio bursts in the era of SKA , Sci. China Phys. Mech. Astron. 66 (2023) 120412 [ 2307.01605 ].
- [83] T.-N. Li, S.-J. Jin, H.-L. Li, J.-F. Zhang and X. Zhang, Prospects for Probing the Interaction between Dark Energy and Dark Matter Using Gravitational-wave Dark Sirens with Neutron Star Tidal Deformation , Astrophys. J. 963 (2024) 52 [ 2310.15879 ].
- [84] J. Yu, Z. Liu, X. Yang, Y. Wang, P. Zhang, X. Zhang et al., Measuring the Hubble Constant of Binary Neutron Star and Neutron Star-Black Hole Coalescences: Bright Sirens and Dark Sirens , Astrophys. J. Suppl. 270 (2024) 24 [ 2311.11588 ].
- [85] LIGO Scientific, Virgo, 1M2H, Dark Energy Camera GW-E, DES, DLT40, Las Cumbres Observatory, VINROUGE, MASTER collaboration, A gravitational-wave standard siren measurement of the Hubble constant , Nature 551 (2017) 85 [ 1710.05835 ].
- [86] LIGO Scientific, Virgo collaboration, GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral , Phys. Rev. Lett. 119 (2017) 161101 [ 1710.05832 ].
- [87] LIGO Scientific, Virgo, Fermi GBM, INTEGRAL, IceCube, AstroSat Cadmium Zinc Telluride Imager Team, IPN, Insight-Hxmt, ANTARES, Swift, AGILE Team, 1M2H Team, Dark Energy Camera GW-EM, DES, DLT40, GRAWITA, Fermi-LAT, ATCA, ASKAP, Las Cumbres Observatory Group, OzGrav, DWF (Deeper Wider Faster Program), AST3, CAASTRO, VINROUGE, MASTER, J-GEM, GROWTH, JAGWAR, CaltechNRAO, TTU-NRAO, NuSTAR, Pan-STARRS, MAXI Team, TZAC Consortium, KU, Nordic Optical Telescope, ePESSTO, GROND, Texas Tech University, SALT Group, TOROS, BOOTES, MWA, CALET, IKI-GW Follow-up, H.E.S.S., LOFAR, LWA, HAWC, Pierre Auger, ALMA, Euro VLBI Team, Pi of Sky, Chandra Team at McGill University, DFN, ATLAS Telescopes, High Time Resolution Universe Survey, RIMAS, RATIR, SKA South Africa/MeerKAT collaboration, Multi-messenger Observations of a Binary Neutron Star Merger , Astrophys. J. Lett. 848 (2017) L12 [ 1710.05833 ].
- [88] W. Del Pozzo, Inference of the cosmological parameters from gravitational waves: application to second generation interferometers , Phys. Rev. D 86 (2012) 043011 [ 1108.1317 ].
- [89] H.-Y. Chen, M. Fishbach and D.E. Holz, A two per cent Hubble constant measurement from standard sirens within five years , Nature 562 (2018) 545 [ 1712.06531 ].
- [90] R. Nair, S. Bose and T.D. Saini, Measuring the Hubble constant: Gravitational wave observations meet galaxy clustering , Phys. Rev. D 98 (2018) 023502 [ 1804.06085 ].
- [91] LIGO Scientific, Virgo collaboration, A Standard Siren Measurement of the Hubble Constant from GW170817 without the Electromagnetic Counterpart , Astrophys. J. Lett. 871 (2019) L13 [ 1807.05667 ].
- [92] R. Gray et al., Cosmological inference using gravitational wave standard sirens: A mock data analysis , Phys. Rev. D 101 (2020) 122001 [ 1908.06050 ].
- [93] DES, LIGO Scientific, Virgo collaboration, First Measurement of the Hubble Constant from a Dark Standard Siren using the Dark Energy Survey Galaxies and the LIGO/Virgo Binary-Black-hole Merger GW170814 , Astrophys. J. Lett. 876 (2019) L7 [ 1901.01540 ].
- [94] LIGO Scientific, Virgo, VIRGO collaboration, A Gravitational-wave Measurement of the

Hubble Constant Following the Second Observing Run of Advanced LIGO and Virgo , Astrophys. J. 909 (2021) 218 [ 1908.06060 ].

| [95]   | J. Yu, Y. Wang, W. Zhao and Y. Lu, Hunting for the host galaxy groups of binary black holes and the application in constraining Hubble constant , Mon. Not. Roy. Astron. Soc. 498 (2020) 1786 [ 2003.06586 ].                                                            |
|--------|--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| [96]   | DES collaboration, A statistical standard siren measurement of the Hubble constant from the LIGO/Virgo gravitational wave compact object merger GW190814 and Dark Energy Survey galaxies , Astrophys. J. Lett. 900 (2020) L33 [ 2006.14961 ].                            |
| [97]   | A. Palmese, C.R. Bom, S. Mucesh and W.G. Hartley, A Standard Siren Measurement of the Hubble Constant Using Gravitational-wave Events from the First Three LIGO/Virgo Observing Runs and the DESI Legacy Survey , Astrophys. J. 943 (2023) 56 [ 2111.06445 ].            |
| [98]   | LIGO Scientific, Virgo, KAGRA collaboration, Constraints on the Cosmic Expansion History from GWTC-3 , Astrophys. J. 949 (2023) 76 [ 2111.03604 ].                                                                                                                       |
| [99]   | H. Leandro, V. Marra and R. Sturani, Measuring the Hubble constant with black sirens , Phys. Rev. D 105 (2022) 023523 [ 2109.07537 ].                                                                                                                                    |
| [100]  | A. Finke, S. Foffa, F. Iacovelli, M. Maggiore and M. Mancarella, Cosmology with LIGO/Virgo dark sirens: Hubble parameter and modified gravitational wave propagation , JCAP 08 (2021) 026 [ 2101.12660 ].                                                                |
| [101]  | J.R. Gair et al., The Hitchhiker's Guide to the Galaxy Catalog Approach for Dark Siren Gravitational-wave Cosmology , Astron. J. 166 (2023) 22 [ 2212.08694 ].                                                                                                           |
| [102]  | T. Yang, R.-G. Cai and H.M. Lee, Space-borne atom interferometric gravitational wave detections. Part III. Eccentricity on dark sirens , JCAP 10 (2022) 061 [ 2208.10998 ].                                                                                              |
| [103]  | N. Muttoni, D. Laghi, N. Tamanini, S. Marsat and D. Izquierdo-Villalba, Dark siren cosmology with binary black holes in the era of third-generation gravitational wave detectors , Phys. Rev. D 108 (2023) 043543 [ 2303.10693 ].                                        |
| [104]  | J.-Y. Song, L.-F. Wang, Y. Li, Z.-W. Zhao, J.-F. Zhang, W. Zhao et al., Synergy between CSST galaxy survey and gravitational-wave observation: Inferring the Hubble constant from dark standard sirens , Sci. China Phys. Mech. Astron. 67 (2024) 230411 [ 2212.00531 ]. |
| [105]  | G. D´lya, a G. Galg´czi, o L. Dobos, Z. Frei, I.S. Heng, R. Macas et al., GLADE: A galaxy catalogue for multimessenger searches in the advanced gravitational-wave detector era , Mon. Not. Roy. Astron. Soc. 479 (2018) 2374 [ 1804.05709 ].                            |
| [106]  | G. D´lya a et al., GLADE+ : an extended galaxy catalogue for multimessenger searches with advanced gravitational-wave detectors , Mon. Not. Roy. Astron. Soc. 514 (2022) 1403 [ 2110.06184 ].                                                                            |
| [107]  | KAGRA collaboration, Detector configuration of KAGRA: The Japanese cryogenic gravitational-wave detector , Class. Quant. Grav. 29 (2012) 124007 [ 1111.7185 ].                                                                                                           |
| [108]  | LIGO Scientific collaboration, Advanced LIGO , Class. Quant. Grav. 32 (2015) 074001 [ 1411.4547 ].                                                                                                                                                                       |
| [109]  | VIRGO collaboration, Advanced Virgo: a second-generation interferometric gravitational wave detector , Class. Quant. Grav. 32 (2015) 024001 [ 1408.3978 ].                                                                                                               |
| [110]  | LIGO Scientific collaboration, Exploring the Sensitivity of Next Generation Gravitational Wave Detectors , Class. Quant. Grav. 34 (2017) 044001 [ 1607.08697 ].                                                                                                          |
| [111]  | M. Punturo et al., The Einstein Telescope: A third-generation gravitational wave observatory , Class. Quant. Grav. 27 (2010) 194002.                                                                                                                                     |

- [112] N. Seto, S. Kawamura and T. Nakamura, Possibility of direct measurement of the acceleration of the universe using 0.1-Hz band laser interferometer gravitational wave antenna in space , Phys. Rev. Lett. 87 (2001) 221103 [ astro-ph/0108011 ].
- [113] S. Kawamura et al., Space gravitational-wave antennas DECIGO and B-DECIGO , Int. J. Mod. Phys. D 28 (2019) 1845001.
- [114] Y.-Y. Dong, J.-Y. Song, S.-J. Jin, J.-F. Zhang and X. Zhang, Enhancing dark siren cosmology through multi-band gravitational wave synergetic observations , arXiv e-prints (2024) [ 2404.18188 ].
- [115] LGWA collaboration, Lunar Gravitational-wave Antenna , Astrophys. J. 910 (2021) 1 [ 2010.13726 ].
- [116] K. Jani and A. Loeb, Gravitational-wave Lunar Observatory for Cosmology , JCAP 2021 (2021) 044 [ 2007.08550 ].
- [117] LISA collaboration, Laser Interferometer Space Antenna , arXiv e-prints (2017) [ 1702.00786 ].
- [118] T. Robson, N.J. Cornish and C. Liu, The construction and use of LISA sensitivity curves , Class. Quant. Grav. 36 (2019) 105011 [ 1803.01944 ].
- [119] LISA Cosmology Working Group collaboration, Cosmology with the Laser Interferometer Space Antenna , Living Rev. Rel. 26 (2023) 5 [ 2204.05434 ].
- [120] Y.-L. Wu, Hyperunified field theory and Taiji program in space for GWD , Int. J. Mod. Phys. A 33 (2018) 1844014 [ 1805.10119 ].
- [121] W.-H. Ruan, Z.-K. Guo, R.-G. Cai and Y.-Z. Zhang, Taiji program: Gravitational-wave sources , Int. J. Mod. Phys. A 35 (2020) 2050075 [ 1807.09495 ].
- [122] W.-R. Hu and Y.-L. Wu, The Taiji Program in Space for gravitational wave physics and the nature of gravity , Natl. Sci. Rev. 4 (2017) 685.
- [123] Z.-W. Zhao, L.-F. Wang, J.-F. Zhang and X. Zhang, Prospects for improving cosmological parameter estimation with gravitational-wave standard sirens from Taiji , Sci. Bull. 65 (2020) 1340 [ 1912.11629 ].
- [124] TianQin collaboration, TianQin: a space-borne gravitational wave detector , Class. Quant. Grav. 33 (2016) 035010 [ 1512.02076 ].
- [125] H.-T. Wang et al., Science with the TianQin observatory: Preliminary results on massive black hole binaries , Phys. Rev. D 100 (2019) 043003 [ 1902.04423 ].
- [126] L.-F. Wang, Z.-W. Zhao, J.-F. Zhang and X. Zhang, A preliminary forecast for cosmological parameter estimation with gravitational-wave standard sirens from TianQin , JCAP 11 (2020) 012 [ 1907.01838 ].
- [127] J. Luo et al., The first round result from the TianQin-1 satellite , Class. Quant. Grav. 37 (2020) 185013 [ 2008.09534 ].
- [128] V.K. Milyukov, TianQin Space-Based Gravitational Wave Detector: Key Technologies and Current State of Implementation , Astron. Rep. 64 (2020) 1067.
- [129] TianQin collaboration, The TianQin project: current progress on science and technology , PTEP 2021 (2021) 05A107 [ 2008.10332 ].
- [130] NANOGrav collaboration, The NANOGrav 15 yr Data Set: Evidence for a Gravitational-wave Background , Astrophys. J. Lett. 951 (2023) L8 [ 2306.16213 ].
- [131] NANOGrav collaboration, The NANOGrav 15 yr Data Set: Observations and Timing of 68 Millisecond Pulsars , Astrophys. J. Lett. 951 (2023) L9 [ 2306.16217 ].

- [132] EPTA collaboration, The second data release from the European Pulsar Timing Array III. Search for gravitational wave signals , Astron. Astrophys. 678 (2023) A50 [ 2306.16214 ].
- [133] EPTA collaboration, The second data release from the European Pulsar Timing Array - I. The dataset and timing analysis , Astron. Astrophys. 678 (2023) A48 [ 2306.16224 ].
- [134] A. Zic et al., The Parkes Pulsar Timing Array third data release , Publ. Astron. Soc. Austral. 40 (2023) e049 [ 2306.16230 ].
- [135] D.J. Reardon et al., Search for an Isotropic Gravitational-wave Background with the Parkes Pulsar Timing Array , Astrophys. J. Lett. 951 (2023) L6 [ 2306.16215 ].
- [136] H. Xu et al., Searching for the Nano-Hertz Stochastic Gravitational Wave Background with the Chinese Pulsar Timing Array Data Release I , Res. Astron. Astrophys. 23 (2023) 075024 [ 2306.16216 ].
- [137] R.w. Hellings and G.s. Downs, UPPER LIMITS ON THE ISOTROPIC GRAVITATIONAL RADIATION BACKGROUND FROM PULSAR TIMING ANALYSIS Astrophys. J. Lett. , 265 (1983) L39.
- [138] A. Sesana, A. Vecchio and M. Volonteri, Gravitational waves from resolvable massive black hole binary systems and observations with Pulsar Timing Arrays , Mon. Not. Roy. Astron. Soc. 394 (2009) 2255 [ 0809.3412 ].
- [139] V. Corbin and N.J. Cornish, Pulsar Timing Array Observations of Massive Black Hole Binaries , arXiv e-prints (2010) [ 1008.1782 ].
- [140] A. Sesana and A. Vecchio, Measuring the parameters of massive black hole binary systems with Pulsar Timing Array observations of gravitational waves , Phys. Rev. D 81 (2010) 104008 [ 1003.0677 ].
- [141] S. Babak and A. Sesana, Resolving multiple supermassive black hole binaries with pulsar timing arrays , Phys. Rev. D 85 (2012) 044034 [ 1112.1075 ].
- [142] K.J. Lee et al., Gravitational wave astronomy of single sources with a pulsar timing array , Mon. Not. Roy. Astron. Soc. 414 (2011) 3251 [ 1103.0115 ].
- [143] J.A. Ellis, X. Siemens and J.D.E. Creighton, Optimal strategies for continuous gravitational wave detection in pulsar timing arrays , Astrophys. J. 756 (2012) 175 [ 1204.4218 ].
- [144] J.A. Ellis, A Bayesian analysis pipeline for continuous GW sources in the PTA band , Class. Quant. Grav. 30 (2013) 224004 [ 1305.0835 ].
- [145] X.-J. Zhu, L. Wen, G. Hobbs, Y. Zhang, Y. Wang, D.R. Madison et al., Detection and localization of single-source gravitational waves with pulsar timing arrays , Mon. Not. Roy. Astron. Soc. 449 (2015) 1650 [ 1502.06001 ].
- [146] Y. Wang and S.D. Mohanty, Pulsar Timing Array Based Search for Supermassive Black Hole Binaries in the Square Kilometer Array Era , Phys. Rev. Lett. 118 (2017) 151104 [ 1611.09440 ].
- [147] Y. Feng, D. Li, Z. Zheng and C.-W. Tsai, Supermassive Binary Black Hole Evolution can be traced by a small SKA Pulsar Timing Array , Phys. Rev. D 102 (2020) 023014 [ 2005.11118 ].
- [148] C. Yan, W. Zhao and Y. Lu, On Using Inspiraling Supermassive Binary Black Holes in the PTA Frequency Band as Standard Sirens to Constrain Dark Energy , The Astrophysical Journal 889 (2020) [ 1912.04103 ].
- [149] L.-F. Wang, Y. Shao, J.-F. Zhang and X. Zhang, Ultra-low-frequency gravitational waves from individual supermassive black hole binaries as standard sirens , arXiv e-prints (2022) [ 2201.00607 ].

| [150]   | C.M.F. Mingarelli, T.J.W. Lazio, A. Sesana, J.E. Greene, J.A. Ellis, C.-P. Ma et al., The Local Nanohertz Gravitational-Wave Landscape From Supermassive Black Hole Binaries , Nature Astron. 1 (2017) 886 [ 1708.03491 ].                                 |
|---------|------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| [151]   | J.A. Casey-Clyde, C.M.F. Mingarelli, J.E. Greene, K. Pardo, M. Na˜ez n and A.D. Goulding, A Quasar-based Supermassive Black Hole Binary Population Model: Implications for the Gravitational Wave Background , Astrophys. J. 924 (2022) 93 [ 2107.11390 ]. |
| [152]   | T.H. Jarrett, T. Chester, R. Cutri, S. Schneider, M. Skrutskie and J.P. Huchra, 2mass extended source catalog: overview and algorithms , Astron. J. 119 (2000) 2498 [ astro-ph/0004318 ].                                                                  |
| [153]   | P.A. Rosado and A. Sesana, Targeting supermassive black hole binaries and gravitational wave sources for the pulsar timing array , Mon. Not. Roy. Astron. Soc. 439 (2014) 3986 [ 1311.0883 ].                                                              |
| [154]   | NANOGrav collaboration, The NANOGrav 11 yr Data Set: Limits on Supermassive Black Hole Binaries in Galaxies within 500 Mpc , Astrophys. J. 914 (2021) 121 [ 2101.02716 ].                                                                                  |
| [155]   | J.P. Huchra et al., The 2MASS Redshift Survey - Description and Data Release , Astrophys. J. Suppl. 199 (2012) 26 [ 1108.0669 ].                                                                                                                           |
| [156]   | A. Sesana, Z. Haiman, B. Kocsis and L.Z. Kelley, Testing the binary hypothesis: pulsar timing constraints on supermassive black hole binary candidates , Astrophys. J. 856 (2018) 42 [ 1703.10611 ].                                                       |
| [157]   | M. Cappellari, Effect of Environment on Galaxies' Mass-Size Distribution: Unveiling the Transition from outside-in to inside-out Evolution , Astrophys. J. Lett. 778 (2013) L2 [ 1309.1136 ].                                                              |
| [158]   | N.J. McConnell and C.-P. Ma, Revisiting the Scaling Relations of Black Hole Masses and Host Galaxy Properties , Astrophys. J. 764 (2013) 184 [ 1211.2816 ].                                                                                                |
| [159]   | J. Binney and S. Tremaine, Galactic Dynamics: Second Edition , Princeton University Press (2008).                                                                                                                                                          |
| [160]   | A. Sesana and F.M. Khan, Scattering experiments meet N-body - I. A practical recipe for the evolution of massive black hole binaries in stellar environments , Mon. Not. Roy. Astron. Soc. 454 (2015) L66 [ 1505.02062 ].                                  |
| [161]   | S. Genel, M. Vogelsberger, V. Springel, D. Sijacki, D. Nelson, G. Snyder et al., Introducing the Illustris Project: the evolution of galaxy populations across cosmic time , Mon. Not. Roy. Astron. Soc. 445 (2014) 175 [ 1405.3749 ].                     |
| [162]   | V. Rodriguez-Gomez et al., The merger rate of galaxies in the Illustris Simulation: a comparison with observations and semi-empirical models , Mon. Not. Roy. Astron. Soc. 449 (2015) 49 [ 1502.01339 ].                                                   |
| [163]   | Planck collaboration, Planck 2015 results. XIII. Cosmological parameters , Astron. Astrophys. 594 (2016) A13 [ 1502.01589 ].                                                                                                                               |
| [164]   | E.S. Phinney, A Practical Theorem on Gravitational Wave Backgrounds , arXiv e-prints (2001) astro.                                                                                                                                                         |
| [165]   | A. Sesana, Systematic investigation of the expected gravitational wave signal from supermassive black hole binaries in the pulsar timing band , Mon. Not. Roy. Astron. Soc. 433 (2013) 1 [ 1211.5375 ].                                                    |
| [166]   | K.S. Thorne, Gravitational radiation. , Three Hundred Years of Gravitation (1987) 330.                                                                                                                                                                     |
| [167]   | P.F. Hopkins, L. Hernquist, T.J. Cox, T. Di Matteo, B. Robertson and V. Springel, A Unified, merger-driven model for the origin of starbursts, quasars, the cosmic x-ray                                                                                   |

background, supermassive black holes and galaxy spheroids , Astrophys. J. Suppl. 163 (2006) 1 [ astro-ph/0506398 ].

- [168] P.F. Hopkins, G.T. Richards and L. Hernquist, An Observational Determination of the Bolometric Quasar Luminosity Function , Astrophys. J. 654 (2007) 731 [ astro-ph/0605678 ].
- [169] S. Vagnozzi, Implications of the NANOGrav results for inflation , Mon. Not. Roy. Astron. Soc. 502 (2021) L11 [ 2009.13432 ].
- [170] R.-G. Cai, Z.-K. Guo, J. Liu, L. Liu and X.-Y. Yang, Primordial black holes and gravitational waves from parametric amplification of curvature perturbations , JCAP 06 (2020) 013 [ 1912.10437 ].
- [171] V. De Luca, G. Franciolini and A. Riotto, NANOGrav Data Hints at Primordial Black Holes as Dark Matter , Phys. Rev. Lett. 126 (2021) 041303 [ 2009.08268 ].
- [172] P. Schwaller, Gravitational Waves from a Dark Phase Transition , Phys. Rev. Lett. 115 (2015) 181101 [ 1504.07263 ].
- [173] NANOGrav collaboration, Searching for Gravitational Waves from Cosmological Phase Transitions with the NANOGrav 12.5-Year Dataset , Phys. Rev. Lett. 127 (2021) 251302 [ 2104.13930 ].
- [174] J. Ellis and M. Lewicki, Cosmic String Interpretation of NANOGrav Pulsar Timing Data , Phys. Rev. Lett. 126 (2021) 041304 [ 2009.06555 ].
- [175] X. Siemens, V. Mandic and J. Creighton, Gravitational wave stochastic background from cosmic (super)strings , Phys. Rev. Lett. 98 (2007) 111101 [ astro-ph/0610920 ].
- [176] R.Z. Ferreira, A. Notari, O. Pujolas and F. Rompineve, Gravitational waves from domain walls in Pulsar Timing Array datasets , JCAP 02 (2023) 001 [ 2204.04228 ].
- [177] P.C. Peters, Gravitational Radiation and the Motion of Two Point Masses , Phys. Rev. 136 (1964) B1224.
- [178] T. Liu, T. Cohen, C. McGrath, P.B. Demorest and S.J. Vigeland, Multi-messenger Approaches to Supermassive Black Hole Binary Detection and Parameter Estimation. II. Optimal Strategies for a Pulsar Timing Array , Astrophys. J. 945 (2023) 78 [ 2301.07135 ].
- [179] A. Albrecht et al., Findings of the Joint Dark Energy Mission Figure of Merit Science Working Group , arXiv e-prints (2009) [ 0901.0721 ].
- [180] C.M. Hirata, D.E. Holz and C. Cutler, Reducing the weak lensing noise for the gravitational wave Hubble diagram using the non-Gaussianity of the magnification distribution , Phys. Rev. D 81 (2010) 124046 [ 1004.3988 ].
- [181] N. Tamanini, C. Caprini, E. Barausse, A. Sesana, A. Klein and A. Petiteau, Science with the space-based interferometer eLISA. III: Probing the expansion of the Universe using gravitational wave standard sirens , JCAP 04 (2016) 002 [ 1601.07112 ].
- [182] R.N. Manchester, G.B. Hobbs, A. Teoh and M. Hobbs, The Australia Telescope National Facility pulsar catalogue , Astron. J. 129 (2005) 1993 [ astro-ph/0412641 ].
- [183] D.W. Hogg, Distance measures in cosmology , arXiv e-prints (1999) [ astro-ph/9905116 ].
- [184] B. Henriques, S. White, G. Lemson, P. Thomas, Q. Guo, G.-D. Marleau et al., Confronting theoretical models with the observed evolution of the galaxy population out to z=4 , Mon. Not. Roy. Astron. Soc. 421 (2012) 2904 [ 1109.3457 ].
- [185] M. Chevallier and D. Polarski, Accelerating universes with scaling dark matter , Int. J. Mod. Phys. D 10 (2001) 213 [ gr-qc/0009008 ].
- [186] E.V. Linder, Exploring the expansion history of the universe , Phys. Rev. Lett. 90 (2003) 091301 [ astro-ph/0208512 ].

- [187] L. Amendola, Coupled quintessence , Phys. Rev. D 62 (2000) 043511 [ astro-ph/9908023 ].
- [188] M.G. Brown, K. Freese and W.H. Kinney, The Phantom bounce: A New oscillating cosmology , JCAP 03 (2008) 002 [ astro-ph/0405353 ].
- [189] X. Zhang, Can black holes be torn up by phantom dark energy in cyclic cosmology? , Eur. Phys. J. C 60 (2009) 661 [ 0708.1408 ].
- [190] A. Ibarra and D. Tran, Decaying Dark Matter and the PAMELA Anomaly , JCAP 02 (2009) 021 [ 0811.1555 ].
- [191] D.-Z. He, X.-J. Bi, S.-J. Lin, P.-F. Yin and X. Zhang, Expected LHAASO sensitivity to decaying dark matter signatures from dwarf galaxies gamma-ray emission , Chin. Phys. C 44 (2020) 085001 [ 1910.05017 ].
- [192] E. Di Valentino, A. Melchiorri and J. Silk, Planck evidence for a closed Universe and a possible crisis for cosmology , Nature Astron. 4 (2019) 196 [ 1911.02087 ].
- [193] J.-Z. Qi, J.-W. Zhao, S. Cao, M. Biesiada and Y. Liu, Measurements of the Hubble constant and cosmic curvature with quasars: ultracompact radio structure and strong gravitational lensing , Mon. Not. Roy. Astron. Soc. 503 (2021) 2179 [ 2011.00713 ].

Figure 7 . The marginalized posterior of H 0 using the mock dark sirens data based on the M17 model (upper panel) and the C21 model (lower panel) in different PTA detection cases. In each panel, PTA detection cases are considered with N p of 20, 100, and 200, observed over ten years with σ t = 20 ns. The colored lines show the posterior probability density of H 0 for each source with colors representing their SNRs, and the thick black line shows the combined posterior probability density of H 0 . The fiducial value of H 0 is shown as a vertical dashed line.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00609v2-figure-b1fee73c-2cf9-4da8-a531-2d3d59db6d76.jpg)
Si-Ren Xiao Yue Shao Ling-Feng Wang Ji-Yu Song Lu Feng Jing-Fei Zhang Xin Zhang
Nanohertz gravitational waves from a quasar-based supermassive black hole binary population model as dark sirens
Recently, several pulsar timing array (PTA) projects have detected evidence
of the existence of a stochastic gravitational wave background (SGWB) in the
nanohertz frequency band, providing confidence in detecting individual
supermassive black hole binaries (SMBHBs) in the future. Nanohertz GWs emitted
by inspiraling SMBHBs encode the luminosity distances of SMBHBs. They can serve
as dark sirens to explore the cosmic expansion history via a statistical method
to obtain the redshift information of GW sources' host galaxies using galaxy
catalogs. The theoretical analysis of the dark siren method relies on the
modeling of the population of SMBHBs. Using a population model consistent with
the latest SGWB observations is essential, as the SGWB provides significant
information about the distribution of SMBHBs. In this work, we employ a
quasar-based model, which can self-consistently account for the SGWB amplitude,
to estimate the population of SMBHBs. We constrain the Hubble constant using
the mock GW data from different detection cases of PTAs in the future. Our
results show that a PTA consisting of 100 pulsars with a white noise level of
20 ns could measure the Hubble constant with a precision close to $1\%$ over a
10-year observation period, and a PTA with 200 pulsars may achieve this goal
over a 5-year observation period. The results indicate that modeling the SMBHB
population significantly influences the analysis of dark sirens, and SMBHB dark
sirens have the potential to be developed as a valuable cosmological probe.