Mass-redshift dependency of Supermassive Black Hole Binaries for the Gravitational Wave Background

Studying how the black hole (BH) - (galaxy) bulge mass relation evolves with redshift provides valuable insights into the co-evolution of supermassive black holes and their host galaxies. However, obtaining accurate measurement of BH masses is challenging due to the bias towards the most massive and luminous galaxies. Instead we focus on the BH and bulge masses as they vary with redshift using the EAGLE, Illustris, TNG100, TNG300, Horizon-AGN and SIMBA large-scale cosmological simulations. We use an analytical astrophysical model with galaxy stellar mass function, pair fraction, merger timescale and BH-bulge mass relation extended to include redshift evolution. The model can predict the intensity of the gravitational wave background (GWB) produced by a population of supermassive black hole binary (SMBHB) as a function of the frequency. This allows us to compare the predictions of this model with the constraints of Pulsar Timing Array observations. Here, we employ Bayesian analysis for the parameter inference. We find that all six simulations are consistent $\leq 3.5\sigma$ with a range of simulated GWB spectra. By fixing the BH-bulge mass parameters to the simulations we analyze the changes in the constraints on the other astrophysical parameters. Furthermore, we also examine the variation in SMBHB merger rate with mass and redshift between these large-scale simulations.


INTRODUCTION
The co-evolution of galaxies and their supermassive black holes (SMBHs), i.e. the relationship between SMBHs and the dark matter halo potential, their role in the stellar formation activity, their local interactions with the stars and gas, and their fate during the history of galaxy mergers, are key ingredients of recent large cosmological simulations and of our understanding of large-scale structure formation and evolution (see e.g., Habouzit et al. 2021, 2022a, andreferences therein).
A PTA uses radio telescopes to time a network of millisecond pulsars (Sazhin 1978;Detweiler 1979).In principle, once the pulsar rotation irregularities, its possible orbital motion, the dispersion and scattering of its radio signal through the interstellar and heliospheric plasma and the systematics due to the Earth's motion in the Solar System are properly modelled and subtracted from the time series of measured pulsations, one expects to be able to extract the GW imprint from the resulting timing residuals.The analysis requires observations of multiple millisecond pulsars at sub µs precision for several decades (up to about 25 years for ongoing programs) in order to extract a GWB from unmodelled noise.There are several PTA consortia, structured at continental levels and collaborating globally: European PTA (EPTA) (Kramer & Champion 2013;Desvignes et al. 2016;Chen et al. 2021), Parkes PTA (PPTA) in Australia (Manchester et al. 2013;Hobbs 2013;Kerr et al. 2020), North American Nanohertz Observatory for Gravitational Waves (NANOGrav) (Arzoumanian et al. 2016(Arzoumanian et al. , 2018(Arzoumanian et al. , 2020)), Indian PTA (InPTA) (Joshi et al. 2018;Tarafdar et al. 2022), Chinese PTA (CPTA) (Lee 2016;Jiang et al. 2019) and MeerTime PTA (MPTA) in South Africa (Bailes et al. 2020;Spiewak et al. 2022).These PTAs form a world wide organisation, the International PTA (IPTA), where they share their data and coordinate their analysis to eventually detect and hopefully characterise the GW signal (Hobbs et al. 2010;Verbiest et al. 2016;Antoniadis et al. 2022).
NANOGrav, PPTA, EPTA and IPTA have reported the detection of a low-frequency common signal in their pulsar datasets (Arzoumanian et al. 2020;Goncharov et al. 2021;Chen et al. 2021;Antoniadis et al. 2022).This marks the first step towards the detection of a GWB.If the common signal is of gravitational wave origin it should also show a characteristic spatial correlation between the pulsars, called the Hellings-Downs correlation (Hellings & Downs 1983), which the above mentioned collaborations and the Indian Pulsar Timing Array have found evidence for (EPTA Collaboration et al. 2023;Agazie et al. 2023a;Reardon et al. 2023;The International Pulsar Timing Array Collaboration et al. 2023).In addition, the Chinese Pulsar Timing Array concurrently also found significant evidence for a Hellings-Downs correlated signal in their dataset (Xu et al. 2023).
If these recently observed spectral signatures are from a population of supermassive black hole binaries (SMBHBs), they favour heavy black hole masses and short merger timescales.Future detections will improve on these constraints and should allow some relations to be ruled out, in particular those with the lowest GWB.This would open new multi-messenger probes to study SMBHs and their host galaxies (e.g., Pol et al. 2021).
By formulating the relative strength of the GWB as a function of SMBHB merger rate and gravitational wave energy spectrum, we can connect them to astrophysical parameters.The SMBHB merger rate is linked to the galaxy merger rate via a mass relation between the SMBH and galaxy bulge.Using the Galaxy Stellar Mass Function (GSMF), a differential pair fraction of galaxy in binaries and a merger timescale one can compute the galaxy merger rate.The gravitational wave energy spectrum depends on the binary orbital eccentricity and the nature of the environment driving their evolution (Chen et al. 2017b(Chen et al. , 2019)).
The mass relation between the SMBH and galaxy bulge, called the BH-bulge mass relation, is widely studied using both observational data and large-scale cosmological simulations.The different values of the BH-bulge mass parameters for our Universe are constrained using observational data.Although there is currently no consensus, several observational samples suggest that the BH-bulge mass relation could evolve with redshift (Merloni et al. 2010;Kormendy & Ho 2013).In these papers, for a fixed galaxy mass, BHs are on average more massive at high redshift compared to those in similar host galaxies at low redshift.
Studying the evolution of the Universe through observations is a challenging task due to a number of technical limitations.The expansion of the Universe causes the light from the galaxies and SMBHs to shift towards longer wavelengths, making it difficult to detect their emission and accurately measure their properties, such as their mass and accretion rate.For example, it can be difficult to study scaling relations at high redshifts beyond z ∼ 2 due to the challenges of disentangling the light from an active galactic nucleus (AGN) and the light from the host galaxy (Ding et al. 2020).The high redshift galaxies are fainter and smaller than nearby galaxies, which makes it challenging to study their structure and dynamics (Kormendy & Ho 2013).It is important to consider the types of systems that are selected for observation, as this can introduce biases, such as a focus on galaxies with AGNs, which are not representative of the overall galaxy population.These technical limitations can make it difficult to obtain detailed and accurate interpretation of the BH-bulge mass relation.
Large-scale cosmological simulations have been successful in reproducing many aspects of the Universe with a high degree of accuracy.One aspect that has been well reproduced is the largescale structure of the Universe, including the distribution and size of galaxies, clusters of galaxies, and cosmic voids (Genel et al. 2018;Pillepich et al. 2018b).These simulations have also been successful in reproducing the observed distribution of matter in the Universe, including the distribution of dark matter, which is difficult to detect directly (Vogelsberger et al. 2020;Angulo & Hahn 2022, and references therein).In particular, we use: EAGLE (Schaye et al. 2015;Crain et al. 2015), Illustris (Genel et al. 2014;Vogelsberger et al. 2014), TNG100, TNG300 (Springel et al. 2018;Naiman et al. 2018;Marinacci et al. 2018;Nelson et al. 2018;Pillepich et al. 2018a,b), Horizon-AGN (Dubois et al. 2014(Dubois et al. , 2016)), and SIMBA (Davé et al. 2019).
Our aim in this work is to setup the methodology to constrain the SMBHB properties using future PTA observations.We concentrate on the BH-bulge mass relation and test for its redshift dependence.Existing formulations of the BH-bulge mass relation as a function of redshift for z < 5 can be improved in light of ,e.g., the recent developments in cosmological simulations and observations from new instruments.Thus, we formulate an equation for the BH-bulge mass relation taking into account the redshift of the system and apply this equation to fit for BH and galaxy stellar mass data from several large-scale cosmological simulations.This BH-bulge mass relation with redshift dependence is then used in an analytical astrophysical model to compute the intensity of the GWB generated by a population of SMBHBs focusing on the PTA frequency range.Bayesian analysis is used to find the posterior of all the parameters of this GWB model.We also fix the BH-bulge mass parameters to those fitted to the cosmological simulations to constrain the posteriors of other parameters.
The paper is organized as follows.Section 2 describes the astrophysical model to compute the GWB formed by the mergers of a population of SMBHBs in a parametric form using the GSMF, pair fraction and merger timescale.Section 3 focusses on the relation between the galaxy bulge and central black hole mass, where we review the redshift independent relation and extend it by fitting to results from large-scale cosmological simulations.In Section 4, the analysis setup, the priors motivated by observations and large-scale cosmological simulations, and the simulation of GWB detections with different strains are described.We present our results in Section 5 for the different GWB strains and also study the impact of using fixed BH-bulge mass parameters fitted to cosmological simulations.Finally, Section 6 outlines the conclusions.

