Observing and measuring the neutron-star equation-of-state in spinning binary neutron star systems

LIGO and Virgo recently observed the first binary neutron star merger, demonstrating that gravitational-waves offer the ability to probe how matter behaves in one of the most extreme environments in the Universe. However, the gravitational-wave signal emitted by an inspiraling binary neutron star system is only weakly dependent on the equation of state and extracting this information is challenging. Previous studies have focused mainly on binary systems where the neutron stars are spinning slowly and the main imprint of neutron star matter in the inspiral signal is due to tidal effects. For binaries with non-negligible neutron-star spin the deformation of the neutron star due to its own rotation introduces additional variations in the emitted gravitational-wave signal. Here we explore whether highly spinning binary neutron-star systems offer a better chance to measure the equation-of-state than weakly spinning binary-neutron star systems. We focus on the dominant adiabatic quadrupolar effects and consider three main questions. First, we show that equation-of-state effects can be significant in the inspiral waveforms, and that the spin-quadrupole effect dominates for rapidly rotating neutron stars. Second, we show that variations in the spin-quadrupole phasing are strongly degenerate with changes in the component masses and spins, and neglecting these terms has a negligible impact on the number of observations with second generation observatories. Finally, we explore the bias in the masses and spins that would be introduced by using incorrect equation-of-state terms. Using a novel method to rapidly evaluate an approximation of the likelihood we show that assuming the incorrect equation-of-state when measuring source parameters can lead to a significant bias. We also find that the ability to measure the equation-of-state is improved when considering spinning systems.

the likelihood we show that assuming the incorrect equation-of-state when measuring source parameters can lead to a significant bias. We also find that the ability to measure the equation-of-state is improved when considering spinning systems.
Keywords: neutron stars, compact binary mergers, equation of state, spin quadropole (Some figures may appear in colour only in the online journal)

Introduction
On August 17, 2017, Advanced LIGO [1] and Advanced Virgo [2] made the first observation of a binary neutron-star merger [3]. This merger was associated with the short gamma-ray burst GRB170817A [4] and prompted a large-scale observing campaign to characterize the accompanying electromagnetic transients from the entire electromagnetic spectrum [4]. This observation firmly established the field of multimessenger astronomy, and demonstrated its potential to directly probe the physics of neutron stars. Neutron stars are exceptional environments where all four fundamental forces are simultaneously important and consist of matter compressed by their strong self-gravity to densities up to several times the density of an atomic nucleus. Despite much recent progress in theory, experiments, and observations, see e.g. [5][6][7][8][9][10][11][12][13][14][15][16], determining the composition and equation of state of neutron star matter remains a major objective at the forefront of fundamental physics [17] and astrophysics [18,19]. In the coming years, gravitational wave detections of many more binaries involving one or two neutron stars are expected, and will offer new ways to explore the internal structure of neutron stars.
One avenue for measuring the equation-of-state of matter in neutron star binaries is to look for deviations in the emitted gravitational waveform due to tidal coupling between the inspiraling neutron stars. The dominant effect is due to an adiabatic linear tidal interaction [20]. Unfortunately, this effect only becomes important when the stars are close to merger, since it scales as ∼k 2 (R/r) 5 , where k 2 is the tidal Love number, R is the neutron star's radius and r the orbital separation. For neutron stars the gravitational-wave frequency f ∼ 2f orbit ∼ GM/r 3 /π at which tidal effects become noticeable is a few hundred Hz, where the sensitivity of gravitational wave interferometers has begun to decrease. However, for spinning neutron stars there is another equation-of-state-dependent effect that scales as ∼QS 2 /r 2 , where S is the spin angular momentum and Q a dimensionless coefficient characterizing the star's rotational deformation away from the value of black holes Q BH = 1. This spin-quadrupole term, first described in [21], arises because a spinning compact object has a non-spherically symmetric mass distribution due to the rotational flattening at the poles, which distorts the gravitational field around the neutron star. This distortion affects the orbital evolution and the gravitational wave emission.
It is important to understand if these different effects of neutron star matter can be measured when using a gravitational-wave observatory like Advanced LIGO. It is also important to understand if neglecting such terms in current searches will reduce the number of observations that might be made in the coming years. This is because binary inspiral searches are based on cross-correlating the data with a bank of waveform templates, and unmodeled physics can potentially lead to a loss of signals due to inadequate templates. A number of works have addressed these topics [22][23][24][25][26][27][28][29][30][31], but in most cases did not include the spin-quadrupole term in the waveform 3 . The main reason why the spin-quadrupole term has largely been neglected is because of the expectation that a residual birth-spin of a neutron star would have decayed away long before entering the band of interest for ground-based gravitational wave observatories [33]. If the neutron-stars are non-spinning (slowly spinning) the quadrupole moment term is not present in (only a small contribution to) the gravitational wave signal and one must therefore rely only on the tidal effects to gain any information about the equation-of-state. However, it is possible that a neutron star can be spun up in a process known as 'recycling'. Recycled neutron stars have been observed as millisecond pulsars with spin frequencies as large as f spin = 716 Hz [34]. This translates into a dimensionless spin of order χ = S/m 2 = 2π(c/G) f spin I/m 2 ∼ 0.4 assuming the pulsar's mass and radius are m ∼ 1.4M and R ∼ 12 km, and with a moment of inertia I ∼ 1.4 × 10 45 g cm 2 , where the value of I and thus the inferred spin χ for a given rotation frequency depend on the equation of state (note also that the mass of the millisecond pulsar is not known). Such rapidly rotating neutron stars have not yet been observed in binary neutron star systems and it is not clear if binary neutron star systems where at least one of the bodies has such large spins will exist. One possibility is that they could form in dense stellar environments such as globular clusters or galactic centers through dynamical interactions [35]. However, if such systems do exist it may be possible to learn about the neutron star's internal physics by measuring the effect of the spin-quadrupole term, in addition to the tidal deformability term on the orbital evolution. Conversely, neglecting this effect might lead to a bias in all measured parameters, and potentially a loss of detected signals.
In this work we investigate whether deviations due to the neutron stars' equation-of-state would be observable in binary-neutron star systems observed with Advanced LIGO, allowing spins to be as large as χ ∼ 0.4. We particularly focus on the effect that the spin-quadrupole term can have and the effects that neglecting this term might cause.
The organization of this paper is as follows. In section 2 we briefly review the waveform model used in our study and the dominant effects of the internal structure of the neutron stars. In section 3 we give a brief introduction of the data analysis techniques used in this article. In section 4 we investigate the similarity between waveforms whose masses and spins are equal, but where the equation of state is allowed to vary. This does not give a complete picture though, because it is possible that a system with one equation-of-state could be modeled by systems with a different equation-of-state but also different masses and/or spins. Therefore, in section 5 we investigate the 'fitting-factor' between spinning binary neutron star waveforms with different equation-of-states to determine if changes in the equation-of-state can be hidden by changes in the system's intrinsic parameters. In section 6 we use a new method to efficiently evaluate the marginalized posterior probability distribution for the source parameters for four different example cases to determine the bias in masses and spins when neglecting the equation-of-state effects or using wrong values, as well as the improvements in equationof-state measurements for rapidly spinning binaries. In section 7 we briefly discuss how the redshift might be measured from equation-of-state terms. Finally we conclude in section 8. Unless otherwise specified we will use geometric units G = c = 1.

