Matteo Martinelli
email: matteo.martinelli@inaf.it

Francesca Scarcella
email: francesca.scarcella@uam.es

Natalie B Hogg
email: natalie.hogg@ipht.fr

Bradley J Kavanagh
email: kavanagh@ifca.unican.es

Daniele Gaggero
email: daniele.gaggero@ific.uv.es

Pierre Fleury
email: pierre.fleury@ipht.fr

Prepared for

Primordial black holes (PBHs) are compact objects proposed to have formed in the early Universe from the collapse of small-scale over-densities. Their existence may be detected from the observation of gravitational waves (GWs) emitted by PBH mergers, if the signals can be distinguished from those produced by the merging of astrophysical black holes. In this work, we forecast the capability of the Einstein Telescope, a proposed third-generation GW observatory, to identify and measure the abundance of a subdominant population of distant PBHs, using the difference in the redshift evolution of the merger rate of the two populations as our discriminant. We carefully model the merger rates and generate realistic mock catalogues of the luminosity distances and errors that would be obtained from GW signals observed by the Einstein Telescope. We use two independent statistical methods to analyse the mock data, finding that, with our more powerful, likelihood-based method, PBH abundances as small as f PBH ≈ 7 × 10 -6 (f PBH ≈ 2 × 10 -6 ) would be distinguishable from f PBH = 0 at the level of 3σ with a one year (ten year) observing run of the Einstein Telescope. Our mock data generation code, darksirens, is fast, easily extendable and publicly available on GitLab ɭ.

Introduction

The idea that a population of primordial black holes (PBHs) may exist in the Universe and constitute a significant portion of the dark matter (DM) has been debated in the literature since the pioneering works of the late 1960s and early 1970s [START_REF] Zel'dovich | The Hypothesis of Cores Retarded during Expansion and the Hot Cosmological Model[END_REF][START_REF] Hawking | Gravitationally collapsed objects of very low mass[END_REF][START_REF] Chapline | Cosmological effects of primordial black holes[END_REF] (see e.g. [START_REF] Green | Primordial Black Holes as a dark matter candidate[END_REF] for a recent review). The discovery of such a population would have profound implications for fundamental physics, even if it were to represent only a subdominant portion of the DM. The formation of PBHs could reveal precious hints about inflationary and early Universe physics [START_REF] Dolgov | Baryon isocurvature fluctuations at small scales and baryonic dark matter[END_REF][START_REF] Jedamzik | Primordial black hole formation during the QCD epoch[END_REF][START_REF] Garcia-Bellido | Density perturbations and black hole formation in hybrid inflation[END_REF][START_REF] Khlopov | Primordial Black Holes[END_REF][START_REF] Garcia-Bellido | Primordial black holes from single field models of inflation[END_REF][START_REF] Ballesteros | Primordial black hole dark matter from single field inflation[END_REF][START_REF] Pattison | Quantum diffusion during inflation and primordial black holes[END_REF], and their subsequent evolution could impact structure formation [START_REF] Inman | Early structure formation in primordial black hole cosmologies[END_REF] and solve long-standing puzzles related to the early formation of supermassive black holes [START_REF] Volonteri | The origins of massive black holes[END_REF]. Moreover, their detection could help to exclude the presence of thermal weakly interacting massive particles and could strongly constrain models of particle physics which invoke new physics at the Weak scale [START_REF] Lacki | Primordial Black Holes as Dark Matter: Almost All or Almost Nothing[END_REF][START_REF] Adamek | WIMPs and stellar-mass primordial black holes are incompatible[END_REF][START_REF] Bertone | Primordial Black Holes as Silver Bullets for New Physics at the Weak Scale[END_REF][START_REF] Carr | Black holes and WIMPs: all or nothing or something else[END_REF].

PBHs whose masses lie between 1 M and 100 M are particularly promising candidates for study, given the now routine detection of gravitational wave (GW) signals from the mergers of compact objects in this mass range. Gravitational waves from the inspiral, merger and ringdown of a pair of BHs were first detected by the LIGO Scientific Collaboration (henceforth LIGO) and the Virgo Collaboration (henceforth Virgo) in 2015 [START_REF] Abbott | Observation of Gravitational Waves from a Binary Black Hole Merger[END_REF] and the detection of a binary neutron star merger with an electromagnetic counterpart signal was made in 2017 by LIGO and Virgo [START_REF] Abbott | GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral[END_REF]. In 2020, the KAGRA observatory in Japan [START_REF] Akutsu | KAGRA: 2.5 Generation Interferometric Gravitational Wave Detector[END_REF] joined the global network of GW detectors, reporting its first observations in conjunction with the GEO600 instrument in March 2022 [START_REF] Abbott | First joint observation by the underground gravitational-wave detector[END_REF].

The analysis of the combined information about the merger rate, mass and spin distribution of low redshift merger events collected by LIGO, Virgo and KAGRA (LVK) since 2015 have allowed upper limits on the PBH abundance f PBH = Ω PBH /Ω DM in the aforementioned mass window to be set [START_REF] Ali-Haïmoud | Merger rate of primordial black-hole binaries[END_REF][START_REF] Kavanagh | Merger rate of a subdominant population of primordial black holes[END_REF]. This has also triggered a debate about the possibility of identifying a subdominant population of PBHs on top of the (likely dominant) contribution of astrophysical black holes (ABHs) [START_REF] Franciolini | Quantifying the evidence for primordial black holes in LIGO/Virgo gravitational-wave data[END_REF][START_REF] Franciolini | How to assess the primordial origin of single gravitational-wave events with mass, spin, eccentricity, and deformability measurements[END_REF][START_REF] Hütsi | Two populations of LIGO-Virgo black holes[END_REF], motivated in particular by the detected substructure in the mass distribution of these objects [START_REF] Abbott | Population Properties of Compact Objects from the Second LIGO-Virgo Gravitational-Wave Transient Catalog[END_REF][START_REF] Abbott | The population of merging compact binaries inferred using gravitational waves through GWTC[END_REF]. Definitively detecting such a population of PBHs will become a concrete possibility with the advent of the third generation of GW observatories, such as the Einstein Telescope (ET) [START_REF] Maggiore | Science Case for the Einstein Telescope[END_REF] and Cosmic Explorer (CE) [START_REF] Reitze | Cosmic Explorer: The U.S. Contribution to Gravitational-Wave Astronomy beyond LIGO[END_REF]. 1Assessing the potential constraining power of the ET on PBHs is the focus of this work.

The most up-to-date configuration of the planned ET facility is known as ET-D [START_REF] Hild | Beyond the Second Generation of Laser-Interferometric Gravitational Wave Observatories[END_REF], which proposes a three-armed observatory consisting of three interferometers arranged in an equilateral triangle. In the ET-D configuration, each detector is in fact made up of a pair of detectors -one sensitive to a lower frequency range and the other to a higher frequency range -thus greatly increasing the overall sensitivity of the instrument with respect to the current generation of GW observatories. Furthermore, the triangular shape of the observatory will enable improved sky localisation of GW events [START_REF] Vitale | Parameter estimation for binary black holes with networks of third generation gravitational-wave detectors[END_REF]. Lastly, the interferometers and detectors will all be constructed underground, in an effort to reduce seismic noise [START_REF] Beker | Improving the sensitivity of future GW observatories in the 1-10 Hz band: Newtonian and seismic noise[END_REF].

The less noise in the detector, the smaller the amplitude of GWs -or strain -that can be detected. This noise is typically quantified by the strain amplitude spectral density (ASD), The strain amplitude spectral densities (ASD) of advanced LIGO (black) and the ET-D configuration of the Einstein Telescope (green), along with the characteristic strain of two PBH merger events, with component masses of 10 M (orange) and 30 M (purple) taking place at z = 1, that is, with a luminosity distance D ≈ 6.8 Gpc.

which we show for advanced LIGO2 (aLIGO) and ET-D 3 in Fig. 1. These ASD curves effectively show the lowest GW strain that can be detected and we show for comparison the characteristic strains of two PBH merger events with the example masses we use in the rest of this work: 10 M and 30 M . From this plot, we can see that ET-D will be sensitive to strains around two orders of magnitude smaller than what aLIGO can currently detect, and will also probe a much broader range of frequencies.

With all these factors taken together, the ET-D design is expected to yield many more observations of GWs at ever-greater cosmic distances (or redshifts) than the current generation of terrestrial detectors.

A crucial difference between ABH mergers and PBH mergers is the redshift evolution of their merger rates. While there is broad consensus on a steeply decreasing rate for ABHs beyond z 2 (see e.g. [START_REF] Dvorkin | Metallicity-constrained merger rates of binary black holes and the stochastic gravitational wave background[END_REF] and references therein), the PBH merger rate is expected instead to be a monotonically increasing function of redshift. In fact, a significant number of PBH binaries are expected to form by gravitational decoupling from the Hubble flow before matterradiation equality. As shown in [START_REF] Ali-Haïmoud | Merger rate of primordial black-hole binaries[END_REF][START_REF] Nakamura | Gravitational waves from coalescing black hole MACHO binaries[END_REF][START_REF] Sasaki | Primordial Black Hole Scenario for the Gravitational-Wave Event GW150914[END_REF], the distribution of the orbital parameters for these primordial binaries peaks at low values of both semi-major axis and angular momentum, and therefore low values of the merger time (compared to the Hubble time scale). Hence, even a small population of PBHs is expected to dominate the merger rate at early times. This idea naturally suggests a focus on high redshift GWs for the purposes of detecting PBHs [START_REF] Koushiappas | Maximum redshift of gravitational wave merger events[END_REF]. They may manifest as resolved individual events [START_REF] Chen | Distinguishing Primordial Black Holes from Astrophysical Black Holes by Einstein Telescope and Cosmic Explorer[END_REF][START_REF] De Luca | Bayesian Evidence for Both Astrophysical and Primordial Black Holes: Mapping the GWTC-2 Catalog to Third-Generation Detectors[END_REF][START_REF] De Luca | The Minimum Testable Abundance of Primordial Black Holes at Future Gravitational-Wave Detectors[END_REF][START_REF] Ng | On the single-event-based identification of primordial black hole mergers at cosmological distances[END_REF], but also as a measurable contribution to the stochastic GW background, as recently pointed out in [START_REF] Atal | Constraining changes in the merger history of (P)BH binaries with the stochastic gravitational wave background[END_REF][START_REF] Braglia | Tracking the origin of black holes with the stochastic gravitational wave background popcorn signal[END_REF]. In this article we focus on the former case: the prospects for detecting an anomalously large number of individual resolved events at large cosmic distances. Specifically, we propose (i) a careful theoretical modelling of both the PBH and ABH merger rates; and (ii) a realistic assessment of the ability of future observatories -specifically, the ET -to disentangle the two populations.

As far as modelling is concerned, a useful starting point is the formalism for the PBH merger rate developed in [START_REF] Ali-Haïmoud | Merger rate of primordial black-hole binaries[END_REF][START_REF] Nakamura | Gravitational waves from coalescing black hole MACHO binaries[END_REF][START_REF] Ioka | Black hole binary formation in the expanding universe: Three body problem approximation[END_REF], based on the assumption of an initially uniform and isotropic distribution of PBHs all having the same mass (in other words, a "monochromatic" mass function). In this scenario, neighbouring PBHs may decouple from the Hubble flow before matter-radiation equality, forming binaries which merge with an ever decreasing rate throughout the age of the Universe. We extend this formalism to include the effects of the early-time formation of PBH clusters [START_REF] Chisholm | Clustering of primordial black holes: basic results[END_REF][START_REF] Chisholm | Clustering of Primordial Black Holes. II. Evolution of Bound Systems[END_REF][START_REF] Raidal | Formation and Evolution of Primordial Black Hole Binaries in the Early Universe[END_REF] on the evolution of the binaries. The "background" to these PBH merger events is provided by the mergers of ABHs. The ABH merger rate traces the binary ABH birth rate, with some delay, which in turn is expected to trace the star formation rate, again with some delay [START_REF] Dvorkin | Metallicity-constrained merger rates of binary black holes and the stochastic gravitational wave background[END_REF]. A careful description of PBH and ABH merger rates as a function of redshift is therefore crucial in understanding how well future GW observatories can distinguish these two populations.

In this work, we make the conservative imposition of monochromatic mass functions for both the ABH and PBH populations -in other words, all the BHs in a given population have the same mass -and focus on the redshift dependence of the merger rates as the only discriminant between the two. We also focus on PBHs in the 1 M -100 M mass range, but we note that several other investigations of the capability of future observatories to detect PBHs selected the sub-solar mass range as an interesting avenue for study [START_REF] Bhatiani | Confirmation of Planet-mass Objects in Extragalactic Systems[END_REF][START_REF] Mukherjee | Prospects of discovering subsolar primordial black holes using the stochastic gravitational wave background from third-generation detectors[END_REF][START_REF] Pujolas | Prospects for probing gravitational waves from primordial black hole binaries[END_REF].

The uncertainty in the measured luminosity distances of events at high redshift is also expected to play an important role in separating the two populations of BHs. The authors of [START_REF] Ng | On the single-event-based identification of primordial black hole mergers at cosmological distances[END_REF] simulate the response of different detector networks located in Europe and the USA, and assess whether a single-event-based PBH identification can be unambiguous, if the redshift is large enough. The authors conclude that the typical redshift measurement is not precise enough to conclude with certainty that a single source is of primordial origin.

Motivated by these results, we present a framework that implements a statistically sound assessment of the capability of ET to:

• Detect an excess of merger events at high redshift with respect to the astrophysical expectation;

• Measure the f PBH associated with the detection, if present, or constrain f PBH to lie within some range.

