Equation-of-state insensitive relations after GW170817

TL;DR

This paper improves the robustness of neutron star property relations derived from gravitational wave data by reducing equation-of-state variability, thereby decreasing systematic uncertainties in future measurements.

Contribution

It introduces a method to reduce equation-of-state variability in neutron star relations by applying priors from GW170817 data, enhancing the reliability of gravitational wave inferences.

Findings

Reduced equation-of-state variability by over 50%.

Improved insensitivity of neutron star relations to the equation of state.

Systematic uncertainties in tidal deformability extraction are decreased.

Abstract

The thermodynamic relation between pressure and density (i.e. the equation of state) of cold supranuclear matter is critical in describing neutron stars, yet it remains one of the largest uncertainties in nuclear physics. The extraction of tidal deformabilities from the gravitational waves emitted in the coalescence of neutron star binaries, such as GW170817, is a promising tool to probe this thermodynamic relation. Equation-of-state insensitive relations between symmetric and antisymmetric combinations of individual tidal deformabilities, the so-called "binary Love relations", have proven important to infer the radius of neutron stars, and thus constrain the equation of state, from such gravitational waves. A similar set of relations between the moment of inertia, the tidal deformability, the quadrupole moment, and the compactness of neutron stars, the so-called `I-Love-Q' and `C-Love' relations, allow for future tests of General Relativity in the extreme gravity regime. But even the most insensitive of such relations still presents some degree of equation-of-state variability that could introduce systematic uncertainties in parameter extraction and in model selection. We here reduce this variability by more than 50% by imposing a prior on the allowed set of equations of state, derived from the posteriors generated from the analysis of GW170817. The resulting increase in insensitivity reduces systematic uncertainties in the extraction of the tidal deformability from future gravitational wave observations, although statistical uncertainties currently dominate the error budget, and will continue to do so until the era of Voyager-class detectors.

Figures

Tables

Table 1. Table 1: Fit parameters for the I-Love-Q and C-Love relations using the constrained set and the fitting functions in Eq. ( 6 ) (top) and in Eq. ( 8 ) (bottom).

$y$	$x$	$\alpha$	$K_{yx}$	$a_1$	$a_2$	$a_3$	$a_4$	$a_5$	$b_1$	$b_2$	$b_3$
$\overline{I}$	$\Lambda$	(–)	(–)	$1.493$	$0.06410$	$0.02085$	$-5.018\times {10}^{-4}$	$3.16\times {10}^{-7}$	(–)	(–)	(–)
$\overline{Q}$	$\Lambda$	(–)	(–)	$0.2093$	$0.07404$	$0.05382$	$-5.018\times {10}^{-3}$	$1.576\times {10}^{-4}$	(–)	(–)	(–)
$\overline{I}$	$\overline{Q}$	(–)	(–)	$1.383$	$0.5931$	$-0.02161$	$0.04190$	$-2.968\times {10}^{-3}$	(–)	(–)	(–)
$\overline{I}$	$\Lambda$	$2/5$	$0.5313$	$1.287$	$0.09888$	$-2.300$	(–)	(–)	$-1.347$	$0.3857$	$-0.02870$
$\overline{Q}$	$\Lambda$	$1/5$	$3.555$	$-2.122$	$2.72$	$-1.491$	(–)	(–)	$0.8644$	$-0.1428$	$-1.397$
$\overline{I}$	$\overline{Q}$	$2$	$0.008921$	$10.59$	$-37.46$	$43.18$	(–)	(–)	$-2.361$	$1.967$	$-0.5678$
$C$	$\Lambda$	$-1/5$	$0.2496$	$-919.6$	$330.3$	$-857.2$	(–)	(–)	$-383.5$	$192.5$	$-811.1$

Table 2. Table 2: Maximum relative and fractional EoS variation in the I-Love-Q, C-Love, R-Love, and binary Love relations using the unconstrained set, the constrained set and variations reported in previous work  Yagi and Yunes ( 2013b , 2016a ) . The maximum absolute EoS variation is also reported in the C-Love, R-Love and binary Love cases in parentheses. The maximum variation in the constrained set case is better than a factor of two smaller than the variability in the unconstrained case and in previous work. The maximum EoS variation in the unconstrained set is slightly larger than that found in previous work because the former is built from a large random sampling of EoSs.

EoS-insensitive	Maximal EoS Variability
Relation	Previous	Unconstrained	Constrained
$\overline{I}-\Lambda$	$0.0059$	$0.0077$	$0.0031$
$\overline{Q}-\Lambda$	$0.010$	$0.013$	$0.0047$
$\overline{I}-\overline{Q}$	$0.012$	$0.015$	$0.0057$
$C-\Lambda$	$0.065$	$0.072$	$0.022$
	(–)	($0.018$)	($0.0066$)
$R-\Lambda$	–	$0.056$	$0.022$
	(–)	($880\text{ m}$)	($360\text{ m}$)
${\Lambda }_a-{\Lambda }_s$	$\sim 0.50$	$0.57$	$0.21$
$q=0.90$	(–)	($190$)	($37$)
${\Lambda }_a-{\Lambda }_s$	$\sim 0.20$	$0.25$	$0.083$
$q=0.75$	(–)	($320$)	($52$)
${\Lambda }_a-{\Lambda }_s$	$\sim 0.025$	$0.038$	$0.018$
$q=0.50$	(–)	($240$)	($29$)

Table 3. Table 3: Maximum relative and fractional EoS variation in the I-Love relation, fitting to three different sets of data: using the constrained set plus hybrid EoSs, using the constrained set alone, and using only the hybrid EoSs. In all 3 cases the hybrid EoSs are EoS-insensitive to ∼ 1 similar-to absent 1 \sim 1 %, which is a slight decrease in universality relative to the hadronic only EoSs.

Fitting	Maximal EoS Variability
Case	Constrained	Hybrid
Combined	$0.0044$	$0.014$
(Case 1)
Constrained only	$0.0031$	$0.017$
(Case 2)
Hybrid only	$0.0084$	$0.010$
(Case 3)

Table 4. Table 4: Similar to Table  3 but for the C-Love relation. For all 3 cases, the hybrid EoSs are only universal up to a minimum of ∼ 1 similar-to absent 1 \sim 1 % (fractional difference from the fit), and the constrained EoSs typically outperform the hybrid ones (other than the third case where the hybrid stars show slightly better agreement to the fit). The second case is also shown in Fig.  4 .

Fitting	Maximal EoS Variability
Case	Constrained	Hybrid
Combined	$0.037$	$0.055$
(Case 1)
Constrained only	$0.022$	$0.072$
(Case 2)
Hybrid only	$0.058$	$0.045$
(Case 3)

Table 5. Table 5: Coefficients to the fits given by Eqs. ( 19 )-( 22 ) for the relative error on Λ a subscript Λ 𝑎 \Lambda_{a} in the improved binary Love EoS-insensitive relations presented in this paper.

${\mu }_1$	$3.509\times {10}^{-3}$	${\sigma }_1$	$-2.074\times {10}^{-7}$
${\mu }_2$	$9.351\times {10}^{-1}$	${\sigma }_2$	$1.492\times {10}^{-3}$
${\mu }_3$	$-18.07$	${\sigma }_3$	$-4.891\times {10}^{-2}$
${\mu }_4$	$27.56$	${\sigma }_4$	$8.207\times {10}^{-1}$
${\mu }_5$	$-10.10$	${\sigma }_5$	$-1.308$
		${\sigma }_6$	$-63.76$
		${\sigma }_7$	$11.14$
		${\sigma }_8$	$75.25$
		${\sigma }_8$	$-23.69$

Table 6. Table 6: Approximate signal-to-noise ratio ρ GW170817 A subscript superscript 𝜌 𝐴 GW170817 \rho^{A}_{\text{GW170817}} and ( 1 ​ σ 1 𝜎 1\sigma ) statistical uncertainty on the extraction of λ 0 subscript 𝜆 0 \lambda_{0} had a single event like GW170817 been observed by future interferometer A 𝐴 A and had interferometer A 𝐴 A observed N A subscript 𝑁 𝐴 N_{A} events in a 1 year observation, using aLIGO, A + , Voyager, CE, and ET. The number of events N A subscript 𝑁 𝐴 N_{A} , and the combined statistical uncertainty depends on the binary NS merger detection rate, and thus we include results assuming an upper, a central and a lower limit on this rate. The statistical uncertainties on λ 0 subscript 𝜆 0 \lambda_{0} becomes comparable with the systematic uncertainty (set to be P 90 = 13.19 subscript 𝑃 90 13.19 P_{90}=13.19 ) from using the improved binary Love relations with detectors of Voyager-class or better.

Detectors (A)	GW170817		Multiple events
	${\rho }_{\text{GW170817}}^A$	${\sigma }_{\text{GW170817}}^A$	$N_A$			${\sigma }_N^A$
			Low	Central	High	Low	Central	High
