© The Authors 2024

1. Introduction

Third-generation detectors of wide-band gravitational waves (GWs) such as the Einstein Telescope (ET; Hild et al. 2011; Punturo et al. 2010) or the Cosmic Explorer (CE; Dwyer et al. 2015; Abbott et al. 2017; Reitze et al. 2019) will be able to probe a much larger volume of the Universe, and will therefore have a much higher detection rate (Maggiore et al. 2020) compared to the current second-generation detectors. ET will have a detection sensitivity down to 1 Hz (Hild et al. 2008; Hild 2012) and will therefore have the ability to detect binary black holes (BBHs) of high mass; that is, in the range of 10² − 10⁴ M_⊙ (Huerta & Gair 2011a,b; Gair et al. 2011; Amaro-Seoane & Santamaría 2010). Assuming ET-D (Hild et al. 2011) design sensitivity, the expected detection rates are ∼10⁵ − 10⁶ BBH detections and ∼7 × 10⁴ binary neutron star (BNS) detections in one year (Regimbau et al. 2012, 2014; Belgacem et al. 2019). Given its increased detection sensitivity, the ET will also be able to detect the coalescence of compact binaries with a total mass of 20–100 M_⊙, typical of black hole–black hole (BH–BH) or black hole–neutron star (BH–NS) binaries, up to redshift z ≈ 20 and even higher. LIGO-Virgo-Kagra (LVK) has already detected 90 gravitational wave events to date. The most recent detections are presented in the GWTC-3 catalogue (Abbott et al. 2023). Given that third-generation gravitational wave detectors will have high detection rates and redshift reach, and will also be able to provide strong constraints on population properties thanks to the smaller uncertainties on the physical parameters of the GW events, (Singh et al. 2022; Yi et al. 2022; Van Den Broeck 2010) in this paper, we investigate the prospects of using the ET as a single instrument to reconstruct the functional form of the star formation rate (SFR) for a redshift range of 0 ≤ z ≤ 20 for different compact binary populations.

The SFR density for Population (Pop) I stars for low redshifts is already very well constrained (Madau & Dickinson 2014), but there are large uncertainties for the first stars in the Universe, also known as Pop III stars. Pop III stars are expected to be the first sources of light and play a crucial role in the early cosmic evolution by producing the very first heavy elements (Bromm & Larson 2004; Bromm et al. 2009). Several observational methods have been used to probe this early stellar population (Komissarov & Barkov 2010; Toma et al. 2011; Mészáros & Rees 2010; Campisi et al. 2011; Ma et al. 2015), but direct observation of these Pop III stars remains to be achieved. The remnants of Pop III stars have been studied as possible sources of gravitational waves (Bond & Carr 1984; Belczynski et al. 2004; Kulczycki et al. 2006; Kowalska et al. 2012; Kinugawa et al. 2014; Hartwig et al. 2016).

There are some major differences between the evolution of Pop III stars and that of metal-polluted Pop I and II compact binaries (see Bromm & Larson 2004 for detailed review). Pop III stars are expected to be more massive than Pop I stars (Hosokawa et al. 2012; Stacy et al. 2012), with masses of 10 − 100 M_⊙, and so the Pop III star binaries are expected to evolve into BBHs. The main differences are as follows. (i) The initial conditions, such as initial mass function, initial binary mass ratio, and initial separations and eccentricities are different (Belczynski et al. 2017). (ii) The wind mass loss is very different. Heavy wind mass loss is seen for Pop I stars, smaller wind mass loss for Pop II stars, and almost no wind mass loss for Pop III stars. This in turn affects the binary separation evolution as wind mass loss widens the separation, meaning that Pop III stars are subject to more frequent binary interactions, such as Roche-lobe overflow (RLOF) and common envelope (CE) interactions. The wind mass loss affects the neurton star and black hole mass, as there is more mass available for Pop III stars to form compact remnants (for a description of the effect of winds on mass, see Belczynski et al. 2010). (iii) There are also differences in radial evolution. Pop III stars have smaller radii than Pop I and II stars (Belczynski et al. 2017). This leads to a smaller number of binary interactions for Pop III stars (RLOF and CE) as compared with Pop I and II stars and thus counter balances the effect of wind loss. As there has not yet been direct detection of Pop III stars, the uncertainties due to different evolutionary parameters are still to be verified. It is therefore necessary to obtain strong constraints on the merger rate densities as a function of redshift in order to have a better understanding of the evolution of these first stars.

Assuming that the Pop III stars are formed by the collapse of dark matter halos, the SFR depends on multiple factors, including halo mass, reionisation history, metal enrichment in the intergalactic medium, and accretion rate. The ability of a primordial gas to cool and condense in the early Universe – which in turn depends on the size of the ‘mini-halos’ – is one of the factors affecting star formation efficiency. The star formation inside mini-halos, which are smaller in mass than the critical halo mass (Yoshida et al. 2003), can be suppressed by ultraviolet background in the Lyman-Werner bands. Metal enrichment in the intergalactic medium – where the main contribution is from supernova explosions – determines when the formation of first stars will terminate. (De Souza et al. 2011 explored many such factors in calculating the SFR of Pop III stars.) Therefore, any observational constraint on the SFR will be a crucial step in constraining formation scenarios.

For any given population of stars, three major questions regarding the SFR are as follows. (i) Firstly, we do not yet know when the star formation of Pop III stars terminated; (ii) secondly, the redshift at which the SFR peaks in unknown; and finally, (iii) the functional form of SFR at high redshift remains to be determined. In this paper, we aim to find the answers to these questions with a given set of detections with single ET. We define the termination redshift as the redshift where the SFR drops to 1% of its peak value. Vitale et al. (2019) specifically showed how detections of GWs from inspiraling BBHs by third-generation detectors can be used to measure the SFR of massive stars with high precision up to redshifts of approximately 10 assuming that all sources come from galactic fields. The authors mention that they assume that the time-delay distribution is the same for all sources at all redshifts and neglected the dependence of this distribution on the mass and spin of the source. In order to estimate the redshift of the inspiralling binaries, these latter authors assume a three-detector 3G network. In this paper, we simulate multiple mock populations for compact binaries originating from Pop I+II (field binaries) and Pop III stars using a more realistic distribution of time delay to construct the mock population. While estimating the SFR we make no assumptions about the originating population of the compact binary or about the functional form of the SFR. We estimate the parameters of the compact binaries and the SFR with ET as a single instrument.

2. Plan of the paper

In Singh & Bulik (2021; SB1 hereafter), we developed an algorithm to break the chirp mass–redshift degeneracy in the detected GW signal from the coalescence of a compact binary with ET as a single instrument, and thus to estimate the parameters of the merging compact binary. We estimated the area of localisation, chirp mass, redshift, and mass ratios by estimating their posterior distribution for short-duration GW signals from inspiraling compact binary systems. In the subsequent work in the series (Singh & Bulik 2022; SB2 hereafter), we further developed the algorithm, taking into account the effect of the rotation of the Earth on the antenna pattern function in order to analyse the long-duration signals from coalescing low-mass compact binary systems. We used this algorithm to further analyse realistic populations of compact binary systems originating from Pop I+II, Pop III, and globular cluster (GC) populations in (Singh et al. 2022; S22 hereafter). In S22, we concluded that ET as a single instrument is capable of detecting and distinguishing different compact binary populations separated in chirp mass–redshift space. We also estimated the merger rate density and found that, although our estimates for Pop I+II and GC populations are in good agreement with the true merger rate density of the respective populations, the deviation from the true value is much larger in the case of Pop III, especially for higher redshifts of z > 7 (see Fig. 7 in S22). In the present work, we use an improved version of the algorithm to estimate the parameters of the merging compact binaries.

In Sect. 3 we describe the population models and the SFRs used to generate our mock populations. We discuss estimations of the parameters of the detected compact binaries in Sect. 4. To obtain an estimate of the SFR of a given population, it is crucial to have a correct assessment of the merger rate density as a function of redshift. In order to estimate the merger rate density accurately, we calculate the detection efficiency. Taking into account the fraction of binaries that do not cross the detection threshold, we reconstruct the true merger rate density of the population. This is described in detail in Sect. 5. Subsequently, assuming a functional form for the delay time, which is the time from the formation of the stars to the merger of the compact binaries formed from those stars, we proceed to reconstruct the SFR of the population. This process is described in detail in Sect. 6. We discuss our conclusions in Sect. 7.

3. Mock source catalogue

