Host redshifts from gravitational-wave observations of binary neutron star mergers

Inspiralling compact binaries as standard sirens will soon become an invaluable tool for cosmology when advanced interferometric gravitational-wave detectors begin their observations in the coming years. However, a degeneracy in the information carried by gravitational waves between the total rest-frame mass $M$ and the redshift $z$ of the source implies that neither can be directly extracted from the signal, but only the combination $M(1+z)$, the redshifted mass. Recent work has shown that for binary neutron star systems, a tidal correction to the gravitational-wave phase in the late-inspiral signal that depends on the rest-frame source mass could be used to break the mass-redshift degeneracy. We propose here to use the signature encoded in the post-merger signal to deduce the redshift to the source. This will allow an accurate extraction of the intrinsic rest-frame mass of the source, in turn permitting the determination of source redshift and luminosity distance solely from gravitational-wave observations. This will herald a new era in precision cosmography and astrophysics. Using numerical simulations of binary neutron star mergers of very slightly different mass, we model gravitational-wave signals at different redshifts and use Bayesian parameter estimation to determine the accuracy with which the redshift can be extracted for a source of known mass. We find that the Einstein Telescope can determine the source redshift to $sim 10$--$20%$ at redshifts of $z<0.04$.


I. INTRODUCTION
The prospects for gravitational-wave (GW) astronomy in the era of advanced detectors are promising, with several detections expected before the end of the decade when Advanced LIGO [1], Advanced Virgo [2] and KAGRA [3] become fully operational. Among the sources of GWs expected to be detected are the inspiral and coalescence of binary neutron stars (BNSs), neutron star-black hole binaries, and binary black holes. Population models suggest that the detection rate of compact binary coalescences for BNSs will be ∼ 10 yr −1 , when Advanced LIGO [4] reaches its design sensitivity.
The inspiral of compact binary systems are also known as standard sirens [5], as their luminosity distance can be extracted from GW observations alone, without the need for any detailed modelling of the source, or of the properties of the media along the GW path. This is because the observed amplitude of GWs during the inspiral phase reaches the detector essentially unaltered and depends on a small number of parameters, which can all be measured using a network of GW detectors. These parameters include the total gravitational mass and mass ratio of the system, the spins of the compact objects, the orientation of the binary's orbital plane with respect to the line of sight, the source's position on the sky and the luminosity distance to the source. GW observations can very accurately measure the signal's phase evolution, which depends only on the total mass and mass ratio of a binary. Simultaneously, a network of detectors can determine the sky position, polarisation amplitudes and the distance to the binary. The observed total mass, however, is not the system's intrinsic mass M (i.e., mass as measured in the rest frame of the source) but the redshifted mass M z ≡ M (1 + z). This is known as the mass-redshift degeneracy.
The mass-redshift degeneracy is detrimental to the application of GW observations for cosmological inference. The relationship of the source's luminosity distance to its redshift on cosmological scales is precisely that which allows us to probe the parameters governing a cosmological model. Breaking the mass-redshift degeneracy requires an electromagnetic identification to tie the source to its host galaxy and thereby extract the source's redshift. It was thought, until recently, that there is no way to infer the source's redshift from GW observations alone.
To use of GW observations to extract information that is necessary for cosmography [e.g., estimation of the Hubble parameter and the dark energy equation of state (EOS)] and astrophysics [e.g., measurement of the masses and radii of neutron stars (NSs) and the EOS of matter at supranuclear densities], requires precision measurements of both the luminosity distance and intrinsic mass of the source. The massredshift degeneracy forces reliance on electromagnetic identification of host galaxies [6][7][8][9][10][11], which may be possible only very rarely. For example, using gamma-ray bursts for identification of the host galaxy greatly reduces the available signal population for cosmography and could potentially lead to observational bias. This is because gamma-ray emission is believed to be strongly beamed along a jet [12][13][14], while GW emission is expected to be approximately isotropic (quadrupolar) and hence only a small fraction (∼ 10 −3 ) of all GW events will have gamma-ray burst counterparts [8]. Determining the electromagnetic counterparts to binary-black hole mergers is also a very active area of research and several simulations have already been performed in this context [15][16][17][18][19] to provide first estimates on the properties and energetics of these emissions.
Some authors have explored other approaches to measure cosmological parameters without the aid of electromagnetic counterparts. These methods either assume that the intrinsic NS mass is constrained to be in a small range around [20][21][22][23] or make use of a statistical approach to measure the cosmological parameters from a population of events [24], as originally proposed by Schutz [5]. For BNS systems, however, there are two signatures in the GW signal from the late-inspiral and post-merger stages that depend on the rest-frame source mass and could potentially provide a measure of the source's intrinsic mass. Such a measurement would break the mass-redshift degeneracy and help return both the source redshift and intrinsic total mass from GW observations alone. The first effect is concerned with the correction in the orbital phase due to tidal effects that appear at order (v/c) 10 beyond the leading order in the post-Newtonian approximation to Einstein's equations [25], where v and c are the orbital and light velocities, respectively. This is a secular effect that becomes important as the two bodies approach each other and was first proposed as a cosmological tool by [26]. The second effect occurs in the post-merger stage and causes a significant departure in the post-Newtonian evolution of the system. Unless the two stars are not very massive, the newly formed object is a hypermassive neutron star (HMNS) [27], which develops a bar-mode deformation, which can survive even for a fraction of a second (cf. Fig. A1 of [28]), delaying the birth of the black hole and emitting GWs in a narrow frequency range. The importance of this effect has not yet been considered for cosmological exploitation.
In this paper, we propose to use the GW signal including the inspiral phase and the HMNS signature to extract the intrinsic gravitational masses and the source redshift. The power spectrum of the HMNS stage of the merger of BNS systems has been shown to contain prominent spectral features that vary smoothly in frequency with the total gravitational mass of the system [29,30]. The inspiral stage of the waveform can be used to obtain a highly accurate measurement of the redshifted total mass of the system M z . This allows us to constrain the true values of the rest mass and redshift to a relatively narrow band spanning the full range of the (z, M ) plane. Independently, the HMNS stage of the waveform allows us to measure the redshifted fundamental frequencies of two prominent spectral features to a reasonable accuracy. Using these and an empirically determined relationship between the total gravitational mass of the binary and the rest-frame fundamental frequencies, we are then able to independently constrain a second region of the (z, M ) plane. The localised intersection of the two regions in the (z, M ) plane allows us to break the massredshift degeneracy present in both measurements and make an estimate of the redshift and gravitational mass of the system. A cartoon of this idea is shown in Fig. 1, which can be directly compared to an example of the results of this analysis shown in Fig. 4.
The detectability of the HMNS part of the waveform is not likely to be high for sources observed in advanced detectors, as these features lie at frequencies significantly higher than the sensitive bands of ground based detectors. We will, therefore, explore how we might use the signature of HMNS in the context of the Einstein Telescope (ET), a third-generation ground-based interferometric GW detector [31].
The rest of the paper is organised as follows: In Section II we describe the numerical waveforms used for this analysis. In Section III we describe our robust, but ad-hoc, parameterisation and modelling of the HMNS power spectrum. In Section IV we describe the analysis methods used to simulate and measure the HMNS spectral features. We then describe the procedure with which these measurements are combined to obtain the redshift and gravitational masses of the source. Finally, in Section V we conclude with discussions of our results and future directions for this research.

