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

Musfar Muhamed Kozhikkal, Siyuan Chen, Gilles Theureau, Melanie Habouzit, Alberto Sesana

► To cite this version:

Musfar Muhamed Kozhikkal, Siyuan Chen, Gilles Theureau, Melanie Habouzit, Alberto Sesana. Mass-redshift dependency of Supermassive Black Hole Binaries for the Gravitational Wave Background. Monthly Notices of the Royal Astronomical Society, 2024, 531 (1), pp.1931-1950. 10.1093/mn-ras/stae1219. hal-04126263

HAL Id: hal-04126263

https://hal.science/hal-04126263v1

Submitted on 26 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Mass-redshift dependency of supermassive black hole binaries for the gravitational wave background

Musfar Muhamed Kozhikkal Ⓔ,<fnref n="1" />⋆ Siyuan Chen Ⓔ,<fnref n="2" />,<fnref n="3" />⋆ Gilles Theureau Ⓔ,<fnref n="1" />,<fnref n="4" />,<fnref n="5" /> Mélanie Habouzit Ⓔ,<fnref n="6" />,<fnref n="7" /> and Alberto Sesana Ⓔ,<fnref n="8" />,<fnref n="9" />,<fnref n="10" />

Accepted 2024 April 30. Received 2024 February 19; in original form 2023 May 31

## ABSTRACT

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 time-scale 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 $\le 3.5\sigma$ with a range of simulated GWB spectra. By fixing the BH–bulge mass parameters to the simulations we analyse 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.

**Key words:** black hole physics – gravitational waves – methods: analytical – galaxies: evolution – galaxies: formation.

## 1 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, and references therein).

Moreover, the SMBH pair formation process in the post-merger galaxy potential and their inspiral to coalescence, produces gravitational waves (GWs) in the low frequency domain, observable either as a stochastic gravitational wave background (GWB) or

as individual continuous gravitational wave sources with pulsar timing array (PTA) experiments (nHz–μHz) (Foster & Backer 1990; Rajagopal & Romani 1995; Jaffe & Backer 2003; Sesana, Vecchio & Colacino 2008), or with the future spatial laser interferometers like the Laser Interferometer Space Antenna (LISA) ($10^{-4}$–$10^{-2}$ Hz) (Amaro-Seoane et al. 2017, 2023).

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 yr for ongoing programmes) in order to extract a GWB from unmodelled noise. There are several PTA consortia, structured at continental levels and collaborating globally:

<footnote n="1">¹Laboratoire de Physique et Chimie de l'Environnement et de l'Espace, Université d'Orléans / CNRS, F-45071 Orléans Cedex 02, France</footnote>

<footnote n="2">²Shanghai Astronomical Observatory, Chinese Academy of Sciences, Shanghai 200030, P. R. China</footnote>

<footnote n="3">³Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, P. R. China</footnote>

<footnote n="4">⁴Observatoire Radioastronomique de Nançay, Observatoire de Paris, Université PSL, Université d'Orléans, CNRS, F-18330 Nançay, France</footnote>

<footnote n="5">⁵Laboratoire Univers et Théories LUTH, Observatoire de Paris, Université PSL, CNRS, Université de Paris, F-92190 Meudon, France</footnote>

<footnote n="6">⁶Zentrum für Astronomie der Universität Heidelberg, ITA, Albert-Ueberle-Str. 2, D-69120 Heidelberg, Germany</footnote>

<footnote n="7">⁷Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany</footnote>

<footnote n="8">⁸Dipartimento di Fisica 'G. Occhialini', Università degli Studi di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy</footnote>

<footnote n="9">⁹INFN, Sezione di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy</footnote>

<footnote n="10">¹⁰INAF – Osservatorio Astronomico di Brera, via Brera 20, I-20121 Milano, Italy</footnote>

