Probing the Mass Relation between Supermassive Black Holes and Dark Matter Halos at High Redshifts by Gravitational Wave Experiments

Numerous observations have shown that almost all galaxies in our Universe host supermassive black holes (SMBHs), but there is still much debate about their formation and evolutionary processes. Recently, gravitational waves (GWs) have been expected to be a new and important informative observation, in particular, in the low-frequency region by making use of the Laser Interferometer Space Antenna (LISA) and Pulsar Timing Arrays (PTAs). As an evolutionary process of the SMBHs, we revisit a dark matter (DM) halo–SMBH coevolution model based on the halo merger tree employing an ansatz for the mass relation between the DM halos and the SMBHs at z = 6. In this model, the mass of SMBHs grows through their mergers associated with the halo mergers, and hence, the evolutionary information must be stored in the GWs emitted at the mergers. We investigate the stochastic gravitational background from the coalescing SMBH binaries, which the PTAs can detect, and also the GW bursts emitted at the mergers, which can be detected by the mHz band observations such as LISA. We also discuss the possibility of probing the mass relation between the DM halos and the SMBHs at high redshift by future GW observations.


INTRODUCTION
Supermassive Black Holes (SMBHs) typically have masses in the range of 10 6 M ⊙ ≲ M BH ≲ 10 10 M ⊙ and are found at the centers of almost all observed galaxies.The masses of SMBHs are strongly correlated with various features of the host galaxy, such as velocity dispersion, bulge mass, and luminosity (see, e.g., McConnell & Ma 2013).These observations suggest that SMBHs and their host galaxies have co-evolved to date.However, the formation and growth processes of SMBHs are still unknown.To clarify this conundrum, gravitational wave astronomy is expected to provide important clues in the near future through observations of GWs in the low-frequency region with the Laser Interferometer Space Antenna (LISA) and Pulsar Timing Arrays (PTAs).There are many theoretical predictions for the GWs from SMBH binaries.The event rate of SMBH coalescence is theoretically estimated by using phenomenological prediction (Jaffe & Backer 2003;Sesana 2013;Ravi et al. 2015;Sesana et al. 2016), semi-analytical models based on the extended Press-Schechter formalism or cosmological N-body simulation (Wyithe & Loeb 2003;Enoki et al. 2004;Sesana et al. 2004Sesana et al. , 2009;;Yang et al. 2019), or cosmological hydrodynamical simulations (Salcido et al. 2016;Kelley et al. 2017;Katz et al. 2020;Volonteri et al. 2020).
In addition to the relations between SMBHs and various features of host galaxies, there is also a strong correlation between the masses of SMBHs and their host dark matter (DM) halos.For instance, Ferrarese (2002) estimated the local M BH -M halo relation at z ∼ 0. Furthermore, Shimasaku & Izumi (2019) estimated M BH -M halo relation at z ∼ 6 using 49 Quasi-Stellar Object (QSO) datasets and compared it with the relation at z ∼ 0. They argued that these results, together with the relation to galaxy properties, provide clues to understanding how SMBHs grow with galaxies, but the relation at z = 6 was limited to QSOs with DM halo masses in the range of 10 11 M ⊙ ≲ M halo ≲ 10 13 M ⊙ .Their results are thus considered biased toward heavier masses than the most populated mass range M halo ≲ 10 12 M ⊙ .
Using a simple DM halo-SMBH coevolution model based on an extended Press-Schechter formalism that only considers the SMBH mass increase due to mergers, Bansal et al. (2023) constrained the M BH -M halo relation characterized by a power law index n in the region below 10 12 M ⊙ and a minimum halo mass M lim that can contain SMBH at z = 6.The purpose of this study is to extend their study further and discuss the possibility of restricting the M BH -M halo relation by extrapolating from the observed mass range, considering its effect on the GW background radiation from SMBH mergers, and considering the detectability of such GWs in future observations.In addition to their model, in this study, we also consider the SMBH coalescence timescale and the GWs generated by the coalescence of the SMBH binaries.The M BH -M halo relation at z = 6 allowed from GW observations enables us to constrain the SMBH formation process and its growth by gas accretion at high redshifts z ≳ 6.In particular, the parameter M lim may be closely related to the origin of the SMBH.For example, in a scenario where we consider the remnant of PopIII as the origin of the SMBH (King et al. 2008;Banik et al. 2019;Griffin et al. 2019;Zubovas & King 2021), we can consider a BH of about 10M ⊙ in a halo of about M lim ∼ 10 6 M ⊙ , and in a direct collapse scenario (Bromm & Loeb 2003;Chon & Omukai 2020), we can consider a BH of about 10 5 M ⊙ in a halo of about M lim ∼ 10 7 M ⊙ .
The structure of this paper is as follows.In Section 2, we explain our DM halo-SMBH coevolution model and calculate the SMBH mass evolution in this model.In Section 3, we introduce two types of GW radiation, namely Stochastic GW Background (SGWB) and GW burst, generated by SMBH binaries and discuss their detectabilities and implications to the model parameters M lim and n.Finally, in Section 4, we conclude our work and discuss future directions.
The aim of this paper is to investigate the impact of the initial DM halo-SMBH relation on the SGWB and GW bursts.As the DM halo-SMBH coevolution model, we take and modify the semianalytical model in Bansal et al. (2023), in which DM halos merge and evolve to satisfy the Press-Schechter mass function and only the mergers of SMBHs associated with the DM halo mergers lead to the mass evolution of SMBHs.Although it is simple, Bansal et al. (2023) has shown that the model can connect the M BH -M halo relations observed between at a high redshift (z ≈ 6) and the local universe (z = 0).In this section, we describe our DM halo-SMBH coevolution model and the initial DM-SMBH relation at z = 6.