II. NUMERICAL SIMULATIONS OF BNSS
All of our calculations have been performed in full general relativity. The evolution of the spacetime is obtained by using the CCATIE code, a finite-differencing code providing the solution of a conformal traceless formulation of the Einstein equations [32], with a "1 + log" slicing condition and a "Gamma-driver" shift condition. The generalrelativistic hydrodynamics equations are solved using the Whisky code [27,33], with the Marquina flux formula and a PPM reconstruction. For the sake of simplicity we model the NS matter as an ideal fluid with a gamma-law EOS, p = (Γ − 1)ρ with Γ = 2, where p is the pressure, ρ the rest-mass density, and specific internal energy (see [34] for details). The grid hierarchy, with a reflection symmetry condition across the z = 0 plane and a π-symmetry condition across the x = 0 plane, is handled by the Carpet mesh refinement driver [35], where we use 6 refinement level and  I. Properties of our initial data of equal-mass BNSs with the initial coordinate separation 45 km. Reported in the various columns are the baryon mass M b of each star, the ADM mass MADM of the system at initial, the gravitational mass M∞ of each star at infinite separation (M = 2M∞), the circumferential radius R∞ of each star at infinite separation, the compactness C ≡ M∞/R∞ and the orbital frequency f orb at the initial separation. the spacing of the finest grid is 0.15 GM c −2 ∼ 0.221 km. We extract the GWs, consisting of a plus and cross polarisation and sampled in time at a rate of ∆t = 1.68 GM c −3 ∼ 8.27×10 −3 ms equivalent to a sampling rate of ∼ 121 kHz, at a distance R 0 = 500 GM c −2 ∼ 738 km. We analyse only the = m = 2 mode of GWs, which is the dominant one. As initial data, we use quasi-equilibrium irrotational BNSs generated by the multi-domain spectral-method code LORENE [36] under the assumption of a conformally flat spacetime metric. We have considered five equal-mass binaries with an initial coordinate separation of the stellar centres of 45 km, and polytropic EOS, p = Kρ Γ with an adiabatic exponent Γ = 2 and polytropic constant K = 123.6 (in units c = G = M = 1); details on the different binaries are collected in Table I. A very important requirement of our sample of BNSs is that they are only very finely separated in total gravitational mass, with differences that are of the order of 2% only. Producing such a sample at a fixed separation is far from trivial and has represented a major numerical difficulty, stretching the capabilities of the LORENE libraries. Once evolved, the stars perform approximately 3.5 orbits before merger.

