Reconstructing the star formation rate for compact binary populations with the Einstein telescope

The Einstein Telescope (ET) is a proposed third-generation, wide-band gravitational wave (GW) detector. Given its improved detection sensitivity in comparison to the second-generation detectors, it will be capable of exploring the Universe with GWs up to very high redshifts. In this paper, we present a population-independent method to infer the functional form of star formation rate density (SFR) for different populations of compact binaries originating in stars from Population (Pop) I+II and Pop III using ET as a single instrument. We use an algorithm to answer three major questions regarding the SFR of different populations of compact binaries. Specifically, these questions refer to the termination redshift of the formation of Pop III stars, the redshift at peak SFR, and the functional form of SFR at high redshift, all of which remain to be elucidated. We show that the reconstruction of SFR as a function of redshift for the different populations of compact binaries is independent of the time-delay distributions up to $z \sim 14,$ and that the accuracy of the reconstruction only strongly depends on this distribution at higher redshifts of $z\gtrsim 14$. We define the termination redshift for Pop III stars as the redshift where the SFR drops to 1\% of its peak value. In this analysis, we constrain 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 for Pop III stars, which have a true separation of at least $\Delta z \sim 2$. The accurate estimation of the termination redshift depends on correctly modelling the tail of the time-delay distribution, which constitutes delay times of $\gtrsim 8$ Gyr.


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 1Hz (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 2 −10 4 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 5 − 10 6 BBH detections and ∼ 7 × 10 4 binary neutron star (BNS) detections in one year (Regimbau et al. 2012(Regimbau et al. , 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 (The LIGO Scientific Collaboration et al. 2021).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 phys-⋆ e-mail: singh@lapth.cnrs.frical 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 − 100M ⊙ , 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 Rochelobe 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 PopI and II stars (Belczynski et al. 2017).This leads to a smaller number of binary interactions for PopIII stars (RLOF and CE) as compared with PopI 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.

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 §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 §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 §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 §6.We discuss our conclusions in §7.

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 6 cm −3 , 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 We refer to these SFRs as SFR1, SFR2, and SFR3, respectively. These SFRs are shown in Figure 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.
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(Belczynski et al. , 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 Figure 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.

Construction of the mock populations
For both M30B and FS1 compact binaries, we use the population available on the StarTrack1 website.For the FS1 model, we only consider the following data for each binary: (i) the masses of the merging compact objects m i 1,2 , and (ii) the delay time between the formation of the binary at zero-age main sequence (ZAMS) and its coalescence, t i del .For compact binaries of the M30B model, we use the following data for each binary: (i) the masses of the merging compact objects m i 1,2 , (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: 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 8 M ⊙ Belczynski et al. (2020) and 3.5 × 10 9 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 R b (z): 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 Equation ( 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.

Estimation of the parameters of compact binaries with ET
For a merging compact binary of chirp mass M, 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 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 1 and m 2 , the chirp mass is M = (m 1 m 2 ) 3/5 /M 1/5 , and the total mass is M = m 1 + m 2 .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 where t 0 is the time of coalescence in the detector frame, and so (t 0 − 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 by substituting the values of the two polarisations from Equation (3) into Equation ( 4): The functions Θ and the phase Φ c in Equation (5) are functions of the antenna response functions F + and F × and the angle of inclination ι, and are defined as 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) where the redshifted chirp mass M z = (1 + z)M and the reference mass M BNS ≈ 1.218M ⊙ is the chirp mass of an equal mass binary with components of 1.4M ⊙ each.The function ζ ( f max ) is defined as 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, where r 0 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 where the effective antenna response function Θ eff is

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 obs gw is calculated using Equation (4.195) in Maggiore (2007): 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 i th segment, the corresponding values f i−1 , f i of f obs gw will be and The S/N for the i th segment in the j th detector can be written using equation ( 8), as where and where F i + and F i × are the antenna response functions for the j th detector in the i th segment.The effective S/N for the i th segment is where and the function Θ i e f f is We assume that the observables, in the case of the detection of a coalescing binary system, are: (a) the three S/Ns ρ i j for the i th segment of the signal, (b) the phase Φ i o, j for j = (1, 2, 3) for each of the three ET detectors in the i th 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 M 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 ρ i j and Φ i j 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;The LIGO Scientific Collaboration et al. 2021).
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 Θ i eff (see Sec IV of SB2): where ρ i 21 ≡ ρ i 2 /ρ i 1 and ρ i 31 ≡ ρ i 3 /ρ i 1 in the i th segment.We then proceed to constrain the source-dependent quantity Λ (defined in Equation 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: 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 Figure 13 of SB2) given the antenna pattern functions of the triangular configuration of ET (see Appendix of SB2 for detail).

Modification in the algorithm used in SB2
By comparing both sides of Equation ( 25), we find a function F (ρ i 21 , ρ i 31 ), such that where the subscript 's' denotes the true value of Λ from the actual source parameters, and Λ med is the median of the estimated posterior distribution of Λ in the i th segment of the signal: We approximate F (ρ i 21 , ρ i 31 ) to a fit of the form where R is the Marr or Mexican hat function, and is given as A more detailed description of this fit is described in §A.We include this function F (ρ i 21 , ρ i 31 ) in Equation 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 i th segment in the j th detector ρ i j > 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 §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 M 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 §4.2.The chirp mass, total mass, and redshift of these detected sources are denoted M 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 M med,det , M med,det , and z med,det , respectively.

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 Figure 2.This figure shows the estimated median values (M med,det , z med,det ) with respect to the actual values of the parameters (M 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 Figure 2a and 2b, can be compared with Figures 3(c,d) of S22 in order to see the effect of using F (ρ i 21 , ρ i 31 ) 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 Figures 2a, 2c, and 2e.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 Figures 2b, 2d, and 2f, 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 F (ρ i 21 , ρ i 31 ) as an additional prior only improves the parameter estimates for those binary sources for which 0.95 ≲ ρ i j ≲ 1.05, where ρ i j is the ratio of S/N generated in the i th and j th detector of ET.This is explained in more detail in §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.

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.

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 where N yr is the number of mergers per year calculated by integrating the merger rate density given by Equation ( 1).The merger rate density R mer for a given population is calculated for the time T mock : where N (z i ,z i+1 ) is the number of mergers in a redshift bin [z i : From Equation (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 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 R mer (z med,det ) R mer (z s,mock ) and R mer (z med,det ) R mer (z s,det ) .These are shown in Figures 3b, 3e, and 3h.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, 3f, and 3i 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 §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 ).

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 populationdenoted with subscript 'sec'-from the detected sources, assuming that the distributions of the chirp mass and redshift are proportional to M med,det and z med,det , that is, p(M sec ) ∝ p(M 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: 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 D as a function of redshift is defined as where [N sec ] (z i ,z i+1 ) is the number of mergers in the secondary mock population in the redshift bin (z i , z i+1 ) and [N sec,det ] (z i ,z i+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 Figure 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.

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 R mer,recon (z i , z i+1 ) = R mer (z med,det ) D The reconstructed merger rate density is calculated for the three mock populations using Equation ( 34) in ( 35).The reconstructed merger rate densities for Mock 1, Mock 2, and Mock 3 are shown in the right panel Figure 4.The red lines are the true merger rate density, while the green lines show the merger rate density calculated using Equation ( 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 M 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 Figure 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.

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 −1 , and (iii) t −2 , where t del,min < t < t H Gyr. The broken power-law distribution is assumed to be: The broken power law time-delay distribution is motivated by the fit to the actual time-delay distribution as shown in Figure 11 of Belczynski et al. (2017) for FS1 compact binaries.The time-delay distribution t −1 is based on the best fit for the Pop I+II compact binaries (M30B time-delay data is available on the StarTrack 2 website).The SFR estimate at a given z is the sum over contributions from each binary 'b', averaged over its delay time: for z b ini ≡ z(t b mer − t b del ).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-rateweighted 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: 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.
For each of the three mock populations, the reconstruction of SFR as a function of redshift is independent of the assumed timedelay 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 timedelay 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 Figure 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 Figure 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 Figure 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 Sect C of the Appendix.
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.
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 timedelay distributions and the error on the estimates of the merger rate density as a function of redshift.
The top three panels in Figure 6a, Figure 6b, and Figure 6c 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 Figure 6c.As seen in Figure 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.

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 the 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 timedelay 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.

Fig. 1 :
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).
) and are thought to be less massive (∼ 40 − 60M ⊙ ) than Pop III.1 stars (∼ 1000M ⊙ ).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, 100km/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.

Fig. 2 :
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 M (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.
Figures 3a, 3d, and 3g 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 :
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.

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

Fig. 5 :
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 t del,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.

Fig. 6 :
Fig.6: Reconstructed SFRs for the three populations with the reconstructed merger rate density R mer,recon .The maximum value of the SFR has been normalised to 1.The left and right panels show the reconstruction assuming that t del,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.

Table D .
1: List of variables used frequently in the text Notation Definition Reference t ini cosmic time of the formation of a ZAMS binary t ini = t obs − t del t obs cosmic time at which the compact binary is observed to merge.t del delay time t del = t evol + t merg t evol time of evolution from ZAMS to the formation of a compact binary system sim total mass of all stars accompanying the stellar evolution, leading to formation of compact object binaries including the binaries as well as the single stars.M s,mock , M s,mock , z s,mock subscript (s, mock) denotes true source parameters of a binary in the mock population M s,det , M s,det , z s,det subscript (s, det) denotes true source parameters of detected sources M med,det , M med,det , z med,det subscript (med, det) denotes the median of the estimated see Fig. 5 in SB2 ′ merger-rate-weighted SFR