Halo merger tree
The halo merger tree can provide the DM halo mass evolution in the hierarchical structure formation described in the Press-Schechter formalism.There are various methods to construct a halo merger tree (e.g., Zhang et al. 2008).In this paper, we adopt the binary merger tree method first introduced by Lacey & Cole (1993).This method is easy to calculate fast and suitable for describing SMBH binary mergers.We implement this method in our model following Somerville & Kolatt (1999).
In the binary merger tree method, all mergers are binary and the resultant halo mass after a merger is the sum of the two halo mass in the binary merger.Now we consider the probability of one progenitor halo with mass M p1 at redshift z = z 1 evolve to the descendent halo with mass M 0 (M p1 < M 0 ) at z = z 0 (z 0 < z 1 ).In the extended Press-Schechter formalism (Bond et al. 1991), such probability is given by Here ∆S ≡ S(M p1 ) − S(M 0 ) where S ≡ σ 2 (M ) is the mass variance of the overdensity field.In Equation (1), ∆δ c ≡ δ c (z 1 ) − δ c (z 0 ) with δ c (z) = δ c /D(z) where δ c = 1.69 is the threshold of spherical collapse and D(z) is the linear growth factor normalized as D(z = 0) = 1.This probability corresponds to the one of the binary merger in which the descendant halo with mass M 0 forms from two progenitor halos with masses M p1 and M p2 = M 0 − M p1 in the redshift range z 0 < z < z 1 .To construct the halo merger tree, we start from one realization of the halo mass distribution following the Press-Schechter mass function at z = 0. We take the redshift steps ∆z up to z = 6.At each redshift step, we consider the binary merger probability in Equation (1) for each existing DM halo and obtain two progenitors if the binary merger happens.We repeat this procedure at the subsequent redshift step until the redshift reaches z = 6.In the construction of the merger tree, we initially prepare 100 DM halos with mass in the range of 10 11 M ⊙ h −1 ≤ M halo ≤ 10 14 M ⊙ h −1 at z = 0 following the Press-Schechter mass function.In the calculation, ∆δ c is often taken to be constant for the reduction of the calculation cost.In this study, we take ∆δ c = 0.02 and we confirmed that the Press-Schechter mass function is realized at each redshift step with this value of ∆δ c .Note that by introducing lower cut-off mass M lim we ignore the existence of any DM halo with mass lower than M lim in our model.If a mass lower than M lim is selected as a progenitor mass in the merger tree construction, we assume that this mass is accreted onto the other halo progenitor as diffuse DM, and we do not take into account it as a DM halo.We will also discuss M lim in Section 2.3.