To this aim, we generate a mock data set associated to the null hypothesis of f PBH = 0, and several data sets associated to different PBH fractions. We analyse these data using two different methods. First, we use an intuitive two-bin approach to determine the minimum abundance of PBHs that can be detected as a significant "excess" with respect to the astrophysical background, identifying the optimal binning in redshift. Next, we develop a parameter estimation pipeline, in which we estimate the posterior distribution of f PBH , by comparing the mock data generated for different fiducial values of f PBH with the theoretical distribution of event distances.

The article is structured as follows: in Sec. 2 we discuss our theoretical and phenomenological models for the redshift evolution of the merger rate of primordial and astrophysical BHs; in Sec. 3 we describe how we generate our mock GW catalogues; in Sec. 4 we present the results from the "cut-and-count" method, which assesses ET's ability to detect PBHs; in Sec. 5 we present the results from the likelihood-based method, which further analyses ET's ability to measure the PBH fraction; in Sec. 6 we discuss the implications and possible shortcomings of these results and finally in Sec. 7 we conclude. Appendices A to C provide more information on the merger rate suppression due to PBH clustering, the method for computing the signal-to-noise ratio of an event in our mock catalogues, and the effect of lensing on GW measurements. The code associated with the article, darksirens,4 is publicly available on GitLab ɭ.

Conventions, notation, cosmology. We assume a spatially flat ΛCDM cosmology throughout, using H 0 = 67.4 km s -1 Mpc -1 and Ω m = 0.315 as reported by the Planck collaboration [START_REF] Aghanim | Planck 2018 results. VI. Cosmological parameters[END_REF]. A three-bar equality sign (≡) indicates a definition; an upper-case P indicates a probability, while a lower-case p indicates a probability density function (PDF); bold symbols (x, θ, A, . . . ) stand for vectors or matrices. The luminosity distance is simply denoted with a D throughout, because we use no other notion of distance. We endeavour with all heroism to keep explicit factors of G and c where relevant.

Modelling the black hole merger rate

This section summarises our assumptions regarding the modelling of the merger rate of both ABHs and PBHs, which will then be used to produce mock catalogues of GWs.

Preliminary definitions

Correctly interpreting a population of observed merger events relies on our expectations for the distribution of mergers. We therefore need to accurately model this quantity. Physical models for the formation and mergers of BHs typically yield the merger rate density R(z); that is, the number of mergers per unit comoving volume and per unit time in the rest frame of the source, as a function of redshift z. However, from an observational point of view, the relevant quantity is rather the number of mergers that would be observed by an ideal detector per unit redshift and per unit time in the rest frame of the detector. This is what we shall refer to as the merger rate [START_REF] Abbott | Supplement: The Rate of Binary Black Hole Mergers Inferred from Advanced LIGO Observations Surrounding GW150914[END_REF][START_REF] Abbott | Search for intermediate mass black hole binaries in the first observing run of Advanced LIGO[END_REF],

R(z) = R(z) 1 + z dV c dz . (2.1)
The factor of (1 + z) -1 serves to convert a source-frame rate to a detector-frame rate, and V c denotes the comoving volume; dV c /dz then represents the comoving volume of a spherical shell between z and z + dz around the detector,

dV c dz = 4πr 2 (z) c H(z) , with r(z) ≡ z 0 dζ c H(ζ) , (2.2)
the comoving distance at redshift z. Both r(z) and H(z) are computed using the public cosmology solver CAMB [START_REF] Lewis | Efficient computation of CMB anisotropies in closed FRW models[END_REF][START_REF] Howlett | CMB power spectrum parameter degeneracies in the era of precision cosmology[END_REF]. 5 Given a merger rate R(z), the total theoretical number of events that would be detected by an infinitely sensitive instrument during an observation time

T obs reads Ntot = T obs zmax z min dz R(z) , (2.3)
where we limit ourselves to mergers in a redshift range z ∈ [z min , z max ]. Their redshift distribution, still in the case of an ideal detector, is given by

p(z) = T obs R(z) Ntot = R(z) zmax z min dζ R(ζ)
.

(2.4)

Importantly, not all of these mergers will have a signal-to-noise ratio (SNR) large enough to be detectable. We therefore define the selection function f det (z), which encodes the fraction of mergers that are detectable at a given redshift. We compute this explicitly for aLIGO and ET-D by sampling over binary orientations and sky positions, following the approach of Ref. [START_REF] Hall | gw-horizon-plot[END_REF], and requiring a threshold of SNR min = 8 for an event to qualify as detected. We illustrate these selection functions in Fig. 2 for a range of BH masses. With this, the rate of detectable events is given by

R det (z) = f det (z) R(z) . (2.5)
The expected number of detected events, Ndet , and their redshift distribution, p det (z), are defined analogously to Eqs. (2.3) and (2.4), respectively, by substituting R(z) with R det (z).

In the following subsections, we outline our computation of the merger rate density R(z) for ABH and PBH mergers. In Fig. 3, we show examples of these merger rate densities (left panel), along with the corresponding redshift distribution p(z) (right panel). These illustrate why high redshift mergers are so crucial to discovering and constraining a PBH population; above z ∼ 30 -40 the ABH rate quickly becomes negligible and mergers from a potential PBH population come to rapidly dominate.

Primordial black holes

PBH binaries are efficiently formed deep in the radiation era, in what is known as the early-time formation channel. Binaries can also be formed in dense environments through (2.4). This corresponds to the distribution of observed event redshifts for an ideal detector, ignoring the effects of a selection function f det (z).

gravitational capture or three-body interactions [START_REF] Bird | Did LIGO detect dark matter?[END_REF][START_REF] Franciolini | Primordial black hole mergers from three-body interactions[END_REF]. These late-time formation channels are expected to become increasingly important for large PBH abundances and low redshifts. In this work, we focus on high redshift observations and our analysis is restricted to f PBH 10 -3 (in agreement with current constraints, see Sec. 2.4). In this case, we can consider the merger rate to be dominated by the contribution of binaries formed in the early Universe. We compute it following the derivation of [START_REF] Ali-Haïmoud | Merger rate of primordial black-hole binaries[END_REF], which builds on several earlier works [START_REF] Nakamura | Gravitational waves from coalescing black hole MACHO binaries[END_REF][START_REF] Sasaki | Primordial Black Hole Scenario for the Gravitational-Wave Event GW150914[END_REF][START_REF] Ioka | Black hole binary formation in the expanding universe: Three body problem approximation[END_REF]. The PBHs are assumed to be initially Poisson distributed (we neglect any possible initial clustering, see e.g. [START_REF] Bringmann | Towards closing the window of primordial black holes as dark matter: The case of large clustering[END_REF][START_REF] Ballesteros | On the merger rate of primordial black holes: effects of nearest neighbours distribution and clustering[END_REF][START_REF] Young | Initial clustering and the primordial black hole merger rate[END_REF]) with negligible initial velocity dispersion. Under mutual gravitational attraction, PBH pairs decouple from the Hubble flow and form binary systems. Head-on collisions are avoided thanks to the small angular momentum provided by the presence of surrounding PBHs and smooth density perturbations. After formation, the binaries slowly lose energy through GW emission and eventually merge. The coalescence time (for small values of j) is given by [START_REF] Peters | Gravitational radiation and the motion of two point masses[END_REF] t merger = 3 170

c 5 G 3 M 3 PBH a 4 j 7 , (2.6)
where M PBH is the PBH mass, a is the binary's initial semi-major axis and j is the initial dimensionless angular momentum, related to the eccentricity e through j ≡ √ 1e 2 . As the binaries which merge within a Hubble time typically form early in the radiation era [START_REF] Ali-Haïmoud | Merger rate of primordial black-hole binaries[END_REF], we identify the coalescence time with the cosmic time at which the merger occurs.

Given Eq. (2.6), we can estimate the merger rate density from the initial distribution of the orbital parameters (a, j) as

R PBH [z(t)] = 1 2 n PBH da dj p(j|a)p(a) δ [t -t merger (a, j)] , (2.7)
where n PBH = f PBH ρ DM,0 /M PBH is the comoving number density of PBHs, and ρ DM,0 is the energy density of cold dark matter at redshift zero. We obtain distributions for the semi-major axis p(a) and for the angular momentum p(j|a) following [START_REF] Ali-Haïmoud | Merger rate of primordial black-hole binaries[END_REF][START_REF] Kavanagh | Merger rate of a subdominant population of primordial black holes[END_REF]. The latter distribution is computed taking into account the torques generated by all the surrounding PBHs and from the density perturbations in the matter fluid. Most binaries are born with very high eccentricities and small semi-major axes, which correspond to short times to merger, see Eq. (2.6). It follows that the PBH merger rate increases with redshift, with a larger number of binaries merging at early times, as shown in the left panel of Fig. 3.

In estimating the initial distribution of orbital parameters, we do not take into account the early-time disruption of binaries by neighbouring PBHs, which ultimately results in a suppression of the merger rate, as shown in [START_REF] Raidal | Formation and Evolution of Primordial Black Hole Binaries in the Early Universe[END_REF]. However, in this work we restrict ourselves to densities f PBH 10 -3 (see Sec. 2.4 for details) for which the effect is negligible [START_REF] Hütsi | Two populations of LIGO-Virgo black holes[END_REF].

The expression for the merger rate given in Eq. (2.7) assumes that the orbital parameters evolve exclusively through GW emission. In fact, a and j can be altered via a number of mechanisms. Apart from the early-time disruption by neighbours discussed above, one of these is accretion, which can modify the PBH mass and shrink the binary [START_REF] Caputo | Gravitational-wave detection and parameter estimation for accreting black-hole binaries and their electromagnetic counterpart[END_REF][START_REF] Ricotti | Effect of Primordial Black Holes on the Cosmic Microwave Background and Cosmological Parameter Estimates[END_REF][START_REF] Hütsi | Small-scale structure of primordial black hole dark matter and its implications for accretion[END_REF][START_REF] De Luca | The evolution of primordial black holes and their final observable spins[END_REF][START_REF] De Luca | Primordial Black Holes Confront LIGO/Virgo data: Current situation[END_REF]. The effect of accretion is particularly relevant for high PBH masses and results in an enhancement of the high mass tail of the mass function: we leave the study of this effect to a future work where extended mass functions will be considered. Furthermore, if PBHs exist side-by-side with particle DM, the latter forms mini-haloes around the PBHs. In this case, another perturbation to the binary can come from the dynamical friction that the mini-haloes induce on the PBHs when they approach each other [START_REF] Kavanagh | Merger rate of a subdominant population of primordial black holes[END_REF]. While the binaries are dramatically perturbed by this interaction, the shrinking of the binary compensates the increase of angular momentum in such way that the coalescence time is almost unaffected.

The simplified picture of isolated binaries which we have presented so far fails to capture an important aspect: under the action of gravity, PBHs form bound clusters, where complex N -body interactions take place. Various works have been dedicated to studying the effect of clustering on the evolution of binaries [START_REF] Raidal | Formation and Evolution of Primordial Black Hole Binaries in the Early Universe[END_REF][START_REF] Belotsky | Clusters of primordial black holes[END_REF][START_REF] Vaskonen | Lower bound on the primordial black hole merger rate[END_REF][START_REF] Jedamzik | Primordial Black Hole Dark Matter and the LIGO/Virgo observations[END_REF][START_REF] De Luca | The clustering evolution of primordial black holes[END_REF], showing that it can have a significant impact on the merger rate. In the most extreme scenarios, binaries can be easily disrupted in these dense environments. But even if the binary is not completely disrupted, it must be extremely eccentric in order to merge within a Hubble time; given the strong dependence of the time of merger on the angular momentum, t merger ∝ j 7 , even a small increase in j due to interactions within the cluster is sufficient to delay the merger beyond our time.

In this work, we quantify the suppression of the merger rate due to clustering following the semi-analytical modelling presented in [START_REF] Vaskonen | Lower bound on the primordial black hole merger rate[END_REF]. This calculation is based on estimating the fraction of clusters and sub-clusters that undergo core collapse following gravo-thermal instability. It is assumed that all binary systems within these structures do not contribute to the merger rate: they end up being perturbed in the high density cores in such a way that their coalescence time exceeds the age of the Universe. In this sense, the calculation can be considered an over-estimation of the effect (while it is very likely that a binary is perturbed in a dense environment, it is not necessarily so; furthermore, new binaries can be created in the cores). However, this may be partially compensated by the fact that the formalism does not consider perturbations to binaries in stable clusters (i.e. those not affected by core collapse). The calculation of the suppression factor is detailed in Appendix A.

In Fig. 4, we plot the PBH merger rate for different values of M PBH and f PBH , with and without the effects of clustering. The suppression factor increases going towards low redshifts, as clusters of increasingly larger size have sufficient time to undergo core collapse. As expected, it also increases with f PBH . As we explain in Appendix A, for small PBH fractions f PBH 10 -2 , the effect becomes negligible, as most clusters in this case form late in the Universe and do not have sufficient time to undergo gravo-thermal collapse before today. As the disruption in clusters becomes greater with time, it causes an overall enhancement of the slope of the function R(z), making small values at z = 0 compatible with larger overall PBH abundances. Given that current bounds from GWs lie in the range f PBH ∼ 10 -3 -10 -2 , where this effect may still be relevant, we update these to incorporate clustering in Sec. 2.4.

Astrophysical black holes

The ABH merger rate depends on the rate of formation of ABHs from stars via supernovae. We note that there are two channels through which ABH binaries can form: directly from binary stellar systems or from pairs of originally isolated black holes. In both scenarios, the ABH merger rate density R ABH can be expressed as [START_REF] Dvorkin | Metallicity-constrained merger rates of binary black holes and the stochastic gravitational wave background[END_REF][START_REF] Mukherjee | Can we distinguish astrophysical from primordial black holes via the stochastic gravitational wave background?[END_REF]

