The sum of Love: Exploring the effective tidal deformability of neutron stars

TL;DR

This paper investigates how the composition of neutron star matter affects the tidal deformability during binary inspiral, providing insights into systematic effects on gravitational-wave signals.

Contribution

It offers the first detailed analysis of composition effects on tidal deformability, enhancing understanding of neutron star physics in gravitational-wave observations.

Findings

Composition impacts tidal deformability estimates

Systematic effects may influence gravitational-wave signal interpretation

Provides a framework for incorporating matter composition into models

Abstract

Finite size effects come into play during the late stages of neutron star binary inspiral, with the tidal deformability of the supranuclear density matter leaving an imprint on the gravitational-wave signal. As demonstrated in the case of GW170817, this leads to a constraint on the neutron star radius (and hence the equation of state). A deeper understanding of the tidal response requires an analysis of both the state and composition of matter. While these aspects may not have dramatic impact, they could lead to systematic effects that need to be kept in mind as the observational data become more precise. As a step in this direction we explore the role of the composition of matter, which is likely to remain "frozen" during the late stages of binary inspiral. We provide the first in-depth analysis of the problem, including estimates of how composition impacts on the effective tidal deformability. The results provide improved insight into how aspects of physics that tend to be "ignored" impact on binary neutron star gravitational-wave signals.

Tables

Table 1. Table 1: Numerical results for barotropic n = 1 𝑛 1 n=1 polytropes (with Γ 1 = 2 subscript Γ 1 2 \Gamma_{1}=2 ) and l = 2 𝑙 2 l=2 .

mode	${\tilde{\omega }}_n$	$\left|{\tilde{Q}}_n\right|$	${\left(V_n/W_n\right)}_R$	$k_l^n$
$p_4$	9.0471	$4.4927\times {10}^{-5}$	$1.2218\times {10}^{-2}$	$-1.7426\times {10}^{-6}$
$p_3$	7.2581	$3.0791\times {10}^{-4}$	$1.8983\times {10}^{-2}$	$-5.9238\times {10}^{-4}$
$p_2$	5.4144	$2.6168\times {10}^{-3}$	$3.4104\times {10}^{-2}$	$1.3425\times {10}^{-3}$
$p_1$	3.4615	$2.6888\times {10}^{-2}$	$8.3845\times {10}^{-2}$	$-1.6411\times {10}^{-2}$
$f$	1.2267	$5.5792\times {10}^{-1}$	$4.4163\times {10}^{-1}$	0.77528

Table 2. Table 2: Numerical results for n = 1 𝑛 1 n=1 polytropes with Γ 1 = 2.05 subscript Γ 1 2.05 \Gamma_{1}=2.05 and l = 2 𝑙 2 l=2 .

mode	${\tilde{\omega }}_n$	$\left|{\tilde{Q}}_n\right|$	${\left(V_n/W_n\right)}_R$	$k_l^n$
$p_2$	5.4949	$2.5459\times {10}^{-3}$	$3.3111\times {10}^{-2}$	$1.6522\times {10}^{-3}$
$p_1$	3.5204	$2.5865\times {10}^{-2}$	$8.1049\times {10}^{-2}$	$-1.9849\times {10}^{-2}$
$f$	1.2274	$5.5796\times {10}^{-1}$	$4.4009\times {10}^{-1}$	0.77055
$g_1$	0.1845	$1.7657\times {10}^{-3}$	27.8841	$2.2458\times {10}^{-3}$
$g_2$	0.1270	$4.264\times {10}^{-4}$	60.9918	$8.2137\times {10}^{-4}$
$g_3$	0.0975	$1.2493\times {10}^{-4}$	104.4819	$3.3549\times {10}^{-4}$
$g_4$	0.0794	$4.1194\times {10}^{-5}$	158.1560	$1.4542\times {10}^{-4}$

Table 3. Table 3: Numerical results for n = 1 𝑛 1 n=1 polytropes with Γ 1 = 7 / 3 subscript Γ 1 7 3 \Gamma_{1}=7/3 and l = 2 𝑙 2 l=2 .

mode	${\tilde{\omega }}_n^2$	$\left|{\tilde{Q}}_n\right|$	${\left(V_n/W_n\right)}_R$	$k_l^n$
$p_2$	5.9316	$2.2079\times {10}^{-3}$	$2.8417\times {10}^{-2}$	$9.6779\times {10}^{-4}$
$p_1$	3.8411	$2.1197\times {10}^{-2}$	$6.8024\times {10}^{-2}$	$-1.0631\times {10}^{-2}$
$f$	1.2291	$5.5805\times {10}^{-1}$	$4.3227\times {10}^{-1}$	0.74685
$g_1$	0.4361	$1.110\times {10}^{-2}$	4.9584	$1.4299\times {10}^{-2}$
$g_2$	0.3035	$2.6048\times {10}^{-3}$	10.6541	$5.0402\times {10}^{-3}$
$g_3$	0.2343	$7.6158\times {10}^{-4}$	18.0985	$2.0379\times {10}^{-3}$
$g_4$	0.1912	$2.4981\times {10}^{-4}$	27.2719	$8.7876\times {10}^{-4}$
$g_5$	0.1618	$8.7649\times {10}^{-5}$	38.1746	$3.9289\times {10}^{-4}$

