Upper limit set by causality on the tidal deformability of a neutron star

TL;DR

This paper establishes a causality-based upper limit on the tidal deformability of neutron stars, constraining gravitational-wave signals and informing the neutron star equation of state.

Contribution

It derives a fundamental upper bound on neutron star tidal deformability based on causality constraints, linking the EOS to observable gravitational-wave features.

Findings

Causality imposes a strict upper limit on tidal deformability.

Upper limit on phase shift in black hole-neutron star mergers.

Constraints remain strong even with high maximum neutron star mass assumptions.

Abstract

A principal goal of gravitational-wave astronomy is to constrain the neutron star equation of state (EOS) by measuring the tidal deformability of neutron stars. The tidally induced departure of the waveform from that of point-particle (or spinless binary black hole (BBH)) increases with the stiffness of the EOS. We show that causality (the requirement that the speed of sound is less than the speed of light for a perfect fluid satisfying a one-parameter equation of state) places an upper bound on tidal deformability as a function of mass. Like the upper mass limit, the limit on deformabity is obtained by using an EOS with $v_{sound} = c$ for high densities and matching to a low density (candidate) EOS at a matching density of order nuclear saturation density. We use these results and those of [B.D. Lackey et al., Phys. Rev. D 89, 043009 (2014)] to estimate the resulting upper limit on the gravitational-wave phase shift of a black hole-neutron star (BHNS) binary relative to a BBH. Even for assumptions weak enough to allow an maximum mass of $4\ M_\odot$ (a match at nuclear saturation density to an unusually stiff low-density candidate EOS), the upper limit on dimensionless tidal deformability is stringent. It leads to a still more stringent estimated upper limit on the maximum tidally induced phase shift prior to merger. We comment in an appendix on the relation between causality, the condition $v_{sound} < c$, and the condition $dp/d\epsilon < 1$ for the effective EOS governing the equilibrium star.

Equations

$$\frac{dp}{dϵ}\le 1.$$

$$p-p_{\mathrm{match}}=ϵ-{ϵ}_{\mathrm{match}},$$

$$p(ϵ)=\{\begin{array}{ll}{p}_{MS1}(ϵ),& ϵ\le {ϵ}_{\mathrm{match}}\\ ϵ-{ϵ}_{\mathrm{match}}+p_{MS1}({ϵ}_{\mathrm{match}}),& ϵ\ge {ϵ}_{\mathrm{match}}.\end{array}$$

$$\left(1-\frac{2m}{r}\right)\frac{dp}{dr}=-\frac{1}{r^2}(ϵ+p)(m+4\pi r^{3}p),$$

$$\frac{dm}{dr}=4\pi r^2ϵ.$$

$$\delta g_{\alpha \beta }={\delta }_{\mathrm{external}}{g}_{\alpha \beta }+{\delta }_{\mathrm{induced}}{g}_{\alpha \beta },$$

$${\delta }_{\mathrm{external}}{g}_{tt}$$

$${\delta }_{\mathrm{induced}}{g}_{tt}$$

$$Q_{ij}=-\lambda {\mathcal{E}}_{ij},$$

$$k_2=\frac{3\lambda }{2R^5}.$$

$$\mathbf{g}=-e^{2\nu }d{t}^2+\frac{1}{1-2m/r}d{r}^2+r^2(d{\theta }^2+{\sin }^2\theta d{\varphi }^2),$$

$$\left(1-\frac{2m}{r}\right)\frac{d\nu }{dr}=\frac{1}{r^2}(m+4\pi r^{3}p),$$

$$\delta \mathbf{g}=$$

$$+r^2(d{\theta }^2+{\sin }^2\theta d{\varphi }^2)K{Y}_{2,m}(\theta ,\varphi ),$$

$$0=$$

$$-H[\frac{6}{r^2}-4\pi \left(5ϵ+9p+\frac{ϵ+p}{dp/dϵ}\right)$$

$$+4\left(1-\frac{2m}{r}\right){\left(\frac{d\nu }{dr}\right)}^2].$$

$$k_2=$$

$$\times \{2C[6-3Y+3C(5Y-8)+2C^2(13-11Y)$$

$$+2C^3(3Y-2)+4C^4(Y+1)]$$

$${+3(1-2C{)}^2[2-Y+2C(Y-1)]\log (1-2C)\}}^{-1},$$

$$y(r)=r\frac{H^{\prime }(r)}{H(r)},$$

$$\frac{dy}{dr}=$$

$$+\frac{6}{r-2m}-\frac{4\pi r^2}{r-2m}\left[5ϵ+9p+\frac{(ϵ+p{)}^2}{ϵdp/dϵ}\right].$$

$$\Delta {\Phi }_{\mathrm{PN}}=-\frac{3\Lambda }{128\eta }(\pi Mf{)}^{5/3}\left[a_0+a_1(\pi Mf{)}^{2/3}\right],$$

$$a_0$$

$$a_1$$

$$-\sqrt{1-4\eta }\left(1+\frac{4243}{234}\eta -\frac{6217}{234}{\eta }^2-\frac{10}{9}{\eta }^3\right)].$$