GWB CHARACTERISTIC STRAIN
For a population of SMBHBs the characteristic spectrum of the GWB was expressed in Phinney (2001) as where f is the frequency, G is the Newton's constant, c is the speed of light, and z is the redshift.The chirp mass M is given as where M 1 , M 2 are the individual SMBH masses in the binary system.Chen et al. (2019), which is extended by a parameter describing the redshift dependent evolution of the BH-bulge relation, see Section 3.

Analytic model and fitting function
Using the formalism of Chen et al. (2017b) we write dE d f r in terms of sum of harmonics at each eccentricity e n at each orbital frequency of the binary as where F(e) = 1 + (73/24)e 2 + (37/96)e 4 (1 − e 2 ) 7/2 , (4) and J n is the first kind of n th Bessel function.
To increase the computational efficiency Chen et al. (2017b) use the characteristic strain spectrum h c,0 ( f ) of a reference SMBHB with e 0 = 0.9 at f 0 = 10 −10 and peak frequency f p,0 .For a generic SMBHB with e t at f t ̸ = f 0 the strain can be computed as with the peak frequency A trial analytic function for the characteristic spectrum for the reference SMBHB with f = f /(10 −8 Hz) can be written as The constants a 0 , a 1 , a 2 , b 0 , b 1 , b 2 , c 0 , c 1 , c 2 are determined by the fit and are given in Chen et al. (2017b).By considering SMBHBs with different redshifts and chirp masses, we get the characteristic spectrum of a population of SMBHBs as

