Gravitational-wave event rates as a new probe for dark matter microphysics

We show that gravitational waves have the potential to unravel the microphysical properties of dark matter due to the dependence of the binary black hole merger rate on cosmic structure formation, which is itself highly dependent on the dark matter scenario. In particular, we demonstrate that suppression of small-scale structure -- such as that caused by interacting, warm, or fuzzy dark matter -- leads to a significant reduction in the rate of binary black hole mergers at redshifts $z\gtrsim5$. This shows that future gravitational-wave observations will provide a new probe of physics beyond the $\Lambda$CDM model.


I. INTRODUCTION
The standard ΛCDM model of cosmology has been posited to explain a range of observations spanning the largest observable scales to the scale of galaxies.Its success at explaining these data relies on two mysterious components: dark energy in the form of a cosmological constant (Λ) and dark matter (DM).Dark matter, in particular, is key to explaining the formation and evolution of structures such as galaxies and clusters of galaxies in the late-time Universe.
Present observations on supergalactic scales are compatible with the hypothesis that the dark matter is coldi.e., the particles have vanishingly small thermal velocities.In this cold dark matter (CDM) model, the particles also do not have significant nongravitational interactions [1][2][3].These requirements can be satisfied by particlelike dark matter of typical mass keV, or wavelike dark matter of mass 10 −22 eV.However, the key to determining the fundamental nature of dark matter lies in subgalactic scales at large redshifts.This is because the onset of nonlinear structure formation can be very sensitive to the microphysics of the dark matter [4][5][6].In turn, any probe that could shed light on this epoch of structure formation would also provide valuable insights into the nature of dark matter [7,8].
There are three classes of phenomenological particlelike and wavelike dark matter scenarios that generically predict small-scale signatures that differ from the predictions of standard CDM.The first, warm dark matter (WDM) scenario usually assumes negligible interactions but has a small DM particle mass in the low keV range that allows these particles to free-stream out of small-scale perturbations [9][10][11][12][13][14][15][16][17][18][19].The second, interacting dark matter (IDM) scenario makes no strong assumption about the particle mass but endows the DM particle with non-negligible interactions (see Sec. II for details).The third, fuzzy dark matter (FDM) scenario comprises a condensate of ultra-light DM particles of mass ∼ 10 −22 -10 −21 eV whose collective behaviour is wavelike [20][21][22][23], although Lymanα data implies a mass > 10 −20 eV [24].As shown in Fig. 1, although the details differ between WDM, IDM, and FDM, all three scenarios predict a cutoff in the linear matter power spectrum at large wave numbers k, a feature that is often invoked as a possible solution to the claimed "small-scale crisis", a long-standing set of discrepancies where simulations predict more small-scale structures than observed [25].
In this article we present a new observational probe of dark matter microphysics, namely the merger rate of stellar-mass binary black holes (BBHs) throughout cosmic time. 1 With the recent advent of gravitational-wave .Linear matter power spectra for several DM-neutrino interacting scenarios of varying interaction strength uνχ, along with predictions of their "equivalent" FDM and thermal WDM scenarios characterised by their respective particle masses, max and m wdm , tuned to give a similar small-scale suppression.The FDM power spectra have been computed using A x i o n C A M B [34,35]; the WDM power spectra are derived from the ΛCDM power spectrum filtered with a transfer function [10].Vertical lines indicate rough k-scales corresponding to the mass resolutions used in g a l f o r m. astronomy, advanced interferometers like LIGO [26] and Virgo [27] have already detected dozens of BBHs [28,29].Next-generation observatories such as Einstein Telescope [30] and Cosmic Explorer [31] will be capable of detecting nearly all stellar-mass binary black holes in the observable Universe [32], extending the reach of gravitational-wave astronomy out to the beginning of cosmic star formation.The space-based interferometer LISA [33] will also provide a great deal of complementary low-frequency information about this population.The abundances of these systems at different epochs encode the star formation, chemical enrichment, and other baryonic processes occurring within dark matter haloes.As we demonstrate below, the BBH merger rate is highly sensitive to the suppression of small-scale structure induced by dark matter microphysics, as this in turn inhibits star formation and reduces the BBH yield, particularly at redshifts z 5.
In the following, we use as a working example a DMneutrino interacting scenario [4,6,[36][37][38][39][40][41][42][43] to demonstrate quantitatively the potential of the BBH merger rate to probe cosmologies with suppressed small-scale structure.A brief description of this IDM scenario is given in Sec.II.Our analysis pipeline, shown in Fig. 2, is as follows.We compute the linear matter power spectrum of an IDM cosmology using c l a s s [42,44,45], which is then fed into the semianalytic model of structure formation g a lf o r m [46,47] to produce synthetic realisations of galactic Explorer) because the physics underpinning their merger rates is much better understood at present than that of the supermassive BBHs probed by GW searches at lower frequencies (principally by pulsar timing arrays).populations, as described in Sec.III.Concurrently, we perform an N -body simulation of the IDM cosmology using the g a d g e t -4 code [48,49], initialised with the same linear input, in order to cross-check the halo mass functions predicted by g a l f o r m.Finally, the star formation rate and stellar metallicities computed by g a l f o r m are fed into the binary population synthesis code c o mpa s [50][51][52], described in Sec.IV, which computes the gravitational-wave event rate.Our results and conclusions are presented in Sec.V.

