# Implications for the Supermassive Black Hole Binaries from the NANOGrav
15-year Data Set

Yan-Chen Bi<EMAIL_ADDRESS>CAS Key Laboratory of Theoretical Physics,
Institute of Theoretical Physics, Chinese Academy of Sciences,Beijing 100190,
China  School of Physical Sciences, University of Chinese Academy of Sciences,
No. 19A Yuquan Road, Beijing 100049, China  Yu-Mei Wu<EMAIL_ADDRESS>School
of Fundamental Physics and Mathematical Sciences, Hangzhou Institute for
Advanced Study, UCAS, Hangzhou 310024, China  School of Physical Sciences,
University of Chinese Academy of Sciences, No. 19A Yuquan Road, Beijing
100049, China  Zu-Cheng Chen<EMAIL_ADDRESS>Department of Astronomy,
Beijing Normal University, Beijing 100875, China  Advanced Institute of
Natural Sciences, Beijing Normal University, Zhuhai 519087, China  Department
of Physics and Synergistic Innovation Center for Quantum Effects and
Applications, Hunan Normal University, Changsha, Hunan 410081, China  Qing-Guo
Huang<EMAIL_ADDRESS>School of Fundamental Physics and Mathematical
Sciences, Hangzhou Institute for Advanced Study, UCAS, Hangzhou 310024, China
<EMAIL_ADDRESS>CAS Key Laboratory of Theoretical Physics, Institute of
Theoretical Physics, Chinese Academy of Sciences,Beijing 100190, China  School
of Physical Sciences, University of Chinese Academy of Sciences, No. 19A
Yuquan Road, Beijing 100049, China

###### Abstract

NANOGrav, EPTA, PPTA, and CPTA have announced the evidence for a stochastic
signal from their latest data sets. Supermassive black hole binaries (SMBHBs)
are supposed to be the most promising gravitational-wave (GW) sources of
pulsar timing arrays. Assuming an astro-informed formation model, we use the
NANOGrav 15-year data set to constrain the gravitational wave background (GWB)
from SMBHBs. Our results prefer a large turn-over eccentricity of the SMBHB
orbit when GWs begin to dominate the SMBHBs evolution. Furthermore, the GWB
spectrum is extrapolated to the space-borne GW detector frequency band by
including inspiral-merge-cutoff phases of SMBHBs and should be detected by
LISA, Taiji and TianQin in the near future.

Introduction. Supermassive black holes (SMBHs), with masses from
$10^{5}M_{\odot}$ to $10^{11}M_{\odot}$, are thought to reside in the centers
of nearly all galaxies (Akiyama  _et al._ , 2019; Kormendy and Ho, 2013a). The
details of their formation, evolution and interaction with host galaxies need
to be clarified. In the scenario of galaxies coalescence, SMBHs hosted in
their nuclei sink to the center of the remnant due to the interaction with the
surrounding ambient and eventually form bound supermassive black hole binaries
(SMBHBs) (Begelman  _et al._ , 1980). The SMBHBs subsequently harden because
of dynamical interaction with the dense background (Dotti  _et al._ , 2012a)
until the gravitational wave (GW) takes over at sub-parsec separation, forming
an abundant population GW sources in around nHz frequencies (Chen  _et al._ ,
2019a). Moreover, the interaction with the stellar environment and scattering
of ambient stars potentially attenuate the GW spectrum at the lower end of PTA
ranges and tend to increase the binary eccentricity (Kocsis and Sesana, 2011;
Sesana, 2013; Ravi  _et al._ , 2014; Kelley  _et al._ , 2017; Chen  _et al._ ,
2017; Rodig  _et al._ , 2011; Quinlan, 1996; Armitage and Natarajan, 2005),
which leaves imprints on the GW emission.

Therefore, detecting or constraining such signals will help to reveal the
nature and properties of SMBHBs (for a review see (Sesana, 2013)). Since
galaxies are observed to merge quite frequently (Bell  _et al._ , 2006; de
Ravel  _et al._ , 2009) and the observable Universe encompasses several
billions of them, a largish cosmological population of SMBHBs is expected to
produce the stochastic gravitational-wave background (SGWB) (Sesana  _et al._
, 2008). The SGWB is going to be captured by pulsar timing arrays (PTAs)
(Foster and Backer, 1990). SGWBs affect the time of arrivals (TOAs) of radio
pulsar from millisecond pulsars (MSPs) in a correlated way (Hellings and
Downs, 1983). PTAs search for SGWB using the Hellings $\&$ Downs curve.
Besides detecting the SGWB at nano-Hz frequency, the final goal of PTAs is to
extract useful astrophysical information from their data.

Recently, the North American Nanohertz Observatory for Gravitational Waves
(NANOGrav (Agazie  _et al._ , 2023a, b; Afzal  _et al._ , 2023)), the European
Pulsar Timing Array (EPTA (Antoniadis  _et al._ , 2023a, b)), the Parkes
Pulsar Timing Array (PPTA (Zic  _et al._ , 2023; Reardon  _et al._ , 2023)),
and the Chinese Pulsar Timing Array (CPTA (Nan  _et al._ , 2011; Xu  _et al._
, 2023)) have announced the evidence for a stochastic signal consistent with a
SGWB, brings us a great opportunity to prove the imprints of SMBHBs (Agazie
_et al._ , 2023c, d; Antoniadis  _et al._ , 2023c).