R ABH (t, M ABH ) = N ABH ∆tmax ∆t min d∆t p(∆t) R birth (t -∆t, M ABH ) , (2.8)
where R birth (t, M ABH ) is the birth rate density of the stellar remnants as a function of redshift and mass, and p(∆t) is the distribution function of the time delay between ABH formation and merger. In general, this time delay depends on the detailed rate of binary formation and on the orbital parameters of the resulting binaries (semi-major axis and eccentricity). In our model we follow the simplified approach of [START_REF] Dvorkin | Metallicity-constrained merger rates of binary black holes and the stochastic gravitational wave background[END_REF] and express this distribution as p(∆t) ∝ 1/∆t for ∆t min ≤ ∆t ≤ ∆t max with ∆t min = 50 Myr and ∆t max = H -1 0 . This expression is motivated by the high-resolution numerical simulations of binary BH formation via the evolution of isolated binary stars presented in [START_REF] Belczynski | The first gravitational-wave source from the isolated evolution of two 40-100 Msun stars[END_REF]. However, there are a number of uncertainties in the rate of ABH mergers [START_REF] Mapelli | The cosmic merger rate of neutron stars and black holes[END_REF][START_REF] Mapelli | The properties of merging black holes and neutron stars across cosmic time[END_REF][START_REF] Graziani | Cosmic archaeology with massive stellar black hole binaries[END_REF][START_REF] Santoliquido | The cosmic merger rate density evolution of compact binaries formed in young star clusters and in isolated binaries[END_REF][START_REF] Mapelli | Mass and Rate of Hierarchical Black Hole Mergers in Young, Globular and Nuclear Star Clusters[END_REF][START_REF] Mapelli | Formation Channels of Single and Binary Stellar-Mass Black Holes[END_REF][START_REF] Sedda | Population synthesis of black hole mergers with B-POP: the impact of dynamics, natal spins, and intermediate-mass black holes on the population of gravitational wave sources[END_REF]. We therefore include an overall normalising factor N ABH , which we fix based on the observed low redshift merger rate, as we will describe in Sec. 2.4.

We assume that the ABH birth rate density R birth (t, M ABH ) is proportional to the star formation rate (SFR) ψ SFR (z). Galaxy surveys performed in the ultraviolet band and in the far-infrared band can directly measure the instantaneous star formation rate density at different redshifts, with the former range of wavelengths being a direct tracer of short-lived massive stars, and the latter being a signature of UV light emitted by the same population of stars and subsequently absorbed and re-emitted by dust. In-depth analyses of these data allow for modelling of the cosmic history of the SFR, unambiguously highlighting a rising trend at low redshift up to a peak of very intense star formation in the range 1 ≤ z ≤ 2, followed by a fall-off at large redshift [START_REF] Madau | The Star Formation History of Field Galaxies[END_REF][START_REF] Madau | Cosmic Star Formation History[END_REF].

However, the high redshift behaviour of the SFR may be affected by significant uncertainties and biases, mainly due to dust obscuration and to the fact that early star formation took place in very faint galaxies which are typically missed in existing surveys [START_REF] Kistler | The Cosmic Star Formation Rate from the Faintest Galaxies in the Unobservable Universe[END_REF]. Therefore we also consider an alternative tracer of the high redshift SFR, via the rate of gamma ray bursts (GRBs). Several studies (see e.g. [START_REF] Vangioni | The Impact of Star Formation and Gamma-Ray Burst Rates at High Redshift on Cosmic Chemical Evolution and Reionization[END_REF] and references therein) based on this observable have indicated a higher SFR at high redshift compared to studies based on galaxy counts. With this method, the main uncertainty is the (model-dependent) relation between the rate of GRBs and the star formation rate itself.

We follow [START_REF] Nagamine | Abundance of damped Lyman-alpha absorbers in cosmological SPH simulations[END_REF] and adopt the following parametrisation for the SFR,

ψ SFR (z) = k a e b(z-zm) a -b + b e a(z-zm) .
(2.9)

The free parameters k, a, b, z m in this expression are usually fitted by comparison to existing catalogues of different tracers of star formation. Aiming to bracket the uncertainties associated with the astrophysical rate, we consider two different choices for these parameters. The reference set of parameters corresponds to the fiducial model quoted in [START_REF] Chen | Distinguishing Primordial Black Holes from Astrophysical Black Holes by Einstein Telescope and Cosmic Explorer[END_REF] based on a fit to observations of bright galaxies. The values are

k = 0.178 M yr -1 Mpc -3 , z m = 2, a = 2.37, b = 1.8.
We also consider a GRB-based fit as a maximal model, with the values taken from [START_REF] Vangioni | The Impact of Star Formation and Gamma-Ray Burst Rates at High Redshift on Cosmic Chemical Evolution and Reionization[END_REF]:

k = 0.146 M yr -1 Mpc -3 , z m = 1.72, a = 2.8, b = 2.46.
We stress that, unlike in the PBH case, the overall normalisation factor N ABH of the astrophysical rate is not determined a priori. Instead, we set this normalisation using the low redshift binary BH merger rate data provided by GW observatories currently in operation. Since in this work we are considering a hybrid scenario where a portion of the events is ascribed to a PBH population, each value of f PBH actually corresponds to a different normalisation factor for the rate associated to the ABH population, so that the integrated rate in the redshift range z ≤ 1 matches the observed one. The procedure to match the data naturally provides an upper limit on f PBH , set by the requirement not to overshoot the rate measurement. The details of the procedure, and the corresponding bound on f PBH , are described next.

Normalisation of the redshift distributions

Given the models for the redshift evolution of the ABH and PBH rate density described above, we now aim to: (i) obtain an expression for the actual merger rate by convolving the theoretical rate estimates with the detector space-time sensitivity; (ii) compare the predicted merger rate for each value of f PBH with the recent estimates provided by LVK. In this way, we will simultaneously set the normalisation of the astrophysical component for each value of f PBH , and obtain an upper limit on this quantity.

We compare the predicted PBH merger rate to the observational data from the second and third Gravitational-Wave Transient Catalogs (GWTC-2 and GWTC-3, produced by LIGO-Virgo and LVK respectively) for different PBH mass intervals and values of f PBH . In particular, we adopt the merger rate density estimates reported in Table 4 of [START_REF] Abbott | Population Properties of Compact Objects from the Second LIGO-Virgo Gravitational-Wave Transient Catalog[END_REF] for GWTC-2 and Table 4 of [START_REF] Abbott | The population of merging compact binaries inferred using gravitational waves through GWTC[END_REF] for GWTC-3. 6 These are reported assuming a non-evolving rate, while our models for the ABH and PBH merger rates evolve with redshift. We therefore compute the average merger rate density R , defined as

R = zmax 0 R det (z) dz zmax 0 f det (z) 1+z dVc dz dz . (2.10)
This corresponds to the non-evolving merger rate density which would give rise to the same number of expected detectable events as our redshift-dependent model, over the redshift range z ∈ [0, z max ]. For f det (z) we use the selection function for aLIGO shown in Fig. 2 and we fix z max = 2. 7 For each mass bin considered in these catalogues, we make the conservative choice to exclude the values of f PBH that would correspond to an average merger rate density R that exceeds the measured one by 2σ. For values of f PBH smaller than this bound, we adjust the ABH merger rate normalisation N ABH so that the sum of the averaged rates corresponds to the reported one, that is

R PBH (f PBH ) + R ABH (N ABH ) = R LVK obs . (2.11)
We have hence defined an allowed region for f PBH and properly set the normalisation of the astrophysical rate for each value of f PBH . We show the resulting new upper bounds we obtain on f PBH in Fig. 5. These bounds are compatible with previously reported bounds from GWTC-2 including the contribution of ABH mergers [START_REF] Hütsi | Two populations of LIGO-Virgo black holes[END_REF]. We find that the suppression due to late-time clustering, described in Sec. 2.2 and Appendix A, does not impact this upper limit. Its effect is relevant only for values of f PBH larger than O(10 -2 ), as shown in Fig. 4. Though the rate is significantly suppressed by clustering for large values of f PBH , the predicted PBH merger rate would still exceed the merger rate observed by LVK, and large values of f PBH remain excluded. Around our reported bound of f PBH 10 -3 , the suppression becomes negligible. Even so, we emphasise that the upper limit has to be taken cum grano salis, given the strong assumption of a monochromatic mass function. We cannot exclude scenarios featuring non-trivial mass functions and nonnegligible initial clustering, which may potentially evade the bound. The study of such scenarios are beyond the scope of the current work.

Generating mock GW catalogues

Having presented our merger rate calculation and the subsequent updated upper bounds on f PBH from current data, we are ready to turn to the core of this work: assessing the ability of the ET to disentangle ABH mergers from PBH mergers using their redshift distributions. This section presents our method to produce mock data for the ET, accounting for the instrumental response and other observational effects such as lensing of the signal. Since we focus on the information contained in the redshift distribution of the GW events -or, more accurately, their luminosity distance distribution which is the actual observable -in the remainder of this article, a data set D will refer to N det luminosity distance measurements D i with their uncertainty σ i , D = {(D i , σ i )} i=1,...,N det . Our mock generation code, darksirens ɭ, is publicly available and could be easily modified to include other observables.

