Footprints of population III stars in the gravitational-wave background

We investigate detection prospects of the gravitational-wave background (GWB) that originates from the merging of compact objects formed by the collapse of population III stars. Younger population I/II stars lead to a GWB in the LIGO/Virgo frequency band at the inspiral phase, while population III stars would likely show up at the later merger and ringdown phases. We show that, using a network of third-generation detectors, we may be able to separate a population I/II signal from a population III one, provided we can subtract individual coalescence events. A detection of a population III GWB could reveal important information, such as the average redshifted total mass.

We investigate detection prospects of the gravitational-wave background (GWB) that originates from the merging of compact objects formed by the collapse of population III stars. Younger population I/II stars lead to a GWB in the LIGO/Virgo frequency band at the inspiral phase, while population III stars would likely show up at the later merger and ringdown phases. We show that, using a network of third-generation detectors, we may be able to separate a population I/II signal from a population III one, provided we can subtract individual coalescence events. A detection of a population III GWB could reveal important information, such as the average redshifted total mass. Introduction-We have witnessed a rapid expansion of gravitational-wave (GW) astrophysics in the last decade due to the success of the Advanced LIGO and Virgo GW detectors [1,2] at uncovering signals from numerous compact binary coalescences (CBCs) [3]. The most recent Advanced LIGO/Virgo observing run, O3, presented us with dozens of new merger events and has significantly expanded the stellar graveyard [4]. Despite the increase in the number of detections, we are yet to observe with confidence an event that would suggest that the progenitor compact objects are remnants of the oldest stars in the Universe [5] -the theoretically-postulated population III (pop III) stars (GW190521 could be a potential candidate [6]). Pop III stars are thought to have formed at high redshifts and as such have low metallicity compared to the more recently formed, population I/II (pop I/II) stars [7][8][9]. These old stars have hitherto evaded sky surveys [10][11][12], and their detection remains an objective for upcoming experiments, such as the James Webb Space Telescope [13].
Pop III stars may solve some of the puzzles in black hole formation, as well as help understanding the early epochs of the Universe such as reionisation and galaxy evolution [14][15][16]. Numerical simulations show that these primordial stars could have led to the formation of supermassive black holes at high redshifts [17][18][19]. Mergers of such heavy remnants would appear in the millihertz frequency range explored by future space-based detectors such as LISA. The scope of this study, however, is detection prospects of terrestrial detector networks, and we therefore focus on models that predict a pop III signal in the LIGO/Virgo frequency range. The contribution to the gravitational-wave background (GWB) 1 from a superposition of unresolved pop III-seeded CBCs has been explored in several studies [21][22][23][24][25]. They show significant deviation of a pop III star signals from a pop I/II stars 1 Often referred to as the stochastic gravitational-wave background [20].
signal due to different mass and redshift distributions [25].
The GWB is comprised of many sources, of astrophysical or cosmological origin, but we expect the CBC signal to be the foreground to all sources [26]. In this study, for the first time, we consider the possibility of separating pop I/II and pop III GWB contributions. Numerous models suggest the total CBC background is dominated by pop I/II. However, pop III can be uncovered using subtraction techniques and studying the residual backgrounds [27,28]. As the sensitivity of detectors increases and GW interferometers see more individual CBC events, a pop III residual background emerges as the dominant signal over pop I/II residual background. We first study how to detect the GWB from pop III stars, and in the case of a successful detection, we explore subsequent implications, namely information about masses and redshifts of the population. The paper is organised as follows: in Sec. II we present the pop III models and their resulting GWB in different detector networks, highlighting characterisation of the residual backgrounds. We then introduce Bayesian analysis used in Sec. III, and discuss search filters we consider for pop III stars. Sec. IV is an implications study in case of a detection of a pop III signal. We use StarTrack (ST) simulation data [29] and apply our detection methods, ultimately showing consistency of our implications analysis with the underlying population. We select the ST data since it is the most recent extensive catalogue of merging binaries from pop I, II and III stars that lead to a GWB in the LIGO/Virgo frequency range [30].
Population III GWB-The GWB is defined as the superposition of GWs from all unresolved sources. It is characterised by the dimensionless parameter Ω GW [31], expressed as the ratio of GW energy density per logarithmic frequency bin dρ GW /d ln(f ), normalised by the critical energy density of the Universe, ρ c = (3H 2 0 c 2 )/(8πG) : with H 0 = 67.9 km s −1 Mpc −1 [32].
Here we concentrate on the CBC contribution to the GWB, namely from pop I/II stars and the theoretical pop III stars. One can express the quantity Ω GW in terms of CBC source parameters θ (masses and spins), as [30] : where p(θ) is the probability distribution of the source parameters, dE GW /df s in the energy density emitted by a single source at a redshift z with parameters θ, f s is the emitted frequency in the source frame f s = f (1 + z) and z up (θ) is the maximal redshift at which a compact binary with parameters θ can form. The factor (1 + z) in the denominator converts the merger rate R(z, θ) from the source to the detector frame, and E z (z) accounts for the considered cosmology, i.e. the expansion history of the Universe, E z (z) = Ω m (1 + z) 3 + Ω Λ , with Ω m = 0.31, Ω Λ = 0.69 [32].
For the total population, when all sources are included in the background, the merger rate R(z; θ) at redshift z for sources with parameters θ is given in the source frame, per unit of comoving volume and time. It is derived from the star formation rate, corrected by the time delay between the birth of the progenitors and the merger of the compact objects [33]. To calculate the residual background when individually detected sources are removed, one has to multiply the total rate by a factor 1 − (z, θ), where the efficiency (z, θ) is the probability for a source at redshift z with parameters θ to be detected, integrated over inclination, polarisation and position in the sky (see [34]).
In the case of binary neutron stars and neutron starblack hole mergers, we only consider the inspiral phase and assume that the emission of GWs stops at the last stable orbit. For binary black holes, we consider the three different regimes of the coalescence (inspiral, merger and ringdown phase) given by the corresponding phenomological waveforms [35] calculated for circular orbits. The energy density is [36]: where f merg , f ring and f cut are the frequencies at the start of merger, start of ringdown and end of emission in the source frame, respectively. The chirp mass of the binary is a combination of the individual masses of the compact objects, M c = (m1m2) 3/5 (m1+m2) 1/5 . L(f s , f ring , σ) is the Lorentzian function centered at f ring , with width σ, and w m , w r are the normalisation constants ensuring the continuity between the three phases. The factors i and α i and the frequencies f merg , f ring and f cut follow from analytical waveforms detailed in [35] and depend on the symmetric mass ratio η = (m 1 m 2 )/(m 1 + m 2 ) 2 of the progenitors' masses, m 1 and m 2 , and the effective spin of the system χ = [(m 1 s 1 + m 2 s 2 )/(m 1 + m 2 )] L/L. Once the spectrum of Ω GW is calculated, we can estimate the corresponding signal-to-noise ratio (SNR) for a given network of N detectors [38]: with T the observational time, P I and P J the one-sided power spectral noise densities of detectors I and J, and γ IJ the normalised isotropic overlap reduction function characterising the distance and the relative orientation between I and J for sources isotropically distributed in the sky.
In this work, we consider the StarTrack model FS1 for pop III [24]. Assumptions about characteristics of the initial binary pop III stars are discussed in [37]. The FS1 model assumes that pop III stars were formed in large gas clouds with a star formation rate that peaks at redshift z ∼ 12, while the star formation rate for pop I/II stars peaks at z ∼ 2 (see Fig. 4 in [24]). Even though pop III stars are less abundant than pop I/II, the ST model FS1 considers a rather optimistic ratio of pop III to pop I/II stars. The corresponding background and its detectability have been calculated [30] using a catalogue of sources rather than the analytical expression in Eq. 2. The residual background catalogue is obtained by subtracting all sources individually detected by the interferometer network. For each source k we calculate the individual SNR ρ k assuming optimal-matched filtering and uncorrelated gaussian noise in the detectors as follows: wherẽ with F +,I and F ×,I the antenna factors of detector I for polarisations + and × that depend on source inclination Θ k and position in the sky ψ k , whileh k + andh k × are the Fourier transforms of the gravitational waveforms of the source k. A residual catalogue is computed by removing all sources with ρ k > 12.
If pop III exists, its signal will be superposed with a pop I/II signal. In the case of a dominant pop I/II signal, the pop III signal will remain hidden underneath it, and one can only place upper limits on the amplitude of the pop III contribution to Ω GW (f ). If, however, a pop III signal is the dominant one, then we could detect deviations from the 2/3 CBC power law, and even getting insight on the mass and redshift distribution of pop III stars. We explore the last scenario by considering two terrestrial networks of third-generation (3G) detectors: (i) Einstein Telescope (ET) at the Virgo site, and (ii) ET at the Virgo site with two Cosmic Explorers (CE) at the LIGO Hanford and Livingston sites.
Estimates of CBC contributions to Ω GW from ST simulations suggest that pop III signal is lost in the pop I/II foreground. For 2G detector networks -even by including LIGO-Hanford, LIGO-Livingston, Virgo, LIGO-India, and KAGRA -pop III is practically invisible and its contribution to the global SNR is negligible, as it is shown in [30]. However, 3G detectors such as ET and CE, may reveal a pop III background. The future detectors will have unprecedented sensitivity and they will be able to discover a great number of individual CBCs, thereby reducing the GWB originating from unresolved CBCs. For ET+2CE, we uncover pop III after the subtraction of individually resolved merger events. This follows because subtraction methods are less efficient to detect the high redshift and low frequency pop III CBCs. Being more difficult to resolve, binaries from pop III persist, resulting in a large contribution to the residual CBC background in 3G detectors.
We compare in Fig. 1 the total and residual background for the two 3G networks: ET (top) and ET+2CE (bottom). It confirms that the pop III contribution in ET has a very small impact on the combined residual background from pop I/II and pop III, while in ET+2CE the pop III residual background clearly dominates for frequencies below ∼ 20Hz. In addition, Fig. 1 shows a change in the shape of the background: The peak frequency of pop III changes slightly while the slope characterising the end of emission decreases dramatically when we remove individually detected sources.
To demonstrate the impact of subtraction of resolved CBCs on the population, we show in Fig. 2 the probability density of the total redshifted mass, M z tot = (1 + z)(m 1 + m 2 ), and the merger rate R(z) as a function of redshift, between the whole catalogue and the residual one for ET+2CE. Clearly, the sources remaining in the residual catalogue are the ones with the highest redshift, affecting the total corrected mass distribution which is in turn responsible for the changes in the GWB spectrum. We will estimate the ET+2CE residual pop III parameters by filtering the corresponding background and performing a Bayesian analysis.
Detection method-The stochastic pipeline takes strain datas I,J from detectors, I, J, and constructs crosscorrelation statistics using optimal filters [38]: with T the duration of the run, and γ IJ (f ) the normalised overlap reduction function as defined in Eq. 4. The estimator is normalised with S 0 (f ) = (3H 2 0 )/(2π 2 f 3 ) leading to Ĉ IJ (f ) = Ω GW (f ). We assume correlated noise not to be a limiting factor to our detector sensitivity and  consider all noise to be gaussian. The variance is Let us construct a gaussian log-likelihood, where Ω GW (f |θ) represents the GWB model with parameters θ. This allows us to estimate the model parameters by finding the best-fit to the cross-correlation data and minimising the likelihood function. Note that we have made the simplifying assumption that the log-likelihood of a detector network is the sum of log-likelihoods of the individual baselines. To compare models and find which ones fit data better, we perform model selection with Bayes factors. Bayes factor, B M1 M2 , is defined as the ratio of evidences of model M 1 to model M 2 , and if large and positive, demonstrates preference for M 1 over M 2 .
Typically, we model a CBC signal as Hz. This is because the CBCs detected so far have low masses that would lead to an inspiral signal in the low-frequency range. This can be seen in Fig. 1 where the total GWB from pop I/II and III in 3G detectors are presented. The pop III spectrum shows clear deviation from a 2/3 power law, because these further away stars will lead to more redshifted frequencies and therefore be detected in their merger and ringdown phases. We test search filters different from a 2/3 power law to investigate if the deviation from pop I/II signal can be identified in a parameter estimation study. Motivated by the shape of the residual pop III signal in Fig. 1, we consider the following filters: • power law with varying spectral index (PL) • broken power law (BPL) • smooth BPL • triple BPL for f (1) peak < f ≤ f (2) peak , k Ω peak (f /f (2) peak ) α3 for f > f (2) peak , (13) where k = (f (2) peak /f (1) peak ) α2 ensures continuity of the piecewise function.
The priors for each model's parameters can be found in the Appendix. If any of the filters above are preferred over a 2/3 filter, this could be an indication of the presence of a pop III signal.
Implications-In the case of a detection, we examine whether we can constrain the mass/redshift distribution from the optimal search parameters. Following the GWB expression (Eq. (2)), we see that the parameters impacting the background shape are the redshift-dependent merger rate and the black holes' mass distribution. To understand how these population characteristics relate to model parameters, such as peak frequency and slope, we generate multiple spectra. We make simplifying assumptions about our progenitors by assuming spinless, equal-mass binaries [39]. We fix the merger rate and vary the intrinsic mass input, observing how the shape of the GWB spectrum changes. The results we find, however, change with a different choice of merger rate, as described in the Appendix. This is because there is a degeneracy between the effects that merger rate and mass distribution have on the GWB [23]. We thus study the dependence of Ω GW on the redshifted total mass of the population, M z tot = (1 + z)(m 1 + m 2 ), which is related to the merger rate, and find a relationship between the mass and the peak frequency of the spectrum.
We generate GW spectra with a merger rate from ST, varying the redshifted total mass, and we find an agreement (within 10%) between redshifted ringdown frequency and the peak of the spectrum, see Table I. We obtain the same agreement if we use the merger rate from [23], suggesting once more that an estimate of the peak frequency can be used to constrain the average redshifted total mass of the population. This relationship, therefore, holds independently of the model used for the evolution of the pop III binaries.  Table I. Variation of the peak of GWB spectra with a change in redshifted total mass. We find agreement between the peak frequency and the redshifted ringdown frrequency.
Results-We simulate one year of observation time with the ET+2CE network, taking the CBC background from the ST catalogue. We find the best-fit models to the residual GWB that remains after subtracting the individual sources. A model selection study shows preference for other filters over a 2/3 PL, see Table II. The models with a broken power law shown in the last 3 rows are clearly favoured over a single power law model. However, we do not observe a great increase in Bayes factor for the smooth and triple BPL over just a BPL. Therefore, we conclude that a BPL filter is sufficient for a pop III GWB  Already with a varying-index PL search we deduce that the 2/3 filter is not appropriate, since the estimated power law index is α = −0.6, as observed in the corner plot in Fig. 3. The more intricate filters, however, fit the Ω GW spectrum well and capture the presence of the peak successfully. In order to understand the redshifted mass distribution of pop III creating the GWB, we investigate the peak frequency of the signal. We obtain a good estimate of the peak frequency using a BPL search filter as shown in Fig. 4, f peak = 15.4 Hz. The redshifted ringdown frequency that matches the peak frequency, f z ring = 15.4 Hz, corresponds to M z tot = 1076M . The ST redshifted mass distribution shown in the top panel of Fig. 2, gives an average redshifted total mass of the residual population, < M z tot >= 1121M . Therefore, our estimate of average M z tot agrees well with the true value. Finally, the estimated M z tot can be depicted as a curve in the redshift-intrinsic total mass plane since M z tot = (1+z)M tot , see Fig. 5. Note that we have included a 10% uncertainty error in matching of the ringdown and the peak frequency for consistency with our findings in Table I.
Conclusions-GWs emitted from CBCs formed by  Instead of the 2/3 power law, the GWB spectrum peaks in the low-frequency LIGO/Virgo range. A model selection study showed that we could observe the peak caused by the unresolved merger and ringdown GWs from pop III stars. With a good estimate of the peak frequency, we could even deduce the redshifted total mass distribution of the residual population. Taking data from the recent ST binary sources catalogue, we have demonstrated the effectiveness of our Bayesian analysis combined with implications of retrieved parameters. One should note that we have derived the relationship between peak frequency and redshifted total mass assuming equal-mass and non-spinning binaries. Future work should involve relaxing these assumptions. Additionally, one could further investigate the connection between the negative slope of the Ω GW spectrum and pop III properties such as merger rate and mass distribution.
For models with a break frequency we use a uniform prior for the first power law index between 2/3 and 5/3, since this represents the inspiral/merger regime of the CBC. Triple BPL has the third spectral index fixed to α 3 = 2/3 since we expect the inspiral phase of pop I/II signal to dominate at higher frequencies.
Intrinsic mass distribution-We study the relation between peak frequency of pop III GWB spectrum and the mass distribution of the sources. We fix the merger rate as a function of redshift to be the one of ST. For total mass, M tot = m 1 + m 2 varying between 10 and 90 M , . StarTrack merger rate evolved, equal-mass binaries. We find a relationship between peak frequency and total intrinsic mass of the merger. This is a model-dependent statement.
we generate Ω GW spectra and record the frequency at which the spectra are maximum. We then find a best-fit curve for the data, Hz, with f 0 = 12.8 Hz see Fig. 6. However, changing the merger rate to the one from [23], we find a different best fit curve, with f 0 = 53.7 Hz, implying that the intrinsic mass may be difficult to extract from the estimate of the spectrum peak. We find more promising results if we study the redshifted total mass and its relation to peak frequency, as described in the main text.