Exploring the multi-band gravitational wave background with a semi-analytic galaxy formation model

An enormous number of compact binary systems, spanning from stellar to supermassive levels, emit substantial gravitational waves during their final evolutionary stages, thereby creating a stochastic gravitational wave background (SGWB). We calculate the merger rates of stellar compact binaries and massive black hole binaries using a semi-analytic galaxy formation model – Galaxy Assembly with Binary Evolution (GABE) in a unified and self-consistent approach, followed by an estimation of the multi-band SGWB contributed by those systems. We find that the amplitudes of the principal peaks of the SGWB energy density are within one order of magnitude Ω 𝐺𝑊 (cid:24) 10 (cid:0) 9 (cid:0) 10 (cid:0) 8 . This SGWB could easily be detected by the Square Kilometre Array (SKA), as well as planned interferometric detectors, such as the Einstein Telescope (ET) and the Laser Interferometer Space Antenna (LISA). The energy density of this background varies as Ω 𝐺𝑊 / 𝑓 2 (cid:157) 3 in SKA band. The shape of the SGWB spectrum in the frequency range (cid:24) » 10 (cid:0) 4 , 1 … Hz could allow the LISA to distinguish the black hole seed models. The amplitude of the SGWB from merging stellar binary black holes (BBHs) at (cid:24) 100 Hz is approximately 10 and 100 times greater than those from merging binary neutron stars (BNSs) and neutron-star-black-hole (NSBH) mergers, respectively. Note that, since the cosmic star formation rate density predicted by GABE is somewhat lower than observational results by (cid:24) 0 . 2 dex at z < (cid:24) 2 , the amplitude of the SGWB in the frequency range (cid:24) » 1 , 10 4 … Hz may be underestimated by a similar factor at most.


INTRODUCTION
The present and forthcoming experiments aimed at detecting gravitational waves (GWs), including the advanced Laser Interferometer Gravitational wave Observatory (aLIGO) (Harry & LIGO Scientific Collaboration 2010), the Einstein Telescope (ET) (Punturo et al. 2010), the Laser Interferometer Space Antenna (LISA) (Amaro-Seoane et al. 2017), and the Square Kilometre Array (SKA) (Dewdney et al. 2009), herald the advent of multi-band GW astronomy.Joint detection by multiple GW experiments holds the potential to realize the goal of exploring the multi-band stochastic GW background (SGWB), which arises from the energy background resulting from the superposition of various resolved and unresolved GW sources spanning a broad range of frequencies.These sources can be classified into two categories, namely, cosmological sources and astrophysical sources.The joint detection of the multi-band SGWB from cosmological sources has been extensively investigated (e.g.Barish et al. 2021;Campeti et al. 2021).Here, we concentrate on the joint ★ E-mail: xilong.fan@whu.edu.cndetection of the multi-band SGWB originating from astrophysical sources.
In the near future, joint detection of the SGWB in different bands could provide an opportunity to assess the reliability of semi-analytic galaxy formation models and binary population synthesis models in a way that is independent of traditional electromagnetic observations.While detectors such as ET and LISA could directly detect loud GW events (e.g.Sesana et al. 2004;Rosado 2011;Périgois et al. 2021;Boco et al. 2021), the shape and amplitude of the total SGWB provide additional information to constrain the population properties of source models.A previous study by Rosado (2011) investigated the multi-band SGWB from astrophysical sources using independently calculated merger rates of massive black hole binaries and stellar binaries.In this paper, we employ a unified and self-consistent approach to calculate the merger rates of compact binary systems ranging from stellar to supermassive levels.We derive the merger rates of stellar compact binaries and massive black hole binaries from a semi-analytic galaxy formation model called Galaxy Assembly with Binary Evolution (GABE) (Jiang et al. 2019), combined with a rapid binary population synthesis model COSMIC v3.3.0 (Breivik et al. 2020).We comprehensively assess the influence of cosmic star formation rate density (SFRD) on the SGWB by incorporating empirical constraints for the SFRD.We also improve the completeness of the sample of intermediate-mass (100 ⊙ <  < 10 6  ⊙ ) black holes by utilizing the Millennium-II simulation (Boylan-Kolchin et al. 2009) in place of the Millennium simulation (Springel et al. 2005) used in Jiang et al. (2019).
The paper is organized as follows: firstly, the basic equations employed for the estimation of the SGWB signal are presented in Section 2. Next, the methodology utilized for generating mock samples of GW sources is demonstrated in Section 3. In Section 4, the SGWB signals are computed from different astrophysical sources which are expected to be distributed across a broad range of frequencies, followed by a presentation of the main findings.Finally, the primary conclusions drawn from the study are summarized in Section 5.