Waveform model including the quadrupolar spin and tidal deformations of neutron stars
We begin by discussing how the gravitational-wave signal emitted during a binary neutronstar merger would differ from that of a binary black-hole merger with otherwise identical source properties. For binary-neutron star mergers, most of the information in gravitationalwaves will come from the inspiral, and not the merger or post-merger, because the merger signal is emitted at frequencies too high to be easily observed with second generation gravitational wave observatories. Therefore we only consider the inspiral of the two bodies here.
In this work we will consider three effects that distinguish binary-neutron star mergers from binary black hole mergers: First, the effect of the rotational deformation of the components in the case that the components' spins are non-zero [21]. Second, the effect of the tidal deformation of the neutron-stars [20], and third the effect of the merger frequency of the binary [36]. Since we are interested in the dominant effects we neglect higher multipoles beyond the quadrupole [23], dynamical tidal interactions [37][38][39][40][41], gravitomagnetic tides [42,43], spin-tidal couplings [44,45], and the presence of a surface rather than an event horizon [46]. We note that there are some scenarios that predict more exotic, or more extreme deviations between a binary-neutron star merger and a binary-black hole merger [47], we do not consider such scenarios here. We also note that the most prominent imprint of neutron star physics on the gravitational waves is in the post-merger signal that can markedly differ from the ringdown of a binary-black hole merger, see e.g. the review articles [48][49][50] for details and references. However, we do not consider the post-merger signals in this work because they occur at frequencies too high to be observable with current facilities, and because while simulations of the post-merger signal do exist, this epoch is generically not yet well understood due to the complexity of the physics that becomes important.

Approximate frequency-domain description of GW signals from binaries
We will model the gravitational-wave signal emitted by two inspiralling neutron stars using the Post-Newtonian approximation [51]. Throughout this work we assume that the objects move on circular orbits and that their spins are collinear with the orbital angular momentum. In the post-Newtonian framework, the gravitational wave phase evolution is computed by imposing that the power radiated in gravitational waves is balanced by the change in binding energy of the binary. A number of different perturbative expansions, in powers of v/c where v ∼ GM/r is the orbital velocity, can be used to compute the phase evolution given the center-of-mass energy, currently known to O(v/c) 8 beyond the Newtonian result [52,53], and the gravitational-wave flux, for which O(v/c) 7 corrections to the quadrupole formula have been computed [54]. In this work we consider the frequency-domain 'TaylorF2' model describing the = |m| = 2 spherical harmonic mode of the waveform; for an overview of the different Post-Newtonian waveform approximants see e.g. [55]. TaylorF2 waveforms can be expressed in analytic form as Here h ( f ) denotes the Fourier transform of h(t), the time-domain gravitational-wave strain, denotes the chirp mass, D L the luminosity distance to the source, and θ x describes the various orientation angles that only affect the amplitude and overall phase of the observed gravitational waveform [56]. The phase Ψ is computed in the stationary phase approximation by integrating where E is the energy of the system and L GW is the gravitational-wave luminosity. The result of solving equation (3) perturbatively for small f and using the Post-Newtonian equations is where t c is the coalescence time and φ c is a constant phase offset. The coefficients λ i,j for nonspinning point-mass binaries are given in [57]. The first term, λ 0,0 , depends only on the chirp mass M but higher order post-Newtonian corrections also involve the symmetric mass ratio Spin effects first enter at O(v/c) 3 and are characterized at that order by the spin-orbit parameter Here, L is the unit vector in the direction of the orbital angular momentum and the χ i are the dimensionless spin parameters of each object where S denotes the spin angular momentum. At second post-Newtonian order, O(v/c) 4 , spinspin interactions start to influence the signal and are parameterized by Spin-dependent contributions to the phase appear again at higher orders in v/c. The explicit results for these contributions, currently known up to O(v/c) 7 , can be found e.g. in [57].