III. FREQUENCY DOMAIN MODELLING OF THE HMNS
In order to perform the parameter estimation described in Sec. IV we must first be able to parametrise and model the HMNS stage of the waveform. Using our five waveforms as a basis, the current state-of-the-art numerical simulations of BNS systems do not yet give us the insight and accuracy required to model the phase evolution of the HMNS waveform as a function of the system's mass. This, coupled with the assumption that there exists a smooth relationship between the total gravitational mass of the system M and the frequencies of prominent spectral features, forces us to model the signal power rather than the complex waveform. Therefore, in our HMNS analysis, we are insensitive to information encoded in the phasing of the waveform. Unless a semi-analytic description of the phase evolution in the HMNS stage is possible, the one adopted here is probably the only approach feasible.
For each numerical waveform we perform the following procedure in order to compute noise-free power-spectrum ref-erence templates. The time series for both the plus and cross polarisations are pre-processed using a fifth-order high-pass Butterworth filter with knee-frequency 1 kHz and a timedomain Tukey window with alpha parameter 0.25. This is done to suppress the leakage of power from the last few cycles of the inspiral and initial merger stage of the waveform. The discrete Fourier transform is then computed for each polarisation from which we construct the reference template where S h (f ) is the noise spectral density of the detector, which we choose to be that of the ET-B [37] design 1 . Visual inspection of these reference templates as a function of frequency, shown in Fig. 2, allows us to clearly identify the two primary spectral features of interest. The first feature, at frequencies ≈ 1.2-1.6 kHz is approximately Gaussian in profile and moves to higher frequencies for higher mass systems. In contrast, the second feature, at frequencies ≈ 1.7-3 kHz, appears to be best described by a sloping trapezoid with rounded shoulders and a central frequency and bandwidth that also grows with increasing system mass. In addition, there appears to be a third power component at low frequencies, i.e., ≤ 2 kHz, that becomes more dominant as the system mass increases. A reasonable approximation to this third feature is a second Gaussian of lower amplitude and greater variance than that used to model the first feature. For the purposes of this work, this third feature is included only to improve the quality of our model fitting. Mathematically, our entire ad-hoc model of the waveform power-spectrum can be expressed as where A 1 and A 3 are amplitude terms with F 1 , F 3 and W 1 , W 3 as central frequencies and half-widths (standard deviations), respectively, of the Gaussian features. The template is whitened using the detector noise spectral density as done for the reference template. We also define which is a linear slope of amplitude A 2a for f = F 2 − W 2 and amplitude A 2b for f = F 2 + W 2 . Finally, the function is the difference between two simple sigmoid functions which serves to bound our model of the second spectral feature component between the frequencies F 2 ± W 2 with a smooth transition from zero to unity over a fixed frequency range controlled by the parameter s. This ad-hoc template is therefore described by the 11-dimensional parameter vector λ ≡ We employ a simple least-squares fitting procedure to obtain our best fit parameters λ (M ) for each system mass. In all cases there were restrictions on the allowed parameter space ensuring that: (a) A 2a > A 2b , such that the slope of the second spectral feature was negative; (b) F 3 < F 2 , such that the third (Gaussian) spectral feature was restricted to the lower frequencies; (c) the smoothing length s < W 2 /5, such that the smoothed transition regions of the second feature account for less that 10% of the total feature width. As can be seen from Fig. 2, our choice of model and parameter restrictions provides a reasonable fit to the numerical data across our range of masses. We note that the quality of fit does begin to deteriorate at higher masses. This is due to the fact that as the mass of the binary increases, the HMNS is further away from a stable equilibrium and its dynamics is much more violent; in particular, the bar-deformed object rapidly spins up via the copious emission of GWs, leading to a very broad spectrum for the F 2 frequency.
The model fit depends on the parameter set λ, but we are only truly interested in a measure of the characteristic frequencies corresponding to the two dominant spectral features. For the lower-frequency Gaussian feature we choose to use the central frequency of the corresponding fit, F 1 , as this measure. For the higher-frequency, less symmetric, second feature we choose to define its characteristic frequency as the average frequency within the bandwidth of the second feature weighted by our best fit power model. This choice was made in an attempt to more robustly track the location of the power of the second feature. Hence, the lower and upper characteristic frequencies are defined as Taking the best fit parameters for each system mass and computing the corresponding f 1 and f 2 values, we indeed validate the initial hypothesis that there is a smooth relationship between these values and the total gravitational mass of the system. This relationship is shown in Fig. 3 where we plot total gravitational mass versus f 1 and f 2 . Also plotted are the following best-fit second and first-order polynomials for f 1 and f 2 , respectively, and whose expressions in Hz are tions represent our empirically determined relationship between the characteristic frequencies of the two spectral features and the total gravitational mass. These phenomenological frequencies essentially track the eigenfrequencies of the HMNS, which has been recently shown to behave as an isolated, self-gravitating system, whose dynamics can be described as a superposition of different oscillation modes (see [38] for details). Some remarks are worth making. First, expression (6) are clearly tuned to our choice of EOS, but equivalent expressions can be derived for any EOS [30]. Second, the functional dependence of the frequencies on the total mass of the system is only weakly dependent on the mass ratio; as a result, although our simulations refer to equal-mass systems, we expect the functions (6) to be a good approximation also for unequal-mass binaries (see also [29] where this is discussed in detail). Finally, as anticipated in the Introduction, the importance of the relations (6) is that they can be employed to break the mass-redshift degeneracy.