Merger rates and SGWB
The merger rate of compact binary systems, which denotes the event number per comoving volume per cosmic time at redshift , can be described as: where N () is the number density (i.e., the number of GW sources in per comoving volume per unit redshift at redshift ).Throughout the paper, we adopt a fiducial cosmological model with Ω m = 0.25, Ω Λ = 0.75 and  0 = 73 km/s/Mpc (WMAP1, Spergel et al. 2003), which is also adopted in the Millennium simulation.Then, one can get the relation between the redshift  and the cosmic time , i.e., d/d The dimensionless energy density of the SGWB, Ω GW , is defined as: where  c = 3 0 2  2 /8 is the critical energy density,  GW is the present-day energy density of GWs.According to Phinney (2001), where d (   ) /d (   ) is the energy spectrum of a single GW source in logarithmic interval, and  r = (1 + )  is the frequency in the cosmic rest frame of the GW source,  is the observed frequency.Following Ajith et al. (2008) and Zhu et al. (2011b), we adopt the GW energy spectrum with where  c = ( 1  2 ) 3/5 ( 1 +  2 ) −1/5 is the chirp mass,  1 and  2 are the masses of the primary and secondary stars in the GW source; ring are constants to make spectrum continuous at the boundary;  merg ,  ring and  cut are characteristic frequencies to divide the different stages of the compact binaries.The parameters  merg ,  ring ,  and  cut can be calculated with the formula Ajith et al. (2008), where  =  1  2 / 2 is the symmetric mass ratio with  =  1 +  2 .The values of the coefficients [ 1 ,  2 ,  3 ] are [2.9740× 10 −1 , 4.4810 × 10 −2 , 9.5560 × 10 −2 ], [5.9411 × 10 −1 , 8.9794 × 10 −2 , 1.9111 × 10 −1 ], [5.0801 × 10 −1 , 7.7515 × 10 −2 , 2.2369 × 10 −2 ], and [8.4845 × 10 −1 , 1.2848 × 10 −1 , 2.7299 × 10 −1 ] for  merg ,  ring ,  and  cut , respectively, which are also presented in the Table I of Ajith et al. (2008).The Eq. ( 4) is developed to describe inspiral-mergerringdown stages for non-spinning merging black hole binaries on circular orbits, but we assume that for BNSs and NSBHs, the spectrum is still suitable, as in previous studies (e.g.Chen et al. 2019;Capurri et al. 2021).
Two other quantities are also used to describe the SGWB besides Ω GW (  ): the one-sided spectral power density  h (  ), which is related to the detection criterion (e.g. the signal to noise ratio) , and the characteristic strain ℎ  (  ).They are related to Ω GW (  ) with the expression: The dimensionless characteristic strain  yr −1 is usually introduced in the pulsar timing array (PTA) experiment, and the relation between where the equality assumes a power-law form for the background and  yr = 1 −1 is the reference frequency.Note that for the SGWB from supermassive ( > 10 6  ⊙ ) binary black holes, in the frequency ∼ 10 −9 − 10 −7 Hz,  ≈ −2/3 (e.g.Zhu et al. 2019).
The merger rate () and the energy density of SGWB Ω GW (  ) are both related to the number density N () as presented in Eqs. ( 1) and ( 3).An analytical method to calculate N () (and then () or Ω GW (  )) is carried out by using an analytical model for SFRD and metallicity evolution history (e.g.Zevin et al. 2020) and combining with the properties of merging compact binaries.Alternatively, these quantities also can be computed by using the semi-analytic galaxy formation model.The later method is employed in this work as discussed in Section 3. In the semi-analytic method, one can compute N () with where Δ is the number of merging stellar compact binaries in the redshift interval Δ, and Δ  is the volume of the simulation box.
Eq. ( 3) can be rewritten as: where 'num' is the total number of the GW sources.Furthermore, Eq. ( 8) can also be calculated in terms of ,  1 , and  2 bins as: where  and  are the bin numbers of  1 and  2 , respectively, which are determined by the division of the  1 ×  2 space,  is the bin number of redshift , which is determined by the division of redshift, Δ is the number of merging stellar compact binaries in the redshift-primary mass-secondary mass interval Δ,  1 ,  2 .In this work, we assume that the values of d (   ,  1 ,  2 ) /d   are equal to their values at the bin center.