As mentioned in the previous section, in S22 we showed that ET as a single instrument is capable of detecting and distinguishing different compact binary populations separated in chirp mass–redshift space. Based on this earlier result, in the present work, we simulate three mock populations for compact binaries originating from Pop I+II and Pop III stars: Two mock populations are generated that consist of only compact binaries of Pop III stars using two realistic star formation densities evolving over redshift. We assume that compact binaries originating from these populations are clearly distinguishable from other populations. One mock population is generated consisting of compact binaries from both Pop I+II and Pop III stars, allowing us to estimate the SFR without assuming the ability to distinguish between these populations.

We use the FS1 model of the population of compact binary systems from the first, metal-free Pop III stars generated by Belczynski et al. (2017). The initial conditions for generating this model are based on the models obtained by Ryu et al. (2016) using N-body simulation, and assuming the formation of Pop III stars from a mini halo of ∼2000 AU in size (Stacy & Bromm 2013). The number density of the gas medium is chosen to be 10⁶ cm⁻³, following the parameters specified by Stacy & Bromm (2013). The details of this population of metal-free binaries, such as initial mass function, mass ratio, orbital separations, and eccentricities, are described in Belczynski et al. (2017).

De Souza et al. (2011) considered two populations of Pop III stars: (i) Pop III.1 stars, which are the first-generation stars formed from initial conditions determined cosmologically, and (ii) Pop III.2 stars, which are zero-metallicity stars that formed from a primordial gas, influenced by an earlier generation of stars. Pop III.2 stars are expected to form in an initially ionised gas (Johnson & Bromm 2006; Yoshida et al. 2007) and are thought to be less massive ( ∼ 40 − 60 M_⊙) than Pop III.1 stars ( ∼ 1000 M_⊙). De Souza et al. (2011) calculated the SFR for Pop III.1 and Pop III.2 using three different values of the galactic wind – namely v_wind = 50, 75, 100 km s⁻¹ – in order to incorporate the effect of metal enrichment by galactic winds, and two different star formation efficiency values, f_⋆ = 0.001 and f_⋆ = 0.1. In the present paper, we consider three SFRs calculated by De Souza et al. (2011) for Pop III.2 stars: (i) v_wind = 50 km s⁻¹ and f_⋆ = 0.001, (ii) v_wind = 100 km s⁻¹ and f_⋆ = 0.001, and (iii) a very optimistic case with v_wind = 50 km s⁻¹ and f_⋆ = 0.1. We refer to these SFRs as SFR1, SFR2, and SFR3, respectively.

These SFRs are shown in Fig. 1 and we see that the major effect of metal enrichment by galactic winds in the case of SFR1 and SFR2 is manifested in the termination (the redshift where the SFR drops to 1% of its peak value) of these SFRs. SFR2 has greater galactic wind velocity as compared to SFR1 and we see that SFR2 terminates at z ∼ 5.3, whereas SFR1 terminates at a much later time, at z ∼ 3.2. Therefore, any constraint on the termination redshift will be helpful in providing information about the formation scenarios of these first stars (for details see Bromm & Larson 2004). While the metal enrichment via galactic winds is the same for SFR1 and SFR3, the increased star formation efficiency, f_⋆, leads to an overall increase in the absolute value of SFR as a function of redshift.

Fig. 1. Star formation rates used to construct the mock populations in this analysis. SFRs for Pop III are adopted from De Souza et al. (2011) and the SFR for Pop I+II is taken from Belczynski et al. (2020).

For the Pop I+II compact binaries, we use the M30B generated by Belczynski et al. (2020) and an upgraded version of the population synthesis code StarTrack (Belczynski et al. 2002, 2008). These authors generated multiple binary stellar evolution models consistent with LVK O1/O2 merger rates for BBH and BNS mergers (Abbott et al. 2019). The input physics for the M30 model is summarised in Table 2 of Belczynski et al. (2020). The extension ‘B’ specifies the models that do not allow a common envelope with Hertzsprung gap donors. The M30B model used in this analysis was generated using the SFR specified in Eq. (16) in Belczynski et al. (2020), taking into account the evolution of metallicity with redshift using Eq. (18) in Belczynski et al. (2020). We refer to this SFR as the SFRL and show it in Fig. 1. We construct the following three mock populations for our analysis:

Mock 1. Population consisting of only Pop III compact binaries, constructed using the model FS1 and SFR1.

Mock 2. Population consisting of only Pop III compact binaries, constructed using the model FS1 and SFR2.

Mock 3. Population consisting of both Pop I+II and Pop III compact binaries, constructed with the model M30B for Pop I+II using the SFRL, and with the model FS1 for Pop III using SFR3. We used SFR3 here to have a non-negligible merger rate of Pop III binaries as compared to Pop I+II compact binaries.

The expected merger rates per year for Mock 1 and Mock 2 are 2429 and 1450, respectively, and as we generated 80 725 for Mock 1 and 68 500 compact binaries for Mock 2, these correspond to an observation time of ∼33 yr and ∼47 yr, respectively. In the case of Mock 3, we generated 162 106 compact binaries, which corresponds to ∼0.5 yr of observation for Mock 3.

3.1. Construction of the mock populations

For both M30B and FS1 compact binaries, we use the population available on the StarTrack¹ website. For the FS1 model, we only consider the following data for each binary: (i) the masses of the merging compact objects $m_{1,2}^i$, and (ii) the delay time between the formation of the binary at zero-age main sequence (ZAMS) and its coalescence, $t_{\text{del}}^i$. For compact binaries of the M30B model, we use the following data for each binary: (i) the masses of the merging compact objects $m_{1,2}^i$, (ii) the merger redshift, and (ii) the merger rate density for each binary in the observer frame of reference.

For a ZAMS binary with masses m_1, 2, delay times t_del, and metallicity Z formed at cosmic time t_ini = t_obs − t_del, corresponding to redshift z_ini = z(t_ini), the delay time t_del = t_evol + t_merg. The t_evol is the time of evolution from ZAMS to the formation of a compact binary system; t_merg is the time till the merger; and t_obs is the cosmic time at which the compact binary is observed to merge. Then, for a binary system ‘b’, the merger rate density per unit redshift as a function of redshift is given as:

$$\begin{array}{c}\hfill {\mathcal{R}}_b(z)=\frac{1}{(1+z)}\frac{\mathrm{d}V}{\mathrm{d}z}\left(\frac{\mathrm{SFR}(z_{i,\mathrm{ini}}),Z_b)}{M_{\mathrm{sim}}}\right),\end{array}$$(1)

where M_sim is the total mass of all stars that must accompany the stellar evolution, leading to formation of compact object binaries. These include the binaries and the single stars. M_sim for M30B and FS1 is 2.8 × 10⁸ M_⊙Belczynski et al. (2020) and 3.5 × 10⁹ M_⊙ (Belczynski et al. 2017), respectively.

Equation (1) gives the redshift dependence of the merger rate density of each binary, and so the probability density of a merger of a type b to happen in the Universe at redshift z is proportional to ℛ_b(z):

$$\begin{array}{c}\hfill P(b,z)\propto {\mathcal{R}}_b(z),\end{array}$$(2)

which is discrete in the index b and continuous in z. The probability distribution can be obtained by normalisation. We construct the mock populations using different SFRs, as mentioned above. To generate a mock population, we sample random binaries from the distribution given by Eq. (2) and to each compact binary system drawn in this manner we assign random values to the four angles: the right ascension α, the angle of declination δ, the polarisation angle ψ, and the inclination angle ι of the binary with respect to the direction of observation. The values of cos δ, α/π, cos ι, and ψ/π are chosen to be uncorrelated and distributed uniformly over the range [ − 1, 1]. In the following section, we describe the method to estimate the parameters of merging compact binaries detected with ET as a single instrument.

4. Estimation of the parameters of compact binaries with ET

For a merging compact binary of chirp mass ℳ, located at a luminosity distance D_L, the two polarisations of GWs at time t < t_c according to the theory of general relativity (Allen et al. 2012) are

$$\begin{array}{cc}& h_{+}(t)=-\frac{1+{cos}^2\iota }{2}\left(\frac{G\mathcal{M}}{c^2{D}_{\mathrm{L}}}\right){\left(\frac{t_c-t}{5G\mathcal{M}/c^3}\right)}^{-1/4}\hfill \\ \hfill & \times cos[2{\mathrm{\Phi }}_c+2\mathrm{\Phi }(t-t_c;M,\mu )],\hfill \end{array}$$(3a)