Table 4. Table 4: The accumulated contribution to the tidal deformability from the different modes, in order of relevance of the contribution, for the three different models we consider. We expect (for the barotropic Γ = 2 Γ 2 \Gamma=2 case) to have k l ≈ 0.259909 subscript 𝑘 𝑙 0.259909 k_{l}\approx 0.259909 . If we add up the contributions from the different modes in each case, the results converges to the expected answer. The mode sum is always dominated by the f-mode with a modest correction from the other modes of the star, but the enhanced importance of the g-modes with increasing stratification is notable.

${\Gamma }_1=2$		${\Gamma }_1=2.05$		${\Gamma }_1=7/3$
mode	$k_l$	mode	$k_l$	mode	$k_l$
$f$	0.27528	$f$	0.27055	$f$	0.24685
$+p_1$	0.25887	$+p_1$	0.25526	$+g_1$	0.26115
$+p_2$	0.26021	$+p_2$	0.25653	$+p_1$	0.25052
$+p_3$	0.26015	$+g_1$	0.25878	$+g_2$	0.25556
		$+g_2$	0.25960	$+p_2$	0.25653
		$+g_3$	0.25993	$+g_3$	0.25856
		$+g_4$	0.26008	$+g_4$	0.25944
				$+g_5$	0.25983
	$9\times {10}^{-4}$		$7\times {10}^{-4}$		$3\times {10}^{-4}$

Equations

$$t_M\sim \frac{2\mathrm{months}}{T_9^6},t_D\sim \frac{20\mathrm{s}}{T_9^4},$$

$${\partial }_t\left(\Delta \rho +\rho {\nabla }_i{\xi }^i\right)=0$$

$$\Delta =\delta +{\mathcal{L}}_{\xi }$$

$$\Delta v^i=\delta v^i={\partial }_t{\xi }^i$$

$${\partial }_t^2{\xi }_i+\frac{1}{\rho }{\nabla }_i\Delta p-\frac{\Delta \rho }{{\rho }^2}{\nabla }_{i}p+{\nabla }_i\Delta \Phi =-{\nabla }_i\chi$$

$${\partial }_t^2{\xi }_i+\frac{1}{\rho }{\nabla }_i\delta p-\frac{1}{{\rho }^2}\delta \rho {\nabla }_{i}p+{\nabla }_i\delta \Phi =-{\nabla }_i\chi$$

$${\nabla }^2\delta \Phi =4\pi G\delta \rho$$

$$\chi =-\frac{G{M}^{\prime }}{|\mathbf{r}-\mathbf{D}(t)|}=-G{M}^{\prime }\sum _{l\ge 2}\sum _{m=-l}^l\frac{W_{lm}{r}^l}{D^{l+1}(t)}{Y}_{lm}{e}^{-im\psi (t)}$$

$$W_{20}=-\sqrt{\pi /5},W_{2\pm 2}=\sqrt{3\pi /10},W_{2\pm 1}=0.$$

$$\hat{\chi }(\omega ,r)={\int }_{-\infty }^{\infty }\chi (t,r)e^{i\omega t}dt=-G{M}^{\prime }\sum _{l\ge 2}\sum _{m=-l}^l{W}_{lm}{r}^l{Y}_{lm}\left[{\int }_{-\infty }^{\infty }\frac{e^{i\omega t-im\psi (t)}}{D^{l+1}(t)}dt\right]$$

$$\hat{\chi }(\omega ,r)=\sum _{l\ge 2}\sum _{m=-l}^l{v}_l{r}^l{f}_{lm}(\omega )Y_{lm}$$

$$\dot{\psi }={\left[\frac{GM}{D^3}\right]}^{1/2}=\Omega ⟶\psi \approx \Omega t$$

$$\hat{\chi }(\omega ,r)=\sum _{l\ge 2}\sum _{m=-l}^l{v}_l{r}^l\delta (\omega -m\Omega )Y_{lm}$$

$$\dot{D}=-\frac{64G^3}{5c^5}\frac{M^{4/3}{\mathcal{M}}^{5/3}}{D^3}$$

$$\mathcal{M}={\mu }^{3/5}{M}^{2/5}$$

$$t_D=\frac{D}{|\dot{D}|}\ll \frac{1}{\Omega }$$

$${\nabla }_i{\hat{\xi }}^i=0$$

$$\delta \rho =\Delta \rho -{\xi }^j{\nabla }_i\rho ⟶\delta \hat{\rho }=0$$

