{'fraction_non_alphanumeric': 0.029923903312444048, 'fraction_numerical': 0.018397493285586392, 'mean_word_length': 4.857498689040377, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Accuracy of inference on the physics of binary evolution from gravitational-wave observations

The properties of the population of merging binary black holes encode some of the uncertain physics of the evolution of massive stars in binaries. The binary black hole merger rate and chirp mass distribution are being measured by ground-based gravitational-wave detectors. We consider isolated binary evolution and explore how accurately the physical model can be constrained with such observations by applying the Fisher information matrix to the merging black hole population simulated with the rapid binary population synthesis code COMPAS. We investigate variations in four COMPAS parameters: common envelope efficiency, kick velocity dispersion, and mass loss rates during the luminous blue variable and Wolf--Rayet stellar evolutionary phases. We find that 1000 observations would constrain these model parameters to a fractional accuracy of a few percent. Given the empirically determined binary black hole merger rate, we can expect gravitational-wave observations alone to place strong constraints on the physics of stellar and binary evolution within a few years.


INTRODUCTION
Gravitational waves from binary black hole coalescences (Abbott et al. 2016a(Abbott et al. , 2017c have recently been observed by the ground-based gravitational-wave detectors of the Advanced Laser Interferometer Gravitational-Wave Observatory (aLIGO; Aasi et al. 2015) and Advanced Virgo (AdV; Acernese et al. 2015). These observations provide a revolutionary insight into the properties of the population of binary black holes. The catalogue of detections will grow rapidly as the instruments continue to improve their sensitivity (Abbott et al. 2018). In this paper, we analyse how such a catalogue will make it possible to infer the physics of binary evolution by performing inference on parametrised population synthesis models.
A number of channels for the formation of binary black holes have been proposed (see, e.g., Abbott et al. 2016c; E-mail: compas@star.sr.bham.ac.uk Miller 2016;Mandel & Farmer 2017, for reviews). In this study, we assume that all merging binary black holes form through classical isolated binary evolution via a commonenvelope phase (Postnov & Yungelson 2014;Belczynski et al. 2016). While all events observed to date are consistent with having formed through this channel (Stevenson et al. 2017a;Eldridge et al. 2017;Giacobbo et al. 2018), a future analysis would need to hierarchically include the possibility of contributions from multiple channels (e.g., Stevenson et al. 2017b;Zevin et al. 2017;Talbot & Thrane 2017).
Previous efforts to explore how stellar and binary population synthesis models could be constrained with gravitational-wave observations (e.g., Bulik et al. 2004;Bulik & Belczyński 2003;Mandel & O'Shaughnessy 2010;Gerosa et al. 2014;Stevenson et al. 2015) have typically focused on a discrete set of models, usually obtained by varying one evolutionary parameter at a time (e.g., Voss & Tauris 2003;Dominik et al. 2012;Mennekens & Vanbeveren 2014). In this paper, we consider the realistic scenario in which the astrophysical model is described by a multi-dimensional set of continuous parameters which may be strongly correlated. We ask how well we could constrain these parameters with a large observational data set.
The main tool we use to tackle this problem is the Fisher (information) matrix. Fundamentally, if we make an observation of a process, and we have a model for that process that depends on some parameters, then the Fisher matrix quantifies how much we can learn about the parameters in our model from the observation we made. We derive an expression for the Fisher matrix for binary-population synthesis. We use this to quantify how much we can learn about the population parameters from observations of binary black holes using the current generation of groundbased gravitational-wave detectors. While we concentrate on gravitational-wave observations here, the method is applicable to other data sets, and the best constraints may come from combining multiple complementary observations.
We use Fisher matrices to demonstrate that it may be possible to precisely measure the population parameters in binary-population synthesis models with ∼ 1000 observations of binary black hole mergers. At the expected rate of gravitational-wave detections (Abbott et al. 2017a), this could be within a few years of the detectors reaching design sensitivity (∼ 2-3 yr at design sensitivity for our fiducial model); the observing schedule for gravitational-wave observatories is given in Abbott et al. (2018).
We first give an introduction to our binary population synthesis model in section 2, together with a description of the model parameters we wish to infer using gravitationalwave observations. In section 3, we demonstrate how we transform the raw outputs of our binary population synthesis model by considering observational selection effects and redshift-and metallicity-dependent star formation rates. In section 4 we introduce the statistical tools used in this paper: (i) the likelihood function representing the probability of an observation given our model, (ii) a method for including measurement uncertainties in observations, and (iii) the Fisher matrix, which quantifies the sensitivity of our model to changes in its underlying parameters. The results of applying this methodology to binary population synthesis models are presented and discussed in section 5, and we discuss our conclusions in section 6.

POPULATION SYNTHESIS OF MASSIVE STELLAR BINARIES
Many of the details of binary evolution are currently uncertain (Postnov & Yungelson 2014;De Marco & Izzard 2017). Population synthesis models efficiently, albeit approximately, simulate the interactions of a large number of binaries in order to capture population wide behaviour and thoroughly explore the space of initial conditions. Uncertainties in the physics underlying isolated binary evolution are captured within population synthesis models through tunable parameters, which we call population parameters.
In this paper we focus on four population parameters which have an impact on binary black hole formation. We use the rapid population synthesis code COMPAS. 1 This uses the stellar evolutionary models of Hurley et al. (2000). Final black hole masses are calculated using the delayed model of . With the exception of the variations to the four population parameters we describe in section 2.1, we employ the Stevenson et al. (2017a) fiducial model throughout this paper.

Supernova kick velocity
The asymmetric ejection of matter (Janka & Mueller 1994;Burrows & Hayes 1996;Janka 2013) or emission of neutrinos (Woosley 1987;Bisnovatyi-Kogan 1993;Socrates et al. 2005) during a supernova can provide a kick to the stellar remnant. This birth kick is on the order of hundreds of km s −1 for neutron stars (Hobbs et al. 2005a). The typical strength of supernova kicks imparted to black holes is not well constrained observationally (Wong et al. 2014;Mandel 2016;Repetto et al. 2017), although they may be reduced relative to neutron star through the accretion of material falling back onto the stellar remnant .
In COMPAS, the strength of supernova kicks is parametrised using the dispersion parameter for a 3dimensional Maxwell-Boltzmann distribution σ kick . A kick velocity v kick is drawn from the distribution Alternative parametrisations for the supernova kick have been considered by Bray & Eldridge (2016), who did not find sufficient evidence to prefer them; here, we only consider continuous variations to model parameters, including the kick velocity in the Maxwell-Boltzmann distribution. The kick is modified to account for mass fallback, so that the final kick imparted to the black hole is where f fb is the fraction of matter that falls back on to the black hole, calculated according to the delayed model of Fryer et al. (2012). For carbon-oxygen core masses greater than 11M , f fb = 1 and so many heavy black holes receive no natal kick in this model (Belczynski et al. 2016;Stevenson et al. 2017a). Whilst observations of the proper motions of isolated Galactic pulsars (Hobbs et al. 2005b) suggest a value of σ kick = 265 km s −1 , we choose a fiducial σ kick = 250 km s −1 to match Stevenson et al. (2017a).

Common-envelope efficiency
When mass transfer is dynamically unstable and initially proceeds on the very short dynamical timescale of the donor, a shared, non co-rotating common envelope is formed around the donor core and the companion (Paczynski 1976). The details of the common-envelope phase are amongst the least well understood across all phases of isolated binary evolution (for a review, see Ivanova et al. 2013).
In COMPAS, the classical energy formalism (Webbink 1984) is employed to parametrise uncertainty in the physics of the common envelope. When a binary begins a commonenvelope phase, each envelope is bound to its core, with a total binding energy approximated by where G is Newton's constant, M core,(1,2) are the core masses of the two stars, M (1,2) and R (1,2) are the stellar masses and radii, respectively, and λ CE(1,2) are the corresponding stellar-structure parameters introduced by de Kool (1990) and are functions of star's evolutionary state (e.g., Dewi & Tauris 2000;Kruckow et al. 2016).
The loss of co-rotation between the orbit of the cores and the common envelope leads to energy dissipation which causes the cores to spiral in. Some of this lost orbital energy may be eventually used to eject the common envelope. The efficiency with which this transfer of energy occurs is uncertain, and is characterised by the free parameter αCE. In order to determine the separation after the common-envelope phase, the classical energy formalism compares the binding energy of the envelope to the energy transferred from the orbit ∆E orbit so that If the binary has sufficient orbital energy to completely expel the envelope, we consider this a successful common-envelope event. Unsuccessful ejections lead to a merger before a binary black hole system is formed. We follow Stevenson et al. (2017a) in assuming that common-envelope phases initiated by main sequence of Hertzsprung gap donors always lead to mergers (cf. the pessimistic model of Dominik et al. 2012). The fiducial choices of the parameters in COMPAS are λCE = 0.1 and αCE = 1.0. We explicitly leave λCE fixed whilst making small perturbations to αCE; however, this is an issue of labelling, since it is the product of these two free parameters which is ultimately of importance to the common-envelope physics ).

Mass-loss multipliers
Throughout their lives, stars lose mass through stellar winds. The wind mass-loss rate depends strongly on the star's luminosity and is generally highest for high mass, high metallicity stars. The dearth of observations of low metallicity environments means wind mass-loss rates are poorly constrained at low metallicities, and at high masses where stars are intrinsically rare. These are precisely the regimes where the progenitors of gravitational-wave sources are likely to form (Belczynski et al. 2016;Eldridge & Stanway 2016;Lamberts et al. 2016;Stevenson et al. 2017a;Giacobbo et al. 2018).
COMPAS employs the approximate wind mass-loss prescriptions detailed in Belczynski et al. (2010). For hot O/Bstars, we employ the wind mass-loss prescription of Vink et al. (2001). Our Wolf-Rayet wind mass-loss rates come from Hamann & Koesterke (1998). For other phases the mass-loss prescriptions from Hurley et al. (2000) are used. Uncertainty in mass-loss rates can have a significant impact on stellar evolution; for example, Renzo et al. (2017) find that there is a ∼ 50 per cent uncertainty in the mapping between initial and final masses when considering different mass-loss prescriptions when modelling solar-metallicity, non-rotating, single stars, with initial masses between 15M and 35M . There are particular phases of stellar evolution where the mass-loss rates lack strong constraints by observations. We parametrise the mass-loss rates in two of these phases with tunable population parameters.
During the luminous blue variable (LBV) phase (Humphreys & Davidson 1994), extremely massive stars undergo a relatively short episode of rapid mass loss which strongly impact the binary's future evolutionary trajectory (e.g., Mennekens & Vanbeveren 2014); observational constraints on the physics of LBV stars are currently uncertain (Smith 2017). 2 Following Belczynski et al. (2010), we parametrise this rate in terms of a multiplicative factor fLBV used to modify the basic prescription, so that the rate becomeṡ our fiducial value for this factor is fLBV = 1.5 (Belczynski et al. 2010).
During the Wolf-Rayet phase, stars have lost their hydrogen envelopes and have high but relatively poorly constrained mass-loss rates (Crowther 2007). We use a multiplicative constant fWR to modify the base rate: where L is the stellar luminosity, Z is the metallicity, Z = 0.02 is approximately the bulk metallicity of our Sun, and m = 0.86 is an empirically determined scaling factor (Vink &

MODEL PREDICTIONS
In this paper we evaluate the impact of the tunable parameters described above on the rate of detections and the measured chirp-mass distribution of binary black holes. The chirp mass M is a particular combination of the component masses M1, M2 which is measured well from the gravitational-wave frequency evolution during the binary inspiral (Cutler & Flanagan 1994;Abbott et al. 2016b), The chirp mass is just one of the parameters measurable through gravitational waves, other observables such as component masses, spins and the distance to the source can also be inferred (Abbott et al. 2016b). For simplicity, we have chosen to focus on chirp mass since it is the best measured. This is a conservative approach, as we have neglected information about other parameters; however, the methods presented here are easily extendible to include other observables.
In order to represent the distribution of chirp masses produced by the population synthesis model, we chose to bin our systems by chirp mass. Throughout this paper, we use 30 bins of equal width, ranging from the lowest to the highest chirp masses present in our dataset. The number of bins is determined by the scale length of variability in the chirp-mass distribution and the chirp-mass measurement uncertainty discussed below; the results are insensitive to halving the number of bins.
The raw output of a population synthesis model is a list of the initial conditions and final outcomes of all the binaries simulated. In order to compare this output to astronomical observations, it is necessary to process the data further, in order to account for the history of star formation in the Universe and the observational selection effects. We describe this processing below.

Cosmic history
In order to focus our computation on black hole progenitors, we only simulate systems with primary masses between 7M < M1 < 100M . We assume that all stars are in binaries with primary masses ranging between 0.01-150M following the initial mass function of Kroupa (2001) with a flat mass-ratio distribution (Sana et al. 2012). At formation, binaries are assumed to have a uniform-in-the-logarithm distribution of orbital separations (Öpik 1924;Abt 1983) and zero orbital eccentricity; for more detailed studies of massratio and orbital distributions, see Duchêne & Kraus (2013); Moe & Di Stefano (2017). COMPAS simulations produce a rate of binary black hole formation per unit star formation mass M form , where τ delay is the delay time, defined as the time from the birth of a binary to its coalescence (Peters 1964). To compute the total rate of binary black hole mergers per unit comoving volume per unit time we need to convolve the COMPAS formation rate with the amount of metallicityspecific star formation per unit volume per unit time at the birth of the binaries. Delay times can range from a few Myr to Gyr, and observations show that both the metallicity and star formation rates in galaxies evolve significantly over these timescales (Madau & Dickinson 2014). We use the star formation rate distribution of Madau & Dickinson (2014) and the metallicity distribution of Langer & Norman (2006). Other distributions have been suggested (e.g., Savaglio et al. 2005;Ma et al. 2016;Vangioni et al. 2015), and the cosmic history of metallicity evolution adds an additional source of uncertainty to our model predictions. Future studies could consider how metallicity evolution could be included with the other model parameters and inferred from binary observations. In figure 1 we provide an illustration of the metallicity-specific star formation rate at redshifts z = 0.5, 1, and 1.5, and also indicate metallicities at which we performed simulations for this study. We use these to The metallicity-specific star formation rate as a function of metallicity at three different redshifts, using the starformation-rate distribution of Madau & Dickinson (2014) and the metallicity distribution of Langer & Norman (2006). The vertical dashed lines indicate the metallicities at which we undertook simulations for this study. Metallicities above Z = 0.02 contribute negligibly to the binary black hole merger rate. translate the star formation rate into the merger rate at redshift z where ts is the time measured in the frame of reference of the merger, Vc is the comoving volume and we use cosmological parameters from Ade et al. (2016). Figure 2 shows the local merger rate at three different redshifts after accounting for changes in star formation rate and cosmology.

Selection effects
Gravitational-wave detectors are not equally sensitive to every source. The distance to the source, its orientation and position relative to the detectors, as well as the physical characteristics of the source all affect how likely it is that the system would be detectable. The detectability of a signal depends upon its signal-to-noise ratio (SNR). The SNR in a single detector is defined as (Finn 1992) where h(f ) is the waveform measured by the detector, S(f ) is the one-sided noise power spectral density, and fmin and fmax are the limits of the frequency range considered. For simplicity, we assume that signals are detected if their single-detector SNR exceeds a threshold value of 8 (Abbott et al. 2018). To model the waveforms, we use the IM-RPhenomPv2 (Hannam et al. 2014;Husa et al. 2016;Khan et al. 2016) and SEOBNRv3 (Pan et al. 2014;Babak et al. 2017) approximants; 3 these include the inspiral, merger and ringdown phases of a binary black hole coalescence, and allow for precession of the black hole spins. We incorporate the effects of cosmological redshift, which manifest as an apparent increase in the system masses, M obs = (1+z)Ms (Krolak & Schutz 1987;Holz & Hughes 2005). We assume a detector sensitivity equal to aLIGO in its design configuration (Aasi et al. 2015;Abbott et al. 2018).
We optimise our computations, reducing the number of waveform calculations required, by exploiting the fact that the parameters extrinsic to the gravitational-wave source, such as its position and orientation, only couple into the overall amplitude of the wave via F× ≡ 1 2 sin(2ψ)[1 + cos 2 (θ)] cos(2φ) where A, DL, i, ψ, θ and φ are the gravitational-wave amplitude, luminosity distance, inclination, polarization, and polar and azimuthal angles of the source location in the detector frame, respectively (Krolak & Schutz 1987;Cutler & Flanagan 1994). Therefore, we need only compute the phase evolution for a given combination of intrinsic binary parameters, such as masses, once, and then marginalize over the extrinsic parameters (with the exception of DL) as described in Finn & Chernoff (1993).
For a system with a given (M1, M2, DL), we determine the fraction of extrinsic parameter realisations for which the observed SNR passes our threshold, and label this as our detection probability P det .
We can use this detection probability to transform the merge rate given in Eq. (9) into a rate of detections. Integrating over the merger redshift gives the total detection 3 We use the implementations publicly available in the LAL suite software package wiki.ligo.org/DASWG/LALSuite. Total rate: 484 yr −1 Figure 3. The rate and chirp-mass distribution of the binary black hole coalescences we expect aLIGO to observe at design sensitivity, taking into account cosmic history and selection effects, for the COMPAS fiducial model as described in Stevenson et al. (2017a).
where ts is time in the source frame and t obs = (1 + z)ts is time in the observer's frame. Figure 3 shows the rate and chirp-mass distribution of binary black hole mergers detected at aLIGO design sensitivity. The mass distribution is shifted to higher masses relative to the intrinsic merger rate plotted in figure 2 because selection effects favour heavier systems which emit louder gravitational-wave signals. Some of the sharp features in this plot are the consequence of simulating systems on a discrete grid of metallicities (cf. Dominik et al. 2013). LBV winds tend to reduce high mass stars to a narrow, metallicitydependent range of black hole masses. We discuss the impact of these features in section 6.

The Fisher information matrix
The Fisher matrix quantifies the amount of information that a set of observable random variables (in our case, the merger rate and chirp-mass distributions) carries about the parameters (in our case, the four tunable parameters described in section 2) of a distribution that models these observables. Specifically, the Fisher matrix F for a set of random variables D (the data) which are dependent on a set of parameters {λ} is defined element-wise as where L is the likelihood function, defined as the probability of acquiring the observed data D given the model parameters, and the angle brackets indicate an expectation over the data realisation. We introduce the likelihood for our problem in the section below.
Under certain conditions, the inverse of the Fisher matrix gives a lower bound (the Crámer-Rao bound) on the covariance matrix for those dependent parameters (Vallisneri 2008); we discuss the regime of validity of the Fisher matrix inverse as an approximation to the covariance matrix in section 5.2. The covariance matrix tells us how sensitive our data is to a change in the model parameters. We can also examine which combinations of dependent parameters are degenerate and which combinations yield the greatest information gain. 4 The Fisher matrix quantifies the sensitivity of predicted observations to model parameters, and provides a bound on the accuracy of parameter inference. This approach assumes that the model is correct. The correctness of the model can be evaluated through other means. For example, model selection can be used to compare distinct models, whether these are different formation channels or different prescriptions for describing the physical processes of binary evolution (e.g., Mandel & O'Shaughnessy 2010;Vitale et al. 2017a;Stevenson et al. 2017b;Zevin et al. 2017;Talbot & Thrane 2017), or model-independent clustering can be used without reference to particular models (e.g., Mandel et al. 2015.

The COMPAS likelihood function
For this study we assume that we have a gravitational-wave catalogue of merging binary black holes, formed via the isolated binary evolution channel, and we focus on two observable characteristics of such a dataset: the rate of detections and the distribution of chirp masses for the observed systems.

(16)
The first term is the likelihood of observing binary black holes at a given rate. We assume that observable binary black holes coalesce in the Universe as a Poisson process with rate parameter µ, which is predicted by our population synthesis model, and total number of observations N obs , accumulated in a time t obs . The Poisson likelihood is log L (N obs |{λ}, t obs ) = N obs log(µt obs ) − µt obs − log(N obs !).
The second term is the likelihood of observing a given chirpmass distribution. As described in section 3, we have chosen to represent our chirp-mass distribution in bins. In this case the correct likelihood is a multinomial distribution (Stevenson et al. 2015) log 4 This is analogous to identifying the chirp mass as being the best measured combination of masses from gravitational-wave observations.
where K is the number of chirp-mass bins, c k is the number of observed systems falling into the k-th bin with k c k = N obs , and p k is the probability predicted by the model that a system falls into the k-th bin. Thus, µ and p k are functions of the tunable model parameters λ, while c k and N obs are observables. Given the likelihood, we can now calculate the Fisher matrix.

Computing the Fisher matrix
In order to compute the Fisher matrix, we need to find the second derivatives of the likelihood with respect to the population parameters and average over the possible observations drawn according to the same likelihood distribution. First differentiating the total-rate log likelihood, Meanwhile, differentiating the chirp-mass portion of the log likelihood yields The expectation value of N obs over this Poisson likelihood with rate parameter µt obs is just N obs = µt obs ; similarly, c k = µt obs p k . Therefore, the Fisher matrix is where we used k p k = 1 to eliminate the second term from Eq. (20). Crucially, this expression contains only first-order derivatives of the observables with respect to the population parameters. These derivatives can be readily and reliably estimated using population synthesis models, as described below.

Evaluating the first derivatives
We have shown in Eq. (21) that the Fisher matrix can be computed using just the first derivatives of the binned rates with respect to the population parameters. To compute derivatives, we simulated binary populations using a suite of variations to the population parameters discussed in section 2.1. We used the same set of random seeds to the random number generator in COMPAS, so that for each variation the initial conditions (i.e. masses and separation) and random effects (i.e. kick directions) remain fixed. This allows us to accurately measure the derivatives by estimating the differential rather than absolute rates, reducing the uncertainty associated with a limited number of simulations. We made six perturbations to the fiducial model for each population parameter (three negative and three positive). The perturbations were chosen to be sufficiently small that we could reliably estimate first derivatives numerically. A full list of the variations we used can be found in table 1. For each of the quantities we are differentiating, we have a set of overconstrained simultaneous equations for the first and second derivatives according to the leading terms in the Taylor series, which we can write in matrix form If we label the three terms in Eq. (22) as y, X and β respectively, then the maximum-likelihood solution for the derivativesβ can be computed directly as (Anton & Rorres 2000, section 9.3) We use this approach to compute all of the derivatives in Eq. (21) and combine them into an estimate of the Fisher matrix. The Fisher matrix can then be inverted to provide the Crámer-Rao lower bound on the covariance matrix of the astrophysical parameters evaluated at the COMPAS fiducial model.

Measurement uncertainty
The measurements of chirp masses will be subject to a certain amount of measurement uncertainty. We use a simplified treatment of this measurement uncertainty based on the methodology of Gair et al. (2010), see their appendix A. We assume that the probability of finding a system in an incorrect bin is given by a Gaussian distribution about the centre of the correct bin into which the system was placed in the simulation. Let fi be the fraction of system predicted by the simulation to lie in the i-th bin, which is centred on chirp mass µi and has left and right edges at chirp masses µ − i and µ + i , respectively. Then the probability pi of observing a system in the i-th bin is where σi is the standard deviation of the measurement in the i-th bin. In the limit of σi tending to zero, we recover perfect measurement accuracy, pi = fi. An illustration of this treatment of the measurement errors is presented in figure 4. The chirp-mass measurement uncertainty depends strongly on the total mass of the source, with the most massive sources spending the fewest inspiral cycles in band, leading to the largest measurement uncertainty (e.g., Abbott et al. 2016a). It also scales inversely with the source SNR. Here, we crudely approximate this as a fixed fractional uncertainty on the chirp mass of 3 per cent (cf., Vitale et al. 2017b). We therefore modify the binned rates according to Eq. (24), using a standard deviation σi = 0.03µi. (1); α CE is the efficiency of common-envelope ejection from Eq. (4); f WR is the multiplier for Wolf-Rayet wind mass loss from Eq. (6), and f LBV is the multiplier for luminous blue variable mass loss described in Eq. (5).
Our fiducial model appears in the top row. For each of these population parameter combinations we also varied metallicity. We used 12 different metallicities, which were evenly spaced in the log between 0.005Z and Z , where we use a solar metallicity Z = 0.02. We therefore had a total of 300 model variations. We simulated 1, 197, 989 binaries for each of these variations.  This method of incorporating measurement errors is a simplification. The formally correct approach would be to incorporate them on a per-system basis, which would involve a modification of the likelihood function. Performing the analysis in this way would correctly account for correlations between bins, whereas in the simplified approach bins are modified independently, losing information and slightly swelling the uncertainty.

Uncertainty quantification
The rate derivatives used to compute the Fisher matrix at the COMPAS fiducial model depend on the particular population realisation used in the calculation. We quantify the impact of simulation realisation noise, due to the finite number of simulated binaries, with bootstrapping. We recompute the Fisher matrix by re-creating data sets of the same size as the original simulated data set by drawing samples from it with replacement.
By repeating this process many times and observing the spread in the results, we can observe how much the point estimates change under different population realisations (different sets of binary initial conditions). Our full dataset consists of 359, 396, 700 binary simulations, which consists of the same set of 1, 197, 989 ZAMS binaries evolved under each of 300 different model variations (the 25 population parameter combinations listed in table 1, each simulated at the 12 different metallicities shown in figure 1). To generate one bootstrap sample Fisher matrix: (i) We randomly choose 1, 197, 989 initial conditions, with replacement, from our original set of initial conditions. (ii) For each of the 25 population parameter combinations in table 1, we find the systems from the bootstrap initial conditions which become merging binary black holes, and calculate their total rate and chirp-mass distribution (taking into account cosmic history, selection effects and measurement uncertainty). (iii) We use Eq. (22) and Eq. (23) to compute the derivatives of the total rate and chirp-mass distribution bin heights, with respect to each population parameter. (iv) We use these derivatives to compute the Fisher matrix, using Eq. (21).
We repeat the above steps 1500 times in order to express the uncertainty coming from the realisation of the initial conditions, i.e. from the simulation statistical fluctuations.
In principle, this model uncertainty could be overcome with more simulations, unlike the more fundamental uncertainties stemming from a finite number of observations and chirpmass measurement uncertainty. We discuss the relative contributions of these sources of uncertainty in section 6.

RESULTS AND DISCUSSION
Using the method described in section 4 we computed the elements of the Crámer-Rao lower bound on the covariance matrix for the population parameters σ kick , αCE, fLBV and fWR. We computed simulation uncertainties on these elements by taking 1500 bootstrap samples from the 1, 197, 989 sets of initial conditions simulated for the binaries, specifically varying the metallicities, initial masses and separations. Using these results we are able to explore what can be learned about these population parameters using gravitational-wave observations of binary black holes. Results are presented for N obs = 1000 observations, a sufficiently large number to ensure the validity of our results; we discuss the effect of changing the number of observations in section 5.2. Figure 5 shows the distribution of standard deviations of each of the population parameters. We see that it will be possible to measure αCE, fLBV and fWR with fractional accuracies of ∼ 2 per cent after 1000 observations. We will be less sensitive to the value of σ kick . This is an expected result, since the natal kicks of black holes are reduced according to Eq. (2), and many of the more massive ones do not get a kick at all.
The fractional uncertainties on all of the parameters are quantities of order N −1/2 obs ≈ 0.03 for N obs = 1000. Varying the parameters by their full dynamic range would change the rate by O(N obs ). For example, reducing αCE from 1 to 0 would make binary black hole formation through a commonenvelope phase impossible, reducing the expected number of detections from N obs to ∼ 0.
The measurement accuracy with which the tunable population parameters can be inferred using 1000 gravitationalwave observations can be alternatively interpreted from the perspective of model selection.For example, the median of the distribution for the standard deviation of αCE is ∼ 0.02. Therefore, if αCE different from the fiducial value by 6 per cent, the fiducial model could be ruled out with a confidence of ∼ 3σ ≈ 99.7 per cent.
We can examine the full multivariate normal behaviour of the population parameters. Figure 6 shows marginalised univariate distributions and bivariate projections of the 90 per cent confidence interval for each of the bootstrap samples. This plot shows that most pairwise correlations between most population parameters are negligible. Figure 7  shows the correlations between αCE and fWR, and between αCE and fLBV. Bootstrapping indicates an 88 per cent confidence that αCE and fWR are anti-correlated. Increasing αCE increases the efficiency with which orbital energy is transferred into the common envelope. An increased efficiency means that there will be less tightening of the binary, so fewer systems will come sufficiently close together to merge within a Hubble time. Losing mass through winds widens the orbit, meaning that increasing the Wolf-Rayet wind massloss rate creates more systems which are too wide to merge within a Hubble time. Increased mass loss also results in the black holes being less massive, therefore increasing the time required for them to merge through gravitational-wave emission from a given initial separation (Peters 1964). These correlations mean that increasing (or decreasing) both αCE and fWR would compound the effect on the rates, so their bivariate distribution (in figure 6) is narrower in this direction. Conversely, the effects of increasing one whilst decreasing the other would partially cancel out, and thus the bivariate distribution is wider in that direction. The confidence in the anti-correlation between αCE and fLBV is only 76 per cent, and there is insufficient evidence for correlation between other parameter pairs.

Information from the total detection rate
To gain further insight into the correlations between the inferred parameters, we now consider what we could learn about the population parameters by considering only the total rate at which gravitational waves are observed. It is impossible to constrain the four-dimensional population parameter vector considered in this paper with a single observable, the binary black hole detection rate. In this case, all that can be learned about the population parameters is the value of some linear combination of them. We construct a detection rate Fisher matrix, using only the total rate log likelihood of Eq. (17), and perform an eigendecomposition. We expect to see that there is only one eigenvector whose eigenvalue is non-zero. We verified that this is true for all 1500 of our bootstrap samples, which provided a useful sanity check of our results. Next, by examining the eigenvector whose eigenvalue is non-zero, we can find the linear combination of population parameters to which we are sensitive. Figure 8 shows a univariate representation of this direction (with its distribution from bootstrapping over simulations). The components of the vector parallel to fLBV and σ kick axes are broadly consistent with zero. Most of the information learned solely from the total detection rate is in the αCE-fWR plane. The fact that both values are simultaneously positive implies that they are correlated; this is the same correlation as was discussed at the beginning of this section.
Whilst we can only measure this specific combination of population parameters using only the total detection rate, we can constrain parameter combinations in the ∼ αCE + fWR direction to within a few per cent from the total rate. Figure 9 shows the standard deviation along the line defined by this combination of population parameters a −1/2 , where a is the principal eigenvalue. This can be interpreted in the same way as the standard deviations in figure 5, and matches the expected value of O(N −1/2 obs ). We see that if this combination of population parameters differed from our fiducial values by more than a few per cent, we would be able to confidently rule our model out after 1000 observations. However, we also see from figure 9 that including the chirp-mass distribution would significantly improve measurements of this parameter combination.

Standard Deviation Along Principal Eigenvector
Rates only Including M Figure 9. Distribution of the standard deviation of the particular linear combination of population parameters corresponding to the principal eigenvector of the total detection rate Fisher matrix. The measurement accuracy is computed using only information from the total rate (blue) and after including information from the chirp-mass distribution (green). The distributions come from considering all 1500 bootstrapped sets.

Number of observations
The expected number of observations only appears as a multiplicative term in Eq. (21), so that the standard deviations in figure 5 simply scale as N −1/2 obs . However, the results presented here are predicated on the assumption that the inverse of the Fisher information matrix is a good approximation to the covariance, and not just a lower bound. This in turn requires the likelihood to be approximately Gaussian, i.e. the linear single approximation (LSA; Vallisneri 2008) should hold. Only if the predicted parameter uncertainties are smaller than the neighbourhood in which the LSA is valid does the Fisher matrix provide a self-consistent estimate of the accuracy of parameter inference. This effectively sets a minimal threshold on the number of observations required for self-consistency in our estimates.
When computing the derivatives, as described in section 4.4, we measure the terms in a Taylor expansion of an observable (binned) rate f as a function of the population parameter λ, In order to verify the validity of the LSA, we need to check that each f is indeed linear when ∆ is of the same order as the computed standard deviations for the population parameters. We require that the linear term is dominant in the Taylor series, so that We find N obs = 1000 to be a sufficient lower limit on the number of observations necessary to ensure the LSA is valid. At 1000 observations, the best measured combination of parameters is constrained at the per cent level, and this will continue to improve as we expand the catalogue of observations.
For smaller numbers of observations, the LSA will break down. The probability distribution for the model parameters may no longer be a multi-dimensional Gaussian so the Fisher matrix is likely to under-estimate the inference uncertainty.

CONCLUSIONS
We have, for the first time, quantitatively analysed how accurately gravitational-wave observations of binary black hole mergers will constrain binary population synthesis models described by a multi-dimensional parametrisation. When ground-based detectors have accumulated 1000 observations of merging binary black holes, we have shown that we will measure binary population synthesis model parameters with an accuracy of a few per cent. Equivalently, we will be able to distinguish models for which the population parameters only differ by a few per cent.
Our analysis accounts for three distinct sources of uncertainty in the inference of population parameters using gravitational-wave observations. The first is due to the finite number of observations. We show when the linear signal approximation holds (section 5.2), the accuracy with which population parameters can be inferred scales with the inverse square root of the number of observations. The second is the chirp-mass measurement uncertainty in individual observations. We only model this approximately (section 4.5) but find that it is unlikely to be limiting factor in inference. The third source of uncertainty is simulation uncertainty: the accuracy in predicted detection rate and chirp-mass distribution is limited by the finite number of COMPAS simulations. This uncertainty, which we quantify with bootstrapping (section 4.6), is only limited by computational cost, and be reduced indefinitely with more simulations or more efficient sampling (e.g., Andrews et al. 2017).
There is, of course, potential systematic uncertainty in the binary evolution models themselves: for example, it is probable that the αCE parameter is not universal, as assumed here, but depends on the binary properties during the common-envelope phase. Model inference techniques such as those described here should be combined with model verification and with weakly modelled inference (e.g., . We show the expected detection rate and chirp-mass distribution of merging binary black holes in figure 3. The sharp features in the chirp-mass distribution are due to only simulating systems at a small number (12) of metallicities, replacing the integral over metallicity in Eq. (9) with a discrete sum. Mass loss, particularly during the luminous blue variable phase, leads to a pile up of black hole masses from the most massive stars at particular metallicity-dependent values. The subsequent discrete sum over metallicities overpopulates some bins in the chirp-mass distribution relative to neighbouring bins (cf. Dominik et al. 2013). This can impact our results, causing us to overstate the accuracy with which we will be able to measure population parameters. This issue can be addressed in the future by interpolating model predictions over metallicity (e.g., using Gaussian process emulators as described by Barrett et al. 2017), producing a smooth set of predictions.
Our primary intention with this paper was to introduce a methodology for evaluating the accuracy with which astrophysical model parameters can be estimated based on the rates and properties of observed transients. We considered a four-dimensional parameter space, but the number of dimensions is limited only by computational cost. It is also straightforward to add more terms than just the chirp-mass distribution to Eq. (16) in order to investigate other observable characteristics of binary populations such as mass ratios and spins (e.g., Stevenson et al. 2017b;Talbot & Thrane 2017;Zevin et al. 2017). Furthermore, this analysis can be used for other populations than observations of binary black hole mergers via gravitational-waves in this paper. Other observed populations, such as Galactic binary pulsars, Xray binaries, Wolf-Rayet stars, short gamma-ray bursts or luminous red novae (for a review, see De Marco & Izzard 2017), can provide orthogonal constraints on the parameters governing binary evolution (cf. figure 9). Over the coming decade, such measurements will help us to accurately determine the physics of massive binary evolution.