10 0 10 1 10 2 M PBH [M ] 10 -4 10 -3 10 -2 10 -1 10 0 f PBH CMB Microlensing X-ray LIGO-Virgo [2019] LIGO-Virgo updated [GWTC-2]
LVK updated [GWTC-3] Figure 5. Updated LIGO-Virgo and LVK bounds on f PBH , based on the observed low redshift merger rate reported in recent GW transient catalogues. The various bins in M PBH correspond to the reported upper limits on the merger rate across different mass bins. For comparison, we also show bounds from the CMB due to PBH accretion [START_REF] Serpico | Cosmic microwave background bounds on primordial black holes including dark matter halo accretion[END_REF], microlensing constraints from the high redshift star Icarus [START_REF] Oguri | Understanding caustic crossings in giant arcs: characteristic scales, event rates, and constraints on compact dark matter[END_REF], and X-ray observations of the Milky Way [START_REF] Manshanden | Multi-wavelength astronomical searches for primordial black holes[END_REF]. We obtained these bounds from the PBHbounds repository [START_REF] Kavanagh | PBHbounds [Code v1.0[END_REF].

Parameters

The parameters that must be set to produce a mock catalogue of GW distance measurements can be divided into four categories:

• Cosmology. We assume a spatially flat homogeneous and isotropic ΛCDM cosmological background. The parameters to be set are the Hubble-Lemaître constant H 0 and the total matter density parameter Ω m .

• Primordial black holes. Following Sec. 2.2, the PBH population is characterised by the masses of the objects, M PBH , the fraction of DM made up of PBHs, f PBH , and whether or not clustering of PBHs is considered (we remark once again that the effect of clustering is negligible within the current implementation).

• Astrophysical black holes. Following the model outlined in Sec. 2.3, the ABH population is characterised by the masses of the objects M ABH and the SFR parameters z m , a and b which enter into Eq. (2.9).

• Specifications. This class of parameters allows us to set the observational time T obs of the survey and the SNR threshold, SNR min , which determines if a candidate event is detected or not. It also allows us to customise details of the mock generation, i.e. to specify the redshift range over which calculations are done (z ∈ [z min , z max ]) and whether or not to include the effect of lensing.

We remind the reader that we assume for simplicity a single mass for PBHs, M PBH , and a single mass for ABHs, M ABH . Recall also that the normalisation of the ABH merger rate N ABH is not set as an external parameter but is fixed by comparison with low redshift GW observations, as detailed in Sec. 2.4.

For the results presented in Secs. 4 and 5, we fix the parameters used to generate the mock data to the fiducial values listed in Table 1, unless otherwise specified in the text. We report the survey specifications assumed for the ET observations in Table 2.

Cosmology

ABH parameters PBH parameters

Ω m H 0 [km s -1 Mpc -1 ] M ABH [M ] z m a b M PBH [M ]

Sketch of the generation algorithm

Once the parameters are chosen, darksirens produces a mock catalogue D of GW distances with their uncertainties as they would be measured by ET. For f PBH = 0, the data is a mix of ABH and PBH mergers, which takes into account the potential clustering of the latter, the effect of lensing on distance measurements, and ET's instrumental uncertainties. We now explain exactly how the mock data are generated. The very first step consists in computing the total number N tot of mergers occurring in the redshift range [z min , z max ]. We draw a random number from a Poisson distribution whose mean, Ntot , is given by Eq. (2.3), in which we consider the total merger rate R(z) = R ABH (z) + R PBH (z) as described in Sec. 2. Each data point i = 1, . . . , N tot is then produced as follows:

1. Randomly draw the "true" redshift z i of the merger from the redshift distribution p(z)

given in Eq. (2.4), once again considering the total merger rate.

2. Convert z i into the "true" unlensed luminosity distance Di ≡ D(z i ) of the merger, using the background cosmological model and CAMB.

3. Compute the unlensed signal-to-noise ratio (SNR) ρi of the event. The position, polarisation and inclination of the event are drawn randomly and the SNR computed based on these quantities along with the specifications of the ET. The full details of this calculation are given in Appendix B.

4. Determine the weak gravitational lensing magnification µ i of the event by randomly drawing it from the theoretical PDF p(µ) -see Appendix C for details on lensing and its statistics. The lensing magnification enhances (or reduces) the SNR as

ρ i = √ µ i ρi ,
and reduces (or increases) the luminosity distance as Di = Di / √ µ i .

5. Determine whether the SNR is large enough for the event to be properly detected: if ρ i < 8, then the event is considered to be too faint to be a true GW candidate and is removed from the catalogue. The choice of discarding events with an SNR smaller than eight follows the approach taken by the LIGO collaboration as an estimate for the detection threshold [START_REF] Abbott | Sensitivity of the Advanced LIGO detectors at the beginning of gravitational wave astronomy[END_REF].

6. Compute the measured luminosity distance, D i = Di + ∆D i , where ∆D i represents the instrumental error on the measurement. The latter is drawn from a Gaussian distribution N (0, σ inst i ), with error inversely proportional to the SNR,

σ inst i = 2 Di /ρ i [94].
7. Compute the total uncertainty σ i on this data point as the quadratic sum of the instrumental and lensing uncertainties,

σ 2 i = σ inst i 2 + σ lens i 2 = 2 Di ρ i 2 + σ 2 κ (z i ) D2 i , (3.1)
where σ 2 κ (z) is the variance of the weak-lensing convergence -see Appendix C for details.

The end product is a catalogue of distance measurements with their uncertainties, D = {(D i , σ i )} i=1,...,N det , where N det ≤ N tot is the number of events that survive the SNR cut of step 5. We have checked that N det agrees on average with the theoretical expectation Ndet . Figure 6 shows an example of a catalogue that can be obtained using this approach. The figure shows the distance D i of the events that survive the SNR cut,8 together with their relative error, for both ABH (blue) and PBH (yellow). The data shown here are obtained setting f PBH = 10 -5 .

Distance uncertainty

Correctly modelling uncertainties in the luminosity distance of the GW events is crucial to disentangling high redshift and low redshift populations of mergers. A full analysis pipeline based on real GW data would provide as an output a probability distribution for the true luminosity distance P ( D) of the merger. Instead, in our catalogues, the distance uncertainty is described by only a single number σ i , given in Eq. (3.1). We therefore model the probability distribution for the luminosity distance as a Gaussian, given explicitly as:

p( D|D i ) = 1 √ 2πσ i exp - ( D -D i ) 2 2σ 2 i . (3.2)
In the statistical modelling which we will describe in Sec. 4 and Sec. 5, it is important to distinguish between P (D i | D) and p( D|D i ). In particular, some of the calculations presented in those sections are simplified if we identify the right hand side of Eq. (3.2) as p(D i | D), rather than p( D|D i ). However, this does not follow an intuitive definition of what is typically meant by the measurement error. We would typically consider that D i is our best estimate of the luminosity distance D, with σ i parametrising our uncertainty on D, matching the definition given in Eq. (3.2). We also find that interpreting the right hand side of Eq. (3.2) as p(D i | D) leads to inconsistencies. Consider the case where both D i and σ i are large, in which case the right hand side of Eq. (3.2) is a relatively flat function of D. If we interpret this as p(D i | D), it would imply that a wide range of D values would all give rise to the same observed D i with similar probabilities. However, we know that the measurement error becomes smaller as we reduce D, in which case large D i values should be very unlikely. The only consistent definition for the error is therefore given by Eq. (3.2).

We will now present two different analysis methods of the mock data we have generated: an intuitive binned approach, and a more complete Bayesian analysis of the data, in Sec. 4 and Sec. 5 respectively. With the first method we aim to estimate the detectability threshold of a PBH population, while with the second method we also assess the capability of the ET to measure the PBH fraction. We refer to these methods as "cut-and-count" and "likelihoodbased" respectively. 1 and Table 2, with f PBH = 10 -5 . The blue dots show the observed ABH events, while the yellow dots represent the observed PBH events. At a fixed distance D, the typical uncertainty is smaller for PBHs than for ABHs due to our choice of a larger PBH mass (M PBH = 10 M , M ABH = 7 M ). The dotted horizontal line shows the level of instrumental uncertainty corresponding to our SNR cut of ρ i = 8. The events that lie above this line still satisfy the SNR condition we impose, but receive an extra contribution to the error from lensing.

Detecting PBHs: the cut-and-count method

In this section, we present the "cut-and-count" method which we use to determine the smallest PBH fraction that could be detected by the ET. The idea is simple and intuitive: since we expect no ABHs to be formed beyond some high redshift (as there is a significant delay between the beginning of the Universe, the birth and death of the first stars and hence the formation and merging of ABHs), any sufficiently high redshift GW event produced by a BH merger should be the result of merging primordial BHs.

In practice, however, large distances are also the most uncertain ones, since the GW signal is typically much fainter than closer events. Hence the question we aim to answer is better phrased as: what is the lowest value of f PBH that would still produce a sufficient number of large-distance events so that they would be statistically distinguishable from f PBH = 0 when observing with the ET? We hence cut the data set into two subsets containing the small-and large-distance events, count the number of events in the large-distance subset, and compare to the expected number for f PBH = 0. This method only exploits a fraction of the available information in the data set, but it is nevertheless a natural starting point.

Description of the method

Let D * ≡ D(z * ) be an arbitrary distance threshold (z * is the corresponding background redshift threshold). Given a data set D, we divide it into two subsets: a small-distance subset D ≤ on the one hand, made of the events with distances beneath the threshold (D ≤ D * ); and a large-distance subset D > on the other hand, with events beyond the threshold (D > D * ). We call N > the cardinal of D > , i.e., the number of events above the threshold D * ; the larger f PBH , the larger the expected N > .

To be more specific, the computation of N > does not simply consist in counting the number of events whose best-fit distance is above the threshold. Due the measurement errors, it may happen that an event truly lies beneath the distance threshold ( Di < D * ) but is actually measured beyond it (D i > D * ), or vice versa. In order to account for this, we calculate the probability that the true distance lies beyond the arbitrary distance threshold as

P i = 1 √ 2πσ i ∞ D * d Di exp - ( Di -D i ) 2 2σ 2 i , (4.1)
using the definition of the uncertainty from Eq. (3.2). We can therefore obtain the value of N > by summing over all the events the probability of falling within the bin being considered,

N > ≡ N det i=1 P i . (4.
2)

The uncertainty on N > is twofold. First, since the data set is discrete, it inevitably comes with a Poisson uncertainty with variance

σ 2 P = N> , (4.3)
where N> is the theoretical expectation value of N > . 9 On the other hand, the number N > of large-distance events may be seen as the sum of N det independent Bernoulli variables, the i th one having a probability P i of being equal to 1, and a probability 1 -P i of being equal to zero. This observational contribution to the variance of N > is the sum of the individual Bernoulli variances,

σ 2 B = N det i=1 P i (1 -P i ) . (4.4)
Summing the two sources of uncertainty in quadrature yields the total error on N > ,

σ > = σ 2 P + σ 2 B . (4.5)
As an example, we consider the null case of a data set D 0 containing no PBH (f PBH → 0). Figure 7 shows the expected value of N > as a function of the redshift threshold z * , to which the N > of an observed data set should be compared. The error bars of that figure are computed following the approach described aobve, and they represent the uncertainty σ > on N > .

Optimal distance threshold and smallest detectable PBH fraction

Let us now determine what is the smallest detectable PBH fraction f PBH that could be detected using the cut-and-count method. This will also be the occasion to determine the optimal choice for the arbitrary distance threshold, that is, the value of D * leading to the maximal sensitivity in f PBH .

For that purpose, we generate mock data sets D f PBH for 100 values of f PBH ranging from 10 -6 to 10 -2 . For each of these data sets, we apply the cut-and-count method for 10 different values of the redshift threshold z * , taken as equispaced in the range [START_REF] Dolgov | Baryon isocurvature fluctuations at small scales and baryonic dark matter[END_REF][START_REF] Raidal | Formation and Evolution of Primordial Black Hole Binaries in the Early Universe[END_REF]. This yields a number N > (D f PBH , z * ) for each mock data set and each value of z * . Those numbers come with their own uncertainty σ > (D f PBH , z * ), computed as described in Sec. 4.1 for the null case D 0 . We estimate the statistical significance of the detection of a non-zero f PBH as

S(D f PBH , z * ) ≡ |N > (D f PBH , z * ) -N > (D 0 , z * )| σ 2 > (D f PBH , z * ) + σ 2 > (D 0 , z * ) , (4.6)
which represents the "number of sigmas" with which we could claim a detection.

In the following, we consider a "significant" detection to be one made at 3σ; in other words, for each value of z * , the smallest detectable f PBH is determined by finding the data set D f PBH such that S(D f PBH , z * ) > 3. In the top left panel of Fig. 8, both curves (for M PBH = 10 M in red and M PBH = 30 M in yellow) show the evolution of the smallest detectable f PBH with z * . The optimal value for the redshift threshold, leading to the best sensitivity in f PBH , is found to be z * ≈ 10, corresponding to a distance of D * ≈ 106 Gpc.

The existence of an optimal value for z * arises from a trade-off between two effects: on the one hand, increasing z * reduces the contamination due to ABH events in the high redshift subset; on the other hand, reducing z * increases the number of events in that high redshift subset and hence the significance of the detection. Our results suggests that choosing z * < 10 implies too much ABH contamination to easily assess whether or not PBHs are present in the data set, while z * > 10 worsens the statistics and only provides a detection for large values of f PBH . These conclusions are valid for both M PBH = 10 M and M PBH = 30 M , although the exact value of the optimal z * does change when the mass of the PBHs is changed.

We expect the SFR model to affect the results, as it determines how fast the probability of ABH mergers vanishes with redshift, see Fig. 3. We verify the impact of this uncertainty by repeating our analysis using the merger rate obtained assuming the GRB SFR model. In the top right panel of Fig. 8, we compare the GRB SFR results with the baseline case (M PBH = 10 M , fiducial SFR and T obs = 1 yr). We can see how in the GRB case, the optimal threshold redshift z * increases due to the ABH merger probability being non-vanishing up to higher redshifts than in the fiducial case. We also explore the impact of the observation time of the survey on these results; we compare our baseline results with an extended survey time for ET, setting T obs = 10 yr. As can be seen in the bottom panel of Fig. 8, increasing the time of the survey generally leads to a lower detection threshold for f PBH , while also decreasing the dependency of such a threshold on the choice of the binning strategy. With a greater survey time, more high redshift events are observed, meaning that a large detection significance is still achievable even as the high redshift bin is moved to larger values of z * .

Overall, we find that a fraction f PBH ≈ 2 × 10 -5 is the lowest value detectable by the ET with this method for T obs = 1 year, with a weak dependency on the mass of the progenitor systems and on the SFR model chosen. In our baseline case, for the optimal z * ≈ 10, such a fraction of PBH corresponds to N > = 16 ± 4.6, while the no PBH case yields N > = 1 ± 1.7, a result that puts this f PBH over the 3σ threshold we consider for detection.

Figure 9 shows more details for the results obtained in our baseline case, when taking z * = 10. The left panel shows in red the results for N > for the sampled values of f PBH , together with their error bars and the uncertainty on N > in the no PBH case (grey band). In this plot, we also show the trend of the detection significance S with f PBH (blue line). The right panel shows the counts in the two redshift bins for f PBH = 2 × 10 -5 in red, i.e. the first value for which the counts in the large-distance bin result in S ≥ 3. Together with the counts we also show the distance and its error for the events contained in this mock data set, as a function of the event redshifts (which would not be observed in reality).

The conclusions we have found for our baseline case using this cut-and-count approach are compatible with those that can be obtained using the approach of [START_REF] Ng | On the single-event-based identification of primordial black hole mergers at cosmological distances[END_REF], where the discriminant for the presence of PBH in the observed data set is the observation of events at z > 30, where the contribution of ABHs is negligible. With our simulated data we indeed find that f PBH ≈ 2 × 10 -5 is the lowest value in the baseline case for which we have at least one event above z = 30 with 99.7% confidence level. 1 andTable 2) (red points). The grey band shows the uncertainty for N > in the no-PBH case. The blue line shows instead the trend of the detection significance for a departure from the no-PBH case as a function of f PBH . Right: Counts in the small-and large-distance bins are shown by the red bar plot for the baseline case, with z * = 10 and f PBH = 2 × 10 -5 , together with the errors we determine on the value of the counts (left y-axis). In blue, the distances and their errors for all the events contained in the data set are shown as a function of the true redshift z true (right y-axis), while the top x-axis shows the redshift z obs that could potentially be inferred from the distance measurements.

5 Quantifying the PBH fraction: a likelihood-based method Having followed our simple cut-and-count method to quantify the lowest value of f PBH that would be statistically distinguishable from f PBH = 0 (i.e. the ability of the ET to detect PBHs), we now utilise a more powerful, likelihood-based method to assess the ET's potential capacity to both detect PBHs and measure the f PBH associated with a detection. The analysis that follows assumes fixed values of the cosmological parameters Ω m and H 0 , considering them as priors set from independent experiments. We thus neglect the potential degeneracies that could exists between f PBH and the cosmological model. This approach is justified by the fact that even future GW experiments are expected to yield much looser cosmological constraints than established probes such as Planck.

Description of the method

We describe an unbinned likelihood-based approach to comparing the data to a given model for the merger rate, which may include a contribution from both ABH and PBH mergers. To do this, we must compute the likelihood L(f PBH ), which is the probability to observe the data set D given the PBH fraction f PBH .

We write the probability distribution for the true luminosity distances of ABH and PBH mergers as p ABH ( D) and p PBH ( D|f PBH ) respectively, where we have indicated explicitly that the distribution for PBH mergers depends on f PBH . The probability that a BH merger has a true luminosity distance in the range [ D, D + d D] is given by

p( D|f PBH ) d D = N ABH Nobs p ABH ( D) d D + N PBH Nobs p PBH ( D|f PBH ) d D , (5.1)
where Nobs (f PBH ) = N ABH + N PBH (f PBH ) is the total expected number GW events and N ABH and N PBH are the expected numbers of ABH and PBH mergers respectively. These are the numbers of events expected to be observed above the SNR threshold, taking into account the detection efficiency of the observatory, as detailed in Sec. 2.1. We do not observe the true distance D but instead an estimate of the luminosity distance D. The probability that we observe an event with estimated luminosity distance D can be written as