Stellar environment
The GWB energy spectral shape is affected by the environmental coupling.A super-efficient inspiral can cause a bend in the GWB spectrum in the PTA frequency range (Sesana 2013a;Ravi et al. 2014;Huerta et al. 2015;Chen et al. 2017b).At short separations the gravitational radiation starts to dominate the binary evolution, after a phase where the energy loss was driven by interactions with stellar or gaseus environment (Sampson et al. 2015).We now consider that in our model stellar hardening dominates at low frequency until it is overtaken by the GW emission at the transition frequency where the chirp mass M 9 = M /(10 9 ) is rescaled, ρ i,100 is the stellar density of the environment within the SMBHB influence radius, the additional multiplicative factor ρ 0 includes all systematic uncertainties while estimating ρ i,100 and σ 200 is the stellar central velocity dispersion in the galaxy, which are given by the stellar density distribution's inner slope is given by γ ∈ (0.5, 2), a is the characteristic radius and M is the total bulge mass of the galaxy, which are expressed as

Merger rate
The merger rate in Eqn (1) can be written in terms of SMBHB mass as M BH can be parameterized using galaxy bulge mass as shown in Section 3.An astrophysical observable based description of the galaxy merger rate is given in Sesana (2013a); Sesana et al. (2016); Chen et al. (2019) as where q is the galaxy binary mass ratio with the primary galaxy mass M, Φ(M, z) = (dn G /d log M) z is the GSMF estimated at redshift z, the differential pair fraction of the galaxy binaries is F (z, M, q) = (d f /dq) z,M and the merger timescale τ(z, M, q) = z z ′ (dt/dz)dz is obtained by integrating over the instantaneous redshift dz between the redshifts at the start z ′ and end z of the galaxy merger.Using a flat ΛCDM model one finds Here we use energy density ratios Ω M = 0.3, Ω k = 0, Ω Λ = 0.7 and Hubble constant H 0 = 70 kmMpc −1 s −1 .

Galaxy Stellar Mass Function
The Galaxy Stellar Mass Function (GSMF) describes the number density of galaxies as a function of their stellar mass.The assembly of stellar mass and the evolution of the stellar formation rate through cosmic time can be traced using the GSMF and is a major estimate of the characteristics of the galaxy population.
This astrophysical observable can be parameterized and fitted in the form of a Schechter function (Conselice et al. 2016).To take into account redshift evolution we can write the GSMF using parameters from Mortlock et al. (2015) as

Pair fraction
The differential pair fraction of the galaxy binaries at M and z with respect to q can be written as (Mundy et al. 2017) with f 0 = f ′ 0 q γ f dq.Integrating over q then gives

Merger timescale
The timescale of the evolution of a binary galaxy from the dynamical friction can be used to approximate the full merger timescale, which can be written using a parameterisation with τ 0 , α τ , β τ , γ τ as where h 0 = 0.7 is the Hubble parameter.Substituting these observables into Eqn (18) gives

ASTROPHYSICS OF SMBH MASS
The final ingredient to describe the SMBHB merger rate in Eqn ( 16) is a relation between the galaxy stellar mass and the central black hole mass.We first express the bulge mass of a galaxy using its total stellar mass, and then use the resulting BH -bulge mass relation to extract the BH mass needed for the computation of the merger rate.
The fraction of the total stellar mass assigned to the bulge mass depends on the galaxy morphology and galaxy mass regime.For this work the phenomenological stellar-bulge mass relation (Bernardi et al. (2014); Sesana et al. (2016)) is used (25) This relation focusses on spherical and elliptical galaxies, which dominate the PTA GWB signal.Higher mass galaxies M ≤ 10 10 M ⊙ have been observed to be correlated with the size of the bulge and disk, while lower mass galaxies do not.

Large-scale cosmological simulations
In this paper we investigate the differences in the BH-bulge mass relations produced in EAGLE, Illustris, TNG100, Horizon-AGN, SIMBA, and TNG300, and quantify the evolution of the relation with redshift.The galaxy stellar mass and the corresponding SMBH mass from the simulations are given in Habouzit et al. (2021).The conversion of the stellar mass of the galaxies into their bulge mass is done using Eqn (25).Thus, the BH-bulge mass relation is connected to the BH-galaxy stellar mass relation.
Cosmological simulations model the dark matter and baryonic contents of the Universe in an expanding space-time.All the simulations studied in this paper have a volume of ⩾ 100 3 cMpc 3 , a dark matter mass resolution of ∼ 5 × 10 6 − 8 × 10 7 M ⊙ , and a spatial resolution of 1 − 2 ckpc.As such, the simulations capture the time evolution of the galaxies with a total stellar mass in the range M = 10 9 − 10 11−12 M ⊙ and their BHs.Baryonic processes taking place at small scales below the galactic scale are modelled as subgrid physics (e.g., supernova (SN) and AGN feedback).Although theoretically based on the same idea, these processes are modelled differently in each simulation.For example, AGN feedback releases energy in the BH surrounding but the implementation in the simulation can rely on the injection of thermal energy only, thermal and kinetic energy or momentum in a given direction to mimic an outflow or jet.The subgrid physics of the simulations impact the evolution of both galaxies and BHs (Habouzit et al. 2021(Habouzit et al. , 2022a)).
There is no consensus on the shape nor on the time evolution of the BH-bulge relation produced by the EAGLE, Illustris, TNG100, Horizon-AGN, SIMBA, and TNG300 simulations (Habouzit et al. 2021(Habouzit et al. , 2022b)).The shape of the BH-bulge relation in the low-mass end (M ⩽ 10 10.5 M ⊙ ) is mainly driven by BH seeding mass, strength of SN feedback and BH accretion modelling.The massive end is affected by the modelling of AGN feedback and BH accretion.Half of the simulations have more massive BHs at high redshift than at z = 0 at fixed galaxy stellar mass.The other simulations follow the opposite trend.On average, the time evolution of the relation depends on whether BHs grow more efficiently than their host galaxies (see summary in Fig. 11 in Habouzit et al. (2021)).

