Gravitational Waves from Supermassive Black Hole Binaries in Ultraluminous Infrared Galaxies

Most massive galaxies in the local universe host supermassive black holes (SMBHs). One fundamental goal in astrophysics is to understand the origins of these SMBHs and their host galaxies. Galaxy mergers play an important role in their evolution on a cosmological timescale (e.g., Kormendy & Ho 2013), in the assembly of massive galaxies, and in fueling gas to nuclear SMBHs. A natural outcome of galaxy mergers containing SMBHs is the formation of binary SMBHs. If the binary SMBHs coalescence within a Hubble time, a significant fraction of their rest-mass energy is emitted as gravitational waves (GWs). The GW emission is a good tool to probe the cosmological evolution of SMBHs in the framework of hierarchical structure formation.

Pulsar timing array (PTA) experiments enable us to directly address the GW emission in the nHz−μHz band. There are three ongoing PTA experiments: the European Pulsar Timing Array (EPTA), the Australian Parkes Pulsar Timing Array (PPTA), and the North American Nanohertz Observatory for Gravitational Waves (NANOGrav). Their measurements have recently provided upper limits on the strength of the GW background (GWB) from binary SMBHs (Lentati et al. 2015; Shannon et al. 2015; Arzoumanian et al. 2016).

The theoretical GWB signal in the PTA band has been predicted by using semi-analytical calculations (e.g., Jaffe & Backer 2003) and cosmological simulations (e.g., Sesana et al. 2009). Those models incorporated various combinations of physical processes that affect the evolution of binary SMBHs in merging galaxies (Begelman et al. 1980): eccentric binary evolution (Enoki & Nagashima 2007), viscous drag from a circumbinary gaseous disk (Kocsis & Sesana 2011), dynamical friction (DF; McWilliams et al. 2014; Kulier et al. 2015), and multibody BH interactions (Bonettiet al. 2018; Ryu et al. 2018). Combined with cosmological hydrodynamical simulations, Kelley et al. (2017) investigated the impact of various environmental processes on the binary evolution. In spite of great efforts in theory, there are still many uncertainties for model parameters and observed empirical relations (Sesana 2013).

In this Letter, we address this issue with a different approach following observational results. We consider ultraluminous infrared galaxies (ULIRGs), which are one of the best tracers of merging galaxies containing SMBHs, as sources of a GWB. Recent observations by Herschel and the Wide-field Infrared Survey Explorer (WISE) have provided a large sample of bright IR galaxies at 0 < z < 4 and allow us to explore their infrared spectral energy distributions (SEDs; Delvecchio et al. 2014). Adopting the luminosity function (LF) of ULIRGs, we study the development of a GWB due to coalescing binary SMBHs driven by merging ULIRGs. Intriguingly, we find a tension with the most stringent PTA upper limits that constrain the fraction of binary SMBHs that coalesce by the present time, depending on the Eddington ratio of BHs in ULIRGs, the typical lifetime of ULIRGs, and the fraction of ULIRGs associated with SMBH binaries.

We consider here ULIRGs with total infrared luminosities of ${L}_{\mathrm{IR},\mathrm{tot}}\gtrsim {10}^{12}\,{\text{}}{L}_{\odot }$ as good tracers of merging galaxies that host at least one accreting SMBH in each system. At low redshifts of z < 1, infrared observations have revealed that the morphologies of ULIRGs are almost exclusively caused by mergers (e.g., Surace et al. 1998) and that the number fraction containing AGNs increases with ${L}_{\mathrm{IR},\mathrm{tot}}$ and reaches almost unity in ULIRGs (e.g., Veilleux et al. 2002; Ichikawa et al. 2014). At higher redshifts of z > 1, the merger fraction is still uncertain since active star formation in gas-rich galaxies alone could produce a similar level of infrared luminosities observed as ULIRGs (e.g., Kartaltepe et al. 2012).