Signal to Noise ratio
We first describe the formalism we used to calculate the signal to noise ratio (SNR) for ground-based and space-based detectors.The expected SNR for a cross-correlation search for an unpolarized and isotropic SGWB is given by Allen & Romano (1999); Thrane & Romano (2013).Once the SGWB is calculated, we can estimate the SNR for the given noise power spectral density  n and overlap reduction function Γ   (  ).Γ   (  ) describes the reduction in sensitivity to the SGWB of two detectors labeled by  and  due to their separation and non-optimal orientations (Flanagan 1993) (when  = , we define the transfer function   (  ) = Γ   (  )).A normalized overlap reduction function: is often used to ensure    (0) = 1 for two identical, co-located and co-aligned detectors, where  is the opening angle of the interferometer.In this work, we estimate the SNR with two detectors for both LIGO and ET.Assuming the two detectors are almost co-located, the (  ) sin 2  for LIGO can be approximated to 1 (Nishizawa et al. 2009;Fan & Zhu 2008) and for ET, it can be approximated to −3/8 (Amalberti et al. 2022) over the frequency range we studied.For ground-based detectors, where  is the observation period of the experiment,  max and  min are the upper and lower sensitivity bounds of the experiment.We adopt the  n of aLIGO and ET from the LIGO document 1 .For 1 https://dcc.ligo.org/LIGO-T1500293/publicLISA, the √ 2 in the Eq. ( 11) can be reduced due to the use of data from only one detector (Thrane & Romano 2013): where   (  ) =   (  )/  (  ) is the strain spectral sensitivity.We calculated   (  ) using the analytical formulas presented in Robson et al. (2019), considering the galactic confusion noise.
For SKA, we use the python package ℎ (Hazboun et al. 2019) to calculate the SNR: where  Nyq = 1/(2Δ), Δ is the minimum time interval between two observations. eff (  ) is the effective strain-noise power spectral density for the whole PTA in which a single distinct pair is labeled by  and : where    are the Hellings and Downs factors, is the individual pulsar strain-noise power spectral density, N −1  (  ) is the inverse-noise-weighted transmission function and   (  ) = 1/(12 2  2 ).

Power-law integrated sensitivity curves
Power-law integrated (PI) sensitivity curves are often used to estimate the detection abilities of various detectors, which takes into account the enhancement in sensitivity that comes from integrating over frequency (Thrane & Romano 2013).Assuming a set of powerlaw SGWBs Ω  (  /  ref )  with indices e.g.,  = {−8, −7, ..., 7, 8} and an arbitrary reference frequency  ref , the PI sensitivity curves are the envelope of these SGWBs.

Merging massive black hole binaries
Utilizing the Millennium (MSI) and Millennium-II (MSII) simulations in conjunction with a semi-analytic galaxy formation model-GABE, we simulate the emergence and growth of galaxies and their associated black holes in a cosmological context.While the MSII simulation can resolve all luminous galaxies and their associated black holes in the universe, it lacks sufficient volume to provide adequate statistics to model extremely large supermassive black holes, which can lead to noisy predictions of ultra-low frequency SGWB.To address this limitation, we leverage the MSI, which has a volume 125 times greater than that of MSII, to model ultra-low frequency background.
As is common in many galaxy formation models, GABE assumes that a seed black hole is responsible for the formation of the central black hole of a galaxy.However, the exact mass of the seed black hole is not well-constrained by observations (see Volonteri 2010;Woods et al. 2019 for reviews).In this study, two different seed models are considered: the light-seed and heavy-seed models.The light-seed model proposes that the seed black holes are formed through the collapse of Population III stars (e.g.Bond et al. 1984;Madau & Rees 2001;Pezzulli et al. 2016), while the heavy-seed model suggests that they are the product of the direct collapse of atomic gas in dark matter halos with a virial temperature of about 10 4 K (e.g.Bromm & Loeb 2003;Begelman et al. 2006;Lodato & Natarajan 2006;Regan et al. 2017).In practice, the light-seed model involves the addition of a 10 2  ⊙ black hole to a galaxy at emergence, while the heavy-seed model adds a 10 4  ⊙ black hole to a galaxy when it first appears.Notably, the choice of seed black hole mass has a negligible effect on subsequent stellar evolution or the final black hole mass in the semi-analytic models utilized in this study, as the seed mass only represents a small fraction of the final black hole mass (≳ 10 6  ⊙ ).Additionally, any initial differences in seed mass are quickly erased by the following rapid growth during the first gas-rich galaxy merger.
The growth model of massive black holes adopted in GABE is similar to that of Croton et al. (2006).In the early stages, seed black holes primarily grow through the process of merging with other black holes and accreting cold disk gas during galaxy mergers, which is commonly referred to as the 'quasar mode'.The growth during a merger event is: where  BH,sat is the central black hole mass of the satellite galaxy,  sat and  cen are the masses (cold gas and stars) of the satellite and the central galaxy, respectively,  cold is the total mass of the cold gas in the satellite and central galaxies,  vir is the virial velocity of the central dark matter halo, and  = 0.03 is a parameter that describes the growth rate of the black hole.Once these black holes reach a sufficient mass (such as > 10 7  ⊙ ), they begin to accrete hot gas from the host galaxy in a steady and continuous manner, releasing an enormous amount of energy into the surrounding circumgalactic medium.This process is known as the 'radio mode' and is associated with the suppression of star formation.Such growth rate via absorbing hot gas is: where  is a parameter that describes the intensity of absorption,  hot is the mass of the hot gas, and  vir and  vir are the virial mass and velocity of the halo, respectively.While the former mode is the main driver of massive black hole growth in the early stages, the latter dominates at lower redshifts.
In GABE, the assumption is made that two massive black holes merge instantaneously when their host galaxies merge.This assumption, however, is an oversimplification, as the actual process of two massive black holes dissipating their orbital energy and sinking to the center of the potential well before coalescence involves several physical processes.These processes include dynamical friction with the stellar and gaseous background, three-body interactions with core stars, gravitational and viscous torques of the circumbinary disk, and the emission of GWs when the separation decreases to ∼ 0.01 pc scale.These dissipative processes can span a few billion years (Yu 2002;Vasiliev et al. 2015;Sesana & Khan 2015).Despite this, the disk and bulge configuration of galaxies considered in GABE are too simple to account for these physical processes.Therefore, as did in Yang et al. (2019) and Curyło & Bulik (2022), no time delay between massive black hole binary mergers and galaxy mergers is considered in our work.The Chabrier (Chabrier 2003) initial mass function (IMF) is adopted in GABE, while the rapid binary population synthesis model COSMIC v3.3.0 adopts the Kroupa IMF (Kroupa 2001).However, the Kroupa IMF and Chabrier IMF are similar for  * > 1 ⊙ , so we assume no calibration is required in the SGWB calculation especially in Eq. ( 21).