O2	$3.2\times {10}^1$	$1.7\times {10}^2$	–	–	–	–	–	–
aLIGO	$9.1\times {10}^1$	$1.1\times {10}^2$	$2.0\times {10}^1$	$9.8\times {10}^1$	$3.0\times {10}^2$	$1.8\times {10}^2$	$8.3\times {10}^1$	$4.7\times {10}^1$
A+	$1.8\times {10}^2$	$4.6\times {10}^1$	$1.6\times {10}^2$	$7.9\times {10}^2$	$2.4\times {10}^3$	$5.9\times {10}^1$	$2.5\times {10}^1$	$1.4\times {10}^1$
Voyager	$4.3\times {10}^2$	$2.5\times {10}^1$	$2.2\times {10}^3$	$1.1\times {10}^4$	$3.2\times {10}^4$	$2.1\times {10}^1$	$9.6\times {10}^0$	$5.3\times {10}^0$
ET-D	$1.4\times {10}^3$	$6.9\times {10}^0$	$7.2\times {10}^4$	$3.4\times {10}^5$	$1.1\times {10}^6$	$3.8\times {10}^0$	$1.7\times {10}^0$	$9.6\times {10}^{-1}$
CE	$2.8\times {10}^3$	$7.7\times {10}^0$	$3.0\times {10}^5$	$1.4\times {10}^6$	$4.4\times {10}^6$	$3.7\times {10}^0$	$1.7\times {10}^0$	$9.0\times {10}^{-1}$

Equations

$$\Gamma (x)=\exp \sum _k^N{\gamma }_k{x}^k,$$

$$\frac{dϵ(p)}{dp}=\frac{ϵ(p)+p}{p\Gamma (p)}.$$

$$Q_{ij}=-\lambda {\epsilon }_{ij}.$$

$$\Phi =-\frac{1+g_{tt}}{2}$$

$$\displaystyle+\frac{1}{2}\varepsilon_{ij}x^{i}x^{j}+\mathcal{O}\Bigg{(}\frac{r^{3}}{R^{3}}\Bigg{)},$$

$${\Lambda }_s\equiv \frac{{\Lambda }_2+{\Lambda }_1}{2},{\Lambda }_a\equiv \frac{{\Lambda }_2-{\Lambda }_1}{2},$$

$$\ln y=a_1+a_2\ln x+a_3(\ln x{)}^2+a_4(\ln x{)}^3+a_5(\ln x{)}^4,$$

$${\bar{I}}^{\text{N}}=K_{\bar{I}\Lambda }{\Lambda }^{2/5},{\bar{Q}}^{\text{N}}=K_{\bar{Q}\Lambda }{\Lambda }^{1/5},{\bar{I}}^{\text{N}}=K_{\bar{I}\bar{Q}}{\bar{Q}}^2,$$

$$y=K_{yx}{x}^{\alpha }\frac{1+{\sum }_{i=1}^3{a}_i{x}^{-i/5}}{1+{\sum }_{i=1}^3{b}_i{x}^{-i/5}},$$

$$C=\sum _{k=0}^2{a}_k(\ln \Lambda {)}^k.$$

$$C^N=K_{C\Lambda }{\Lambda }^{-1/5}.$$

$${\Lambda }_a=F_n(q)\frac{1+{\sum }_{i=1}^3{\sum }_{j=1}^2{b}_{ij}{q}^j{\Lambda }_s^{-i/5}}{1+{\sum }_{i=1}^3{\sum }_{j=1}^2{c}_{ij}{q}^j{\Lambda }_s^{-i/5}}{\Lambda }_s^{\alpha },$$

$$F_n(q)\equiv \frac{1-q^{10/(3-n)}}{1+q^{10/(3-n)}}.$$

$$b_{ij}=\left[\begin{array}{cc}-14.40& 14.45\\ 31.36& -32.25\\ -22.44& 20.35\end{array}\right]c_{ij}=\left[\begin{array}{cc}-15.25& 15.37\\ 37.33& -43.20\\ -29.93& 35.18\end{array}\right],$$

$$\Lambda ={\lambda }_0+{\lambda }_1\left(1-\frac{M}{M_0}\right)+\mathcal{O}\left[{\left(1-\frac{M_0}{M}\right)}^2\right].$$

$$r({\Lambda }_s,q)\equiv {\Lambda }_a^{\text{relation}}({\Lambda }_s,q)-{\Lambda }_a^{\text{true}},$$

$${\Lambda }_a^{\text{new relation}}={\Lambda }_a^{\text{relation}}({\Lambda }_s,q)+\mathcal{N}({\mu }_r({\Lambda }_s,q),{\sigma }_r({\Lambda }_s,q)),$$

$${\mu }_r({\Lambda }_s,q)$$

$${\sigma }_r({\Lambda }_s,q)$$

$${\mu }_r({\Lambda }_s)$$

$${\mu }_r(q)$$

$${\sigma }_r({\Lambda }_s)$$

$${\sigma }_r(q)$$

$$\mathcal{R}({\Lambda }_s,q)={\lambda }_0^{\text{relation}}({\Lambda }_s,q)-{\lambda }_0^{\text{true}}.$$

$${\theta }^a=(\ln A,{\varphi }_c,t_c,\ln \mathcal{M},\ln \eta ,{\chi }_s,{\chi }_a,{\lambda }_0,{\lambda }_1),$$

$$\Delta\theta^{a}=\sqrt{\Big{(}\tilde{\Gamma}^{-1}\Big{)}^{aa}},$$

$$\tilde{\Gamma}_{ab}\equiv\Big{(}\frac{\partial h}{\partial\theta^{b}}\Big{|}\frac{\partial h}{\partial\theta^{a}}\Big{)}+\frac{1}{\sigma_{\theta^{a}}^{2}}\delta_{ab}\,.$$

$$(a|b)\equiv 2{\int }_0^{\infty }\frac{{\tilde{a}}^{\ast }\tilde{b}+{\tilde{b}}^{\ast }\tilde{a}}{S_n(f)}df,$$

$$f(\rho )=3{\rho }_{\text{th}}^3/{\rho }^4$$

$$N_A=\Delta {\tau }_0\underset{0}{\overset{z_h}{\int }}4\pi [a_0{r}_1(z){]}^2\mathcal{R}r(z)\frac{d\tau }{dz}dz.$$

$$a_0{r}_1(z)=\frac{1}{H_0}\underset{0}{\overset{z}{\int }}\frac{d{z}^{\prime }}{\sqrt{(1-{\Omega }_{\Lambda })(1+z^{\prime }{)}^3+{\Omega }_{\Lambda }}},$$

$$\frac{d\tau }{dz}=\frac{1}{H_0}\frac{1}{1+z}\frac{1}{\sqrt{(1-{\Omega }_{\Lambda })(1+z^{\prime }{)}^3+{\Omega }_{\Lambda }}},$$

