Impact of Population III stars on the astrophysical gravitational-wave background

We probe the astrophysical gravitational-wave background resulting from compact binary coalescences, focusing on Population III binary black holes. We exploit results of state-of-the-art simulations on the evolution of Population I-II and III binaries, considering a variety of initial condition and star formation rate models for the latter. The contribution from Population III binary black holes is found to be very small, with no effect on the gravitational-wave spectrum. A network of third-generation detectors will detect easier individual Population III binaries, due to their significantly higher masses, hence decreasing even further their residual contribution.


I. INTRODUCTION
During its first three observing runs, the LIGO-Virgo-KAGRA (LVK) Collaboration has detected almost a hundred gravitational wave (GW) signals from compact binary coalescences (CBCs) [1][2][3], while at the end of the ongoing O4 run this number is expected to increase considerably [4].This multitude of detections allow us to continuously refine our understanding of the stellar population of compact binaries and constrain astrophysical models.Apart from the transient signals, there is a continuous effort to search for the gravitational-wave background, a random signal made up of the superposition of numerous GWs throughout the history of the Universe [5,6].
While a gravitational-wave background may result from various cosmological mechanisms [7], the astrophysical GW background (AGWB) made up from CBCs mergers is expected to be the dominant contributor [8,9].To estimate the AGWB, only binaries produced by Population I-II stars are usually considered.There is, however, an older population of stars -typically addressed in the literature as Population III stars -which could result in a sizeable contribution to the AGWB [10][11][12][13].
The increasing number of studies focusing on Population III stars in recent years [14][15][16][17][18][19][20][21][22][23][24][25][26][27][28], combined with the updated information gained on Population I & II stars from the LVK Collaboration detections, motivates us to revisit the AGWB from both Population I-II and Population III binaries.We compute the expected total AGWB from Population I-II & III binaries for a network of 3G detectors, as well as the unresolved AGWB made up of binaries that are not individually detectable.We show that the contribution of GW signals from Population III mergers to the total AGWB is negligible, as already found in previous studies [13].We then show that since compact binaries from Population III stars are expected to be much more massive than those of Population I-II stars, hence easier to be detected individually in the astrophysical GW foreground, their contribution to the AGWB is completely lost after foreground subtraction.
This paper is organised as follows: In Section II, we revisit the standard formalism to study the AGWB, focusing on an approach with catalogues of sources.We then describe the simulated population models and the corresponding compact binary parameters.In Section III, we outline the process to remove the astrophysical GW foreground detected by a 3G detector network.In Section IV we present our findings for the total and unresolved AGWB.We also discuss statistics for Population I-II and III binaries in terms of the signal-to-noise ratio.In Section V, we summarise our conclusions and briefly comment on the differences between our and related previous studies.

II. ASTROPHYSICAL GRAVITATIONAL-WAVE BACKGROUND
We briefly outline the formalism to estimate the AGWB.We then summarise our method to build catalogues of compact binaries and highlight our selected binary parameters.

A. Formalism
We describe the AGWB by the normalised GW energy spectrum [29] where ρ c = 3H 2 0 /(8πG) is the critical energy density, with H 0 the Hubble's constant and G the Newton's gravitational constant, c the speed of light, and f the GW frequency in the detector frame.The integrated flux density [30], depends on the merger rate, R(z), the comoving distance, d c (z), and the GW energy spectrum in the source frame, dE GW /df s , with f s = f (1 + z) the GW frequency in the source frame.We express the merger rate R(z) in terms of the merger rate per unit comoving volume, R V (z), as where We adopt a flat ΛCDM model with Ω m = 0.31, Ω Λ = 0.69 and H 0 = 67.66km/s/Mpc[31].The (1 + z) factor accounts for cosmic expansion.
The GW energy spectrum for a single source reads [32] where h+ , h× are the Fourier transforms of the two polarisation modes of the gravitational waveform in the frequency domain.
Given a catalogue of sources, the GW energy spectrum, Eq. ( 4), can be computed for each source individually.The integrated flux density, Eq. ( 2), can then be replaced by a sum over all sources [33]: where h+ , h× are now taken in the detector frame, where i denotes a particular source.The total number of sources, N , is related to the total time of observation, T , as: In what follows, we use catalogues of simulated GW sources, which have been stored in redshift bins.We thus perform a discrete summation to calculate N , instead of computing the integral expression given in Eq. (6).The number of sources at each redshift bin is then given by with total number of sources We opt for a redshift resolution ∆z i = 0.1, ∀i and fix N = 10 5 for statistical convergence (as in Ref. [34]).
The SEVN binary population synthesis code, provides 3 of the binary parameters, namely the two compact binary masses and the eccentricity of the orbit.
The high number of considered IC and SFR models aims to account for the high uncertainty that characterises the properties of Population III stars.For a brief description of the different models, see Appendix A.
Once we have the Star Formation Rate (SFR), we use the (publicly available) code CosmoRate [43,44] to getfor a given redshift bin -the Merger Rate Density (MRD) and a set of binaries, from which we randomly select N i , Eq. ( 7), to produce a catalogue.
The aforementioned steps to produce a catalogue are summarised in the flowchart, Fig. 1.
As mentioned previously, SEVN provides only the compact binary masses in the source frame, m 1(source) and m 2(source) (where 1 stands for the primary and 2 for the secondary compact object), and the eccentricity, e.We have to choose 8 additional parameters: 2 intrinsic parameters (spins) and 6 extrinsic ones (polarisation angle, coalescence phase, right ascension, declination, inclination angle, and coalescence time).We sample as follows (see, Table I): • The spins in the z-direction, χ 1(z) and χ 2(z) , are sampled from uniform distributions in the range [−0.75, 0.75] for black holes, and [−0.05, 0.05] for neutron stars [47].We set the spins in the x-and y-direction to 0.
• The coalescence time in the detector data segment, t c , is set at 0 for all signals.
As for the gravitational waveform we use to model the signals, we have opted for IMRPhenomXHM [48].This choice has been motivated by a preference for probing the harmonic modes of the signals to compute the AGWB (as considered in Ref. [12]), rather than precession.