Merging stellar compact binaries
We have utilized the rapid binary population synthesis model COS-MIC v3.3.0 (Breivik et al. 2020) to derive properties of simple stellar populations (SSPs).The model is based on the binary star evolution code presented in Hurley et al. (2002), and has been substantially enhanced to incorporate prescriptions for massive star evolution and binary interactions.The primary objective of COSMIC is to simulate compact binaries and their progenitors, making it particularly well-suited to our research goals.
For every SSP with metallicity , we calculate event numbers of stellar compact binary mergers per solar mass Δ SSP, (,  delay ,  1 ,  2 ).Here,  represents categories of merging stellar compact binaries: BNSs, NSBHs, and BBHs;  1 and  2 denote the masses of the primary and secondary stars, respectively.We divide  1 and  2 into equal-width bins in the ranges of 1.24 ⊙ − 3 ⊙ for neutron stars and 3 ⊙ − 45 ⊙ for black holes.The delay time  delay is the duration between the merger time and the formation time of their progenitor zero age main sequence binaries.We consider the cases of  = 0.0001, 0.0003, 0.001, 0.004, 0.01, 0.02, and 0.03 in our work.
The COSMIC initialization is configured based on the study by Zevin et al. (2020), wherein Table 1 lists various options for the initialization, and we adopt the [Initial conditions, CE survival, CE efficiency, Remnant mass, Natal kicks] = [S+2012, Pessimistic, 1.0, Delayed, Bimodal] option for this investigation.Fig. 1 presents the delay time distributions of the BNSs, NSBHs, and BBHs mergers in the context of SSPs.For all types of merging compact stellar binaries and metallicities, the event rate  SSP follows an almost power law distribution ∝  −1 delay after reaching its peak, as several studies have previously pointed out (e.g., Dominik et al. 2012).
To self-consistently derive the distribution of merging stellar compact binaries in the universe, we apply our SSPs to galaxy catalogues generated by GABE.In GABE, a star formation activity is triggered if the gaseous disk is too massive to be stable or a galaxy merger occurs.Galaxies exhibit star formation events triggered by the disk instability that occurs when the cold gas within their disks surpasses a critical mass (Croton et al. 2006): where  max represents the maximum circular velocity of the dark matter halo, and  gas,d denotes the scale radius of the cold gas disk.Subsequently, the cold gas within the disk undergoes conversion into stars at a rate denoted by  * : where  = 0.02 is a parameter which describes the star forming efficiency and  dyn = 3 gas,d / max is the characteristic timescale at the edge of the star-forming disk.When galaxies undergo mergers, multi-band gravitational wave background 5 2.5 2.0 1.5 1.0 0.5 0.0 0.5 1.0 1.5 the process can lead to star formation events known as star bursts.
During these episodes, the cold gas present in the galaxies undergoes conversion into stars at a specific mass ratio (Somerville et al. 2001): where  sat and  cen are the total mass (cold gas mass and star mass) of the satellite and central galaxies, respectively.After a star formation activity is triggered, the subsequent stellar evolution is assumed to be completed in an instantaneous way in GABE.Stars evolve to their final states (e.g. age ∼ 10 Gyr) right after their birth.During this process, 43% of the stellar mass is recycled back to the circumstance through stellar winds and supernovae, while metals are produced with a yield of 3% of the initial stellar mass.These newly created metals are thoroughly mixed with the pre-existing cold gas, and spread to other galaxy components through baryonic cycle in galaxies.The resulting metal enrichment of galaxies in GABE exhibits a reasonable agreement with observational mass-metallicity relations (see Fig. 2 and section 2.2.7 of Jiang et al. 2019).Note that the instantaneous assumption is only used to calculate the recycled mass and metals.The detailed evolution of compact binaries and the time delay of each merger have been specifically considered, as described below.
Each star formation event is treated as the birth of a SSP and the corresponding information of the SSP, such as the formation time, mass, and metallicity, is recorded.This method allows us to generate star formation histories for galaxies of various types.Subsequently, the total number of mergers within the simulation volume in the redshift interval Δ is calculated by: where  represents species of merging stellar compact binaries: BNSs, NSBHs or BBHs;   and   are mass and metallicity of the -th SSP respectively, determined self-consistently in GABE;  SSP is the total number of the SSPs;  delay, is -th delay time bin of kind  object of the -th SSP; and  max is the number of delay time bins in -th SSP needed to cover the redshift interval Δ.The equation reveals that Δ  in the redshift interval Δ is contributed by the stellar compact binaries that will merge in this interval (the secondary sum over ) from all SSPs (the first sum over ).See more details in Appendix.Thus, by utilizing the aforementioned approach, we are able to obtain the numbers and rates of merging events for various types of stellar compact binaries.The findings from this analysis are subsequently presented and discussed in the following section.
It is worth noting that, according to the binary population synthesis model COSMIC v3.3.0, the mass of the heaviest stellar black holes is estimated to be around 45 ⊙ , while the lightest massive black holes have a mass of approximately 100 ⊙ with the lightseed model (MSII).It is important to emphasize that no channels have been introduced in this work to generate black holes within the pair-instability mass gap.Nevertheless, merging BBHs have been observed within this mass gap, such as the case of GW190521 (Abbott et al. 2020), which may involve a primary black hole of 85 ⊙ and a secondary black hole of 66 ⊙ .It is suggested that events similar to GW190521 may have a dynamical origin (Romero-Shaw et al. 2020), or even a cosmological origin (Clesse & Garcia-Bellido 2020), but this goes beyond the scope of our current investigation.