Recent observations have discovered extremely bright ULIRGs with ${L}_{\mathrm{IR},\mathrm{tot}}\simeq {10}^{13}\mbox{--}{10}^{14}\,{\text{}}{L}_{\odot }$, the so-called hyperluminous infrared galaxies (HyLIRGs; Assef et al. 2015; Tsai et al. 2015). This population shows relatively clear signatures of galaxy mergers even at higher redshifts. The merger fraction has been estimated by quantifying their morphologies as ${f}_{{\rm{m}},\mathrm{gal}}\sim 62\pm 14 \%$ for hot, dust-obscured galaxies at z ∼ 3 (Fan et al. 2016) and ∼75% for HyLIRGs at 1.8 < z < 2.5 (Farrah et al. 2017). Moreover, the merger fraction tends to increase with bolometric luminosity and becomes ∼80% at the brightest end of ${L}_{\mathrm{bol}}\sim (1\mbox{--}5)\times {10}^{14}\,{\text{}}{L}_{\odot }$ (Glikman et al. 2015). Those measurements give a lower limit of the fraction because such merger signatures would smooth out after several dynamical timescales.

The enormous power of U/HyLIRGs is produced by deeply buried and rapidly accreting SMBHs. Most of the AGN radiation is absorbed by surrounding dust and is re-emitted at mid-infrared wavelengths. By decomposing the hot dust emission from their SEDs, the AGN contribution to the total (IR) luminosity increases with ${L}_{\mathrm{IR},\mathrm{tot}}$ (e.g., Murphy et al. 2011). Namely, the luminosity ratio is estimated as ${L}_{\mathrm{IR},\mathrm{AGN}}/{L}_{\mathrm{IR},\mathrm{tot}}\simeq 0.2\mbox{--}0.3$ for ULIRGs (e.g., Ichikawa et al. 2014) and reaches ≃0.7–1.0 for HyLIRGs (e.g., Jones et al. 2014; Farrah et al. 2017).

Using the SED decomposition technique for a sample of Herschel-selected galaxies within the Great Observatories Origins Deep Survey (GOODS)-South and the Cosmic Evolution Survey (COSMOS) fields, Delvecchio et al. (2014) have reconstructed the AGN bolometric LF at 0.1 < z < 3.8. The bolometric correction factor is estimated by solving the radiative transfer equation for a smooth dusty structure irradiated by the AGN accretion disk, instead of adopting the bolometric correction shown in Hopkins et al. (2007). The AGN bolometric LF for ULIRGs is well fit by

$$\begin{eqnarray}&&\displaystyle \frac{d\phi (L,z)}{d\mathrm{log}L}={\phi }_{\star }{\left(\displaystyle \frac{L}{{L}_{\star }}\right)}^{\beta }\exp \left[-\displaystyle \frac{\{\mathrm{log}(1+\tfrac{L}{{L}_{\star }})\}{}^{2}}{2{\sigma }^{2}}-\displaystyle \frac{L}{{L}_{{\rm{c}}}}\right],\end{eqnarray} \tag{ 1 }$$

where L is the AGN bolometric luminosity, and the values of fitting parameters (ϕ_⋆, L_⋆, β, and σ) are listed in Table 1 of Delvecchio et al. (2014). Here, we set an exponential cutoff above a critical luminosity ${L}_{{\rm{c}}}(\equiv {10}^{14}\,{\text{}}{L}_{\odot })$ because there is no detection for bright AGNs with $L\gt {10}^{14}\,{\text{}}{L}_{\odot }$ due to the lack of such bright ones or due to the limited observation volume. In addition, we set a minimum value of the AGN luminosity driven by galaxy major mergers to ${L}_{\min }={10}^{12}\,{\text{}}{L}_{\odot }$. Considering $L/{L}_{\mathrm{IR},\mathrm{AGN}}\simeq 3$ (Delvecchio et al. 2014) and ${L}_{\mathrm{IR},\mathrm{AGN}}/{L}_{\mathrm{IR},\mathrm{tot}}\simeq 0.2\mbox{--}0.3$ (Ichikawa et al. 2014), ${L}_{\min }={10}^{12}\,{\text{}}{L}_{\odot }$ corresponds to ${L}_{\mathrm{IR},\mathrm{tot}}\,\simeq (1\mbox{--}2)\times {10}^{12}\,{\text{}}{L}_{\odot }$. Since the merger fraction of IR galaxies with ${L}_{\mathrm{IR},\mathrm{tot}}\lt {10}^{12}\,{\text{}}{L}_{\odot }$ may be significantly below unity at high redshift, we do not consider such galaxies as sources of a GWB. Figure 1 (left panel) shows the luminosity function of ULIRGs for different redshifts.

Zoom In Zoom Out Reset image size