$$-{\omega }^2{\hat{\xi }}_i+\frac{1}{\rho }{\nabla }_i\delta \hat{p}+{\nabla }_i\delta \hat{\Phi }=-{\nabla }_i\hat{\chi }$$

$${\nabla }^2\delta \hat{\Phi }=0$$

$${\nabla }^2\delta \hat{p}=0$$

$$\delta \hat{\Phi }=\sum _l{\Phi }_l{Y}_{lm}$$

$${\Phi }_l=c{r}^l+\frac{d}{r^{l+1}}$$

$${\Phi }_l=c{r}^l$$

$$p_l=b{r}^l$$

$${\chi }_l=v{r}^l$$

$$-{\omega }^2{\xi }_l+\frac{1}{\rho }{\partial }_r{p}_l+{\partial }_r{\Phi }_l=-{\partial }_r{\chi }_l$$

$${\omega }^2{\xi }_l=\left(\frac{1}{\rho }b+c+v\right)l{r}^{l-1}$$

$$\Delta p=\delta p+{\xi }^j{\nabla }_{j}p=0⟶b{r}^l+{\xi }_l{p}^{\prime }=0\text{mbox}atr=R$$

$$p^{\prime }=-\rho {\Phi }^{\prime }=-\frac{G{M}_{\star }\rho }{R^3}r$$

Full text

I Introduction

The breakthrough detection of signals from inspiralling and merging black-hole binaries gw150914 ; catalog and the spectacular neutron star merger event GW170817 gw170817 demonstrate the promise of gravitational-wave astronomy. As the sensitivity of the detectors improves we can expect other exciting discoveries. In particular, further events involving neutron stars should allow us to put tighter constraints on the physics of matter under extreme conditions.

A typical neutron star binary system may spend as much as 15 minutes in the sensitivity band of advanced ground-based interferometers (above 10 Hz). The detection of, and extraction of parameters from, such signals is of great importance for both astrophysics and nuclear physics. From the astrophysics point of view observed event rates should lead to insight into the formation channel(s) for these systems, while the nuclear physics aspects relate to the poorly constrained equation of state for matter at supranuclear densities (see for example LOFT ). Focussing on the nuclear physics, binary neutron star signals constrain the equation of state of supranuclear matter in two ways. First of all, finite size effects impact on the inspiral signal. The deformability of the stellar fluid leaves an imprint on the late-time chirp, an effect encoded in the tidal Love numbers hind1 ; hind2 ; hind3 ; agat . Meanwhile, oscillations of the merger remnant, which depend on the hot equation of state, are expected to leave a robust signature baus1 ; baus2 ; baus3 ; bern ; luc1 ; luc2 . However, the high-frequency nature of the merger signal makes it difficult to detect with the current generation of detectors (even at design sensitivity) clark ; ligo_pm .

In this paper we focus on the tidal deformability. Our aim is simple; we want to explore to what extent the composition of the neutron star matter enters the problem. The motivation for this is clear. As the star is deformed by the tidal interaction, matter is driven out of equilibrium and it is easy to argue that the the relevant nuclear reactions are too slow to re-establish equilibrium on the timescale of inspiral111Note that the situation is different for the tide in an ordinary stellar binary, where the orbital evolution is excruciatingly slow. In that context, the assumption of a “static” tide, encoding the deformability in the Love numbers, does not need to be questioned.. This is clear from, for example, the estimates in haensel . At the simplest level, the relevant equilibration timescales are

[TABLE]

for the modified and direct Urca reaction, respectively. The temperature is scaled to hot systems, $T_{9}=T/10^{9}$ K, but inspiralling neutron stars are old and cold, typically in the range $T_{9}\leq 0.01$, which would make both $t_{M}$ and $t_{D}$ much longer than the time it takes a given system to move through the sensitivity band of a ground-based interferometer (minutes). Moreover, at the expected temperature the star’s interior should be superfluid, in which case reactions are exponentially suppressed. However, superfluidity brings additional aspects to the problem (which warrant a separate discussion) so we will focus on the (simpler) case of non-superfluid matter in the following.

The upshot of the equilibration argument is that the equation of state is no longer barotropic, as has been assumed in virtually every previous analysis of the tidal problem (most of which follow the steps laid out in hind2 ; hind3 ). Instead, we want to establish to what extent a frozen matter composition leads to a noticeable effect on, for example, the Love numbers222The issue also arises in connection with the so-called I-Love-Q relations ilq , but we will not consider that aspect here., and if this in turn impacts on the extraction of neutron star parameters from an observed signal. The effect may well be small enough that we can ignore it, given the anticipated observational “errors”, but we need to make sure that this is the case. In essence, we want to quantify the systematic “errors” associated with the current models.

II The adiabatic problem

In order to set the stage for the discussion, it is useful to remind ourselves of the context. The tide raised by a binary companion (here treated as a point particle, which should be a good enough approximation for our purposes) induces a linear response in the primary. In order to quantify this response we solve the linearised fluid equations in Newtonian gravity. Assuming that the star is non-rotating, we first of all have the perturbed continuity equation