Event numbers and merger rates
We present the event number distributions of merging massive black hole binaries in the MSII simulation over chirp mass and redshift, where both light-seed and heavy-seed models are taken into account.The resulting distributions are illustrated in the upper panels of Fig. 2. Our analysis reveals that the redshift-chirp mass distribution of both seed models is quite similar above the chirp mass   = 10 6  ⊙ , whereas it becomes increasingly distinct with decreasing chirp mass.Additionally, we also display the event number distributions of merging stellar compact binaries, which are found to be independent of the seed models, in the lower panels of the same figure.In the lower right panel of Fig. 2, we notice that most massive BBH mergers occur at low redshift ( < 2).These merging BBHs of high mass actually form in metal-poor SSPs ( ≲ 0.004) at high redshift, while they exhibit significant delay times, ranging from ∼ 1 Gyr to the age of the universe.As a result, these high-mass BBHs do not undergo immediate merging after their birth; instead, their mergers occur at low redshifts.
The upper panel of Fig. 3 presents the merger rates of different types of stellar binaries.Specifically, the local merger rates of BNSs and BBHs are consistent with those reported in the third Gravitational-Wave Transient Catalog (GWTC-3) (The LIGO Scientific Collaboration et al. 2021), at 436.8 and 37.0 Gpc −3 yr −1 , respectively.However, the local merger rate of NSBHs is found to be lower, at 3.0Gpc −3 yr −1 , compared to that reported in GWTC-3.Possible sources of this discrepancy include uncertainties in the semi-analytic galaxy formation model and binary population synthesis model, as well as the formation channels of binary compact objects.In particular, the event rate may be influenced by the uncertainties in the SFRD and metallicity evolution history in the semi-analytic models, and by the unclear physical processes of binary evolution in the binary population synthesis model.Additionally, our calculation only considers isolated field binaries, whereas merging stellar compact binaries may also form through other channels, such as dynamical encounters and cosmological processes (e.g.Bavera et al. 2022).
The redshift evolution of merger rates of BBHs, BNSs, and NS-BHs above shows very broad, flat peaks around  ∼ 3.7,  ∼ 2.9, and  ∼ 3.4, spanning a significant redshift range, respectively.These peaks are located at higher redshifts compared to most other studies (e.g., Fig. 6 in Santoliquido et al. 2022), primarily because the SFRD derived by GABE is slightly different from that in other studies.We display the SFRD from GABE and compare it with that from Madau & Fragos (2017), in the lower panel of Fig. 3.We also present various observational estimates compiled by Loiacono et al. (2021), including, the results from UV surveys (Schiminovich et al. 2005;Wyder et al. 2005;Dahlen et al. 2007;Reddy & Steidel 2009;Robotham & Driver 2011;Bouwens et al. 2012a,b;Cucciati et al. 2012;Schenker et al. 2013;Bouwens et al. 2015), the results from infrared, mm, and radio-selected galaxies (Sanders et al. 2003;Takeuchi et al. 2003;Karim et al. 2011;Magnelli et al. 2011Magnelli et al. , 2013;;Gruppioni et al. 2013;Rowan-Robinson et al. 2016;Dunlop et al. 2017;Novak et al. 2017), the result from an optical-NIR observation (Driver et al. 2018), and the result from gamma-ray bursts (Kistler et al. 2009).In addition, we include estimates from the ALPINE collaboration (Gruppioni et al. 2020;Loiacono et al. 2021;Khusanova et al. 2021).Below z ∼ 2, the SFRD from GABE is lower than the median of the observational results by ∼ 0.2 dex (see the dashed line plotted in the lower panel of Fig. 3, which is boosted by this factor.).At higher z, there are large scatters in the observational results of the SFRD, the SFRD from GABE is consistent with some observations (e.g., Kistler et al. 2009;Rowan-Robinson et al. 2016;Gruppioni et al. 2020;Loiacono et al. 2021).The noteworthy feature of the SFRD from GABE is its peak at higher redshifts (z∼ 4), which distinguishes it from the SFRD presented in Madau & Fragos (2017).Notably, some other semi-analytic models also exhibit broad peaks, albeit occurring at lower redshifts (e.g.see Fig. 11 of Henriques et al. 2015).
The lower SFRD at  < 2 in our model certainly underestimates the merger rates of the stellar compact binaries in the redshift range by about 0.2 dex at most.This upper limit is derived assuming an extreme case where the SFRD is underestimated throughout all redshift ranges.We should keep this as a caveat to the results relevant to the stellar compact binaries in the following sections.
Fig. 4 shows the merger rates of massive black hole binaries, The merger rates of massive black hole binaries in the universe calculated using GABE.The total merger rates of the MSI and MSII models are shown in solid lines.Three sub-populations are divided based on the chirp mass   of these binaries: sub1,   < 8.7 × 10 2  ⊙ for light-seed and   < 8.7 × 10 4  ⊙ for heavy-seed model; sub2, 8.7 × 10 2 ≤   < 8.7 × 10 6  ⊙ for light-seed and 8.7 × 10 4 ≤   < 8.7 × 10 6  ⊙ for heavyseed; and sub3,   ≥ 8.7 × 10 6  ⊙ including the light-seed (MSI) model and both the light-and heavyseed (MSII) models, represented by solid lines in different colors.The MSI model exhibits approximately 1.5 − 2 magnitudes lower merger rates than the MSII models, attributed to its 125 times worse particle mass resolution.Note that the two solid lines of the MSII models are exactly the same and overlap with each other, as they actually share the same dark matter halo merger trees and only the embedded seed mass is different.We also divide the mergers into different sub-populations according to their chirp masses: sub1 (  < 8.7 × 10 2  ⊙ for lightseed and   < 8.7 × 10 4  ⊙ for heavy-seed); sub2 (8.7 × 10 2 ≤   < 8.7×10 6  ⊙ for light-seed and 8.7×10 4 ≤   < 8.7×10 6  ⊙ for heavy-seed); and sub3 (  ≥ 8.7 × 10 6  ⊙ ).Agreed with the hierarchical clustering of ΛCDM cosmology, massive black hole binary mergers are mostly contributed by sub1 at the beginning and are then generally dominated by sub2.Though sub3 has the highest chirp mass, its event rate is still the lowest even at  ∼ 0.