III. UNRESOLVED AGWB
An AGWB present in the data will be dominated by a foreground of loud detected signals; removing them from the dataset would allow to estimate the AGWB truly composed of individually unresolvable signals.
The unresolved contribution to the Ω GW is made up of all signals below some SNR threshold.For the i−th binary, we define the network optimal SNR as where is the detector response, P I the power spectral density of detector I, and M the total number of considered detectors.We denote by F lp the location-phase factor, and by F + , F × the antenna pattern functions.We consider a three 3G detector network: an L-shaped Cosmic Explorer (CE) detector with 20km-long arms and post-merger optimised PSD located at the site of LIGO Hanford, an L-shaped CE detector with 40km-long arms and compact-binary optimised PSD located at the site of LIGO Livingston [49], and a triangular (that is, composed of three distinct interferometers) Einstein Telescope (ET) detector with 10km-long arms located at the site of Virgo [50].We set the minimum frequency at f min = 5 Hz and the maximum frequency at f max = 1024 Hz [51], as contributions to the SNR outside this frequency range are negligible.To compute the SNR we use the GWBENCH package [52].
We compute the optimal SNR of every signal in our catalogue.The calculation of the unresolved AGWB follows the same procedure as outlined for the total AGWB, with the sole difference that this time we only consider the signals in our catalogues that pass the SNR threshold.We consider two choices for the SNR threshold: ρ thr = 12 (conservative), and ρ thr = 8.

IV. RESULTS
In what follows, we present the total AGWB produced separately by mergers of Population I-II stars and Population III stars, and discuss the characteristics of the latter versus the former.We then discuss the difference in the unresolved AGWB for Population I-II versus Population III, noting that this depends on the detector sensitivities.In our study, we have used a three 3G detector network and fixed the SNR threshold.

A. Total background
We plot in Fig. 2 the total Ω GW from mergers of all three binary classes (BBHs, BHNSs, and BNSs) produced by Population I-II stars, as well as the individual contributions.We also plot the Ω GW resulting from Population III binaries, considering the most optimistic (J19, LAR1 ) and most pessimistic (SW20, TOP5 ) cases (we refer the reader to Appendix B, where we present results for all Population III cases).We also show the total AGWB from Population I-II, and III binaries, where for the latter we selected the most optimistic model.
We observe that for Population I-II, the BHNS contribution is the weakest.The AGWB seems to be dominated by the BBH mergers in the whole frequency range except for the highest frequencies, above f ∼ 900 Hz, where the BNS contribution prevails.The contribution from Population III binaries peaks at ∼ 10 −13 or ∼ 6 × 10 −11 depending on the assumed SFR and IC model.In all cases, we observe a plateau which in the most optimistic scenario is in the region f ∈ [20,200] Hz.Clearly, even in this most optimistic case, the resulting increase in the total AGWB after including BBHs from Population III stars is only minimal, hence insufficient to lead to significant deviation from the expected ∝ f 2/3 spectrum [32].Considering the contribution from Population I-II and III binaries, the Ω GW peaks just above 10 −9 .
We next investigate the AGWB from Population I-II and Population III binaries when considering only unresolved sources.