$$\begin{array}{cc}& h_{\times }(t)=-cos\iota \left(\frac{G\mathcal{M}}{c^2{D}_{\mathrm{L}}}\right){\left(\frac{t_c-t}{5G\mathcal{M}/c^3}\right)}^{-1/4}\hfill \\ \hfill & \times sin[2{\mathrm{\Phi }}_c+2\mathrm{\Phi }(t-t_c;M,\mu )],\hfill \end{array}$$(3b)

where c is the speed of light, G is the gravitational constant, μ is the reduced mass of the binary system, ι is the angle of inclination of the orbital plane of the binary system with respect to the observer, and Φ(t−t_c;M,μ) is the orbital phase of the binary system. For a binary system with component masses m₁ and m₂, the chirp mass is ℳ = (m₁m₂)^3/5/M^1/5, and the total mass is M = m₁ + m₂. The quantities t_c and Φ_c are the time and phase, respectively, at the termination of the waveform (Allen et al. 2012). The strain h(t) generated in a detector due to this waveform is

$$\begin{array}{c}\hfill h(t)=F_{+}{h}_{+}(t+t_c-t_0)+F_{\times }{h}_{\times }(t+t_c-t_0),\end{array}$$(4)

where t₀ is the time of coalescence in the detector frame, and so (t₀ − t_c) is the travel time from the source to the detector. F₊, F_× are the antenna response functions of one of the three detectors in ET. The strain can be rewritten as

$$\begin{array}{cc}& h(t)=-\left(\frac{G\mathcal{M}}{c^2}\right)\left(\frac{\mathrm{\Theta }}{4D_L}\right){\left(\frac{t_0-t}{5G\mathcal{M}/c^3}\right)}^{-1/4}\hfill \\ \hfill & \times cos[2{\mathrm{\Phi }}_0+2\mathrm{\Phi }(t-t_c;M,\mu )]\hfill \end{array}$$(5)

by substituting the values of the two polarisations from Eq. (3) into Eq. (4): The functions Θ and the phase Φ_c in Eq. (5) are functions of the antenna response functions F₊ and F_× and the angle of inclination ι, and are defined as

$$\begin{array}{c}\hfill \mathrm{\Theta }\equiv 2{[F_{+}^2{(1+{cos}^2\iota )}^2+4F_{\times }^2{cos}^2\iota ]}^{1/2},\end{array}$$(6)

such that 0 < Θ < 4 and

$$\begin{array}{c}\hfill 2{\mathrm{\Phi }}_0=2{\mathrm{\Phi }}_c-arctan\left(\frac{2F_{\times }cos\iota }{F_{+}(1+{cos}^2\iota )}\right).\end{array}$$(7)

Assuming that the three ET detectors have identical noise, the signal to noise ratio (S/N), ρ_j, for j = (1, 2, 3) in each of the three ET detectors is given as (Taylor & Gair 2012)

$$\begin{array}{c}\hfill {\rho }_j\approx 8{\mathrm{\Theta }}_j\frac{r_0}{D_{\mathrm{L}}}{\left(\frac{{\mathcal{M}}_z}{{\mathcal{M}}_{\mathrm{BNS}}}\right)}^{5/6}\sqrt{\zeta \left(f_{\max }\right)},\end{array}$$(8)

where the redshifted chirp mass ℳ_z = (1 + z)ℳ and the reference mass ℳ_BNS ≈ 1.218 M_⊙ is the chirp mass of an equal mass binary with components of 1.4 M_⊙ each. The function ζ(f_max) is defined as

$$\begin{array}{c}\hfill \zeta \left(f_{\max }\right)=\frac{1}{x_{7/3}}{\displaystyle {\int }_1^{2f_{max}}}\frac{\mathrm{d}f{\left(\pi M_{\odot }\right)}^2}{{\left(\pi f{M}_{\odot }\right)}^{7/3}{S}_h\left(f\right)},\end{array}$$(9)

where S_h(f) is the power spectral density (PSD) with the ET-D noise curve (Hild et al. 2011) for the ET-D configuration and,

$$\begin{array}{cc}& x_{7/3}={\displaystyle {\int }_1^{\infty }}\frac{\mathrm{d}f{\left(\pi M_{\odot }\right)}^2}{{\left(\pi f{M}_{\odot }\right)}^{7/3}{S}_h\left(f\right)},\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill & r_0^2=\frac{5}{192\pi }{\left(\frac{3G}{20}\right)}^{5/3}{x}_{7/3}\frac{M_{\odot }^2}{c^3},\hfill \end{array}$$(11)

$$\begin{array}{cc}\hfill & f_{\max }=785\left(\frac{M_{\mathrm{BNS}}}{M(1+z)}\right)\mathrm{Hz},\hfill \end{array}$$(12)

where r₀ is the characteristic distance sensitivity and f_max is the frequency at the end of the inspiral phase. For the ET triangular configuration consisting of three detectors, the combined effective S/N is given as

$$\begin{array}{c}\hfill {\rho }_{\mathrm{eff}}=8{\mathrm{\Theta }}_{\mathrm{eff}}\frac{r_0}{D_{\mathrm{L}}}{\left(\frac{{\mathcal{M}}_z}{1.2M_{\odot }}\right)}^{5/6}\sqrt{\zeta \left(f_{\max }\right)},\end{array}$$(13)

where the effective antenna response function Θ_eff is

$$\begin{array}{c}\hfill {\mathrm{\Theta }}_{\mathrm{eff}}={({\mathrm{\Theta }}_1^2+{\mathrm{\Theta }}_2^2+{\mathrm{\Theta }}_3^2)}^{1/2}.\end{array}$$(14)

4.1. Summary of the algorithm used in SB2

In SB2, we estimate the parameters of the merging compact binaries, taking into account the change in the antenna pattern with the rotation of the Earth in order to analyse the long-duration GW signals. We assume the location of the ET detector to be the Virgo site, and we analyse the signal every 5 minutes. The observed GW frequency $f_{\mathrm{gw}}^{\mathrm{obs}}$ is calculated using Eq. (4.195) in Maggiore (2007):

$$\begin{array}{c}\hfill f_{\mathrm{gw}}^{\mathrm{obs}}=\frac{1}{\pi }{\left(\frac{5}{256}\frac{1}{{\tau }_{\mathrm{obs}}}\right)}^{3/8}{\left(\frac{G{\mathcal{M}}_z}{c^3}\right)}^{-5/8},\end{array}$$(15)

where τ_obs is the time to coalescence, which is measured in the observer’s frame. The minimum frequency given the detection sensitivity of the detector and the frequency f_max determines the limit on τ_obs in the detection band. For τ_i − 1 and τ_i, which are the initial and the final values of τ_obs, respectively, for the ith segment, the corresponding values f_i − 1, f_i of $f_{\mathrm{gw}}^{\mathrm{obs}}$ will be

$$\begin{array}{c}\hfill f_{i-1}=\frac{1}{\pi }{\left(\frac{5}{256}\frac{1}{{\tau }_{i-1}}\right)}^{3/8}{\left(\frac{G{\mathcal{M}}_z}{c^3}\right)}^{-5/8},\end{array}$$(16)

and

$$\begin{array}{c}\hfill f_i=\frac{1}{\pi }{\left(\frac{5}{256}\frac{1}{{\tau }_i}\right)}^{3/8}{\left(\frac{G{\mathcal{M}}_z}{c^3}\right)}^{-5/8}.\end{array}$$(17)

The S/N for the ith segment in the jth detector can be written using Eq. (8), as

$$\begin{array}{c}\hfill {\rho }_j^i\approx 8{\mathrm{\Theta }}_j^i\frac{r_0}{D_{\mathrm{L}}}{\left(\frac{{\mathcal{M}}_z}{{\mathcal{M}}_{\mathrm{BNS}}}\right)}^{5/6}\sqrt{{\zeta }^i(f_{i-1},f_i),}\end{array}$$(18)

where

$$\begin{array}{c}\hfill {\zeta }^i(f_{i-1},f_i)=\frac{1}{x_{7/3}}{\displaystyle {\int }_{f_{i-1}}^{f_i}}\frac{\mathrm{d}f{\left(\pi M_{\odot }\right)}^2}{{\left(\pi f{M}_{\odot }\right)}^{7/3}{S}_h\left(f\right)},\end{array}$$(19)

and

$$\begin{array}{c}\hfill {\mathrm{\Theta }}_j^i\equiv 2{[{(F_{+}^i)}^2{(1+{cos}^2\iota )}^2+4{(F_{\times }^i)}^2{cos}^2\iota ]}_j^{1/2},\end{array}$$(20)