Initial mass relation between SMBHs and DM halos
Shimasaku & Izumi (2019) has reported the mass relation between the DM halos and SMBHs around z = 6 using QSO datasets.Motivated by this work, Bansal et al. (2023) suggested the power-law type of the mass relation between the DM halos and SMBHs at z = 6, where (M 0 , M 1 ) = (3.5 × 10 12 M ⊙ , 1.5 × 10 9 M ⊙ ).We adopt this power-law type of the mass relation characterized by a power-law index n as the initial mass distribution of SMBHs in our calculation.
As introduced in the previous section, M lim is the lowest mass of DM halos and, in our merger tree, the increment of mass lower than M lim is regarded as the accretion of ambient DM, which does not host an SMBH.Here, as the fiducial value, we set M lim at z = 6 to the virial mass with the virial temperature T vir = 10 5 K at z = 6 as adopted in Wyithe & Loeb (2003), M lim = M halo (T vir = 10 5 K) ≃ 3.5 × 10 9 M ⊙ h −1 .

DM halo-SMBH coevolution
Based on the halo merger history from z = 6 to z = 0 and the initial condition above, we consider the evolution of SMBH via SMBH mergers accompanied by their host DM halo mergers.Generally, the SMBH coalescence is delayed from the mergers of DM halos.In order to take into account the delay, we adopt the SMBH merger scenario following Volonteri et al. (2003).
At a merging time t mg , a lighter DM halo is captured into a heavier one.Then the lighter DM halo sinks to the center of the heavier DM halo due to dynamical friction.Although the shape of the lighter DM halo decays during the dynamical friction process, the lighter DM halo can keep its central SMBH.Therefore, when the lighter DM halo reaches the center of the heavier DM halo, we can assume that the coalescence between SMBHs that each DM halo hosts occurs.Hence we can take this dynamical friction timescale t df as the SMBH merger timescale 1 .We calculate t df as where P = M L /M H is the mass ratio of the lighter-heavier DM halos (M L ≤ M H ), r i is the initial radius, ϵ is the circularity, V circ is the circular velocity of the new halo with mass M L + M H , r vir is the virial radius, and ln Λ ≃ ln(1 + P ) is the Coulomb logarithm.We calculate V circ and r vir following Barkana & Loeb (2001), where We choose the set of constants in Equation ( 3), Θ = (r i /r vir ) 2 ϵ α = 0.3, following Volonteri et al. (2003).Then the coalescence of the SMBH binary occurs at time t mg + t df after the DM halo merger at time t mg .The mass of the final SMBH equals the total mass of two SMBHs.New System (Halo Mass:   +   ) Figure 1.Illustration for the merger of two systems (DM halos with SMBHs) in our model.The heavier system consists of a DM halo with mass M H (a dashed circle in red) and two SMBHs (black-filled circles with red outlines).One SMBH is at the center of the DM halos and the other SMBH is falling to the center.The lighter system includes a DM halo with mass M L (a dashed circle in blue) and an SMBH at the center (black-filled circles with blue outlines).Through the DM halo merger, a new system has a DM halo with mass M H + M L (a dashed circle in navy) and two SMBHs coming from the heavier system (one is at the center and the other is falling to the center).If P < 0.3, we neglect the SMBH from the lighter system in the new system (the left branch).Otherwise, we assume that the lighter system is falling into the center of the new system with the dynamical friction time scale determined by P and leads to the SMBH coalescence.