IV. ANALYSIS OF THE DATA
We now describe how we simulate the occurrence and subsequent measurement of BNS signals in third-generation GW detectors. We separate this process into two parts. The first is the simulation and measurement of the redshifted characteristic frequencies of the HMNS spectral features in the presence of Gaussian detector noise. The second is the independent measurement of redshifted mass parameters using the inspiral stage of a BNS waveform. We then show how these measure- ments can be combined to infer the gravitational mass and redshift of a source using the results of the previous section.

A. Measurement of the HMNS phase
The numerical waveforms comprise time series of the plus and cross polarisations of the GW signal at a distance of R 0 = 500 M which we label h  × (t) respectively. Using this data, we are able to simulate waveforms from BNS systems with arbitrary orientations, sky position, GW polarisation, phase and redshift.
We first compute the redshifted and distance-scaled polarisationsh where D L is the luminosity distance and is a function of the redshift z and a choice of cosmological parameters. We use a standard cosmological model described by the parameters: h 0 = 0.71, Ω m = 0.27, and Ω Λ = 0.73. In order to evaluate these waveforms at arbitrary redshifted frequencies we use a cubic spline interpolation scheme on the over-resolved discrete Fourier transform of the original time series data. To simulate GWs with arbitrary sky positions and nuisance parameters we recognise that the input polarisations are con-sistent with a circularly polarised wave, so that the signal at the detector is where φ, θ, ψ are the source right ascension, declination and polarisation angles, respectively. The cosine of the inclination angle is given by µ and we apply an arbitrary constant phase factor ϕ. Explicit definitions of the antenna patterns F (+/×) i , for the ET detector (where i indexes the three ET interferometers) can be found in [10]. Independent Gaussian noiseñ i (f ) of spectral amplitude matching the ET-B noise curve is then added to the frequency domain waveforms corresponding to each ET interferometer. The simulated data is theñ where i indexes the ET interferometers and k the discrete frequency grid on which the waveform has been evaluated. This grid is defined by the length and sampling time of the original numerical waveforms such that the frequency resolution ∆f = (N ∆t) −1 = 31.2118 Hz and only the 86 frequencies between 800 Hz and 3.5 kHz are included. At this point we apply a fifth-order high-pass Butterworth filter to the data with knee frequency 1 kHz and a time-domain Tukey window with alpha parameter 0.25. We therefore treat the data in exactly the same way as done in the template generation procedure and for the same reasons, namely, to minimise contributions to the signal power from the last few cycles of the inspiral and merger. It is also important to perform this filtering after the waveform has been redshifted and noise has been added, since in this way any spectral features resulting from this pre-processing will not show any artificial cosmological dependence.
We define the detector noise weighted power as At any given frequency, this power is governed by a noncentral χ 2 distribution with six degrees of freedom. The noncentrality parameter of this distribution is given by the squared optimal signal-to-noise ratio (SNR) in that frequency bin which is approximately equal to our power-spectrum template defined by Eq.
(2). We note that the freedom in our choice of amplitude parameters in the template allows us to use an equality rather than a proportionality in the above relationship. It then follows that the likelihood can be written as A1 A2a A 2b A3 F1,z F2,z F3,z W1,z W2,z W3,z min 0 0 0 0 1100 1700 800 30 150 30 max 100 100 100 100 1600 2700 1500 150 500 500 TABLE II. The range of priors of the model parameters used in the calculation of the posterior distributions on λz. Additional prior constraints defined by A1,z > A3,z, A2a > A 2b , and sz < W2,z/5 are also applied.
where I 2 is the second-order modified Bessel function of the first kind and L is the total number of frequency bins. We note that our template will now be sensitive to the redshifted frequencies present in the data and not the the intrinsic frequencies. Therefore, we define a redshifted parameter vector λ z containing the same amplitude parameters as λ, but with the frequencies F j,z ≡ F j /(1 + z), W j,z ≡ W j /(1 + z) and s z ≡ s/(1 + z). Consequently our parameters of interest become the redshifted characteristic frequencies f j,z ≡ f j /(1 + z).
We then apply a nested-sampling algorithm [39] to the data to provide samples from the posterior distributions of the λ z parameter set. We assume uniform priors on all parameters and use output posterior samples to generate posterior samples of the redshifted frequencies f 1,z and f 2,z . These samples represent the posterior probability distribution function (PDF) p(f j,z |{P }). We keep the prior parameter ranges of λ z constant for all simulations and therefore treat each simulation identically, independent of system mass and redshift. The range of our priors are given in Table II, where we also define the additional constraints that the amplitude of the first feature must be greater than that of the third feature; the slope of the second feature must be negative; and finally the smoothing length s z must be less than 10% of the total width of the second feature.
The resultant typical measurement uncertainties in f 1 and f 2 for each simulated system as a function of redshift are given in Table III. These intermediate results show, as expected, that the accuracy of measurement (reported in brackets) is O(few %) and decreases with increasing distance. It is also clear that there is a mild trend towards higher percentage uncertainties in higher-mass systems and that the percentage errors in frequency are comparable between the two spectral features.