Empirical BH-bulge mass relation
The BH-bulge mass relation is a key quantity for our understanding of the co-evolution of galaxies and their central black holes.The redshift independent BH-bulge mass relation (Kormendy & Ho 2013) which is usually used in the literature is given by where M BH is the mass of the SMBH at the centre of the galaxy with bulge mass M bulge and N denotes a Gaussian distribution.α and β are the BH-bulge mass parameters that determine the slope and normalization of the relation respectively.On the logarithmic scale the relation becomes a straight line with scattering ε where the parameters can be deduced from a least squares fit.Reviews on different models and parameters of this relation can be found in Sesana (2013b) and Schutte et al. (2019).

BH-bulge mass relation with redshift dependence
Our goal is to formulate a parametric redshift dependent BH-bulge mass relation for z ≤ 5 since the GWB should be detectable with PTA in this redshift range.We propose a relation given by where we consider an additional BH-bulge mass parameter γ * , which determines the extent of the evolution of the black hole mass with redshift.Positive γ * values result in larger black hole masses as the redshift increases, while negative γ * values have the opposite effect.This relation is based on the assumption that the BH-bulge mass relation evolves only through the normalization parameter β * , while the slope α * remains constant with redshift.This assumption is based on the observation that the correlation between the mass of the black hole and the mass of the bulge is largely set by the processes that lead to the formation of the bulge, which happen early in the galaxy's history.These processes are not expected to change significantly over cosmic time and so the slope parameter is expected to remain relatively constant.However, the normalization parameter is expected to change with redshift because the growth of the black hole and the bulge are linked through complex feedback processes.These feedback processes are expected to change over time as the galaxy evolves, and so the normalization parameter is expected to evolve with redshift (Kormendy & Ho 2013).
Eqn ( 29) is fitted to the SMBH and galactic bulge masses for each of the six cosmological simulations separately.For a given simulation α * is the slope in the logarithmic scale given by the linear least squares fit over all redshifts z ≤ 5. β * is the intercept at M bulge = 10 11 M ⊙ of the least squares fit at z = 0.The intercepts at different redshifts is used to compute γ * .The amount of scattering of the SMBH mass from the phenomenological fit in Eqn ( 29) is denoted by ε.Table 1 lists the BH-bulge mass relation parameters for these cosmological simulations.The variation of the masses for these simulations are plotted in Figure 1 as they evolve with redshift.Figure 2 Fitting the black hole and bulge masses to this relation, we obtain large variability of the parameters and higher scattering values for most of the simulations we have used.Thus, in this work we use the previous relation that produces stable parameter values at different redshifts and lower scattering values consistent with the simulations.
We note that our relation (Eqn ( 29)) is a simple approximation to the simulations extending the redshift independent BH-bulge relation.Therefore, our best fit values may not be fully representative of the results from the simulations.This caveat should be kept in mind with the results presented in Section 5.
In order to find the redshift volume that PTA can probe for galaxy and SMBH mergers we consider z m , which is the maximum redshift that is used to compute the volume, as an additional parameter to see the change in GWB characteristic strain if the volume is larger or smaller.
The effect of each of the 20 astrophysical parameters on the GWB is shown in Figure A1 for a fiducial choice of values.Using the corresponding values from Table 1 for the six large-scale simulations the differences in the GWB spectra can be seen in Figure 3.
With this parametric model in hand we can set up the Bayesian analysis to use simulated PTA detections to infer what posterior constraints can be achieved for each parameter.Figure 3. GWB characteristic strain spectra in the PTA range from a set of fiducial values (see Figure A1) showing the differences when using the BHbulge mass parameters for the EAGLE, Illustris, TNG100, Horizon-AGN, SIMBA, and TNG300 simulations respectively.An analytic sensitivity curve from the IPTA DR2 (Antoniadis et al. 2022) is plotted to guide the eye.

Simulated GWB detections
Different values for the 20 parameters within the prior ranges give different GWB characteristic strain.Depending on the values of the 20 parameters, we can simulate a straight line or a curve bending down at low frequency for the GWB characteristic strain in the frequency range of 10 −9 − 10 −6 .In our model the straight line and curve spectra are associated with circular and eccentric SMBHB populations respectively.We created datasets for these two different shapes of spectrum for strain values of 0.5 , 1 , 2 , 3 , 4 × 10 −15 at the reference frequency of f = 1/1year ( f ≈ 10 −7.5 Hz) as shown in Figure 4. PTAs typically search at frequencies that are multiples of 1/T span , where T span is the total observation time span of the PTA dataset.For simplicity and computational efficiency, we use the five lowest bins with T span = 25 years.The values of the parameters used to create the different simulated spectra are chosen by hand and given in Table B1.These sets of parameters are non-unique and thus not necessarily representative for the given GWB spectrum.

Likelihood function
To (32) The Parallel Tempering Markov Chain Monte-Carlo (PTM-CMC) Sampler (Ellis & van Haasteren 2017) is used with log 10 A det ( f ) taken from the simulated GWB datasets and σ det = 0.09 for the Bayesian analysis.

