Cosmological inference using only gravitational wave observations of binary neutron stars

[Abridged] This study presents the first Bayesian investigation of the accuracy with which the cosmological parameters can be measured using information coming \emph{only} from the gravitational wave observations of binary neutron star systems by Einstein Telescope. We find, by direct simulation of $10^3$ detections of binary neutron stars, that, within our simplifying assumptions, $H_0,\Omega_m,\Omega_\Lambda,w_0$ and $w_1$ can be measured at the $95\%$ level with an accuracy of $\sim 8\%,65\%,39\%,80\%$ and $90\%$, respectively. We also find, by extrapolation, that a measurement accuracy comparable with current measurements by Planck is possible if the number of gravitational wave events observed is $O(10^{6-7})$. We conclude that, while not competitive with electro-magnetic missions in terms of significant digits, gravitational wave alone are capable of providing a complementary determination of the dynamics of the Universe.

Detectors beyond the second generation are already being envisaged. For instance, the Einstein gravitational * walter.delpozzo@unipi.it wave Telescope (ET) [24] is a proposed underground detector consisting of three 10 km arm-long Michelson interferometers in a triangular topology with opening angles of 60 degrees [25]. The strain sensitivity is estimated as factor 10 better than second generation detectors, down to frequencies of 1 − 3 Hz depending on the actual configuration of the instrument [26]. The high sensitivity promises the detection of a very large number of gravitational waves (GW) signals with large signal-to-noise ratios (SNR), thus allowing for unprecedented population studies as well as extremely accurate measurements of the physical parameters of coalescing binary systems [24].