In order to convert the LF to the BH mass function (MF), it is necessary to obtain the Eddington ratio (${\lambda }_{\mathrm{Edd}}=L/{L}_{\mathrm{Edd}}$). Since the broad-line regions of ULIRGs are completely obscured in the optical, it is almost impossible to estimate their BH masses and thus the Eddington ratio using the optical spectra. However, a well-defined sample of quasars obtained from the Sloan Digital Sky Surveys (SDSS) catalog suggests that the typical Eddington ratio for those quasars is ≃0.3 (Kollmeier et al. 2006). On the other hand, the largest SMBHs have a maximum mass limit at ${M}_{\mathrm{BH},\max }\sim {\rm{a}}\,\mathrm{few}\times {10}^{10}\,{\text{}}{M}_{\odot }$, which is nearly independent of redshift (Netzer 2003; Wu et al. 2015). The radiation luminosity from quasars hosting the most massive BHs is estimated as $L\simeq {10}^{15}{\lambda }_{\mathrm{Edd}}\,{\text{}}{L}_{\odot }({M}_{\mathrm{BH},\max }/3\times {10}^{10}\,{\text{}}{M}_{\odot })$. Since this luminosity would be ≳L_c, we obtain ${\lambda }_{\mathrm{Edd}}\gtrsim 0.1{({M}_{\mathrm{BH},\max }/3\times {10}^{10}{\text{}}{M}_{\odot })}^{-1}$. Thus, we adopt λ_Edd = 0.3 as our fiducial value. Figure 1 (right panel) shows the BH MF in ULIRGs for different redshifts. This BH MF of ULIRGs is consistent with that of SDSS quasars (QSOs) obtained by Kelly & Shen (2013), where the BH masses are estimated by using the width of the broad emission lines and the AGN continuum luminosity.

Luminous QSOs and ULIRGs are much rarer than normal galaxies, which is expected since those luminous phases have a lifetime shorter than a Hubble time. The lifetime is one of the most fundamental quantities for estimating the intrinsic number density of those luminous objects. The QSO lifetime can be observationally constrained by several methods (Martini 2004 and references therein). Overall, the QSO lifetime lies in the range of $1\,\mathrm{Myr}\lesssim {t}_{\mathrm{life}}\lesssim 100\,\mathrm{Myr}$. Using galaxy merger simulations, Hopkins et al. (2006) demonstrated that the lifetime tends to decrease with increasing luminosity, namely, ${t}_{\mathrm{life}}\simeq 10\mbox{--}50\,\mathrm{Myr}$ for ${L}_{\mathrm{bol}}\gt {10}^{13}\,{\text{}}{L}_{\odot }$. This shorter lifetime is consistent with $1\,\mathrm{Myr}\lesssim {t}_{\mathrm{life}}\lesssim 20\,\mathrm{Myr}$ obtained from observations of extended Lyα emission near luminous QSOs at 2.5 ≲ z ≲ 2.9 with ultraviolet luminosities of ${L}_{\mathrm{UV}}\sim {10}^{14}\,{\text{}}{L}_{\odot }$ (Trainor & Steidel 2013). In this Letter, we adopt a conservative value of ${t}_{\mathrm{life}}\simeq 30\,\mathrm{Myr}$ as our fiducial case.

Following Phinney (2001), we estimate the GW energy density per logarithmic frequency interval³ as

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{GW}}(f)=\displaystyle \frac{1}{{\rho }_{{\rm{c}}}{c}^{2}}\int \displaystyle \frac{{d}^{2}N}{d{{ \mathcal M }}_{{\rm{c}}}{dz}}\,\displaystyle \frac{1}{1+z}\,{f}_{r}\displaystyle \frac{{{dE}}_{\mathrm{gw}}}{{{df}}_{r}}\,d{{ \mathcal M }}_{{\rm{c}}}{dz},\end{eqnarray} \tag{ 2 }$$

where ρ_c is the critical mass density of the universe at z = 0, f_r is the GW frequency in the source's cosmic rest frame, $f={f}_{r}/(1+z)$ is the observed GW frequency, z is the redshift when the GWs are produced, ${d}^{2}N/d{{ \mathcal M }}_{{\rm{c}}}{dz}$ is the comoving number density of GW events with chirp masses of [${{ \mathcal M }}_{{\rm{c}}},{{ \mathcal M }}_{{\rm{c}}}+d{{ \mathcal M }}_{{\rm{c}}}]$ that occurs at cosmic times corresponding to the redshift range between z and z + dz,

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}N}{d{{ \mathcal M }}_{{\rm{c}}}{dz}}\,d{{ \mathcal M }}_{{\rm{c}}}{dz}=\displaystyle \frac{{f}_{{\rm{m}},\mathrm{gal}}}{{f}_{\mathrm{duty}}}\,\displaystyle \frac{d\phi (L,z)}{d\mathrm{log}L}\,d\mathrm{log}L\,{dz},\end{eqnarray} \tag{ 3 }$$