p(D|f PBH ) = p(D| D) p( D|f PBH ) d D .
(5.2)

The distribution of "true" luminosity distance p( D|f PBH ) is given by the theoretical expectation in Eq. (5.1). The term p(D| D) can be obtained from the measurement uncertainty in Eq. (3.2) via Bayes' theorem [START_REF] Bayes | An essay toward solving a problem in the doctrine of chances[END_REF],

p(D| D) = p( D|D) p( D) p(D) . (5.3)
Here, p(D) and p( D) are the overall probability distributions of D and D, where by overall we mean that they are marginalised over all theory parameters (which in our case consists of only f PBH ). The term p(D) enters as an overall normalisation which can be pulled out of the integral in Eq. (5.2) and does not depend on the theory parameters. This can therefore be safely neglected in the likelihood. The final task is then to compute p( D), for which we need a prior on the PBH fraction Pr(f PBH ):

p( D) = 1 0 p( D|f PBH ) Pr(f PBH ) df PBH .
(5.4)

We assume an uninformative log-flat prior, Pr(f

PBH ) ∝ 1/f PBH , with f PBH ∈ [10 -9 , 10 -3 ].
For a sample of N obs observed merger events, the likelihood can then be written as

L(D|f PBH ) = Nobs (f PBH ) N obs e -Nobs (f PBH ) N obs ! × i=1,N obs p(D i |f PBH ) . (5.5)
The first term is the Poisson probability to observe N obs merger events, given that we expect to observe Nobs (f PBH ). The second term accounts for the contribution of each observed event to the likelihood, where D i are the estimated luminosity distances. Combining the results above, the probability of observing a merger at an estimated distance D i is given by

p(D i |f PBH ) ∝ p( Di |D i ) p( Di ) p( Di |f PBH ) d Di , (5.6)
which depends on the uncertainty σ i through p( Di |D i ) in Eq. (3.2).

We then adopt a Bayesian approach in order to construct the posterior distribution function p(f PBH |D) and thus determine projected constraints on the PBH fraction,

p(f PBH |D) ∝ L(D|f PBH )Pr(f PBH ) .
(5.7)

We sample log 10 f PBH from a uniform prior log 10 f PBH ∈ [-9, -3] using the public cosmological sampling code Cobaya [START_REF] Torrado | Cobaya: Code for Bayesian Analysis of hierarchical physical models[END_REF], while fixing all other parameters to the values used to generate mock data. We analyse the results using GetDist [START_REF] Lewis | GetDist: a Python package for analysing Monte Carlo samples[END_REF] to obtain the bounds achievable on f PBH using this approach.

Future constraints on f PBH

In this section we discuss the bounds on f PBH that ET can potentially obtain, as a function of the "true" value of this parameter, using the likelihood-based method. In order to estimate the constraints, we generate 10 mock data sets with fiducial values f fid PBH logarithmically distributed in the range 10 -6 , 10 -4 . For each fiducial value of f fid PBH , we compute the likelihood associated to each data set, and obtain the posterior distribution function applying Eq. (5.7). We define the mean of the distribution as the "measured value" f meas PBH . We show the results of this analysis in Fig. 10. In the top left panel, we visualise the 68% and 99.7% confidence level bounds (i.e. 1σ and 3σ) on the measured PBH fraction for each of the chosen f fid PBH in our baseline settings, together with the mean value obtained for f meas PBH . We interpolate between the results obtained for our ten values of f fid PBH , in order to visualise the trend of these bounds.

We notice how the qualitative behaviour of these results is similar to that of the previous method; the analysis highlights how for low fiducial values (f PBH 10 -5 ) the ET is not able to detect the presence of PBHs, and only an upper bound can be placed on f PBH . Choosing a threshold of 3σ for a non-vanishing value of f meas PBH to be considered a detection, as we did in Sec. 4, we find that ET will be able to detect the presence of PBHs for f fid PBH ≈ 7 × 10 -6 . Such a value is the result of interpolating between the 99.7% confidence level lower bounds obtained for the different data sets, and then finding for which value of f fid PBH this function would exclude the no PBHs case (set to be at f PBH = 10 -7 ). Such a value is roughly three times lower than what we found with the cut-and-count method, highlighting how using the full amount of information present in the data set helps to boost the survey sensitivity.

For higher values of f fid PBH , a detection is possible and the constraining power of ET increases with higher fiducial values. In the right panel of Fig. 10, we show the posterior distribution of f meas PBH for three different f fid PBH ; here we highlight how for low f fid PBH , the posterior distribution is extremely flat at low f PBH , as the ET cannot distinguish between fractions of PBHs that produce a very low number of events. Moving towards higher f fid PBH , the posterior becomes increasingly peaked, showing how these PBH abundances could be measured with very high precision.

These results highlight how, should the Universe contain a high enough number of PBHs, the method proposed in this work will be able to constrain the value of f PBH . For example, considering f fid PBH = 1.3 × 10 -5 , i.e. the first generated data set above the detection threshold of f PBH ≈ 7 × 10 -6 found above, we find that one can obtain a measurement of f meas PBH with a precision of 35 % at 68 % confidence level.

Finally, we compare the baseline results with what can be achieved increasing the observation time to T obs = 10 yrs. This comparison is depicted in the bottom panel of Fig. 10, where the baseline bound is shown in red, while the result for the extended observation time is shown in navy. Our results confirm that extending the survey time allows the detection of lower fractions of PBH; specifically, the smallest detectable f fid PBH in the ten year case is about five times smaller than in the one year case. For a given detectable value of f PBH , the uncertainty on its measurement is also reduced; for instance, with f fid PBH = 1.3 × 10 -5 we find that f meas PBH can be obtained with a precision of 11 % with 68 % confidence in the ten year case.

Discussion

The rationale of the two statistical frameworks presented here is to present complementary approaches, and emphasise different aspects of the analysis of GWs from BH merger events as a tool to detect and measure a PBH population. In particular, the cut-and-count method is especially designed to intuitively capture and visualise the essential aspect of the problem, namely the potential discovery of an anomalous excess of high redshift events beyond some reference distance. On the other hand, the likelihood-based approach exploits the capability of a Bayesian framework while using all the available information encoded in the redshift distribution of the events. In both cases, the characterisation of the instrument response plays a crucial role, and the main results we obtained revealed the capacity of the ET to disentangle astrophysical and primordial BH merger events on solid statistical grounds. However, when searching for signatures of new physics against a highly uncertain astrophysical background, a key question is whether such signatures remain undetected simply due to the limitations of instrumental sensitivity or instead due to poor understanding of either the background or the signal itself, and the systematic uncertainties associated with each. In the search for PBHs we are considering here, the modelling of both the background and the signal poses a clear challenge.

Beginning with the signal, our model relies on a set of simplifying assumptions which includes: (i) a statistically uniform initial spatial distribution of PBHs (i.e. no initial clustering); and (ii) a monochromatic mass function (i.e. all the PBHs have the same mass). Regarding (i), as pointed out in [START_REF] De Luca | The Minimum Testable Abundance of Primordial Black Holes at Future Gravitational-Wave Detectors[END_REF], initial clustering may significantly change the picture: the authors of that study aimed to find the minimum PBH abundance testable by future detectors and concluded that abundances as small as f PBH ∼ 10 -10 can be probed, if PBHs are highly clustered at formation; in the case of no initial clustering, they estimated that values of f PBH ∼ 10 -5 can provide at least one event per year at the ET. Our results based on the mock data generation and careful statistical analysis are compatible with the latter estimate. We leave the assessment of the role of the initial clustering, as well as a more refined analysis of the late-time clustering, to future work.

Regarding (ii), realistic PBH production scenarios generally predict an extended mass function rather than a monochromatic one as assumed here. Assessing the role of the mass function requires us to understand the role of the PBH mass itself in the merger rate. On the one hand, the larger the chirp mass of a PBH binary, the higher the amplitude of the corresponding GW signal, which may result in a higher SNR, and hence an increased detectability of the PBH population. However, changing the chirp mass also affects the frequency range of the GW emitted by the merger. Our choice for the PBH mass, M PBH = O(10 M ), roughly coincides with an optimal detectability of high-redshift events, given the ET noise curve. 10 On the other hand, for a fixed f PBH , the larger the PBH mass, the fewer PBHs, and hence the rarer the merger events. Therefore, if we consider a broad mass function, we expect that a high-mass tail would be associated with rare but (potentially) high-SNR events, while a low-mass tail would correspond to abundant but low-SNR events. Overall, for a relatively narrow mass function we do not expect our main conclusions to dramatically change, but we leave the study of more complex mass functions for future work.

Related to the above point, we chose to keep the masses of the ABH and PBH populations fixed in the statistical analysis. We have not used any information about the inferred masses of the mergers in this work, instead focusing on how the distance information can disentangle the two populations. Of course, including mass information is likely to enhance our ability to detect a PBH population, since detecting a deviation from the conventional ABH mass distribution model would be a hint at a possible population of PBHs, but that would require further assumptions about the ABH mass distribution.

The possibility to detect a subdominant population of PBHs by analysing the mass function of binary BH mergers has progressively gained momentum with the detection of many BH merger events with masses around 30 M [START_REF] Abbott | Population Properties of Compact Objects from the Second LIGO-Virgo Gravitational-Wave Transient Catalog[END_REF][START_REF] Abbott | The population of merging compact binaries inferred using gravitational waves through GWTC[END_REF]. For example, the authors of [START_REF] De Luca | Bayesian Evidence for Both Astrophysical and Primordial Black Holes: Mapping the GWTC-2 Catalog to Third-Generation Detectors[END_REF] claim that there is decisive evidence in favour of a two-population model from these data, with a second, subdominant, population peaked around 30 M . This population displays a significant high-mass tail, hence its presence is further supported by the detection of a BH with mass in the pair-instability gap region [START_REF] Abbott | GW190521: A Binary Black Hole Merger with a Total Mass of 150M[END_REF] (several astrophysical models of BH formation in the pair-instability mass gap also exist, see e.g. [START_REF] Abbott | Properties and Astrophysical Implications of the 150 M Binary Black Hole Merger GW190521[END_REF]).

Motivated by these considerations, we emphasise that an exploration of the impact of broad mass distributions for the PBH population (possibly with multiple peaks on different mass scales [START_REF] Carr | Cosmic conundra explained by thermal history and primordial black holes[END_REF]) on the PBH detection prospects with the ET is an important future goal in this research line. This would require a more careful treatment of late-time clustering than that presented in this work, given the complex phenomenology that may arise (for instance, mass segregation in clusters [START_REF] Trashorras | The clustering dynamics of primordial black boles in N -body simulations[END_REF]).

We have assumed throughout this work that PBH binaries formed in the early Universe dominate the merger rate. Dense environments, such as PBH clusters, can potentially lead to the disruption of early Universe binaries. We have considered the effect of late-time clustering (see Appendix A), which seems to be negligible for the small PBH abundances considered in this work; the impact of initial PBH clustering is left for future work. PBH binaries could also be disrupted by interactions with other objects; while it is unlikely that they will be disrupted by stellar encounters in the disks of typical galaxies [START_REF] Sasaki | Primordial Black Hole Scenario for the Gravitational-Wave Event GW150914[END_REF], it seems feasible that they may be disrupted in sufficiently dense environments, such as close to the centre of the host galaxy. This is unlikely to affect the high redshift signals we study here, with the first haloes having little effect on early-formed binaries [START_REF] Ali-Haïmoud | Merger rate of primordial black-hole binaries[END_REF]. On the other hand, dense environments can lead to the formation of a population of late-time PBH binaries. This is typically expected to be subdominant [START_REF] Ali-Haïmoud | Merger rate of primordial black-hole binaries[END_REF][START_REF] Bird | Did LIGO detect dark matter?[END_REF], in particular for the small values of f PBH considered in this work. These late-time binaries may be more difficult to distinguish from ABH ones, as they are expected to form only at low redshift, in which case alternative approaches are necessary, such as a measurements of the clustering bias [START_REF] Herrera | Learning How to Surf: Reconstructing the Propagation and Origin of Gravitational Waves with Gaussian Processes[END_REF]. However, as long as the population of early-time binaries remains relatively unperturbed, high redshift observations would be able to probe it, and hence detect the presence of PBHs, even if this population were to be subdominant with respect to the late-time one. Regarding the measurability discussion, a sizable population of late-time binaries merging at relatively high redshifts could lead us to overestimate the abundance of PBHs. In neglecting this, we rely on the expectation that the merger rate due to such binaries is only relevant at low redshifts and for large PBH abundances.

There are also a number of ways in which the modelling of the astrophysical background could alter the results presented here. In particular, we emphasise that our model for the ABH evolution, characterised by a peak at z ∼ 2 and a monotonic decrease at higher redshifts, is designed to capture the rate associated to binary systems made of second-and third-generation stars (usually called Population I and II). These systems may form via different channels, such as binary stellar evolution in galactic fields, or dynamical formation through multi-body interactions in star clusters.

However, a significant contribution to the binary BH merger rate, including a peak at high redshift, z ∼ O [START_REF] Ballesteros | Primordial black hole dark matter from single field inflation[END_REF], can be expected from the first-generation (Population III) stars, which formed out of the pristine gas left over after cosmological nucleosynthesis, and generated the first heavy elements in the Universe [START_REF] Vangioni | The Impact of Star Formation and Gamma-Ray Burst Rates at High Redshift on Cosmic Chemical Evolution and Reionization[END_REF][START_REF] De Souza | Populations III.1 and III.2 gamma-ray bursts: constraints on the event rate for future radio and X-ray surveys[END_REF]. Taking into account this additional component would certainly require, again, to address the mass dependence in more detail. In fact, this population would be modelled as a high (or intermediate) mass model of star formation at high redshift, to be added to the background model [START_REF] Ng | Probing multiple populations of compact binaries with third-generation gravitational-wave detectors[END_REF]. An even more detailed treatment would imply a marginalisation over the parameters that describe the metallicity dependence of the SFR, treated as nuisance parameters. We leave this more detailed treatment to future studies, and emphasise that the numerical and statistical tools presented here are the ideal framework to consistently address uncertainties in both the signal and background models.