II. INTERACTING DARK MATTER
We consider the case in which the dark matter scatters elastically with (massive) neutrinos, via a constant, velocity-independent cross section σ 0 as described in [42].This kind of interaction can be realised in some particle physics theories beyond the Standard Model [5,[95][96][97].From the phenomenological perspective, however, the cosmological dynamics of this interaction can be captured by a single parameter, the interaction strength, normally parameterised in terms of the dimensionless quantity [39,40,42,65] where m χ is the DM mass, σ 0 the interaction cross section, and σ Th ≈ 6.65×10 −29 m 2 is the Thomson scattering cross section.The case of DM-massive neutrino interaction was previously implemented in the Einstein-Boltzmann solver c l a s s [44,45] in [42], which we use in this work to generate the linear matter power spectra for our IDM scenarios.
As noted earlier, the interaction suppresses the matter power spectrum below a certain scale, similar to the predictions of WDM scenarios.As collisions occur between dark matter and neutrinos, the former is prevented from collapsing to form structures, leading to a suppressed linear matter power through collisional damping.This type of damping is physically distinct from free-streaming in WDM scenarios, where the large initial velocities result in the WDM particles escaping from initial overdensities.It is also distinct from the damping seen in FDM scenarios, which arises from quantum pressure hindering gravitational collapse.In general, the stronger the coupling, the later the DM decouples from the neutrinos and hence the larger the scales affected (see Fig. 1).Observe also that IDM scenarios predict in addition prominent acoustic oscillations at higher wave numbers.The origin of these acoustic oscillations is analogous to the baryon acoustic oscillations (BAO) from the coupling of photons and baryons at early times, and the scale at which the oscillations appear is again governed by the dark matter and neutrino decoupling times.
Cosmic microwave background (CMB) anisotropy and BAO measurements currently constrain the DM-neutrino interaction strength to u νχ 10 −4 (95% confidence) [42].Including Lyman-α forest measurements strengthens the bound to u νχ 10 −5 (95% confidence), albeit with an apparent ∼ 3σ preference for a nonzero value centred around u νχ ∼ 5 × 10 −5 [43] (although other studies claim to rule out WDM with similar suppression [98,99]).Given these constraints on u νχ , the interaction only occurs at appreciable rates at redshifts z 1000.At z 1000, the DM particles are effectively collisionless and cold as in standard CDM, and hence amenable to standard N -body simulations of nonlinear structure formation.

III. SIMULATING STRUCTURE FORMATION
Our modelling of nonlinear structure is based on a combination of semianalytic computations and N -body simulations.We use the semianalytic model of structure formation g a l f o r m, first introduced in [46], which provides a flexible and computationally inexpensive way to produce synthetic realisations of galactic populations in a cosmological setting.The latest version [47]   Figure 3. Galaxy luminosity functions in the bJ band for several IDM scenarios of different interaction strengths.Observational data are from the 2dF survey [107].It is evident that while the more weakly interaction IDM scenarios and ΛCDM are roughly consistent with observations, the uνχ ∼ 10 −6 case deviates significantly, and even in the case with feedback turned off, cannot reproduce the low-luminosity tail, despite overpredicting at high luminosities.This corresponds to a slightly weaker interaction than the preferred value of uνχ = 5.5 +2.6 1.1 × 10 −6 from [43].We do not show even more strongly interacting IDM scenarios, as these deviate from observations even further.
We generate merger trees using the Monte Carlo technique described in [102] (which is itself based on the extended Press-Schechter (EPS) theory, and is calibrated to the results of N -body simulations).In models in which the linear power spectrum P (k) has a cutoff, as in our IDM scenario, a small correction is required to the EPS formalism: to obtain the variance of the density field σ(M h ), P (k) needs to be convolved with a sharp k-space filter rather than with the real-space top-hat filter used for CDM [104].Using our Monte Carlo technique rather than N -body simulations to generate merger trees has the advantage that IDM scenarios can be studied at minimum computational expense while avoiding the complication of spurious fragmentation in filaments that occurs in N -body simulations with a resolved cutoff in P (k) (e.g., [105,106]).
We calibrate the astrophysical parameters in g a lf o r m to obtain the best possible fit to observational data for ΛCDM.In particular, we adjust parameters controlling the strength of supernova feedback to achieve the best possible match to the z = 0 galaxy luminosity functions in the b J and K-bands, using observational data from [107,108].We find that for IDM scenarios with u νχ < 10 −6 , the same parameters used for ΛCDM provide an equally good fit to the observational data.For u νχ 10 −6 , however, the feedback strength (as a function of halo mass) has to be reduced significantly to obtain the best possible fit.As shown in Fig. 3, this is still a worse fit than for the more weakly interaction scenarios, particularly for galaxies fainter than 19th mag-  .The suppression of the halo mass function, as described by Eq. 2, for a range of interaction strengths.The dotted lines show the mass resolution of our g a l f o r m merger trees for both the standard and high-resolution settings, while the dashed line shows the atomic hydrogen cooling limit, below which star formation is not expected [109].
nitude.We find that it is impossible for scenarios with u νχ 10 −5 to produce realistic galaxy distributions, since even with all astrophysical feedback mechanisms suppressing galaxy formation turned off, not enough galaxies are produced.The failure to reproduce the observed galaxy luminosity functions at z = 0 therefore already poses a stronger constraint on u νχ than the previous tightest bounds from Lyman-α and rules out the preferred value of u νχ = 5.5 +2.6 1.1 × 10 −6 found in [43].To ensure the validity of our g a l f o r m results, we compute the merger trees with two different mass resolution settings, which determine the smallest tracked progenitor halos.As shown in Fig. 4, our "standard" resolution is well below the minimum halo mass expected to form in IDM models with interaction strength u νχ > 10 −6 .For the models with weaker interactions (u νχ ≤ 10 −7 ), the resolution is not high enough to capture the formation of the smallest halos, though it does capture a large part of the suppression in the u νχ ∼ 10 −7 case.This motivates our choice of the "high" resolution halo merger trees.Conveniently, this corresponds nearly exactly with the atomic hydrogen cooling limit, which marks the lowest mass haloes expected to form stars [109].This means that, although even our high-resolution halo merger trees don't resolve the smallest haloes expected to form in the models with u νχ < 10 −7 , we are still able to capture the majority of the star formation expected to occur in these models.

A. N -body simulations and the halo mass function
To ensure that g a l f o r m accurately predicts the halo mass function of the IDM scenarios, we cross-check it by performing N -body simulations using a version of and the Sheth-Tormen formalism (dashed lines).The choice of uνχ = 3.4 × 10 −4 has been motivated by the CMB+BAO constraint found in Ref. [42].It is evident that the standard Sheth-Tormen halo mass function does not provide a good description of our simulation results.
g a d g e t -4 [48] modified to include the effects of massive neutrinos [49] and initialised with the matter power spectra of the IDM scenarios of Fig. 1.We find that the g a l f o r m and N -body results are compatible with each other to the ∼ O(10%) level in the mass range over which we generate merger trees using the Monte Carlo technique.
Our simulations are performed using a modified version of g a d g e t -4 [48], which includes massive neutrinos via the S u p e r E a s y linear response method of [49].g a d g e t -4 is a massively parallel N -body code supporting both collisionless (gravity-only) simulations and hydrodynamical simulations for baryonic matter and star formation.We restrict ourselves to collisionless simulations, which suffice for the purpose of cross-checking the g a l f o r m halo mass function output.Initial conditions are generated at z = 49 using the built-in n g e n i c implementation, with linear matter power spectra from our modified version of c l a s s as input.For interaction strengths allowed by current cosmological data, dark matter-neutrino decoupling happens at redshifts z 1000, well before the simulation initialisation time.Thus, the interaction affects nonlinear structure formation only through initial conditions, and there is no need to modify the main body of g a d g e t -4 itself.
To probe the effects of our IDM scenario at a range of scales, we run our simulations using N = 512 3 particles in several different box sizes, (50 Mpc) 3 , (100Mpc) 3 , (500Mpc) 3 , and (1Gpc) 3 , for different values of the interaction strength u νχ , including a noninteracting (u νχ = 0) set.All other cosmological parameters are kept fixed to the CMB+BAO best-fit values found in [42].We also perform a set of simulations using the best-fit parameters found in [43], which is phenomenologically similar to our  2); dotted and dot-dash lines are two fitting functions from Refs.[110,111].The bottom panel of each plot shows the fractional residuals of the three fitting functions.
We use the friends-of-friends algorithm to identify haloes in our simulations.For halo masses we use the convention M 200c , i.e., the mass contained within a sphere with mean density 200 times the critical density of the Universe.The haloes from each simulation are sorted into 30 evenlyspaced logarithmic mass bins to obtain the halo number density per mass bin, i.e., the halo mass function.As is well-described in the literature (see, e.g., [105,106,[112][113][114]), the halo mass function of cosmologies with a smallscale cutoff exhibits unphysical discreteness effects on small scales (relative to the box size).We therefore conservatively choose to disregard the low-mass end of each halo mass function thus constructed, cutting it off just above the scale at which discreteness effects become evident.Then, to produce a halo mass function that spans a range of mass scales, we simply stitch together the pieces obtained from our different box-size runs.
We have attempted to describe the IDM halo mass functions semianalytically using the Sheth-Tormen mass function, with the IDM linear matter power spectra as input and a real-space spherical tophat filter for smoothing.However, as is evident in Fig. 5, the Sheth-Tormen mass function thus constructed fails to capture the low-mass suppression relative to the CDM case seen in our Nbody simulation results.We therefore forego the Sheth-Tormen approach completely, opting instead to use a fitting function to describe the suppression relative to the CDM case.
We find the following fitting function to serve the purpose well: where M β ≡ (4π)ρ m /3k −3 β is the mass corresponding to the wave number k β at which the linear matter power spectrum of an IDM scenario is suppressed by β% relative to the CDM case, ρm is the comoving mean matter density, and the parameters α and β are to be determined by fitting (2) to our N -body results.The fit is achieved by forming and then minimising a figure-of-merit (FoM) of the form We find the best-fit values to be α = 0.9 and β = 10.The fit is displayed against simulation results in Fig. 6.We also give in Table I the M 10 values corresponding to several interacting strengths u νχ .To confirm the robustness of the fit, we test the fitting function (2) against simulations at several low u νχ values outside of the calibration range, and find the fit to be accurate to better than ∼ 10%.
Finally, we remark that extended Press-Schechter formalism used in g a l f o r m does not exactly use our fitting function, but produces a halo mass function that matches uνχ M10/M k10/h Mpc −1 3.4 × 10 −4 1.0 × 10 14 0.17 Table I.IDM interacting strengths and the corresponding suppression mass scale M β that appears in the fitting function (2) for the halo mass function.We find β = 10 to be the bestfit value against simulation results.k10 is the corresponding k-scale.
it to the same or better precision as the fitting function matches the simulations.

IV. GRAVITATIONAL-WAVE EVENT RATES
Our key observable is the redshift-dependent BBH merger rate density R BBH (z).In principle we could also use our simulations to study the merger rates of binary neutron stars and black hole-neutron star binaries; however, these are only detectable at lower redshifts due to their smaller masses and correspondingly weaker gravitationalwave signals.They are thus less sensitive as a probe of the high-redshift suppression effects we are interested in.
Stellar-mass BBH mergers are thought to form primarily through isolated binary evolution, in which massive binary stars end their lives as black holes and radiate orbital energy through gravitational-wave emission until they eventually merge, generating an observable gravitationalwave event [115].The BBH merger rate is thus essentially a delayed tracer of star formation, whose normalisation depends on the efficiency with which massive binary stars are converted into BBHs.This efficiency is mostly determined by the stellar metallicity (i.e., the fraction of the stellar mass that is in elements heavier than helium).
We model the merger rate using the binary population synthesis code c o m pa s [50][51][52], which generates a synthetic BBH population by sampling individual stellar binaries from initial mass and metallicity distributions and evolving them over cosmic time.c o m pa s tracks the stellar and orbital evolution of each binary, accounting for a range of physical processes such as stellar winds, mass transfer, BH formation at the end of each star's lifetime, and subsequent shrinking of the orbit via gravitationalwave radiation reaction.We use a c o m pa s dataset of 20 million evolved binaries (resulting in ≈ 0.7 million BBHs) presented in [116], which is publicly available at [117].This gives us the BBH formation efficiency as a function of initial mass and metallicity, as well as the delay time between star formation and BBH merger.By combining this with a model for the star formation rate density and metallicity distribution as functions of redshift, we can use the c o m pa s "cosmic integration" module [118] to average over the synthetic population and obtain the cosmic BBH merger rate.These inputs are typically phenomenological models chosen to fit low-redshift observational data.In our case, we instead use the star formation rates and metallicity distributions generated by g a l f o r m, allowing us to model how the BBH rate changes as we vary the underlying dark matter model.
Our key results for the BBH merger rate are shown in Fig. 7.In addition to the baseline ΛCDM case, we consider a range of IDM models with u νχ = 3.4 × {10 −6 , 10 −7 , 10 −8 } (recall that models with stronger DMν interactions than this are unable to reproduce the galaxy luminosity function data shown in Fig. 3).The factor of 3.4 comes from the CMB+BAO limits in [42]; we have chosen interaction strength values that are whole factors  [119] is shown as a green band.Curves marked "hr" are high resolution, meaning that the g a l f o r m merger trees are tracked to smaller progenitor masses.The curve marked "G.C." is the contribution from globular clusters in the noninteracting scenario. of 10 smaller.The rates that we obtain depend on the mass resolution of our g a l f o r m merger trees, with finer resolution allowing us to access additional star formation taking place in smaller haloes, resulting in a greater total number of BBHs.Most of the results shown in Fig. 7 use a benchmark mass resolution of M h ≥ 10 9.5 M , but we additionally show results for CDM and for u νχ = 3.4 × 10 −8 at a finer resolution of M h ≥ 10 8.5 M .At this mass resolution, we are able to resolve the majority of the star formation in the Universe (under the assumption that it takes place through the cooling of atomic hydrogen), except for that from the ultrafaint dwarf population.We do not anticipate these galaxies to contribute significantly to the overall BBH merger rate.
One limitation of our approach is that c o m pa s only accounts for BBHs formed through isolated binary evolution, neglecting other formation channels such as dynamical capture in dense stellar environments [115,120] (globular clusters [121], nuclear star clusters [122], etc.) or in AGN disks [123], as well as the possible presence of primordial BBHs [124].However, these additional merger channels can be safely neglected in a first analysis, as they are expected to be subdominant compared to the contribution from isolated binary evolution.For example, Fig. 7 shows the contribution from dynamical BBH assembly in globular clusters, as calculated in [125] by combining cluster N -body simulations [126] with a semianalytical cosmological model of globular cluster abundances [127].We see that this contribution is subdominant at all redshifts, justifying our focus here on isolated binary evolution.In principle, however, one could incorporate subdominant channels such as these into our analysis with additional modelling, e.g., so that we capture the suppression of globular cluster abundances due to IDM.

A. BBH merger rate as a function of host halo mass
The advantage of our proposed method for constraining dark matter models is that it allows us to target the high-redshift, low-mass haloes that are most strongly suppressed in these models.This is highlighted in Fig. 8, where we break the total redshift-dependent BBH merger rate down into bins in halo mass for CDM and for an IDM model with u νχ ∼ 10 −7 .We see that the strongest contribution to the observed merger rate, across all redshifts, comes from light haloes, M h 10 10 M .These are precisely the haloes that are suppressed for DM-ν interactions on the order of u νχ ∼ 10 −7 (cf.Fig. 4), which explains why the BBH merger rate is a useful probe of this suppression.
In particular, we see in Fig. 8 that the suppression is most apparent at high redshifts, z 7.This again highlights the advantages of our method, as galaxies in these light haloes are challenging to observe directly at such high redshift, even with JWST.For instance, at z 7 the observational limit for JWST assuming a Hubble Ultra Deep Field-like survey is around M ∼ 10 7.5 M , or M h ∼ 10 9−10 M [128,129] depending on the stellar-tohalo mass relation assumed.Even under the assumption of this optimistic scenario, the BBH merger rate technique allows us, in principle, to probe halo masses one to two orders of magnitude lower.

B. Future gravitational-wave detection forecasts
In order to convince ourselves that differences in the high-z BBH merger rate will be detectable with future gravitational-wave observatories, we compute forecasts for the expected number of BBHs, N > (z * ), detected above a given threshold redshift z * , as well as the expected uncertainty in this quantity due to Poisson fluctuations in the number of BBHs at each redshift, as well as redshift measurement errors.(This is based on the "cut-and-count" method discussed in [130].)We assume one year of observations with a third-generation interferometer network consisting of Einstein Telescope [131] and two Cosmic Explorers [132], accounting for the detection efficiency and redshift uncertainty associated with this network.As shown in Fig. 9, we can clearly distinguish between the different N > (z * ) predictions, allowing us to confirm or rule out a small-scale suppression of the scale caused by DM-neutrino interactions down to the level of u νχ ∼ 10 −7 .
Figure 9 clearly demonstrates that the main limitation to the cosmological information we can extract from future gravitational-wave detections is not statistical uncertainty, but systematics associated with modelling choices.In particular, since the observational study of BBHs is still at an early stage, there are numerous astrophysical uncertainties associated with the stellar physics that governs isolated binary evolution.We return to this issue later in this section, where we investigate the impact of these uncertainties and show that they are not degenerate with the DM suppression effect.
Our key observable quantity is constructed as follows: we choose some redshift threshold z * , and define N > (z * ) as the number of detected BBHs whose best-fit inferred redshift falls above this threshold, ẑ > z * .Our results show that different dark matter scenarios will predict different values of N > .We can predict whether these models are distinguishable in a given observational scenario by forecasting the corresponding uncertainty on N > .It is straightforward to write down the mean value of this quantity, Here, T obs is the total observing time, V is the comoving volume as a function of redshift, R is the mean BBH  merger rate density as a function of redshift, θ is a vector of BBH parameters (masses, spins, etc.) with population distribution p pop (θ|z) (which we allow to evolve over redshift), P det (z, θ) is the probability of detecting a BBH with parameters θ at redshift z, and P > (z, θ|z * ) is the probability of inferring a best-fit redshift greater than z * for the same BBH, the width of which encapsulates redshift measurement uncertainties for a given GW detector network.We describe how to calculate each of these probabilities below.
In order to calculate the variance of N > , we assume that the number, n, of BBH signals emitted from a comoving volume V in a source-frame time T = T obs /(1 + z) is a Poisson random variable with mean n = V T R. We further assume that Poissonian fluctuations in the BBH rate at different redshifts are uncorrelated with each other.This allows us to write the covariance matrix between these BBH number counts as We thus find that which shows that N > is not itself Poissonian, since its mean and variance are unequal so long as the θ integral is not equal to unity.We assume the BBH signals are detected using a matched-filter search, in which the detection statistic is a noise-weighted inner product of the data d with the signal template h, with S(f ) the noise power spectral density of the detector.In Gaussian noise, the signal-to-noise ratio (SNR) of the search is then approximately a Gaussian random variable with unit variance and mean ¯ = (h, h); so for example, having ¯ = 5 would constitute a 5σ detection in this idealised scenario.In reality, the noise is not Gaussian, and we require a larger SNR in order to confidently discriminate a BBH from various non-Gaussian noise transients.A commonly adopted detection threshold is ¯ ≥ 8.We can therefore write the detection probability as which interpolates between (almost) zero for faint signals and (almost) one for strong signals, with a characteristic width set by the randomness of the observed SNR.
We approximate the redshift inference using a Fisher forecast.The Fisher matrix for parameters θ i (including z) of the signal model h(θ i ) is given by In the strong-signal limit, the inferred redshift ẑ is a Gaussian random variable with mean equal to the true redshift z, and standard deviation given by the inverse Fisher matrix, We therefore have which interpolates between (almost) zero for nearby BBHs and (almost) one for distant BBHs, with a width set by σ z .In order to evaluate the expected SNR ¯ and the redshift uncertainty σ z , we need a signal model h(θ, z).For this, we use the phenomenological hybrid waveform model of [133].For simplicity, we assume that the BBHs are all equal-mass and nonspinning (as is approximately true for most of the BBHs detected by LIGO/Virgo thus far [119]), and we average out all of the extrinsic parameters (sky location, polarisation angle, inclination, etc.).This reduces the number of relevant parameters to just two: the redshift z and the total mass M .In terms of these parameters, the frequency-domain waveform model can be written as where f merge , f ring , f cut , f width , ϕ 0 , and {ψ k } are all constants given in [133].This expression is appropriate for z → 0; in order to generalise to arbitrary redshift we simply replace r → d L (z) (with d L the luminosity distance), M → (1 + z)M , and f → f /(1 + z) for all of the frequencies.
We use this waveform to compute the expected SNR ¯ and the redshift uncertainty σ z (the latter of which involves computing and inverting the 2 × 2 Fisher matrix for the mass and redshift), as shown in Fig. 10.The final missing ingredient to calculate N > and its variance is the population distribution of the BBH masses, for which we use the best-fit broken-power-law model from [134], which transitions from a shallow m −1.58 power law at low masses (5.9 M < m < 41 M ) to a steeper one, m −5.59 , at large masses (41 M < m < 87 M ).Here m is the source-frame mass of each individual black hole, so M = 2(1 + z)m.

C. Varying astrophysical model parameters
While our results show that the high-z BBH merger rate is sensitive to DM microphysics, it undoubtedly also depends on the astrophysical modelling choices we have made in our pipeline.In order to investigate the extent to which variations in these choices are degenerate with the DM suppression effect, we use the c o m pa s simulation data from [135] (publicly available at [136]), consisting of 19 alternative models in which the binary stellar evolution modelling is systematically modified, isolating the effects of mass transfer, common envelopes, supernovae, and stellar winds.We focus on modifications to c o m pa s modelling choices rather than those in g a l f o r m, as the latter are currently much more tightly constrained by fitting to observational data.In fact, the predictive power of the g a l f o r m model is enhanced by the fact that the same set of parameters that show consistency with low redshift observational data also predict reasonably accurate galaxy abundances at high redshift (relevant to the epochs of interest in this paper) across multiple wavelengths [137][138][139].
In Fig. 11 we show the BBH formation rate and merger rate in 21 different scenarios: our fiducial CDM model ("A"), 19 alternative CDM models ("B"-"T") with alternative modelling of binary stellar evolution, and a fiducial IDM model with u νχ = 3.4 × 10 −7 .We see that the IDM model is clearly distinguished from the other 20 models in terms of its BBH formation rate, with a characteristic high-z suppression relative to these other models.However, the distinction between IDM and the other alternative models is blurred when considering the BBH merger rate, which we can understand as being due to the way these models modify the distribution of delay times between formation and merger.Nonetheless, a more detailed analysis reveals that IDM is still distinguishable from the alternative models, even if we only have access to the BBH merger rate.To demonstrate this, we define a parameterised merger rate model, where the index i labels models IDM, B, C, . . ., T. The case where θ i = 0 for all i corresponds to our fiducial model A, while having R equal to a particular alternative model corresponds to θ i = 1 for that case and = 0 otherwise.Between these special points in parameter space, (13) linearly interpolates between the 21 different models.We can therefore investigate the degeneracy between the different models by carrying out a Fisher analysis on the θ i parameters.We stress that this is not intended to represent a realistic model of the BBH merger rate over the full space of astrophysical model parameters (this would require us to carry out simulations for a large number of points sampled from the 20+ dimensional parameter space, which is prohibitively expensive), but rather to provide an approximate tool for assessing the degeneracy of IDM suppression with astrophysical modelling choices.
To perform our Fisher analysis for the parameters {θ i }, we imagine sorting a large number of observed BBHs into redshift bins, and counting the number N (z) in each bin. 2  We model these counts with a Gaussian likelihood, where the mean count in each bin, N (z), is calculated from (4) using the parameterised merger rate (13), and 2 These number counts are related to the "cut-and-count" quantity N> from the previous section by N (z) = N>(z min ) − N>(zmax), where zmax, z min are the upper and lower edges of bin z.
the variance σ 2 z is calculated from ( 6) using the fiducial rate R A .(Recall that this variance includes redshift measurement errors.)The Fisher matrix is then In Fig. 12 we show the covariance matrix, Cov(θ i , θ j ) = (F −1 ) ij , computed from (15) using 38 redshift bins of width ∆z = 0.25 over the range z ∈ [0, 9.25], assuming one year of observations with a third-generation GW interferometer network as described in the main text.We are particularly interested in the on-diagonal elements of this matrix, which indicate the measurement variance for each of the θ i ; roughly speaking, a variance significantly smaller than unity for a given model indicates that the BBH merger rate in that model can be distinguished from all other models with one year of data.We find that Var(θ IDM ) ≈ 0.18, so that the DM suppression effect is not degenerate with the other models considered here, at least for an interaction strength of u νχ ∼ 10 −7 ; in fact, IDM is the only model in which this is the case, indicating that the alternative models B-T are all degenerate with each other.This variance corresponds to a ∼ 2.3σ detection of the IDM suppression after one year of observations.Since the Fisher information grows linearly with observation time (see ( 4) and ( 6)), this implies a 5σ detection with ∼ 4.6 yr of data.

V. SUMMARY AND CONCLUSION
We have investigated the suppression of the merger rate of binary black holes in scenarios where dark matter scatters off of neutrinos, taking this as an example of a broader family of dark matter models in which structure is suppressed on small scales.Our predictions are based on a simulation pipeline (cf.Fig. 2) that captures all of the crucial physical ingredients of the problem, from the initial (linear) matter power spectra, through the formation and evolution of galaxies, to a population of merging BBHs.Our key results are shown in Figs.7  and 9, which respectively show the suppression of the BBH merger rate and demonstrate that this suppression will be detectable with next-generation gravitational-wave observatories.
All of the dark matter scenarios we consider are consistent with the local (z = 0) merger rate inferred by LIGO/Virgo, R BBH (z = 0) = [16-61] Gpc −3 yr −1 [119], with only a small suppression ( 20%) of this local rate in the u νχ ≤ 3.4 × 10 −7 IDM scenarios compared to ΛCDM.At higher redshifts, however, we see that their differences become increasingly significant.Physically this reflects that most star formation takes place in smaller haloes at these early epochs, and these smaller haloes are where the effects of DM interactions are most pronounced.Even for interactions as weak as u νχ = 3.4 × 10 −8 there is a small but visible difference from CDM in the predictions from our high-resolution simulations.Decreasing the interaction strength further down to u νχ = 3.4 × 10 −9 , however, makes the suppression scale so small that there are very few star-forming haloes that are missing compared to ΛCDM, and consequently the BBH merger rate for this scenario is essentially indistinguishable from that in ΛCDM.
To confirm that the different dark matter models can indeed be distinguished in this way, we have carried out a Fisher forecast for the redshift-dependent BBH number count N > (z * ), as shown in Fig 9 .We find that future observations with Einstein Telescope and Cosmic Explorer will allow us to confirm or rule out suppression caused by DM-neutrino interactions down to the level of u νχ ∼ 10 −7 .
The main challenge for this method is the significant astrophysical uncertainties in the modelling.We stress however that there is good reason to be optimistic regarding these systematic uncertainties: any changes to the binary evolution model will have an impact at all redshifts and should therefore have minimal degeneracy with the specific high-redshift suppression we are targeting.We envisage that the very large numbers of BBHs detected by third-generation interferometers at low redshifts will be able to pin down the astrophysical modelling, allowing us to use the high-redshift tail as a sensitive probe of structure formation and dark matter.We demonstrate this explicitly by investigating 19 alternative c o m pa s models [135]; we find that the suppression effect from IDM (at the level of u νχ ∼ 10 −7 ) can be confidently distinguished from each of these alternatives with next-generation GW observations.Similarly, we anticipate that uncertainties associated with the g a l f o r m modelling governing star formation and chemical enrichment will be reduced by comparing against future low-redshift observational data for faint galaxies, as well as high-redshift observations with, e.g., the James Webb Space Telescope [140], leaving little degeneracy with DM microphysics.
To conclude, we have shown that the binary black hole merger rate has the potential to become an important probe of deviations from the ΛCDM model, particularly the suppression of structure from dark matter interactions.As a by-product of our analysis, we have shown that dark matter-neutrino interactions with u νχ 10 −6 are already strongly in tension with the observed abundance of faint galaxies.With the next generation of gravitationalwave observatories probing the gravitational-wave event rate out to the earliest epochs of star formation, it will be possible to strengthen these constraints by another order of magnitude or more, giving us highly sensitive information about the earliest stages of nonlinear structure formation.

AUTHOR CONTRIBUTIONS
C.B., M.S. and A.C.J. proposed the idea for this project.M.M. performed the N -body simulations, halofinding, and halo mass function analysis with input from Y.W., S.B. computed the galaxy populations, and A.C.J. made the gravitational-wave predictions.All authors contributed to every stage of the paper, including the writing.

3 P 3 ]
Figure1.Linear matter power spectra for several DM-neutrino interacting scenarios of varying interaction strength uνχ, along with predictions of their "equivalent" FDM and thermal WDM scenarios characterised by their respective particle masses, max and m wdm , tuned to give a similar small-scale suppression.The FDM power spectra have been computed usingA x i o n C A M B[34,35]; the WDM power spectra are derived from the ΛCDM power spectrum filtered with a transfer function[10].Vertical lines indicate rough k-scales corresponding to the mass resolutions used in g a l f o r m.