where ${f}_{{\rm{m}},\mathrm{gal}}$ is the merger fraction of galaxies inferred from the morphologies of U/HyLIRGs (we set ${f}_{{\rm{m}},\mathrm{gal}}\simeq 1;$ see Section 2), and ${f}_{\mathrm{duty}}\equiv {t}_{\mathrm{life}}\cdot ({dz}/{{dt}}_{r})$ is the duty cycle of ULIRGs. Since a constant ${t}_{\mathrm{life}}\simeq 30\,\mathrm{Myr}$ is adopted, the duty cycle is independent of the luminosity.

The GW emission spectrum from a merging binary in the rest frame is given by

$$\begin{eqnarray}&&\displaystyle \frac{{{dE}}_{\mathrm{gw}}}{{{df}}_{r}}=\displaystyle \frac{{(\pi G)}^{2/3}{{ \mathcal M }}_{{\rm{c}}}^{5/3}}{3{f}_{r}^{1/3}},\end{eqnarray} \tag{ 4 }$$

where E_gw is the energy of the GW. The chirp mass is written as ${{ \mathcal M }}_{{\rm{c}}}^{5/3}\equiv {{qM}}_{\mathrm{BH}}^{5/3}/{(1+q)}^{1/3}$, where M_BH is the mass of the primary SMBH and q(<1) is the mass ratio of two SMBHs. We suppose that the primary BH follows the MF in Figure 1. Thus, we implicitly assume that the primary SMBH is located at the center of a ULIRG after the galaxy mergers and is responsible for the ULIRG activity in a lifetime of t_life, while the secondary BH is still located off center in a lower-density region. Since the secondary BH would decay its orbit via DF on a timescale of t_DF ∼ 100 Myr (e.g., Yu 2002), the binary formation would occur after the ULIRG phase (i.e., t_DF ≳ t_life). This assumption would be plausible because the number fraction of AGNs that are dual SMBHs is as small as ∼10% at z < 1 (Comerford et al. 2013).

As discussed in Section 2, most ULIRGs are triggered by major mergers of galaxies. We set here a minimum value of the BH mass ratio to q_min ≃ 0.1. The mass-ratio distribution of SMBHs, Φ(q) at q_min ≤ q ≤ 1, is uncertain. However, the chirp mass averaged over q is less uncertain, namely, $\langle {{ \mathcal M }}_{{\rm{c}}}^{5/3}\rangle /{M}_{\mathrm{BH}}^{5/3}\simeq 0.47$, 0.34, and 0.3 for ${\rm{\Phi }}(q)=\mathrm{const}.$, ${\rm{\Phi }}(q)\propto 1/q$, and ${\rm{\Phi }}(q)\propto \delta (q-1/3)$, respectively. We adopt the last one for a more conservative estimate.

In the course of a galaxy merger embedding two SMBHs, a variety of physical processes affect the binary evolution to the coalescence in a certain timescale (e.g., Begelman et al. 1980; Merritt 2013). Since the delay time between formation and coalescence of BH binaries is still uncertain, rather than attempting to model this delay, we assume a uniform delay time of t_delay. To consider the delay effect on BH mergers, we evaluate the LF and duty cycle in Equation (3) at $\tilde{z}\equiv z+{\rm{\Delta }}z$, where Δz is determined by solving ${t}_{\mathrm{delay}}={\int }_{z+{\rm{\Delta }}z}^{z}\tfrac{{{dt}}_{r}}{{dz}}{dz}$. The delay time tends to be shorter than a Hubble timescale for SMBH binaries with higher total masses (${M}_{\mathrm{BH}}\gt {10}^{8}\,{\text{}}{M}_{\odot }$) and mass ratios (q > 0.2) (Khan et al. 2016; Kelley et al. 2017). Within their model uncertainties, the delay time is estimated as ≃0.35–6.9 Gyr (see Table B1 in Kelley et al. 2017), which is significantly longer than the SMBH binary formation (t_DF ∼ 100 Myr) and the ULIRG's lifetime (t_life < 100 Myr).