SMBHBs are also one of the most promising GW sources for space-borne GW
detectors, such as LISA (Amaro-Seoane  _et al._ , 2017), Taiji (Luo  _et al._
, 2020) and TianQin (Luo  _et al._ , 2016), in the frequency band around
$10^{-4}\sim 10^{-1}$ Hz. However, the expected event rate of SMBHB sources in
space-borne detectors is still unknown. During the SMBHBs evolution, the
inspiral phase falls into the range of PTA while the merger phase in the
sensitive interval of space-borne GW observatories (Chen  _et al._ , 2019a;
Amaro-Seoane  _et al._ , 2017). This insight indicates that the observation of
PTA and space-borne GW detector are complementary (Steinle  _et al._ , 2023).
With the joint observation, the complete description of SMBHBs will be well
constrained.

In this letter, we use the NANOGrav 15-year data set to constrain the GWB from
SMBHBs, and find that the SMBHBs orbits have a large eccentricity when GWs
begin to dominate the SMBHBs evolution. In addition, we extrapolate the GWB
spectrum from $10^{-9}$ Hz to $10^{-1}$ Hz (see Fig. 1), providing a promising
GW source for the space-borne GW detectors.

Figure 1:  The expected SGWB spectrum $\Omega_{\rm{GW}}(f)$ from the
astrophysical SMBHBs in the frequency range $[10^{-10},10^{-1}]\rm{Hz}$. The
light-blue violin plots show the posterior of the first 14-frequency bins for
the NANOGrav 15-yr data set. The green line shows the medium spectrum obtained
from the NANOGrav analysis based on the astro-informed model, and the light-
green shaded region indicates the $90\%$ credible interval. We also show the
expected power-law integrated (PI) curves of LISA, TianQin and Taiji which
represent their detection abilities with solid, dashed, dash-dotted gray
lines, respectively.

SGWB from SMBHB. SGWB from SMBHBs is the most promising target of PTA. In
general SMBHB is supposed to form following the merge of two galaxies. After
galaxies merging, their central SMBHs sink into the center of the merger
remnant and form a bound binary. Initially the binary orbit shrinks due to the
energy and angular momentum exchanging with surrounding stars and cold gas.
Then GW radiation will dominate the evolution at turn-over frequency $f_{t}$
(corresponding to initial eccentricity $e_{0}$), bringing the binary to final
coalescence (Dotti  _et al._ , 2012b). The spectrum of SGWB is composed from
the sum of all SMBHBs emitting at a given observed frequency $f$. The present-
day energy density of SGWB $\Omega_{\mathrm{GW}}(f)$ are given by (Phinney,
2001)

$\Omega_{\mathrm{GW}}(f)=\frac{8\pi Gf}{3H_{0}^{2}c^{2}}\int
dzd\mathcal{M}\frac{d^{2}n}{dzd\mathcal{M}}\frac{dE_{\rm GW}}{df_{r}}\ ,$ (1)

where $H_{0}=67.4\ {\rm km\ sec^{-1}\ Mpc^{-1}}$ (Aghanim  _et al._ , 2020) is
the Hubble constant, $f_{r}=(1+z)f$ is the source-frame GW frequency.
$\mathcal{M}=Mq^{3/5}/(1+q)^{1/5}$ is the chirp mass of SMBHB, where $M$ is
the primary SMBH mass and $q$ is the SMBHB mass ratio. Here,
$d^{2}n/dzd\mathcal{M}$ is the SMBHB population and ${dE_{\rm GW}}/{df_{r}}$
is the energy spectrum of single SMBHB. The relevant ranges in the integrals
used here are $0\leq z\leq 5$ and $10^{5}M_{\odot}\leq\mathcal{M}\leq
10^{11}M_{\odot}$.

The SMBHB population ${d^{2}n}/{dzd\mathcal{M}}$ in Eq. (1) are estimated from
astrophysical observation. It is widely acknowledged that nearly every galaxy
has a SMBH at their centers. This relation can be used as an astrophysical
model to describe the merger rate. In this astro-informed formation model, the
astronomical surveys serve relatively strong constrains naturally and SMBHB
population is allowed to interact with the environment. The population can be
expressed as (Sesana  _et al._ , 2009)

$\displaystyle\frac{d^{2}n}{dz^{\prime}d\mathcal{M}}=\int\frac{d^{3}n_{G}}{dz^{\prime}dM_{G}dq_{G}}\frac{dM_{G}}{dM}\frac{dq_{G}}{dq}\frac{dM}{d\mathcal{M}}dq\
,$ (2)

where $M_{G}$ is the primary galaxy mass and $q_{G}$ is the mass ratio between
the two paired galaxy, $q=q_{G}^{\alpha_{*}}$ is the SMBHB mass ratio
mentioned before. The relations used for transforming the galaxy mass $M_{G}$
into the primary black hole mass $M$ are given as (Sesana, 2013; Kormendy and
Ho, 2013b)