Figure 2 .
Figure 2.An illustration of our pipeline.Linear power spectra of IDM scenarios are computed using c l a s s, which are then fed into g a d g e t and g a l f o r m as initial conditions.The g a l f o r m output is cross-checked with the halo mass function from g a d g e t and fed into c o m pa s, which computes the gravitational-wave event rate. includes

Figure 4
Figure 4.The suppression of the halo mass function, as described by Eq. 2, for a range of interaction strengths.The dotted lines show the mass resolution of our g a l f o r m merger trees for both the standard and high-resolution settings, while the dashed line shows the atomic hydrogen cooling limit, below which star formation is not expected[109].

Figure 6 .
Figure 6.IDM halo mass function n idm normalised to a reference CDM mass function n cdm for three interacting strengths uνχ.Solid lines denote results from N -body simulations; dashed lines represent our fitting function (2); dotted and dot-dash lines are two fitting functions from Refs.[110,111].The bottom panel of each plot shows the fractional residuals of the three fitting functions.

uνχ = 3 . 4 • 10 − 8 uνχ = 3 . 4 • 7 uνχ = 3 . 4 • 10 Figure 7 .
Figure 7.The binary black hole merger rate density over cosmic time, as predicted by our pipeline for ΛCDM and a range of IDM scenarios.The local z = 0 rate (90% C.L.) inferred from LIGO/Virgo observations[119] is shown as a green band.Curves marked "hr" are high resolution, meaning that the g a l f o r m merger trees are tracked to smaller progenitor masses.The curve marked "G.C." is the contribution from globular clusters in the noninteracting scenario.