By using Equations (2)–(4), the energy spectrum of the total GW emission due to merging SMBHs in ULIRGs is calculated as

$$\begin{eqnarray}{{\rm{\Omega }}}_{\mathrm{GW}}(f) & = & 1.53\times {10}^{-9}\,{{ \mathcal F }}_{\mathrm{GW}}\\ & & \times \,{\left(\displaystyle \frac{{\lambda }_{\mathrm{Edd}}}{0.3}\right)}^{-5/3}{\left(\displaystyle \frac{{t}_{\mathrm{life}}}{30\mathrm{Myr}}\right)}^{-1}{\left(\displaystyle \frac{f}{10\mathrm{nHz}}\right)}^{2/3},\end{eqnarray} \tag{ 5 }$$

and the characteristic strain is estimated as

$$\begin{eqnarray}{h}_{{\rm{c}}}(f) & = & 3.32\times {10}^{-15}\,{{ \mathcal F }}_{\mathrm{GW}}^{1/2}\\ & & \times \,{\left(\displaystyle \frac{{\lambda }_{\mathrm{Edd}}}{0.3}\right)}^{-5/6}{\left(\displaystyle \frac{{t}_{\mathrm{life}}}{30\mathrm{Myr}}\right)}^{-1/2}{\left(\displaystyle \frac{f}{10\mathrm{nHz}}\right)}^{-2/3},\end{eqnarray} \tag{ 6 }$$

where ${{ \mathcal F }}_{\mathrm{GW}}$ is the ratio of the GW amplitude with an assumed value of t_delay to that without delay. We also define the number fraction of SMBHs that coalesce within a Hubble time as ${{ \mathcal F }}_{\mathrm{coal}}$. As shown in Table 1, the coalescence fraction decreases monotonically with t_delay and drops sharply at t_delay ≳ 7 Gyr. On the other hand, the GW amplitude slightly increases with t_delay (i.e., ${{ \mathcal F }}_{\mathrm{GW}}\gt 1$) and decreases at t_delay ≳ 10 Gyr significantly (i.e., ${{ \mathcal F }}_{\mathrm{GW}}\lt 1$). This is because a short delay time barely reduces the number of merger events occurring in a Hubble time, but induces BH mergers at lower redshift. In Figure 2, we plot the total GW spectrum of interest (solid). The upper limits from the PTA are presented by triangle symbols at frequencies where the limit becomes the most stringent: EPA (Lentati et al. 2015), NANOGrav (Arzoumanian et al. 2016), and PPTA (Shannon et al. 2015). Without the delay effect, merging SMBHs in ULIRGs would overproduce a GWB. The GWB has also been overproduced based on the observed (McWilliams et al. 2014) and simulated (Kulier et al. 2015) abundance of massive galaxies, and based on the observed periodic quasar candidates (Sesana et al. 2018), assuming that these objects all host SMBH binary mergers. With these PTA constraints, we obtain an upper limit for ${{ \mathcal F }}_{\mathrm{GW}}$ as

$$\begin{eqnarray}&&{{ \mathcal F }}_{\mathrm{GW}}\lesssim 0.43{\left(\displaystyle \frac{{\lambda }_{\mathrm{Edd}}}{0.3}\right)}^{5/3}\left(\displaystyle \frac{{t}_{\mathrm{life}}}{30\,\mathrm{Myr}}\right),\end{eqnarray} \tag{ 7 }$$

which corresponds to t_delay ≳ 10.9 Gyr and ${{ \mathcal F }}_{\mathrm{coal}}\lesssim 0.16$. We also plot the sensitivity that is achievable with the Square Kilometre Array (SKA) with 50 pulsars for a T_obs = 10 years observation. If a GWB will not be detected even by such a planned detector, it would imply a strong constraint on ${{ \mathcal F }}_{\mathrm{coal}}\ll 1 \%$.

Zoom In Zoom Out Reset image size

Table 1.  The Relation between the Delay Time, Coalescence Fraction, and GWB Amplitude

tdelay(Gyr)	${{ \mathcal F }}_{\mathrm{coal}}$	${{ \mathcal F }}_{\mathrm{GW}}$	${h}_{c}\times {10}^{15}$
0	1.0	1.0	3.32
3	1.0	1.16	3.57
7	0.75	1.22	3.67
10	0.31	0.83	3.03
11	0.14	0.39	2.06
12	0.022	0.06	0.81