As we were finalising this work, Ng et al. [START_REF] Ng | Constraining high-redshift stellar-mass primordial black holes with next-generation ground-based gravitational-wave detectors[END_REF] proposed an independent analysis on the same topic. Reference [START_REF] Ng | Constraining high-redshift stellar-mass primordial black holes with next-generation ground-based gravitational-wave detectors[END_REF] explores how a network consisting of two third generation detectors (ET and CE) can be used to distinguish between ABH and PBH populations using high redshift mergers. Their analysis differs from ours in that they consider exclusively population III ABHs and restrict their analysis to events above z > 8, as well as working directly with measurement errors on the event redshift, rather than on the luminosity distance, as we do here. Their results are broadly consistent with ours, highlighting that a PBH population with f PBH ∼ O(10 -5 ) should be within reach of third generation detectors. We find a one year sensitivity estimate for f PBH roughly a factor of three stronger, which may be partly due to differences in the distance uncertainties and population models assumed. This emphasises the importance of careful modelling of these details in quantifying the PBH detection and measurement potential of future GW detectors.

Conclusions

In this article we have assessed the capability of the forthcoming Einstein Telescope (ET) gravitational wave (GW) observatory to detect and measure a subdominant population of primordial black hole (PBH) mergers using a novel statistical framework.

We have described a procedure that computes the redshift evolution of both the expected background associated with the astrophysical black hole (ABH) merger events, and the signal under consideration associated with the PBH population (assumed to be characterised by masses of the same order of magnitude). In our modelling, we paid particular attention to the impact of late-time clustering of PBHs. We set the relative normalisation of the two contributions by comparison with the latest (low redshift) data released by the LIGO, Virgo and KAGRA collaborations. This procedure naturally provided an updated upper limit on the fraction of dark matter in the form of PBHs (quantified using the parameter f PBH ), based on the third Gravitational Wave Transient Catalog (GWTC-3) -see Fig. 5.

The key feature of our merger rate models is the different behaviour at high redshift of signal (i.e. PBH mergers) and background (i.e. ABH mergers), with the former monotonically increasing with increasing distance, and the latter steadily decreasing. Motivated by this qualitative aspect, and taking into account the expected high sensitivity of ET, we presented two methods to assess the potential of discovering a distant PBH population, which both accurately take into account the experimental sensitivity. Both methods are based on the generation and analysis of a set of mock data catalogues, associated to the null hypothesis (no PBHs) and to different values of the fraction of DM in the form of PBHs.

The first method (dubbed "cut-and-count") relies on a two-redshift-bin data analysis strategy that highlights the PBH contribution as a significant excess in the high redshift bin.

A Merger rate suppression due to PBH clustering

In this appendix, we discuss the suppression of the PBH merger rate due to PBH clustering. We follow the numerical implementation given in Ref. [START_REF] Vaskonen | Lower bound on the primordial black hole merger rate[END_REF] but here we provide additional details which are useful for understanding why the suppression of the merger rate is negligible for the ranges of f PBH allowed by present bounds.

A.1 PBH cluster formation

Even if PBHs are uniformly distributed in the early Universe (i.e. according to a Poisson distribution), their discrete nature means that the mean PBH density in a given region can be subject to large fluctuations. Regions with overdensities of PBHs begin to collapse early, forming PBH clusters [START_REF] Inman | Early structure formation in primordial black hole cosmologies[END_REF][START_REF] Chisholm | Clustering of primordial black holes: basic results[END_REF][START_REF] Chisholm | Clustering of Primordial Black Holes. II. Evolution of Bound Systems[END_REF]. Here, we briefly review this mechanism, before commenting on the relevance for the survival of PBH binaries.

At matter-radiation equality, t i = t eq , we consider a spherical region of radius r i , containing an overdensity δ i = (ρ iρ eq )/ρ eq . This overdense region will evolve as an isolated matter-dominated Universe with effective curvature [START_REF] Gunn | On the Infall of Matter into Clusters of Galaxies and Some Effects on Their Evolution[END_REF][START_REF] Fillmore | Self-similar gravitational collapse in an expanding universe[END_REF][START_REF] Bertschinger | Self-similar secondary infall and accretion in an Einstein-de Sitter universe[END_REF],

K = 8πGρ eq δ i r 2 i 3 = 2G Mi δ i r i , (A.1)
where Mi = (4π/3)Gρ eq r 3 i would be the initial mass of the region in the unperturbed

The mean density of a cluster will be

ρ ≡ 3M 4πR 3 ≈ 18π 2 ρ c a -3 = 18π 2 ρ c a -3 eq f 3 PBH N 3/2 , (A.9)
where the second equality follows from the theory of spherical collapse [START_REF] Peebles | The large-scale structure of the universe[END_REF] and from setting a = 1/z coll (the scale factor at which the cluster forms). Here, ρ c is the critical density. Next, we set the PBH velocity v in the cluster equal to some typical velocity dispersion v ≈ σ v ≈ GM/R. Using the definition of ρ in Eq. (A.9), we then have

σ 3 v = GM R 3/2 = 4π 3 G 3/2 M √ ρ . (A.10)
Assuming that there is no substantial growth or evaporation of the cluster after formation, the relaxation timescale can then be written as Here, we have defined ρ m,eq = ρ c a -3 eq (the matter density at matter-radiation equality). We have also taken the numerical values ρ c = 2.78 × 10 11 h 2 M Mpc -3 [START_REF] Zyla | Review of Particle Physics[END_REF] and H 0 = h × 100 km s -1 Mpc -1 , with h = 0.673 [START_REF] Aghanim | Planck 2018 results. VI. Cosmological parameters[END_REF]. This value for the relaxation time t r matches the result given in Eq. ( 7) of [START_REF] Vaskonen | Lower bound on the primordial black hole merger rate[END_REF].

A.3 Cluster collapse

In order to undergo core collapse before redshift z, a PBH cluster must be formed with fewer than a critical number of members N c (z). This critical number N c (z) is obtained by equating the core-collapse timescale to the time between formation and collapse,

18 t r (N c ) = t(z) -t(z c ) , (A.13)
where z c is the redshift at which clusters of size N c are formed. Here, t r (N c ) is the relaxation time for clusters of size N c , defined in Eq. (A.12). Recall also from the previous sections that

z c ≈ z eq f PBH √ N c . (A.14)
With these definitions, we can calculate N c (z) from Eq. (A.13). This is illustrated in Fig. 11. There, we see that for f PBH = 1, the critical size of clusters is N c ∼ 5000 at z = 0. This decreases with decreasing f PBH . Below f PBH ∼ 0.01, core-collapse becomes more or less irrelevant, as the critical number of PBHs tends to one.

A.4 Merger rate suppression factor

Following Ref. [START_REF] Vaskonen | Lower bound on the primordial black hole merger rate[END_REF], the suppression factor is expressed as the probability of a binary not belonging to a cluster that undergoes gravo-thermal collapse. Estimating this factor requires knowledge of the PBH halo mass function, which we describe using the analytical model discussed in [START_REF] Inman | Early structure formation in primordial black hole cosmologies[END_REF]. We indicate the probability of a PBH belonging to a cluster of N elements at given z and f PBH with P (N |z, f PBH ). We thus have

P (N |z, f PBH ) ∝ exp[-N/N * (z, f PBH )],
where N * (z, f PBH ) is the characteristic number of PBHs in clusters forming at redshift z.

The probability of a binary being disrupted by redshift z is given, according to [START_REF] Vaskonen | Lower bound on the primordial black hole merger rate[END_REF], by the sum of two terms:

1. The probability of the binary belonging to a cluster which has reached instability before redshift z, i.e. a cluster with N ≤ N c (z), P where the probabilities of binaries belonging to a cluster P bin (N ) are taken to be ∝ P (N ) but normalised to exclude P (N = 1):

∞ N =2 P bin (N ) = 1. (A.16)
2. The probability of a binary belonging to a sub-cluster of N ≤ N c within a larger cluster (N > N c ),

P (2) p = N >Nc Nc N =3 P sub (N |z c ) P bin (N |z c ) , (A.17)
where the probabilities of belonging to a sub-cluster P sub (N ) are taken to be proportional to P (N ), but normalised summing up to the size of the cluster that contains them

N N =2 P sub (N ) = 1. (A.18)
All the probabilities in the expressions above are evaluated at z = z c , which is the typical redshift of formation of the clusters with N c elements. Then z c is obtained by requiring that N * (z c ) = N c . Finally, the probability of a binary not being perturbed is given by P np = 1 -P

p -P

p . This suppression factor is shown in Fig. 12 for different values of z and f PBH .

B Obtaining the signal-to-noise ratio of an event

In this appendix, we provide a detailed description of how our mock catalogues are constructed. We start by introducing GW signals in Appendix B.1, move to a discussion of the luminosity distance from GWs in Appendix B.2 and finally the computation of the uncertainty on the luminosity distance in Appendix B.3.

B.1 Gravitational waves

The merging of massive compact objects results in the emission of GWs. These can be detected by the strain they produce in GW interferometers, housed in observatories like LIGO and Virgo. Let us note the different terms we use in the following discussion: the term GW detector refers to the equipment used to measure the strain itself; the term GW interferometer refers to the combination of two separate vacuum tubes or "arms" housing the laser and mirrors with a detector at the junction of the arms; and the term GW observatory refers to an entire facility, which may contain multiple detectors, or as in the case of LIGO, multiple interferometers at different physical locations (LIGO Livingston and LIGO Hanford).

The GW strain is given in the transverse traceless gauge as

h(t) = F + (θ, φ, ψ) h + (t) + F × (θ, φ, ψ) h × (t), (B.1)
where h +,× are the plus and cross polarisations of the metric perturbation h αβ , F +,× are the corresponding antenna pattern functions which describe the angular dependence of the sensitivity of the interferometers. The angles (θ, φ) are the spherical polar coordinates of the wave source in the celestial sphere, θ = π/2 being the detector plane; ψ denotes the angle between the two-dimensional basis of the plane orthogonal to the line of sight with respect to which the plus and cross polarisations are defined, and the natural basis (e θ , e φ ) associated with the celestial coordinates (θ, φ) [START_REF] Maggiore | Gravitational Waves: Volume 1: Theory and Experiments[END_REF]. More simply put, ψ is the polarisation angle.

For the results presented in this work, we generate mock catalogues of GW events based on the configuration and sensitivity of the ET, a proposed third-generation GW observatory. Specifically, we consider the ET-D configuration, which involves three interferometers arranged in an equilateral triangle shape and two detectors at each vertex of the triangle, one sensitive to higher frequencies and one to lower, i.e. six detectors in total. The antenna patterns for a single interferometer with a 60 degree opening angle are given by

F + = √ 3 2 1 2 (1 + cos 2 θ) cos(2φ) cos(2ψ) -cos θ sin(2φ) sin(2ψ) , (B.2) F × = √ 3 2 1 2 (1 + cos 2 θ) cos(2φ) sin(2ψ) + cos θ sin(2φ) cos(2ψ) . (B.3)
The antenna patterns for the two other interferometers in the ET-D configuration are then given by F +,× (θ, φ + 2π/3, ψ) and F +,× (θ, φ + 4π/3, ψ).

We assume that the basic information available from the detection of GWs from merging compact objects by the ET will be the luminosity distance to the merger event, D, and the uncertainty on that luminosity distance, σ, which has contributions from the instrumental noise and from the uncertainty due to weak lensing of the GW.

B.2 Luminosity distance

To obtain the luminosity distance for a given simulated event, we firstly draw a redshift from a probability distribution defined based on the merger rate of the progenitor system being considered (binary PBHs or binary PBHs), which we previously described in Section 2.

However, unlike binary neutron star mergers which produce an electromagnetic counterpart, the redshift of an individual11 binary BH merger is always unknown -hence the name dark siren. We therefore convert the simulated redshift to a luminosity distance in the standard ΛCDM cosmology using CAMB. All of the analysis presented in the main text is based solely on the luminosity distance information, rather than the redshift information, in order to accurately simulate the actual analysis that could be done with a catalogue which contains only binary BHs.

Finally, we rescale the luminosity distance to account for weak lensing of the signal. The full computation is shown in Appendix C.

B.3 Uncertainty on the lensed luminosity distance

The second step in the generation of the mock catalogue is then to compute the uncertainty on the lensed luminosity distance for a given event. As we mentioned, this has two components: the instrumental noise and the noise due to weak lensing of the GW. The instrumental noise depends on the signal-to-noise ratio for a given event.

The unlensed signal-to-noise ratio for a single detector is given by

ρ det = ∞ -∞ df h(f ) K * (f ) ∞ -∞ df 1 2 S n (f ) |K(f )| 2 , (B.4)
where h(f ) is the waveform of the GW (the strain as a function of the frequency f ), K(f ) is the filter used in the matched filtering process which maximises the signal-to-noise ratio and S n is the noise power spectral density [START_REF] Li | Extracting Physics from Gravitational Waves: Testing the Strong-field Dynamics of General Relativity and Inferring the Large-scale Structure of the Universe[END_REF]. Asterisks denote the complex conjugate. The optimal unlensed signal-to-noise ratio is obtained using a Wiener filter, and is given by

ρ opt = 4 fupper f lower df h(f )h * (f ) S n (f ) 1 2
, (B.5)