A. Cosmological inference with gravitational waves
When estimating the parameters of GW sources, and in particular the coalescences of binary neutron stars and black holes, the luminosity distance can be observed directly [27,28]. This makes GW an ideal laboratory to place samples in the Hubble diagram in a manner that is arXiv:1506.06590v2 [gr-qc] 11 Jan 2017 free from the potential systematics affecting electromagnetic (EM) methods. Unfortunately, in the vast majority of the cases, the redshift cannot be measured from GW alone and this piece of information needs to be extracted by means other than GW.
In the recent years, various proposals have been put forward to overcome this difficulty. For instance, one can assume that the coalescences of compact binaries are the progenitors of short gamma ray bursts (sGRB) [29]. In this case, coincident observations of a GW event and the correspondent sGRB would allow the measurement of the luminosity distance from GW and the redshift from spectroscopy of the host galaxy [17,[30][31][32] 1 . For second generation interferometers, this method indicates a relative accuracy on the measurement of the Hubble constant H 0 of few percent in the case where 10-15 such events are detected. However, whether the coalescences of compact binaries are the progenitors of sGRBs is still a matter of debate. Also, the fraction of GW events also observable as sGRBs might be as low as 10 −3 [35] due to sGRB beaming effects.
An alternative approach, following broadly the argument first given in Ref. [27], would be to statistically identify the possible host galaxies of a GW event to obtain a distribution of possible redshifts associated with each GW detection. This method should yield ∼ 5% percent accuracy on H 0 using 20-50 events [18] observed by Advanced LIGO/Virgo. A similar methodology has also been applied to space-based detectors [36,37].
A few methods aim at extracting the redshift using GW observations alone. For example, one can use the knowledge of the (rest frame) mass function of NS and the measured (redshifted) mass to infer the redshift of the source [19,38]. In this framework, second generation interferometers should infer the Hubble constant H 0 with ∼ 10% accuracy using about 100 events [19].
The results of advanced interferometers can be greatly improved by third generation instruments such as ET. In fact, ET can probe regions of the Universe where the effects of the dark energy will be substantial, thus allowing an independent sampling of the cosmic history.
The potentialities of ET have already been investigated by various groups [30,31,39], concluding that, when only a limited set of cosmological parameters is considered, the accuracy of the inference is comparable to that of current EM measurements.

B. Outline
In this paper, we will expand on the approach proposed by Messenger and Read [40] in which if one of the 1 Kilonovae are also expected EM counterparts to BNS coalescences [33]. However their utility as cosmological probes is yet unclear due to their intrinsically faint luminosities [e.g. 34] which limits the distance at which they can be confidently detected.
two components is a NS, information about the equation of state (EOS) allows a direct measurement of the rest-frame masses and thus of the source redshift [40]. Using Fisher matrix formalism, the authors estimate the accuracy with which z can be measured to be ∼ 8 − 40%, depending on the EOS and on the distance to the source. Recently, a similar investigation was carried out in [41] using a more realistic Monte Carlo data analysis method. The authors concluded that the average uncertainty is closer to 40% for a hard EOS and essentially independent of redshift. Nevertheless, given the large number of sources that can be observed by ET and the possibility of combining information across them, even the large uncertainty reported in Ref. [41] could be sufficient to obtain interesting indications on the accuracy with which ET will measure the cosmological parameters. In this paper, we explore this idea in a simplified scenario and conclude that ET can indeed set bounds that are comparable to current EM measurement. We are interested in the cosmological information that can be inferred exclusively from the observation of gravitational waves. We will thus not discuss the potential of coincident GW-EM detections which are presented elsewhere [31,39]. We note here that, because of the co-location of the three ET interferometers and because of its topology, its expected sky resolution is extremely poor. Consequently, the probability of a successful EM-GW association is a priori very small. Note that at the time of ET, second generation detectors are expected to be operational with improved sensitivities [42]. For a substantial fraction of the loudest GW events, the sky localisation from a network made of ET and advanced detectors will be vastly improved compared to ET alone. In this case, some of the aforementioned EM+GW methods might become feasible and used to yield constraints on the cosmological parameters.
The rest of the paper is organised as follows. In footnote 1 we cast the problem in a Bayesian framework, and identify the necessary components to arrive at the cosmological inference. In Eq. (10) we describe the procedures of simulating GW events and the detector noise, and the implementation of the analysis. In Fig. 3 we present the results of our simulations and finally in table III we summarise and discuss our results. The mathematical solution to the problem of the inference of the cosmological parameters in the presence of a detection threshold is given for completeness in Appendix A.

II. METHOD
In this section we present a Bayesian solution to the problem of computing posterior probability density functions for a set of cosmological parameters from GW data. We broadly follow the presentation in [18]. Note that the treatment is not specific to the case of ET.
A. Inference of the cosmological parameters in the absence of a detection threshold Consider a catalogue of GW events E ≡ { 1 , . . . , N }. Each event is defined as a stretch of data d i (t) given by the sum of noise n i (t) and a gravitational wave signal where Θ i indicates the set of all parameters of the signal i.
The noise is taken to be a stationary Gaussian process with a zero mean and covariance defined by its one-sided spectral density S n (f ) such that where I represents all the relevant information for the inference problem, a tilde represents the Fourier transform, and where we have introduced a scalar product between two real functions A(t) and B(t) as The likelihood of observing the event i is then given by where S is the signal model that relates the signal parameters Θ i to a gravitational wave signal h. Moreover, the signal-to-noise ratio (SNR) ρ can be succinctly written as The posterior distribution for any parameter in our signal model S is related to the likelihood in Eq. (4) through the application of Bayes' theorem where p( Θ i |S, I) is the prior probability distribution for the parameters Θ i . When multiple independent detectors are included in the analysis, the likelihood (Eq. (4)) generalises to For this work, we are only interested in the posterior probability for a subset of parameters Ω ≡ {H 0 , Ω m , Ω Λ , . . .}. Therefore, we marginalise over the remaining subset of parameters θ i , i.e.
Where we have introduced the so-called "quasilikelihood" [43] Finally, the posterior for Ω given an ensemble of events E can be shown to be Therefore, in order to obtain the posterior for Ω, we need to perform a multi-dimensional integral in Eq. (9) for each of the GW events. The description of this procedure and the generation of data follow in Eq. (10).

III. ANALYSIS
In this section we describe the simulation that was performed. Firstly, we outline the generation of the data, consisting of the GW signal model, the astrophysical and cosmological assumptions regarding the source population, and the simulation of the detector noise. Secondly, we show the data analysis implementation with which the simulated data was analysed. In particular, we describe the construction of the quasi-likelihood, and it subsequent use to arrive at our cosmological inference. The GW signals, and the detector noise have been generated using the LIGO Analysis library (LAL) [44].

A. Astrophysical and cosmological assumptions
The NS masses are distributed according to a Gaussian distribution with a mean of 1.35M and a standard deviation of 0.15M [45] which is assumed constant throughout the cosmic history. For the NS equation of state we consider three cases; a hard EOS, a medium and a soft EOS. They are labelled as MS1 [46], H4 [47] and SQM3 [48]. We investigate these three cases since in [40] it was shown that the accuracy with which the redshift can be measured depends on the magnitude of the physical effects related to the details of the EOS. One can think of these three cases as an optimistic, a realistic and a pessimistic one, respectively.
The events are distributed uniformly in the cosine of the inclination, polarisation and time of arrivals. The events are also uniformly distributed in comoving volume. Their redshifts are sampled from the probability density given by [49] p(z| Ω) = dR(z) dz where R(z) is the cosmic coalescence rate. It is worth nothing that p(z| Ω) is an explicit function of Ω. The differential cosmic coalescence rate is equal to where r 0 is the local rate, e(z) is the cosmic star formation rate and V is the comoving volume. In a Friedmann-Robertson-Walker-Lemaitre (FRWL) universe, the rate of change of V with z is given by where we have introduced the Hubble parameter and the luminosity distance [50] H 0 is the Hubble constant, Ω m is the matter fractional density, Ω Λ is the fractional energy density of dark energy, is a convenient parametrisation to capture the effects of the redshift evolution of dark energy [51]. For Ω we chose fiducial values of where h = H 0 /100km·Mpc −1 ·s −1 . Even though ET horizon distance is 37 Gpc (z 4.15 for our fiducial cosmology), we limit our analysis to z max = 2 as this corresponds approximately the sky averaged horizon distance of 13 Gpc for BNS systems [52]. For simplicity, we decided to assume a star formation rate e(z) that does not change with redshift and is therefore irrelevant for our problem.
We simulated 1,000 binary NS events as observed by ET. The parameters θ i of each individual source have been generated according to the assumptions described in Eq. (10). The corresponding waveformh(f, Θ i ) was then added to Gaussian noise which is coloured according to the amplitude spectral density shown in Fig. 2. In Fig. 1 we show the network SNR distribution in the ET detector computed using Eq. (5).
Differently from most existing literature, we do not filter our sources with any SNR threshold. If we were to do so, we would be introducing a selection bias [53]. Note that, due to the potentially large number of sources observed by ET and their distribution in co-moving volume, the vast majority of them will in practice not be detected in a search which uses the SNR as decision statistics. It is known that ignoring these unregistered sources leads to a significant bias in the estimation of "global" parameters, see [53].
The main reason for the emergence of biases is intimately linked to the functional form for the prior on z Eq. (11). Since Eq. (11) quantifies the prior expectation regarding the distribution of sources in the comoving volume, it is an explicit function of Ω. The quasilikelihoods for the majority of our simulated events are almost uniform in Ω, see Section III C 2, therefore our inference is greatly influenced by the prior distribution: if one were to analyse sources that are louder than some threshold SNR the overall population of events would appear on average closer than the actual population. At the same time, the observed distribution of z would follow the actual cosmological distribution. Since p(z| Ω) Eq. (11) relates D L and z via Ω, this leads for estimates of Ω m → 1 and h → 0. Similarly, if one were to consider only events that are quieter than a given SNR threshold, Ω m → 0 and h → 1, thus Ω Λ → 1. In Appendix A, we give a mathematical solution to the problem of inferring Ω that accounts for sets of unobserved events. However, the solution in A is not computationally treatable with current techniques and for a large number of events, therefore we opted for an SNR selection threshold of 0 and analysed all simulated events. Furthermore, our choice relies on the capacity of distinguishing between low SNR GW signals and low SNR background events due to noise in the detector. A discussion of this problem can be found in [53]. The triangular configuration of ET provides an additional tool to study the distribution of signal and noise low SNR events. Thanks to its topology, ET admits the construction of a null stream which is devoid of any GW signal as the sum of the outputs of the individual Michelson detectors [25]. Being a pure noise process, the analysis of the null stream can be used to understand the SNR distribution of noise events which can then be used to infer the SNR distribution of quiet sources.

B. The signal model
In the previous paragraph we introduced the signal model S without specifying its properties. In this section, we lay out the assumptions that go in the construction of S.
In modelling the GW from a binary system, we limit our analysis to the inspiral phase of the coalescence process. We model the inspiral using an analytical frequency domain 3.5 post-Newtonian waveform in which we ignore amplitude corrections and the effects of spins. This is not a big limitation as NS are expected to be slow rotators [54]. In particular, we use the so-called TaylorF2 approximant [55], which can be written as where the waveform is written in terms of the amplitude A( Θ) and the phase Φ( Θ; f ).
The amplitude of the waveform A( Θ) is given by where we have introduced the chirp mass M = m 3/5 1 m 3/5 2 /(m 1 + m 2 ) 1/5 , D L is the luminosity distance defined in Eq. (15), (α, δ) signify the sky position of the source, and (ι, ψ) give the orientation of the binary with respect to the line of sight [55].
The wave phase can be written in the form where the ψ n are the so-called post-Newtonian coefficients (see e.g. [56]), which are functions of the component masses m 1 and m 2 , and (t c , φ c ) are the time and phase of coalescence. Note that all masses are defined in the observer frame, and the rest frame mass m rest is related to the observed mass through where z is the redshift of the GW source.
The description of the phase in Eq. (20) assumes that the object is a point particle, and thus cannot be tidally deformed. However, since we consider all of our events are binary NS coalescences, we modify the gravitational wave phase in Eq. (20) by including the finite-size contributions to the phase. These in turn depend on the tidal deformability λ(m rest ) [21] of the star which is a function of its equation of state and its rest frame mass. The finite-size contributions to the GW phase, as a function of observed masses, is given by where the sum is over the components of the binary, χ a = m a /M , λ a = λ(m a ) where m a are the component masses, M is the total mass, and η = m 1 m 2 /M 2 .
Knowledge of the EOS and using information encoded in the GW tidal phase contribution allows to measure the redshift of the source [40]. While the EOS is not known yet, various studies have shown that it could be possible to infer it from observations of BNS with second generation detectors [20-23, 57, 58]. In what follows, we will assume that the nature of the NS interior is known.

C. Data analysis
For our analysis, we assumed a noise curve for ET corresponding to the "B" configuration [59], corresponding to the projected sensitivity achievable with the current technologies Fig. 2. Given the anticipated rates of compact binary coalescences [13], the detection rates of binary NS systems in ET are expected to lie in the range 10 3 − 10 7 yr −1 [24].
The parameters θ of our signal model are the component masses m 1 and m 2 , inclination ι, polarisation ψ, right ascension α and declination δ, the time of coalescence t c , the phase of coalescence φ c , luminosity distance D L and the redshift z. In our analysis we ignored the presence of spins, as it is believed to be small in binary NS systems [54]. We analysed our ensemble of sources assuming that the EOS of the NS is known, thus accounting to a total of three analysis runs (one for each of our predefined hard, medium and soft EOSs). To obtain the posterior probability distribution on the cosmological parameters Ω we proceeded in three steps. Firstly, we analysed each source to compute a quasi-likelihood as a function of the redshift z and the luminosity distance D L . Secondly, these quasi-likelihoods are then converted into quasi-likelihoods as a function of the cosmological parameters Ω as shown in Eq. (9). Finally, the posterior probability function for the cosmological parameters given an ensemble of events are computed from Eq. (10).

Obtaining the quasi-likelihood
For each event i , we compute the quantity that is a partially marginalised quasi-likelihood, where the marginalisation is done on all parameters that are not relevant to the inference of Ω. These are λ ≡ (m 1 , m 2 , ψ, ι, φ c , t c , α, δ). The further marginalisation over z and D L will be described later on. For the time being, let us describe the details of the analysis for the computation of Eq. (23). The above integral was computed using a Nested Sampling algorithm [60] implemented similarly to what described in [61]. For each of the three analysis runs, we chose the same prior probability distributions for all parameters, with the exception of the component masses. For the common parameters we used uniform probability distributions on the 2-sphere for sky position (α, δ) and orientation (ψ, ι) and uniform in the time of coalescence t c with a width of 0.1 seconds around the actual coalescence time. For the first marginalisation, we choose uniform sampling distributions for both D L and z in the intervals [1, 10 5 ] Mpc and [0,2], respectively.
For the component masses, the priors were different across the different runs; each EOS in fact predicts not only the functional form of the tidal deformability λ(m) that enters in the phase of the GW waveform, but also the maximum permitted mass of the NS itself. Therefore, for the three EOSs under consideration we used maximum expected rest frame mass M max of 2.8, 2.0, 2.0M for MS1, H4 and SQM3 respectively. The prior probability distribution for the component masses was then chosen to be uniform between 1M and M max .

Cosmological inference
The marginalisation over the redshift and luminosity distance was then performed as follows: once a cosmological model is introduced, z and D L are not independent parameters anymore, but they are related unequivocally by Eq. (15), thus -after some algebra which can be found in [18] -we are left with the following integral to compute where p(z| Ω, SI) is given in Eq. (11) and we chose, consistently with the sources generation, z max = 2.
One of the problems we needed to overcome in order to perform the integral in Eq. (24) was how to represent L( i , D L (z), z, Ω) in a tractable way. In fact, one of the outputs of the Nested Sampling algorithm is a set of samples drawn from the integrand in Eq. (23) which is difficult to manipulate -in particular difficult to integrate -without making any assumptions about the underlying probability distribution.
A possible treatment of the problem would be to use the samples from Eq. (23) and approximate it using a normalised histogram. This procedure was successfully used in other unrelated studies [14], however, for our purposes it is not accurate enough. In fact, an histogram representation is dependent on a parameter, the bin size (or equivalently the number of bins once the range is specified), which cannot be inferred from the data but has to be chosen arbitrarily. The majority of the quasi-likelihoods in Eq. (24) tend to be very uniform over the cosmological parameter space for individual sources and, as noted in Section III A, the inference is strongly dominated by the prior on z. Most sources are close or below the detection threshold of the detector, thus a single source, in general, yields very little information about the underlying cosmology, therefore any small fluctuation in the histogram approximation due to the random variation of the number of samples in any specific bin would be amplified and would eventually lead to a biased estimate of the posterior probability density for Ω. As an example, compare the panels in Fig. 3. The left panel shows samples from Eq. (24) for a source having a network SNR=25. In this case, the isoprobability contours are almost consistent with a normal distribution. The right panel shows instead a source with network SNR=4. In this case the samples are almost uniformly distributed, thus if one were to approximate Eq. (24) with a two dimensional histogram, different choices of the bin size would result in different approximations, which would yield to very different inference of Ω. Instead, we decided to follow a different course of action. Given a set of samples, for any partition of the parameter space, the resulting probability distribution of the observations is always a multinomial distribution, therefore the "probability distribution" of the occurrences in each bin is a Dirichlet distribution. The above property defines a Dirichlet Process [62] which can be used to define an analytical representation of the underlying probability distribution of which we have only a finite number of samples available. For the mathematical details and definitions, the reader is referred to the original paper by Ferguson [62] or to the more recent discussions in [63]. We used the approximate variational algorithm [64] as implemented in [65] to find the Dirichlet Process Gaussian Mixture Model to represent the integrand in Eq. (24). The output of this procedure is an analytical representation of the target probability distribution as an infinite mixture of Gaussian distributions which is analytical and continuous. This form can then be used as the integrand in Eq. (24) and the integral itself can be evaluated using another Nested Sampling algorithm whose output can then be combined to compute the posteriors for Ω, Eq. (10).
For this last marginalisation, we assume priors on Ω that are uniform in all parameters. In particular, h ∈

IV. RESULTS
For each of the equations of state assumed we considered various scenarios within the set of (h, Ω m , Ω Λ , w 0 , w 1 ). We show posterior distributions for the cosmological parameters obtained from the joint analysis of 1, 000 events. We also report confidence intervals from number of events greater than 1, 000 obtained via extrapolation.
Moreover, since all EOS yield very similar results, we choose to report only posteriors for the cosmological parameters obtained from the MS1 equation of state.
In the following subsections we report on the values of the various cosmological parameters in each different cosmological model we considered. We note here that, as expected, the accuracy of the cosmological parameters measurement is better for the flat case and get progressively worse with the increasing dimensionality of the model under consideration. Also, all uncertainties we report are at the 95% confidence level.
Posterior distributions on the parameters of all cosmological models from the analysis of 10 3 BNS events are reported in Figs. 4, 5 and 6.
Depending on the actual astrophysical rate, ET will observe between 10 3 and 10 7 BNS events per year [24]. It is computationally unfeasible at the moment to simulate and analyse in a realistic way such a large number of events. We therefore extrapolated the results we obtained for 1,000 sources to the maximum expected number of events. We assumed that the central limit theorem holds, in other words that our posteriors for 10 3 events are approximately Gaussian, and simply scaled the variance of the one dimensional posteriors by the number of sources N . Tables I,II and III show the extrapolation of the 95% widths (2σ) to a number of events ranging from 10 4 to 10 7 for all relevant parameters for each of the cosmological models we considered in our analysis.

A. H0
We find that 10 3 BNS observations yield the following results: for the model (i) the accuracy is 0.05(7%) which remains approximately constant in the general case (ii) and worsens to 0.08(11%) in the general case (iii). In comparison with other GW studies, we find our measurements to be significantly worse. For instance [18] reports a 95% accuracy of ∼ 10% with second generation interferometers and similarly so do [32] and [38]. In comparison to traditional methods, the most constraining measurement from the Planck experiment in conjuction with other methods reports an accuracy of ∼ 0.5% [66]. So with 10 3 GW sources current measurements are far more accurate than those obtained with our method using ET. However, when extrapolating to the potential number of observable BNS systems, we find the Plancklike accuracy is reached with ∼ 10 5 BNS observations. A further order of magnitude improvement is seen for 10 7 BNS observations, see Table I.

B. Ωm and ΩΛ
We find that with 10 3 BNS events Ω m can be measured with an accuracy of 0.125(47.5%), 0.135(45%) and 0.19(65.5%) for the models (i), (ii) and (iii) respectively. The above numbers compare very unfavorably with the ∼ 2% yield by Planck [67]. However, if more sources are observed the current accuracy is reached with 10 6−7 sources, depending on the actual cosmological model.
The situation is similar for the measurement of the cosmological constant Ω Λ . With 10 3 BNS observations we find an accuracy of 0.23(37.5%) and 0.275(39%) for the models (ii) and (iii) respectively. As a comparison Planck reports ∼ 0.9% [66]. However, a similar uncertanty is reached with 10 6−7 sources. The results are summerised in Table II.

C. Dark Energy parameters
In the case in which the D L − z relation is modified to allow for a time varying cosmological constant, see Eq. (16), the parameters w 0 and w 1 were included in the model. From the analysis of 10 3 events, at 95% confidence, we find ∆w 0 = 0.8 and ∆w 1 = 0.9. In comparison, assuming 10 3 BNS events with optical counter parts and electromagnetic priors on the remaining cosmological parameters [31] finds ∆w 0 ≈ 0.1 and ∆w 1 = 0.3. Most recent determination by Planck [67] of the parameters w 0 and w 1 reports ∆w 0 ≈ 0.2 and ∆w 1 ≈ 0.5. Extrapolation to the potential number of sources observable in 1 year shows that the accuracy on w 0 and w 1 can be improved by 2 orders of magnitude, thus 1 order of magnitude better than the current best estimates, see Table III. It is worth noting that the posterior distributions for w 0 and w 1 are not very Gaussian, therefore the extrapolations to a large number of events might not be as reliable as for the other cosmological parameters.

V. DISCUSSION
In this study we investigated the potentialities of BNS observations with ET as cosmological probes. In particular, we quantified the cosmological information that can be extracted from pure GW observations of BNS. The ingredient that allows the measurement of the redshift is the knowledge of the NS equation of state and thus of the NS tidal deformability.
We simulated 1,000 events and relied on extrapolation to the expected 10 4 −10 7 events per year. The main result of this study is that from GW alone, ET could measure all cosmological parameters with an accuracy that is comparable with current state-of-the-art measurements from EM missions. This is the very first study of this kind, and therefore it should be regarded as a proof-of-principle. Our analysis General FRW+DE 0.9 × 10 0 2.9 × 10 −1 0.9 × 10 −1 2.9 × 10 −2 0.9 × 10 −2 relies on a set of simplifying assumptions which should be progressively relaxed for a comprehensive investigation. We neglected the NS spins and eventual precession of the orbital plane. We do not expect this to be a serious limitation, since NS are expected to be slow rotators [54], however it is not clear how a small degree of precession would impact the analysis of ET data. There is evidence that for binary black holes and NS black hole binaries accounting for precession of the orbital plane leads to more accurate distance measurements [68]. If a redshift measurement can be associated to these classes of sources, the determination of Ω would improve substantially, also thanks to the considerably larger volume that is observable by ET.
Due to its low frequency sensitivity, BNS signals in the local Universe (z < 1) will be in the ET sensitive band for a time scale of hours to days, depending on the observer frame masses. It is clear that this problem cannot currently be investigated using a fully realistic simulation since the generation of the most accurate waveforms and cutting edge data analysis algorithms are not yet sufficiently fast. A further complication arises from the number of signals itself; a detection rate of 10 7 events per year implies an average time delay between signals of ∼ 3 seconds. Given the duration of the signals in band, it follows that several signals would be present simultaneously in the ET data stream at any given time. No systematic study nor algorithm yet exists to investigate this problem.
We further assumed perfect knowledge of the NS EOS. While there are indications that second generation interferometers could measure the EOS [20][21][22][23], it is unlikely that we would know it without any form of uncertainty. However, we do not consider this a serious limitation as long as the error bar on any given value of the NS mass m and its tidal deformability λ(m) is sufficiently small to avoid confusion between different EOS. Moreover, we did not include in our analysis the merger part of the waveform, which, especially for the most distant sources, would be in the sensitive band of ET. The inclusion of the merger in the analysis would yield more information about the BNS mass and spin parameters as well as introducing a possible further constraint on its redshift [69], allowing for more precise measurements.
A further element that deserves future investigation is the effect of detection thresholds. We give a formal solution to the problem of the inference of the cosmological parameters in Appendix A. However, we did not explore the details of its practical implementation, which we defer to a future study. We ignored the effects of the detector calibration uncertainties over the inference of the GW event parameters as well their impact over the global inference of Ω. At the end of the first observing run of Advanced LIGO, typical amplitude uncertainties (which is relevant for the determination of D L ) and phase uncertainties (relevant for the estimation of z) were estimated at ∼ 10% and 10 degrees, respectively [70]. Simulations indicate that ignoring the presence of such calibration errors does not lead to significant bias in the estimation of the GW parameters, as long as the SNR is not very large [71]. However, data analysis models for GW analysis that marginalise over calibration uncertainties are now available [72] and routinely utilised for the actual analysis [10]. The additional calibration uncertainty increase the statistical uncertainty of the inferred parameters by a similar amount. In our case, the largest source of uncertainty would come from the amplitude calibration and thus on the determination of D L for the GW sources. Assuming, naively, an ET uncertainty budget of ∼ 10%, we estimate a similar degradation of our inference over Ω.
We also ignored the effects of weak lensing. Weak lensing is a zero mean process [73], thus, when averaging over thousands of sources, it will not induce an overall bias in the estimate of Ω. A proper account of the lensing uncertainty, would lead to similar conclusions as the detector calibration uncertainty.
Even with the caveats discussed above, our study shows that even considering exclusively information coming from GW alone (with no input from any EM association) ET is capable to fully probe the evolution of the Universe and determine the value of Ω with reasonable accuracy. We emphasise that our results apply to a pure GW-based inference of Ω. A more accurate determination of Ω from GW-EM joint detections may be possible, thus the results presented in this study should be regarded as a lower-limit to what the actual potentiality of ET is as a cosmological probe. Nevertheless, we showed that GW alone can be a feasible complementary and cross-validating route to probe the dynamics of the Universe.
For a set of M non-detected events, the likelihood will be given by The number of non-detected events M ( Ω, z) is a nuisance parameters which depends on Ω, the rate R( Ω, z), the detection threshold ρ th , the observation time T and the observed volume V ( Ω, ρ th ) as It is interesting to verify that Eq. (A7) reduces to Eq. (10) in the limit of ρ th → 0. In this limit we have Assuming that the rate R( Ω, z) is given by the integral of Eq. (12), we recover the form of the likelihood (10) which we used in our study.