B. Measurement of the inspiral phase
We treat the measurement of the inspiral stage of the waveform separately from the numerical simulations of the HMNS stage. At the redshifts relevant for the ET detector, the SNR from a BNS inspiral signal is after averaging over sky position, polarisation and inclination angles. We ignore the use of tidal information in the inspiral stage for redshift inference (we discuss this in Sec. V). Measurement of just the inspiral phase of the signal allows us to infer the redshifted total mass M z , but not separately We show pairs of results in units of Hz for f1 (upper rows) and f2 (lower rows) for each redshift and for each of the five total gravitational masses. Each value represents half of the span of the 68% confidence region on the frequency measurement averaged over 100 different noise realisations, source and sky orientations.
the gravitational mass or redshift. Given that the SNR is going to be very high, the redshifted-mass parameters will be very well constrained from the inspiral phase alone. To quantify the accuracy of this measurement, we use a Fishermatrix approach, which is a good approximation when the SNR is large and identical to that used in previous work on GW parameter estimation [40,41]. We consider as our signal model the frequency domain stationary phase approximation of a non-spinning BNS inspiral signal, correct to 3.5 post-Newtonian order. It follows that the redshifted chirp mass M z ≡ M z η 3/5 and symmetric mass ratio η ≡ m 1 m 2 /M 2 , where m 1 and m 2 are the component masses, has fractional measurement uncertainties of ∆M z /M z ≈ 5 × 10 −6 (z/0.1) and ∆η/η ≈ 10 −3 (z/0.1) for z < 1. This corresponds to a fractional uncertainty in the redshifted total mass of which is dominated by the uncertainty in the symmetric mass ratio of the system. In subsequent use of this expression we will not assume that the system consists of equal component masses. We also include a constant scale factor α and take a conservative approach by setting it equal to 10, therefore overestimating the measurement uncertainty of the redshifted total mass. In contrast to the HMNS stage, our analysis of the inspiral stage is quite simplistic. Given one of our original set of five numerically evolved binary systems and a choice of redshift, we assume a Gaussian uncertainty in the measurement of M z based on our Fisher matrix result. We then draw a random Gaussian variable M z , with mean equal to the original redshifted total mass and a standard deviation governed by Eq. (14) with α = 10, where a prime stands for the measured mass, which includes the measurement uncertainty, as opposed to the true redshifted total mass M z . The corresponding likelihood function is then a Gaussian with the same standard deviation and mean equal to the randomly drawn M z value.