[TABLE]

where $\xi^{i}$ is the displacement vector associated with the Lagrangian perturbation fs78a

[TABLE]

(with $\delta$ the corresponding Eulerian variation and $\mathcal{L}_{\xi}$ the Lie derivative along $\xi^{i}$) such that the perturbed velocity is given by

[TABLE]

The perturbed Euler equation is

[TABLE]

or, equivalently, provided that the background is in hydrostatic equilibrium,

[TABLE]

(where it is useful to keep in mind that the Lagrangian variation commutes with the gradient $\nabla_{i}$).

We also have the Poisson equation for the gravitational potential

[TABLE]

while the tidal potential due to the presence of the binary partner (which generates the fluid perturbation) is given by a solution to $\nabla^{2}\chi=0$.

In a coordinate system centred on the primary, which we will take to have mass $M_{\star}$, we have holai

[TABLE]

where $M^{\prime}$ is the mass of the secondary. The orbit of the companion is taken to be in the plane $[D(t),\pi/2,\psi(t)]$ where $D$ is the binary separation and $\psi$ is the orbital phase. For $l=2$ (which leads to the main contribution to the gravitational-wave signal) we have pt77

[TABLE]

The last result follows from symmetry; $W_{lm}$ must vanish for all odd $l+m$.

We need to work out the deformation of the star due to the companion’s tidal field. That is, we are looking for a solution to (6) given the specific form for $\chi$. The problem is usually explored in the frequency domain. Simplistically, this is often taken to mean that the time-dependence of perturbed quantities is harmonic, proportional to $e^{i\omega t}$. In many situations this is sufficient, but in the present case we need to be a little bit more careful. The starting point would be the text-book Fourier transform (or, perhaps rather a Laplace transform, as we are dealing with an initial-value problem, noting that this would also involve changing the lower limit in the integral below). Formally, we need

[TABLE]

using hats to indicate frequency domain quantities. In general, we have

[TABLE]

and in order to proceed we need to account for the evolution of the orbital separation $D$ and the phase $\psi$. The standard approach is to assume that the evolution is adiabatic and make use of the stationary phase approximation. In the simplest case one would progress by connecting a sequence of circular orbits with fixed separation l94 . Then we have

[TABLE]

and it follows straight away that

[TABLE]

As we are dealing with a linear problem, the delta function in frequency is inherited by all perturbed quantities and hence the solution will (at the end of the day) only have support at distinct frequencies $\omega=m\Omega$. In particular, there will be a distinction between the $m=0$ contribution, which is time independent and leads to a static tide, and the dynamical contributions from the $m\neq 0$ terms (usually discussed in terms of resonances associated with the star’s oscillation modes l94 ; ks ; wynn ).

The problem becomes more complicated when we account for the orbital evolution. The frequency support in (11) becomes less obvious and the inversion to the time domain more involved. However, this is the situation we are dealing with when we consider neutron star binaries sweeping through the sensitivity band of a gravitational-wave interferometer. To leading order, the evolution timescale due to gravitational radiation reaction is given by

[TABLE]

where the dot is a time derivative and we have introduced the chirp mass

[TABLE]

with the total mass $M=M_{\star}+M^{\prime}$ and the reduced mass $\mu=M_{\star}M^{\prime}/M$, as usual. The adiabatic approximation should be reasonable as long as

[TABLE]

and as long as we work at this level we can make progress in accounting for the shrinkage of the orbit (see, for example, the discussion in l94 ; ks ). However, in order to reach the precision required for gravitational-wave searches (and parameter extraction) we need to do better. The problem is messy as the orbital evolution requires a higher order post-Newtonian analysis. Moreover, this may still not be sufficient, as the post-Newtonian scheme breaks down as the system approaches merger, leading to nonlinear simulations becoming necessary. Attempts to match approximate solutions and simulation results set the (still developing) state of the art bern2 ; diet2 ; kawa ; diet3 ; diet4 ; fouc . However, the main point we need to appreciate is simple: A realistic model will not lead to (13). Rather, we end up with (11) which means that we have to account for a frequency dependent factor ($f$) in order to complete a dynamical description of the tidal response. This is a key part of the problem, but we will not attempt to solve it here. Instead, we take the view that the relevant frequency dependence (i.e. the function $f(\omega$)) is “known” and focus on the response of the stellar fluid. Admittedly, this leaves the model incomplete but we nevertheless believe that our analysis adds valuable insight.

III A simple homogeneous model

Let us turn to the issue of solving (6) in the frequency domain. As a first example, we consider a homogeneous, incompressible star. The calculation is straightforward and turns out to be instructive. We have $\Delta\hat{\rho}=0$ and $\nabla_{i}\rho=0$, which means that the continuity equation reduces to

[TABLE]

Moreover

[TABLE]

The Euler equation then simplifies to

[TABLE]