<footnote>*E-mail: musfarmuhamed@gmail.com (MMK); siyuan.chen@shao.ac.cn (SC)</footnote>
European PTA (EPTA) (Kramer & Champion 2013; Desvignes et al. 2016; Chen et al. 2021), Parkes PTA (PPTA) in Australia (Hobbs 2013; Manchester et al. 2013; Kerr et al. 2020), North American Nanohertz Observatory for Gravitational Waves (NANOGrav) (Arzoumanian et al. 2016, 2018, 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 organization, the International PTA (IPTA), where they share their data and coordinate their analysis to eventually detect and hopefully characterize 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 data sets (Arzoumanian et al. 2020; Chen et al. 2021; Goncharov 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 InPTA have found evidence for (EPTA Collaboration 2023; Agazie et al. 2023a; Reardon et al. 2023; The International Pulsar Timing Array Collaboration 2024). In addition, the CPTA concurrently also found significant evidence for a Hellings-Downs correlated signal in their data set (Xu et al. 2023).

If these recently observed spectral signatures are from a population of SMBH binaries (SMBHBs), they favour heavy BH masses and short merger time-scales. 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 multimessenger 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 time-scale 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, Sesana & Del Pozzo 2017b; Chen, Sesana & Conselice 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 \sim 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 large-scale 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 (Crain et al. 2015; Schaye et al. 2015), Illustris (Genel et al. 2014; Vogelsberger et al. 2014), TNG100, TNG300 (Marinacci et al. 2018; Naiman et al. 2018; Nelson et al. 2018; Springel et al. 2018; Pillepich et al. 2018a, b), Horizon-AGN (Dubois et al. 2014, 2016), and SIMBA (Davé et al. 2019).

Our aim in this work is to set-up 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 time-scale. Section 3 focusses on the relation between the galaxy bulge and central BH 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 set-up, 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.

## 2 GWB CHARACTERISTIC STRAIN

For a population of SMBHBs the characteristic spectrum of the GWB was expressed in Pinney (2001) as

$$h_c^2(f) = \frac{4G}{\pi c^2 f} \int_0^\infty dz \int_0^\infty dM \frac{dE}{dM} \frac{d^2n}{dzdM}, \quad (1)$$
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

$$M = \frac{(M_1 M_2)^{3/5}}{(M_1 + M_2)^{1/5}}, \qquad (2)$$

where $M_1, M_2$ are the individual SMBH masses in the binary system. The amount of energy emitted as GWs by each individual binary $\frac{dE}{d\dot{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³ of SMBHB mergers) per unit redshift and chirp mass $\frac{d^2n}{d\dot{z}dM}$ can be derived from astrophysical observables or from a phenomenological function.

Below we summarize the parametric model from Chen et al. (2017b, 2019), which is extended by a parameter describing the redshift dependent evolution of the BH–bulge relation, see Section 3.

## 2.1 Individual binary

### 2.1.1 Analytic model and fitting function

Using the formalism of Chen et al. (2017b) we write $\frac{dE}{d\dot{f}_r}$ in terms of sum of harmonics at each eccentricity $e_n$ at each orbital frequency of the binary as

$$\frac{dE}{d\dot{f}_r} = \frac{M^{5/3}(\pi G)^{2/3}}{3(1+z)f^{1/3}} \sum_{n=1}^{\infty} \frac{g_n(e_n)}{F(e_n)(n/2)^{2/3}}, \qquad (3)$$

where

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

$$g_n(e) = \frac{n^4}{32} \left[ \begin{aligned} & (J_{n-2}(ne) - 2eJ_{n-1}(ne) + \frac{2}{n}J_n(ne) + 2eJ_{n+1}(ne) \\ & -J_{n+2}(ne))^2 + (1 - e^2)(J_{n-2}(ne) - 2J_n(ne) \\ & +J_{n+2}(ne))^2 + \frac{4}{3n^2}J_n^2(ne) \end{aligned} \right] \qquad (5)$$

and $J_n$ is the first kind of nth 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 \neq f_0$ the strain can be computed as

$$h_c(f) = h_{c,0} \left( f \frac{f_{p,0}}{f_{p,t}} \right) \left( \frac{f_{p,t}}{f_{p,0}} \right)^{-2/3}, \qquad (6)$$

with the peak frequency

$$f_p = \frac{1293 f}{181} \left[ \frac{e^{12/19}}{1 - e^2} \left( 1 + \frac{121 e^2}{304} \right)^{870/2299} \right]^{3/2} . \qquad (7)$$

A trial analytic function for the characteristic spectrum for the reference SMBHB with $\bar{f} = f/(10^{-8} \text{ Hz})$ can be written as

$$h_{c,\text{fit}}(f) = a_0 \bar{f}^{a_1} e^{-a_2 \bar{f}} + b_0 \bar{f}^{b_1} e^{-b_2 \bar{f}} + c_0 \bar{f}^{c_1} e^{-c_2 \bar{f}}. \qquad (8)$$

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

$$h_c^2(f) = \int_0^\infty dz \int_0^\infty dM \frac{d^2n}{d\dot{z}dM} h_{c,\text{fit}}^2 \left( f \frac{f_{p,0}}{f_{p,t}} \right) \left( \frac{f_{p,t}}{f_{p,0}} \right)^{-4/3} \\ \times \left( \frac{1+z}{1+z_0} \right)^{-1/3} \left( \frac{M}{M_0} \right)^{5/3}. \qquad (9)$$

### 2.1.2 Stellar environment

The GWB energy spectral shape is affected by the environmental coupling. A superefficient 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 gaseous environment (Sampson, Cornish & McWilliams 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

$$f_t = 0.356 \text{ nHz} \left( \frac{1}{F(e)} \frac{\rho_{i,100}}{\sigma_{200}} \rho_0 \right)^{3/10} M_9^{-2/5}, \qquad (10)$$

where the chirp mass $M_9 = M/(10^9)$ is re-scaled, $\rho_{i,100}$ is the stellar density of the environment within the SMBHB influence radius, the additional multiplicative factor $\rho_0$ includes all systematic uncertainties while estimating $\rho_{i,100}$ and $\sigma_{200}$ is the stellar central velocity dispersion in the galaxy, which are given by

$$\rho_{i,100} = \frac{\rho_i}{100} \approx \left( \frac{2M_{\text{BH}}}{M} \right)^{\gamma/(\gamma-3)} \frac{(3-\gamma)M}{400\pi a^3}, \qquad (11)$$

$$\sigma_{200} = \frac{\sigma}{200} = \frac{261}{200} \left( \frac{M_{\text{BH}}}{10^9} \right)^{0.228}. \qquad (12)$$

The stellar density distribution's inner slope is given by $\gamma \in (0.5, 2)$, *a* is the characteristic radius and *M* is the total bulge mass of the galaxy, which are expressed as

$$a = 239(2^{1/(3-\gamma)} - 1) \left(\frac{M}{10^9}\right)^{0.596}, \qquad (13)$$

$$M = 1.84 \times 10^{11} \left( \frac{M_{\text{BH}}}{10^9} \right)^{0.862}. \qquad (14)$$

## 2.2 Merger rate

The merger rate in equation (1) can be written in terms of SMBHB mass as

$$h_c^2(f) = \frac{4G}{\pi c^2 f} \int_0^\infty dz \int_0^\infty dM_{\text{BH}} \frac{dE}{d\dot{f}_r} \int_0^1 \frac{d^3n}{dz dM_{\text{BH}} dq_{\text{BH}}} dq_{\text{BH}}, \qquad (15)$$

where SMBHB merger rate and BH mass are

$$\frac{d3n}{dz dM_{\text{BH}} dq_{\text{BH}}} = \frac{d^3 n_G}{dz dM dq} \frac{dM}{dM_{\text{BH}}} \frac{dq}{dq_{\text{BH}}}, \qquad (16)$$

$$M_{\text{BH}} = \frac{M(1 + q_{\text{BH}})^{1/5}}{q_{\text{BH}}^{3/5}}. \qquad (17)$$

$M_{\text{BH}}$ can be parametrized 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), and Chen et al. (2019) as

$$\frac{d^3 n_G}{dz' dM dq} = \frac{\Phi(M, z) \mathcal{F}(z, M, q) dt_r}{M \ln(10) \tau(z, M, q) dz}, \qquad (18)$$

where *q* is the galaxy binary mass ratio with the primary galaxy mass *M*, $\Phi(M, z) = (dn_G/d\log M)_z$ is the GSMF estimated at redshift *z*, the differential pair fraction of the galaxy binaries is $\mathcal{F}(z, M, q) = (d\dot{f}/dq)_{z,M}$ and the merger time-scale $\tau(z, M, q) = \int_z^\infty (dt/d\bar{z})d\bar{z}$ is obtained by integrating over the instantaneous redshift $d\bar{z}$ between the redshifts at the start $z'$ and end $z$ of the galaxy merger. Using a
flat lambda cold dark matter (ΛCDM) model one finds

$$ \frac{dt}{dz} = \frac{1}{H_0(1+z)\sqrt{\Omega_M(1+z)^3 + \Omega_k(1+z)^2 + \Omega_\Lambda}} \quad (19) $$

Here, we use energy density ratios $\Omega_M = 0.3$, $\Omega_k = 0$, $\Omega_\Lambda = 0.7$ and Hubble constant $H_0 = 70$ km Mpc$^{-1}$ s$^{-1}$.

### 2.2.1 Galaxy stellar mass function

The 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 parametrized 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

$$ \Phi(M, z) = \ln(10)10^{\Phi_0 + z\Phi_I} \left(\frac{M}{M_0}\right)^{1+\alpha_0+z\alpha_I} \exp\left(-\frac{M}{M_0}\right). \quad (20) $$

### 2.2.2 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)

$$ \mathcal{F}(z, M, q) = f'_0 \left(\frac{M}{10^{11}}\right)^{\alpha_f} (1+z)^{\beta_f} q^{\gamma_f} = \frac{df_{\text{pair}}(z, M)}{dq} \quad (21) $$

with $f_0 = f'_0 \int q^{\gamma_f} dq$. Integrating over q then gives

$$ f_{\text{pair}} = f_0 \left( \frac{M}{10^{11}} \right)^{\alpha_f} (1+z)^{\beta_f}. \quad (22) $$

### 2.2.3 Merger time-scale

The time-scale of the evolution of a binary galaxy from the dynamical friction can be used to approximate the full merger time-scale, which can be written using a parametrization with $\tau_0$, $\alpha_\tau$, $\beta_\tau$, $\gamma_\tau$ as

$$ \tau(z, M, q) = \tau_0 \left( \frac{M h_0}{0.4 \times 10^{11}} \right)^{\alpha_\tau} (1+z)^{\beta_\tau} q^{\gamma_\tau}, \quad (23) $$

where $h_0 = 0.7$ is the Hubble parameter.

Substituting these observables into equation (18) gives

$$ \frac{d^3 n_G}{dz'dM dq} = \frac{10^{\Phi_0 + z\Phi_I} f'_0 \left(\frac{0.4}{h_0}\right)^{\alpha_\tau} \left(\frac{M}{10^{11}}\right)^{\alpha_f - \alpha_\tau} \left(\frac{M}{M_0}\right)^{\alpha_0 + z\alpha_I}}{M_0 \tau_0} \times e^{-M/M_0} (1+z)^{\beta_f - \beta_\tau} q^{\gamma_f - \gamma_\tau} \frac{dt}{dz}. \quad (24) $$

## 3 ASTROPHYSICS OF SMBH MASS

The final ingredient to describe the SMBHB merger rate in equation (16) is a relation between the galaxy stellar mass and the central BH 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

$$ M_{\text{bulge}} = \begin{cases} \left( \frac{\sqrt{6.9}}{(\log M - 10)^{1.5}} \exp \frac{-3.45}{\log M - 10} + 0.615 \right) M & \text{if } \log M > 10 \\ 0.615 M & \text{if } \log M \le 10. \end{cases} \quad (25) $$

This relation focusses on spherical and elliptical galaxies, which dominate the PTA GWB signal. Higher mass galaxies $M \le 10^{10} M_\odot$ have been observed to be correlated with the size of the bulge and disc, while lower mass galaxies do not.

## 3.1 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 equation (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 $\ge 100^3$ cMpc$^3$, a dark matter mass resolution of $\sim 5 \times 10^6 - 8 \times 10^7 M_\odot$, 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_\odot$ and their BHs. Baryonic processes taking place at small scales below the galactic scale are modelled as sub-grid 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 sub-grid physics of the simulations impact the evolution of both galaxies and BHs (Habouzit et al. 2021, 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, 2022b). The shape of the BH–bulge relation in the low-mass end ($M \le 10^{10.5} M_\odot$) 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).

## 3.2 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 BHs. The redshift independent BH–bulge mass relation (Kormendy & Ho 2013) which is usually used in the literature is given by

$$ M_{\text{BH}} \sim \mathcal{N} \left\{ \left( \frac{M_{\text{bulge}}}{10^{11} M_{\odot}} \right)^{\alpha} 10^{\beta}, \varepsilon \right\} \quad (26) $$

$$ \log_{10} M_{\text{BH}} = \alpha \log_{10} \left( \frac{M_{\text{bulge}}}{10^{11} M_{\odot}} \right) + \beta, \quad (27) $$
<table><caption>**Table 1.** Best-fitting parameters with uncertainties for the redshift dependent BH–bulge relation from equation (29) using BH and bulge masses from the six large-scale cosmological simulations.</caption><thead><tr><th>Simulation</th><th>α<sub>*</sub></th><th>β<sub>*</sub></th><th>γ<sub>*</sub></th><th>ε</th></tr></thead><tbody><tr><td>EAGLE</td><td>1.39 ± 0.027</td><td>8.23 ± 0.039</td><td>0.01<sup>+0.022</sup><sub>-0.035</sub></td><td>0.21<sup>+0.079</sup><sub>-0.076</sub></td></tr><tr><td>Illustris</td><td>1.28 ± 0.040</td><td>8.38 ± 0.088</td><td>0.18<sup>+0.046</sup><sub>-0.063</sub></td><td>0.08<sup>+0.144</sup><sub>-0.058</sub></td></tr><tr><td>TNG100</td><td>1.23 ± 0.022</td><td>8.91 ± 0.074</td><td>-0.02<sup>+0.025</sup><sub>-0.014</sub></td><td>0.16<sup>+0.078</sup><sub>-0.047</sub></td></tr><tr><td>HorizonAGN</td><td>1.03 ± 0.026</td><td>8.50 ± 0.036</td><td>0.07<sup>+0.008</sup><sub>-0.020</sub></td><td>0.08<sup>+0.032</sup><sub>-0.048</sub></td></tr><tr><td>SIMBA</td><td>1.24 ± 0.046</td><td>8.78 ± 0.063</td><td>-0.15<sup>+0.080</sup><sub>-0.064</sub></td><td>0.28<sup>+0.055</sup><sub>-0.050</sub></td></tr><tr><td>TNG300</td><td>1.29 ± 0.019</td><td>8.91 ± 0.050</td><td>-0.02<sup>+0.007</sup><sub>-0.007</sub></td><td>0.26<sup>+0.256</sup><sub>-0.115</sub></td></tr></tbody></table>

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. $\alpha$ and $\beta$ 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 $\varepsilon$ 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, Reines & Greene (2019).

### 3.3 BH–bulge mass relation with redshift dependence

Our goal is to formulate a parametric redshift dependent BH–bulge mass relation for $z \le 5$ since the GWB should be detectable with PTA in this redshift range. We propose a relation given by

$$M_{BH} \sim N \left\{ \left( \frac{M_{\text{bulge}}}{10^{11} M_{\odot}} \right)^{\alpha_{*}} 10^{\beta_{*} + \gamma_{*} z}, \varepsilon \right\}, \quad (28)$$

$$\log_{10} M_{BH} = \alpha_{*} \log_{10} \left( \frac{M_{\text{bulge}}}{10^{11} M_{\odot}} \right) + \beta_{*} + \gamma_{*} z, \quad (29)$$

where, we consider an additional BH–bulge mass parameter $\gamma_*$, which determines the extent of the evolution of the BH mass with redshift. Positive $\gamma_*$ values result in larger BH masses as the redshift increases, while negative $\gamma_*$ values have the opposite effect.

This relation is based on the assumption that the BH–bulge mass relation evolves only through the normalization parameter $\beta_*$, while the slope $\alpha_*$ remains constant with redshift. This assumption is based on the observation that the correlation between the mass of the BH 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 BH 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).

Equation (29) is fitted to the SMBH and galactic bulge masses for each of the six cosmological simulations separately. For a given simulation $\alpha_*$ is the slope in the logarithmic scale given by the linear least-squares fit over all redshifts $z \le 5$. $\beta_*$ is the intercept at $M_{\text{bulge}} = 10^{11} M_{\odot}$ of the least squares fit at $z = 0$. The intercepts on different redshifts is used to compute $\gamma_*$ The amount of scattering of the SMBH mass from the phenomenological fit in equation (29) is denoted by $\varepsilon$. Table 1 lists the BH–bulge mass relation parameters for these cosmological simulations. The variation of the masses for

these simulations are plotted in Fig. 1 as they evolve with redshift. Fig. 2 shows the best-fitting values of $\gamma_*$ and $\varepsilon$ for the simulations at different redshifts. The values are approximately constant across the 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

$$M_{BH} \sim N \left\{ \left( \frac{M_{\text{bulge}}}{10^{11} M_{\odot}} \right)^{\alpha_{*}} 10^{\beta_{*} (1+z)^{\gamma_{*}}} \varepsilon \right\}, \quad (30)$$

$$\log_{10} M_{BH} = \alpha_{*} \log_{10} \left( \frac{M_{\text{bulge}}}{10^{11} M_{\odot}} \right) + \beta_{*} + \gamma_{*} \log_{10}(1+z). \quad (31)$$

Fitting the BH 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 (equation 29) is a simple approximation to the simulations extending the redshift independent BH–bulge relation. Therefore, our best-fitting 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.

### 4 BAYESIAN ANALYSIS SET-UP

Using all the parts described above the characteristic spectrum can be computed as a function (equation 15) with 19 parameters which can be estimated from astrophysical observables. These parameters are five GSMF parameters $\Phi_0, \Phi_f, M_0, \alpha_0, \alpha_f$, four pair fraction parameters $f_0, \alpha_f, \beta_f, \gamma_f$, four merger time-scale parameters $\tau_0, \alpha_r, \beta_r, \gamma_r$, four BH–bulge mass parameters $\alpha_*, \beta_*, \gamma_*, \varepsilon$, and two parameters $e_0, \rho_0$ related to the individual binary GW emission.

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 Fig. 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 Fig. 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.

### 4.1 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 data sets for these two different shapes of spectrum for strain values of 0.5, 1, 2, 3, 4 $\times 10^{-15}$ at the reference frequency of $f = 1/1$ yr ($f \approx 10^{-7.5}$ Hz) as shown in Fig. 4. PTAs typically search at frequencies that are multiples of $1/T_{\text{span}}$, where $T_{\text{span}}$ is the total observation time span of the PTA data set. For simplicity and computational efficiency, we use the five lowest bins with $T_{\text{span}} = 25$ yr. The values of the
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Best-fitting 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.</figcaption></figure>

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.

## 4.2 Likelihood function

To simulate a detection of the GWB, we assume at each frequency a Gaussian distribution of central logarithmic amplitude $\log_{10} A_{\det}(f)$
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Variation of best-fitting $\gamma_*$ and $\epsilon$ values with redshift for the EAGLE, Illustris, TNG100, Horizon-AGN, SIMBA, and TNG300 simulations. The values are approximately constant across the redshift range for a given simulation, thus, one set of parameters ($\alpha_*, \beta_*, \gamma_*, \epsilon$) (see Table 1) can be used to represent the corresponding simulation.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** GWB characteristic strain spectra in the PTA range from a set of fiducial values (see Fig. A1) showing the differences when using the BH–bulge 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.</figcaption></figure>

<caption>**Figure 4.** Simulated GWB detections with different characteristic spectra in the PTA range used in the Bayesian analysis. Table B1 provides a non-unique set of parameter values for these spectra. To guide the eye the analytic sensitivity curve from the IPTA DR2 (Antoniadis et al. 2022).</caption>

and width $\sigma_{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 following Chen et al. (2017a) and Middleton et al. (2018) can be written as

$$p_{\text{det}}(d | A_{\text{trial}}(f)) \propto \exp \left\{ - \frac{(\log_{10} A_{\text{trial}}(f) - \log_{10} A_{\text{det}}(f))^2}{2\sigma_{\text{det}}(f)^2} \right\}. \quad (32)$$

The parallel tempering Markov Chain Monte Carlo (PTMCMC) sampler (Ellis & van Haasteren 2017) is used with $\log_{10} A_{det}(f)$ taken from the simulated GWB data sets and $\sigma_{det} = 0.09$ for the Bayesian analysis.

## 4.3 Prior choice

The prior for the BH–bulge mass relation is constrained using all possible masses from the six different simulations for $z \le 5$ to set the allowed range as shown in Fig. 5 and the initial test values are given the Table 2. Only combinations of $\alpha_*, \beta_*, \gamma_*$ 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 \in [0.1, 5]$, which includes the above range.
<figure><img src="image_5.png" /><figcaption>**Figure 5.** The BH masses of the EAGLE, Illustris, TNG100, Horizon-AGN, SIMBA, and TNG300 simulations for all redshift $z \le 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.</figcaption></figure>

<table><caption>**Table 2.** Prior choice for the parameters of the redshift dependent BH-bulge mass relation.</caption><thead><tr><th>Description</th><th>Parameter</th><th>Range</th></tr></thead><tbody><tr><td>BH-bulge mass relation average slope</td><td>α<sub>*</sub></td><td>[1,1.5]</td></tr><tr><td>BH-bulge mass relation norm at z = 0</td><td>β<sub>*</sub></td><td>[8,9]</td></tr><tr><td>BH-bulge mass relation norm redshift evolution</td><td>γ<sub>*</sub></td><td>[-0.5,0.5]</td></tr><tr><td>BH-bulge mass relation scatter</td><td>ε</td><td>[0.05,0.5]</td></tr><tr><td>Maximum redshift</td><td>z<sub>m</sub></td><td>[0.1,5]</td></tr></tbody></table>

<table><caption>**Table 3.** Prior choice for the parameters of the other astrophysical observables.</caption><thead><tr><th>Description</th><th>Parameter</th><th>Range</th></tr></thead><tbody><tr><td>GSMF norm</td><td>Φ<sub>0</sub></td><td>[-3.4, -2.4]</td></tr><tr><td>GSMF norm redshift evolution</td><td>Φ<sub>I</sub></td><td>[-0.6,0.2]</td></tr><tr><td>GSMF scaling mass</td><td>log<sub>10</sub> M<sub>0</sub></td><td>[11,11.5]</td></tr><tr><td>GSMF mass slope</td><td>α<sub>0</sub></td><td>[-1.5, -1]</td></tr><tr><td>GSMF mass slope redshift evolution</td><td>α<sub>I</sub></td><td>[-0.2,0.2]</td></tr><tr><td>Pair fraction norm</td><td>f<sub>0</sub></td><td>[0.01,0.05]</td></tr><tr><td>Pair fraction mass slope</td><td>α<sub>f</sub></td><td>[-0.5,0.5]</td></tr><tr><td>Pair fraction redshift slope</td><td>β<sub>f</sub></td><td>[0,2]</td></tr><tr><td>Pair fraction mass ratio slope</td><td>γ<sub>f</sub></td><td>[-0.2,0.2]</td></tr><tr><td>Merger time norm</td><td>τ<sub>0</sub></td><td>[0.1,10.0]</td></tr><tr><td>Merger time mass slope</td><td>α<sub>τ</sub></td><td>[-0.5,0.5]</td></tr><tr><td>Merger time redshift slope</td><td>β<sub>τ</sub></td><td>[-3,1]</td></tr><tr><td>Merger time mass ratio slope</td><td>γ<sub>τ</sub></td><td>[-0.2,0.2]</td></tr><tr><td>Binary eccentricity</td><td>e<sub>0</sub></td><td>[0.01,0.99]</td></tr><tr><td>Stellar density factor</td><td>log<sub>10</sub> ρ<sub>0</sub></td><td>[-2,2]</td></tr></tbody></table>

For the other parameters, we adopt the same prior choice as Chen et al. (2019) shown in Table 3.

# 5 RESULTS

## 5.1 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 pa-

rameters in our model. We first investigate the consistency of the fitting values that are an approximate representation of the complex simulations with the different shapes and strains of the simulated PTA detections.

Fig. 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 fitting 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 Fig. 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 Fig. 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.

## 5.2 Constraints on astrophysical observables

### 5.2.1 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 BH–galaxy bulge mass redshift dependent relation, which is why we focus on the parameters $α_*$, $β_*$, $ε$ and $z_m$ in the following. Figs 8 and 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 Fig. 8, there is little difference between the posteriors (black) and the priors (green; described in Section 4.3). All simulation fitting values lie within the allowed region. A detected amplitude of $h_c = 0.5 \times 10^{-15}$ provides little extra information. As the amplitude increases to $h_c = 4 \times 10^{-15}$, certain regions of the parameter space are ruled out. Noticeably, $β_*$ and $ε$ both show a tendency for larger values. As high-redshift BHs 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 Fig. 9, only marginally changes the findings from circular populations. The eccentricity and the
<figure><img src="image_6.png" /><figcaption>**Figure 6.** *P*-values from Kolmogorov–Smirnov tests for the fitting 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.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 7.** Mahalanobis distances between the fitting values from the simulations and the median values from the posterior constraints for all simulated data sets and both curved (circles) and straight line (crosses) spectra.</figcaption></figure>

environment of the BHs 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 Fig. 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 BH–bulge mass relation.

Given the PTA detections of a common signal of amplitude ~2.5 × 10⁻¹⁵ 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⁻¹⁵. However, we use a 25 yr observation time span, compared to the ~15 yr 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.

The corner plots for the complete 20 parameters with amplitudes $h_c = 0.5 \times 10^{-15}, 1 \times 10^{-15}, 2 \times 10^{-15}, 3 \times 10^{-15}$, and $4 \times 10^{-15}$ for both circular and eccentric population of SMBHBs are presented in the online supplementary material.

### 5.2.2 SMBHB merger rate

An interesting quantity that can be computed from our model is the merger rate of the SMBHBs from equation (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
<figure><img src="image_8.png" /><figcaption>**Figure 8.** Posterior distributions of $\alpha_*$, $\beta_*$, $\gamma_*$, $\epsilon$, and maximum redshift $z_m$ for two different straight line characteristic strain spectra of 0.5 and $4 \times 10^{-15}$ at $f = 1/1$ yr are shown as dark contours on the left and right panels, respectively. The prior distributions are denoted by light lines. The values of the parameters for the six large-scale cosmological simulation are also shown as crosses, which can be used to judge the consistency between a simulation and the posteriors obtained from a simulated PTA GWB detection.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 9.** Same style as Fig. 8, but for two different curved characteristic strain spectra of 0.5 and $4 \times 10^{-15}$ at $f = 1/1$ yr.</figcaption></figure>

the mass and redshift to get Figs 11 and 12 showing $dn/dM$ and $dn/dz$, respectively.

The merger rates with respect to the SMBHB mass in Fig. 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 \times 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
<figure><img src="image_10.png" /><figcaption>**Figure 10.** Posterior distributions of all the BH-bulge mass parameters including redshift variation (in the rows) as the detected characteristic strain value increases (in the columns) for both straight line (solid) and curved (dashed) spectra. For comparison the prior distributions are shown by the shaded areas.</figcaption></figure>

for smaller mass binaries in the 2σ 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 close by binaries. As such, we have binned the posterior samples by their maximum redshift. Within each bin, we plot the merger rate within a common range of redshifts in Fig. 12. For example 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 \times 10^{-15}$ the drop of the
<figure><img src="image_11.png" /><figcaption>**Figure 11.** 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 dark, medium and light lines, respectively.</figcaption></figure>

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 right most column in Fig. 12, where the 2σ lower bound drop of the merger rate moves from $z \approx 1$ to $z \approx 4$.

## 5.3 Constraints from simulations

By fixing the BH-bulge mass parameters to the best-fitting values from the simulations given in Table 1, we can see how the constraints on the other parameters are affected. For computational cost reasons we only analyse the $h_c = 0.5, 2, 4 \times 10^{-15}$ detections for both straight line and curved GWB spectra. Fig. 13 shows the median values and central 68 per cent 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 ($\Phi_0, \Phi_I, M_0, \alpha_0, \alpha_I$) are already well-constrained from observations. Parameters that play only a sub-dominant role, like those for the pair fraction ($f_0, \alpha_f, \beta_f, \gamma_f$) and maximum redshift $z_m$ are only slightly constrained towards larger values for stronger GWB strains. The eccentricity $\epsilon_0$ and stellar density $\rho_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 SMBHBs are fixed to results from cosmological simulations. The (second column, third row) panel in Fig. 13 on the merger time norm $\tau_0$ shows for all six simulations this decrease of the median values as well as the shrinkage of central 68 per cent credible regions. The other three parameters describing the merger

time ($\alpha_\tau, \beta_\tau, \gamma_\tau$) 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 \times 10^{-15}, 2 \times 10^{-15}$, and $4 \times 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.

### 5.3.1 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 Fig. 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 \times 10^{-15}$ and (2) $\beta_f$ at $h_c = 4 \times 10^{-15}$. At a characteristic spectrum value of $h_c = 0.5 \times 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 $\alpha_f$ for both Illustris and SIMBA shows a preference for low values at $h_c = 0.5 \times 10^{-15}$, followed by no preference at $h_c = 2 \times 10^{-15}$, and then higher values at $h_c = 4 \times 10^{-15}$ for both the curved and straight line spectra. The posterior distributions of the pair fraction redshift slope $\beta_f$ exhibit similar evolution with amplitude as in the case of $\alpha_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
<figure><img src="image_12.png" /><figcaption>**Figure 12.** 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 dark, medium and light 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$.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 13.** Median values and central 68 per cent credible regions of the 1D marginalized posteriors from the 16 parameters using the best-fitting values for the BH-bulge relation for the six simulations and three different strain values and shapes. The prior distributions are indicated by the left most point in each panel. Crosses indicate straight line spectra, while circles show curved spectra. A set of six points represent the results from one simulated detection case for all six simulations.</figcaption></figure>

larger amplitudes for both circular and eccentric populations. More massive and distant galaxy pairs are required to produce the GWB 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 \times 10^{-15}$ in SIMBA reveal a correlation between higher posterior values of $\alpha_\tau$ and higher values of eccentricity $e_0$ and $\rho_0$. This is in contrast to the general behaviour that increasing characteristic spectrum values lead to lower values of $\alpha_\tau$, $\beta_\tau$, and $\tau_0$. The posteriors change from

being in broad agreement with the priors at $h_c = 0.5 \times 10^{-15}$ for both curved and straight line spectra and both simulations to trending very clearly towards lower merger times at $h_c = 4 \times 10^{-15}$. The main difference between Illustris and SIMBA seems to be that the merger time redshift slope $\beta_\tau$ is more constrained for Illustris, whereas it is the merger time norm $\tau_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 $\rho_0$ parameters. The straight line spectra lead to almost no constraints at all strains. On the other
<figure><img src="image_14.png" /><figcaption>**Figure 14.** 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 shaded areas indicate the prior distributions for comparison.</figcaption></figure>

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 $\rho_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 Fig. 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 $\gamma_*$ 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.

### 5.3.2 Merger rate 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 Figs 15 and 16. The most prominent feature of $dn/dM$ in Fig. 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 Fig. 14. The other feature that is visible from Fig. 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.
<figure><img src="image_15.png" /><figcaption>**Figure 15.** Same style as Fig. 11, but for merger rates with respect to the chirp mass of the SMBHB using values fixed to Illustris and SIMBA in the Bayesian analysis.</figcaption></figure>

The merger rate with respect to the redshift is shown in Fig. 16 for Illustris and SIMBA with the panels defined in the same way as in the previous section and Fig. 12. An important feature here is that no GWB strain amplitude of $h_c = 4 \times 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 Fig. 16 are very similar to those in Fig. 12, the merger rates for Illustris become nearly constant across all redshifts at amplitude of $h_c = 2 \times 10^{-15}$ already.

## 6 DISCUSSION AND CONCLUSIONS

The parametric astrophysical model presented in this work describes the intensity of the GWB 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 BH, 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 are 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-fitting parameters

serve as representative values for a Bayesian analysis. In general, all six simulations are consistent within $\le 3.5\sigma$ 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 $\gamma_*$ 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 $\gamma_*$ becomes more restricted. We find the tightest constraints for $\beta_*$ from a GWB detection in the PTA range, while $\alpha_*$ does not change much from the prior.

Varying the maximum redshift parameter in the model seen in Fig. A1 shows that the dominant fraction of the SMBHB population can be found within $z_m \sim 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 $\alpha_*$, $\beta_*$, $\gamma_*$, and $\varepsilon$ for redshifts $z \le 5$. The results depend on the specific parametric function and thus
<figure><img src="image_16.png" /><figcaption>**Figure 16.** Same style as Fig. 12, but for merger rates with respect to redshift using values fixed to Illustris and SIMBA in the Bayesian analysis. The left three columns are for Illustris, whereas the right three columns are for SIMBA.</figcaption></figure>
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 us-
ing observational data. Further studies to find and test the op-
timal 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 BHs.

ACKNOWLEDGEMENTS

We would like to thank Marta Volonteri for helpful discussions. We
acknowledge the support of our colleagues in the European Pulsar
Timing Array. MH acknowleges support from the Gliese fellowship
of the Zentrum für Astronomie and the MPIA fellowship of the
Max-Planck-Institut für Astronomie. AS acknowledges the financial
support provided under the European Union’s H2020 ERC Consol-
idator Grant ‘Binary Massive Black Hole Astrophysics’ (B Massive,
Grant Agreement: 818691). We acknowledge financial support
from ‘Programme National de Cosmologie and Galaxies’ (PNCG),
and ‘Programme National Hautes Energies’ (PNHE) funded by
CNRS/INSU-IN2P3-INP, CEA and CNES, France. We acknowledge
financial support from Agence Nationale de la Recherche (ANR-18-
CE31-0015), France.

DATA AVAILABILITY

The data used in this article shall be shared on reasonable request
to the corresponding authors. Supplementary material can be found
online at doi:10.5281/zenodo.11118728.

REFERENCES

Agazie G. et al., 2023a, ApJ, 951, L8

Agazie G. et al., 2023b, ApJ, 952, L37

Amaro-Seoane P. et al., 2017, preprint (arXiv:1702.00786)

Amaro-Seoane P. et al., 2023, Living Rev. Relativ., 26, 2

Angulo R. E., Hahn O., 2022, Living Rev. Comput. Astrophys., 8, 1

Antoniadis J. et al., 2022, MNRAS, 510, 4873

Arzoumanian Z. et al., 2016, ApJ, 821, 13

Arzoumanian Z. et al., 2018, ApJ, 859, 47

Arzoumanian Z. et al., 2020, ApJ, 905, L34

Bailes M. et al., 2020, PASA, 37, e028

Bernardi M., Meert A., Vikram V., Huertas-Company M., Mei S., Shankar
F., Sheth R. K., 2014, MNRAS, 443, 874

Chen S., Middleton H., Sesana A., Del Pozzo W., Vecchio A., 2017a,
MNRAS, 468, 404

Chen S., Sesana A., Del Pozzo W., 2017b, MNRAS, 470, 1738

Chen S., Sesana A., Conselice C. J., 2019, MNRAS, 488, 401

Chen S. et al., 2021, MNRAS, 508, 4970

Conselice C. J., Wilkinson A., Duncan K., Mortlock A., 2016, ApJ, 830, 83

Crain R. A. et al., 2015, MNRAS, 450, 1937

Davé R., Anglés-Alcázar D., Narayanan D., Li Q., Rafieferantsoa M. H.,
Appleby S., 2019, MNRAS, 486, 2827

Desvignes G. et al., 2016, MNRAS, 458, 3341

Detweiler S., 1979, ApJ, 234, 1100

Ding X. et al., 2020, ApJ, 888, 37

Dubois Y. et al., 2014, MNRAS, 444, 1453

Dubois Y., Peirani S., Pichon C., Devriendt J., Gavazzi R., Welker C.,
Volonteri M., 2016, MNRAS, 463, 3948

Ellis J., van Haasteren R., 2017, jellis18/PTMCMCMSampler: Official Release.
Zenodo, available at: https://zenodo.org/records/1037579

EPTA Collaboration, InPTA Collaboration, 2023, A&A, 678, A50

EPTA Collaboration, InPTA Collaboration, 2024, A&A, 685, A94

Foster R. S., Backer D. C., 1990, ApJ, 361, 300

Genel S. et al., 2014, MNRAS, 445, 175

Genel S. et al., 2018, MNRAS, 474, 3976

Goncharov B. et al., 2021, ApJ, 917, L19

Graham A. W., 2012, ApJ, 746, 113

Habouzit M. et al., 2021, MNRAS, 503, 1940

Habouzit M. et al., 2022a, MNRAS, 509, 3015

Habouzit M. et al., 2022b, MNRAS, 511, 3751

Hellings R. W., Downs G. S., 1983, ApJ, 265, L39

Hobbs G., 2013, Class. Quant. Grav., 30, 224007

Hobbs G. et al., 2010, Class. Quant. Grav., 27, 084013

Huerta E. A., McWilliams S. T., Gair J. R., Taylor S. R., 2015, Phys. Rev. D,
92, 063010

Jaffe A. H., Backer D. C., 2003, ApJ, 583, 616

Jiang P. et al., 2019, Sci. China Phys. Mech. Astron., 62, 959502

Joshi B. C. et al., 2018, J. Astrophys. Astron., 39, 51

Kerr M. et al., 2020, PASA, 37, e020

Kormendy J., Ho L. C., 2013, ARA&A, 51, 511

Kramer M., Champion D. J., 2013, Class. Quant. Grav., 30, 224009

Lee K. J., 2016, in Qian L., Li D., eds, ASP Conf. Ser. Vol. 502, Frontiers in
Radio Astronomy and FAST Early Sciences Symposium 2015. Astron.
Soc. Pac., San Francisco, p.19

Manchester R. N. et al., 2013, PASA, 30, e017

Marinacci F. et al., 2018, MNRAS, 480, 5113

Merloni A. et al., 2010, ApJ, 708, 137

Middleton H., Chen S., Del Pozzo W., Sesana A., Vecchio A., 2018, Nat.
Commun., 9, 573

Mortlock A. et al., 2015, MNRAS, 447, 2

Mundy C. J., Conselice C. J., Duncan K. J., Almaini O., Häußler B., Hartley
W. G., 2017, MNRAS, 470, 3507

Naiman J. P. et al., 2018, MNRAS, 477, 1206

Nelson D. et al., 2018, MNRAS, 475, 624

Phinney E. S., 2001, preprint (arXiv:astro-ph/0108028)

Pillepich A. et al., 2018a, MNRAS, 473, 4077

Pillepich A. et al., 2018b, MNRAS, 475, 648

Pol N. S. et al., 2021, ApJ, 911, L34

Rajagopal M., Romani R. W., 1995, ApJ, 446, 543

Ravi V., Wyithe J. S. B., Shannon R. M., Hobbs G., Manchester R. N., 2014,
MNRAS, 442, 56

Reardon D. J. et al., 2023, ApJ, 951, L6

Sampson L., Cornish N. J., McWilliams S. T., 2015, Phys. Rev. D, 91, 084055

Sazhin M. V., 1978, Soviet Ast., 22, 36

Schaye J. et al., 2015, MNRAS, 446, 521

Schutte Z., Reines A. E., Greene J. E., 2019, ApJ, 887, 245

Sesana A., 2013a, Class. Quant. Grav., 30, 224014

Sesana A., 2013b, MNRAS, 433, L1

Sesana A., Vecchio A., Colacino C. N., 2008, MNRAS, 390, 192

Sesana A., Shankar F., Bernardi M., Sheth R. K., 2016, MNRAS, 463, L6

Spiewak R. et al., 2022, PASA, 39, e027

Springel V. et al., 2018, MNRAS, 475, 676

Tarafdar P. et al., 2022, PASA, 39, e053

The International Pulsar Timing Array Collaboration, 2024, ApJ, 966, 105

Venemans B. P., Walter F., Zschaechner L., Decarli R., De Rosa G., Findlay
J. R., McMahon R. G., Sutherland W. J., 2016, ApJ, 816, 37

Verbiest J. P. W. et al., 2016, MNRAS, 458, 1267

Vogelsberger M. et al., 2014, MNRAS, 444, 1518

Vogelsberger M., Marinacci F., Torrey P., Puchwein E., 2020, Nat. Rev. Phys.,
2, 42

Xu H. et al., 2023, Res. Astron. Astrophys., 23, 075024
APPENDIX A: GWB CHARACTERISTIC SPECTRA FROM VARIATIONS OF THE ASTROPHYSICAL PARAMETERS

<figure><img src="image_17.png" /><figcaption>**Figure A1.** Effect on the GWB spectrum from variations of each of the 20 astrophysical parameters within the range given in Tables 2 and 3. The default values for the parameters are set as: $\Phi_0 = -2.6$, $\Phi_I = -0.45$, $M_0 = 11.25$, $\alpha_0 = -1.15$, $\alpha_I = -0.1$, $f_0 = 0.02$, $\alpha_f = 0.1$, $\beta_f = 0.8$, $\gamma_f = 0.1$, $\tau_0 = 0.8$, $\alpha_\tau = -0.1$, $\beta_\tau = -2$, $\gamma_\tau = -0.1$, $\alpha_*= 1.1$, $\beta_* = 8.2$, $\gamma_* = -0.2$, $\varepsilon = 0.3$, $e_0 = 0.9$, $\log_{10}\rho_0 = 0.1$, $z_m = 2.0$. Each panel shows the change in the GWB spectrum by varying only one parameter.</figcaption></figure>
APPENDIX B: PARAMETER VALUES FOR THE SIMULATED DATA SETS

<table><caption>**Table B1.** Values of astrophysical parameters used to create the different simulated detections shown in Fig. 4. Note that these are just one possible set for each spectrum and are neither unique nor necessarily representative.</caption><thead><tr><th rowspan="2">Parameter<br>h<sub>c</sub> at f = 1/1yr</th><th colspan="8">Line</th><th colspan="3">Curve</th></tr><tr><th>0.5 x 10<sup>-15</sup></th><th>1 x 10<sup>-15</sup></th><th>2 x 10<sup>-15</sup></th><th>3 x 10<sup>-15</sup></th><th>4 x 10<sup>-15</sup></th><th>0.5 x 10<sup>-15</sup></th><th>1 x 10<sup>-15</sup></th><th>2 x 10<sup>-15</sup></th><th>3 x 10<sup>-15</sup></th><th>4 x 10<sup>-15</sup></th></tr></thead><tbody><tr><td>Φ<sub>0</sub></td><td>-2.9</td><td>-2.6</td><td>-2.9</td><td>-2.9</td><td>-2.9</td><td>-2.55</td><td>-2.5</td><td>-2.5</td><td>-2.5</td><td>-2.6</td></tr><tr><td>Φ<sub>1</sub></td><td>-0.45</td><td>-0.45</td><td>-0.45</td><td>-0.45</td><td>-0.45</td><td>-0.45</td><td>-0.255</td><td>-0.1</td><td>0.08</td><td>0.095</td></tr><tr><td>M<sub>0</sub></td><td>11.25</td><td>11.25</td><td>11.25</td><td>11.25</td><td>11.3</td><td>11.35</td><td>11.35</td><td>11.35</td><td>11.35</td><td>11.2</td></tr><tr><td>α<sub>0</sub></td><td>-1.15</td><td>-1.15</td><td>-1.15</td><td>-1.15</td><td>-1.15</td><td>-1.1</td><td>-1.1</td><td>-1.1</td><td>-1.1</td><td>-1.1</td></tr><tr><td>α<sub>1</sub></td><td>-0.1</td><td>-0.1</td><td>-0.1</td><td>-0.1</td><td>-0.1</td><td>-0.12</td><td>-0.12</td><td>-0.12</td><td>-0.12</td><td>-0.12</td></tr><tr><td>f<sub>0</sub></td><td>0.015</td><td>0.02</td><td>0.03</td><td>0.03</td><td>0.035</td><td>0.022</td><td>0.03</td><td>0.03</td><td>0.03</td><td>0.03</td></tr><tr><td>α<sub>f</sub></td><td>0.3</td><td>0.3</td><td>0.3</td><td>0.3</td><td>0.3</td><td>-0.15</td><td>-0.15</td><td>-0.15</td><td>-0.15</td><td>-0.15</td></tr><tr><td>β<sub>f</sub></td><td>0.6</td><td>0.6</td><td>0.6</td><td>0.6</td><td>0.6</td><td>0.8</td><td>1.2</td><td>1.48</td><td>1.3</td><td>1.7</td></tr><tr><td>γ<sub>f</sub></td><td>0.1</td><td>0.1</td><td>0.1</td><td>0.1</td><td>0.1</td><td>0.1</td><td>0.1</td><td>0.1</td><td>0.1</td><td>0.1</td></tr><tr><td>τ<sub>0</sub></td><td>2.</td><td>1.8</td><td>2.</td><td>2.</td><td>2.</td><td>0.8</td><td>0.8</td><td>0.8</td><td>0.8</td><td>0.8</td></tr><tr><td>α<sub>τ</sub></td><td>0.1</td><td>0.1</td><td>0.1</td><td>0.1</td><td>0.1</td><td>0.2</td><td>0.2</td><td>0.2</td><td>0.2</td><td>-0.2</td></tr><tr><td>β<sub>τ</sub></td><td>-1.8</td><td>-2.</td><td>-2.1</td><td>-2.3</td><td>-2.5</td><td>-2.</td><td>-2.</td><td>-2.</td><td>-2.1</td><td>-2.1</td></tr><tr><td>γ<sub>τ</sub></td><td>-0.1</td><td>-0.1</td><td>-0.1</td><td>-0.1</td><td>-0.1</td><td>-0.1</td><td>-0.1</td><td>-0.1</td><td>-0.1</td><td>0.1</td></tr><tr><td>α<sub>*</sub></td><td>1.1</td><td>1.1</td><td>1.1</td><td>1.1</td><td>1.1</td><td>1.</td><td>1.</td><td>1.</td><td>1.</td><td>1.</td></tr><tr><td>β<sub>*</sub></td><td>8.5</td><td>8.5</td><td>8.65</td><td>8.8</td><td>8.8</td><td>8.0</td><td>8.0</td><td>8.0</td><td>8.0</td><td>8.1</td></tr><tr><td>γ<sub>*</sub></td><td>0.1</td><td>0.1</td><td>0.3</td><td>0.3</td><td>0.3</td><td>0.1</td><td>0.1</td><td>0.1</td><td>0.1</td><td>0.1</td></tr><tr><td>ε</td><td>0.35</td><td>0.35</td><td>0.35</td><td>0.35</td><td>0.35</td><td>0.2</td><td>0.2</td><td>0.2</td><td>0.2</td><td>0.2</td></tr><tr><td>e<sub>0</sub></td><td>0.5</td><td>0.5</td><td>0.5</td><td>0.5</td><td>0.5</td><td>0.9</td><td>0.9</td><td>0.9</td><td>0.9</td><td>0.9</td></tr><tr><td>log<sub>10</sub> ρ<sub>0</sub></td><td>-0.1</td><td>-0.1</td><td>-0.1</td><td>-0.1</td><td>-0.1</td><td>-0.1</td><td>-0.1</td><td>-0.1</td><td>-0.1</td><td>-0.1</td></tr><tr><td>z</td><td>1.5</td><td>2.0</td><td>3.0</td><td>3.0</td><td>3.0</td><td>1.4</td><td>2.0</td><td>3.0</td><td>3.0</td><td>3.0</td></tr></tbody></table>

This paper has been typeset from a TeX/LATEX file prepared by the author.