The multi-band SGWB
By utilizing the merger rates and GW energy spectra, the SGWB contributed by massive black hole binaries and stellar binaries can be obtained via Eqs.( 8) and ( 9), respectively.The middle panel of Fig. 5 illustrates the SGWB spectra produced from different sources.Moreover, the total SGWB generated from all considered sources is presented in the lower panel of Fig. 5 in the frequency range spanning from nHz to kHz.To investigate the influence of the seed black hole models, we calculate the SGWB with two different seed models.Note that we only present the results from the MSI with the light-seed model, because we only use ultra-low frequency band of the MSI results which are unaffected by seed models.In general, the amplitudes of the principal peaks of the SGWB energy density are within one   4).Middle panel: the SGWB from merging massive black hole binaries with different seed models, as well as merging BNSs, NSBHs, and BBHs.The inspiral phases are plotted separately.Lower panel: the total multi-band SGWB.In addition, the sensitivity curves of aLIGO, ET, LISA, and SKA are shown in both the middle and lower panels.order of magnitude in all two MSII cases under consideration, i.e., Ω GW ∼ 10 −9 to 10 −8 .Furthermore, to illustrate the impact of GW energy spectra on SGWB spectra, we calculate the GW energy spectra for binaries with several typical masses using Eq.( 4), and present them in the upper panel of Fig. 5.
To predict the detectability of the SGWB, the PI sensitivity curves of current and forthcoming detectors are plotted as guidelines in the middle and lower panels of Fig. 5.The threshold SNR is taken as SNR th = 3, and the observation periods for ground-based detectors (i.e., aLIGO and ET) and LISA are assumed to be  = 1 year and  = 4 years, respectively.For the SKA, we adopt the number of pulsars as  p = 200, the rms timing residual  = 50 ns with no red noise, the observing period of  = 10 years, and the average observation cadence of 1 per week (Campeti et al. 2021;Mingarelli et al. 2019).
In the ultra-low frequency band (∼ 10 −9 − 10 −7 Hz) which are mainly contributed by merging supermassive black hole binaries, the SGWB amplitude predicted by the MSII is almost one order of magnitude higher than that of the MSI.This is simply because the inadequate volume of the MSII , which results in poor statistics for extremely large supermassive black holes, leading to noisy predictions in this frequency band (see section 3.1 for more details).We report the dimensionless characteristic strain  yr −1 as 5.1 × 10 −16 , which falls within the same order of magnitude as those predicted by previous theoretical studies and observed by previous experiments2 The sensitivity curve of SKA shows its capability of detecting the SGWB across the ultra-low frequency range from ∼ 10 −9 to 10 −7 Hz, following the typical power law Ω  ∝  2/3 .
In the low frequency band (∼ 10 −4 − 1Hz), the SGWB spectra obtained from the MSII have much flatter shapes and higher amplitudes compared to those of the MSI.These differences simply reflect the fact that a larger number of intermediate-mass massive black holes are not resolved in the MSI.These suggest that the Ω GW signal can be detected by LISA in this frequency band, irrespective of the seed model employed.Furthermore, the middle panel of Fig. 5 indicates that the merging massive black hole binaries contributes substantially to the SGWB in this frequency range.Within the frequency range of approximately 10 −4 to 10 −2 Hz, LISA has the ability to detect the SGWB generated by merging BBHs and BNSs.Interestingly, the SGWB signal predicted by the heavy-seed model exhibits a higher and flatter spectrum than that of the light-seed model in the same frequency band.This is due to the higher event number of massive black holes in the range of ∼ 10 4 − 10 6  ⊙ in the heavy-seed model, as demonstrated in the upper panel of Fig. 2.These sources dominate the GW energy spectra in the LISA band, as shown in the upper panel of Fig. 5, hence enabling LISA to potentially discriminate between black hole seed models.
In the high frequency band (∼ 1 − 10 4 Hz), the SGWB is primarily dominated by the contribution from merging BBHs ( ∼ 1 to 100Hz) and merging BNSs (∼ 10 3 to 10 4 Hz).Specifically, at ∼ 100Hz, the Ω GW contribution from merging BBHs is ∼ 10 and ∼ 100 times larger than that from merging BNSs and NSBHs, respectively.It is worth noting that the contribution from merging NSBHs is negligible, primarily due to their low event rates, as shown in Fig. 3.In the frequency range of ∼ 100 to 10 3 Hz, the total SGWB spectra exhibit small variations, with amplitudes of Ω GW ∼ 2 × 10 −9 .It should be noted that this outcome is dependent on the adopted GW energy spectrum of BNSs presented in Eq. (4).Typically, the energy spectra of BNSs are truncated at the frequency of innermost stable circular orbit (∼ 800Hz) for simplicity, and the post-merger phases of BNSs are not considered (e.g.Abbott et al. 2018).As depicted in the middle panel of Fig. 5, in the frequency range of ground-based detectors, the SGWB spectra are primarily dominated by the inspiral stages of the binary systems, consistent with the findings of previous research (e.g.Abbott et al. 2016Abbott et al. , 2017;;Marassi et al. 2011).In addition, the ET has the capability to detect the SGWB emanating from not only stellar binaries but also massive black hole binaries.The Ω  from merging massive black hole binaries in the light seed scenario displays double peaks.The peaks at ∼ 10 Hz are mostly caused by the first generation of mergers involving massive black holes with masses close to their initial seed mass (e.g.≤ 2 times the light seed mass).These massive black holes do not undergo a substantial accretion phase.The SGWB spectra obtained in this study exhibit similar shapes to those from previous works, as reported in (e.g., Chen et al. 2019;Capurri et al. 2021).We also note, similar to the merger rates in our model, the estimation of the SGWB could also be underestimated by ∼ 0.2 dex at most.