Prior choice
The prior for the BH-Bulge mass relation is constrained using all possible masses from the six different simulations for z ≤ 5 to set the allowed range as shown in Figure 5 and the initial test values are given the Table 2.Only combinations of α * , β * ,γ * and z = (0, z max ) that give relations within the boundaries are accepted.This ensures that the BH-bulge relations are compatible with those from the simulations for all redshifts between 0 and z max .
It is assumed that the redshift volume that PTAs are sensitive to is between 1.5 and 2.5.To study the effects of evolution with redshift and to test this assumption, we consider an extended range of z m ∈ [0.1, 5], which includes the above range.
For the other parameters we adopt the same prior choice as Chen et al. (2019) shown in Table 3.

Consistency of the cosmological simulations with GWB detections
Simulated PTA GWB detections are used to perform the Bayesian analysis to find the posterior constraints on the astrophysical parameters in our model.We first investigate the consistency of the fitted values that are an approximate representation of the complex simulations with the different shapes and strains of the simulated PTA detections.
Figure 6 shows the p-values from Kolmogorov-Smirnov (KS) tests on whether the parameters from the six simulations can be consistent with being drawn from the underlying posterior distributions.Since we have four parameters in the redshift dependent BH-bulge relation, each parameter is investigated independently.In general, for most of the simulations and simulated detections, the p-value are well above 0.1, indicating that the fitted values are possible draws from the posterior distributions.For α * , γ * and ε the p-values do not vary much for each simulations across the different strains and shapes of the GWB spectrum.The main changes can be seen for β * , this could be due to dominant role β * plays in determining the overall GWB strain level.We can very broadly see two trends in the p-values: 1. where they tend to grow as the GWB strain increases and 2. where they behave in the opposite way.Looking at Figure 6 the simulations can be separated by the two trends into two groups: 1. TNG100, TNG300 and SIMBA, following the first trend and 2. EAGLE, Illustris and HorizonAGN, which behave by the second trend.
As the KS tests are performed on marginalized 1D distributions and do not take covariances into account, we also employ the Mahalanobis distance to give another quantity for the consistency between a simulation and a simulated PTA detection.All simulations give distances between about 1 to 3.5 for all GWB strains and spectral shape, see Figure 7.The same two groups of simulations can be found to follow the same trend, where the first (TNG100, TNG300 and SIMBA) have decreasing distances and the second (EAGLE, Illustris and HorizonAGN) become less consistent.In general, the first look to be more consistent with simulated PTA detection than the second group.