Equation-of-state effects in the gravitational-wave signals
For the purpose of our study we will consider two physical effects that lead to an imprint of the equation-of-state on the gravitational waves: spin-and tidally-induced deformations. For each of these we will focus on the dominant quadrupolar effect. As discussed above, the spin of the objects first enters equation (4) as an order (v/c) 3 correction to the leading order λ 0,0 -term, due to the coupling between the orbital angular momentum and the components' spin encoded in β given in equation (6). However, this term does not depend on the deformation of the objects and is thus independent of the nature of the object when expressed in this way 4 .
Finite size effects that depend on the underlying equation-of-state first enter the gravitational-wave signal as an order (v/c) 4 correction through the quadrupole-monopole interaction. This effect arises because a rotating neutron star is not spherically symmetric, which is physically a purely Newtonian effect despite the Post-Newtonian-like scaling with the frequency. The star's rotation causes a centrifugal flattening of its mass distribution into an oblate shape, which in turn creates a distortion in the gravitational field it generates. At large distances from the star, the leading order deviation of the Newtonian gravitational potential away from that of a nonspinning body is characterized by a quadrupole moment scalar where χ = |χ| is the magnitude of the dimensionless spin defined in (7) and Q is a dimensionless parameter characterizing the quadrupole deformation [21,58]. These Newtonian notions can be generalized to general relativity by considering the spacetime of a rotating neutron star and identifying the quadrupole moment from the asymptotic fall-off behavior of the metric potentials at large distances from the neutron star, or equivalently by using a more formal definition of multipole moments [58]. For a black hole, the quadrupole is given by the exact relation Q BH = −χ 2 m 3 , in accordance with the no-hair property. The difference between the black hole and neutron star quadrupole moments affects the orbital motion and rate of inspiral. This results in the following leading order contribution to the phase where The spin-induced deformation term enters again at higher orders. In our analysis we include the O(v/c) 6 term that is given for the case where both spins aligned with L by [59] λ QM 6,0 = The equation-of-state imprint in the gravitational-wave signals that has received the most attention in recent years is due to tidal effects. As the quadrupole-monopole term, these are Newtonian effects but they scale with the orbital velocity v as O(v/c) 5 and therefore become important later in the inspiral. The neutron star deforms in response to the companion's nonuniform gravitational potential across its mass distribution. Similar to the rotationally-induced buldge, the tidal bulges distort the object's exterior gravitational field, which in turn affects the orbital motion and gravitational wave emission. The dominant effect is characterized by a tidally induced quadrupole scalar of the form where E is the companion's tidal field. In Newtonian gravity, E ∼ −m comp /r 3 but this is generalized for relativistic systems to a definition in terms of the Riemann tensor characterizing the spacetime curvature produced by the companion. The coefficient Λ is the dimensionless tidal deformability parameter, which vanishes for black holes, Λ BH = 0. The adiabatic quadrupolar tidal effects give the following contribution to the TaylorF2 phasing [60]: λ 12,0 = 5(πMf 0 ) 7/3 512η 12/5 Here, Λ and δΛ are combinations of the individual tidal parameters given bỹ We note that for neutron stars, it was shown that the dimensionless parameters Q and Λ, characterizing the spin-and tidally-induced quadrupolar deformations of the neutron star's exterior spacetime, encode similar equation-of-state information and can be related in an approximately equation-of-state independent way [61,62].

Characteristic parameters
To illustrate the features of the parameters Q and Λ we consider representative examples of proposed equation-of-state models. The equation-of-state of neutron star matter has long remained a scientific challenge, despite theoretical advances and improved constraints from nuclear experiments and astrophysics. While most models largely agree up to densities around nuclear density, the large extrapolations required to apply known nuclear physics to the extreme conditions in neutron star interiors result in a wide range of possible equationsof-state. Among the numerous candidate equation-of-state models we consider two cases: one where matter is more compressible and the neutron star is compact (SLy, [63,64]) and a model where matter is stiff and the neutron star thus has a large radius for a given mass (MS1b [65]). We note that the latter model is already disfavored by the results from analysis of the first binary neutron star observation [3], but is not yet confidently ruled out, and for our purposes will serve as an upper bound on the size of the matter effects. The two different equations of state lead to different global parameters of the neutron star, as shown in figure 1 for the radius, rotational quadrupole parameter, and tidal deformability, as a function of the neutron star's mass. For comparison, we also included a few other equation-of-state models that assume a different composition such as hyperons or quark matter at high density, different values of parameters, or use different methods of calculation. As seen from the plot, the examples SLy and MS1b bracket a range of plausible equations-of-state.
An illustration of where the effects discussed above become important for signals observable by Advanced LIGO is shown in figure 2. The plot depicts the normalized accumulation of information about parameters per logarithmic frequency interval versus the gravitationalwave frequency for an equal-mass binary. The quantity shown in the plot is the normalized are the intrinsic source parameters and S n is the noise power spectral density, for which Advanced LIGO's zero-detuned high-power configuration [66] was used. The significance of this quantity for measuring the parameters will be explained in detail in the subsequent sections. It is the integrand for the diagonal elements in the Fisher information matrix, up to the factor of 1/f which converts to a logarithmic frequency interval. Each curve is normalized to its individual maximum value, except for the tidal parameter which is normalized by its value at a reference frequency of 1 kHz. We observe that for the mass-ratio and spin parameters (η, β, σ) the major contribution to the information comes from similar frequency ranges, while information about the tidal parameter accumulates at much higher frequencies. This is an important feature that we will return to when discussing our results. Note that in the Post-Newtonian waveform the symmetric mass ratio η first enters at a lower order than the spin parameters so that one might expect the information about η to be concentrated at lower frequencies. But because η also enters at all higher post-Newtonian orders, the distribution in the plot is shifted to higher frequencies than the spin parameters, for which only the leading order effect was included to generate this plot.
Additional information about the equation-of-state can come from the frequency at which the merger or tidal disruption occurs, where the latter is mainly relevant for mixed neutron star-black hole binaries. The TaylorF2 waveforms considered here describe only the inspiral portion of the signal and are usually terminated at the frequency corresponding to the inner-most stable circular orbit (ISCO) of a nonspinning system of a test particle orbiting a Schwarzschild black hole with the given masses [56]. We will also test the impact of not using the 'ISCO'-criterion but instead a fit from numerical relativity simulations for the merger frequency to terminate the inspiral signal [36]. 8 10 12 14 16 Radius (km)