Detection ability
In addition to estimating from PI sensitivity curves, we have calculated the SNR of the SGWB, as defined in Section 2.2, to evaluate the detectability of upcoming detectors.We present the results in Table 1.The obtained SNRs for aLIGO, ET, and SKA in the given observing periods are approximately 5, 700, and 42.84, respectively (note that the obtained SNRs for aLIGO and ET might be underestimated up to ∼ 0.2 dex.In any case, the SNR is much higher than a detection SNR threshold.).However, the SNR from LISA exhibits significant dependence on the adopted black hole seed model, yielding values of approximately 6.7 × 10 4 and 1.5 × 10 5 for the light-seed and heavy-seed models, respectively.This indicates that LISA has the potential to differentiate between black hole seed models.It is of significance to estimate the time required for detectors to collect sufficient data and announce detections, which is equivalent to determining the observing periods necessary for a detector to exceed a given threshold SNR.In this study, we have evaluated the observing periods required for the detectors under consideration by setting the threshold SNR as SNR th = 3.Our findings are presented in Table 2.The results indicate that the SGWB is expected to be detected by aLIGO in a relatively short timescale of a few months, while SKA may require a duration of approximately two years.Notably, regardless of the adopted black hole seed model, the ET and LISA detectors may detect the SGWB immediately after the collection of a sufficient duration of data to calculate the SNR.As a result of the underestimated SGWB in the high frequency band in our model, the above estimation of time for detection by aLIGO and ET is conservative.