Note. (1) Delay time, (2) coalescence fraction of BHs, (3) effective coalescence fraction, (4) total GW amplitude, and (5) stochastic GWB amplitude at $f={10}^{-8}\,\mathrm{Hz}$. We adopt ${\lambda }_{\mathrm{Edd}}=0.3$ and ${t}_{\mathrm{life}}=30\,\mathrm{Myr}$.

Download table as:  ASCIITypeset image

In Figure 3, we present the evolution of the total stochastic GWB (black), and the GWB due to SMBHs with ${M}_{\mathrm{BH}}\geqslant {10}^{9}\,{\text{}}{M}_{\odot }$ (red) and $\lt {10}^{9}\,{\text{}}{M}_{\odot }$ (blue). Without the delay, half of the stochastic GWB energy is produced by merging SMBHs with ${M}_{\mathrm{BH}}\gtrsim {10}^{9}\,{\text{}}{M}_{\odot }$ at z > 1.5, while others are due to less massive ones at z < 1.5. This result reflects the shape and redshift-evolution of the LF of ULIRGs and the MF of SMBHs. In fact, a GWB from higher-mass BHs dominates at a higher redshift. The GWB amplitude decreases significantly at t_delay ≳ 10 Gyr because a larger fraction of SMBHs in ULIRGs do not merge within a Hubble time. Our results are qualitatively consistent with previous work (e.g., Sesana et al. 2008), concluding that a GBW in the PTA band is dominated by nearby and massive binary SMBHs (z < 2 and ${M}_{\mathrm{BH}}\gt {10}^{8}\,{\text{}}{M}_{\odot }$). However, it is worth emphasizing that coalescing binary SMBHs, even in a rare population of high-z ULIRGs associated with gas-rich major mergers, which are quite different from multiple dry mergers occurring at low redshift, can produce a GWB close to the present-day upper limit. Since other types of galaxies unlike ULIRGs may have additional SMBH mergers and may contribute to the GWBs (see, e.g., McWilliams et al. 2014), our result gives a lower limit on the total GWB in the PTA band.

Zoom In Zoom Out Reset image size

We briefly discuss possible biased estimates of the BH masses caused by a scatter of the Eddington ratio distribution (Shen et al. 2008). For the SDSS quasar sample, the Eddington ratio approximately follows a log-normal distribution with a dispersion of σ_E ≃ 0.25 dex for a fixed BH mass (see Figure 12 in Kelly & Shen 2013). We estimate the mass bias as ${\gamma }_{{\rm{L}}}{\sigma }_{{\rm{E}}}^{2}\mathrm{ln}10$ dex, assuming a power-law shape for the underlying true MF and a symmetric Gaussian scatter in $\mathrm{log}{M}_{\mathrm{BH}}$ around a mass-independent mean value, where γ_L(<0) is the slope of the AGN bolometric LF, and we find that the bias effect reduces the GWB amplitude by half. However, we note that the mean Eddington ratio for SMBHs of interest is significantly lower than our fiducial value (Kelly & Shen 2013). Therefore, the GWB amplitude would be rather enhanced (${{\rm{\Omega }}}_{\mathrm{GW}}\propto {\lambda }_{\mathrm{Edd}}^{-5/3}$), and the bias effect would be canceled out.

The growth of SMBHs in galactic nuclei can be constrained by observations of present-day BH remnants. Adopting the LF of ULIRGs and a radiative efficiency ${\epsilon }_{r}$, we can estimate the BH mass density accreted during the ULIRG phases, and we can compare it to that observed in the local universe (Soltan 1982; Yu & Tremaine 2002). The estimated BH mass density is given by

$$\begin{eqnarray}&&{\rho }_{\mathrm{BH}}\simeq 3.5\times {10}^{4}\,\left(\displaystyle \frac{1-{\epsilon }_{r}}{{\epsilon }_{r}}\right)\,{\text{}}{M}_{\odot }\,{\mathrm{cMpc}}^{-3}\end{eqnarray} \tag{ 8 }$$

(Delvecchio et al. 2014). In order not to exceed the value observed in the local universe, ${\rho }_{\mathrm{BH},\mathrm{obs}}\simeq {4.2}_{-1.0}^{+1.2}\times {10}^{5}\,{\text{}}{M}_{\odot }\,{\mathrm{cMpc}}^{-3}$ (e.g., Shankar et al. 2009), the radiative efficiency is required to be ${\epsilon }_{{\rm{r}}}\gtrsim {0.076}_{-0.018}^{+0.023}$. This efficiency is consistent with similar arguments for bright QSOs in the optical and X-ray bands (e.g., Yu & Tremaine 2002; Hopkins et al. 2007; Ueda et al. 2014).