where $F_{+}^i$ and $F_{\times }^i$ are the antenna response functions for the jth detector in the ith segment. The effective S/N for the ith segment is

$$\begin{array}{c}\hfill {\rho }_{\mathrm{eff}}^i=8{\mathrm{\Theta }}_{\mathrm{eff}}^i\frac{r_0}{D_{\mathrm{L}}}{\left(\frac{{\mathcal{M}}_z}{1.2M_{\odot }}\right)}^{5/6}\sqrt{{\zeta }^i(f_{i-1},f_i),}\end{array}$$(21)

where

$$\begin{array}{c}\hfill {({\rho }_{\mathrm{eff}}^i)}^2={({\rho }_1^i)}^2+{({\rho }_2^i)}^2+({\rho }_3^i){,}^2\end{array}$$(22)

and the function ${\mathrm{\Theta }}_{\text{eff}}^i$ is

$$\begin{array}{c}\hfill {({\mathrm{\Theta }}_{\mathrm{eff}}^i)}^2={({\mathrm{\Theta }}_1^i)}^2+{({\mathrm{\Theta }}_2^i)}^2+{({\mathrm{\Theta }}_3^i)}^2.\end{array}$$(23)

We assume that the observables, in the case of the detection of a coalescing binary system, are: (a) the three S/Ns ${\rho }_j^i$ for the ith segment of the signal, (b) the phase ${\mathrm{\Phi }}_{o,j}^i$ for j = (1, 2, 3) for each of the three ET detectors in the ith segment of the signal, (c) the GW frequency at the start and end of each segment of the detected signal, (d) the redshifted chirp mass ℳ_z, and (e) the frequency at the end of the inspiral, corresponding to the innermost stable circular orbit, f_max. We assume the measurement errors on the S/Ns to be Gaussian, such that the standard deviations for ${\rho }_j^i$ and ${\mathrm{\Phi }}_j^i$ are σ_ρ = 1 and σ_Φ = π/ρ, respectively. This is a conservative assumption as compared to the errors on the S/Ns of the GW detections mentioned in GWTC-2 and GWTC-3 (Abbott et al. 2021, 2023).

One of the main building blocks of the algorithm in SB2 is that we use the ratios of S/N in each segment in order to constrain ${\mathrm{\Theta }}_{\text{eff}}^i$ (see Sect. IV of SB2):

$$\begin{array}{c}\hfill {\rho }_{21}^i=\frac{{\mathrm{\Theta }}_2^i}{{\mathrm{\Theta }}_1^i}\equiv {\mathrm{\Theta }}_{21}^i\mathrm{and}{\rho }_{31}^i=\frac{{\mathrm{\Theta }}_3^i}{{\mathrm{\Theta }}_1^i}\equiv {\mathrm{\Theta }}_{31}^i,\end{array}$$(24)

where ${\rho }_{21}^i\equiv {\rho }_2^i/{\rho }_1^i$ and ${\rho }_{31}^i\equiv {\rho }_3^i/{\rho }_1^i$ in the ith segment. We then proceed to constrain the source-dependent quantity Λ (defined in Eq. (37) of SB2) for each segment in order to constrain the binary parameters such as chirp mass, total mass, mass ratio, redshift, and luminosity distance:

$$\begin{array}{c}\hfill \mathrm{\Lambda }\equiv {\left(\frac{8r_0}{D_{\mathrm{L}}}{\left(\frac{{\mathcal{M}}_z}{{\mathcal{M}}_{\mathrm{BNS}}}\right)}^{5/6}\right)}^{-1}\approx \frac{{\mathrm{\Theta }}_{\mathrm{eff}}^i\sqrt{{\zeta }^i(f_{i-1},f_i)}}{{\rho }_{\mathrm{eff}}^i}.\end{array}$$(25)

With this algorithm, we find that the chirp masses are overestimated, while the redshift is underestimated. This was also seen in the estimates of the binary parameters for Pop III binaries carried out in S22. The origin of this ‘bias’ is the probability distribution of Θ (see Fig. 13 of SB2) given the antenna pattern functions of the triangular configuration of ET (see Appendix of SB2 for detail).

4.2. Modification in the algorithm used in SB2

By comparing both sides of Eq. (25), we find a function $\mathcal{F}({\rho }_{21}^i,{\rho }_{31}^i)$, such that

$$\begin{array}{c}\hfill \mathcal{F}({\rho }_{21}^i,{\rho }_{31}^i)=\frac{{\mathrm{\Lambda }}_{\mathrm{s}}}{{\mathrm{\Lambda }}_{\mathrm{med}}},\end{array}$$(26)

where the subscript ‘s’ denotes the true value of Λ from the actual source parameters,

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_{\mathrm{s}}={\left(\frac{8r_0}{D_{\mathrm{L}}}{\left(\frac{{\mathcal{M}}_z}{{\mathcal{M}}_{\mathrm{BNS}}}\right)}^{5/6}\right)}^{-1},\end{array}$$(27)

and Λ_med is the median of the estimated posterior distribution of Λ in the ith segment of the signal:

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_{\mathrm{med}}={\left(\frac{{\mathrm{\Theta }}_{\mathrm{eff}}^i\sqrt{{\zeta }^i(f_{i-1},f_i)}}{{\rho }_{\mathrm{eff}}^i}\right)}_{\mathrm{median}}.\end{array}$$(28)

We approximate $\mathcal{F}({\rho }_{21}^i,{\rho }_{31}^i)$ to a fit of the form

$$\begin{array}{c}\hfill \mathcal{F}({\rho }_{21}^i,{\rho }_{31}^i)\approx 1+\frac{\mathcal{R}({\rho }_{21}^i,{\rho }_{31}^i,1,1,0.022)}{1.36\times {10}^6}+\frac{\mathcal{R}({\rho }_{21}^i,{\rho }_{31}^i,1,1,0.7),}{8}\end{array}$$(29)

where ℛ is the Marr or Mexican hat function, and is given as

$$\begin{array}{cc}& \mathcal{R}(x,y,x_0,y_0,\sigma )=\frac{1}{\pi {\sigma }^4}(1-\frac{1}{2}\left(\frac{{(x-x_0)}^2+{(y-y_0)}^2}{{\sigma }^2}\right)).\hfill \\ \hfill & \times e^{-\frac{{(x-x_0)}^2+{(y-y_0)}^2}{2{\sigma }^2}}\hfill \end{array}$$(30)

A more detailed description of this fit is described in Sect. A. We include this function $\mathcal{F}({\rho }_{21}^i,{\rho }_{31}^i)$ in Eq. (39) in SB2 as prior information, and continue the rest of the analysis to estimate the parameters of a given binary as described in SB2. In this work, we fix the detection at a threshold value of accumulated effective S/N ρ_eff > 8, and the S/N for ith segment in the jth detector ${\rho }_j^i$ > 3 in at least one segment for j = (1, 2, 3), corresponding to the three ET detectors comprising the single ET. We then use this improved algorithm to estimate the chirp mass, total mass, and redshift of a detected compact binary system from the three mock populations.

We consider three populations of compact binaries originating from Pop III stars and merging within a Hubble time based on three different SFRs. The construction of all three populations is described in Sect. 3.1. For a given population of compact binaries, the true values of the chirp mass, total mass, and redshift of these ‘sources’ are represented as ℳ_s, mock, M_s, mock, and z_s, mock, respectively. A binary source is considered as ‘detected’ if it crosses a detection threshold set on the S/N as mentioned in Sect. 4.2. The chirp mass, total mass, and redshift of these detected sources are denoted ℳ_s, det, M_s, det, and z_s, det, respectively. The posterior probability distribution for the chirp mass and redshift for each of these detected sources is estimated using the algorithm described above (see Fig. 5 in SB2 for an example). The median values of these estimated posterior probability distributions of chirp mass, total mass, and redshift for each detected compact binary source are represented as ℳ_med, det, M_med, det, and z_med, det, respectively.

4.3. Chirp mass and redshift estimates