SUMMARY AND DISCUSSION
In this study, we explore the potential for the joint detection of the multi-band SGWB from astrophysical sources by future groundbased, space-based detectors and pulsar timing arrays.To this end, we have developed a self-consistent methodology that combines the MSI and MSII simulations with a semi-analytic galaxy formation model-GABE and a binary population synthesis model (COSMIC v3.3.0) to model galaxy and binary evolution in the universe.By treating binary and galaxy evolution in a self-consistent manner, we obtain a comprehensive model of the population of stellar compact binaries (BNSs, NSBHs, BBHs) and massive black hole binaries in the universe.
Our results show that the total multi-band SGWB from merging massive black hole binaries and stellar compact binaries has main peaks with amplitudes of ∼ a few 10 −9 (Fig. 5), which could be easily detected by future GW detection experiments at different sensitivity bands (e.g.ET, LISA and SKA).The shape of the SGWB spectrum in the low frequency band may be influenced by the assumed seed black hole models, which could be determined by LISA.Furthermore, in the high frequency band, the SGWB from merging BBHs is ∼ 10 and ∼ 100 times higher than the SGWB from merging BNSs and merging NSBHs at ∼ 100 Hz.The SGWB is nearly flat in the ∼ 100−1000 Hz frequency range, with an amplitude of ∼ 2 × 10 9 .This result strongly depends on the GW energy spectrum adopted when the frequency is over ∼ 800 Hz.We note that the GW energy spectra of BNSs are closely related to the post-merger physics of neutron stars, and thus the SGWB over ∼ 800 Hz may have the potential to test the equation of state of BNSs.
Note that the SFRD derived from GABE is slightly lower than the observational results by ∼ 0.2 dex at low redshifts (z<∼ 2).Consequently, our prediction of the SGWB produced by stellar compact binaries might be underestimated by a similar factor at most.Given uncertainties in other effects, like the time-delay function, metallicity in the binary population synthesis model, our results on high frequency band (∼ [1 , 10 4 ] Hz) of the SGWB are still conservative.

Figure 1 .
Figure 1.The delay time distribution  SSP of merging BNSs, NSBHs, and BBHs is shown from left to right for seven different SSPs with varying metallicities.The age of the SSPs is presented along the x axis.

Figure 2 .
Figure 2. Upper panel: the event number distribution  of merging massive black hole binaries over chirp mass and redshift, computed using the MSII simulation in two scenarios of light-seed and heavy-seed models.The values of event numbers are indicated by color bars in the unit of log 10 (d 2  /(ddlog 10   ) ). Lower panel: the event number  distribution of merging stellar compact binaries over chirp mass and redshift.The values of event numbers are indicated by color bars in the unit of log 10 (d 2  /(dd  ) ).

Figure 3 .
Figure 3. Upper panel: The merger rates of BNSs, NSBHs and BBHs in the universe calculated from GABE.The blue, red and green bars are the estimation range of local merger rates (The LIGO Scientific Collaboration et al. 2021) according to GWTC-3.Lower panel: The SFRD in GABE, together with the SFRD inMadau & Fragos (2017).The empirical constraints for the SFRD available from the literature based on multiwavelength observations are also shown.All the data points and SFRD values correspond to a Chabrier IMF(Chabrier 2003).The peaks of the merger rates and GABEderived SFRD are marked by the intersection of two dotted lines.

Figure 5 .
Figure5.Upper panel: the GW energy spectra of various merging systems, including BNSs, NSBHs, BBHs, and massive black hole binaries.The spectra are calculated using Eq.(4).Middle panel: the SGWB from merging massive black hole binaries with different seed models, as well as merging BNSs, NSBHs, and BBHs.The inspiral phases are plotted separately.Lower panel: the total multi-band SGWB.In addition, the sensitivity curves of aLIGO, ET, LISA, and SKA are shown in both the middle and lower panels.

Figure A1 .
Figure A1.Schematic illustration of Eq. (21).The two arrays shown in blue represent two SSPs formed at  0 respectively, while the green array represents Δ  () within a redshift interval Δ.The thick lines and scales below each array depict the flow of cosmological time.

Table 1 .
The SNR for aLIGO, ET, LISA and SKA with light-seed (MSI), light-seed (MSII) and heavy-seed (MSII) model.The observing periods are assumed as 1, 4 and 10 year(s) for ground-based detectors, LISA and SKA, respectively.For SKA we adopt  p = 200 pulsars, the rms timing residual  = 50 ns with no red noise, and the average observation cadence of 1 per week.

Table 2 .
The observing time required for aLIGO, ET, LISA, and SKA to achieve a threshold SNR of 3, with consideration of the light-seed (MSI), light-seed (MSII), and heavy-seed (MSII) models.