DM Halo
The timescale given by Equation ( 3) is determined by P .Here we introduce the critical value P * with which the dynamical friction time scale is roughly the same as the Hubble time scale at time t mg .If P of the merging DM halos is smaller than P * , the dynamical time scale exceeds the cosmological time scale and the coalescence of SMBHs does not happen in the cosmological time scale.Therefore, we neglect the SMBH coalescence with P < P * .If P is larger than P * , the dynamical time scale is short enough to induce the SMBH coalescence in the cosmological time.Therefore, we assume that the SMBH coalescence with P ≥ P * leads to the mass growth of the central SMBH of the heavier DM halo at time t mg + t df .In our redshift range, we confirmed that the dynamical friction time scale is always less than the Hubble time scale with P ≥ 0.3.The DM halo merger with P ≥ 0.3 corresponds to the so-called major merger.Through our simulations, we adopt the constant critical value P * = 0.3.In our model, the SMBH coalescence is delayed from the DM merger.Therefore, there is the possibility that, before the satellite SMBH coalesces with the central SMBH, the host DM halo has a merger with another DM halo.As shown in Figure 1, let us suppose that the heavier DM halo (red dashed circle) has one central SMBH and one satellite SMBH, while the lighter DM halo (blue dashed circle) has one central SMBH.After the DM halo merger at time t mg , the merged DM halo (navy) succeeds the central and satellite SMBHs in the heavier DM halo of the progenitors (red).If P = M L /M H ≥ 0.3, we regard the central SMBH in the lighter DM halo as the satellite SMBH in the merged DM halo which eventually coalesces with the central SMBH at t mg + t df (the right branch in Figure 1)2 .While, if P = M L /M H < 0.3, we neglect the merger between the central SMBHs of the two DM halos because they do not coalesce in the cosmological timescale (the left branch in Figure 1).Now we show the results of the SMBH mass growth in our model.Figure 2 represents the initial mass relation between the masses of DM halos and SMBHs at z = 6 with different power-law indexes in our calculations.The mass distribution of DM halos is obtained from the procedure we have described above.Then we assigned SMBHs to DM halos at z = 6, following Equation (2).For comparison, we also plot the observed mass relation between DM halos and SMBHs at z = 6 and z = 0 as in red and gray dots, respectively.We also represent the local mass relation reported in Ferrarese (2002) as the black dotted line.Some of the DM halo masses among observed DM halo-SMBH systems are larger than 10 12 M ⊙ at z = 6.Such heavier DM halos are rare objects in the Press-Schechter mass function and are not found in our simulation.The main mass range of DM halos in our simulations is less than 10 12 M ⊙ .
Figure 3 tells us the resultant mass relation between DM halos and SMBHs at z = 0 from the initial relation in Figure 2. In the left panel of Figure 3, we show the dependence on the power-law index n with the fiducial M lim .Although all of DM halo mass is below 10 12 M ⊙ at z = 6 initially, they evolve through the mergers and accretions in our simulations.And the resultant mass range expands to 10 14 M ⊙ at z = 0. SMBHs also evolve only through the major mergers of DM halos with P > 0.3.Relatively, larger DM halos have a larger possibility to experience a major merger than smaller DM halos and larger DM halos tend to host larger SMBHs.As a result, larger SMBHs exclusively have mass growth, while the mass of smaller SMBHs hardly grows.Accordingly the mass relation between DM halos and SMBHs becomes steeper as the universe evolves.We find out that n = 1.0, 1.4 and 1.8 at z = 6 shift to 1.3, 1.6 and 1.9 at z = 0, respectively.In our models, a smaller power-law index leads to the effective mass growth of SMBHs on the large mass side.
In the lower panel of Figure 3, we plot the mass relations for different M lim .In contrast to n, the difference of M lim does not affect the evolution of the mass relation.Although changing M lim just modifies the total SMBH mass, it does not change the mass dependence of the SMBH mass growth.
At the last of this section, it is worth mentioning about the SMBH growth via gas accretions.Although we omit it in our model for simplicity, the gas accretion onto an SMBH is another major contribution for the mass growth of SMBHs.Therefore, we underestimate the final SMBH mass at z = 0. Including the mass growth via the gas accretion shifts up the SMBH mass z = 0 in the all of mass range.Although there exists the gap between the observed local relations and the relations in our simulations in Figure 3, the gas accretion onto SMBHs can decrease this gap.

GRAVITATIONAL WAVES FROM SMBH BINARIES IN THE COEVOLUTION MODEL
Based on the DM halo-SMBH coevolution model introduced in the previous section, we investigate GWs from SMBH binaries in the process of coalescing evolutions.First, we calculate the SGWB that could be mainly produced by the superposition of a large number of inspiraling phases and discuss the constraint from the PTA experiments on the power-law index, n, which has been introduced in Equation (2).Then, we also estimate the expected event rate of SMBH mergers in the LISA experiment.