A brief recap of binary-neutron star data analysis techniques
In this section we provide a brief recap of a number of the data analysis techniques that we will use in later sections in this work. For a more complete introduction to these topics we refer the reader to [67,68]. Consider a stretch of data s(t), recorded by a gravitational-wave observatory. This data is assumed to consist of colored, Gaussian noise n(t) with the possible presence of a gravitational-wave signal h(t). The noise is described by the one-sided noise power-spectral density S n ( f ), defined by where E[·] denotes the expectation value over independent noise realizations. We denote these assumptions of the noise properties with I. When evaluating the likelihood of a signal h(t) being present in the detector data, one can determine the probability of obtaining the given data realization if no signal is present, P(s|n, I), compared to the probability of obtaining the same data if a signal is present, P(s|h + n, I). These probabilities can be calculated according to [67,68] which reduces to P(s|n) in the case that h = 0. Here a|b defines a noise-weighted inner product according to where ã represents the Fourier transform of a. Then the relative probability of the two hypotheses is given by for S n the zero-detuned high power configuration of Advanced LIGO and each curve normalized to its maximum value.
This can be maximized over the unknown amplitude of the signal, A, to give the matched-filter signal-to-noise ratio that is routinely used in searches for compact binary mergers [56,69] When attempting to measure the parameters of a known signal in the data, a slightly different question is asked. We discussed in the previous section how a gravitational-wave signal from a binary neutron star merger depends on a variety of parameters, which we will collectively denote ξ i . When attempting to measure the parameters of a known signal one wishes to know what is the probability of the signal having a specified set of parameters, P(ξ i |s, I). According to Bayes' theorem this is given by The term P(s, I) in this application acts only as a normalization factor such that P(ξ i |s, I)dξ = 1. The quantity P(ξ i |I) represents the prior belief that some values of the signal parameters ξ i are expected to be more likely than others. Finally P(s|ξ i , I) represents the probability of obtaining a specific noise realisation s given a specific choice of parameters and our stated assumptions I. If one can evaluate P(ξ i |s, I) at all possible values of ξ i one has a direct measurement of the probability of different parameters. In general the parameter space is too large to allow a direct measurement and instead techniques to draw samples from the underlying probability distribution are employed [70]. An alternative, and much quicker, method to compute the expected bias in the peak of P(s|ξ i , I) due to underlying noise is to use the Fisher Information Matrix [71]. This is very quick to evaluate, but one must be careful when using this as it provides an estimation of the matched-filter between two waveforms (h(ξ i )|h(ξ i + δξ i )), which is only valid when δξ i tends to 0. For this to be valid for small, but non-negligible values of δξ i , as might be expected when estimating the parameters of gravitational-wave signals, the underlying parameter space metric must not vary strongly in the parameters used to evaluate the Fisher information matrix [71]. In this work we will instead attempt to measure P(ξ i |s, I) directly, making some assumptions to reduce the dimensionality of the parameter space, as we will discuss later in section 6.