C. Inference in the (z, M ) plane
Starting with the measurement of the inspiral phase described in the previous section, we can write the joint posterior PDF of the total gravitational mass and redshift as In the (z, M ) plane this defines a narrow "stripe" of probability spanning a region that extends from high masses at low redshift, to low masses at high redshift (see the blue curve in Fig. 1). We have assumed flat prior distributions for M and z.
For the HMNS measurement, the following procedure is applied to each spectral feature indexed using j. We note that any given redshift and total gravitational mass of a BNS system defines a specific point in the (z, M ) plane. From Eq. (6) we can determine the intrinsic characteristic frequency of the spectral feature in question at this value of M . We can also calculate the corresponding redshifted characteristic frequency using the relation f j,z = f j /(1 + z). The probability density associated with obtaining any particular value of f j,z , hence that particular (z, M ) pair, is given by the marginalised posterior PDF, p(f 1,z |{P }), obtained from our analysis of the HMNS data described in Sec. IV A. We can, therefore, write This function will describe an arc of probability in the (z, M ) plane that sweeps almost orthogonally to that obtained from the measurement of the inspiral phase. An increase in the observed redshifted frequency can be obtained via either an increase in the total mass of the binary or a decrease in redshift of the source. Consequently, a high-redshift high-mass system will generally have similar redshifted characteristic frequencies to a low-redshift low-mass system. Combining the information from both spectral features and the inspiral phase of the signal and assuming statistical independence, the final joint posterior distribution of z and M is simply the product of all three distributions We show examples of joint posterior probabilities of z and M for all five systems studied, and for multiple redshift values, in Fig. 4. These examples mirror the original conceptual sketch in Fig. 1, where we can see how, in practice, the massredshift degeneracy is broken through the use of the spectral properties of the HMNS GW signal. For sources at redshifts z = 0.01-0.04 the uncertainty in the measurement of the redshift is ∆z ∼ 10%-20%, over the full range of simulated system masses. In addition, the gravitational mass can be measured with fractional errors of < 1% in all cases.
In Tables IV and V we give representative values for the fractional uncertainty on the measured redshift and total gravitational mass, as a function of their true simulated values. These uncertainty estimates are obtained by marginalising the joint distribution p(z, M |M z , {P }) over the total mass to obtain p(z|M z , {P }), and over the the redshift to obtain p(M |M z , {P }). From these marginalised distributions we compute our representative uncertainties as half of the minimum interval to contain 68% of the total probability (analogous to the 1-sigma uncertainty for a Gaussian distribution).