Stochasic gravitational wave background
To describe the energy density of the SGWB per logarithmic frequency interval, dρ GW /d ln f , we use the characteristic strain amplitude h c (f ) which is given as (see, e.g., Maggiore 2018) and the SGWB from SMBH binaries can be estimated as Here, d 3 n c /dzdM 1 dM 2 is the comoving number density of the coalescing SMBH binaries with masses M 1 ∼ M 1 + dM 1 and M 2 ∼ M 2 + dM 2 at z ∼ z + dz, which can be calculated in the DM halo-SMBH coevolution model discussed in the previous section, and dE GW (f, z, M 1 , M 2 )/d ln f is the energy spectrum of the GW radiation from a coalescing SMBH binary with masses of M 1 and M 2 at the redshift z.As a characteristic frequency of the GWs from such a BH binary, f ISCO , which is the GW frequency corresponding to the orbital angular frequency of the innermost stable circular orbit (ISCO), can be given by where M tot = M 1 + M 2 .For the lower frequencies than the f ISCO the coalescing BH binaries can be considered to be in the inspiral phase, and the most sensitive frequency range of the PTA experiments is about 10 −9 ∼ 10 −8 Hz.As we have shown in the previous section, the mass range of the SMBHs in our setup is typically 10 4 ∼ 10 10 M ⊙ (see Figures. 2 and 3).Thus, for the PTA experiments, we can use the energy spectrum of the GW radiation in the inspiraling phase, which is given by (see, e.g., Maggiore 2018) where is the chirp mass of the binary.For simplicity, we have assumed the circular orbit and dE GW /d ln f = 0 for f > f ISCO .Then, from the above formula, we can calculate the characteristic strain amplitude of the SGWB from SMBH binaries in our DM halo-SMBH coevolution model and it basically depends on two model parameters, n and M lim , which have been introduced in the ansatz given by Equation (2).
Figure 4 shows the characteristic strain amplitudes for the cases with n = 1.0 (left) and 1.4 (right), respectively.The upper panels show the contribution from each mass bin in the amplitude, and the lower ones show the contribution from each redshift bin in the amplitude.In each panel, the total amplitude is shown by a black line.From this figure, one can find that the contribution from the coalescing SMBHs with masses in the range of 10 7 − 10 9 M ⊙ at 0 ≤ z ≤ 1 is dominant.This is simply because the larger system can emit the larger amplitude of GWs, as shown by Equation ( 7).While there is a contribution from the SMBH binaries with M BH ≳ 10 9 M ⊙ in the case with n = 1.0, it is small.This is due to the small number of such heavier SMBHs, which can be seen in Figure 3.As for the dependence on M lim , as shown in Figure 3 the population of the SMBHs at the lower redshift is almost independent of M lim , and Figure 4 shows that the dominant contribution in the SGWB is coming from the coalescing SMBHs at the redshifts 0 ≤ z ≤ 1.Thus, from these facts, the amplitude of the SGWB from the coalescing SMBHs in the frequency band of the PTA experiments (f ∼ 10 −9 − 10 −8 Hz) is almost independent of M lim .
In Figure 5, we show the comparison of the theoretical predictions in our coevolution model and the PTA observational results.The left panel shows the theoretical predictions of the characteristic strain in the cases with n = 1.0 (blue), 1.4 (orange), and 1.8 (green) for M lim = 3.5 × 10 9 M ⊙ h −1 , and the upper bound obtained from the PPTA project (Shannon et al. 2015) is represented by a triangle.We also include the possible signal of SGWB corresponding to the common-spectrum process reported from the NANOGrav experiment (Arzoumanian et al. 2020) represented by the red band and the sensitivity curve expected for the future PTA experiment by SKA with the black dashed-dotted line (see, e.g., Campeti et al. 2021).The right panel shows the contours of the amplitude of the characteristic strain at f = 1yr −1 on the n-M lim plane.First, from the left panel, one can find that the amplitude of SGWB becomes smaller as n increases, which can be also seen in Figure 4.This is simply because the typical mass of SMBHs in the lower redshift becomes smaller for larger n as shown in Figure 3.As the possible signal of SGWB reported by the NANOGrav experiment, which is represented by the red band, we employ  with A = 1.92 +0.75 −0.55 .The right panel tells us that the current observational results obtained from the PTA experiments imply n ≳ 1.2 and this implication is almost independent of the value of M lim .This means that the PTA experiments can provide information about the mass relation between the DM halos and the SMBHs at z = 6.