$$A=\{\begin{array}{ll}{A}_{\mathrm{PN}},& Mf\le .01\\ A_{\mathrm{PN}}{e}^{-\eta \Lambda B(\Lambda ,\eta ,{\chi }_{\mathrm{BH}})(Mf-.01{)}^3},& Mf>.01\end{array}$$

$$\Delta \Phi =\{\begin{array}{ll}\Delta {\Phi }_{\mathrm{PN}}(Mf),& Mf\le .02\\ -\eta \Lambda E(\eta ,{\chi }_{\mathrm{BH}})(Mf-.02{)}^{5/3}& \\ +\Delta {\Phi }_{\mathrm{PN}}(.02)& \\ +(Mf-.02)\Delta {\Phi }_{\mathrm{PN}}^{\prime }(.02),& Mf>.02.\end{array}$$

$$B=e^{b_0+b_1\eta +b_2{\chi }_{\mathrm{BH}}}+\Lambda e^{c_0+c_1\eta +c_2{\chi }_{\mathrm{BH}}},$$

$$E=e^{g_0+g_1\eta +g_2{\chi }_{\mathrm{BH}}+g_3\eta {\chi }_{\mathrm{BH}}},$$

$$p=ϵ-{ϵ}_{\mathrm{S}},$$

$$M_{\max }\propto {ϵ}_{\mathrm{S}}^{-1/2},R_{\max }\propto {ϵ}_{\mathrm{S}}^{-1/2},$$

$${\lambda }_{\max }\propto {ϵ}_{\mathrm{S}}^{-5/2}.$$

Full text

I Introduction

As the second generation of ground-based interferometric graviational wave detectors (Advanced LIGO, Advanced Virgo, KAGRA, and LIGO-India) approach design sensitivity, we are likely to detect each year the inspiral and coalescence of several compact binary systems that include neutron stars, both BHNS and binary neutron-star (BNS) systems. These observations can constrain the neutron-star equation of state (EOS), which gives the pressure $p$ in terms of the energy density $\epsilon$. A stiffer EOS, where the pressure increases rapidly with density, yields stars with larger radii and larger tidal effects on the waveform, governed by the star’s tidal deformablity. In particular, tidal distortion during inspiral increases with the stiffness of the EOSs. Because energy is lost both to gravitational waves and to the work needed to distort the stars, the inspiral proceeds more rapidly for stars with greater tidal deformability. The result is a waveform in which the increase in frequency is more rapid and in which coalescence occurs sooner – at lower frequency.

Beginning with work by Kochanek Kochanek (1992) and Lai and Wiseman Lai and Wiseman (1996), a number of authors have studied the effect of tides on inspiral waveforms. Simulations Pannarale et al. (2015); Chatziioannou et al. (2015); Zhuge et al. (1996); Uryū et al. (2000); Faber et al. (2002, 2004); Bejger et al. (2005); Gondek-Rosińska et al. (2007); Shibata et al. (2005); Shibata and Uryū (2002); Oechslin and Janka (2007); Taniguchi and Gourgoulhon (2003); Shibata (1999); Shibata and Uryū (2000); Shibata and Taniguchi (2006); Shibata et al. (2003); Miller et al. (2004); Marronetti et al. (2004); Liu et al. (2008) of BHNS and BNS systems and analytic approximations in the context of post-Newtonian theory Vines et al. (2011) and the Effective-One-Body (EOB) formalism Hinderer et al. (2016); Bernuzzi et al. (2015a, b) are nearing the precision needed to extract neutron-star deformability from observations with the projected sensitivity of Advanced LIGO. Recent estimates of the measurability of tidal effects and the ability of these observatories to constrain the EOS with signals from BHNS and BNS systems are given in Hotokezaka et al. (2016); Damour et al. (2012); Del Pozzo et al. (2013); Agathos et al. (2015); Wade et al. (2014) and references therein.

In this work, we obtain the upper limit imposed by causality on the tidal deformability of neutron stars and estimate the resulting constraint on the maximum departure of the waveform of a BHNS inspiral from a corresponding spinless BBH inspiral.111After this paper was posted to arXiv, Moustakidis Moustakidis et al. (2016) pointed out a preprint by him and his coauthors that also obtains upper limits on neutron star mass and tidal deformability imposed by bounds on the speed of sound, including $v_{\mathrm{sound}}\leq c$. However, they use a matching density (described in the next section) 50% higher than ours, giving less conservative results. Furthermore, they do not consider tidal effects during late inspiral, whereas we apply the results of Lackey et al. (2014) to do so. The limit is analogous to the upper limits on neutron-star mass $M_{NS}$ Rhoades Jr and Ruffini (1974); Brecher and Caporaso (1976) and radius $R$ Lattimer (2012). In each case, one assumes an EOS of the form $p=p(\epsilon)$ that is known below an energy density $\epsilon_{\mathrm{match}}$, and one obtains a limit on $M$ and $R$ by requiring that the EOS be causal for $\epsilon>\epsilon_{\mathrm{match}}$ in the sense that the sound speed, given by $\sqrt{dp/d\epsilon}$, must be less than the speed of light. Because the sound speed is a measure of the stiffness of the EOS, this is a constraint on the stiffness. An upper limit on tidal deformability then implies an upper limit on the departure of gravitational wave phase shifts from corresponding waveforms of BBH inspiral.

We use metric signature $-+++$ and gravitational units with $G=c=1$.

II Method

II.1 Causal EOS

For a perfect fluid with a one-parameter EOS $p=p(\epsilon)$, causality implies that the speed of sound, $\sqrt{dp/d\epsilon}$, is less than the speed of light. That is, the dynamical equations describing the evolution of fluid and metric are hyperbolic, with characteristics associated with fluid degrees of freedom lying outside the light cone unless

[TABLE]

There is some inaccuracy in using the one-parameter EOS that governs the equilibrium star to define the characteristic velocities of the fluid, because fluid oscillations with the highest velocities have frequencies too high for the temperature of a fluid element and the relative density $Y_{i}$ of each species of particle to reach their values for the background fluid star. Nevertheless, using a result of Geroch and Lindblom Geroch and Lindblom (1991), we show in Appendix A that causality implies the equilibrium inequality (1) for locally stable relativistic fluids satisfying a two-parameter EOS $p=p(\epsilon,s)$, where $s$ is the entropy per baryon. For the multi-parameter equation of state $p=p(\epsilon,s,Y_{i})$, with $Y_{i}$ the relative density of each species of particle, one must assume without proof that causality implies $v_{\rm sound}<1$; the equilibrium inequality (1) again follows from local stability.

The speed of sound is a measure of the stiffness of the EOS. The well-known upper limit on the mass of a neutron star and a corresponding upper limit on its radius are obtained by using the stiffest EOS consistent with causality and with an assumed known form at low density. That is, above a density $\epsilon_{\rm match}$, the EOS is given by

[TABLE]

where $p_{\rm match}$ is fixed by continuity to be the value of $p$ at $\epsilon_{\rm match}$ for the assumed low-density EOS. The upper limits on mass and radius are then found as functions of the matching density $\epsilon_{\rm match}$.

In this work, we again use an EOS of this form to find an upper limit on neutron-star deformability. To be conservative, as our low-density EOS we choose the MS1 EOS Mueller and Serot (1996), which is among the stiffest candidate equations of state. Our matched causal EOS is then given by

[TABLE]

In computing the deformability, we consider only irrotational neutron stars; and in estimating the effect of tides on the inspiral phase, we neglect resonant coupling of tides to neutron-star modes. Tidal deformation of slowly rotating relativistic stars is treated by Pani et al. Pani et al. (2015); and Essick et al. Essick et al. (2016) argue that tidal excitation of coupled modes may alter the waveform in BNS systems.

II.2 Static, Spherical Stars

We next construct the sequence of static spherical stars based on the causal EOS (3). We numerically integrate the Tolman-Oppenheimer-Volkoff (TOV) equation Oppenheimer and Volkoff (1939),

[TABLE]

where $m(r)$ is the total mass-energy inside radius $r$, related to $\epsilon$ by

[TABLE]

A member of the sequence is specified by its central density $\epsilon_{c}$. Its circumferential radius $R$ is the value of the Schwarzschild coordinate $r$ at which $p(r)=0$, and its gravitational mass is $M=m(R)$.

II.3 Calculating the Tidal Deformability

The departure of the inspiral of a BHNS binary from point-particle (or spinless BBH) inspiral depends on the tidal deformation of the neutron star due to the tidal field of its companion. For large binary separation, the metric near the neutron star can be written as a linear perturbation of the Schwarzschild metric of the unperturbed star that has two parts: The tidal field of the companion, expressed in Schwarzschild coordinates about the center of mass of the neutron star, has the form of an external quadrupole field; and the induced quadrupole distortion of the neutron star gives a second quadrupole contribution to the perturbed metric. That is, outside the support of the star, the quadrupole perturbation is a sum,

[TABLE]

of two time-independent solutions to the field equations linearized about a vacuum Schwarzschild geometry. In a gauge associated with asymptotically Cartesian and mass centered coordinates, the contributions to the perturbed metric have the form

[TABLE]

with no $r^{-3}$ contribution, and

[TABLE]

Here $n^{i}=x^{i}/r$ is an outward-pointing unit vector, $\mathcal{E}_{ij}$ is the tracefree tidal field from the black hole, and $Q_{ij}$ is the neutron star’s induced quadrupole moment. The quadrupole moment tensor $Q_{ij}$ is proportional to $\mathcal{E}_{ij}$,

[TABLE]

and the constant of proportionality $\lambda$ is the tidal deformability of the neutron star. It measures the magnitude of the quadrupole moment induced by an external tidal field and is proportional to the (dimensionless) $\ell=2$ tidal Love number Flanagan and Hinderer (2008)

[TABLE]

After constructing the one-parameter family of spherical stars satisfying Eqs. (3), (4), and (5), we tidally perturb them, compute $k_{2}$ and the radius $R$ of each star, and then find the tidal deformability $\lambda$ from Eq. (10). To calculate $k_{2}$, we use the method described by Hinderer Hinderer (2008): A perturbation of the spherically symmetric background metric

[TABLE]

with $\nu(r)$ determined by

[TABLE]

is found in the Regge-Wheeler gauge Regge and Wheeler (1957), with $\delta\mathbf{g}$ a linear, quadrupolar, static, polar-parity perturbation given by Hinderer (2008) 222Note that, because this gauge does not conform to the constraints of an asymptotically Cartesian and mass-centered chart, there are additional terms in the expansion of the asymptotic metric.

[TABLE]

where $H$ and $K$ are both functions of $r$. The perturbed Einstein equation gives a differential equation for $H$ Hinderer (2008):

[TABLE]

In vacuum, $H$ can be written as a linear combination of $P_{2}^{2}(r/M-1)$ and $Q_{2}^{2}(r/M-1)$, where $P_{2}^{2}$ and $Q_{2}^{2}$ are the $\ell=m=2$ associated Legendre functions. When expanded in powers of $M/r$ at infinity, $P_{2}^{2}(r/M-1)=\mathcal{O}(M/r)^{3}$ and $Q_{2}^{2}(r/M-1)=\mathcal{O}(r/M)^{2}$. The coefficient of $P_{2}^{2}$ is therefore related to the quadrupole moment of the star, and the coefficient of $Q_{2}^{2}$ is related to the tidal field applied by the black hole. By matching $H(r)$ and its derivative across the surface of the star, one can show Hinderer (2008)

[TABLE]

where $C=M_{\mathrm{NS}}/R$ is the compactness of the star, $Y=RH^{\prime}(R)/H(R)$, and $R$ is the radius of the star. Since $k_{2}$ depends on $Y$ and not $H$ or $H^{\prime}$ individually, Postinkov, Prakash, and Lattimer Postnikov et al. (2010) and Lindblom and Indik Lindblom and Indik (2014) define

[TABLE]

which gives rise to the first-order differential equation

[TABLE]

To find $Y=y(R)$, we numerically integrate Eq. (14) and evaluate $y$ at the surface of the star.

Despite appearances, the expression in curly braces in Eq. (13) is $\mathcal{O}(C^{5})$ due to cancellations of terms in curly brackets that are polynomial in $C$ with terms from the expansion of $\log(1-2C)$. For stars of small compactness, calculating $k_{2}$ directly from Eq. (13) is difficult because it requires that both the numerator and denominator of the right side are accurately calculated to a large number of decimal places. As a result, we expand $k_{2}$ to 20 orders in $C$. Since $k_{2}$ is $\mathcal{O}(C^{0})$, this allows for much more accurate results for small $C$. The compactness has a maximum value of $1/2$, so this expansion converges for all stars.

II.4 Estimating the Gravitational Wave Phase Shift due to Tidal Deformability

The tidal deformability $\lambda$ defined in the last section accurately describes the actual deformation of a neutron star in a binary system only when the neutron star is far from the other compact object. This is for several reasons: As the neutron star approaches the other object, linear perturbation theory and the assumption of a static spacetime used to define $\lambda$ break down; higher-order multipoles in the metric become important; as the star spirals in, its orbital angular velocity increases and becomes comparable to the frequencies of the star’s normal modes, and this enhances the star’s response to the tidal perturbation Hinderer et al. (2016); and, finally, the neutron star may be tidally disrupted before merger. Nevertheless, the tidal deformability turns out essentially to determine the departure of gravitational waveforms from those spinless BBH inspiral in numerical simulations Lackey et al. (2014, 2012); Read et al. (2013).

A post-Newtonian expansion Vines et al. (2011) describes the effect of tidal deformability on the phase of the gravitational waveform to linear order in $\lambda$:

[TABLE]

where $\Delta\Phi$ is the difference in gravitational wave phase between a spinless BBH and a BHNS binary, $\Lambda=\lambda/M_{\mathrm{NS}}^{5}$ is the dimensionless tidal deformability, $M=M_{\mathrm{BH}}+M_{\mathrm{NS}}$ is the total mass of the binary system, $\eta=M_{\mathrm{BH}}M_{\mathrm{NS}}/M^{2}$ is the symmetric mass ratio, $f$ is the linear frequency of the gravitational radiation, and $a_{0}$ and $a_{1}$ are functions of $\eta$:

[TABLE]

Where Eq. (15) is valid, in the early inspiral when the frequency $f$ is low, it allows us to easily compute the phase change (the amplitude of the waveform is also affected by tidal deformability, but in this regime the difference in amplitude is small). However, tidal effects are largest during late inspiral when the frequency is high.

To extend the analytic computation to late inspiral, Lackey et al. Lackey et al. (2014) fit the amplitude and phase of the gravitational waveforms of neutron star-black hole inspirals to the results of numerical simulations, for black hole spins $\chi_{\mathrm{BH}}$ between -.5 and .75, and mass ratio $M_{\mathrm{BH}}/M_{\mathrm{NS}}$ in the range 2 to 5. The resulting expressions (below) depend on post-Newtonian theory for low frequencies, when the neutron star is still far from the black hole. At high frequencies, the fits to numerical results take over:

[TABLE]

[TABLE]

Here the subscript $PN$ indicates the corresponding result from post-Newtonian theory; $B$ is a function of $\Lambda$, $\eta$, and $\chi_{\mathrm{BH}}$; and $E$ is a function of $\eta$ and $\chi_{\mathrm{BH}}$. The parameters of $B$ and $E$ were determined by the numerical fit. In particular,

[TABLE]

with $\{b_{0},\ b_{1},\ b_{2}\}=\{-64.985,\ \!-2521.8,\ \!555.17\}$ and $\{c_{0},\ c_{1},\ c_{2}\}=\{-8.8093,\ \!30.533,\ \!0.6496\}$ as the fitting parameters. Similarly,

[TABLE]

with

$\{g_{0},\ g_{1},\ g_{2},\ g_{3}\}$ = $\{-1.9051,\ \!15.564,\ \!-0.41109,\ \!5.7044\}$. While high tidal deformabilities increase $|\Delta\Phi|$ relative to a point-particle waveform at a given frequency $f$, they also cause stars to be tidally disrupted earlier in the inspiral, damping the resulting gravitational waves. We define the cutoff frequency $f_{\rm cutoff}$ to be the frequency at which effects from tidal deformation dampen the amplitude by a factor of $e$ relative to the post-Newtonian waveforms. To estimate the total effect of tidal deformability on the phase of the waveform throughout the inspiral, we chose to evaluate $\Delta\Phi$ at $f_{\rm cutoff}$.

The errors in the fitting parameters reported in Lackey et al. (2014) correspond to errors in $\Delta\Phi(Mf_{\mathrm{cutoff}})$ of $\sim\!\!15\%$ for typical binary parameters. The $\Delta\Phi$-values reported below should therefore not be taken as accurate predictions of the tidally-induced phase shift. Still, we expect that applying this fit to the matched causal EOS yields an upper limit on $|\Delta\Phi|$ with roughly the same error, especially considering the emphasis in Lackey et al. (2014) on avoiding over-fitting and the lower errors reported for larger $\Lambda$-values. A more accurate calculation of the phase shift from BHNS or BNS systems with our causal EOS requires numerical simulations (now in progress for BNS systems Markakis et al. (tion)) or use of the EOB formalism.

III Results

Most neutron stars observed by gravitational waves in binary inspiral are likely to have masses in or near the $1.25\ M_{\odot}$ to $1.45\ M_{\odot}$ range seen in binary neutron star systems, a range consistent with formation from an initial binary of two high-mass stars. We will see that the causal limit on the dimensionless deformability $\Lambda$ is a monotonically decreasing function of $M$ and is therefore more stringent for higher mass stars. On the other hand, the fraction of matter above nuclear density is smaller in a low-mass neutron star, and that fact limits the effect of a causal EOS above nuclear density. The net result is that the limit on $\Lambda$ set by causality is close to the values of $\Lambda$ associated with candidate neutron-star EOSs for matching densities near nuclear density.

III.1 Effect of Matching Density on Constraints

To understand the results we present in this section, it is helpful first to consider models for which the causal form (2) extends to the surface of the star, where $p=0$. (Here we follow Brecher and Caporaso Brecher and Caporaso (1976) and Lattimer Lattimer (2012).) That is, we consider models based on the EOS

[TABLE]

and having finite energy density $\epsilon_{S}$ at the surface of the star. Because the only dimensionful constant is $\epsilon_{S}$, having (in gravitational units) dimension length-2, and the mass $M$ and radius $R$ each have dimension length, we have the exact relations

[TABLE]

where $R_{\rm max}$ is the maximum radius among models with central density greater than $\epsilon_{\rm nuc}$ (low-density models have larger radii). Because the deformability $\lambda$ has dimension ${\it length}^{5}$, we similarly have

[TABLE]

Using this truncated causal EOS is equivalent to taking $\epsilon_{\mathrm{S}}=\epsilon_{\rm match}$ in the matched causal EOS and discarding the envelope of the star below $\epsilon_{\rm match}$.

We emphasize that the truncated EOS (18) is used only heuristically, to explain the near power-law dependence on $\epsilon_{\rm match}$ of the maximum mass, radius, and deformability. (The exact dependence of the maximum mass, radius, and deformability on $\epsilon_{\mathrm{match}}$ is reported below.) Because the truncated EOS sets the pressure to zero below $\epsilon_{\rm match}$, it underestimates the maximum radius and deformability. As noted earlier, to obtain a conservatively large upper limit on maximum deformability, we use the matched causal EOS (3), which has a stiff candidate EOS for $\epsilon<\epsilon_{\rm match}$.333 There is something paradoxical in using the truncated EOS (18) as an approximation to the EOS that gives the largest possible neutron stars: As Lattimer Lattimer (2012) points out (following Koranda et al. Koranda et al. (1997)), this same EOS gives maximally compact neutron stars, stars with the smallest possible radius for a given mass, among all EOSs consistent with a maximum mass at or above a largest observed value, $M_{\rm observed}$. In these and other papers Haensel and Zdunik (1989); Lattimer and Prakash (2011), Eq. (18) is chosen so that the softest possible EOS (namely $p=0$) is used up to high density; the stiff causal EOS above that density then allows $M_{\rm max}\geq M_{\rm observed}$. The resolution is this: For a fixed maximum mass, Eq. (18) yields neutron stars with the smallest possible radii. On the other hand, for a fixed $\epsilon_{\mathrm{match}}$ (i.e. for a given density up to which we assume a known EOS), Eq. (3) gives neutron stars with the largest possible radii; and, for low matching and surface densities, the difference between the matched causal EOS (3) and the truncated EOS (18) becomes negligible. Equivalently, $\epsilon_{\mathrm{match}}\rightarrow 0$ corresponds to $M_{\mathrm{max}}\rightarrow\infty$, and the difference between the softest and the stiffest possible EOSs vanishes as $M_{\mathrm{m}ax}\rightarrow\infty$. Physically, this happens because a stiff EOS is required to support large masses.

For $\epsilon_{\rm match}\lesssim\epsilon_{\rm nuc}$, where

[TABLE]

is nuclear saturation density (the central density of large nuclei), the contribution of the envelope to mass and radius is small enough that the dependence on $\epsilon_{\rm match}$ is very nearly the dependence on $\epsilon_{\mathrm{S}}$ in the truncated star: $M_{\rm max}$ and $R_{\rm max}$ are each nearly proportional to $\epsilon_{\rm match}^{-1/2}$, and $\lambda_{\rm max}$ is nearly proportional to $\epsilon_{\rm match}^{-5/2}$, where $M_{\rm max}$ is the maximum neutron-star mass consistent with causality and with a low density EOS below $\epsilon_{\rm match}$; and $R_{\rm max}$ and $\lambda_{\rm max}$ are again the corresponding maximum radius and deformability among models with central density greater than $\epsilon_{\rm nuc}$. This behavior can be seen in Figs. \IfInteger01a\StrGobbleRight01(\StrRight01) and \IfInteger01b\StrGobbleRight01(\StrRight01), where linear least-squares fits to the leftmost 10 data points in each plot satisfy

[TABLE]

The rightmost data points in each plot diverge from the line because, at higher matching densities, a larger envelope obeys the low-density (MS1) EOS.

Of greater astrophysical relevance than the upper limit on $\lambda$, however, is the constraint on the dimensionless tidal deformability, $\Lambda=\lambda/M_{\mathrm{NS}}^{5}=\frac{2}{3}k_{2}R^{5}/M^{5}$, that governs the waveform of a binary inspiral. As we will see below, because of the factor $M_{\mathrm{NS}}^{-5}$, $\Lambda$ is monotonically decreasing with increasing mass for central density above $\epsilon_{\rm match}$. The physically interesting constraint on $\Lambda$ is then a constraint at known mass: Inspiral waveforms detected with a high enough signal-to-noise ratio to measure their tidal departure from point-particle inspiral will also have the most accurately measured neutron-star masses. The dependence of $\Lambda$ on $\epsilon_{\rm match}$ for fixed mass cannot be found from the previous dimensional analysis, but it is easy to see that $\Lambda(M,\epsilon_{\rm match})$ is a monotonically decreasing function of $\epsilon_{\rm match}$: As $\epsilon_{\rm match}$ increases and less of the star is governed by the stiffer causal EOS, the star becomes more compact: $R$ decreases at fixed $M$. In addition, as the density profile becomes more centrally condensed, the tidal Love number $k_{2}$ decreases, because, for a given radius, the external tidal force has less effect on a more centrally condensed star. Decreasing $R$ and $k_{2}$ gives a sharp decrease in $\Lambda$, as shown in \IfInteger0Fig. 2Fig. \StrGobbleRight01(\StrRight01) for a $1.4M_{\odot}$ star. For $.33\epsilon_{\rm nuc}<\epsilon_{\rm match}<1.2\epsilon_{\rm nuc}$ we find a near power-law dependence,

[TABLE]

III.2 Comparison between Constraint and Results from Candidate EOSs

We begin by displaying the limit set by causality on the dimensionful tidal deformability $\lambda$ as a function of mass, with $\epsilon_{\rm match}$ taken to be $\epsilon_{\rm nuc}$. There is remaining uncertainty in the equation of state at $\epsilon_{\rm nuc}$, and we obtain a conservative upper limit by matching to the MS1 EOS Mueller and Serot (1996), which is particularly stiff for $\epsilon\lesssim\epsilon_{\mathrm{nuc}}$.

The mass-radius relation for the family of neutron stars obeying the matched causal EOS is indicated by “Matched Causal” in Fig. \IfInteger03a\StrGobbleRight01(\StrRight01). As we saw in Eqs. (22), matching to MS1 below $\epsilon_{\rm nuc}$ is a weak constraint, giving $M_{\rm max}=4.1M_{\odot}$ and $R_{\rm max}>18\ \mathrm{km}$, both significantly larger than their values for any of the candidate EOSs shown. These candidate EOSs include SLy Douchin and Haensel (2001), which is one of the softest EOSs that allow for $2M_{\odot}$ neutron stars, MPA1 Müther et al. (1987), which is slightly stiffer, H4 Lackey et al. (2006), which is stiff at low densities and soft at high densities, and MS1 Mueller and Serot (1996), which is particularly stiff at all densities. The maximum masses allowed by these EOSs are all between $2$ and $2.8\ M_{\odot}$, and the radii are all between $10$ and $15\ \mathrm{km}$.

In Figure \IfInteger03b\StrGobbleRight01(\StrRight01), the top curve displays an upper limit on $\lambda$ as a function of neutron-star mass obtained from the matched causal EOS. The comparison $\lambda(M)$ curves for the same candidate EOSs of Fig. \IfInteger03a\StrGobbleRight01(\StrRight01) show the decreasing deformability associated with stars of decreasing stiffness and radius. Note, however, that the maximum value of $\lambda$ for each EOS occurs at a smaller mass than that of the model with maximum radius. This is due to the increase in central condensation as the mass increases, resulting in an decrease in $k_{2}$. The maximum of the $\lambda(M)$ curve for the matched causal EOS gives the mass-independent upper limit $\lambda<1.5\times 10^{37}\ \mathrm{g\ cm^{2}\ s^{2}}$, for $\epsilon_{\rm match}=\epsilon_{\rm nuc}$, with dependence on $\epsilon_{\rm match}$ given by Eq. (22c) for smaller matching density.

The corresponding upper limit $\Lambda_{\rm max}(M)$ on the dimensionless deformability is given by the top curve in \IfInteger0Fig. 4Fig. \StrGobbleRight01(\StrRight01), for $\epsilon_{\rm match}=\epsilon_{\rm nuc}$. (The dependence on $\epsilon_{\rm match}$ was shown in Fig. 2 for a representative $1.4M_{\odot}$ star.) Since $\Lambda\propto C^{-5}$, $\Lambda$ is large for small masses and relatively small for larger masses. As a result, it is not meaningful to speak of a mass-independent maximum of $\Lambda$, but it is meaningful to compare $\Lambda$-values at constant mass. The most striking feature of \IfInteger0Fig. 4Fig. \StrGobbleRight01(\StrRight01) is how close the curve $\Lambda_{\rm max}(M)$ is to the range of $\Lambda$ allowed by current candidate EOSs. This stringent constraint on $\Lambda$ is in sharp contrast to the larger departures of the curves giving $R_{\rm max}(M)$ and $\lambda_{\rm max}(M)$ in \IfInteger0Fig. 3Fig. \StrGobbleRight01(\StrRight01) from the corresponding curves for candidate EOSs. For $1.4\ M_{\odot}$ stars, for example, it places the constraint that $\Lambda\leq 2300$. For comparison, $1.4\ M_{\odot}$ stars resulting from the SLy, MPA1, H4, and MS1 EOSs have $\Lambda$-values of 300, 490, 900, and 1400, respectively.

One might naively expect $|\Delta\Phi(f_{\rm cutoff})|$ to increase monotonically with $\Lambda$ and therefore to decrease monotonically with the mass $M_{\mathrm{NS}}$ of the neutron star (note that, although $\Delta\Phi$ is positive when evaluated at a given time, it is negative when evaluated at a given frequency). This is not the case, because while $|\Delta\Phi|$ increases with $\Lambda$ (and decreases with $M_{\mathrm{NS}}$) when evaluated at a fixed frequency, $f_{\rm cutoff}$ decreases monotonically with $\Lambda$ (and increases monotonically with $M_{\mathrm{NS}}$). That is, stars with high dimensionless tidal deformability are tidally disrupted at a larger distance from the black hole, corresponding to a smaller orbital (and gravitational wave) frequency. A neutron star with high tidal deformability therefore has fewer cylces during which to accumulate phase relative to a point-particle. As a result, the effect of $M_{\mathrm{NS}}$ on $|\Delta\Phi(f_{\rm cutoff})|$ is complicated, and depends on EOS and the parameters of the binary.

Nevertheless, stiffer EOSs result in larger values of $|\Delta\Phi|$ for given neutron star masses or mass ratios. As can be seen in \IfInteger0Fig. 5Fig. \StrGobbleRight01(\StrRight01), $|\Delta\Phi(f_{\rm cutoff})|$ decreases with mass ratio for all EOSs. On the other hand, $|\Delta\Phi(f_{\rm cutoff})|$ has complicated behavior with respect to neutron star mass for all EOSs when the spin of the companion black hole is zero (\IfInteger0Fig. 6Fig. \StrGobbleRight01(\StrRight01)). In addition, one can see in \IfInteger0Fig. 5Fig. \StrGobbleRight01(\StrRight01) and \IfInteger0Fig. 6Fig. \StrGobbleRight01(\StrRight01) that $|\Delta\Phi(f_{\rm cutoff})|$ increases with the stiffness of the EOS, and is largest for our EOS, but only by a few radians at most. Here, based on our estimate of $\Delta\Phi$, the constraint set by causality is remarkably strong, stronger than the already stringent constraint on $\Lambda$: $\Delta\Phi_{\rm max}(M)$ differs from its value for the stiffest candidate equation of state by less than 14%. The strength of the causal constraint is due to (a) the fact that $\Lambda$ is largest at small mass, where the causal EOS governs the smallest fraction of the star, and (b) a smaller cutoff frequency for the stiffest EOSs that reduces the time over which the phase can accumulate.

As shown in \IfInteger0Fig. 7Fig. \StrGobbleRight01(\StrRight01), for a given black hole mass $M_{\mathrm{BH}}$ and zero black hole spin $\chi_{\mathrm{BH}}$, $|\Delta\Phi(f_{\rm cutoff})|$ increases with $M_{\mathrm{NS}}$ for the Matched Causal EOS. In addition, for a given $M_{\mathrm{NS}}$, $|\Delta\Phi(f_{\rm cutoff})|$ decreases with $M_{\mathrm{BH}}$. Changing $\chi_{\mathrm{BH}}$ can change the qualitative behavior of $|\Delta\Phi(f_{\rm cutoff})|$, as can be seen in \IfInteger0Fig. 8Fig. \StrGobbleRight01(\StrRight01). In particular, a corotating companion black hole tends to make $|\Delta\Phi(f_{\rm cutoff})|$ increase with mass, while antirotating companions tend to make $|\Delta\Phi(f_{\rm cutoff})|$ decrease with mass. For a given $M_{\mathrm{NS}}$, higher (corotating) spins result in smaller $|\Delta\Phi(f_{\rm cutoff})|$, but the effect decreases with increasing $M_{\mathrm{NS}}$.

Figure \IfInteger09\StrGobbleRight01(\StrRight01) shows how $|\Delta\Phi(f_{\rm cutoff})|$ varies with mass ratio for several neutron star masses and 0 black hole spin. For a given $M_{\mathrm{NS}}$, $|\Delta\Phi(f_{\rm cutoff})|$ decreases with increasing mass ratio. For a given mass ratio, $|\Delta\Phi(f_{\rm cutoff})|$ decreases with neutron star mass. The effect decreases in magnitude as the mass ratio increases.

Finally, \IfInteger0Fig. 10Fig. \StrGobbleRight01(\StrRight01) shows $|\Delta\Phi(f_{\rm cutoff})|$ varying with mass ratio for several black hole spins and $M_{\mathrm{NS}}=1.4\ M_{\odot}$. $|\Delta\Phi(f_{\rm cutoff})|$ decreases with mass ratio regardless of the value of $\chi_{\mathrm{BH}}$, but for a given mass ratio, $|\Delta\Phi(f_{\rm cutoff})|$ decreases with $\chi_{\mathrm{BH}}$; it is smallest for corotating black holes, and largest for antirotating black holes.

IV Conclusion

By using a stiffest causal EOS consistent with causality at high density, matched to the MS1 Mueller and Serot (1996) EOS below a density $\epsilon_{\rm match}$, we have set upper limits on the quadrupole tidal deformability $\lambda$ and on the dimensionless tidal deformability $\Lambda$ as a function of neutron star mass. The limit on $\Lambda$, given by Eq. (23) for a 1.4 $M_{\odot}$ neutron star, is conservative, because we have matched to an EOS (MS1) that is stiff below nuclear density: With this low-density EOS and a match at $\epsilon_{\rm nuc}$, the corresponding upper mass limit is $4.1M_{\odot}$. Using the constraint on dimensionless tidal deformability and the Lackey et al. analytic fit to numerical data Lackey et al. (2014), we then estimated the induced phase shift of a BHNS inspiral and merger waveform.

The implied upper limit on the accumulated phase shift $|\Delta\Phi|$ at merger depends on the parameters of the binary, but it is surprisingly close to the range of phase shifts seen in candidate EOSs. Assuming one can neglect resonant interactions of the tidal field with neutron-star modes, we think this conclusion is secure. We emphasize, however, that our upper limits on $|\Delta\Phi|$ rely on an analytic expression based on full numerical simulations for models with a set of EOSs significantly less stiff than the matched causal EOS. Work has begun on numerical simulations to obtain an upper limit on the departure of double neutron star inspiral waveforms from the point-particle (or spinless BBH) case.

Acknowledgements.

We would like to thank Benjamin Lackey and Lee Lindblom for sharing data with us, Lindblom for a useful discussion regarding multiple-parameter EOSs, and James Lattimer for helpful comments on the paper. This work was supported in part by NSF Grant Nos. PHY-1307429 and PHY-1607585.

Appendix A Comments on causality

With the assumption that the equilibrium equation of state of the neutron star and its perturbations are governed by the same one-parameter equation of state, causality implies $dp/d\epsilon<1$. That is, as mentioned in the text, the time-evolution of a barotropic fluid is described by a hyperbolic system whose characteristics lie within the light cone precisely when $dp/d\epsilon<1$ Lichnerowicz (1965). The frequencies of stellar perturbations, however, are too high for the temperature of a fluid element and the relative density $Y_{i}$ of each species of particle to reach their values for the background fluid at the same pressure: Heat flow and nuclear reactions are incomplete.

Because of this, one cannot precisely identify the maximum speed of signal propagation in the fluid with the equilibrium value

[TABLE]

If short wavelength, high frequency perturbations are too rapid for heat flow and for nuclear reactions to proceed, their speed of propagation is

[TABLE]

One therefore expects causality to imply

[TABLE]

This is known to be true for a relativistic fluid with a two-parameter EOS of the form $p=p(\epsilon,s)$: Its dynamical evolution then involves heat flow and is governed by the equations of a dissipative relativistic fluid. Causal theories of this kind were first introduced by Israel and Stewart Israel (1976); Israel and Stewart (1979) and by Liu et al. Liu et al. (1986). The general class of such theories was analyzed by Geroch and Lindblom Geroch and Lindblom (1991), who pointed out that, for dissipative fluids obeying $p=p(\epsilon,s)$, causality implies the inequality (25),

[TABLE]

Now a star is unstable to convection if

[TABLE]

Thus, for a locally stable spherical star (a self-gravitating equilibrium configuration of a relativistic dissipative fluid) based on a two-parameter EOS $p=p(\epsilon,s)$, causality implies

[TABLE]

Thus, at least for two-parameter dissipative fluids, one can rule out the possibility that dispersion in a dissipative fluid could lead to a group velocity smaller than the phase velocity (see, for example Bludman and Ruderman Bludman and Ruderman (1970)) $v_{\rm sound}$ and thereby allow $v_{\rm sound}>1$ without superluminal signal propagation.

For a dissipative fluid obeying a multi-parameter EOS of the form $p=p(\epsilon,s,Y_{i})$, we are not aware of a general proof that causality implies the inequality (25). One has only the weaker statement, for a locally stable spherical star based on an EOS equation of state $p=p(\epsilon,s)$, $v_{\rm sound}<1$ implies the equilibrium inequality (28). There is one additional caveat: The core of a neutron star is likely to be a superfluid, and taking that into account could lead to small corrections in the speed of sound.

Finally, we note that for candidate EOSs, although the inequality $v_{\rm sound}<1$ is stronger than the the equilibrium inequality (28) used to place upper limits on mass, radius and, in the present paper, on deformability, the difference is small. The fractional difference

[TABLE]

is primarily due to composition (to the constant values of $Y_{i}$), and it is less than 5%. (It is approximately half the fractional difference between the adiabatic index $\gamma=\Gamma_{1}$ and the index $\Gamma$ governing the equilibrium configuration; the difference determines the Brunt-Vaïsälä frequency, a characteristic frequency of $g$-modes, and an estimate can be found, for example, in Ref. Reisenegger and Goldreich (1992).)