B. Unresolved background
We plot in Fig. 3 the unresolved AGWB for our two choices of ρ thr .First, we consider the Population I-II contribution.Imposing an SNR threshold of 12 seems to drop the Ω GW by ∼ 1 order of magnitude in the whole frequency range except for the highest frequencies, where the BNS mergers contribute most.This is expected, as GW signals from BNSs are weaker than those from BBHs and dominate after foreground subtraction.
We next consider Population III for the most optimistic case (J19, LAR1 ).One would expect BBH mergers from Population III stars to be overall harder to detect than those of Population I-II stars, given that the merger rate density for Population III binaries peaks at higher redshifts (see Appendix B for a detailed figure of the computed MRD).We observe, however, that the corresponding unresolved Ω GW , for an SNR threshold of 12, lies ∼ 2 orders of magnitude below the total Ω GW for f ≤ 30 Hz, where it follows the ∝ f 2/3 dependence.Notice that Ω GW has lost it characteristic plateau, dropping abruptly for frequencies above f ∼ 30 Hz.Such frequencies correspond to BBHs with very low redshifts, thus high SNRs, therefore removed by the foreground subtraction.
To further quantify the difference in the foreground Note that these results assume a perfect foreground subtraction.Even though this is not the case [34,53], we are not considering any errors that would decrease the efficiency of the subtraction, since the contribution from Population III binaries to the total AGWB is quite small, and lost after foreground subtraction.

V. DISCUSSION & CONCLUSIONS
We have studied the contribution from Population III compact binaries mergers to the astrophysical gravitational-wave background and compared to that from Population I-II.For Population I-II binaries, we have considered all three types of mergers (BBH, BNS, BHNS).For Population III, we have considered a variety of models with different initial conditions for the binary population synthesis code and different star formation rate.
Our analysis has determined that including GW signals from Population III binaries has a very small impact on the total Ω GW , as already shown in Ref. [12].Our study, however, has yielded results that contradict expectations on the unresolved AGWB.Specifically, we find Population III binaries to be on average easier to detect than Population I-II binaries, regardless of the assumed initial conditions and star formation rate of the former.For a given SNR threshold, foreground subtraction has a greater impact in the case of Population III binaries, resulting in their contribution to the total AGWB being lost.Hence the AGWB is characterised and distinguished from the cosmological GWB by its ∝ f 2/3 spectrum.Loud individual GW signals from Population III binaries could be present in the data and identifiable from their higher redshifts compared to Population I-II FIG. 4. Histogram with the network optimal SNRs of all BBHs for Population I-II & Population III (J19, LAR1 ).Despite the higher redshifts typically associated with Population III BBHs compared to Population I-II, the distribution of the former is shifted towards higher SNRs due their significantly higher masses.

binaries.
Let us note that the difference in our results with respect to previous ones [12,13] is due mainly to the different merger rates and binary masses of the considered catalogues.To build our catalogues, we have employed results from the SEVN binary population synthesis code, whereas the study of Ref. [12,13] is based on simulations with StarTrack [54].We refer the reader to models M10 and FS1 for Population I-II and III, respectively, of Ref. [18].In our models the MRD peaks at z ≃ 4.5 and z ≥ 7.5 for Population I-II and Population III, respectively, whereas for the StarTrack models the MRD peaks correspondingly at z ≃ 2 and z ≃ 12.For our models, the average BBH total source mass is M source = 16.8M⊙ for Population I-II, and ranges from M source = 49.6M⊙ (H22, KRO1 ) to M source = 64.0M⊙ for Population III binaries (see Appendix B for a detailed figure of the total source mass distribution), whereas for the StarTrack models the corresponding masses (for mergers within z < 2) are M source = 29.7M⊙ and M source = 63.4M⊙ .Finally, we note that the BHNS and BNS waveform adopted in Ref. [12,13] considered only the inspiral phase.LIGO Laboratory, which is a major facility fully funded by the National Science Foundation.
The authors acknowledge computational resources provided by the LIGO Laboratory and supported by NSF Grants PHY-0757058 and PHY-0823459.
We are indebted to Filippo Santoliquido for many useful discussions on Population III stars.We also thank Ssohrab Borhanian for suggestions on how to run GWBENCH.NK thanks Ansh Gupta, Claire Rigouzzo, Michelle Gurevich, and MS thanks Tania Regimbau for discussions.We thank Nelson Christensen for carefully reviewing this work as a part of the LVK collaboration's internal review process.
NK is supported by King's College London through an NMES Funded Studentship.MS acknowledges support from the Science and Technology Facility Council (STFC), UK, under the research grant ST/X000753/1.This manuscript was assigned LIGO-Document number P2400127.