GW burst
Let us discuss the detectability of the GW bursts from the individual SMBH mergers in our coevolution model by using future planned GW experiments such as LISA.After the inspiraling phase, which has been discussed above, the SMBH binaries finally enter the merging phase and emit strong burst-like GWs.The waveform of GWs emitted during this phase is difficult to estimate analytically, and is, in general, estimated by calculations based on numerical relativity.Here, we employ IMRPhenomenD (Khan et al. 2016) using PyCBC3 to generate the Fourier waveform of a burst.
The criteria for whether a GW burst is detected is closely related to the signal-to-noise ratio (SNR) which is defined as (Babak et al. 2021) where A is the amplitude of the burst in Fourier space and S SciRD (f ) is the official required sensitivity of LISA introduced in the LISA Science Requirements Document, given as (Babak et al. 2021) We adopt ρ ≥ 5 as the detection criteria.Thus, the number of detectable GW bursts per year coming from z ∼ z + dz can be calculated as where d(z) is the comoving distance to the redshift z and Φ(z, M 1 , M 2 ) equals to 1 when ρ ≥ 5 , otherwise 0.
Figure 6 shows the distribution of the expected number of the GW bursts which would be detected in the LISA with respect to the peak frequency.The left panel shows the result with changing the power-law index as n = 1.0 (blue), 1.4 (orange), and 1.8 (green).In the right panel, we change the limiting mass, M lim from 1.0 × 10 8 M ⊙ h −1 to 1.0 × 10 10 M ⊙ h −1 .As we have discussed in Sec. 2, in our DM-halo and SMBH coevolution model, we assume the power-law relation given by Equation 2 at z = 6 as the initial condition, and then our calculation can not trace the halo merger history before z = 6 and the evolutionary history of the SMBHs.Based on our criteria for the coalescing SMBHs discussed in the previous section, the coalescing SMBH binaries associated with the mergers of halos which occur at z ≥ 6 are expected to merge at z ≥ 3. Thus, to obtain a conservative estimation, here, we only count the burst events which occur at z ≤ 3.In our coevolution model, as can be seen in Figure 2, as n increases small halos host the SMBHs with smaller mass.In general, smaller halos exist in larger numbers, and then for the larger n, the typical mass of the SMBHs in our Universe becomes smaller.Since the peak frequency of the burst is inversely proportional to the mass of the binary system, which is roughly given by f peak ≃ 1.9 × 10 −4 Hz 1 1 + z the distribution of the burst events would shift toward the higher frequency as n increases and eventually goes outside the LISA band.As can be seen in the left panel of Figure 6, therefore, the total number of detectable events decreases as n increases.As for the dependence on the limiting mass shown in the right panel, unlike the case of SGWB, the expected number of events also depends on the limiting mass.Basically, the increase in the limiting mass means that the number of halos hosting the SMBHs decreases, and then the number of burst events should also decrease.
In Figure 7, we also plot the contours of the expected number of burst events in the LISA experiment in the n-M lim plane.Basically, from this figure, once the number of burst events from z ≤ 3 is specified, we can obtain the constraint on the mass relation between the DM halos and SMBHs parameterized by n and M lim .However, note that in our model, we only take the mergers into account for the mass evolution of SMBHs, neglecting the gas accretion onto the BHs.Naively, by including the gas accretion effect, the SMBHs are expected to become more massive, and then, in terms of the GW bursts in the LISA band, the expected number of events could be larger compared with that showed in the Figure 7, which is estimated in the case of mass growth only through mergers.In this sense, our model basically gives a conservative estimation for the number of burst events, and it means that we can exclude the region in this figure where the expected number is larger than the observed one.As an example, this figure shows that if we would detect 36 burst events by LISA we could exclude the region where n ≲ 1.6 and M lim ≲ 6.0 × 10 8 M ⊙ h −1 .Unlike the case of SGWB, the LISA experiment can give us information about the limiting mass and it should be a key of unveiling the formation process of the SMBHs at high redshifts.