where f upper and f lower are the cutoff frequencies for the strain, beyond which it is assumed to be zero. The waveform h(f ) is the Fourier transform of the strain h(t) given by Eq. (B.1). We randomly draw the angles φ and ψ from a uniform distribution between 0 and 2π. However, the angles θ and ι (the inclination of the event i.e. the angle between the source plane and the detector plane) should have their cosine uniformly distributed, meaning that for these quantities we draw their values randomly from a uniform distribution between -1 and 1 and then take the arccosine of the result. The inclination of the event enters into the function used to weight the final signal-to-noise ratio sum, as we will see in a moment. We use the publicly available package PyCBC12 to generate a mock waveform h(f ) for a given event, using the IMRPhenomD waveform model, and by inputting the masses and simulated luminosity distance of the event. We keep the spins fixed to zero. We also generate the frequency range, and hence f upper and f lower , using PyCBC.

Assuming that the noise can be described by a Gaussian stochastic process, the power spectral density is given by is the autocorrelation of the noise n(t) [START_REF] Li | Extracting Physics from Gravitational Waves: Testing the Strong-field Dynamics of General Relativity and Inferring the Large-scale Structure of the Universe[END_REF]. In practice, we use the publicly available power spectral density data for ET-D 13 . With the waveform and power spectral density in hand, the signal-to-noise ratio for an event seen in a single interferometer can be computed. The total signal-to-noise ratio for an event i seen by the three-armed ET-D observatory is obtained by summing the squares of the individual interferometer signal-to-noise ratios, ρi = w ρ 2 i opt, 1 st arm + ρ 2 i opt, 2 nd arm + ρ 2 i opt, 3 rd arm , (B.8)

S n (f ) = 2
where the weighting factor w is given by w = F 2 + (θ, φ, ψ)(1 + cos 2 ι) 2 + 4F 2 × (θ, φ, ψ) cos 2 ι, (B.9) Events that are above the signal-to-noise ratio threshold (black dashed line) are considered to be properly observed GW events (red for ABHs and yellow for PBHs), while the black and pink points are respectively the ABHs and PBHs with signal-to-noise ratios that make them too faint to be observed. in the linear regime, harmonic transverse-traceless gauge, and assuming that the GW's wavelength is much shorter than the typical curvature radii of the background spacetime geometry, whose covariant derivative is denoted with ∇ ρ in Eq. (C.1).

Introducing the wave ansatz h µν = H µν exp(iw), where H µν and w respectively denote the amplitude and the phase of the GW, Eq. (C.1) implies (i) that the wave follows null geodesics of the background spacetime; and (ii) that its amplitude satisfies

(1 + z) -1 DH µν = const., (C.2)
where z is the source's redshift and D is the electromagnetic luminosity distance -see, e.g., [START_REF] Dalang | Horndeski gravity and standard sirens[END_REF] for further details. It follows that H µν ∝ D -1 , that is, the energy of a GW "dilutes" just like the energy of an electromagnetic wave as it propagates. Since this is valid in any background spacetime, we conclude that GWs experience the same gravitational lensing effects as light.

In the following, we only consider weak-lensing effects, in the sense that we neglect the possibility that a GW source might be multiply imaged. We also assume that the luminosity distance can be computed from the propagation of an infinitesimal beam of null geodesics. In that framework, the distortions of the beam with respect to the homogeneous-isotropic FLRW case are customarily encoded in the so-called distortion matrix. If θ is the observed incoming direction of GWs and β the direction in which they would be observed in FLRW, then the distortion matrix is defined as the Jacobian matrix of θ → β(θ),

A ≡ dβ dθ =   1 -κ -γ 1 -γ 2 + ω -γ 2 -ω 1 -κ + γ 1   , (C.3)
where κ is called the convergence of the beam, γ = γ 1 + iγ 2 its shear distortion, and ω = O(γ 2 ) its rotation. The correction to the observed luminosity distance, D, relative to the FLRW case, D, is quantified by the magnification µ, which is related to the distortion matrix as follows,

µ -1 ≡ D D 2 = d 2 β d 2 θ = det A = (1 -κ) 2 -|γ| 2 + ω 2 . (C.4)
In particular, at lowest order in κ, γ, ω, the magnification only depends on the convergence, µ ≈ 1+2κ. We shall use this approximation for the estimate of the variance of the magnification -see next subsection. It follows from Eq. (C.4) that the absolute uncertainty on the luminosity distance due to lensing may be estimated as

σ lens = 1 2 D σ µ , (C.5)
where σ µ denotes the standard deviation of the magnification.

C.2 Variance of the magnification in weak lensing

We aim to model the PDF of the magnification of the observed GW signals in a realistic inhomogeneous universe. As a first step, we shall estimate the variance of that distribution. Previous analyses, both in the context of supernova cosmology [START_REF] Betoule | Improved cosmological constraints from a joint analysis of the SDSS-II and SNLS supernova samples[END_REF] or GW cosmology [START_REF] Sathyaprakash | Cosmography with the Einstein Telescope[END_REF][START_REF] Zhao | Determination of Dark Energy by the Einstein Telescope: Comparing with CMB, BAO and SNIa Observations[END_REF][START_REF] Cai | Estimating cosmological parameters by the simulated data of gravitational waves from the einstein telescope[END_REF][START_REF] Du | Future constraints on dynamical dark-energy using gravitational-wave standard sirens[END_REF][START_REF] Jin | Forecast for cosmological parameter estimation with gravitational-wave standard siren observation from the Cosmic Explorer[END_REF][START_REF] Hogg | Constraints on the distance duality relation with standard sirens[END_REF], considered a dispersion of magnitude due to that grows linearly with redshift [START_REF] Jönsson | Constraining dark matter halo properties using lensed Supernova Legacy Survey supernovae[END_REF], ∆m = 0.055 z, which corresponds to a magnification dispersion of σ µ = 4 ln(10)∆m/5 = 0.10 z.

While this linear approximation may be valid at low redshift, we expect it to fail at the high redshifts (z ∼ 10 to 100) considered here. The intuition is that at high redshift the Universe is increasingly homogeneous, thereby reducing the growth of the lensing dispersion.

In order to get a more accurate estimate of σ µ across a wide range of redshifts, we shall use perturbation theory at second order. In that framework, we first notice that

σ 2 µ ≡ µ 2 -µ 2 = 4σ 2 κ + O(κ 4 ) , (C.6)
where σ 2 κ is the variance of the convergence, because κ = 0 at linear order and ω ∼ κ 2 ∼ γ 2 . Hence, we may estimate σ µ at second order from the linear-order results on σ κ .

In the flat-sky and Limber approximations, the variance of the convergence for a source at redshift z is related to the convergence power spectrum P κ ( , z) as where η 0 denotes today's conformal time, and r = r(z) is the comoving distance to the source. We compute P W and the associated quantities using CAMB. For a ΛCDM cosmology with H 0 = 67.4 km s -1 Mpc -1 , Ω m = 0.315, Ω b = 0.05, A s = 2 × 10 -9 , n s = 0.965, we find that the standard deviation of the magnification as a function of redshift is very well fit (with sub-percent accuracy) by

σ κ (z) = a arctan (1 + b z c ) d -1 (C.9)
with a = 0.116, b = 1.26, c = 1.46, d = 0.268. 14 The result is depicted in Fig. 14, where we also indicate the linear ansatz of [START_REF] Jönsson | Constraining dark matter halo properties using lensed Supernova Legacy Survey supernovae[END_REF] for comparison. The latter is a reasonably good approximation of σ κ (z) up to z ≈ 3, but it highly overestimates it at high z. At z = 100, we find σ κ ≈ 16%, which is more than an order of magnitude below the prediction of the linear ansatz. This emphasises the importance of a careful modelling of lensing at high redshift.

Figure 1 .

 1 Figure1. The strain amplitude spectral densities (ASD) of advanced LIGO (black) and the ET-D configuration of the Einstein Telescope (green), along with the characteristic strain of two PBH merger events, with component masses of 10 M (orange) and 30 M (purple) taking place at z = 1, that is, with a luminosity distance D ≈ 6.8 Gpc.

Figure 2 .

 2 Figure 2. The selection functions f det (z) for aLIGO (black) and ET-D (green) for three different total masses M : 10 M (dot-dashed lines), 20 M (dashed lines) and 30 M (solid lines). The selection function is defined as the fraction of mergers which are detectable at a given redshift.

Figure 4 .

 4 Figure 4. Average merger rate density of PBH binaries formed in the early Universe, see Eq. (2.10), as a function of f PBH , for M PBH = 10 M , M PBH = 30 M and M PBH = 100 M . The dashed lines are obtained considering the disruption in clusters as described in Sec. 2.2. The coloured band indicatesfor comparison the latest measure of the total binary BH merger rate as reported in[START_REF] Abbott | The population of merging compact binaries inferred using gravitational waves through GWTC[END_REF].

Figure 6 .

 6 Figure 6. Data set generated using darksirens ɭ with the baseline parameters of Table1and Table2, with f PBH = 10 -5 . The blue dots show the observed ABH events, while the yellow dots represent the observed PBH events. At a fixed distance D, the typical uncertainty is smaller for PBHs than for ABHs due to our choice of a larger PBH mass (M PBH = 10 M , M ABH = 7 M ). The dotted horizontal line shows the level of instrumental uncertainty corresponding to our SNR cut of ρ i = 8. The events that lie above this line still satisfy the SNR condition we impose, but receive an extra contribution to the error from lensing.

Figure 7 .

 7 Figure 7. Value of N > as a function of the redshift threshold z * , for a data set with a negligible number of PBH. The error bars represent the uncertainty σ > on N > .

5 TFigure 8 .

 58 Figure 8. Top left: Lowest 3σ-detectable f PBH as a function of the redshift z * discriminating between the low redshift and high redshift subsets, for M PBH = 10 M (red lines) and M PBH = 30 M (yellow lines). Top right: Lowest 3σ-detectable f PBH as a function of z * , for M PBH = 10 M for the fiducial ABH merger rate (red line) and for the one obtained using the GRB SFR (pink line). Bottom: Lowest 3σ-detectable f PBH as a function of the chosen z * , for M PBH = 10 M for an observation time T obs = 1 yr (red line) and for T obs = 10 yr (blue line).

Figure 9 .

 9 Figure 9. Left: Number N > of events in the high redshift subset, z > z * = 10, for each value of f PBH in the baseline case (see Table1 and Table 2) (red points). The grey band shows the uncertainty for N > in the no-PBH case. The blue line shows instead the trend of the detection significance for a departure from the no-PBH case as a function of f PBH . Right: Counts in the small-and large-distance bins are shown by the red bar plot for the baseline case, with z * = 10 and f PBH = 2 × 10 -5 , together with the errors we determine on the value of the counts (left y-axis). In blue, the distances and their errors for all the events contained in the data set are shown as a function of the true redshift z true (right y-axis), while the top x-axis shows the redshift z obs that could potentially be inferred from the distance measurements.

Figure 10 .

 10 Figure 10. Left: recovered mean value f meas PBH (solid line) and 68% and 99.7% confidence level limits (red bands) as a function of the fiducial f fid PBH (black dots). Right: posterior distributions obtained on f meas PBH for low (green), intermediate (blue) and high (navy) values of f fid PBH , with the vertical dotted line showing the value of f fid PBH corresponding to each color. Bottom: Comparison between the 99.7% (3σ) confidence regions for the baseline case (outer band) and for an observation time of T obs = 10 yrs (inner band).

Figure 11 .

 11 Figure 11. Left: Time-scales associated with PBH clusters. Solid lines show the time between the redshift of cluster formation z c and a given redshift z, for clusters containing different numbers of PBHs N PBH . Horizontal dashed lines show the core-collapse timescale for clusters with N PBH members.Where the solid and dashed lines (of a given colour) cross then N PBH = N c (z). Right: Critical number of PBHs per cluster N c as a function of redshift z. Clusters with N < N c (z) will be disrupted due to core-collapse before redshift z.

P

  bin (N |z c ) , (A.[START_REF] Adamek | WIMPs and stellar-mass primordial black holes are incompatible[END_REF]

Figure 12 .

 12 Figure 12. Merger rate suppression factor due to interactions in gravo-thermally unstable clusters, computed following[START_REF] Vaskonen | Lower bound on the primordial black hole merger rate[END_REF]. Left: Probability of a binary not being disrupted as a function of f PBH for different redshifts. Notice that, depending on the value of z, there is a threshold value of f PBH below which no binaries are perturbed. Right: Same probability as a function of redshift, for different values of f PBH .

  τ ) e i2πf τ , (B.6)whereR(τ ) = n(t + τ ) n(t) (B.7)

Figure 13 .

 13 Figure[START_REF] Volonteri | The origins of massive black holes[END_REF]. Trend of the signal-to-noise ratio with event redshift for ABHs (crosses) and PBHs (dots). Events that are above the signal-to-noise ratio threshold (black dashed line) are considered to be properly observed GW events (red for ABHs and yellow for PBHs), while the black and pink points are respectively the ABHs and PBHs with signal-to-noise ratios that make them too faint to be observed.

  related to the power spectrum P W of the Weyl potential W = (Φ + Ψ)/2, Φ and Ψ being the Bardeen potentials, via P κ ( , z) = 2 ( + 1)

  Figure 3. Distribution of merger events as a function of redshift for ABHs and PBHs. We show results for our fiducial ABH model (solid blue), the more optimistic ABH model (dashed blue) based on GRB data (Sec. 2.3), and our PBH model (Sec. 2.2) with f PBH = 10 -5 (yellow). Left: Comoving merger rate density per unit volume per unit time R(z). Right: Redshift probability distribution p(z), defined in Eq.

	10 0
	p(z)