BH-bulge mass relation
Looking more closely at the posterior constraints on the parameters of the model we note that most of them are very similar to their priors, indicating that they are either already well constrained by other observations or they play only a mild role in the amplitude of the strain values.One of the two main constrained observables is the merger time.It depends on the strength of the GWB, where a higher amplitude leads to shorter merger times and a lower amplitude allows for longer merger times.
The other constrained observable is the black hole -galaxy bulge mass redshift dependent relation, which is why we focus on the parameters α * , β * , ε and z m in the following.Figure 8 and Figure 9 show their 2D and 1D posterior distributions for the cases of circular and eccentric populations respectively.We show only the cases for the smallest and largest amplitudes from our simulated detections.
First, looking at circular population in Figure 8, there is little difference between the posteriors (black) and the priors (green) (described in Section 4. holes tend to be heavier, a trend for faraway SMBHBs also starts to emerge.
Introducing a bend at the lowest frequencies from eccentric population of SMBHBs, shown in Figure 9, only marginally changes the findings from circular populations.The eccentricity and the environment of the black holes can only have an effect, if the bend is more prominent in the PTA frequency band.
The evolution of the parameter constraints with amplitude can be found in Figure 10.As the characteristic strain amplitude becomes higher most parameters α * , β * , ε and z m prefer higher values and γ * becomes closer to zero.This suggests that a PTA detection can put constraints on the redshift evolution of the BHbulge mass relation.
Given the PTA detections of a common signal of amplitude ∼ 2.5 × 10 −15 and the recent evidence for the GW origin, the constraints on our model will be between the two closest matching simulated detections at 2 and 3 × 10 −15 .However, we use a 25 year observation time span, compared to the ∼ 15 years of the most recent PTA data sets.Our model also takes into account for the possibility that the BH-bulge relation could evolve with redshift and samples for the maximum redshift, which is equivalent to the volume of space, that PTAs can constrain.Extensive astrophysical interpretation was performed by EPTA (Antoniadis et al. 2023) and NANOGrav (Agazie et al. 2023b), which are consistent with our findings.

SMBHB merger rate
An interesting quantity that can be computed from our model is the merger rate of the SMBHBs from Eqn (16) given the constraints on the parameters from the simulated GWB detections.Following Chen et al. (2019), we first integrate over the mass ratio, leaving a merger rate by redshift and chirp mass.Next, we can integrate over the mass and redshift to get Figure 11 and Figure 12 showing dn/dM and dn/dz respectively.
The merger rates with respect to the SMBHB mass in Figure 11 are very similar between the circular and eccentric populations at most investigated strain amplitudes.This indicates that environmental effects are not strongly covariant with the population properties.Only at the lowest amplitude h = 0.5 × 10 −15 differences become noticeable with the eccentric population having a larger number of low mass binaries and the high mass drop-off at lower masses, compared to the circular population.This general trend persists through increasing amplitudes, but becomes less significant.In general, with larger amplitude the rate of massive binary mergers also grows.The median merger rate moves towards a drop-off at higher mass.Additionally, one can see an increase of the merger rate for smaller mass binaries in the 2-sigma range, especially in the circular population.
As we introduced the maximum redshift as a free parameter, the merger rates with respect to the redshift drop to zero at different maximum redshifts.This mimics the expectation that the GWB that PTAs are sensitive to will be dominated by closeby binaries.As such, we have binned the posterior samples by their maxmimum redshift.Within each bin we plot the merger rate within a common range of redshifts in Figure 12.E.g., in the left most column we selected all the posterior samples with a maximum redshift of z < 1 and plot the integrated merger rate between 0.1 and 0.5.
In general, the median merger rate as a function of redshift is nearly constant across most redshifts, amplitudes and different populations.A small raise as the detected amplitude increases can be seen.The difference between the two populations is very small with the eccentric population requiring an overall larger number of mergers.As in the mass dependent merger rates, the main differences can be seen at the lowest amplitude.At h = 0.5 × 10 −15 the drop of the lower bounds of the merger rate at high redshift is clearly visible.This is consistent with the prior assumption of a possible decreasing number of SMBHBs at high redshifts contributing to the GWB.Consequently, the number of samples for the highest maximum redshift is also low.If a detection favours high amplitudes, more binaries even at large distances are required to produce the GWB.This can be seen most prominently in the rightmost column in Figure 12, where the 2-sigma lower bound drop of the merger rate moves from z ≈ 1 to z ≈ 4.

Constraints from simulations
By fixing the BH-bulge mass parameters to the best fit values from the simulations given in Table 1, we can see how the constrains on the other parameters are affected.For computational cost reasons we only analyze the h c = 0.5 , 2 , 4 × 10 −15 detections for both straight line and curved GWB spectra.Figure 13 shows the median values and central 68% of the 1D marginalized posterior distributions for all six simulations.
In general, most parameters have similar posterior compared to the prior constraints (in light green).The five parameters related to the GSMF (Φ 0 , Φ I , M 0 , α 0 , α I ) are already well constrained from observations.Parameters that play only a subdominant role, like those for the pair fraction ( f 0 , α f , β f , γ f ) and maximum redshift z m are only slightly constrained towards larger values for stronger GWB strains.The eccentricity e 0 and stellar density ρ 0 parameters are degenerate.However, we can see that straight line spectra result in low eccentric binaries in low stellar dense environments, while a curved spectrum indicates the need for either eccentricity or dense stellar environments of the binaries.Lastly, the most important observable when using a fixed BH-bulge relation is the merger time.A short merger time is needed to produce a stronger GWB, especially if the masses of the SMBHs are fixed to results from cosmological simulations.The (second column, third row) panel in Figure 13 on the merger time norm τ 0 shows for all six simulations this decrease of the median values as well as the shrinkage of central 68% credible regions.The other three parameters describing the merger time (α τ , β τ , γ τ ) play a minor role and are thus not much more constrained compared to the prior.The corner plots for the 16 parameters with amplitudes h c = 0.5×10 −15 , 2×10 −15 , and 4×10 −15 for both circular and eccentric population of SMBHBs using the fitted BH-bulge mass parameters from the simulations can be found in the online supplementary material.

Parameter constraints from Illustris and SIMBA
Below, we focus on constraints from the Illustris and SIMBA simulations.These two cosmological simulations are chosen since Illustris shows positive while SIMBA shows negative evolution of the SMBH mass with redshift, and thus are the extreme two cases in these six large-scale cosmological simulations.The distribution of all the astrophysical parameters for both curved and straight line characteristic spectra of the GWB and with fixed BH-bulge mass parameter values to match Illustris and SIMBA are given in the Figure 14, where Illustris is shown in orange, SIMBA in purple and the prior of the parameters in the green shaded regions.
The evolution of the pair fraction parameters displays similarities in both the curved and straight line characteristic spectra for different amplitudes in Illustris and SIMBA, with two noticeable differences: 1. f 0 at h c = 0.5 × 10 −15 and 2. β f at h c = 4 × 10 −15 .At a characteristic spectrum value of h c = 0.5 × 10 −15 , both Illustris and SIMBA exhibit posteriors similar to the prior for the pair fraction norm f 0 .However, at this strain amplitude, while Illustris trends towards larger values, SIMBA behaves in the opposite way.As the strain value increases, both Illustris and SIMBA start to display posterior distributions that prefer larger values of f 0 .The pair fraction mass slope α f for both Illustris and SIMBA shows a preference for low values at h c = 0.5 × 10 −15 , followed by no preference at h c = 2 × 10 −15 , and then higher values at h c = 4 × 10 −15 for both the curved and straight line spectra.The posterior distributions of the pair fraction redshift slope β f exhibit similar evolution with amplitude as in the case of α f for both simulations, except at the largest strain value, where the trend is more pronounced in Illustris compared to SIMBA.In conclusion, the pair fraction increases with larger amplitudes for both circular and eccentric populations.More massive and distant galaxy pairs are required to produce the gravitational wave background at higher strains.Illustris tends to require more pairs than SIMBA for the same amplitude.The curved characteristic spectra with values of h c = 0.5 , 2 × 10 −15 in SIMBA reveal a correlation between higher posterior values of α τ and higher values of eccentricity e 0 and ρ 0 .This is in contrast to the general behaviour that increasing characteristic spectrum values lead to lower values of α τ , β τ and τ 0 .The posteriors change from being in broad agreement with the priors at h c = 0.5 × 10 −15 for both curved and straight line spectra and both simulations to trending very clearly towards lower merger times at h c = 4 × 10 −15 .The main difference between Illustris and SIMBA seems to be that the merger time redshift slope β τ is more constrained for Illustris, whereas it is the merger time norm τ 0 for SIMBA.We can see that the curved and straight line spectra at the same amplitude mostly impact the eccentricity e 0 and stellar density ρ 0 parameters.The straight line spectra lead to almost no constraints at all strains.On the other hand, the curved line spectra show the correlation between these two parameters in creating a bend at low GW frequencies.It should be noted that the posteriors look to be less well constrained for e 0 and ρ 0 with larger amplitudes in the curved spectra case.This could be from the difficulties of PTA detections to accurately measure a bend in the GWB spectrum, especially for our simulated detections with limited frequency coverage.
Finally, the inclusion of the maximum redshift z m parameter allows to gauge where the most dominant SMBHBs can be found for a simulated PTA GWB detection and a chosen cosmological simulation.The last row in Figure 14 shows that in general larger strain values require binaries to be concentrated at higher redshifts.For both the curved and straight line spectra Illustris constraints more strongly to large maximum redshifts, while SIMBA only shows a weak trend in the same direction.This could be the effect of the γ * parameter that describes the BH-bulge relation, where positive values, like in Illustris, produce more massive BHs at higher redshifts.Whereas the negative value in SIMBA leads to the most massive BHs being at smaller redshifts.

Mergerrate constraints of Illustris and SIMBA
It is interesting to look at dn/dM and dn/dz as in the previous section for the Illustris and SIMBA simulations as shown in Figure 15 and Figure 16.The most prominent feature of dn/dM in Figure 15 is a shift towards larger mass from Illustris to SIMBA for the same GWB strain.This is consistent with the prediction that SIMBA produces lower mass binaries than Illustris and thus need more binaries to match the emitted GWB strain.The free parameters are adjusted to get same amplitude as explained above from Figure 14.The other feature that is visible from Figure 15 is the variation between circular and eccentric binaries producing a straight and curved GWB strain spectrum respectively.There is no difference in the median of the merger rates with respect to the SMBHB mass in Illustris for the circular and eccentric binaries with same amplitude, however we can see small differences at the lower 2-sigma boundaries.SIMBA clearly shows a slight variation in binary chirp mass between circular and eccentric binaries.
The merger rate with respect to the redshift is shown in Figure 16 for Illustris and SIMBA with the panels defined in the same way as in the previous section and Figure 12.An important feature here is that no GWB strain amplitude of h c = 4 × 10 −15 can be obtained from a circular population of binaries with redshifts z < 1.0 and only very few eccentric populations could produce such amplitude in our sampling.Within the small number statistical uncertainties it seems that such a large amplitude is rarely achieved by any simulation within z < 1.0.While in general the results from Illustris and SIMBA in Figure 16 are very similar to those in Fig- ure 12, the merger rates for Illustris become nearly constant across all redshifts at amplitude of h c = 2 × 10 −15 already.

DISCUSSION AND CONCLUSIONS
The parametric astrophysical model presented in this work describes the intensity of the gravitational wave background as a function of the frequency.The focus was on the redshift dependent BH-bulge mass relation.By understanding the processes and relationships concerning the formation and co-evolution of galaxies and their central black holes, we have used an analytical expression in order to refine current astrophysical models.This allowed us to compare the predictions of this model with the constraints from PTA observations.Large-scale cosmological simulations help us to study the evolution of the Universe since observational unbiased data is hard to produce.We have fitted our redshift dependent BH-bulge relation to a suite of six simulations: EAGLE, Illustris, TNG100, TNG300, Horizon-AGN and SIMBA.The obtained best fit parameters serve as representative values for a Bayesian analysis.In general, all six simulations are consistent within ≤ 3.5σ with the range of shapes and strains of our simulated PTA GWB detections.The simulations can be broadly separated into two groups: 1. TNG100, TNG300 and SIMBA, which become more consistent with PTA detection as the GWB increases in amplitude and 2. EAGLE, Illustris and HorizonAGN, which behave in the opposite way.This separation coincidentally also follows the sign of the fitted γ * values of these simulations.
We simulated PTA detections to see how much they can help to constrain the posteriors of the parameters of the redshift dependent BH-bulge mass relation.As the redshift increases the value of γ * becomes more restricted.We find the tightest constraints for β * from a GWB detection in the PTA range, while α * does not change much from the prior.
Varying the maximum redshift parameter in the model seen in Figure A1 shows that the dominant fraction of the SMBHB population can be found withing z m ∼ 1.5 − 2.5 with the SMBHBs at higher redshifts only contributing a small, but not negligible, amount to the GWB.The study of higher redshift galaxies will be useful to determine the redshift evolution of BH-bulge mass relation.There still are difficulties to observe higher redshift galaxies.
Our proposed BH-bulge mass relation is a first-order extension of the standard redshift independent linear scaling relation.It can fit the masses from the simulations while maintaining approximately constant values of α * , β * , γ * and ε for redshifts z ≤ 5.The results depend on the specific parametric function and thus can only approximate the complexity of the masses given by the simulations.
Additionally, (Graham 2012) propose a double power law for the redshift independent BH-bulge mass relation of galaxies using observational data.Further studies to find and test the optimal shape of a redshift dependent BH-bulge mass relation are required.
Another interesting area for further improvement of the model is the galaxy stellar-bulge mass relation.The phenomenological stellar-bulge mass relation we have used is more suitable for elliptical and spheroidal galaxies, so a relation containing spiral galaxies including a degree of the spirality will be ideal to study a wide range of galaxies and their central black holes.Φ 0 -2.9 -2.6 -2.9 -2.9 -2.9 -2.55 -2.5 -2.5 -2.0.5 0.5 0.5 0.5 0.5 0.9 0.9 0.9 0.9 0.9 log 10 ρ 0 -0.1 -0.1 -0.1 -0.1 -0.1 -0.1 -0.1 -0.1 -0.1 -0.1 z 1.5 2.0 3.0 3.0 3.0 1.4 2.0 3.0 3.0 3.0 Table B1.Values of astrophysical parameters used to create the different simulated detections shown in Figure 4.Note that these are just one possible set for each spectrum and are neither unique nor necessarily representative.

Figure 1 .
Figure1.Best fit BH-bulge relations for the EAGLE, Illustris, TNG100, Horizon-AGN, SIMBA, and TNG300 simulations as they evolve with redshift.The BH and bulge masses with uncertainties from the simulations are consistent with the BH-bulge relations at different redshifts.

Figure 4 .
Figure 4. Simulated GWB detections with different characteristic spectra in the PTA range used in the Bayesian analysis.TableB1provides a nonunique set of parameter values for these spectra.To guide the eye the analytic sensitivity curve from the IPTA DR2(Antoniadis et al. 2022).
simulate a detection of the GWB we assume at each frequency a Gaussian distribution of central logarithmic amplitude log 10 A det ( f ) and width σ det ( f ), which are the detection measurement errors.With the GWB computed from a trial parameter set log 10 A trial ( f ) the likelihood function followingChen et al. (2017a) and Middleton

Figure 5 .
Figure5.The BH masses of the EAGLE, Illustris, TNG100, Horizon-AGN, SIMBA, and TNG300 simulations for all redshift z ≤ 5 as a function of the galaxy bulge mass.The maximum and minimum BH mass at the corresponding bulge mass are used to construct an allowed range for the BH-bulge relation, which is shown by the dashed lines.

Figure 6 .Figure 7 .
Figure 6.P-values from Kolmogorov-Smirnov tests for the fitted values from the simulations on the 1D marginalized posterior distributions for (α * , β * , γ * , ε).Crosses and circles indicate p-values from the straight line and curved spectra simulated detections respectively.

Figure 8 :Figure 9 :Figure 10 .
Figure 8: Posterior distributions of α * , β * , γ * , ε and maximum redshift z m for two different straight line characteristic strain spectra of 0.5 and 4 × 10 −15 at f = 1/1year are shown as black contours on the left and right panel respectively.The prior distributions are denoted by light green lines.The values of the parameters for the six large-scale cosmological simulation are also shown as coloured crosses, which can be used to judge the consistency between a simulation and the posteriors obtained from a simulated PTA GWB detection.

Figure 11 .
Figure11.Merger rates with respect to the chirp mass of the SMBHB for increasing characteristic strain values for both straight line (solid) and curved (dashed) spectra computed from the posterior distributions of the Bayesian analysis.The median, central 1 and 2 σ ranges are indicated by black, dark red and light green lines respectively.

Figure 12 .Figure 13 .
Figure12.Merger rates with respect to redshift for increasing characteristic strains values for both straight line (solid) and curved (dashed) spectra computed from the posterior distributions of the Bayesian analysis.The median, central 1 and 2 σ ranges are indicated by black, dark red and light green lines respectively.Each column represents a GWB detection at a strain amplitude.While each row represents the selection of SMBHBs within a given maximum redshift z m .

Figure 14 .
Figure14.Posterior distributions of selected astrophysical parameters for both straight line (solid) and curved (dashed) characteristic spectra of increasing strain with BH-bulge mass parameter values fixed to those representative of Illustris and SIMBA.The green shaded areas indicate the prior distributions for comparison.

Figure 15 .
Figure15.Same style as Figure11, but for merger rates with respect to the chirp mass of the SMBHB using values fixed to Illustris and SIMBA in the Bayesian analysis.
The amount of energy emitted as GWs by each individual binary dE d f r is dependent on the GW frequency in the source rest frame ( f r = (1+z) f ).The SMBHB merger rate (comoving number density in Mpc 3 of SMBHB mergers) per unit redshift and chirp mass

Table 1 .
shows the best fit values of γ * and ε for the simulations at different redshifts.The values are approximately constant across the Best fit parameters with uncertainties for the redshift dependent BH-bulge relation from Eqn (29) using BH and bulge masses from the six large-scale cosmological simulations.different redshifts, thus allowing us to use the average as a set of parameters that can approximately reproduce the masses from the simulations at all redshifts.An alternative redshift dependent BH-bulge mass relation used by e.g., Venemans et al. (2016) is written as

Table 2 .
Prior choice for the parameters of the redshift dependent BH -bulge mass relation.

Table 3 .
Prior choice for the parameters of the other astrophysical observables.
et al. (2018) can be written as