Mass ratio & secondary mass
The secondary mass distribution was either identical to the primary mass distribution, or determined assuming a power law distribution [63,76] for the mass ratio.In both cases, the final mass ratio distribution also depends on the primary mass distribution.

Orbital period
The orbital period distribution is either a power law favouring shorter periods [76], or a Gaussian distribution [63].

Eccentricity
The eccentricity distribution is a power law with either positive [15,16,20] or negative [76] exponent, favouring either higher or lower eccentricities, respectively.

Star formation rate
All considered models for the star formation rate are consistent with the Thomson scattering optical depth value estimated by the Planck Collaboration [77].Their peak varies significantly, from z ≃ 8 (J19 ) to z ≃ 20 (SW20 ), depending on different physical assumptions.First, H22 is a semi-analytic model that samples and traces individual stars, based on dark matter halo merger trees and calibrated to reproduce observables in the Universe [78,79].Next, J19 was obtained using the hydrodynamical/N -body code GIZMO [80], considering both the chemical and radiative feedback of core-collapse and pairinstability supernovae.Likewise, LB20 was the result of simulations with GIZMO, extrapolated to lower redshifts and following a Madau-Dickinson form [81]. Finally, SW20 was constructed from hydro-dynamical cosmological simulations that ran on the adaptive mesh refinement code ENZO [82].

Total background
We plot in Fig. 5 the Ω GW from Population III BBH mergers, for the 4 SFR models discussed in Section II B. Each of the 4 panels considers a single SFR model (listed at the top left corner), and 11 different IC models.Note that the individual curves appear in general wobbly in the region f ≥ 200 Hz.This is because high frequencies in the Ω GW correspond to signals that are emitted by low redshift sources, which are inevitably few in our catalogues.Thus, Ω GW at these frequencies is characterised by higher statistical uncertainties.Notice that in all cases there is a plateau, as discussed in Section IV A.

Merger rate density
We show in Fig. 6 the MRD over redshift for Population I-II and Population III BBHs.Even in the most optimistic case (J19, LAR1 ), the Population III MRD peak value (70.6 Gpc −3 yr −1 ) is significantly lower than the one for Population I-II BBH MRD peak value (172.2Gpc −3 yr −1 ), and lies at a much higher redshift (4.55 and 7.65, respectively).In all other cases, the MRD peaks at z ≥ 7.65.

Mass distribution
We plot in Fig. 7 the total source mass distribution for Population I-II and Population III BBHs.This is ex-pected to be generally higher for the latter because of the corresponding extremely low metallicity, which implies that (a) the initial mass function is more top-heavy as compared to metal-rich stars and (b) there is virtually no mass lost in stellar winds [83,84].Indeed, the M source distribution takes typically higher values for all Population III models compared to Population I-II.Interestingly, we notice that (J19, LAR1 ), which appears rather pessimistic in terms of the expected M source , turns out to be the most optimistic case in terms of the Ω GW (as seen in Fig. 5) as a result of its high MRD.Note that the fast drop around M source ∼ 100 can be associated to the pair-instability mass gap, which for Population III BBHs simulated with SEVN has a lower edge at M source = 86M ⊙ [27].

FIG. 2 .
FIG. 2. The AGWB resulting from mergers of BBHs, BHNSs, and BNSs produced by Population I-II stars, as well as their total.Including the BBH mergers from Population III stars has a negligible impact even when considering the most optimistic model.The AGWB for the most pessimistic Population III model is also provided for comparison.

FIG. 3 .
FIG.3.The unresolved AGWB resulting from mergers of binaries produced by Population I-II & III stars, for an SNR threshold ρ thr ∈ {8, 12}.The AGWB before the foreground subtraction is also plotted for comparison.Removing the foreground of unresolvable sources from our catalogue seems to reduce the AGWB more in the case of Population III binaries.Finally, we provide the total AGWB from all three Populations before and after the foreground subtraction.

FIG. 5 .
FIG. 5.The GWB resulting from mergers of BBHs produced by Population III stars, for 4 different SFR models (4 subfigures) and 11 different initial condition models.All models have been taken from [26, 27].

FIG. 6 .
FIG. 6.The BBH merger rate density (MRD) for Population I-II and the J19, LAR1 Population III model, as computed using CosmoRate.The MRD for every other Population III model is also plotted.

FIG. 7 .
FIG. 7. The BBH total source mass distribution for Population I-II and the J19, LAR1 Population III model.The latter is clearly characterised by more massive, on average, BBHs.The Msource distribution for every other Population III model is also plotted.

TABLE I .
Summary of the distributions we use to sample in all compact binary parameters.