Waveform mismatch with known masses and component spins
In this section we begin our exploration of the effect that equation-of-state dependent terms in the waveform model can have on gravitational-wave searches with a simple question. If we assume that all binary neutron star systems have the same-known-values of component masses and component spins, would it be possible to observe the difference between binary neutron-star systems with different equations-of-state, or a binary-black hole merger with the same component masses and spins. If the answer to this question were 'no' then equation-ofstate terms would not be possible to measure with observatories like Advanced LIGO.
To answer this question we note that the matched-filter signal-to-noise ratio between a filter waveform h and a stretch of data s containing noise n and a signal g is a linear sum Assuming the noise is Gaussian and stationary the average value of n|h over multiple noise realizations is 0, so one can consider a 'zero-noise' realization where n|h = 0 and neglect the noise contribution. Then the normalized matched-filter, or overlap, between the two waveforms gives the fraction of the optimal signal-to-noise ratio that is recovered when searching for a signal g using h as the waveform filter. Finally, we wish to maximize this quantity over the unknown source orientation and sky location. With the waveform model we are using, these parameters enter as a combination of an overall phase shift, overall amplitude shift and overall time-shift [69,72]. Therefore we define the 'match' as the overlap maximized over a phase shift and a time-shift, which is easily computed as described in [56,69] The value of this match at which signals would be distinguishable depends upon the signalto-noise ratio of the signal, as well as the geometry of the parameter space being considered and any strong priors being used. A simple rule-of-thumb, described in [73] argues that signals can be distinguished if the signal-to-noise ratio squared is reduced by an absolute value of 1 when searching for h 1 using h 2 as a filter, compared to the optimal signal-to-noise ratio where h 1 is used as the filter. For a signal-to-noise ratio of 15 this corresponds to waveforms being indistinguishable if the match is larger than 0.9978, or for a signal-to-noise ratio of 25, if the match is larger than 0.9992. We first use this match to determine whether second generation gravitational wave observatories will be sensitive to variations in the value of the spin quadrupole moment scalar, Q. This is done by generating two waveforms where the only difference is in the value of Q and computing the match between them. By repeating this procedure over a range of masses and spins we can evaluate where in the parameter space it might be possible to distinguish differences in waveforms due to variations in the spin quadrupole moment scalar. Here we set the tidal deformability parameter, Λ, of both bodies to 0 and for all waveforms use a termination frequency corresponding to the black hole 'ISCO'-criterion as described in section 2.
In figure 3, we show the match, as a function of the component spins, between waveforms with Q = 1 and waveforms modelled with either Q = 4 or Q = 12. These fiducial values are chosen as examples from the range exhibited in figure 1, e.g. Q ∼ 4 for the SLy model at m ∼ 1.6M , and Q ∼ 12 for the MS1b model at m ∼ 1.1M . The component masses here are chosen to be 1.35M for both bodies although we note that the plots look qualitatively similar when using different component masses and the same values of Q. From these plots we can see that the effect of the neutron-star self-spin deformation will have a negligible effect on the emitted gravitational-wave signal if the dimensionless spins of both bodies are less than χ = 0.05, as would be expected for non-recycled neutron stars. However, if we consider neutron-star systems with spins as large as χ = 0.4, as might be possible with recycled neutron stars, then the self-spin deformation causes very large mismatch between waveforms. Exotic compact objects can have much larger Q than neutron stars (see e.g. [74] for the case of boson stars), which would give a much more noticeable effect.
In figure 4 we show the match, as a function of the component spins, between waveforms modelled assuming that both bodies are binary black holes and waveforms modelled assuming both bodies are neutron stars described by a given equation of state. Here we use both the MS1b and SLy equations of state, described in section 2. These waveforms include not only the effect of the components' spin-quadrupole term, but also tidal terms and a difference in the termination frequency. In the top panels of figure  The difference in termination frequency reduces the match by only 0.0006 in the case of MS1b, which is much smaller than the contribution from the tidal deformation. However, the Λ and Q terms are strongly mass dependent, as seen from figure 1, and will be larger at lower masses, and smaller at higher masses. In the bottom panels of figure 4 we have computed matches for signals chosen with component masses uniformly distributed between 1 and 3 solar masses, and with component dimensionless spins distributed uniformly between -0.4 and 0.4. The results from this complete 4-dimensional parameter space is then plotted as a smoothed projection into the two dimensional parameter space of total mass and a mass-weighted spin term. There is some small variation of the match in the two dimensions projected away, which causes some noisiness in the smoothed plot, but the general trend is clear. Here we can clearly see that for increasing values of total mass and increasing values of the component spins the equation-ofstate dependent terms become increasingly important. Points close to 0 on the x-axis will have very little contribution from the spin-quadrupole terms, and the mismatch here is mainly due to the tidal deformation. The decreasing mismatch as the component spins increase is due to the increasingly important contribution of the spin-quadrupole terms.
These results demonstrate that the commonly-considered tidal terms are the dominant effect arising from the neutron-star matter when the component spins are small. However, if the spins are large, as would be expected for recycled neutron stars, the spin-quadrupole terms are the dominant equation-of-state related term and cause significant mismatches between otherwise identical waveforms. Therefore the spin-quadrupole terms must not be neglected when considering neutron-star systems with large spins.