The total present-day energy density in GW radiation is estimated from

$$\begin{eqnarray}&&{{ \mathcal E }}_{\mathrm{GW}}\equiv {\int }_{0}^{\infty }{\rho }_{{\rm{c}}}{c}^{2}{{\rm{\Omega }}}_{\mathrm{GW}}(f)\displaystyle \frac{{df}}{f},\end{eqnarray} \tag{ 9 }$$

and ${{ \mathcal E }}_{\mathrm{GW}}\simeq 4.3\times {10}^{-16}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-3}$ for ${t}_{\mathrm{delay}}=11\,\mathrm{Gyr}$, ${\lambda }_{\mathrm{Edd}}\,=0.3$, and ${t}_{\mathrm{life}}=30\,\mathrm{Myr}$. Note that ${{ \mathcal E }}_{\mathrm{GW}}\propto {\lambda }_{\mathrm{Edd}}^{-0.98}{t}_{\mathrm{life}}^{-1}$. As a result, the ratio of the total GW energy to the present-day SMBH rest-mass energy is estimated as

$$\begin{eqnarray}&&\displaystyle \frac{{{ \mathcal E }}_{\mathrm{GW}}}{{\rho }_{\mathrm{BH}}{c}^{2}}\lesssim \displaystyle \frac{0.20\,{\epsilon }_{r}}{1-{\epsilon }_{r}},\end{eqnarray} \tag{ 10 }$$

because ${t}_{\mathrm{delay}}\gtrsim 11\,\mathrm{Gyr}$ is required from the PTA observations. We note that for unequal-mass binaries, the GW radiative efficiency from a BH with a mass of qM, which is gravitationally captured in a circular orbit by a BH with mass of M, is given by _gw ≃ (0.057 + 0.4448η + 0.522η²)/(1 + q) ≃ 0.1192, where $\eta =q/{(1+q)}^{2}$ and q = 1/3 is set (Lousto et al. 2010). Approximating ${\epsilon }_{\mathrm{gw}}\simeq {\epsilon }_{r}$, we therefore obtain the interesting constraint that the contribution of BH mergers to the present-day BH mass density is less than $\lesssim 20\langle 1+z\rangle \% ;$ see Equation (7) in Phinney (2001).

The brightest U/HyLIRGs that have experienced active star formation at high redshift would be observed as massive elliptical galaxies in the local universe. An important consequence from Equation (7) is that ≳80% of the binary SMBHs formed in ULIRGs neither coalesce within a Hubble time nor contribute to the GWB. For SMBH binaries with mass ratios of q > 0.1, the DF caused by surrounding stars with velocity dispersion σ_⋆ would carry the binary separation down to

$$\begin{eqnarray}&&{a}_{{\rm{h}}}\sim \displaystyle \frac{12q}{1+q}\left(\displaystyle \frac{{M}_{\mathrm{BH}}}{{10}^{9}\,{\text{}}{M}_{\odot }}\right){\left(\displaystyle \frac{{\sigma }_{\star }}{300\mathrm{km}{{\rm{s}}}^{-1}}\right)}^{-2}\,\mathrm{pc},\end{eqnarray} \tag{ 11 }$$

where the binary becomes hard, ejects stars, and stalls the orbital decay (e.g., Merritt 2013). The formation timescale of an SMBH hard binary is typically ∼100 Myr (e.g., Yu 2002), which is much shorter than both a Hubble time and the delay time t_delay of interest. This suggests that a remnant population of O(1–10) pc binaries would be left at the centers of nearby massive ellipticals. Although no such binaries have been detected by PTAs to date, this nondetection has already yielded interesting constraints on their mass ratios (Schutz & Ma 2016) and anisotropy in the GWB (Mingarelli et al. 2017).

We thank Jeremiah Ostriker, Alberto Sesana, and Yoshiki Toba for useful discussions and comments. This work is partially supported by the Simons Foundation through the Simons Society of Fellows (K.I.), by JSPS KAKENHI 18K13584 and JST grant "Building of Consortia for the Development of Human Resources in Science and Technology" K.I.), and by NASA grant NNX17AL82G and NSF grant 1715661 (Z.H.).