while the Poisson equation becomes

[TABLE]

Since we also know that $\nabla^{2}\hat{\chi}=0$, we see from (19) that

[TABLE]

Expanding in spherical harmonics (here, and in the following, leaving out the implied, frequency dependent, factor of $f$ from (11) in all perturbation expressions), it is easy to see that the radial problem is degenerate (explicitly depending only on $l$, with the source term determining the $m$ dependence). Introducing notation such that (omitting the $m$ label on the different quantities for clarity)

[TABLE]

and similar for other perturbed quantities, the solution to the radial part of Laplace’s equation for a given value of $l$ can be written (suppressing the indices on the coefficients, which should not cause confusion)

[TABLE]

That is, the (regular) interior solution is

[TABLE]

Similarly, we have

[TABLE]

while (11) leads to

[TABLE]

The radial component of (19) leads to (after using the orthogonality of the spherical harmonics)

[TABLE]

(here and in the following steps we take $\xi_{l}$ to represent the radial component of the displacement) or

[TABLE]

We also know that the Lagrangian variation of the pressure must vanish at the surface, so

[TABLE]

Using hydrostatic equilibrium

[TABLE]

(with primes representing derivatives with respect to $r$) we have

[TABLE]

We also require continuity of the perturbed gravitational potential across the surface. This leads to

[TABLE]

Meanwhile, the derivative of the potential must satisfy (remembering that the density is not continuous at the surface in the case under consideration)

[TABLE]

Combining the two conditions, we see that the interior potential should satisfy

[TABLE]

That is, we have

[TABLE]

In the absence of a tidal interaction, the problem reduces to that of free oscillations of the body. That is, the solution should provide the normal modes of the star. Thus, setting $v=0$ we can sanity check the calculation against the standard f-mode result for uniform density stars (see, for example, nagw ). From (28) we get

[TABLE]

What changes when we add the tidal potential? In essence, we need $\delta\hat{\Phi}\to\delta\hat{\Phi}+\hat{\chi}$ as the total potential and its derivative have to be continuous at the surface. However, the tidal potential is already continuous, so the condition (34) remains unchanged. This means that we have the two relations

[TABLE]

(from before) and

[TABLE]

Combining these, we arrive at

[TABLE]

Introducing the dimensionless frequency

[TABLE]

and making use of the f-mode result, we arrive at the final relation

[TABLE]

Let us now make contact with the Love number and the tidal deformability. In general, the matching of the gravitational potential at the star’s surface provides the multipole moment, $I_{l}$, of the body according to

[TABLE]

If we also use

[TABLE]

then the Love numbers $k_{l}$ are defined as

[TABLE]

For the homogeneous model we then have (cf. the discussion in ogilvie )

[TABLE]

The model is admittedly simplistic, but it has led us to a useful result. The final relation (45) links the tidal deformability, expressed in terms of the Love number, to the frequency of the fundamental mode of the star. We can now make a number of observations: First of all, in the limit $\tilde{\omega}\ll\tilde{\omega}_{f}$ we have

[TABLE]

This is the expected result for the “equilibrium” tide, as we arrive at the same answer by setting $\omega=0$ from the outset willpo . Secondly, we see how the Love number depends on the f-mode frequency. This helps explain phenomenological relations like those developed in chan . Next, we note that, in the case of a fixed orbital distance, the time-domain result splits into the expected static contribution (from $m=0$) and a time dependent component with the f-mode resonances at $\omega_{f}=\pm m\Omega$ (from the $m\neq 0$ terms) holai ; l94 ; ks ; wynn . However, the difference is entirely due to the frequency dependence of $f(\omega)$ which comes into play when we invert the transform to the time domain. Finally, and perhaps most notably, the relation (45) provides the frequency dependence of the star’s tidal response. This is a key ingredient for a model involving evolving orbits and, as we will see, crucial for an understanding of the role the matter composition plays in the problem.

IV The effective Love number

Inspired by the homogeneous model, let us now consider the problem for more realistic stellar models. In essence, we need to solve

[TABLE]

for some given equation of state.

First of all, recalling the usual argument for the equilibrium tide, we assume that the matter is in chemical equilibrium (we will relax this assumption later). That is, we take the equation of state is to be barotropic, $p=p(\rho)$, such that

[TABLE]

with $c_{s}^{2}$ the speed of sound (calculated for matter in beta equilibrium). We also know that the unperturbed background configuration is such that

[TABLE]

where (as before) a prime indicates a radial derivative and we have introduced the gravitational acceleration $g$.

Expanding in spherical harmonics (as before) and introducing

[TABLE]

we have the radial component of the Euler equation (in the $\omega\to 0$ limit)

[TABLE]

and the angular part

[TABLE]

Meanwhile, the perturbed Poisson equation becomes

[TABLE]

where we have made use of the fact that $\chi_{l}$ solves the corresponding homogeneous equation.