$$r(z)=\left\{\begin{array}[]{ll}1+2z&z\leq 1\\ \frac{3}{4}(5-z)&1\leq z\leq 5\\ 0&z\geq 5\\ \end{array}\,,\right.$$

Full text

I Introduction

The thermodynamic relation between pressure and density in cold, supranuclear matter, the so-called equation of state (EoS), remains an unsolved problem in nuclear physics. The EoS is critical to our understanding of neutron stars (NSs) because it determines many NS observables, such as their mass, radius, moment of inertia ($I$), quadrupole moment ($Q$) and tidal deformability (or Love number). Unfortunately, terrestrial experiments can only probe the EoS to around nuclear saturation density ($\rho_{\text{sat}}\approx 2.5\times 10^{14}\text{ g/cm}^{3}$) Li et al. (2008); Tsang et al. (2009); Centelles et al. (2009); Li and Chen (2005); Chen et al. (2005). Although some temperature-dependent heavy-ion collision experiments can probe higher densities Danielewicz et al. (2002), astrophysical observations of NSs remain ideal for constraining the EoS of cold and ultra dense, nuclear matter.

Independent measurements of NS observables can be used to constrain the nuclear EoS. For example, electromagnetic observations of the mass and radius of certain NSs have been used to place confidence limits in the mass-radius plane, and thus constrain the EoS Guver and Ozel (2013); Ozel et al. (2010); Steiner et al. (2010); Lattimer and Steiner (2014); Ozel and Freire (2016). These observations, however, may potentially suffer from large systematic errors Miller and Lamb (2016); Miller due to uncertainties in the astrophysical modeling of X-ray bursts. The gravitational waves (GWs) emitted in the coalescence of NS binaries may be a cleaner probe of nuclear physics. During the early inspiral, the orbital separation is large enough that the tidal fields are negligible; but as the orbital separation decreases due to GW emission, tidal forces grow and the NSs respond by developing deformations determined by their nuclear EoS. These deformations source additional multipole radiation as well as affect the orbital trajectory of the binary, thus altering the GWs emitted, encoding within the latter the NS EoS Hinderer (2008); Flanagan and Hinderer (2008).

The GWs emitted by binary NSs in the late inspiral must then depend on the tidal deformabilities $\Lambda_{1}$ and $\Lambda_{2}$111Throughout this work we use the convention that binary quantities are indexed as $1,2$ where 1 always correspond to the more massive object., which control the linear response of the star’s quadrupole deformation to the (electric-type, quadrupole) tidal field of the companion (to leading order in a post-Newtonian expansion Blanchet (2014)) Flanagan and Hinderer (2008); Vines et al. (2011). These parameters, however, enter the GW waveform model at the same post-Newtonian order, making them degenerate, and thus, very difficult to estimate independently with current GW data Wade et al. (2014a). Instead, one can extract certain combinations of the tidal deformabilities, such as a certain mass-weighted tidal deformability $\tilde{\Lambda}$ Favata (2014); Wade et al. (2014a), or one can extract the (mass-independent) coefficients $(\lambda_{0},\lambda_{1},\ldots)$ of a Taylor expansion of the tidal deformabilities about some fiducial mass Messenger and Read (2012); Del Pozzo et al. (2013); Yagi and Yunes (2016a). Current detectors are not sensitive enough to accurately measure any of these coefficients, but future detectors will, and the information from multiple events could then be combined, since the Taylor expansion coefficients should be common to all events.

Lacking enough sensitivity in current GW observations, one is forced to only estimate the mass-weighted tidal deformability, but this prevents the independent extraction of $\Lambda_{1}$ and $\Lambda_{2}$. Yagi and Yunes Yagi and Yunes (2016b, a) proposed a solution to this problem by finding “approximately universal” or “EoS-insensitive” relations between the symmetric and anti-symmetric combinations of the tidal deformabilities $\Lambda_{s,a}=\frac{1}{2}(\Lambda_{1}\pm\Lambda_{2})$, the so-called “binary Love relations.” These relations can be used to analytically express $\Lambda_{s}$ in terms of $\Lambda_{a}$ (or vice-versa), making the mass-weighted tidal deformability a function of only $\Lambda_{a}$. A measurement of the mass-weighted tidal deformability then implies a measurement of $\Lambda_{a}$, and through the use of the binary Love relations, also a measurement of $\Lambda_{s}$, which then allows for the inference of the individual tidal deformabilities $\Lambda_{1}$ and $\Lambda_{2}$ Yagi and Yunes (2016b, a). With those at hand, one can further use EoS-insensitive relations between the tidal deformabilities and the compactness, the so-called “C-Love relations” Yagi and Yunes (2013a, b); Maselli et al. (2013); Yagi and Yunes (2017), to infer the radii of the NSs, and thus, to place two constraints in the mass-radius plane, one for each star in the binary. This idea was recently implemented for GW170817 Abbott et al. (2017a), allowing EoS-independent constraints on the mass-radius curve using GW data Chatziioannou et al. (2018); Abbott et al. (2018a).

EoS-insensitive relations can in fact be used for more than just measuring the nuclear EoS. For years, the theoretical physics community considered the possibility of using measurements of NS properties, such as the mass, the radius and the moment of inertia, to constrain deviations from General Relativity in the strong-field regime. Certain modified theories of gravity, such as scalar tensor theories with spontaneous scalarization Damour and Esposito-Farese (1996), Einstein-Æther and Hořava gravity Eling et al. (2007); Yagi et al. (2014a, b), dynamical Chern-Simons gravity Yunes et al. (2010); Yagi et al. (2013); Gupta et al. (2018), beyond Horndesky theories Babichev et al. (2016); Sakstein et al. (2017), modify such NS observables, but unfortunately, these modifications are typically degenerate with the nuclear EoS Pani and Berti (2014); Minamitsuji and Silva (2016); Maselli et al. (2016). Yagi and Yunes proposed to solve this problem by finding EoS-insensitive relations between the moment of inertia, the tidal deformability (or Love number) and the quadrupole moment, the so-called “I-Love-Q” relations Yagi and Yunes (2013a, b). Given a measurement of the Love number for a given NS, for example through GW observations, the I-Love-Q relations can be used to infer the moment of inertia or the quadrupole moment. A second independent measurement of either of these two quantities, for example through binary pulsar observations Lattimer and Schutz (2005) or observations with the Neutron star Interior Composition ExploreR (NICER) Ozel et al. (2016), then allows an EoS-insensitive test of General Relativity in the strong field regime Yagi and Yunes (2013a, b); Gupta et al. (2018); Doneva and Pappas (2018).

The implementation of the EoS-insensitive relations in data analysis has to somehow contend with the fact that these relations are in fact not exactly universal, but rather present different (albeit small) levels of EoS variability. In the case of the binary Love and C-Love relations to infer the radii of NSs with GW170817, the problem is solved by marginalizing over the EoS variability Chatziioannou et al. (2018). Presumably, this same procedure can be applied in the future when carrying out tests of General Relativity with the I-Love-Q relations, using a combination of binary pulsar, NICER and GW data. But this marginalization procedure may not always be as important; as constraints in the mass-radius plane become more stringent with future GW observations, the allowed space of EoSs will shrink, which in turn must naturally decrease the degree of EoS variability in all EoS-insensitive relations. This is the main focus of this paper.

I.1 Executive Summary

We begin by investigating the increase in EoS insensitivity due to the constraints placed by GW170817 on the allowed space of EoSs. We use the posterior probability distribution on the pressure-density plane Carney et al. (2018); Abbott et al. (2018a) obtained from GW170817 Abbott et al. (2017a) to inform the set of EoSs that is compatible with this observation. We then generate two large samples of spectral EoSs Lindblom (2018), one in which the EoSs are directly sampled from the posterior probability distribution (the “constrained” sample) and another in which this constraint not enforced (the “unconstrained” sample). We repeat the analysis done by Yagi and Yunes Yagi and Yunes (2016a, 2013b) on both sets of EoSs and find that the EoS-insensitive relations present less EoS variability with the constrained set. In particular, the EoS-insensitivity increases by a factor of $\sim 60$% in the binary Love relations (for stars with mass ratio larger than $0.75$), by a factor of $\sim 70\%$ in the C-Love relations, and by factors of $\sim 50$% in the I-Love-Q relations (see Table 2 and Fig. 3 for more details).

With this study at hand, we then carry out additional related studies on EoS-insensitive relations that go beyond the work in Yagi and Yunes (2016a, 2013b). First, we investigate the relation between the NS radius and its tidal deformability, the R-Love relations, for both sets of EoSs, as these are critical in order to place constraints on the radius from a measurement of the Love number. We use the C-Love relations to construct the R-Love relations and find that the maximum EoS variability drops from $\sim 880$ m in the unconstrained case to $\sim 360$ m in the constrained case (see Sec. III.2 for more details). Second, we study the EoS-universality of hybrid stars, which experience strong first-order phase transitions from hadronic to quark matter in the core Paschalidis et al. (2018); Most et al. (2019); Burgio et al. (2018); Montana et al. (2019). We find that the I-Love-Q and C-Love relations remain EoS-insensitive for these hybrid stars, although the EoS variability increases slightly (see Sec. III.1.2 and III.2.2 for more details). However, we also find that the binary Love relations are not EoS-insensitive for a mixed binary with a (massive) hybrid star and a (low-mass) hadronic star, due to the large separation in mass-weighted tidal deformability between the constituent stars (see Sec. III.3 for more details).

Last but not least, we study the importance of using the improved EoS-insensitive relations in future GW observations. The use of EoS-insensitive relations introduces systematic uncertainties in parameter estimation because of the intrinsic non-zero EoS-variability in these relations, which one must marginalize over. These uncertainties are currently irrelevant because statistical uncertainties in parameter estimation are much larger with current detectors. But as the detector sensitivity is improved, the signal-to-noise ratio and the number of events that will be detected will increase, therefore decreasing the statistical uncertainties below the systematic one due to EoS variability. We carry out Fisher analyses to estimate when the statistical uncertainties become comparable to the systematic uncertainties due to EoS variability and find that this occurs for Voyager-class detectors.

Figure 1 shows this result in more detail. We present the (Fisher-estimated) statistical uncertainty in the measurement of $\lambda_{0}$ with various detectors (LIGO O2 aLI , Advanced LIGO at design sensitivity (aLIGO) aLI , LIGO A+ (A+)Ap_ , Voyager Ap_ , the Einstein Telescope (ET) ET , and the Cosmic Explorer (CE) Ap_ ) for an event similar to GW170817. The x-axis shows the signal-to-noise ratio for a GW170817 event detected with each of these detectors. We also present the combined statistical uncertainty $\sigma^{A}_{N}$ after $N$ detections with each of these instruments, with the top and bottom of the region representing the most optimistic and pessimistic expectation for the number of detections expected from the binary NS merger rate Abbott et al. (2017b) (see also Table 6). These statistical uncertainties should be compared to the systematic uncertainty in $\lambda_{0}$ due to EoS-variability improved with the constrained set. We find that the statistical and the systematic uncertainties cross for Voyager-class detectors.

The remainder of this paper presents the details of the results summarized above and it is organized as follows. We begin with complementary background and theory material in Sec. II. We continue in Sec. III by finding new and improved binary Love, I-Love-Q, and C-Love relations, and considering how well hybrid star EoSs agree with the EoS-insensitive relations. We next examine these improved relations and question whether or not they are useful for future interferometers in Sec. IV. We conclude in Sec. V by discussing our results and mentioning avenues of future work. Throughout this paper, we have adopted geometric units of $G=1=c$, unless otherwise stated.

II Background and theory

In this section we review how the EoS can be represented analytically through a spectral decomposition, and how the observation of GW170817 constrains the space of possible EoSs. We then proceed to discuss how one computes the tidal deformabilities and the EoS-insensitive relations.

II.1 Spectral representations of NS EoSs

The structure of a NS and its tidal interactions in a binary system rely heavily on the underlying state function (or equation of state - EoS) describing the relationship between the pressure ($p$) and energy density ($\epsilon$) of nuclear matter. Given that all currently proposed EoSs utilize certain approximations Oertel et al. (2017); Baym et al. (2018), one method to study a wide range of physically realizable EoSs is to parameterize them such that any realistic EoS can be represented with a small number of parameters. Spectral representations Lindblom (2010); Lindblom and Indik (2012a, 2014); Lindblom (2018); Abbott et al. (2018b) parameterize EoSs by performing spectral expansions on the adiabatic index $\Gamma(p)$222Another way of parameterizing EoSs is through a piecewise polytropic formulation Read et al. (2009); Lackey and Wade (2015); Carney et al. (2018).:

[TABLE]

where $x\equiv\log{(p/p_{0})}$ for a minimum pressure $p_{0}$. The equation of state is then determined by an integration of the differential equation:

[TABLE]

Using this formalism, any valid EoS can be approximated through the choice of $N$ spectral coefficients $\gamma_{k}$, and we here choose $N=4$, tabulated for several common EoSs in Table 1 of Lindblom (2018).

We here wish to consider EoS-insensitive relations using two sets of EoSs: a “constrained set” consistent with the observation of GW170817 and an “unconstrained set” that does not impose this prior. In both cases, we model the EoS with a $C^{0}$-piecewise function that equals the low-density crust EoS of SLy333SLy represents the unified equation of state based on the effective Skyrme Lyon nucleon-nucleon interactions, developed by Douchin and Haensel, 2001 Douchin and Haensel (2001). Douchin and Haensel (2001) below half nuclear saturation density $\rho_{\text{stitch}}=1.3\times 10^{14}\text{ g/cm}^{3}$ Read et al. (2009), and equals the spectral decomposition described above outside the crust. For the latter, we restrict the spectral coefficients to the ranges $\gamma_{0}\in[0.2,2]$, $\gamma_{1}\in[-1.6,1.7]$, $\gamma_{2}\in[-0.6,0.6]$, $\gamma_{3}\in[-0.02,0.02]$, and the adiabatic index is further restricted to $\Gamma\in[0.6,4.5]$ Lindblom and Indik (2012b). Moreover, we impose the following two restrictions: (i) causality within 10%, i.e. that the speed of sound of the fluid be less than the speed of light to 10% (following the analysis of Ref. Abbott et al. (2018a)), and (ii) a high maximum mass, i.e. that the resulting EoS supports NSs with masses at least as high as $1.97\text{ M}_{\odot}$, consistent with astrophysical observations Demorest et al. (2010); Antoniadis et al. (2013); Zhao (2015).

The unconstrained set is then defined by drawing random samples in the spectral coefficients within their allowed prior ranges, and then eliminating any EoS that either leads to an adiabatic index $\Gamma$ outside the allowed range, breaks the causality restriction, or breaks the maximum mass restriction. To obtain the constrained set of EoSs we analyze the publicly available data for GW170817 from the Gravitational Wave Open Science Center Gravitational Wave Open Science Center () (GWOSC); Vallisneri et al. (2015). We use the publicly available software library LALInference Veitch et al. (2015); LIGO Scientific Collaboration and Virgo Collaboration (2017) to sample the EoS posterior, similar to Ref Abbott et al. (2018a). Our analysis uses the same settings and prior choice as those of Abbott et al. (2018a); Carney et al. (2018). In both cases, each set consists of 100 members, with a subset of these shown in the left panel of Fig. 2.

In addition, we investigate 10 transitional quark-hadron matter stars, which undergo strong first-order phase transitions at a pressure $P_{\text{tr}}$, leading to the hadronic branch departing into a quark-matter branch at a given transitional mass Paschalidis et al. (2018); Alford and Sedrakian (2017); Seidov (1971); Zdunik and Haensel (2013); Alford et al. (2013). In particular, we focus on the ACS and ACB models (corresponding to speed of sound, and piecewise polytropic representations of the EoS beyond the phase transition, respectively) described in Paschalidis et al. (2018), and also shown in the right panel of Fig. 2. These result in two distinct types of NSs, based on their mass: (i) massive ($M\geq\text{ M}_{\text{tr}}$) hybrid stars which have quark-matter inner cores and nuclear matter elsewhere (henceforth, we denote such stars as hybrid stars (HSs)), and (ii) low-mass ($M\leq\text{ M}_{\text{tr}}$) hadronic stars with no internal transition to quark matter (henceforth, we denote these stars as simply NSs).

II.2 NS Tidal deformability

The GWs emitted in the coalescence of binary NSs, such as GW170817, provide valuable insight into their internal structure. As the stars spiral into each other, they tidally deform due to their companion’s tidal field. The tidal field can be decomposed into multipole moments, with the leading-order field being even parity and quadrupolar, and thus leading to an even parity and quadrupolar deformation in response. The constant of proportionality that determines how quadrupolarly deformed a star becomes in the presence of a quadrupolar tidal field is called the tidal deformability and denoted $\lambda$. The magnitude of $\lambda$ depends strongly on the underlying structure, as well as the fluid nature of the NS Char and Datta (2018). Therefore, this parameter has the largest impact on the GW phase, and thus, it is encoded in GW observations of binary NS coalescence.

Consider a NS of mass $M$ in the presence of the (spatial) tidal tensor $\varepsilon_{ij}$ of its companion. In response to this tidal field, the NS will deform away from sphericity and acquire a (spatial) quadrupole moment tensor $Q_{ij}$. The tidal deformability $\lambda$ quantifies this linear response via Flanagan and Hinderer (2008); Hinderer (2008); Yagi and Yunes (2013b):

[TABLE]

One can extract $Q_{ij}$ and $\varepsilon_{ij}$ – and thus $\lambda$, or its dimensionless form $\Lambda\equiv\lambda/M^{5}$ – via a double asymptotic expansion of the gravitational potential in a buffer zone, defined by ${\cal{L}}\gg r\gg R$:

[TABLE]

where $R$ is the radius of the star and ${\cal{L}}$ is the length scale associated with the radius of curvature induced by the companion.

The calculation of the tidal deformability therefore requires the calculation of the metric component $g_{tt}$. One typically starts by first constructing a spherically symmetric, non-spinning background solution, whose stellar radius and mass are determined from $p(R)=0$ and $M=(1-g_{rr}(R)^{-1})R/2$. One then perturbatively introduces a tidal deformation and solves the perturbed Einstein equations in the NS interior, and matches this solution to an exterior solution at the surface, modulo two integration constants. The requirement that the metric be at least $C^{1}$ fixes the ratio between the two constants of integration, yielding the tidal deformability, as discussed in more detail in Hinderer (2008).

Similar to tidal perturbations, one can consider rotational perturbations. At first order in rotation, the stellar moment of inertia can be extracted by looking at the asymptotic behavior of the $g_{t\phi}$ metric component. At second order in rotation, the quadrupole moment can be extracted by studying the asymptotic behavior of the $g_{tt}$ metric component.

Consider now a binary NS system, as that which produced the GW170817 event, where each star individually experiences the tidal field of the other star. Thus, each star possesses tidal deformabilities $\Lambda_{1}$ and $\Lambda_{2}$ that enter the GW phase and encode information about each star’s EoS. Due to strong correlations between these parameters, individual extraction is very difficult with current interferometer sensitivity limitations. However, the waveform templates can be strategically reparameterized to instead include independent linear combinations of $\Lambda_{1}$ and $\Lambda_{2}$ that mitigate correlations. One common reparameterization is through the introduction of the mass-weighted tidal deformability $\tilde{\Lambda}=\tilde{\Lambda}(\Lambda_{1},\Lambda_{2})$ and the parameter $\delta\tilde{\Lambda}=\delta\tilde{\Lambda}(\Lambda_{1},\Lambda_{2})$ Favata (2014); Wade et al. (2014a). Since these parameters enter the GW phase at 5 and 6PN orders respectively444A term of $N$PN order is proportional to $v^{2N}$ relative to the leading-order term in the expression, where $v$ represents the relative velocity of binary constituents., they partially break the degeneracies between $\Lambda_{1}$ and $\Lambda_{2}$.

II.3 EoS-insensitive relations

Current GW interferometry is not yet sensitive enough to accurately extract both tidal parameters $\tilde{\Lambda}$ and $\delta\tilde{\Lambda}$. In a search to remedy this, Yagi and Yunes Yagi and Yunes (2016a) found that symmetric and antisymmetric combinations of the tidal deformabilities

[TABLE]

display EoS-insensitive properties to a high degree, showing EoS variations of at most 20% for binaries with masses less than $1.7\text{ M}_{\odot}$ and using a representative sample of 11 EoSs. These “binary Love relations” allow one to analytically break degeneracies between the tidal parameters: one can substitute $\Lambda_{a}=\Lambda_{a}(\Lambda_{s})$ in the GW model, thus completely eliminating $\Lambda_{a}$ from the parameter list (or vice-versa). This is important for two reasons: (i) the new model allows for the more accurate extraction of $\Lambda_{s}$ (or $\Lambda_{a}$ if $\Lambda_{s}$ is eliminated), and (ii) the relations allow for the inference of $\Lambda_{a}$ given a measurement of $\Lambda_{s}$ (or vice-versa), and from this for the inference of the individual tidal deformabilities $\Lambda_{1}$ and $\Lambda_{2}$. A simple Fisher analysis has shown that the binary Love relations improve parameter estimation of $\tilde{\Lambda}$ by up to an order of magnitude Yagi and Yunes (2016b, a).

Similar EoS-insensitive relations have been found between individual NS observables: the moment of inertia ($I$), the tidal deformability (Love), the quadrupole moment ($Q$), and the compactness ($C$), known as the “I-Love-Q” and “C-Love” relations Yagi and Yunes (2013a, b); Maselli et al. (2013). These relations are EoS-insensitive to better than 1% and 6% respectively, and they have been important applications in both GW astrophysics Kumar and Landry (2019) and experimental relativity Yagi and Yunes (2013a, b); Gupta et al. (2018); Doneva and Pappas (2018). For example, these relations and the measurement of the tidal deformabilities allow for the inference of several other stellar properties, such as the moment of inertia, the compactness, the spin ($\chi$), and the radius of NSs Kumar and Landry (2019). Analyses such as that of Kumar and Landry (2019) could benefit from the improvement of such EoS-insensitive relations, which we derive in this paper.

III EoS-insensitive relations

In this section, we repeat and improve the analyses of Yagi and Yunes (2016a, 2013b) through the use of the two sets of EoSs (the constrained and the unconstrained set) and through new fitting functions that properly limit Newtonian results.

III.1 I-Love-Q relations

Here we present our results on the I-Love-Q universality, comparing when possible to Fig. 1 of Yagi and Yunes (2013b). In particular, we consider two distinct classes of NSs: nuclear matter EoSs and hybrid quark-hadron star EoSs as described in Sec. II. We begin in Sec. III.1.1 by fitting the new I-Love-Q relations using the constrained set of EoSs. This is followed in Sec. III.1.2 by an analysis and discussion of how well the hybrid star EoSs agree with the improved binary Love relations.

III.1.1 Nuclear matter stars

Following Yagi and Yunes (2013b), we first fit the data for each EoS-insensitive relation to the following function:

[TABLE]

where $y$ and $x$ correspond to NS observables $\bar{I}$, $\bar{Q}$, and $\Lambda$, and the updated coefficients are given in the top of Table 1. This fitting function, however, does not limit properly to the Newtonian results Yagi and Yunes (2013b):

[TABLE]

where $\Lambda^{-1/5}\sim C$ when $C\ll 1$. One can thus improve the fitting function to

[TABLE]

where $\alpha$ is either $2/5$, $1/5$, or $2$ for the $\bar{I}-\Lambda$, $\bar{Q}-\Lambda$, and $\bar{I}-\bar{Q}$ relations respectively, and where the new fitting coefficients are presented in the bottom of Table 1. While the two fits result in similar $R^{2}$ values555$R^{2}$ is the coefficient of determination, defined as $\sum_{i}(f_{i}-\bar{y})^{2}/\sum_{i}(y_{i}-\bar{y})^{2}$, where $\bar{y}$ is the mean data value, and $f_{i}$, $y_{i}$ are the modeled and actual data values of $\sim 0.999995$, the fit in Eq. (8) has the advantage that it properly limits to the Newtonian result as $\Lambda\gg 1$ Yagi and Yunes (2016a).

Figure 3 shows the improved I-Love-Q relations between the dimensionless moment of inertia $\bar{I}\equiv I/M^{3}$, the dimensionless quadrupole moment $\bar{Q}\equiv Q/M^{3}$, and the dimensionless tidal deformability $\Lambda$. The fits are done using the new fitting function in Eq. (8) and over either the constrained set or the unconstrained set separately. The bottom panels show the relative fractional difference between the fit and the data for each set of EoSs. The new fits to the constrained set shows considerably more EoS-insensitivity than the fit to the unconstrained set. These results are summarized in Table 2, which tabulates the maximum EoS variation in each fit. Clearly then, the EoS-insensitive relations can be made more universal by restricting the EoSs through the use of observations. Such improvements can be beneficial to future studies in GW astrophysics, such as those of Kumar and Landry (2019), and in experimental relativity, such as those of Gupta et al. (2018).

III.1.2 Hybrid quark-hadron stars

Let us now focus on the I-Love-Q relations of hybrid stars, and their compatibility with their nuclear matter counterparts. For concreteness, we consider three different sets of EoSs in the hybrid case:

the complete set of 100 constrained EoSs combined with the 10 hybrid star EoSs, 2. 2.

the complete set of 100 constrained EoSs alone, 3. 3.

the complete set of 10 hybrid star EoSs alone.

For each of these cases, we compute the I-Love-Q relations, we fit the data to Eq. (8) and we compute the relative fractional difference.

The fractional differences of hybrid stars from the fit for the second case (fit to only the constrained EoSs) is shown with dashed green lines in Fig. 3, while the maximum EoS variation for the I-Love relation is shown in Table 3 for both the constrained and hybrid star cases. The hybrid star EoSs obey the I-Love-Q relations in each case to better than $\sim 1.7$%, a variability that is slightly higher than that found when using only nuclear matter EoSs in previous works Yagi and Yunes (2013a, b), and is consistent with Paschalidis et al. (2018). The universality cannot be improved much through the introduction of new fits, only bringing the maximal EoS variation down to $\sim 1$% for the fits constructed with only hybrid star EoSs. From this study, we conclude that hybrid star EoSs do obey the traditional nuclear I-Love-Q relations computed with nuclear EoS data, albeit with a slight decrease in universality to $\sim 1.7$%. We find this decrease in universality to be somewhat consistent with similar works Paschalidis et al. (2018); Bandyopadhyay et al. (2018); Lau et al. (2017); Han and Steiner (2018). In such investigations, departures from universality of up to $\sim 1\%-2\%$ were found with various hybrid star EoSs, compatible with our result of $\sim 1.7$%.

III.2 C-Love Relations

In this subsection we focus on the EoS-insensitive C-Love relations, as introduced in Yagi and Yunes (2013a, b); Maselli et al. (2013), as well as on the R-Love relations, since these play a key role in GW data analysis for the extraction of M-R credible intervals. As in the previous subsection, we consider both nuclear matter EoSs, as well as hybrid quark-hadron star EoSs.

III.2.1 Nuclear matter stars

Following Ref. Maselli et al. (2013), we begin by fitting the data for each set of EoSs to the simple curve

[TABLE]

Doing so, yields $a_{0}=0.3617$, $a_{1}=-0.03548$, and $a_{2}=0.0006194$ for the constrained set of EoSs, similar to what was found in Maselli et al. (2013). As in the I-Love-Q case, however, the above fitting function does not limit properly to the Newtonian result

[TABLE]

We thus repeat the fit using Eq. (8) as the fitting function, with the new fitting coefficients presented in Table 1.

The left panel of Fig. 4 shows the C-Love relations for both the constrained and unconstrained sets, along with the corresponding absolute differences666We present absolute differences instead of fractional differences since the former is what matters directly to the GW data analysis. from the fits (instead of the fractional differences as done back in Fig. 3). The fit to the constrained EoSs suppresses the EoS variability compared to the fit to the unconstrained set, as well as that of previous work Maselli et al. (2013). The maximal EoS variation is compared between these three cases in Table 2.

From the C-Love relations, we can also compute directly the R-Love relation using $R(\Lambda)=M/C(\Lambda)$ for NSs of mass $M$. The right panel of Fig. 4 shows the R-Love relations for the constrained and unconstrained sets, with the bottom panel showing the absolute difference of the data and the fits. The C-Love relations allow us to infer the NS radius to better than $\sim 350\textrm{m}$ in the constrained case, while the error goes up to $1,000\textrm{m}$ in the unconstrained case, as also tabulated in Table 2. The systematic uncertainty in the radius using the constrained fit is thus comparable to the $\sim 140\textrm{m}$ systematic uncertainty in the radius due to the choice of EoS to model the crust Gamba et al. (2019). See also Refs. Annala et al. (2018); Raithel et al. (2018) for related work on the R-Love relations.

III.2.2 Hybrid quark-hadron stars

Let us now focus on whether the C-Love relations hold for hybrid stars. Much like in Sec. III.1.2, we perform 3 separate fits: one to the constrained set plus 10 hybrid EoSs, another to the constrained set only, and a third to the 10 hybrid EoSs only. We then compare the EoS insensitivity in each case.

In the top left panel of Fig. 4, we show the C-Love relations for hybrid stars (dashed green curves), while in the bottom panel, we show the absolute difference of the relations for such stars from the fit constructed only from the constrained set of EoSs (dashed green curves). The C-Love relations for hybrid stars remain EoS-insensitive, but the degree of universality is not as high as in the case of hadronic NSs. Table 4 compares the maximal EoS variability for the constrained and hybrid star EoSs fitting to the three data sets described above. As in Sec. III.1.2, the maximal EoS variation for hybrid stars fluctuates only slightly ($\sim 4.5\%-7\%$) in each case. From this, we conclude that hybrid stars do obey the C-Love relation derived for nuclear matter stars, with the caveat that the maximum universality increases to $\sim 7.1\%$. For completeness, we also show the R-Love relations for hybrid stars in Fig. 4. The absolute differences from the (constrained EoS) fit are similarly displayed in the bottom panel, showing a maximum EoS variability of $\sim 50\text{ m}$, consistent with the uncertainties found in the unconstrained relations.

III.3 Binary love relations

Let us now consider the binary Love relations. As for other relations, we consider nuclear matter stars and hybrid stars separately.

III.3.1 Nuclear matter stars

Following Ref. Yagi and Yunes (2016a), we begin by fitting the binary Love relations to the constrained and unconstrained sets using the two-dimensional curve:

[TABLE]

where $q\equiv m_{2}/m_{1}$ is the mass ratio with $m_{2}\leq m_{1}$, and $F_{n}(q)$ is the Newtonian-limiting controlling factor, given by

[TABLE]

The coefficients of the fit to the constrained set are

[TABLE]

where we set $n=0.743$ and $\alpha=1$ as was done in Ref. Yagi and Yunes (2016a). Unlike the individual I-Love-Q relations, the binary Love relations depend also on the mass ratio.

Figure 5 shows the improved binary Love relation for 3 different mass ratios ($q=0.9$, 0.75 and 0.5) for both the constrained and unconstrained sets of EoSs. Once again we find that the constrained set shows a considerable increase in EoS-insensitivity relative to both the unconstrained set and previous works Yagi and Yunes (2016a). As before, the fit to the unconstrained set shows similar, yet slightly larger EoS variation than that found in previous work Yagi and Yunes (2016a) due to the random sampling used in our analysis. The absolute error between the fit and the data is approximately independent of the mass ratio $q$, with a subtle maximum at $q=0.75$, which mimics the behavior of the controlling factor $F_{n}(q)$. As before, maximum EoS variation is tabulated in Table 2 for each value of mass ratio considered here.

We can extract two additional conclusions from Fig. 5. First, as the mass ratio becomes more extreme (i.e. away from unity), the ranges of values that $\Lambda_{s}$ and $\Lambda_{a}$ can take decrease. This is because as the separation between $m_{1}$ and $m_{2}$ increases, the effective number of possible $\Lambda_{1}(m_{1})$-$\Lambda_{2}(m_{2})$ configurations decrease, and thus, the allowed parameter space shrinks. Second, the fit to the constrained set does not extend to as high values of $\Lambda_{s}$ as the fit to the unconstrained set does. This is because of how the sets of EoSs were generated in our analysis, with the unconstrained set including some EoSs that are stiffer777Typically, the “stiffness” of an EoS is determined by the amount of pressure gained given an increase in density. Stiff EoSs have steep pressures and predict larger maximum mass NSs, while soft EoSs have shallow pressures and predict smaller maximum mass NSs. and others that are softer than those constrained by the GW170817 event. This results in the unconstrained set containing a larger range of tidal deformabilities $\Lambda_{1,2}$ (from both above and below) than the constrained set, which ultimately leads to a larger range of $\Lambda_{s}=\frac{1}{2}(\Lambda_{1}+\Lambda_{2})$ values.

III.3.2 Hybrid quark-hadron stars

Let us now consider the binary Love relations for hybrid stars. As described in Sec. II, transitional quark-hadron matter stars undergo strong first-order phase transitions at a transitional pressure $P_{\text{tr}}$, where the hadronic branch departs into a quark-matter branch at the corresponding transitional mass $M_{\text{tr}}$. These transitions result in two distinct types of NSs, based on their observed mass: (i) massive ($M\geq M_{\text{tr}}$) hybrid stars that have quark-matter inner cores and nuclear matter elsewhere (which recall we denote “HS”); and (ii) low-mass ($M\leq M_{\text{tr}}$) hadronic stars with no internal transition to quark matter (which recall we denote “NS”). We therefore expect the binary Love relations for hybrid stars to present behavior identical to their purely hadronic counterparts below the transitional mass (or correspondingly above a transitional $\Lambda_{s,\textrm{tr}}$), and very different behavior for higher masses (or correspondingly for smaller $\Lambda_{s,\textrm{tr}}$).

This behavior is exactly what we observe in Fig. 5 for the binary Love relations of hybrid stars. Indeed, contrary to the case of EoS-insensitive relations for isolated hadronic stars (I-Love-Q, C-Love, and R-Love) which remain moderately EoS-insensitive for all hybrid stars, the binary Love relations depart from their hadronic counterparts below some transitional $\Lambda_{s}$. This is because at low pressures ($\Lambda_{s}>\Lambda_{s,\textrm{tr}}$) the binary is purely of NS/NS type (with the EoS for both stars a member of the constrained set), but once the critical pressure is reached ($\Lambda_{s}<\Lambda_{s,\textrm{tr}}$), one or both stars transition into the hybrid branch, resulting in large reductions in tidal deformability, as shown in the left panel of Fig. 4. When one star lies on the hadronic-matter branch while the other is on the quark-matter branch, the difference between tidal deformabilities becomes larger than expected for pure hadronic-matter stars. This results in large deviations in the sums and differences of tidal deformabilities $(\Lambda_{2}\pm\Lambda_{1})$, disrupting the overall universality and generating the large “bump” in the binary Love relations for hybrid stars888This discrepancy is not present in the I-Love-Q and C-Love relations due to their single-star nature. seen in Fig 5.

From this analysis, we conclude that binaries containing one or more hybrid stars do not satisfy the same binary Love relations as binaries that contain only traditional nuclear matter stars. This does not imply, however, that more sophisticated binary Love relations cannot be constructed that will remain EoS-insensitive and able to model both types of binary systems. Such relations, however, would necessarily have to include (at least) one new parameter that determines the transition between the hadronic and the quark branches, such that the “bump” in the binary Love relations can be properly modeled. Work along these lines is outside the scope of this paper.

IV Impact on future observations

We have now shown that binary NS merger observations can help improve the EoS-insensitive relations, but the question remains: is it worth it? Current interferometer sensitivities are not yet high enough to accurately constrain the dominant tidal parameter $\tilde{\Lambda}$. For example, GW170817 was detected by the second LIGO observing run (“O2”) aLI and Virgo Acernese et al. (2015), and was able to constrain $\tilde{\Lambda}$ to a $90\%$ credible interval centered at $\mu_{\tilde{\Lambda}}=395$ and with a width of 325 Abbott et al. (2019) (or $\sigma_{\tilde{\Lambda}}\approx 198$). This corresponds to statistical uncertainties of ${\cal{O}}(80\%)$, which dominates the error budget compared to the small systematic uncertainties picked up by EoS variation in the EoS-insensitive relations. This implies that currently, the use of improved EoS-insensitive relations will only make a negligible difference on the extraction of tidal parameters.

In this section, we further explore this question and study when in the future the new set of improved relations will become important as current detectors are improved and new ones are built. In Sec. IV.1, we first estimate the systematic uncertainties introduced by using the improved binary Love relations. In Sec. IV.2, we estimate the statistical uncertainties on the extraction of tidal parameters, and compare them to the above-mentioned systematic uncertainties. This is repeated for 5 future detectors, where multiple detections become important, but to combine constraints from multiple events we cannot use the $\tilde{\Lambda}$ parameterization, as this depends on the masses of the binary constituents. To remedy this, we re-parameterize the waveform in terms of the coefficients $\lambda_{0}$ (sometimes called $\Lambda_{1.4}$) and $\lambda_{1}$, obtained from a Taylor expansion of the dimensionless tidal deformability $\Lambda$ about a “canonical” reference mass $M_{0}=1.4\text{ M}_{\odot}$ Del Pozzo et al. (2013); Yagi and Yunes (2016a):

[TABLE]

The Taylor coefficients $\lambda_{0}$ and $\lambda_{1}$ are mass independent, and thus they are identical in all binary NS observations, and their posteriors may be combined.

IV.1 Error Marginalization

Although the improved binary Love relation $\Lambda_{a}^{\text{relation}}(\Lambda_{s},q)$ shows a high degree of universality, any residual EoS dependence in the relation could introduce a systematic bias and lead to incorrect inferences about the correct EoS. Let us then discuss a method to marginalize over the residual EoS-dependence of the binary Love relations. The residual is here defined as

[TABLE]

where $\Lambda_{a}^{\text{relation}}(\Lambda_{s},q)$ is given by the binary Love fit found in Sec. III.3, while $\Lambda_{a}^{\text{true}}$ is the true value predicted by choosing a particular EoS in the constrained set and solving for the tidal deformabilities numerically.

Following the proposal in Chatziioannou et al. (2018), one can model the residual EoS-sensitivity by enhancing the binary Love relations through

[TABLE]

where $\mathcal{N}(\mu_{r},\sigma_{r})$ is a normal distribution with mean and variance $\mu_{r}$ and $\sigma_{r}^{2}$. Let us further assume that the residuals $r$ obey a Gaussian distribution with mean and standard deviation that can be decomposed as

[TABLE]

where $\mu_{r}(q)$ and $\mu_{r}(\Lambda_{s})$ are the means marginalized over $\Lambda_{s}$ and $q$ respectively, while $\sigma_{r}(q)$ and $\sigma_{r}(\Lambda_{s})$ are the standard deviations marginalized over $\Lambda_{s}$ and $q$ respectively.

Clearly then, to account for the residual EoS-sensitivity, we must first find the marginalized mean and standard deviation of the residual function, which we accomplish as follows. We first generate $\Lambda_{s}^{\textrm{true}}$ data by sampling through various values of mass ratio $q\in[0.36,1]$ and symmetric tidal deformability $\Lambda_{s}\in[4,4600]$ over all 100 elements of the constrained set of EoSs. For each value of $(q,\Lambda_{s})$ we then compute $\Lambda_{a}^{\textrm{relation}}$ from the binary Love relation, and then compute the residuals through Eq. (15). We then proceed to marginalize over $q$ by binning the residuals in $\Lambda_{s}$, which returns a distribution of residuals that is only a function of the binned $\Lambda_{s}$, and from which we can compute the mean $\mu_{s}(\Lambda_{s})$ and the standard deviation $\sigma_{r}(\Lambda_{s})$. Repeating this procedure by marginalizing over $\Lambda_{s}$ then allows us to compute the mean $\mu_{s}(q)$ and the standard deviation $\sigma_{r}(q)$. Finally, we fit the means and standard deviations to the functions Chatziioannou et al. (2018)

[TABLE]

where the fitting parameters $\mu_{i}$ and $\sigma_{i}$ are tabulated in Table 5.

In what follows, however, we will be interested in comparing an estimate of the systematic uncertainties in $\lambda_{0}$ due to the residual EoS-sensitivity in the binary Love relations to the statistical error in the extraction of this parameter. To estimate the former, we first numerically calculate $\lambda_{0}^{\textrm{true}}$ by sampling again in $(q,\Lambda_{s})$ within the same ranges as before and over all 100 elements of the constrained and the unconstrained sets of EoSs. We then use the binary Love relation to calculate $\lambda_{0}^{\text{relation}}$ over the $(q,\Lambda_{s})$ sampled points, and from this we compute the residual.

[TABLE]

Figure 6 shows a histogram of ${\cal{R}}$ for both the constrained and unconstrained set of EoSs. The standard deviations of the two histograms, $\sigma=9.764$ and $\sigma=78.28$, show a large decrease in EoS-variability from the unconstrained to the constrained sets of EoSs. From this distribution, we can additionally find the 90%, 99%, and 100% credible interval on $\lambda_{0}$, namely $P_{90}=13.19$, $P_{99}=42.38$, and $P_{100}=111.32$ for the constrained EoSs.

Figure 7 shows the standard deviations of the $\lambda_{0}$ residuals binned in $q$ (top panel) and $\Lambda_{s}$ (bottom panel), as was done previously for the $\Lambda_{a}$ residuals. The EoS-variability in $\lambda_{0}$ is dominated by the region around $q\sim 0.5$ and $\Lambda_{s}\sim 2000$. Henceforth, in order to take into account these high-error regions of parameter space which may get averaged out, we approximate the EoS-variability in the binary Love relations with the 90% credible interval (rather than the $1\sigma$ uncertainty), which for the $\lambda_{0}$ parameter is $P_{90}=13.19$ (indicated by the dashed green line in Fig. 7 and the dashed indigo line in Fig. 1). We can similarly derive a conservative estimate of the systematic uncertainty due to EoS-variability in the binary Love relations in terms of $\tilde{\Lambda}_{a}$ to find $12.01$, which is what we use later in Fig. 9. Here we note that these results rely on the assumption of a purely hadronic EoS, as the constrained set of EoSs do not contain any hybrid stars. This assumption was also made in Ref. Abbott et al. (2018a) upon the derivation of the posterior probability distribution.

IV.2 Future Observations

We now estimate the feasibility of using the improved EoS-insensitive relations in future GW observations of the coalescence of binary NSs. We estimate the statistical accuracy to which parameters can be extracted through a simple Fisher analysis Finn (1992); Cutler and Flanagan (1994), assuming sufficiently high signal-to-noise ratio and Gaussian noise Cutler and Flanagan (1994); Berti et al. (2005); Poisson and Will (1995).

Consider then a waveform template $h$ parameterized in terms of

[TABLE]

where $\eta\equiv m_{1}m_{2}/m^{2}$ is the symmetric mass ratio with $m_{1,2}$ and $m$ being the individual and total masses, $\mathcal{M}=m\eta^{3/5}$ is the chirp mass, $A\equiv{\mathcal{M}^{5/6}}/({\sqrt{30}\pi^{2/3}D_{L}})$ is a normalized amplitude factor with $D_{L}$ the luminosity distance, and $\chi_{s,a}=\frac{1}{2}(\chi_{1}\pm\chi_{2})$ are the symmetric and antisymmetric dimensionless non-precessing spin parameters. Given a sufficiently loud observation and Gaussian noise, the resulting posterior distribution on the estimated parameters is Gaussian with a root-mean-square error given by

[TABLE]

where the Fisher information matrix $\tilde{\Gamma}$ is

[TABLE]

Here $\sigma_{\theta^{a}}$ are the parameters’ prior root-mean-square estimate, and the inner product is defined by:

[TABLE]

with the overhead tilde and the star representing the Fourier transform and complex conjugation respectively.

The inner product, and thus the root-mean-square error, depends on the spectral noise density $S_{n}(f)$ and the waveform template $\tilde{h}$. We consider here 6 different interferometer designs, O2 aLI , aLIGO at design sensitivity aLI , A+ Ap_ , Voyager Ap_ , CE ET , ET-D Ap_ shown in Fig. 8, in order to compare the statistical uncertainties accrued on parameter extraction using future upgraded LIGO detectors, as well as third generation detectors. The inner product also depends on the particular waveform model we use for the template. We here consider a “PhenomD” (IMRD) waveform template Khan et al. (2016); Husa et al. (2016) modified by a 6PN tidal correction Wade et al. (2014a) (IMRD + 6PN), as well as a “PhenomD” template modified with a NRTidal correction Samajdar and Dietrich (2018) (IMRD + NRTidal). Considering two template models will allow us to estimate systematic uncertainties due to mismodeling of the GW signal.

We begin by testing our Fisher analysis (with an IMRD + 6PN waveform injection) against a simulated event identical to GW170817 with O2 detector sensitivity aLI . Because only 1 event was detected, we use $\tilde{\Lambda}$ and $\delta\tilde{\Lambda}$ as the tidal parameters for comparison purposes. Further, we scale the luminosity distance such that the signal-to-noise-ratio ($SNR\equiv\rho$) is fixed to $\rho=32.4$, as found in GW170817. We used fiducial template parameter values of $\mathcal{M}=1.22\text{ M}_{\odot}$ for the chirp mass, $\eta=0.249$ for the symmetric mass ratio, $\tilde{\Lambda}=395$ (Corresponding to GW170817), and [math] for the remaining parameters. We also assume low spin priors $|\chi|\leq 0.05$, as well as $\tilde{\Lambda}\leq 3000$ and $|\delta\tilde{\Lambda}|\leq 500$ Wade et al. (2014b). The resulting $90\%$ credible region of the posterior distribution on $\tilde{\Lambda}$ has a range of $\pm 276.99$, which is in close agreement to that found by LIGO Abbott et al. (2017a, 2019)

Next, we consider events similar to GW170817 detected on upgraded detectors and future detectors, as well as the combined statistical uncertainties of $N_{A}$ events detected over a 1 year observation with each detector. The latter is calculated by integrating the local binary NS merger rate over redshifts up to the horizon redshift of each detector. Similar to before, we used fiducial template parameter values of $\mathcal{M}=1.22\text{ M}_{\odot}$ for the chirp mass, $\eta=0.249$ for the symmetric mass ratio, $\lambda_{0}=150$, $\lambda_{1}=-213$ (corresponding to the GW170817), and [math] for the remaining parameters. Also as before, for each Fisher calculation we assume low spin priors $|\chi|\leq 0.05$, as well as $0\leq\lambda_{0}\leq 3207$ and $-4490\leq\lambda_{1}\leq 0$ Del Pozzo et al. (2013)999These are converted from the dimensional forms found in Del Pozzo et al. (2013) into their corresponding dimensionless forms.. The process we use to compute single and combined statistical uncertainties for each detector sensitivity $S_{n}^{A}(f)$ is detailed in App. A. From this, we determine if and when the statistical uncertainties associated with the parameter extraction of $\lambda_{0}$ drop below the systematic EoS variation errors from using the binary Love relations.

Figures 1 and 9, and Table 6 summarize our results graphically using an IMRD + 6PN template. In Fig. 9, one can see the statistical accuracy to which $\tilde{\Lambda}$ can be estimated as a function of $\rho^{A}_{\text{GW170817}}$ for single GW170817-like events, i.e. the signal-to-noise ratio that future detector $A$ would have measured for a GW170817 event. The systematic uncertainty in $\tilde{\Lambda}$ (the horizontal dashed line) becomes comparable to the statistical uncertainty for Voyager-class detectors or better. Figure 1 presents similar results but for the parameter $\lambda_{0}$, for which one can combine posteriors and obtained a reduced combined statistical error. This figure also shows that the statistical and systematic uncertainties become comparable for Voyager-class detectors or better. The careful reader will notice that Figs. 1 and 9 show that the single-event statistical uncertainties using CE is higher than using ET despite the latter having a larger SNR for a GW170817-like event. This is because the 3-detector geometry of ET enables it to have higher sensitivity than CE at frequencies above 300Hz, which is precisely where the tidal deformabilities are encoded, but ET has lower sensitivity at lower frequencies, where a lot of the signal-to-noise ratio accumulates.

Is the statistical error calculated here robust to mismodeling systematics in the template? We repeat the analysis above using an IMRD + NRTidal template model and find results consistent with those presented above. Indeed, Fig. 1 shows the individual and combined statistical error on $\lambda_{0}$ when using this template model (maroon circles and dashed region). The accuracy to which $\lambda_{0}$ can be measured with this template model is systematically better than when using the IMRD+6PN model (by roughly a factor of two). This implies that more accurate template models will indeed be required in the third-generation detector era, as previously pointed out in Samajdar and Dietrich (2018). However, our conclusion that for Voyager-class detectors or better the statistical and binary Love systematic uncertainties become comparable seems robust.

V Conclusion and Discussion

The recent GW observation of a binary NS coalescence, GW170817, placed constraints on the supranuclear matter EoS for NSs. We used this observation to generate a restricted set of spectral EoSs that agree with it, in order to reduce the uncertainties upon the extraction of tidal parameters from future GW events. Previous work by Yagi and Yunes Yagi and Yunes (2013a, b, 2016b, 2016a) had found EoS-insensitive relations between symmetric and antisymmetric combinations of NS tidal deformabilities, which aid in the extraction of said tidal parameters. We here found that the GW170817-constrained set of EoSs are more EoS-insensitive by a factor of $\sim 60$% for stars with mass ratios of $0.75$ relative to previous work. Similarly, we find an increase in EoS-insensitivity in the C-Love and I-Love-Q relations by factors of $\sim 75$% and $\sim 50$% respectively. The former further allowed us to improve the R-Love relation leading to uncertainties below $400$ meters in the entire parameter space. We also studied the EoS-insensitive relations of hybrid stars and found that for the most part the relations remain insensitive for isolated stars, albeit with slightly higher EoS variability. The binary Love relations, however, do not satisfy the same EoS-insensitive relations as in the case of purely hadronic star, when the binary contains at least one hybrid star.

The second half of this paper focused on when the improvements in the EoS-insensitive relations would become necessary in future detectors. Current detectors are not yet sensitive enough that the systematic uncertainties in the EoS-insensitive relations make a large difference. However, we did find that for Voyager-class detectors or better, the systematic uncertainties in the EoS-insensitive relations due to EoS variability become comparable to statistical uncertainties in the estimation of $\tilde{\Lambda}$ and $\lambda_{0}$. We also considered the effect of waveform mismodeling and found that the above conclusion remains robust, but that more accurate waveform models will be necessary to take full advantage of the improved sensitivity of future detectors.

Future work on this subject could entail an investigation into the improvement of alternative EoS-insensitive relations, such as the multipole Love relations between various $\ell$-th order electric, magnetic, and shape tidal deformabilities, as discussed in Yagi (2014). Lackey et al. Lackey et al. (2018, 2017) presented surrogate models of non-spinning effective-one-body waveforms with the use of universal relations. By reducing the number of waveform model parameters, surrogate models aid in the extraction of NS observables from GW detections. The improvement in the multipole Love relations can then be used to increase the accuracy of such surrogate models. Another possible avenue for future research includes a more comprehensive analysis into the intricacies of new hybrid star binary Love relations.

Acknowledgments

KY acknowledges support from NSF Award PHY-1806776. K.Y. would like to also acknowledge networking support by the COST Action GWverse CA16104. N.Y. acknowledges support from NSF grant PHY-1759615 and NASA grants NNX16AB98G and 80NSSC17M0041. This research has made use of data, software and/or web tools obtained from the Gravitational Wave Open Science Center (https://www.gw-openscience.org), a service of LIGO Laboratory, the LIGO Scientific Collaboration and the Virgo Collaboration. LIGO is funded by the U.S. National Science Foundation. Virgo is funded by the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale della Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions by Polish and Hungarian institutes. The authors are grateful for computational resources provided by the LIGO Laboratory and supported by National Science Foundation Grants PHY-0757058 and PHY-0823459.

Appendix A Computation of statistical uncertainties

In this appendix, we detail the process used to compute the statistical uncertainties on the extraction of $\lambda_{0}$ from the gravitational waveform. This is accomplished with a simple Fisher analysis described in Sec. IV.2, where we first consider single-events similar to GW170817 as if they were detected on future detectors $A\equiv($ O2 aLI , aLIGO aLI , A+ Ap_ , Voyager Ap_ , ET ET , CE Ap_ $)$. We conclude by simulating a population of $N_{A}$ events for each interferometer, approximating the number of events detected on interferometer $A$ over an observing period of one year. This is used to combine the statistical uncertainties, resulting in an approximation on the overall measurement accuracy of $\lambda_{0}$, which is compared to the systematic uncertainties computed in Sec. IV.1. The process used to achieve this is outlined below:

1. (i)

Perform a Fisher analysis as outlined in Sec. IV.2 using detector sensitivity $S_{n}^{A}(f)$, while restricting the luminosity distance $D_{L}$ such that an SNR of $\rho^{\text{O2}}_{\text{GW170817}}=32.4$ would be achieved on O2 sensitivity $S_{n}^{\text{O2}}(f)$. Here we assume low spin priors $|\chi_{1,2}|\leq 0.05$, as well as $0\leq\lambda_{0}\leq 3207$ and $-4490\leq\lambda_{1}\leq 0$ Del Pozzo et al. (2013) (these are converted from their corresponding dimensional forms). This results in an SNR $\rho^{A}_{\text{GW170817}}$ and a single-event statistical uncertainty $\sigma_{\text{GW170817}}^{A}$ accrued in the extraction of $\lambda_{0}$ on detector $A$. 2. (ii)

Generate a population of $N_{A}$ events corresponding to the expected binary NS merger detection rate for detector $A$, following the probability distribution Schutz (2011); Chen and Holz (2014):

[TABLE]

with a network SNR threshold of $\rho_{\text{th}}=8$. The number of events $N_{A}$ is calculated by taking into account the BNS merger rate history throughout all redshift values within detector As’ horizon redshift $z_{h}$, as shown by Eq. (10) of Ref. Cutler and Harms (2006):

[TABLE]

Here, $a_{0}r_{1}(z)$, $\frac{d\tau}{dz}$, and $r(z)$ for our chosen cosmology are given by:

[TABLE]

[TABLE]

[TABLE]

where $H_{0}=70\text{km s}^{-1}\text{Mpc}^{-1}$ is the local Hubble constant and $\Omega_{\Lambda}=0.67$ is the universe’s vacuum energy density. Here we choose an observing period of $\Delta\tau_{0}=1$ year, and calculate the detection rate for the upper, central, and lower limits of the local binary NS coalescence rate density $\mathcal{R}=1540^{+3200}_{-1220}\text{ Gpc}^{-3}\text{yr}^{-1}$ Abbott et al. (2017b), giving the rates $N_{A}$ shown in the second column of Table 6. 3. (iii)

Compute the combined population standard deviation $\sigma_{N_{A}}$, taking into account sources at varying redshifts as was done in Eq. (3) of Ref. Yagi and Tanaka (2010):

[TABLE]

Here, we compute the populations’ single-event uncertainties $\sigma_{i}^{A}(z)$ via the simple SNR scaling factor:

[TABLE]

where $\rho_{i}^{A}$ is the SNR of simulated event $i$ computed from Eq. (28), and $\rho_{\text{GW170817}}^{A}$ and $\sigma_{\text{GW170817}}^{A}$ are the known SNRs and uncertainties of the GW170817 event from step (i). This results in the combined statistical uncertainties for the lower, central, and upper limits of the local binary NS coalescence rate density $\mathcal{R}$, depicted by the cyan shaded region in Fig. 1.