V. CONCLUSIONS
A well-known problem of the detection of GWs from compact-object binaries at cosmological distances is the socalled mass-redshift degeneracy, namely, that GW measurements allow the determination of the redshifted mass M z = M (1 + z), but not the gravitational mass M or the redshift z separately. GW observations allow the measurement of the luminosity distance but this degeneracy restricts the cosmological application of GW observations, since it is the relation between the source's luminosity distance and its redshift that allows us to probe cosmological models. Until recently it was thought that coincident electromagnetic (EM) and GW observations would be required to break the mass-redshift degeneracy, as EM observations would allow the host galaxy to be identified and hence the extraction of the redshift. In this paper we have described a novel approach to this problem that does not require an electromagnetic counterpart and exploits instead the information encoded in the HMNS stage of a BNS waveform to break the mass-redshift degeneracy.
We have described how, with the use of five numerically generated BNS waveforms of very slightly differing mass, we have been able to construct frequency-domain powerspectrum reference templates. The templates were designed to capture the evolution of two primary spectral features in the HMNS stage of the waveforms, as a function of the total gravitational mass. The characteristic frequencies of these spectral features were then fitted to polynomial functions of mass providing us with an ad-hoc approximation to the characteristic frequencies for any mass. A Bayesian inference method was used to test the ability of the ET to measure the characteristic frequencies in the HMNS stage of the signal. These frequency measurements were coupled with our precomputed, empirical, frequency-mass relation and a measurement of the redshifted mass from the inspiral phase of the signal, allowing us to determine both the redshift and gravitational mass separately.
We have shown that in an analysis based on the signal's power spectrum, and ignoring all phase information within the HMNS stage, the measurement uncertainties in the redshift of sources at z = 0.01-0.04 is ∼ 10%-20%, over the full range of simulated system masses. In addition we find that the gravitational mass can be measured with fractional accuracies of < 1% in all cases.
We have specifically ignored the tidal effects in the late insprial phase that have been previously shown to be useful in redshift measurements. This choice was made to simplify our  4. The joint posterior distributions on the redshift and total gravitational mass of the BNS system for single representative realisations of noise and system parameters. Each row of plots represents a simulated signal of one of the five system masses (see Table I) ranging from low (bottom row) to high (top row) mass. Columns represent different simulated redshifts ranging from 0.01 (left) to 0.04 (right) in steps of 0.01. The green, blue and red contours represent the posterior contributions from the inspiral measurement, the first HMNS spectral feature and the second HMNS spectral feature respectively. The black contours represent the final posterior distribution combining all measurements and the black dots indicate the true simulated redshift and total mass values. In all cases the contours enclose 68% of the probability.
analysis and clearly identify the potential of this new technique. It is encouraging to find that the two approaches both have comparable redshift sensitivities of ∼ O(10%), implying that a combination of their results will improve the overall redshift estimate.
Under the hypothesis that there is a single universal NS EOS, it is highly likely that by the time of ET the NS EOS will be tightly constrained via various observations including direct GW detections from the advanced GW detectors. We have limited our study to a single EOS, but based on previous studies [29,42] would expect that the general result holds for all realistic possibilities. Of primary interest here is the general concept that there exists an additional "matter-effect" found in the HMNS stage of the waveform that can provide frequency markers from which redshift information can be obtained.
We should stress that this analysis is one of the first attempts to perform parameter estimation on the HMNS stage of BNS signals. Whilst the numerically generated waveforms and Bayesian parameter estimation techniques used here represent the state-of-the-art, our analytic signal model approximation and mass-frequency fitting is necessarily simplistic. The percentage measurement uncertainties on the redshift of a BNS source. We show results for each of our five different mass systems and for each of four different redshifts. We give three fractional redshift uncertainties for each combination. The first and second correspond to results using the inspiral measurement plus the first and second spectral features respectively. The third result (in boldface) is from a combination of the inspiral measurement and both spectral features. We have performed analyses of 100 different noise realisations, source and sky orientations for each redshift and mass combination. For each realisation we compute a quantity equal to half of the span of the 68% confidence region on the redshift measurement. The quoted value is the median of this quantity over the 100 realisations.  We expect that prior to the era of third generation GW detectors, the understanding of BNS mergers through numerical relativity and direct GW detections will enable us to significantly enhance our ad-hoc models. In the future, a significant improvement in the accuracies of redshift measurements using the HMNS stage could become possible if a realistic model of the phase evolution were constructed. In such a scenario, an analysis of the type presented here may be applicable to signals found in the advanced GW detectors.