Let us now pause to note that the problem appears to be over determined. We seem to have too many equations for the number of variables. However, taking a radial derivative of (52) we get

[TABLE]

Using this in (51) we have

[TABLE]

which is consistent with (48) as long as we have a barotropic equation of state. This identity reduces the number of equations, so the problem is well posed, after all. In essence, this is the Newtonian version of the result discussed in penner .

Now we have

[TABLE]

which is easily solved by integrating from the centre to the surface of the star. At the surface we match to the exterior potential. This matching provides the multipole moment, $I_{l}$, and the Love number, $k_{l}$, reproducing the steps from (42) to (44).

However, as we have already indicated, the calculation of the tidal deformability and the Love number may not be (strictly) valid for realistic neutron star binaries close to merger. Basically, the nuclear reactions required to establish chemical equilibrium are too slow to act on the inspiral timescale. If this is the case we can no longer assume that the equation of state for the perturbations is barotropic. Instead, it would be reasonable to assume that the composition is held frozen as the system sweeps through the sensitivity band of a ground based detector. This changes the response of the stellar fluid to the tidal driving which, in turn, allows us to estimate the impact the matter composition has on the problem.

Taking the matter composition into account (assuming a simple model composed of neutrons, protons and electrons), we have a two-parameter equation of state $p=p(\rho,x_{\rm p})$ (say), where $x_{\rm p}=n_{\rm p}/n$ is the proton fraction. In the Newtonian context the mass density is simply $\rho=m_{B}n$ where $n$ is the baryon number density. Hence, we can think of $\rho$ as a proxy for the number density. Moreover, the continuity equation (2) remains unchanged (although, strictly speaking, it now represents baryon number conservation). As an illustration of the range of values we may need to consider, and the corresponding variation with density, we illustrate the proton fraction for three equation of state parametrisations from the Brussel-Montreal collaboration Bsk1 ; Bsk2 in figure 1. The BSk models have the particular advantage that one can readily work out the different thermodynamical derivatives we require.

In order to account for nuclear reactions, we follow gmode and introduce a new dependent variable $\beta=\mu_{\rm n}-\mu_{\rm p}-\mu_{\rm e}$ which encodes the deviation from chemical equilibrium (with $\mu_{\rm x}$, ${\rm x}={\rm n},{\rm p},{\rm e}$ being the chemical potentials). We then have

[TABLE]

with

[TABLE]

and the thermodynamical derivatives

[TABLE]

The relevant reaction rate is encoded in $\gamma$, assuming a small deviation from equilibrium.

Let us consider the relevant timescales. Introducing a characteristic reaction time

[TABLE]

we see that, if the reactions are fast compared to the dynamics (assumed to take place on a timescale $\sim 1/\omega$) then we have $\omega t_{R}\ll 1$ and if we consider the evolution in the adiabatic limit (where $\omega\to 0$) then

[TABLE]

The system remains in beta-equilibrium and the standard (barotropic) analysis of the tidal deformability holds.

However, we are interested in the opposite limit – where reactions are slow. Then we have $\omega t_{R}\gg 1$ (see gmode ) and we can Taylor expand (62) to get

[TABLE]

However, we can no longer (at least not meaningfully) take the $\omega\to 0$ limit. If we insist on doing this, the problem inevitably becomes barotropic as, for any fixed $t_{R}$, we must cross over into the fast reaction regime. The argument is analogous to that for the g-modes, discussed in gmode . In essence, if we want to consider the impact of a frozen composition we cannot work in the static limit – we have to consider dynamical aspects of the tide and solve the frequency dependent problem.

Motivated by this conclusion, let us see how the analysis changes if we consider a general (compressible and possibly stratified) model. First of all, we need to quantify the difference between the barotropic and the frozen-composition cases. Effectively, this can be done by noting that we have

[TABLE]

with

[TABLE]

in the barotropic case, but with $\Gamma_{1}$ distinct from $\Gamma$, in general. An illustration, again for three BSk models, of the difference between the two adiabatic indices (in the two limits we have discussed) is provided in figure 2. Meanwhile, the relative difference between $\Gamma_{1}$ and $\Gamma$ is shown in figure 3. In the numerical examples considered later, we will take $\Gamma_{1}$ to be constant, with values spanning the range suggested by figure 2. This is, perhaps, not particularly realistic, but it serves as a useful starting point.

IV.1 Mode expansion

Taking a lead from the results for the homogeneous model, we aim to express the driven response of the fluid in terms of a set of normal modes333Note that, while it is common to assume that the modes form a complete set, it is not clear that this is actually the case for more realistic neutron star models., corresponding to solutions $\xi_{n}$ (where $n$ is a label that identifies the modes, say in terms of the number of nodes in the radial eigenfunction and the corresponding spherical harmonics). Letting the (real) mode-frequency be $\omega_{n}$ we have (leaving out the hats for frequency domain quantities in order to avoid the notation becoming cluttered)

[TABLE]

and each individual mode (labelled by $n$) satisfies