Waveform mismatch with unknown masses and component spins
From our results for the observability of effects in the gravitational-wave signals discussed above and displayed in figures 3 and 4, we already drew three conclusions: equation-ofstate dependent effects can have a significant affect on neutron-star waveforms, for slowlyspinning binaries the tidal term is the dominant equation-of-state effect, and for rapidly spinning neutron stars the spin quadrupole term will have a much larger effect than tides. However, these results do not allow us to make the conclusion that the presence of the equation-of-state dependent terms will enable us to measure the equation-of-state in an observation of a binaryneutron-star system. The reason for this is that in general the masses and spins of a binary neutron star system are not known a priori, and so these must also be measured in combination with the equation-of-state related terms. In this section we assume that the component masses and spins are not known and ask if there would be a loss in the optimal signal-to-noise ratio if searching for binary neutron star systems using waveforms with incorrect equation-of-state parameters when the signalto-noise ratio is maximized over the unknown component masses and spins. This measure is important for two reasons, to quantify if there is a significant reduction in the obtained signalto-noise ratio if using the wrong equation-of-state parameters, and if searching for neutron star binary mergers with black hole templates will lead to a reduction in the overall number of binary neutron star mergers that would be observed. If it is not possible to find a waveform with a match close to unity when using an incorrect set of equation-of-state parameters and after maximizing over the unknown component masses and spins it would imply that we would be able to distinguish between the two equations-of state. In this section our focus will be on the question of how using wrong equation-of-state parameters could cause a reduction in the number of observations being made. The impact of neglecting equation-of-state dependent terms on detection rate has been studied in previous works in the context of non-spinning neutron-star systems [31]. We address here, for the first time, this question in the context of binary neutron-star systems with component spins as large as 0.4, where the spin-quadrupole and the tidal deformation terms are considered.
To be able to answer this question we need to calculate what fraction of the signal power is lost after maximizing over the mass and spin parameters. When searching for compact binary mergers a discrete set of waveforms, b i , is used [75][76][77]. The 'fitting factor' (as first defined in [78]) is the maximum overlap between the set of waveform filters and a potential signal waveform h (27) Normally in a search one creates the set of waveforms, b i , to fulfill the criterion that a signal anywhere in the parameter space being covered would have a fitting factor of at least 0.97. However, if signals are not contained within the parameter space being considered, for example if they contain equation-of-state terms not included in b i , the obtained fitting factor can be lower than the expected minimum. A standard practice to evaluate this ( [78][79][80] for example) is to compute the fitting factor for a population of signals and plot the distribution. However, this can sometimes be misleading as often the signals with the lowest values of fitting factors are also ones whose observable gravitational-wave strain is smallest. Therefore we also define, following [81,82], the 'signal recovery fraction' between a population of signals h j and a discrete set of filter waveforms b i where σ(h i ) = h i |h i is proportional to the observable signal power. In this way the signal recovery fraction gives the fraction of signals, described by the population h j , that would be recovered above an arbitrary signal-to-noise ratio threshold using the given set of filter waveforms b i , compared to a theoretical search that includes all possible values of h j in the set of filter waveforms. Or, in short, it quantifies what fraction of signals would be missed because of imperfect coverage of the filter waveforms.
To evaluate fitting factors and signal recovery fractions in this section we first create a set of filter waveforms b i . Here we create a set of filter waveforms using the methods described in [83]. These filter waveforms are constructed assuming that both bodies are black holes (which during the inspiral means spinning point-masses), with component masses between 1 and 3 M and component dimensionless spins χ ∈ [−0.4, 0.4]. The set of filter waveforms is constructed such that the maximal loss in signal-to-noise ratio for any waveform in this parameter space due to the discreteness of the bank is 1% 5 . This number is smaller than the commonly used value of 3% because, as mentioned above, we choose a smaller value here to emphasize the effect of the equation-of-state terms.
In figure 5 we show the signal recovery fraction between a population of black holes in this parameter space and our set of filter waveforms. As before the black holes have Λ = 0 and Q = 1 and use a termination frequency corresponding to the black hole 'ISCO'-criterion. We choose a distribution of component spins uniform in component spin magnitude between −0.4 and 0.4 (reminding the reader that we are restricting to only considering aligned-spin systems in this work). The sky location and orientation of the sources are chosen isotropically, and the signal-recovery fraction measure already assumes a uniform-in-volume distribution. The signal-recovery fraction is evaluated and plotted as a function of the two component massesthat is we choose a distinct set of signals and calculate signal recovery fraction for every point in the component mass space shown in the plots. In all cases in figure 5 we see signal recovery fractions larger than 0.98, which is expected as the waveforms are contained within the parameter space being considered-we are not yet including equation-of-state effects. We also show the fitting factor as a function of the two component spins when the component mass of both bodies is 1.35 M , illustrating that the fitting factor does vary because of the discrete nature of the bank. This figure provides the benchmark for the other plots in this section. There, any reduction in signal recovery fraction or fitting factor from that shown in figure 5 is due entirely to the effects of neutron-star physics that can not be recovered using black-hole waveforms. This would also correspond to the signal loss that is present in current searches for binary neutron stars in LIGO and Virgo data, where such terms are not currently included.
We first measure the signal-recovery fraction with a set of signals, which deviate from binary-black hole waveforms only by the inclusion of the quadrupole-monopole term. We perform two sets of simulations, one where Q is set to a value of 4 and one where it is set to a value of 12. In figure 6 we measure the signal-recovery fraction using the same distribution of signals as figure 5. We can see that in most regions in the component mass parameter space the signal recovery fraction is not noticeably lower than when using low-mass binary-black hole  waveforms. This tells us that in these regions of parameter space the spin-quadrupole term deviations in the waveform are almost completely degenerate with changes in the component spins and masses, as is also expected based on the discussion of figure 2. Current Advanced LIGO and Virgo searches would be able to observe signals that produce waveforms matching the ones used here. At the corners of the parameter space we do notice a significant drop in the signal recovery fraction. This is because the bank of waveform filters we used does not extend past the component mass and spin limits quoted above, in this case the signals here would match well with systems outside of the parameter space (i.e. with component masses <1 or >3M or component spin magnitudes >0.4). We also show the fitting factor as a function of component spins for systems with both component masses equal to 1.35 M and modeled with Q = 12. Here we notice a small reduction in the fitting factor only when the sum of the two component spins is large. We then measure the signal-recovery fraction using a set of signals modeled using the MS1b and SLy equations of state. These include terms for the spin-quadrupole, tidal terms and conditions on the termination frequency as discussed in section 2. We again use the same distribution of component spins, source orientation and sky location as in figure 5. The results of this simulation are shown in figure 7. We can see here that, especially at low masses, ignoring equation-of-state effects results in a signal-recovery fraction as low as 86% if assuming the MS1b equation of state, and as low as 93% if assuming the SLy equation of state. The biggest  To try to identify whether the loss in signal-recovery fraction shown in figure 7 comes primarily from the tidal deformation terms or from the spin-quadrupole terms we perform two additional runs with the MS1b equation of state. In one case we do not include the tidal terms  in the waveform, and in the second case we do not include the self-spin terms. These runs are shown in the middle panels of figure 7. We can clearly see from these plots that the drop in signal-recovery fraction when neglecting equation-of-state terms comes primarily from the tidal terms. This is expected given the degeneracy between the spin-quadrupole terms and the component spins and mass ratio observed in figure 6. From these results we therefore conclude that although the spin-quadrupole term can have a much larger effect on an emitted gravitational waveform than tidal terms, variations in the spin-quadrupole term are strongly degenerate with changes in the masses and component spins. Therefore the presence of spin-quadrupole terms for highly spinning binary neutron star merger waveforms is unlikely to cause a reduction in the number of binary neutron star signals that can be observed with Advanced LIGO. As already shown elsewhere [31], the presence of tidal terms can cause a small reduction in the number of observed signals. These conclusions are consistent with what is expected: tidal terms enter at much higher frequency and therefore cannot be easily mimicked by variations in the masses or spins, as illustrated in figure 2.

Parameter recovery with equation-of-state dependent effects
We have demonstrated in section 4 that the spin-quadrupole terms, often ignored in gravitational-wave data analysis, can have a significant effect on the emitted gravitational-wave signal for systems containing rapidly spinning neutron stars. In section 5 we demonstrated that while this effect is large, it is degenerate with changes in the mass ratio and component spins and one would be able to observe such systems well using only waveforms modeling both bodies as Kerr black holes. We find that the effect of tidal deformation, which only becomes relevant at higher frequencies, is a larger problem when thinking of observing such systems than the spin-quadrupole term.
However, when observing a binary neutron-star system one will also want to measure the parameters of the system, not only the neutron star equation-of-state but also the component masses and spins. Doing so will allow a much better understanding of how neutron stars form, how binary systems evolve, and provide a census of the properties of the astrophysical binary neutron star population. In this section we try to answer two questions. First, if the equationof-state terms are not included in the waveform model being used to estimate the system's parameters, or if an incorrect equation-of-state is assumed, by how much will the values of the parameters that are measured be biased? Second, is it possible to measure the equation-ofstate terms or to test if a specific observation is more compatible with one equation-of-state compared to another.