DISCUSSION AND CONCLUSION
In this paper, we have investigated the potential of GW observations to probe the mass relation between DM halos and SMBHs.The mass relation in high redshifts could be a key to revealing the DM halo-SMBH coevolution and the physics of the SMBH mass growth.
We have constructed the DM halo-SMBH coevolution model according to the merger tree based on the extended Press-Schechter formalism.In our model, the SMBH mass grows through the SMBH coalescence triggered by the DM halo mergers.The initial mass relation of SMBHs and DM halos is given as a power-law type relation with two model parameters; one of the model parameters is the power-law index of the relation at z = 6, n, and the other is the lowest mass of DM halos that can host SMBHs, M lim .Our model can roughly connect the mass relation observed in high redshift (z ∼ 6) to that in the present epoch.
To account for the dynamical friction, we consider the time delay of the SMBH coalescence after the DM halo mergers in our model.This method allows us to treat multiple SMBHs and coalescences in a single DM halo, although the binary merger tree method is used for the DM halo mergers.
Using our SMBH-DM halo coevolution model, we have evaluated two types of GW signals; the GW background in the ∼ nHz frequency range and the GW burst.The former can be probed by PTA experiments and the latter can be detected by the future space gravitational observatory, LISA.
We have shown that the amplitude of the GW background strongly depends on the power-law index of the initial mass relation, n, while it is almost immune to M lim .The current PTA experiments provide the upper limit on the amplitude of the GW background.By combining the initial power law mass relation with the recently observed mass relation around z ∼ 6, we find that the PTA constraint on the GW background gives a limit to the power-law index in the initial mass relation as n ≳ 1.2.We have also demonstrated that near future PTA experiments such as SKA can provide a stronger constraint on n.
For the GW bursts from SMBH coalescences, we have calculated the number of detectable bursts per year with LISA sensitivity.Our results show that the number of detectable GW bursts depends on both n and M lim .Therefore, once we obtain an upper limit or constraint on the number of GW bursts per year, we can impose the limits on n and M lim from the LISA observation, which can then be used to constrain the formation and evolution processes of SMBHs at high redshifts.
Future galaxy and QSO surveys can provide constraints on the mass relation between SMBHs and DM halos in high redshifts.However, in comparison with these observations, GW observations can constrain the mass relation even for the smaller mass region of DM halos and SMBHs which galaxy and QSO surveys can not probe.On the other hand, it is important to construct a precise coevolution model of DM halos and SMBHs to obtain the robust constraint on the mass relation in high redshifts.
In our coevolution model, we do not take into account the gas accretion on SMBHs, which is one of the important processes for SMBH mass growth.Therefore, in our model, the final masses of the SMBHs are underestimated.Since the SMBH mass function shows that the smaller mass BHs are more abundant, including a gas accretion in our model leads to enhancing the number of the SMBH coalescence which contributes to the detectable GW signals.Accordingly, the constraint on the initial conditions of the mass relation in our model could be underestimated.Although it is difficult to model the gas accretion in high redshifts, the detailed observation of the luminosity function of distant QSOs can provide us the information on the efficiency of the gas accretion in high redshifts.We will address these issues in our future work.
In this work, we do not take into account the detailed parameters of the SMBH merger, e.g., the spins and eccentricity of BH binaries.The spins of the BH binary are related to the final BH recoil.The strong final BH recoil in the SMBH coalescence can eject one of the two SMBHs from the host DM halo.The eccentricity of the BH binary affects the spectrum of the GW radiation.Larger eccentricity leads to a decrease in the amplitude of the GW radiation in the PTA frequency range (Enoki & Nagashima 2007;Ravi et al. 2014).Therefore, such parameters of the SMBH merger affect the amplitude of the GW background or the event rate of the GW bursts.Furthermore, if triple BH interaction occurs in a system, one of three BHs will be ejected, and then the radius of the other two BHs will shrink and lead to their coalescence (Volonteri et al. 2003;Bonetti et al. 2018).Thus, triple BH interactions also have a possibility to affect the SMBH merger rate.The orbital evolution via the three-body interaction between the SMBH binary and stars in the system affects the spectrum in the low-frequency (Sesana et al. 2004).As future work, we should consider how these affect SMBH evolution and future GW observations.where a is the separation of two coalescing BHs, M * is the total stellar mass of the host galaxy, σ * = (0.25GM * /R eff ) 1/2 is the velocity dispersion, R eff = 0.1r vir , M BH,L is the mass of the lighter BH, and Λ = ln(1 + M * /M BH,L ).
Employing this timescale for the SMBH mergers, in particular, the merger timescale for lighter BHs would be longer, compared with the timescale given by Equation (3), and some binaries with small masses could not merge within cosmic time.Thus, the resultant GWs in the coevolution model can be relatively sensitive to the masses of SMBHs and also the power-law index n.The results with the timescale given by Equation (A1) are shown in Figure 8.Here, we assume that M * = (Ω b /Ω m )M halo , where M halo = M H + M L is the total halo mass of the new system (see Figure 1), and a = R eff for simplicity.In this figure, the left panel shows the contours of the amplitude of SGWB in n-M lim plane.The amplitude of the SGWB is dominated by the contribution from the large mass BHs, and no matter which time scale is adopted, they always coalesce within the time scale considered, and therefore, the results are not significantly affected.The right panel of Figure 8 shows the event rate of GW bursts for each n and M lim .For the GW bursts, the contribution from small mass BHs dominates, and hence the change in timescale seems to have an impact, especially when n is large and small mass BHs are more abundant.However, in fact, when n is large, the expected number of GW busts was small, and therefore we find that the change in timescale does not have a large impact on the results.