Figure 8 .
Figure 8. BBH merger rate as a function of redshift and host halo mass, in both the (standard-resolution) CDM case and the uνχ ∼ 10 −7 case.

Figure 9 .
Figure 9. Expected number of BBH mergers observed with estimated redshift larger than z * , using one year of observations with a third-generation interferometer network (Einstein Telescope plus two Cosmic Explorers), taking into account the forecast redshift uncertainty and signal-to-noise for the expected BBH population.The shaded regions show the statistical uncertainty C.L.) due to Poisson fluctuations in the number of BBHs, and 'hr' means high resolution, signifying that the g a l f o r m merger trees are tracked to smaller progenitor masses.

Figure 10 .
Figure 10.Left panel: mean signal-to-noise ratio ¯ = (h|h) of equal-mass nonspinning BBHs detected with our thirdgeneration interferometer network, as a function of redshift z and total source-frame mass M .Right panel: forecast redshift uncertainty σz = (F −1 ) zz for the same parameter space.The mass range in both panels is truncated at the boundaries of our broken-power-law mass distribution (12 M < M < 174 M ).

1 IDM 1 IDMFigure 11 .
Figure 11.Formation rate (left panel) and merger rate (right panel) of BBHs in 21 different scenarios: a fiducial CDM model, an IDM model with uνχ ∼ 10 −7 and fiducial astrophysical settings, and 19 alternative CDM models in which various c o m pa s parameters are varied.

Figure 12 .
Figure 12.Covariance of the θi parameters defined in (13) based on the Gaussian Fisher matrix (15), assuming one year of observations with a third-generation GW interferometer network.