Methodology
To answer these questions we wish to evaluate which defines the probability of a signal being present in the data with parameters given by ξ i . When evaluated over all ξ i this defines the probability-density function over the whole parameter space being considered. An introduction of these concepts were given in section 3.
In this work we use a novel technique for evaluating P(ξ i |s, I). For non-precessing binary neutron star waveforms, there is only a weak coupling between the 'extrinsic' parameters of the system-the sky-location, distance, orientation and polarization phase-and the 'intrinsic' parameters-the component masses, spins and the underlying equation of state. Therefore we can make the approximation that analytically maximizing over the unknown extrinsic parameters of the system is equivalent to marginalizing over these parameters. The validity of this approximation is demonstrated in [84]. It is then possible to randomly pick a very large number of points, for each of the simulations described for this work we use 2.5 × 10 12 points, from the underlying distribution given by P(ξ i , I)-restricted to only intrinsic parametersand calculate P(ξ i |s, I) for all points assuming the 'zero-noise' realisation as motivated in section 3. The Fisher information matrix is used to predict the match between each of the points being considered and the parameters corresponding to the true source, using the implementation discussed in [83]. For all points where the Fisher Information Matrix predicts that the match is not negligibly small, P(ξ i |s, I) is calculated numerically using the PyCBC software package [85][86][87]. At all other points P(ξ i |s, I) is assumed to be 0. In this way we can rapidly evaluate P(ξ i |s, I) for the 'zero-noise' realization.

Results
We begin by exploring the parameter bias that occurs if searching for binary neutron star systems using waveforms where the value of the spin-quadrupole term differs from the signal we are looking for. We do this for systems that have a number of different values of masses and component spins. In all cases the signal we are looking for is assumed to have a spin-quadrupole value of Q = 8 on both bodies, and we try to recover this signal assuming Q = 1, 4, 8 or 12 on both neutron stars. This allows us to understand how the bias that will be present in measuring parameters varies as the error on the spin-quadrupole value changes. For this simulation we neglect tidal terms and use a termination frequency corresponding to the binary black hole 'ISCO'-criterion. The results of this are shown in figure 8. Here we show results for four different systems, with the details of those systems given in table 1. We use a dimensionless spin of 0.35 to model the signal in many cases in figure 8. While the bias is largest for systems when the binary neutron star spins are large, there is a visible bias even when the source has no spin on either body. The reason for this is that the signal from a non-spinning binary-neutron star system can match well with a system with non-zero spins. As the signal from spinning systems is altered by the value of the spin-quadrupole term we still observe a bias if we allow high-spinning systems in our prior. We also show the posteriors marginalized over the full parameter space for all cases in table 1. Here the majority of information is coming from the priors, and boundaries that we have placed on the parameter space, rather than from the data. For example we notice that in all cases with component spins of 0.35 the Q = 12 case is strongly disfavored. This is because we have assumed that the prior probability on component spins is flat up to a value of 0.4 and then drops immediately to 0. For Q = 12 we require spins larger than 0.4 to match well to the assumed signal model, which are not permitted, and therefore strongly disfavored. Likewise, in many cases values of Q = 4 are favored above the correct value Q = 8, again this is because the 2-dimensional probability plots shown in figure 8 intersect the boundary of the parameter space to a greater degree with Q = 8 than with Q = 4. In short we are not able to measure the value of Q in any of these simulations.
An additional test of the parameter bias is shown in figure 9 and table 2. In contrast to figure 8, here we use the fits to the various equation of states described in section 2 to evaluate the parameter bias that would be present if we search for a binary-neutron system described by one equation-of-state using waveforms modelled by another. In addition to the spin-quadrupole term, these waveforms include equation-of-state and mass-dependent terms describing the tidal deformation and the termination frequency. In these results we model the astrophysical The blue shaded region denotes the 99% and 99.99% confidence region when using the correct value of the spin-quadrupole term, which is assumed to be Q = 8 here for all simulations. The red shaded regions denote the 99% and 99.99% when assuming incorrect values of the spin-quadrupole term, shown here are results for Q = 1, 4 and 12. As simulation ID 3 is modelled with both component spins of 0, the bias here is much smaller than for other cases. Here, for clarity, we only show the 99% confidence regions for Q = 8 and Q = 1. Table 1. Parameters, and marginalized likelihood, for the runs plotted in figure 8. m 1 and m 2 denote the two component masses in units of component masses. χ 1 and χ 2 denote the two component masses. ρ denotes the signal-to-noise ratio of the signal. L Q=x gives the marginalized likelihood when searching for the signal assuming the spin-quadrupole term, Q, is x for both bodies. In all cases the signal is modelled using Q = 8 and the marginalized likelihoods are normalized so that L Q=8 = 1. Then the ratios of these values give the relative posterior likelihood between the various models. All simulations here assume the Advanced LIGO zero-detuned, high-power noise curve.  Figure 9. Two-dimensional marginalized posterior probability distribution of effective spin and the symmetric mass ratio, η, for the four simulations listed in  Other than that we use the same priors as described above for figure 9. As with figure 9 we see that assuming the incorrect equation-of-state when measuring source parameters can lead to a significant bias. In all cases we can see that the permitted range of mass ratios is much more constrained when using the correct MS1b equation-of-state than when using SLy or assuming two Kerr black holes. This is because the tidal terms become increasingly large on the smaller body as the mass ratio increases for MS1b and this breaks the mass-ratio and spin degeneracy much more easily. We also show the marginalized likelihoods for the various cases in table 2. We can see that for the two high-spin, high signal-tonoise-ratio cases the correct equation of state is strongly favoured over SLy or the Kerr black hole model. For the non-spinning, high signal-to-noise ratio case the discrimination power is much weaker than for the spinning cases. This implies that measuring the equation-of-state for a spinning binary neutron star will be much easier than for a non-spinning binary neutron star. Finally, when reducing the signal-to-noise ratio to 12, as might be expected for many of the binary-neutron star merger detections, we find that it is not possible to distinguish between the various equations of states.