Figure 2 .
Figure 2. The mass relation between DM halo and SMBHs at z = 6.We plot our initial mass relation obtained from the halo merger tree and Equation (2) with M lim = 3.5 × 10 9 M ⊙ h −1 as colored triangles.The different color represents the different power-law index n in Equation (2); n = 1.0 in blue, n = 1.4 in orange and n = 1.8 in green.The faint black and red small circles represent the observed mass relation at z = 0 and z = 6 (Shimasaku & Izumi 2019), respectively.We also plot the local relation given by Ferrarese (2002) in a black dotted line.

Figure 3 .
Figure 3.The mass relation between DM halos and SMBHs at z = 0.The upper three panels show the dependence on the power-law index n with the fiducial M lim = 3.5 × 10 9 M ⊙ h −1 and the lower three panels show the dependence on the minimum halo mass M lim with n = 1.4.In each panel, colored circles represent the mass relations in the simulation and the colored dashed line is the fitting power-law line for the simulation result at z = 0.The faint small black circle and the black dotted line are the same as in Figure 2. Also, in the upper panels, we plot the initial mass relation at z = 6 with the colored triangles.

Figure 4 .
Figure 4.The characteristic strain amplitude of the SGWB from SMBH binaries for the cases with n = 1.0 (left) and 1.4 (right).M lim = 3.5 × 10 9 M ⊙ h −1 is fixed.The upper and lower panels show the contribution from each mass bin and each redshift bin, respectively.

Figure 5 .
Figure5.Comparison of the theoretical predictions in our coevolution model and the PTA observational results.The left panel shows the theoretical predictions of the characteristic strain in the cases with n = 1.0 (blue), 1.4 (orange), and 1.8 (green) for M lim = 3.5 × 10 9 M ⊙ h −1 , and the upper bound obtained from the PPTA project(Shannon et al. 2015) is represented by a triangle.We also include the possible signal of SGWB corresponding to the common-spectrum process reported from the NANOGrav experiment(Arzoumanian et al. 2020) represented by the red band and the sensitivity curve expected for the future PTA experiment by SKA with the black dashed-dotted line (see, e.g.,Campeti et al. 2021).The right panel shows the contours of the amplitude of the characteristic strain at f = 1yr −1 on the n-M lim plane.

Figure 6 .
Figure 6.The distribution of the expected number of the burst events which would be detected in the LISA with respect to the peak frequency.The left panel shows the dependence on the power-law index in our coevolution model with fixing M lim = 1.0 × 10 8 M ⊙ h −1 .The right one shows the result with changing the limiting mass for the case with n = 1.0.

Figure 7 .
Figure 7.The contours of the expected number of burst events in the LISA experiment in the n-M lim plane.The red line denotes the boundary of the excluded region for 36 burst events per year detected by LISA.

Figure 8 .
Figure 8. Theoretical predictions of SGWB and the event rate of GW bursts for the different timescale given by Equation (A1).