The density distribution of the estimated median values with respect to the true values of the parameters for the three populations estimated using the algorithm described above are shown in Fig. 2. This figure shows the estimated median values (ℳ_med, det, z_med, det) with respect to the actual values of the parameters (ℳ_s, det, z_s, det) for each detected compact binary source in three mock populations. The blue contour encloses the 90% probability region of all the detected sources. The plots in the top panel, namely Figs. 2a and 2b, can be compared with Figs. 3c,d of S22 in order to see the effect of using $\mathcal{F}({\rho }_{21}^i,{\rho }_{31}^i)$ as an additional prior on Λ. The estimates of the chirp masses as compared to their true values for the compact binaries in all three populations are shown in Figs. 2a, c, and e. We see that chirp mass is overestimated for half of the population and underestimate for the other half. Using the current algorithm, the median values of the redshift of 90% of the detected sources, shown in Figs. 2b, d, and f, are in agreement with the true values of redshift for Mock 1 and Mock 3, while the redshifts are slightly underestimated for Mock 2. The reason for this it that including $\mathcal{F}({\rho }_{21}^i,{\rho }_{31}^i)$ as an additional prior only improves the parameter estimates for those binary sources for which 0.95 ≲ ρ_ij ≲ 1.05, where ρ_ij is the ratio of S/N generated in the ith and jth detector of ET. This is explained in more detail in Appendices A and B. We quantify the error in the estimate of the redshift in the following section by calculating the merger rate density for the mock populations.

Fig. 2. Estimation of parameters. Density distribution of the median of the estimated posterior with respect to the true values of the parameters: Chirp mass ℳ (left) and redshift z (right) of the detected compact binary sources. Top panel: Mock 1. Middle panel: Mock 2. Bottom panel: Mock 3. The magenta line is a reference for equal values of true and estimated parameters. The blue contour encloses the 90% probability region.

5. Estimation of the merger rate density

In this section, we describe a population-independent method of estimation of the merger rate density as a function of redshift. Using the parameters estimated in the previous section, we calculate the merger rate density for the sources we detect, for a given detection threshold. Then, using the probability distribution of the population parameters of these detected sources, we asses the detection efficiency as a function of redshift. Taking into account this detection efficiency, we then reconstruct the merger rate density for the population of coalescing compact binaries.

5.1. Merger rate density

For a given population of compact binaries, we simulate N_mock number of binaries, of which N_det are detected based on the chosen detection threshold. The expected time taken for these binaries in the mock population to merge is

$$\begin{array}{c}\hfill T_{\mathrm{mock}}=\frac{N_{\mathrm{mock}}}{N_{\mathrm{yr}}}\mathrm{yr},\end{array}$$(31)

where N_yr is the number of mergers per year calculated by integrating the merger rate density given by Eq. (1). The merger rate density R_mer for a given population is calculated for the time T_mock:

$$\begin{array}{c}\hfill R_{\mathrm{mer}}(z_i,z_{i+1})=\frac{1+z_{i+1}}{{\int }_{z_i}^{z_{i+1}}\frac{\mathrm{d}V}{\mathrm{d}z}\mathrm{d}z}\left(\frac{N_{(z_i,z_{i+1})}}{T_{\mathrm{mock}}}\right),\end{array}$$(32)

where N_{(zi, zi + 1)} is the number of mergers in a redshift bin [z_i : z_i + 1].

From Eq. (32), we obtain three sets of merger rate densities for a given population: (i) R_mer(z_s, mock), (ii) R_mer(z_s, det), and (iii) R_mer(z_med, det). These merger rate densities are shown in Figs. 3a, d, and g for Mock 1, Mock 2, and Mock 3, respectively. The shaded region represents the Poisson error. For Mock 1 and Mock 2, we estimated the Poisson error for yearly data sets while for Mock 3 we estimated the Poisson error for monthly data sets.

Fig. 3. Estimate of the merger rate density. Left: merger rate densities for Mock 1, Mock 2, and Mock 3, respectively. The shaded region represents Poisson error. Middle: Relative merger rate density. Right: cumulative probability distribution of the redshifts for Mock 1, Mock 2, and Mock 3. A description of the acronyms is given in Table D.1.

It can be seen that, for all three populations, R_mer(z_med, det)≈R_mer(z_s, det), but R_mer(z_med, det)≪R_mer(z_s, mock). This is because most of the merging compact binary systems at higher redshifts do not cross the detection threshold we have chosen. In order to further compare the merger rate densities, we calculated the relative merger rate densities $\frac{R_{\mathrm{mer}}(z_{\mathrm{med},\det })}{R_{\mathrm{mer}}(z_{\mathrm{s},\mathrm{mock}})}$ and $\frac{R_{\mathrm{mer}}(z_{\mathrm{med},\det })}{R_{\mathrm{mer}}(z_{\mathrm{s},\det })}$. These are shown in Figs. 3b, e, and h. This calculation shows that R_mer(z_med, det)≈R_mer(z_s, det) for Mock 3, while the in cases of Mock 1 and Mock 2, R_mer(z_med, det) is slightly larger than R_mer(z_s, det) in the specific redshift ranges by a factor of ∼3 and ∼4, respectively.

Figures 3c, f, and i show the cumulative probability distribution of the redshifts for Mock 1, Mock 2, and Mock 3 and we see that < 1% of the sources in each of the three mock populations are in the redshift range where the relative merger rate densities are > 1.5. The reason for these few spikes is further explained in Appendix B. As R_mer(z_med, det)≪R_mer(z_s, mock), we calculated the detection efficiency in order to reconstruct the actual merger rate density from the R_mer(z_med, det).

5.2. Detection efficiency

Detection of a coalescing compact binary with a gravitation wave detector is ideally defined by the threshold we choose to set for the S/N of such an event. The higher the S/N threshold value set for detection, the more sources will be left undetected. It is important to note here that, for our analysis, we use the design sensitivity noise curve of ET-D to estimate the S/N and set the threshold for detection, whereas in the case of detections with real noise, one has to take into account the glitches, which will further introduce a selection bias in the detection of events. We neglect the presence of glitches in our analysis. In order to get an accurate estimate of the merger rate density for a population of compact binaries, it is necessary to gauge the number of compact binaries that are not detected as a function of redshift. To this end, we make an assumption that the detected population of the compact binaries truly represents the redshift and chirp mass distribution of the whole population. For a given population, we generate a secondary mock population – denoted with subscript ‘sec’ – from the detected sources, assuming that the distributions of the chirp mass and redshift are proportional to ℳ_med, det and z_med, det, that is, p(ℳ_sec)∝p(ℳ_med, det) and p(z_sec)∝p(z_med, det). We assume that the mass ratio q_sec is uniformly distributed in the range [0,1] with a constraint on total mass M_sec such that (M_med, det)_min ≤ M_sec ≤ (M_med, det)_max, where M_sec is defined as:

$$\begin{array}{c}\hfill M_{\sec }={\mathcal{M}}_{\sec }{\left[\frac{q_{\sec }}{{(1+q_{\sec })}^2}\right]}^{-3/5}.\end{array}$$(33)

To each compact binary of this generated secondary population, we assign random values to the four angular parameters: the right ascension α, the angle of declination δ, the polarisation angle ψ, and the inclination angle ι of the binary with respect to the direction of observation. The values of cos δ, α/π, cos ι, and ψ/π are chosen to be uncorrelated and distributed uniformly over the range [ − 1, 1]. We now use this secondary population to estimate the detection efficiency. Given the detection threshold we chose, the detection efficiency 𝒟 as a function of redshift is defined as

$$\begin{array}{c}\hfill \mathcal{D}(z_i,z_{i+1})={\left[\frac{N_{\sec ,\det }}{N_{\sec }}\right]}_{(z_i,z_{i+1})},\end{array}$$(34)

where [N_sec]_{(zi, zi + 1)} is the number of mergers in the secondary mock population in the redshift bin (z_i, z_i + 1) and [N_sec, det]_{(zi, zi + 1)} is the mergers in this bin that crossed the detection threshold.

The cumulative probability distributions of the secondary mock populations for each of the three mock populations, that is, Mock 1, Mock 2, and Mock 3, are shown in the left panel of Fig. 4 in red, while the cumulative probability distributions of the detected sources from these secondary mock populations are shown in blue. For a given mock population, we use the ratio of the gradient of these curves – taking into account the number of sources in each set of populations – to quantify the detection efficiency.

Fig. 4. Reconstruction of merger rate density. Left: cumulative probability distribution of the redshift of the secondary population (red) and of the detected sources from these secondary mock populations (blue) for Mock 1, Mock 2, and Mock 3. The number in parentheses is the number of sources in each population. Right: reconstructed merger rate density for Mock 1, Mock 2, and Mock 3. The shaded region represents the Poisson error.

5.3. Reconstructed merger rate density

Now we can reconstruct the merger rate density taking into account the detection efficiency calculated in the previous section. The reconstructed merger rate density R_mer, recon is then given as