[TABLE]

where the $C$ operator is messy but we do not need an explicit expression here.

Making use of the inner product from fs78b

[TABLE]

(where the asterisk indicates the complex conjugate) it is easy to show that the modes are orthogonal (at least as long as the frequencies are real). That is, we have (keeping the normalisation of the modes explicit for the moment)

[TABLE]

We use the orthogonality to rewrite (47) as an equation for the mode amplitudes:

[TABLE]

Finally, making use of the perturbed continuity equation

[TABLE]

and integrating by parts, we have (assuming that the density vanishes at the surface of the star)

[TABLE]

In general, e.g. when the star is spinning, it may be practical to express the stellar perturbations with respect to a different set of spherical harmonics holai , but we will not worry about this here.

Making use of the given expression for the tidal potential (11) we have an equation for the driven modes (for each $l$, as the different values of $m$ are still degenerate)

[TABLE]

where we have introduced the “overlap integral”

[TABLE]

In effect, we have a driven set of modes with amplitude

[TABLE]

which may become resonant during a binary inspiral. This is a well-known result l94 ; ks ; pt77 .

Turning to the matching at the surface, the perturbed gravitational potential satisfies

[TABLE]

For the overlap integral we need

[TABLE]

Integrating by parts, we get

[TABLE]

However, we know that the solution should satisfy (since the density now vanishes as $r\to R$)

[TABLE]

so (77) leads to

[TABLE]

and we have

[TABLE]

We recognize $I_{n}$ as the contribution each mode makes to the mass multipole moment.

IV.2 The effective tidal deformability

We now want to connect the mode expansion to the tidal deformability and the effective Love number. The purpose is to discuss what happens far away (above and below) from a given resonance. In order to do this, we need a representation of the perturbed gravitational potential associated with the mode expansion.

In order to obtain the desired result, we need to connect the mode expansion for the displacement to the gravitational potential. It is then useful to consider the components of the displacement vector

[TABLE]

and similarly for the contribution from each mode

[TABLE]

from which it follows that the normalisation constant is given by

[TABLE]

We now see that the tidal problem leads to a fluid displacement of form

[TABLE]

Going back to the Euler equation, the $\theta$-component leads to

[TABLE]

Moreover, at the surface we have $\Delta p=0$ so

[TABLE]

and it follows that

[TABLE]

That is, we have

[TABLE]

Making use of the mode expansion, this becomes

[TABLE]

In the low-frequency limit (with $\omega$ much smaller than the lowest frequency mode444This is a somewhat subtle issue, but the discussion of the g-mode spectrum in gmode ensures that the limiting procedure makes sense.), we have

[TABLE]

where we have used the dimensionless frequency from before.

Now recall that

[TABLE]

and use

[TABLE]

to get

[TABLE]

This leads to the final result

[TABLE]

Alternatively, we can rewrite (93) as

[TABLE]

Introducing the dimensionless overlap integral ks

[TABLE]

we have

[TABLE]

and

[TABLE]

and if we (finally) normalise the modes in such a way that

[TABLE]

we arrive at the expression

[TABLE]

As a quantitative test of (100) we compare results for three models corresponding to a background configuration with $\Gamma=2$ (i.e. a standard $n=1$ polytrope). Our reference model is barotropic, $\Gamma_{1}=\Gamma$, and we compare it to two stratified models, with $\Gamma_{1}=2.05$ and $7/3$, respectively. The results from figure 3 suggest that these cases span the range of “reasonable” parameter values. Moreover, mode results and overlap integrals for the same values of $\Gamma_{1}$ are already available from ks , which means that we have an independent test of the numerics. However, we need to go one step further and add the ratio of the two eigenfunctions at the surface. The numerical results, listed in tables 1-3 (see also figure 4), demonstrate the relative importance of the g-modes for strongly stratified models, evident from the larger values of the overlap integrals $\tilde{Q}_{n}$ and the relative contributions ($k_{l}^{n}$) to the overall Love number.

The numerical results are interesting. In each case, the mode sum converges to the expected value for the Love number. For a barotropic model with $\Gamma=2$ the results should be $k_{l}\approx 0.259909$ willpo and the results from table 1 do, indeed, converge towards this number. This tells us that the sum over the star’s different oscillation modes provides an alternative representation for the Love number. It is important to note that this is true also for the stratified models, see table 4. The fundamental f-mode provides the dominant contribution in all cases (as expected given that this mode most closely resembles the tidal driving force), but in order to have a precise representation we need to account for both pressure and gravity modes. In fact, for the (somewhat extreme) $\Gamma_{1}=7/3$ case the first g-mode is more important than the p-modes.