Redshift effects
It has been pointed out in previous works that measuring the tidal deformability terms Λ i precisely offers a way to directly measure the object's masses and hence the redshift of the source [88]. The reason for this is that gravitational-wave signals are redshifted in an analogous way to electromagnetic signals. In the case of binary-black hole mergers the redshift effect is entirely degenerate with the total mass of the system. That is to say that the gravitational-wave signal observed from a redshifted source will appear exactly the same as that from a non-redshifted source with larger masses by a factor of (1 + z), where z is the redshift. However, in the case of binary neutron stars the tidal deformability breaks this degeneracy. The reason is that Λ i depends on the source-frame mass as Λ i ∼ m −5 i , where a more detailed approx imation involves a linear expansion around a reference value Λ(m = 1.4M ). Therefore one can simultaneously measure the redshifted mass and the tidal deformability of a system, and if the relationship between neutron-star mass and tidal deformability is well understood, one can then determine the redshift of the system because the tidal terms will involve an explicit factor of (1 + z) 5 .
Similarly, the presence of the quadrupole-monopole term also breaks the degeneracy between total mass and redshift, and a precise measurement of the quadrupole-monopole term could likewise allow a determination of a source's redshift. The quadrupole parameter scales with the source-frame mass as Q i ∼ m −1 i , and an analogous reasoning as for the tidal terms applies. However, the results in this paper demonstrate that the quadrupole-monopole term is not measurable to the same accuracy as the tidal deformation term, as changes in the quadrupole-monopole term are degenerate with changes in the component spins and mass ratio. Therefore the ability to measure the redshift of a binary-neutron star system will still depend on the ability to measure the tidal deformation of inspiraling neutron stars.

Conclusion
In this work we have explored the effects of the neutron-star equation-of-state in binary neutron star observations with a particular focus on spinning binary neutron stars and the spin-quadrupole term that is often ignored, as it was in [3]. We have explored the overall distinguishability of waveforms as a function of equation-of-state related terms and found that the spin-quadrupole term has a much larger effect than the tidal-deformability for highly spinning neutron star systems, although both terms are potentially measurable if all other parameters about the system are known. We have explored whether the equation-of-state would have any effect on our ability to observe binary neutron-star mergers, where it is commonly assumed that the compact objects are both black holes. We found that the tidal deformability can lead to a small reduction in the number of observed binary neutron-star systems, as reported in [31], but found that the spin-quadrupole term is largely degenerate with the component spins and mass ratio. Therefore, ignoring this will not result in a reduction in the detection rate other than at the boundaries of the searched parameter space. We have explored the bias in recovered source parameters that can be expected if making incorrect assumptions about the neutron-star equation-of-state and have found that the recovered values of the component spins and mass ratio can be significantly biased. We have also explored the measurability of various equation-of-state terms, finding that the spin-quadrupole term alone cannot be easily measured with Advanced LIGO, but that combined with the effects of the tidal deformability it will be possible to rule out specific equations-of-state, especially in the case that the neutron stars' component spins are large.
As always with investigations of this type, our results are only as good as the waveform model used. The Post-Newtonian waveform model we used has two caveats, first that there may be uncontrolled systematic errors in the phasing model because the Post-Newtonian approximation breaks down in the relativistic regime close to the merger, and second because the waveform stops abruptly at a given frequency, whereas a real binary neutron star source might be expected to have a complex post-merger structure. However, for the case of exactly equal mass systems the TaylorF2 inspiral waveforms give similar results as more sophisticated waveform models [89], and the post-merger signals are expected to be at much larger frequencies than those that Advanced LIGO can observe. Thus, we expect our broad conclusions to remain unchanged in more realistic contexts. While work to improve waveform models in the binary-neutron star regime will be important to remove the potential for systematic errors in gravitational-wave observations, we believe that the results and conclusions that we derive in this manuscript will remain largely unaltered.
With many additional binary neutron-star mergers being expected to be observed by Advanced LIGO and Advanced Virgo in the coming years it will be interesting to apply these methods also on real data and try to measure directly the equation-of-state of neutron stars, as already explored in [3]. However, as such works are allowing for the possibility of very high spin neutron-star systems, neglecting the spin-quadrupole terms will lead to biased results, so we emphasize the importance of including this term both in waveform modelling and in analysis runs performed on real data in the future. The conclusions of this work are also very relevant for proposed third-generation gravitational-wave observatories, such as the Einstein Telescope [90] and Cosmic Explorer [91]. With these instruments the rate of binary neutron star observations will increase by ∼1000 over that expected with second-generation observatories, and allowing occasional observations with signal-to-noise ratio as large as ∼1000.
With such large-SNR observations, ignoring the effect of the rotational deformation of the neutron stars will lead to parameter biases that are much larger than the statistical uncertainty. However, these instruments offer the ability to measure the equation-of-state of neutron-stars with considerably greater accuracy. We leave the problem of extending the results of this manuscript to third-generation observatories to future work.