$\frac{M_{\rm
bulge}}{M_{G}}=\left\\{\begin{array}[]{lr}\frac{\sqrt{6.9}}{(\log
M_{G}-10)^{1.5}}\exp\left(\frac{-3.45}{\log M_{G}-10}\right)-0.615&\quad{\rm
if}\log M_{G}\geq 10\\\ 0.615&\quad{\rm if}\log M_{G}\leq
10\end{array}\right.$ (3) $\displaystyle
M=\mathcal{N}\left\\{M_{*}\left(\frac{M_{\rm
bulge}}{10^{11}M_{\odot}}\right),\epsilon\right\\}\ ,$ (4)

where $\mathcal{N}\\{x,y\\}$ is a log normal distribution with mean $x$ and
standard deviation $y$. $\\{M_{*},\alpha_{*},\epsilon\\}$ are the model
parameters for translation.

The galaxy differential merger rate per unit redshift, galaxy mass and mass
ratio ${d^{3}n_{G}}/{dz^{\prime}dM_{G}dq_{G}}$ can be written as (Sesana  _et
al._ , 2009)

$\frac{d^{3}n_{G}}{dz^{\prime}dM_{G}dq_{G}}=\frac{\Phi(M_{G},z)}{M_{G}\ln
10}\frac{\mathcal{F}(M_{G},z,q_{G})}{\tau(M_{G},z,q_{G})}\frac{dt_{r}}{dz}\ .$
(5)

where $z$ stands for the redshift of galaxy pair and $z^{\prime}$ is the
redshift at galaxy pair merging . The quantitative relation between $z$ and
$z^{\prime}$ is addressed in (Chen  _et al._ , 2019a). $dt_{r}/dz$ is the
relationship between time and redshift assuming a $\Lambda$CDM flat Universe
as follow

$\displaystyle\frac{dt_{r}}{dz}=\frac{1}{H_{0}(1+z)\sqrt{(\Omega_{M}(1+z)^{3}+\Omega_{k}(1+z)^{2}+\Omega_{\Lambda})}}\
.$

$\Phi(M_{G},z)={dn_{G}}/{d\log_{10}M_{G}}$ is the galaxy mass function (GSMF)
measured at redshift of galaxy pair $z$, the explicit expression is (Mortlock
_et al._ , 2015; Conselice  _et al._ , 2016)

$\displaystyle\Phi(M_{G},z)=10^{\Phi(z)}\ln
10\left(\frac{M_{G}}{M_{G0}}\right)^{1+\alpha(z)}\exp{\left(-\frac{M_{G}}{M_{G0}}\right)},$
(7)

where the parameters are
$\Phi(z)=\Phi_{0}+z\Phi_{I},\alpha(z)=\alpha_{0}+z\alpha_{I}$.
$\\{\Phi_{0},\Phi_{I},\alpha_{0},\alpha_{I},M_{G0}\\}$ are the five model
parameters for GSMF. $\mathcal{F}(M_{G},z,q_{G})$ is the differential pair
fraction with respect to the mass ratio of galaxy pair $q_{G}$ and is written
as (Chen  _et al._ , 2019a)

$\displaystyle\mathcal{F}(M_{G},z,q_{G})=\frac{df_{\rm
pair}}{dq_{G}}=f_{0}^{\prime}\left(\frac{M_{G}}{aM_{G0}}\right)^{\alpha_{f}}\left(1+z\right)^{\beta_{f}}q_{G}^{\gamma_{f}}\
.$ (8)

where $aM_{G0}=10^{11}M_{\odot}$ is an arbitrary reference mass (Chen  _et
al._ , 2019a). $\\{f_{0}^{\prime},\alpha_{f},\beta_{f},\gamma_{f}\\}$ are the
four model parameters for pair function. $\tau(M_{G},z,q_{G})$ is the merger
timescale of the pair galaxy and can be expressed as (Snyder  _et al._ , 2017)

$\displaystyle\tau(M_{G},z,q_{G})=\tau_{0}\left(\frac{M_{G}}{bM_{G0}}\right)^{\alpha_{\tau}}\left(1+z\right)^{\beta_{\tau}}q_{G}^{\gamma_{\tau}}\
.$ (9)

where $bM_{G0}=0.4/h_{0}\times 10^{11}M_{\odot}$ is an arbitrary reference
mass (Chen  _et al._ , 2019a).
$\\{\tau_{0},\alpha_{\tau},\beta_{\tau},\gamma_{\tau}\\}$ are the four model
parameters for merger timescale.

To sum up, the present-day energy density $\Omega_{\mathrm{GW}}(f)$ of SGWB
from SMBHB in galaxy model is therefore fully specified by a set of eighteen
model parameters: $\\{\Phi_{0},\Phi_{I},M_{G0},\alpha_{0},\alpha_{I}\\}$ for
the GSMF, $\\{f_{0}^{\prime},\alpha_{f},\beta_{f},\gamma_{f}\\}$ for the pair
fraction, $\\{\tau_{0},\alpha_{\tau},\beta_{\tau},\gamma_{\tau}\\}$ for the
merger timescale, $\\{M_{*},\alpha_{*},\epsilon\\}$ for galaxy-SMBH
transforming relation and $\\{e_{0},\zeta_{0}\\}$ mentioned below for the
SMBHB energy spectrum. The detailed descriptions of parameters are addressed
in Table.3.

The energy spectrum of single SMBHB $dE_{\rm GW}/df_{r}$ are calculated using
its self-similarity. So spectrum in any configuration can be obtained from a
reference spectrum by shift and re-scaling. The fiducial redshift and chirp
mass of reference spectrum are set as $z_{0}=0.02$ and
$\mathcal{M}_{0}=4.16\times 10^{8}M_{\odot}$, respectively. The basic idea of
so-called shift and re-scaling could be found in (Chen  _et al._ , 2017).

Table 1:  Parameters used in Eq. (10), Eq. (11) and Eq. (Implications for the Supermassive Black Hole Binaries from the NANOGrav 15-year Data Set) for a complete description of the energy spectrum (Chen  _et al._ , 2019b). Parameter | $a_{k}$ | $b_{k}$ | $c_{k}$
---|---|---|---
$\nu_{1}$ | $2.9740\times 10^{-1}$ | $4.4810\times 10^{-2}$ | $9.5560\times 10^{-2}$
$\nu_{2}$ | $5.9411\times 10^{-1}$ | $8.9794\times 10^{-2}$ | $1.9111\times 10^{-1}$
$\nu_{3}$ | $8.4845\times 10^{-1}$ | $1.2848\times 10^{-1}$ | $2.7299\times 10^{-1}$
$\sigma_{\nu}$ | $5.0801\times 10^{-1}$ | $7.7515\times 10^{-2}$ | $2.2369\times 10^{-2}$

Past calculations of SMBHBs SGWB spectrum do not consider full spectrum. In
most scenarios, only circular inspiral of SMBHBs is seriously taken into
account. In this letter, the inspiral-merger-cutoff phases are jointed
smoothly following methods proposed by (Poisson, 1993; Finn and Chernoff,
1993; Zhu  _et al._ , 2011; Chen  _et al._ , 2017). We express the complete
description of energy spectrum $dE_{\rm GW}/{df_{r}}$ as follow

$\displaystyle\frac{dE_{\rm GW}}{df_{r}}(f_{r}<\nu_{1})$ $\displaystyle=$
$\displaystyle\frac{\pi c^{2}f}{4G}h^{2}_{\rm c,fit}\left(f\frac{f_{\rm
p,0}}{f_{\rm p,t}}\right)\left(\frac{f_{\rm p,0}}{f_{\rm
p,t}}\right)^{-\frac{4}{3}}$ (10)
$\displaystyle\times\left(\frac{\mathcal{M}}{\mathcal{M}_{0}}\right)^{\frac{5}{3}}\left(\frac{1+z}{1+z_{0}}\right)^{-\frac{1}{3}}$
$\displaystyle\frac{dE_{\rm GW}}{df_{r}}(f_{r}\in[\nu_{1},\nu_{2}))$
$\displaystyle=$
$\displaystyle\frac{(G\pi)^{2/3}\mathcal{M}^{5/3}}{3}\omega_{1}f_{r}^{2/3}$
(11) $\displaystyle\frac{dE_{\rm GW}}{df_{r}}(f_{r}\in[\nu_{2},\nu_{3}))$
$\displaystyle=$
$\displaystyle\frac{(G\pi)^{2/3}\mathcal{M}^{5/3}}{3}\omega_{2}\left[\frac{f_{r}}{1+(\frac{f_{r}-\nu_{2}}{\sigma/2})^{2}}\right]^{2}$

where $\omega_{1}=\nu_{1}^{-1}$, $\omega_{2}=\nu_{1}^{-1}\nu_{2}^{-4/3}$. The
set of parameters $(\nu_{1},\nu_{2},\sigma,\nu_{3})$ can be determined by two
physical parameters: the total mass $M_{\rm total}=M(1+q)$ and the symmetric
mass ratio $\eta=(qM^{2}/M_{\rm total}^{2})$ in terms of
$(a\eta^{2}+b\eta+c)/(\pi GM_{\rm total}/c^{3})$, with coefficients $a,b,c$
given by Table.1. The ratio ${f_{\rm p,0}}/{f_{\rm p,t}}$ for shift is
illustrated as (Chen  _et al._ , 2017)

$\displaystyle\frac{f_{\rm p,0}}{f_{\rm
p,t}}=\frac{f_{0}}{f_{t}}\left[\left(\frac{e_{\rm
ref}}{e_{0}}\right)^{\frac{12}{19}}\frac{1-e_{0}^{2}}{1-e_{\rm
ref}^{2}}\left(\frac{304+{121e_{\rm
ref}^{2}}}{304+{121e_{0}^{2}}}\right)^{\frac{870}{2299}}\right]^{\frac{3}{2}}\
,$

here $f_{0}=10^{-10}{\rm Hz}$ and $e_{\rm ref}=0.9$ are the reference
frequency and ecctricity, respectively. $e_{0}$ is the initial eccentricity we
choose. Once $e_{0}$ is selected, turn-over frequency $f_{t}$ is obtained by

$\displaystyle f_{t}=0.356{\rm
nHz}\left(\frac{\rho_{i,100}}{F(e)\sigma_{200}}\zeta_{0}\right)^{\frac{3}{10}}\mathcal{M}_{9}^{-\frac{2}{5}}\
,$ (14)

where

$\displaystyle F(e)=\frac{1+(73/24)e^{2}+(37/96)e^{4}}{(1-e^{2})^{7/2}}\ ,$
(15)

$\mathcal{M}_{9}=\mathcal{M}/(10^{9}M_{\odot})$ is the rescaled chirp mass,
$\rho_{i,100}=\rho_{i}/(100M_{\odot}\ {\rm pc^{-3}})$ is the velocity
diepersion of stars in the galaxy and $\zeta_{0}$ describes the density of the
stellar environment. $h^{2}_{\rm c,fit}$ is the analytical fitting function of
reference spectrum and takes the form as (Chen  _et al._ , 2017)

$h_{\rm
c,fit}=a_{0}\bar{f}^{a_{1}}e^{-a_{2}\bar{f}}+b_{0}\bar{f}^{b_{1}}e^{-b_{2}\bar{f}}+c_{0}\bar{f}^{-c_{1}}e^{-c_{2}/\bar{f}}\
,$ (16)

where $a_{i},b_{i},c_{i}$, listed in Table.2, are constants determined by
fitting and $\bar{f}=f/(10^{-8}){\rm Hz}$.

Table 2: Constant parameters in the analytic functions of Eq. (16). These parameters are obtained from fitting functions to reference spectrum (Chen  _et al._ , 2017).  parameter | $i=0$ |  $i=1$ |  $i=2$
---|---|---|---
$a_{i}$ | $7.27\times 10^{-14}$ | 0.254 | 0.807
$b_{i}$ | $1.853\times 10^{-12}$ | 1.77 | 3.7
$c_{i}$ | $1.12\times 10^{-13}$ | 0.676 | 0.6

Data and results. The NANOGrav collaboration has performed an analysis on the
15-yr data set by employing a free spectrum that enables independent
variations in the amplitude of the GW spectrum across different frequency
bins. In the analyses, we use the posterior data from NANOGrav (Agazie  _et
al._ , 2023b; Afzal  _et al._ , 2023), and the PTMCMCSampler (Ellis and van
Haasteren, 2017) package to perform the Markov Chain Monte Carlo sampling. All
of the parameters and their prior distributions are listed in Table 3. We note
that these constrained prior distributions are based on those presented in
(Chen  _et al._ , 2019a) and are derived from observational and theoretical
works on the measurement of the GSMF, galaxy pair fraction, merger timescale
and SMBH-host galaxy scaling relations.

The resulting posterior distribution is illustrated in Fig. 2. Our analysis
reveals that the detected stochastic signal provides some new insights into
the differential pair fraction $\mathcal{F}(M_{G},z,q_{G})$, merger timescale
$\tau({M_{G}},z,q_{G})$, galaxy-SMBH mass scaling $M_{G}-M_{\rm{bulge}}-M$
relation, and the initial states of the SMBHBs when GW emission takes over,
such as the eccentricity and the transition frequency $f_{t}$, compared to
other astrophysical observations. Specifically, the preference for a higher
value of the parameter $f_{0}$ and the positive-skewed parameters $\alpha_{f}$
and $\beta_{f}$ suggest larger differential pair fractions in more massive
galaxies, while the preference for a lower value of the parameter $\tau_{0}$
and the negative-skewed parameters $\alpha_{\tau}$ and $\beta_{\tau}$ indicate
shorter merger timescales in more massive galaxies, and the preference for a
higher value of $M_{*}$ corresponds to a higher normalization between the
galaxy bulge mass $M_{\text{bulge}}$ and SMBH mass $M$. The above parameters
entirely contribute to the observed relatively high amplitude of the SGWB
spectrum ($\Omega_{\rm{GW}}=0.93^{+1.17}_{-0.41}\times 10^{-8}$ at
$\text{yr}^{-1}$). Note that the posterior distributions of these parameters
are very similar to those reported in (Middleton  _et al._ , 2021), where the
NANOGrav 12.5-yr data set was used. This is because the spectrum amplitudes in
both data sets are statistically consistent. On the other hand, the parameters
$e_{0}$ and $\zeta_{0}$, which determine the shape of the SGWB spectrum,
display sharp contrasts in the posteriors obtained from the two data sets. For
the NANOGrav 15-yr data set, the distribution of $e_{0}$ indicates that SMBHs
exhibit a large initial eccentricity when transitioning into the GW-emission
dominated process, while the larger value of the parameter $\zeta_{0}$ implies
that massive galaxies have, on average, higher densities than what is
suggested by a standard Dehnen profile (Dehnen, 1993).

Table 3: List of the 18 parameters and their prior in the astro-informed formation model.  parameter |  description |  prior
---|---|---
GSMF
$\Phi_{0}$ |  GSMF norm |  $-2.77_{-0.29}^{+0.27}$
$\Phi_{I}$ |  GSMF norm redshift evolution |  $-0.27_{-0.21}^{+0.23}$
$\log_{10}M_{G0}$ |  GSMF scaling mass |  $11.24_{-0.17}^{+0.20}$
$\alpha_{0}$ |  GSMF mass slope |  $-1.24_{-0.16}^{+0.16}$
$\alpha_{I}$ |  GSMF mass slope redshift evolution |  $-0.03_{-0.14}^{+0.16}$
Pair function
$f_{0}^{\prime}$ |  pair fraction norm |  Uniform $[0.01,0.05]$
$\alpha_{f}$ |  pair fraction mass slope |  Uniform $[-0.5,0.5]$
$\beta_{f}$ |  pair fraction redshift slope |  Uniform $[0,2]$
$\gamma_{f}$ |  pair fraction mass ratio slope |  Uniform $[-0.2,0.2]$
Galaxy merger timescale
$\tau_{0}$ |  merger time norm |  Uniform $[0.1,10]$
$\alpha_{\tau}$ |  merger time mass slope |  Uniform $[-0.5,0.5]$
$\beta_{\tau}$ |  merger time redshift slope |  Uniform $[-3,1]$
$\gamma_{\tau}$ |  merger time mass ratio slope |  Uniform $[-0.2,0.2]$
$M_{\rm bulge}$-$M_{\rm BH}$ relation
$\log_{10}M_{*}$ |  $M_{\rm bulge}$-$M_{\rm BH}$ relation norm |  $8.17_{-0.32}^{+0.35}$
$\alpha_{*}$ |  $M_{\rm bulge}$-$M_{\rm BH}$ relation slope |  $1.01_{-0.10}^{+0.08}$
$\epsilon$ |  $M_{\rm bulge}$-$M_{\rm BH}$ relation scatter |  Uniform $[0.2,0.5]$
Stellar condition
$e_{0}$ |  initial eccentricity |  Uniform $[0.01,0.99]$
$\log_{10}\zeta_{0}$ |  stellar density factor |  Uniform $[-2,2]$
Figure 2:  One- and two-dimensional prior distribution (in the orange color)
and marginalised posterior distribution (in the blue color) for the 18
parameters in the astro-informed model. We show the $1\sigma$ and $2\sigma$
contours in the two-dimensional plots.

Implication for GWs in the space-borne detector frequency band. In the PTA
frequency band, SMBHBs are in the inspiral phase, and their radiation power
can be calculated using Eq. (10). After a prolonged period of mutual inspiral,
these binaries gradually transition to more circular orbits and enter the
merge and ringdown phase, characterized by GW radiation described by Eq. (11)
and Eq. (Implications for the Supermassive Black Hole Binaries from the
NANOGrav 15-year Data Set). Some of these black holes undergo merger and final
coalescence at higher frequencies, entering the space-borne GW detector
frequency band. Now we can deduce the properties of the supermassive black
hole binary population from the PTA results, and further combine Eq. (1) with
Eqs.(10-Implications for the Supermassive Black Hole Binaries from the
NANOGrav 15-year Data Set) to obtain the complete SGWB spectrum generated by
MBHBs spanning both the PTA and LISA/Taiji/TianQin frequency bands, as
depicted in Fig. 1.

We need to emphasize that the general consensus for the GW detection of the
cosmic history of SMBHBs is that PTAs primarily detect the SGWB from the
ensemble of the SMBHB population, while LISA/Taiji/TianQin is expected to
detect the final coalescence stage of individual systems. However, during the
initial stages of the detector in operation, we cannot directly resolve
individual sources, and it is reasonable to consider these sources as
constituting a SGWB. In fact, as depicted in Fig. 1, the spectrum of the SGWB
is sufficiently strong that it is very likely to be detected very soon once
the detectors are in operation.

The relation between merger rate and SMBHB chirp mass is depicted in Fig. 3.
The merger rate of SMBHBs with chirp masses in the range
$10^{5}-10^{11}M_{\odot}$ is
$0.049\,\mathrm{yr}^{-1}\lesssim\mathcal{R}\lesssim 30.58\,\mathrm{yr}^{-1}$ .

Figure 3:  Posterior distributions of the SMBHB merger rate per year per unit
logarithmic chirp mass over the range
$\mathcal{M}\in[10^{5},10^{11}]M_{\odot}$ where the redshift is integrated
over $[0,5]$. The solid line is the median and the shaded region represents
the $90\%$ credible region.

Summary. In this letter, we use the NANOGrav 15-year data set to constrain the
SGWB from SMBHBs, implying that the SMBHBs tend to have a large initial
eccentricity in the transition phase between interaction domination and GW
domination. The SGWB spectrum from SMBHBs is extrapolated from the PTA
frequency band to the space-borne GW detector frequency band. Our results
indicate that such a SGWB from SMBHBs should be also detected by
LISA/Taiji/TianQin in the near future.

Acknowledgements. We acknowledge the use of HPC Cluster of ITP-CAS. QGH is
supported by the grants from NSFC (Grant No. 12250010, 11975019, 11991052,
12047503), Key Research Program of Frontier Sciences, CAS, Grant No. ZDBS-
LY-7009, CAS Project for Young Scientists in Basic Research YSBR-006, the Key
Research Program of the Chinese Academy of Sciences (Grant No. XDPB15). ZCC is
supported by the National Natural Science Foundation of China (Grant No.
12247176 and No. 12247112) and the China Postdoctoral Science Foundation
Fellowship No. 2022M710429.

## References

  * Akiyama  _et al._ (2019) Kazunori Akiyama _et al._ (Event Horizon Telescope), “First M87 Event Horizon Telescope Results. I. The Shadow of the Supermassive Black Hole,” Astrophys. J. Lett. 875, L1 (2019), arXiv:1906.11238 [astro-ph.GA] .
  * Kormendy and Ho (2013a) John Kormendy and Luis C. Ho, “Coevolution (Or Not) of Supermassive Black Holes and Host Galaxies,” Ann. Rev. Astron. Astrophys. 51, 511–653 (2013a), arXiv:1304.7762 [astro-ph.CO] .
  * Begelman  _et al._ (1980) M. C. Begelman, R. D. Blandford, and M. J. Rees, “Massive black hole binaries in active galactic nuclei,” Nature 287, 307–309 (1980).
  * Dotti  _et al._ (2012a) M. Dotti, A. Sesana, and R. Decarli, “Massive black hole binaries: dynamical evolution and observational signatures,” Adv. Astron. 2012, 940568 (2012a), arXiv:1111.0664 [astro-ph.CO] .
  * Chen  _et al._ (2019a) Siyuan Chen, Alberto Sesana, and Christopher J. Conselice, “Constraining astrophysical observables of Galaxy and Supermassive Black Hole Binary Mergers using Pulsar Timing Arrays,” Mon. Not. Roy. Astron. Soc. 488, 401–418 (2019a), arXiv:1810.04184 [astro-ph.GA] .
  * Kocsis and Sesana (2011) Bence Kocsis and Alberto Sesana, “Gas driven massive black hole binaries: signatures in the nHz gravitational wave background,” Mon. Not. Roy. Astron. Soc. 411, 1467 (2011), arXiv:1002.0584 [astro-ph.CO] .
  * Sesana (2013) A. Sesana, “Insights into the astrophysics of supermassive black hole binaries from pulsar timing observations,” Class. Quant. Grav. 30, 224014 (2013), arXiv:1307.2600 [astro-ph.CO] .
  * Ravi  _et al._ (2014) V. Ravi, J. S. B. Wyithe, R. M. Shannon, G. Hobbs, and R. N. Manchester, “Binary supermassive black hole environments diminish the gravitational wave signal in the pulsar timing band,” Mon. Not. Roy. Astron. Soc. 442, 56–68 (2014), arXiv:1404.5183 [astro-ph.CO] .
  * Kelley  _et al._ (2017) Luke Zoltan Kelley, Laura Blecha, and Lars Hernquist, “Massive Black Hole Binary Mergers in Dynamical Galactic Environments,” Mon. Not. Roy. Astron. Soc. 464, 3131–3157 (2017), arXiv:1606.01900 [astro-ph.HE] .
  * Chen  _et al._ (2017) Siyuan Chen, Alberto Sesana, and Walter Del Pozzo, “Efficient computation of the gravitational wave spectrum emitted by eccentric massive black hole binaries in stellar environments,” Mon. Not. Roy. Astron. Soc. 470, 1738–1749 (2017), arXiv:1612.00455 [astro-ph.CO] .
  * Rodig  _et al._ (2011) Constanze Rodig, Massimo Dotti, Alberto Sesana, Jorge Cuadra, and Monica Colpi, “Limiting eccentricity of sub-parsec massive black hole binaries surrounded by self-gravitating gas discs,” Mon. Not. Roy. Astron. Soc. 415, 3033 (2011), arXiv:1104.3868 [astro-ph.CO] .
  * Quinlan (1996) Gerald D. Quinlan, “The dynamical evolution of massive black hole binaries - I. hardening in a fixed stellar background,” New Astron. 1, 35–56 (1996), arXiv:astro-ph/9601092 .
  * Armitage and Natarajan (2005) Philip J. Armitage and Priyamvada Natarajan, “Eccentricity of supermassive black hole binaries coalescing from gas rich mergers,” Astrophys. J. 634, 921–928 (2005), arXiv:astro-ph/0508493 .
  * Bell  _et al._ (2006) Eric F. Bell, Stefanie Phleps, Rachel S. Somerville, Christian Wolf, Andrea Borch, and Klaus Meisenheimer, “The merger rate of massive galaxies,” Astrophys. J. 652, 270–276 (2006), arXiv:astro-ph/0602038 .
  * de Ravel  _et al._ (2009) L. de Ravel _et al._ , “The VIMOS VLT Deep Survey :Evolution of the major merger rate since z~1 from spectroscopicaly confirmed galaxy pairs,” Astron. Astrophys. 498, 379 (2009), arXiv:0807.2578 [astro-ph] .
  * Sesana  _et al._ (2008) Alberto Sesana, Alberto Vecchio, and Carlo Nicola Colacino, “The stochastic gravitational-wave background from massive black hole binary systems: implications for observations with Pulsar Timing Arrays,” Mon. Not. Roy. Astron. Soc. 390, 192 (2008), arXiv:0804.4476 [astro-ph] .
  * Foster and Backer (1990) R. S. Foster and D. C. Backer, “Constructing a Pulsar Timing Array,” Astrophys. J. 361, 300 (1990).
  * Hellings and Downs (1983) R. w. Hellings and G. s. Downs, “UPPER LIMITS ON THE ISOTROPIC GRAVITATIONAL RADIATION BACKGROUND FROM PULSAR TIMING ANALYSIS,” Astrophys. J. Lett. 265, L39–L42 (1983).
  * Agazie  _et al._ (2023a) Gabriella Agazie _et al._ (NANOGrav), “The NANOGrav 15-year Data Set: Observations and Timing of 68 Millisecond Pulsars,” Astrophys. J. Lett. 951 (2023a), 10.3847/2041-8213/acda9a, arXiv:2306.16217 [astro-ph.HE] .
  * Agazie  _et al._ (2023b) Gabriella Agazie _et al._ (NANOGrav), “The NANOGrav 15-year Data Set: Evidence for a Gravitational-Wave Background,”  (2023b), 10.3847/2041-8213/acdac6, arXiv:2306.16213 [astro-ph.HE] .
  * Afzal  _et al._ (2023) Adeela Afzal _et al._ (NANOGrav), “The NANOGrav 15-year Data Set: Search for Signals from New Physics,” Astrophys. J. Lett. 951 (2023), 10.3847/2041-8213/acdc91, arXiv:2306.16219 [astro-ph.HE] .
  * Antoniadis  _et al._ (2023a) J. Antoniadis _et al._ , “The second data release from the European Pulsar Timing Array I. The dataset and timing analysis,”  (2023a), 10.1051/0004-6361/202346841, arXiv:2306.16224 [astro-ph.HE] .
  * Antoniadis  _et al._ (2023b) J. Antoniadis _et al._ , “The second data release from the European Pulsar Timing Array III. Search for gravitational wave signals,” (2023b), arXiv:2306.16214 [astro-ph.HE] .
  * Zic  _et al._ (2023) Andrew Zic _et al._ , “The Parkes Pulsar Timing Array Third Data Release,” (2023), arXiv:2306.16230 [astro-ph.HE] .
  * Reardon  _et al._ (2023) Daniel J. Reardon _et al._ , “Search for an isotropic gravitational-wave background with the Parkes Pulsar Timing Array,” Astrophys. J. Lett. 951 (2023), 10.3847/2041-8213/acdd02, arXiv:2306.16215 [astro-ph.HE] .
  * Nan  _et al._ (2011) Rendong Nan, Di Li, Chengjin Jin, Qiming Wang, Lichun Zhu, Wenbai Zhu, Haiyan Zhang, Youling Yue, and Lei Qian, “The Five-Hundred-Meter Aperture Spherical Radio Telescope (FAST) Project,” Int. J. Mod. Phys. D 20, 989–1024 (2011), arXiv:1105.3794 [astro-ph.IM] .
  * Xu  _et al._ (2023) Heng Xu _et al._ , “Searching for the nano-Hertz stochastic gravitational wave background with the Chinese Pulsar Timing Array Data Release I,”  (2023), 10.1088/1674-4527/acdfa5, arXiv:2306.16216 [astro-ph.HE] .
  * Agazie  _et al._ (2023c) Gabriella Agazie _et al._ (NANOGrav), “The NANOGrav 15-year Data Set: Constraints on Supermassive Black Hole Binaries from the Gravitational Wave Background,” (2023c), arXiv:2306.16220 [astro-ph.HE] .
  * Agazie  _et al._ (2023d) Gabriella Agazie _et al._ (NANOGrav), “The NANOGrav 15-year Data Set: Bayesian Limits on Gravitational Waves from Individual Supermassive Black Hole Binaries,” (2023d), arXiv:2306.16222 [astro-ph.HE] .
  * Antoniadis  _et al._ (2023c) J. Antoniadis _et al._ , “The second data release from the European Pulsar Timing Array: V. Implications for massive black holes, dark matter and the early Universe,” (2023c), arXiv:2306.16227 [astro-ph.CO] .
  * Amaro-Seoane  _et al._ (2017) Pau Amaro-Seoane _et al._ (LISA), “Laser Interferometer Space Antenna,” (2017), arXiv:1702.00786 [astro-ph.IM] .
  * Luo  _et al._ (2020) Ziren Luo, Yan Wang, Yueliang Wu, Wenrui Hu, and Gang Jin, “The Taiji program: A concise overview,” Progress of Theoretical and Experimental Physics 2021 (2020), 10.1093/ptep/ptaa083.
  * Luo  _et al._ (2016) Jun Luo _et al._ (TianQin), “TianQin: a space-borne gravitational wave detector,” Class. Quant. Grav. 33, 035010 (2016), arXiv:1512.02076 [astro-ph.IM] .
  * Steinle  _et al._ (2023) Nathan Steinle, Hannah Middleton, Christopher J. Moore, Siyuan Chen, Antoine Klein, Geraint Pratten, Riccardo Buscicchio, Eliot Finch, and Alberto Vecchio, “Implications of pulsar timing array observations for LISA detections of massive black hole binaries,” (2023), arXiv:2305.05955 [astro-ph.HE] .
  * Dotti  _et al._ (2012b) M. Dotti, A. Sesana, and R. Decarli, “Massive black hole binaries: dynamical evolution and observational signatures,” Adv. Astron. 2012, 940568 (2012b), arXiv:1111.0664 [astro-ph.CO] .
  * Phinney (2001) E. S. Phinney, “A Practical theorem on gravitational wave backgrounds,” (2001), arXiv:astro-ph/0108028 .
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Sesana  _et al._ (2009) 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, 2255 (2009), arXiv:0809.3412 [astro-ph] .
  * Kormendy and Ho (2013b) John Kormendy and Luis C. Ho, “Coevolution (Or Not) of Supermassive Black Holes and Host Galaxies,” Ann. Rev. Astron. Astrophys. 51, 511–653 (2013b), arXiv:1304.7762 [astro-ph.CO] .
  * Mortlock  _et al._ (2015) Alice Mortlock _et al._ , “Deconstructing the Galaxy Stellar Mass Function with UKIDSS and CANDELS: the Impact of Colour, Structure and Environment,” Mon. Not. Roy. Astron. Soc. 447, 2–24 (2015), arXiv:1411.3339 [astro-ph.GA] .
  * Conselice  _et al._ (2016) Christopher J. Conselice, Aaron Wilkinson, Kenneth Duncan, and Alice Mortlock, “THE EVOLUTION OF GALAXY NUMBER DENSITY ATz$<$ 8 AND ITS IMPLICATIONS,” Astrophys. J. 830, 83 (2016).
  * Snyder  _et al._ (2017) Gregory F. Snyder, Jennifer M. Lotz, Vicente Rodriguez-Gomez, Renato da Silva Guimarães, Paul Torrey, and Lars Hernquist, “Massive close pairs measure rapid galaxy assembly in mergers at high redshift,” Monthly Notices of the Royal Astronomical Society 468, 207–216 (2017).
  * Chen  _et al._ (2019b) Zu-Cheng Chen, Fan Huang, and Qing-Guo Huang, “Stochastic Gravitational-wave Background from Binary Black Holes and Binary Neutron Stars and Implications for LISA,” Astrophys. J. 871, 97 (2019b), arXiv:1809.10360 [gr-qc] .
  * Poisson (1993) Eric Poisson, “Gravitational radiation from a particle in circular orbit around a black hole. 1: Analytical results for the nonrotating case,” Phys. Rev. D 47, 1497–1510 (1993).
  * Finn and Chernoff (1993) Lee Samuel Finn and David F. Chernoff, “Observing binary inspiral in gravitational radiation: One interferometer,” Phys. Rev. D 47, 2198–2219 (1993), arXiv:gr-qc/9301003 .
  * Zhu  _et al._ (2011) Xing-Jiang Zhu, E. Howell, T. Regimbau, D. Blair, and Zong-Hong Zhu, “Stochastic Gravitational Wave Background from Coalescing Binary Black Holes,” Astrophys. J. 739, 86 (2011), arXiv:1104.3565 [gr-qc] .
  * Ellis and van Haasteren (2017) Justin Ellis and Rutger van Haasteren, “jellis18/ptmcmcsampler: Official release,”  (2017).
  * Middleton  _et al._ (2021) H. Middleton, A. Sesana, S. Chen, A. Vecchio, W. Del Pozzo, and P. A. Rosado, “Massive black hole binary systems and the NANOGrav 12.5 yr results,” Mon. Not. Roy. Astron. Soc. 502, L99–L103 (2021), arXiv:2011.01246 [astro-ph.HE] .
  * Dehnen (1993) W. Dehnen, “A Family of Potential-Density Pairs for Spherical Galaxies and Bulges,” Mon. Not. Roy. Astron. Soc. 265, 250 (1993).