Should we be surprised about these results? Probably not. As long as the modes of the star are complete, the mode-sum representation of the Love number must converge to the barotropic answer (after all, we can determine $k_{l}$ in the static limit). The contribution of individual modes changes with $\Gamma_{1}$, but the final answer is always the same. Nevertheless, the results are important. They provide the first actual demonstration that the dynamical response of the star, usually discussed in terms of mode resonances l94 ; ks ; wynn , has as its static limit the tidal deformability. In essence, we have an explicit link between the tidal deformability and asteroseismology. The results also point us in the direction we have to go if we want to quantify the impact of matter composition on the tidal response. We need to consider the dynamics encoded in the frequency dependence of (89).

IV.3 The dynamical tide

Having demonstrated that the sum over the star’s oscillation modes provides a precise description of the tidal response in the static limit, let us turn to the dynamical response associated with finite frequencies. The static result gives us confidence that this is a sensible exercise, obviously closely related to previous discussions of resonant mode excitation in binary systems l94 ; ks but at the same time adding insight into the role of the (no longer static) $m=0$ contribution. Taking (89) as our starting point and introducing dimensionless frequencies, as before, we have

[TABLE]

Next we use (97) and normalise the modes to arrive at

[TABLE]

This is the final result. It provides a closed expression for the frequency dependent tidal response (encoded in $k_{l}$) and, given the numerical results from tables 1-3 it is straightforward to obtain the desired dynamical behaviour. This, in turn, allows us to quantify the level at which each individual mode contributes to the overall result. Results for the three different values of $\Gamma_{1}$ we have considered are presented in figure 5. The different panels show the relative contributions to the tidal deformability (compared to that of the f-mode alone). The resonances associated with each mode, which occur when $\tilde{\omega}=\tilde{\omega}_{n}$, are easily distinguishable in each case and the resonance associated with the f-mode leads to a common feature (at $\tilde{\omega}\approx 1.2$) in all panels. The results tell us that modes other than the f-mode contribute to the overall result at the few percent level. Moreover, the results for $\Gamma_{1}=7/3$ bring out the fact that, in this case the leading g-mode dominates over the first p-mode for all frequencies. However, for the (likely) more realistic case with $\Gamma_{1}=2.05$ the g-mode contribution is almost an order of magnitude weaker than that of the first p-mode. These results are important as they provide the first (we believe) demonstration of the level at which frozen matter composition impacts on the tidal response across the range of frequencies relevant for a binary inspiral. Of course, the notion that resonant mode-excitation comes into play is not at all new. However, the discussion tends to focus on $m\neq 0$ modes (which would make a dynamical contribution for systems with a fixed orbital separation). The fact that the $m=0$ contribution has similar features (which should come into play for evolving orbits, when the function $f(\omega)$ is non-trivial, see the discussion in section 2) appears to not have been appreciated previously.

Finally, the overall mode contribution to the tidal response (including only the modes listed in tables 1-3) is illustrated in figure 6. This figure requires some additional explanation. We know from table 4 that the modes we include lead to a truncation error at the $10^{-3}$ level. This sets the level for the horizontal part of the curves towards the left of the figure. The inclusion of further modes should bring this level down (up to numerical precision). Instead, we should focus on the behaviour at frequencies beyond the first mode resonance. This illustrates the difference between the barotropic models and the stratified cases, again suggesting that matter composition may affect the result at the percent level.

V Implications

We have discussed the tidal response of a neutron star during the late stages of neutron star binary inspiral. In particular, we have focussed on the role of the matter composition. This issue has previously been ignored as studies have almost exclusively focussed of barotropic fluid models. However, it is natural to argue (given the timescale involved) that the matter composition should remain “frozen” during the late stages of binary inspiral, leading to a stratified perturbation problem (where the adiabatic index of the perturbation is different from that of the equilibrium background). This connects with previous work on tidal resonances, which has quantified the role of the g-modes (which obviously rely on stratification for their existence). Similarly, the proposed nonlinear pg-instability arras ; wein ; essick ; ligo_pg relies on the coupling between p-modes and g-modes. Of course, one may argue that the impact of stratification on the tidal response should be small enough that it can safely be ignored. Indeed, our numerical results indicate that the difference is at (or below) the level of a few percent. However, it is nevertheless important to quantify this contribution. We need to do this in order to understand systematic “errors” associated with the assumed physics, which ultimately determines the accuracy with which we can hope to extract stellar parameters like the radius from observations. Today’s gravitational-wave detectors are not at a level where a change of a percent in the tidal response makes much difference, but one might want to keep an eye on these issues for future reference. The problem is also important from the physics point of view. Having quantified the level at which the matter composition enters the discussion we can compare to (for example) the role of the elastic crust penner and superfluid components yu . Finally, and perhaps most importantly, our discussion suggests a “new” approach to dynamical tides, putting the mode excitation in focus, not only for resonances associated with the $m\neq 0$ modes but also for the $m=0$ contribution. This in interesting as it may allow us to develop new phenomenological models for the all-important gravitational-wave phasing (as an alternative to the effective Love number prescription from kleff ; stein , see also the recent effort in schmidt ). A first step in this direction is outlined in nap .