Table 1 .

 1 Cosmological and BH related parameters accessible in darksirens ɭ. The values shown here represent our baseline settings used throughout the article. We do not report a value of f PBH as this will be changed depending on the analysis done with the data.

	clustering

Specifications

T obs [yrs] SNR min lensing z min z max

Table 2 .

 2 Parameters

available in darksirens ɭ to specify the characteristics of the survey and the redshift range over which the mock is constructed. The values shown here represent our baseline settings used throughout the article.

  and Red Consolider MultiDark FPA2017-90566-REDC. FS has received financial support through la Caixa Banking Foundation (grant n. LCF/BQ/LI18/11630014) during the early stages of the project. NBH is supported by a postdoctoral position previously funded through two "la Caixa" Foundation fellowships (ID00010434), with codes LCF/BQ/PI19/11690015 and LCF/BQ/PI19/11690018 respectively, and currently funded by the French Commissariat à l'énergie atomique et aux énergies alternatives (CEA). BJK thanks the Spanish Agencia Estatal de Investigación (AEI, MICIU) for the support to the Unidad de Excelencia María de Maeztu Instituto de Física de Cantabria, ref. MDM-2017-0765. DG has received financial support through the Postdoctoral Junior Leader Fellowship Programme from la Caixa Banking Foundation (grant n. LCF/BQ/LI18/11630014) during the early stage of the project. DG was also supported by the Spanish Agencia Estatal de Investigación through the grants PGC2018-095161-B-I00, IFT Centro de Excelencia Severo Ochoa SEV-2016-0597, and Red Consolider MultiDark FPA2017-90566-REDC during the early stages of the project. DG acknowledges funding from the "Department of Excellence" grant awarded by the Italian Ministry of Education, University and Research (MIUR) in October-December 2021. DG also acknowledges support from the INFN grant "LINDARK," and the project "Theoretical Astroparticle Physics (TAsP)" funded by the INFN in October-December 2021. DG acknowledges support from Generalitat Valenciana through the plan GenT program (CIDE-GENT/2021/017) starting from 01/01/2022. In the early stages of this work, MM and PF received the support of a fellowship from "la Caixa" Foundation (ID 100010434). The fellowship codes are LCF/BQ/PI19/11690015 and LCF/BQ/PI19/11690018 for MM and PF respectively.

Other proposed third-generation observatories include LISA[START_REF] Amaro-Seoane | Laser Interferometer Space Antenna[END_REF], DECIGO[START_REF] Kawamura | Current status of space gravitational wave antenna DECIGO and B-DECIGO[END_REF], TianQin and Taiji[START_REF] Gong | Concepts and status of Chinese space gravitational wave detection projects[END_REF], which will all be space-based rather than terrestrial like the ET and CE, allowing for far longer interferometer arms than is possible on Earth. These instruments will hence be sensitive to a very different GW frequency range to that of ET and CE, which means that the latter will have far better prospects for detecting PBHs in the mass range 1-100 M .

https://dcc.ligo.org/LIGO-T1800042/public.

http://www.et-gw.eu/index.php/etsensitivities.

https://darksirens.readthedocs.io

https://camb.info.

For the results presented in Sec. 4 and Sec. 5, we fix the normalisation specifically using GWTC-3, from which the best estimate of merger rate density of BH binaries is R = 22 Gpc -3 yr -1 .

In practice, this integral is largely insensitive to the precise value of zmax, as the selection function drops rapidly above z ∼ 0.1.

A plot including the events which would not survive the SNR cut is shown in Fig.13.

In practice, we estimate N> by generating a number of mock data sets and taking the mean of the N> values obtained for each.

Note however that the optimal PBH mass depends on the precise value of the event's redshift.

There have been proposals to cross-correlate GW events with galaxy catalogues (see e.g.[START_REF] Herrera | Cross-correlation of the astrophysical gravitational-wave background with galaxy clustering[END_REF][START_REF] Mukherjee | Accurate precision cosmology with redshift unknown gravitational wave sources[END_REF]), allowing the redshift of the event to be estimated, but this introduces a good deal of uncertainty, especially with the current small number of events and relatively poor sky localisation.

https://pycbc.org/.

http://www.et-gw.eu/index.php/etsensitivities.

These numerical values depend on the cosmology; our code includes a notebook ɭ capable of generating such a fitting function for other values of the cosmological parameters.

Acknowledgments

We are very grateful to Evan Hall for making the gw-horizon code public and open. We thank Tjonnie Li, Michele Maggiore and Bangalore Sathyaprakash for useful discussions, and Cristina Fernández Suárez for contributions in the early stages of this work.

MM acknowledges funding by the Agenzia Spaziale Italiana (ASI) under agreement n. 2018-23-HH.0. FS was supported by the Spanish Agencia Estatal de Investigación through the grants IFT Centro de Excelencia Severo Ochoa CEX2020-001007-S and PGC2018-095161-B-I00, and, during the early stages of this work, through the grants Severo Ochoa SEV-2016-0597

We demonstrated that a PBH fraction as low as f PBH ≈ 2 × 10 -5 can be detected with a 3σ significance, given a proper choice of the cut position -see Fig. 8.

The second method takes a Bayesian approach to the problem, and fully exploits the information on the redshift dependence of both the signal and background. Within this framework, we provided an assessment of the posterior probability density function for different fiducial values of f PBH . We found that fractions of DM in the form of PBHs as low as O(10 -5 ) can be measured by ET with one year of data taking, with a precision of ∼ 35 % at 68 % confidence level -see Fig. 10. These results demonstrate a well-defined avenue towards a discovery or, in the absence of a discovery, the prospect of setting a significantly improved upper limit on the existence of PBHs in this mass window.

Furthermore, in the process of this work, we have developed the darksirens ɭ code, which allows for the fast generation of realistic mock GW catalogues sourced by ABH and PBH merger events. We have made the code publicly available for community use, with extensive scope for development and application to other investigations.

The final message of this work is clear: the possibility of detecting and measuring a distant population of PBHs is well within the capabilities of future GW observatories such as the ET. The exciting prospect that these exotic objects could constitute a fraction of the dark matter, a substance that despite its fundamental importance in the Universe we still know so little about, makes the pursuit of this goal nothing short of vital. In the rigorous analysis we have here presented, based on the twin pillars of robust statistics and a thorough treatment of the modelling of both signal and noise, we have shown that the ET will be able to measure a fraction of PBHs as low as O(10 -5 ), through observations of their luminosity distances alone. If a high redshift population of PBHs exists in our Universe, the detection of these distant dark sirens is therefore something we no longer need to hope for, but can begin to expect. background. The evolution of the radius of this region r(t) can be solved parametrically as

CRedIT statement

The region collapses then at

or, using the fact that z ∝ t -2/3 in a matter-dominated Universe:

We can relate t eq and r i by assuming that the size of the region corresponds to its comoving size at matter-radiation equality. With this, we find the redshift of collapse as

for δ i 1. The redshift at which the expansion of the region turns around is a little larger, z ta ≈ 0.56z eq δ i .

Consider now a region of the Universe which we expect to contain, on average, N PBHs with total PBH mass M = N M PBH . If a cluster forms from this collection of PBHs, it is likely to form from a typical overdensity in PBHs of δ PBH = √ N /N , due to Poisson fluctuations. If PBHs make up only a fraction f PBH of the total DM density, then the total overdensity will be δ i = f PBH / √ N (assuming that the DM fluctuations are subdominant). Dropping some order one factors, we therefore expect that the PBH cluster to form at a redshift

While this expression is a good estimate of collapse time for large values of N , it fails for smaller values. In our numerical analysis, we rely on the more complete approach discussed in [START_REF] Inman | Early structure formation in primordial black hole cosmologies[END_REF]. However, in the following sections we will make use of Eq. (A.7) to estimate the typical size of clusters which are disrupted at late times.

A.2 Relaxation time

We will assume that a PBH binary is completely disrupted if the PBH cluster it resides in undergoes core collapse due to gravo-thermal instability and that all other PBH binaries are unperturbed. We take the characteristic core-collapse timescale to be t cc ∼ 18 t r [START_REF] Quinlan | The Time scale for core collapse in spherical star clusters[END_REF], where t r is the relaxation time.

The relaxation time t r for a PBH cluster containing N PBHs. For a system with density ρ and component masses m, this relaxation time can be estimated as [START_REF] Quinlan | The Time scale for core collapse in spherical star clusters[END_REF] t r = 0.065

where ln Λ ≈ ln(N/f PBH ) is the Coulomb logarithm associated with interactions between component masses. Assuming that the fraction of PBHs in clusters matches the mass fraction in the Universe, then the total mass of each cluster is M = mN/f PBH .

where F + , F × are the antenna patterns and ι is the inclination. Lastly, following [START_REF] Li | Extracting Physics from Gravitational Waves: Testing the Strong-field Dynamics of General Relativity and Inferring the Large-scale Structure of the Universe[END_REF], we approximate the instrumental uncertainty on the lensed luminosity distance Di for a given event i as

Here, ρi is the unlensed signal-to-noise ratio for the event determined by the observatory in question and µ i is the magnification of the signal-to-noise ratio due to weak lensing. The factor 2 in Eq. (B.10) accounts for the contribution of the uncertainty on the inclination ι to the instrumental noise [START_REF] Li | Extracting Physics from Gravitational Waves: Testing the Strong-field Dynamics of General Relativity and Inferring the Large-scale Structure of the Universe[END_REF]. Once the rescaling of the unlensed signal-to-noise ratio by the lensing magnification has been carried out, it is then trivial to compute the instrumental uncertainty for the event from Eq. (B.10). However, we note that a precise modelling of the uncertainty on the distance measurement would require access to the full parameter estimation pipeline of the experiment. Furthermore, the uncertainty actually depends on the full network of observatories that may be operating in collaboration with the ET, see [START_REF] Vitale | Parameter estimation for binary black holes with networks of third generation gravitational-wave detectors[END_REF]. Extending the uncertainty computation in this way is beyond the scope of this work.

Furthermore, our analysis also rests on the assumption that the expression for the instrumental uncertainty given in Eq. (B.10) is correct. The signal-to-noise ratio in this expression depends on the noise being stationary and Gaussian. In a real detector this is not exactly the case. Depending on the specific template used in the matched filtering process, so-called glitches -artifacts with very high signal-to-noise ratios -can be seen in the detector [START_REF] Nuttall | Characterizing transient noise in the LIGO detectors[END_REF]. The identification and removal of glitches and other noise artifacts through proper characterisation of the detectors is an important part of current gravitational wave analysis pipelines [START_REF] Davis | LIGO detector characterization in the second and third observing runs[END_REF]. However, since we are working solely with mock data, we generate our catalogues using the assumption of stationary, Gaussian noise. This is effectively equivalent to using a real catalogue of events with glitches removed.

Finally, as mentioned in Sec. 3, after computing the lensed signal-to-noise ratio for all the events, we remove all the events with a lensed signal-to-noise ratio of less than eight from the catalogue. This serves to exclude events which may or may not be true detections of GWs. The signal-to-noise ratio threshold that we use follows that of the LIGO collaboration [START_REF] Abbott | Sensitivity of the Advanced LIGO detectors at the beginning of gravitational wave astronomy[END_REF]. In Fig. 13 we show the computed signal-to-noise ratio as a function of redshift for ABH and PBH events, highlighting the significant amount of events which are discarded due to the cut in the signal-to-noise ratio.

C Lensing

This appendix provides details about our modelling of the weak lensing of GW signals. Appendix C.1 is a short theoretical reminder where we state our assumptions and define the relevant quantities to be used. Appendix C.2 indicates how the variance of the magnification is estimated. Appendix C.3 deals with our model for the magnification PDF in the mocks.

C.1 Magnification of gravitational waves

Consider a GW modelled as a small perturbation h µν of an arbitrary background geometry g µν . If the wave is propagating in vacuum, in the sense that we can neglect its direct interaction with matter, then its equation of motion reads [START_REF] Straumann | Graduate Texts in Physics[END_REF] ∇ ρ ∇ ρ h µν = 0 , (C.1)

C.3 Lognormal probability distribution of the magnification

We now turn to the full PDF p(µ) of the magnification. This distribution must satisfy three physical requirements:

1. It must vanish for µ < µ min , where µ min corresponds to the (de)magnification of Zel'dovich's empty beam [START_REF] Zel | Observations in a Universe Homogeneous in the Mean[END_REF]. This is due to the fact that a bundle of null geodesics cannot be less focused than a beam propagating through pure vacuum. In that case, the angular diameter distance coincides with the affine parameter λ along the bundle; the minimum magnification thus reads

The evolution of µ min with redshift z is depicted in the left panel of Fig. 15.

2. The magnification averaged over sources must be unity [START_REF] Wucknitz | From planes to spheres: about gravitational lens magnifications[END_REF]

This magnification theorem assumes that sources are homogeneously distributed in space. Note that the underlying averaging procedure -over sources rather than over random directions in the sky -is essential to this result. See [START_REF] Breton | Theoretical and numerical perspectives on cosmic distance averages[END_REF] and references therein for detailed discussions.

3. The variance of the magnification should coincide with the one evaluated in Appendix C.2,

There are, of course, many models that would satisfy the above three requirements. We adopt, for simplicity, the following shifted lognormal model, The resulting magnification PDF is depicted in the right panel of Fig. 15 for various values of the source redshift. Despite its simplicity, the shifted lognormal model mimics important properties of the expected magnification distribution in the inhomogeneous Universe. On the one hand, p(µ) peaks at µ < 1, thereby encoding the fact that most lines of sight have a lower column density than average, because voids occupy more volume in the Universe. On the other hand, p(µ) exhibits a longer tail towards high magnifications, which encodes that, albeit rare, overdensities can produce large magnifications. Note however that our model tends to underestimate the probability of those high magnifications compared to what is obtained with ray tracing in N -body simulations.