$$\begin{array}{c}\hfill R_{\mathrm{mer},\mathrm{recon}}(z_i,z_{i+1})={\left[\frac{R_{\mathrm{mer}}(z_{\mathrm{med},\det })}{\mathcal{D}}\right]}_{(z_i,z_{i+1})}.\end{array}$$(35)

The reconstructed merger rate density is calculated for the three mock populations using Eq. (34) in Eq. (35). The reconstructed merger rate densities for Mock 1, Mock 2, and Mock 3 are shown in the right panel Fig. 4. The red lines are the true merger rate density, while the green lines show the merger rate density calculated using Eq. (35). The shaded region represents the Poisson error. It can be seen that the merger rate density is reconstructed accurately up to redshift z ∼ 15 for the Mock 1 population, and up to redshift z ∼ 14 and Mock 2.

In the case of Mock 3 binaries, the reconstructed merger rate density at z ∼ 2 is a factor of ∼1.3 smaller that the true merger rate density. This is so because we calculate the detection efficiency by generating a secondary mock population, assuming that the probability distribution of the chirp mass and redshift is proportional to that of ℳ_med, det and z_med, det, respectively. However, in the case of Mock 3 binaries, the merger rate density at z ≲ 2 is dominated by low-mass binaries (see Fig. 12 in Belczynski et al. 2020), and as we do not detect the bulk of these low-mass binaries given the chosen detection threshold (see Fig. in S22), the mass distribution for these objects is not truly represented in the secondary population constructed for Mock 3. The reconstructed merger rate density for Mock 3 is also lower than the true value for z > 8. As seen in Fig. 2e, Pop III binaries constitute a small percentage of Mock 3 binaries, and so they are also under-represented in the secondary population. In order to have an accurate estimate of the detection efficiency in a mixed population of binaries where the merger rate densities of the two individual populations differ by a large factor, it is essential to have a larger observational data set in order for the underlying populations to be accurately represented.

6. Star formation rates

We now proceed to estimate the SFR using our estimate of the merger rate densities. We calculate the SFR using R_mer, recon, assuming three different time-delay t_del distributions: (i) a broken power law, (ii) t⁻¹, and (iii) t⁻², where t_del, min < t < t_H Gyr. The broken power-law distribution is assumed to be:

$$\begin{array}{c}\hfill p(t_{\mathrm{del}})\propto \{\begin{array}{cc}{t}^{-3}\hfill & \mathrm{if}t<1\mathrm{Gyr}\hfill \\ t^{-2}\hfill & \mathrm{if}t⩾1\mathrm{Gyr}.\hfill \end{array}\end{array}$$(36)

The broken power law time-delay distribution is motivated by the fit to the actual time-delay distribution as shown in Fig. 11 of Belczynski et al. (2017) for FS1 compact binaries. The time-delay distribution t⁻¹ is based on the best fit for the Pop I+II compact binaries (M30B time-delay data are available on the StarTrack² website). The SFR estimate at a given z is the sum over contributions from each binary ‘b’, averaged over its delay time:

$$\begin{array}{c}\hfill \mathrm{SFR}(z)=M_{\mathrm{sim}}{\displaystyle \sum _b\left({\int }_{t_{\mathrm{del},\min }}^{t_{\mathrm{H}}}{R}_{\mathrm{mer},\mathrm{recon}}^b\frac{\mathrm{d}{z}_{\mathrm{mer}}}{\mathrm{d}{t}_{\mathrm{del}}}p(t_{\mathrm{del}})\mathrm{d}{t}_{\mathrm{del}}\right),}\end{array}$$(37)

for $z_{\text{ini}}^b\equiv z(t_{\text{mer}}^b-t_{\text{del}}^b)$. M_sim is a constant, and is dependent on the formation scenario. Thus, M_sim provides a normalisation for a formation channel. We assume two different values for minimum time delay: t_del, min = 0.03 and 0.05 Gyr. We denote the Hubble time t_H. In this analysis, we normalise the SFR to one, because the goal is to estimate the functional form of the SFR with redshift. We therefore remove any dependence on the formation channel by removing the dependence on M_sim in the case of Mock 1 and Mock 2 populations, where all the binaries are assumed to be from a single population. In the case of Mock 3, where we have a mixed population, our algorithm estimates a merger-rate-weighted SFR as a function of redshift instead of the true sum of two individual SFRs. We denote this merger-rate-weighted SFR as SFR′.

It should also be noted that while the three mock populations constructed in this analysis are built with the time-delay distribution that was the output of the StarTrack code and so encode the uncertainties for the formation channel parameters, we reconstruct the SFR assuming the three different time-delay t_del distributions described above. We compare the estimated values of the SFRs with the true value of the SFR by calculating the relative entropy (RE), otherwise known as the Kullback–Leibler divergence (D_KL; Cover 1999; Kullback & Leibler 1951; Shlens 2014). RE, or D_KL, quantifies how close a given probability distribution p = p_i is to a given model distribution q = q_i. It can also be said that the D_KL(p||q) is a measure of the inaccuracy of the assumption that the distribution is q when the true distribution is p:

$$\begin{array}{c}\hfill \mathrm{RE}\equiv D_{\mathrm{KL}}(p||q)={\displaystyle \sum _i}{p}_{i}log\left(\frac{p_i}{q_i}\right).\end{array}$$(38)

RE is always non-negative, and it is zero only under the condition that p = q (Cover 1999). We assume that for RE > 1, the distributions cannot be considered to be similar.

To understand the errors in the reconstructed SFR due to the different time-delay distributions and due to the error in the estimated merger rate density, we proceed as follows. As a first step, we estimate the SFR assuming that merger rate densities for each of the three populations R_mer(z_s, mock) are known with high accuracy, so that we can see the effect of only the variation of true time-delay distribution – as compared to the assumed time-delay distributions – on the reconstructed SFR. Then, as a second step, we reconstruct the SFR with R_mer, recon for each of the three mock populations in order to estimate the errors on the SFR due to both the error on the estimated merger rate density and the variation of the true time-delay distribution in comparison to the assumed time-delay distribution.

Figure 5 shows the SFR reconstructed for Mock 1 (top), Mock 2 (middle), and Mock 3 (bottom), assuming an accurate merger rate density with redshift R_mer(z_s, mock). In each of these three figures, the left panel shows the reconstruction assuming that t_del, min = 0.03 Gyr, and the right panel shows the results for t_del, min = 0.05 Gyr. The top three panels in all three plots show the RE for a redshift bin of Δz = 0.5 in width.

Fig. 5. Reconstructed SFRs for the three populations, assuming that the merger rate density is highly accurate. The maximum value of the SFR has been normalised to 1. The left and right panels show the reconstruction assuming that tdel, min = 0.03 and 0.05 Gyr, respectively. The top three panels in all three plots show RE or the KL divergence for a redshift bin of Δz = 0.5 in width. The red curve denotes the true SFR. We note that for the Mock 3 population, we estimate the merger-rate-weighted SFR.

For each of the three mock populations, the reconstruction of SFR as a function of redshift is independent of the assumed time-delay distribution throughout almost the entire redshift range. The comparison of the reconstructed SFR with the true SFR using KL divergence shows that RE < 1 from nearly the termination redshift up to z ∼ 14. In a few redshift bins, RE > 1 due to the fact that the assumed time-delay distributions are not exact representations of the true time-delay distributions of these mock populations. It should be noted that the assumed time-delay distributions in this analysis do not correctly model the longer time-delay distribution t_del > 8 Gyr. This is evident in the case of Mock 2 in Fig. 5b. While it is crucial to have a better model for long time-delay distribution in order to accurately estimate the termination redshift of the SFR, accurate information about the minimum time delay t_del, min is essential in order to estimate the SFR at redshifts beyond z ∼ 14. For Mock 1 and Mock 2, the true values of the termination redshift (which we define as the redshift where the SFR is 1% of its peak value) are ∼3.2 and ∼5.3. As seen from Fig. 5, the errors on the estimates of the termination redshifts due to incorrect modelling of the tail of the time-delay distribution are Δz/z ≲ 7% and Δz/z ≲ 32% for Mock 1 and Mock 2, respectively.

Regarding the reconstruction of the SFR for Mock 3 shown in Fig. 5c, an important point to note is that the estimated SFR is a rate-weighted SFR, and so SFR′≠SFRL + SFR3. The reconstructed SFR clearly shows the presence of two different populations, with one formation peaking at z ∼ 2 and another peaking at z ∼ 10. As in the case of Mock 1 and Mock 2, the RE values in the case of Mock 3 show that the reconstruction is independent of the time-delay distributions up to z ∼ 14 and the accuracy of the reconstruction of SFR strongly depends on the time-delay distribution only at higher redshifts of z ≳ 14. We further verify this conclusion by assuming a few more extreme time-delay distributions. A discussion about the reconstruction of the SFR with these extreme time-delay distributions is presented in Appendix C.

We now proceed to reconstruct the SFR with R_mer, recon for each of the three mock populations. Figure 6 shows the estimated SFRs for Mock 1, Mock 2, and Mock 3. The true SFR is shown in red for each of these three populations. The left panel in each figure shows the reconstructed SFR assuming that t_del, min = 0.03 Gyr, and the right panel shows the results for t_del, min = 0.05 Gyr. The shaded region represents Poisson error. The RE comparing the estimated values of the SFRs with the true value of the SFR in redshift bins of Δz = 0.5 in width is shown in the top three panels for each of the three populations.

Fig. 6. Reconstructed SFRs for the three populations with the reconstructed merger rate density Rmer, recon. The maximum value of the SFR has been normalised to 1. The left and right panels show the reconstruction assuming that tdel, min = 0.03 and 0.05 Gyr, respectively. The top three panels in all three plots show RE or the KL divergence for a redshift bin of Δz = 0.5 in width. The red curve denotes the true SFR. We note that for the Mock 3 population, we estimate the merger-rate-weighted SFR.

As before, we see that for each of the three mock populations, the reconstruction of the SFR as a function of redshift is independent of the assumed time-delay distribution and the comparison of the reconstructed SFR with the true SFR using KL divergence shows that RE < 1 from nearly the termination redshift up to z ∼ 14. The combined error on the termination redshift is due to the error from the improper modelling of the long time-delay distributions and the error on the estimates of the merger rate density as a function of redshift.

The top three panels in Figs. 6a, b, and c show that RE < 1 for 4 ≲ z ≲ 14 for Mock 1, that RE < 1 for 6 ≲ z ≲ 14 for Mock 2, and that RE < 1 for 4 ≲ z ≲ 14 for Mock 3, respectively, irrespective of the assumed time-delay distribution. The errors on the estimate of the termination redshift for Mock 1 and Mock 2 are Δz/z ≲ 22% and Δz/z ≲ 37%, respectively.

For Mock 3, which is a mixed population of Pop I+II+III, we show the reconstructed SFR in Fig. 6c. As seen in Fig. 4f, the R_mer, recon < R_mer(z_s, mock) for z ∼ 2 and z > 8 because of the underestimation of the detection efficiency. As a result of this, RE < 1 up to z ≲ 14, except in the range 4 ≲ z ≲ 8. The reconstructed SFR for this mixed population clearly shows the presence of the peaks of the two different populations, with one formation peaking at z ∼ 2 and another peaking at z ∼ 10. This estimate can be further improved with a larger observational data set, because this would provide a better representation of the underlying populations.

As we only use the inspiral part in this analysis to estimate the parameters of the compact binaries, as done previously in SB1, SB2, and S22, the S/N generated in the ET detectors is underestimated. The estimates on the parameters are therefore conservative, and so the errors estimated on the reconstructed SFR are the upper bounds.

7. Conclusion

In this paper, we used an updated version of the SB2 algorithm to estimate parameters such as chirp mass, redshift, total mass, and mass ratio for compact binaries. We reconstructed the merger rate density and SFR for three mock population models using single ET and assuming a triangular configuration and ET-D design sensitivity. We constructed the mock populations for compact binaries originating in stars from Population (Pop) I+II and Pop III, assuming different SFRs and realistic time-delay distributions.

For a given population, as a first step, we estimated the chirp mass, redshift, and total mass of each detected compact binary. We then estimated the detection efficiency for each population and thus reconstructed the merger rate density taking into account the fraction of binaries that do not cross the detection threshold. We then reconstructed the SFRs using the estimated merger rate density assuming three different functional forms for t_del and two different values of t_del, min. The variable names mentioned in the text are summarised in Table D.1.

For Mock 1 and Mock 2, the true values of termination redshifts are z ∼ 3.2 and z ∼ 5.3, and the errors on the estimation of these termination redshifts due to incorrect modelling of tail of the time-delay distribution are Δz/z ≲ 7% and Δz/z ≲ 32% for Mock 1 and Mock 2, respectively. Taking into account the error on the estimates of the merger rate density in addition to the incorrect modelling of the tail of time-delay distribution, we estimate that the errors on the termination redshift are Δz/z ≲ 22% and Δz/z ≲ 37% for Mock 1 and Mock 2, respectively. We conclude that the farther the true termination redshift is, the larger the error on the estimate will be given the inaccurate modelling of the tail of the time-delay distribution.

For Mock 3, which is a mixed population of Pop I+II+III, the reconstructed merger rate density at z ∼ 2 is a factor of ∼1.3 smaller that the true merger rate density, because the detection efficiency calculated by generating a secondary mock population assumes that the probability distributions of chirp mass and redshift of the underlying population are represented by the distributions of the chirp mass and redshift of the detected population, respectively. However, in the case of Mock 3 binaries, the merger rate density at z ≲ 2 is dominated by low-mass binaries, and given that we do not detect the bulk of these objects, with the chosen detection threshold, the mass distribution for these low-mass binaries is not truly represented in the detected population for Mock 3. The reconstructed merger rate density for Mock 3 is also lower than the true value for z > 8. As Pop III binaries constitute a small percentage of Mock 3 binaries, they are under-represented in the secondary population. This estimate can be further improved with a larger observational data set, because it will provide a better representation of the underlying populations and thus improve our estimate of the detection efficiency.

In conclusion, we provide a method to reconstruct the functional form of the SFR for populations of compact binaries with ET. The SFR as a function of redshift is accurately reconstructed up to redshift z ∼ 14. For all three of our mock populations, we show that the reconstruction of SFR is independent of the time-delay distributions up to z ∼ 14. The accuracy of the reconstruction of the SFR beyond z ∼ 14 strongly depends on the minimum value of the time delay t_del, min. The assumed time-delay distributions in this analysis do not correctly model the longer time-delay distribution t_del > 8 Gyr. While we accurately reconstruct the SFR as a function of redshift for the bulk of each mock population, a better model for the long time-delay distribution is needed in order to accurately estimate the termination redshift of the SFR. We therefore constrained the peak of the SFR as a function of redshift, and show that ET as a single instrument can distinguish the termination redshifts of different SFRs if they have a true separation of at least Δz ∼ 2.

While we used only one population-evolution model for Pop III (FS1) and Pop I+II (M30B), using different population models will not effect the recovery of the peak of the SFRs unless the peak is beyond z ∼ 14; beyond this redshift, we need accurate time-delay distributions to estimate the SFR. Using different population models will also not effect the recovery of the termination redshift because the error comes only from the inaccurate modelling of the tail of the time-delay distribution. For any given population, if the SFR terminates at some redshift of z ∼ 6 for example and if we still detect the binaries from this population at z ∼ 0 − 1, this means these are the binaries with long time delays. The farther the true termination redshift is, the greater the probability that the binaries we detect today will be those with long time delays. This is the reason for the larger error on the estimate of termination redshift for Mock 2 as compared to Mock 1.

Acknowledgments

We thank the anonymous referee for very helpful comments and suggestions. NS is supported by the “Agence Nationale de la Recherche”, grant n. ANR-19-CE31-0005-01 (PI: F. Calore) and is thankful to Astronomical Observatory, University of Warsaw for providing access to computing resources. TB and MC are supported by the grant “AstroCeNT: Particle Astrophysics Science and Technology Centre” (MAB/2018/7) carried out within the International Research Agendas programme of the Foundation for Polish Science (FNP) financed by the European Union under the European Regional Development Fund. TB acknowledges support of NCN through the Harmonia grant UMO-2017/26/M/ST9/00978. This document has been assigned LAPTh document number LAPTH-016/23.

References

Appendix A: Construction of ℱ

In SB2, we found that the estimates for the parameters of the compact binaries were biased. The chirp masses were overestimated, while the redshift was underestimated. The origin of this bias is described in the Appendix of SB2 (see Figure 13 of SB2). In order to find a fitting function, we proceeded as follows. We constructed a mock population of 1500 low-mass compact binaries according to the description given in Sec V of SB2.

We chose to generate this mock population of low-mass compact binaries so as to have a wider range of S/Ns given the longer duration signals from these binaries. We chose the detection threshold such that the accumulated effective S/N ρ_eff > 8, and the S/N for the i^th segment in the j^th detector ${\rho }_j^i>3$ in at least one segment of 5 minutes in duration, for j = (1, 2, 3) corresponding to the three ET detectors comprising the single ET. Then, for each segment, and for each of the detected binaries, we estimated $\frac{{\mathrm{\Lambda }}_{\mathrm{s}}}{{\mathrm{\Lambda }}_{\mathrm{med}}}$ (defined in eq. (27) and eq. (28)) assuming that the measurement errors on the S/Ns are Gaussian, such that the standard deviations for ρ is σ_ρ = 1.

Figure A.1a shows $\frac{{\mathrm{\Lambda }}_{\mathrm{s}}}{{\mathrm{\Lambda }}_{\mathrm{med}}}$ as a function of the S/Ns in the three ET detectors, ρ₂₁ and ρ₃₁ in the left and right panels, respectively. It can be clearly seen that for the bulk of the sources, $1\lesssim \frac{{\mathrm{\Lambda }}_{\mathrm{s}}}{{\mathrm{\Lambda }}_{\mathrm{med}}}\lesssim 2$ for a very narrow range of 0.9 ≲ ρ₂₁ ≲ 1.1 and 0.9 ≲ ρ₃₁ ≲ 1.1. We therefore assume that $\frac{{\mathrm{\Lambda }}_{\mathrm{s}}}{{\mathrm{\Lambda }}_{\mathrm{med}}}=\mathcal{F}({\rho }_{21}^i,{\rho }_{31}^i)$, where $\mathcal{F}({\rho }_{21}^i,{\rho }_{31}^i)$ is as defined in eq. (29).

Fig. A.1. Representation of ℱ. (a) Λs/Λmed as a function of ρ21, ρ31. (b) 3D representation of the function ℱ, which we use as an approximation for Λs/Λmed. (c) Zoomed-in range of (b).

In Figure A.1b and A.1c, we show the function ℱ(ρ₂₁, ρ₃₁) for all the segments of all the binaries that cross the detection threshold. Figure A.1b shows the functional form, while Figure A.1c shows the zoomed-in range to show the clear range of variation of ℱ with ρ₂₁, ρ₃₁. It can be seen that including this function $\mathcal{F}({\rho }_{21}^i,{\rho }_{31}^i)$ in Equation 39 in SB2 as a prior information only affects a very narrow range where 0.95 ≲ ρ₂₁ ≲ 1.05 and 0.95 ≲ ρ₃₁ ≲ 1.05, which is the bulk of population as seen in Figure A.1a, and has negligible effect outside this range.

Appendix B: Effect of ℱ on mock populations

In this paper, we constructed multiple different mock populations as described in §3. For each of these mock populations, we show the variation of ρ₂₁ and ρ₃₁ over the whole range of redshift in Figure B.1. As mentioned in section §A, including ℱ as prior information only affects a very narrow range where 0.95 ≲ ρ₂₁ ≲ 1.05 and 0.95 ≲ ρ₃₁ ≲ 1.05. As can be seen from Figure B.1, there are ≪1% binaries at z ∼ 3 for Mock 1, in the range 4 ≲ z ∼ ≲6 for Mock 2, and at z ∼ 13 for Mock 3 and most of these generate S/Ns outside the range 0.95 ≲ ρ_ij ≲ 1.05 where i,j are two of the three ET detectors. Therefore, including ℱ as prior information leads to no improvement in the bias for these redshift ranges in the respective mock populations.

Fig. B.1. S/N as a function of redshift. Here, we show ρ21 and ρ31 as a function of redshift for compact binaries of (a) Mock 1, (b) Mock 2, and (c) Mock 3 where the subscript denotes one of the three ET detectors.

Appendix C: Reconstruction with extreme time-delay distributions

We concluded in §6, based on Figure 5, that our reconstruction of the SFR is essentially independent of the time-delay distributions up to z ∼ 14 and that the accuracy of the reconstruction of SFR strongly depends on the time-delay distribution only at higher redshifts beyond z ≳ 14. In order to further prove this point, we assume three more time-delay distributions: (i) t⁻³, (ii) t⁻⁴, and (iii) t⁻⁵, where t_del, min < t < t_H Gyr. The results for reconstructed SFR assuming these time delays and assuming that we know the merger rate as a function of redshift with high accuracy are shown in Figure C.1. We can see that for these time-delay distributions, the accuracy of the reconstruction of SFR also strongly depends on the time-delay distribution only at higher redshifts of z ≳ 14.

Fig. C.1. Reconstructed SFRs for the three populations assuming that the merger rate density is known with high accuracy. The maximum value of the SFR is normalised to 1. The left and right panels show the reconstruction assuming that tdel, min = 0.03 and 0.05 Gyr, respectively. The top three panels in all three plots show RE, or the KL divergence, for a redshift bin of Δz = 0.5 in width. The red curve denotes the true SFR. We note that for the Mock 3 population, we estimate the merger-rate-weighted SFR.

Appendix D: Acronyms used in the text

Table D.1 lists the various names of the variables used throughout the text in this paper.

All Tables

All Figures

	Fig. 1. Star formation rates used to construct the mock populations in this analysis. SFRs for Pop III are adopted from De Souza et al. (2011) and the SFR for Pop I+II is taken from Belczynski et al. (2020).
In the text

Fig. 2. Estimation of parameters. Density distribution of the median of the estimated posterior with respect to the true values of the parameters: Chirp mass ℳ (left) and redshift z (right) of the detected compact binary sources. Top panel: Mock 1. Middle panel: Mock 2. Bottom panel: Mock 3. The magenta line is a reference for equal values of true and estimated parameters. The blue contour encloses the 90% probability region.

In the text

	Fig. 3. Estimate of the merger rate density. Left: merger rate densities for Mock 1, Mock 2, and Mock 3, respectively. The shaded region represents Poisson error. Middle: Relative merger rate density. Right: cumulative probability distribution of the redshifts for Mock 1, Mock 2, and Mock 3. A description of the acronyms is given in Table D.1.
In the text

Fig. 4. Reconstruction of merger rate density. Left: cumulative probability distribution of the redshift of the secondary population (red) and of the detected sources from these secondary mock populations (blue) for Mock 1, Mock 2, and Mock 3. The number in parentheses is the number of sources in each population. Right: reconstructed merger rate density for Mock 1, Mock 2, and Mock 3. The shaded region represents the Poisson error.

In the text

Fig. 5. Reconstructed SFRs for the three populations, assuming that the merger rate density is highly accurate. The maximum value of the SFR has been normalised to 1. The left and right panels show the reconstruction assuming that tdel, min = 0.03 and 0.05 Gyr, respectively. The top three panels in all three plots show RE or the KL divergence for a redshift bin of Δz = 0.5 in width. The red curve denotes the true SFR. We note that for the Mock 3 population, we estimate the merger-rate-weighted SFR.

In the text

Fig. 6. Reconstructed SFRs for the three populations with the reconstructed merger rate density Rmer, recon. The maximum value of the SFR has been normalised to 1. The left and right panels show the reconstruction assuming that tdel, min = 0.03 and 0.05 Gyr, respectively. The top three panels in all three plots show RE or the KL divergence for a redshift bin of Δz = 0.5 in width. The red curve denotes the true SFR. We note that for the Mock 3 population, we estimate the merger-rate-weighted SFR.

In the text

	Fig. A.1. Representation of ℱ. (a) Λs/Λmed as a function of ρ21, ρ31. (b) 3D representation of the function ℱ, which we use as an approximation for Λs/Λmed. (c) Zoomed-in range of (b).
In the text

	Fig. B.1. S/N as a function of redshift. Here, we show ρ21 and ρ31 as a function of redshift for compact binaries of (a) Mock 1, (b) Mock 2, and (c) Mock 3 where the subscript denotes one of the three ET detectors.
In the text

Fig. C.1. Reconstructed SFRs for the three populations assuming that the merger rate density is known with high accuracy. The maximum value of the SFR is normalised to 1. The left and right panels show the reconstruction assuming that tdel, min = 0.03 and 0.05 Gyr, respectively. The top three panels in all three plots show RE, or the KL divergence, for a redshift bin of Δz = 0.5 in width. The red curve denotes the true SFR. We note that for the Mock 3 population, we estimate the merger-rate-weighted SFR.

In the text