Role of crustal physics in the tidal deformation of a neutron star

TL;DR

This paper investigates how realistic crustal physics, including shear modulus, influences the tidal deformability of neutron stars during inspiral, finding small but potentially detectable corrections to gravitational wave signals.

Contribution

It introduces a fully general relativistic model incorporating realistic crust equations of state and shear modulus, quantifying their impact on neutron star tidal deformability.

Findings

Realistic crust EOSs cause about 1% correction in the second Love number.

Crustal shear effects on magnetic-type Love number are negligible.

Corrections are subdominant for current detectors but relevant for third-generation observatories.

Abstract

In the late inspiral phase, gravitational waves from binary neutron star mergers carry the imprint of the equation of state due to the tidally deformed structure of the components. If the stars contain solid crusts, then their shear modulus can affect the deformability of the star and, thereby, modify the emitted signal. Here, we investigate the effect of realistic equations of state (EOSs) of the crustal matter, with a realistic model for the shear modulus of the stellar crust in a fully general relativistic framework. This allows us to systematically study the deviations that are expected from fluid models. In particular, we use unified EOSs, both relativistic and non-relativistic, in our calculations. We find that realistic EOSs of crusts cause a small correction, of $\sim 1\%$, in the second Love number. This correction will likely be subdominant to the statistical error expected in LIGO-Virgo observations at their respective advanced design sensitivities, but rival that error in third generation detectors. For completeness, we also study the effect of crustal shear on the magnetic-type Love number and find it to be much smaller.

Tables

Table 1. Table 1: Properties of unified EOSs are tabulated here. Specifically, n 0 subscript 𝑛 0 n_{0} is the saturation density, n t subscript 𝑛 𝑡 n_{t} is the crust-core transition density, K 𝐾 K is the incompressibility, J 𝐽 J and L 𝐿 L are symmetry energy and its slope at saturation density, respectively.

EOS	$n_0$(fm-3)	$K$(MeV)	$J$(MeV)	$L$(MeV)	$n_t$(fm-3)	$M_{\max }/M_{\odot }$	${\mu }_{n=n_t}\left({\mathrm{gcm}}^{-1}{s}^{-2}\right)$
SLy4 Douchin:2001sv	0.159	230.0	32.0	46.0	0.0800	2.05	2.34 $\times {10}^{30}$
KDE0V1 Agrawal:2005ix	0.165	227.5	34.6	54.7	0.0480	1.97	$4.43\times {10}^{29}$
SkI4 Reinhard:1995zz	0.160	248.0	29.5	60.4	0.0359	2.18	$2.24\times {10}^{29}$
NL3 Lalazissis:1996rd ; Grill:2014aea	0.148	270.7	37.3	118.3	0.0548	2.77	$2.20\times {10}^{29}$
NL3$\omega \rho$ Horowitz:2002mb ; Grill:2014aea	0.148	272.0	31.7	55.3	0.0835	2.75	$3.42\times {10}^{30}$
DDME2 Lalazissis:2005de ; Grill:2014aea	0.152	250.9	32.3	51.2	0.0735	2.48	$2.21\times {10}^{30}$
NL3 matched	0.148	270.7	37.3	118.3	0.0740	2.77	$2.73\times {10}^{30}$
NL3$\omega \rho$ matched	0.148	272.0	31.7	55.3	0.0740	2.75	$2.73\times {10}^{30}$
DDME2 matched	0.152	250.9	32.3	51.2	0.0760	2.48	$2.89\times {10}^{30}$

Equations

$$d{s}^2=g_{\alpha \beta }d{x}^{\alpha }d{x}^{\beta }=e^{\nu (r)}d{t}^2-e^{\lambda (r)}d{r}^2-r^{2}d{\theta }^2-r^2{\sin }^2\theta d{\varphi }^2,$$

$$e^{\lambda (r)}={\left[1-\frac{2m(r)}{r}\right]}^{-1}.$$

$$T_{\alpha \beta }=(\rho +P)u_{\alpha }{u}_{\beta }-P{g}_{\alpha \beta },$$

$$\frac{dP(r)}{dr}=-\frac{\left[\rho (r)+P(r)\right]\left[m(r)+4\pi r^{3}P(r)\right]}{r\left[r-2m(r)\right]}$$

$$\frac{d\nu (r)}{dr}=-\frac{1}{\rho (r)+P(r)}\frac{dP(r)}{dr}$$

$$\frac{dm(r)}{dr}=4\pi r^2\rho (r).$$

$$h_{\alpha\beta}(r)=\left(\begin{array}[]{cccc}H_{0}(r)e^{\nu}&0&0&0\\ 0&H_{2}(r)e^{\lambda}&0&0\\ 0&0&r^{2}K(r)&0\\ 0&0&0&r^{2}\sin^{2}\theta K(r)\end{array}\right)P_{2}(\cos\theta).$$

$${\xi }^r=\frac{e^{-\lambda /2}}{r^2}W{P}_2(\cos \theta ),$$

$${\xi }^{\theta }=-\frac{V}{r^2}{\partial }_{\theta }{P}_2(\cos \theta ).$$

$$W^{{}^{\prime }}=\frac{r^2{e}^{\lambda /2}}{2}(-K+H_0)+\frac{3W}{r}+3V{e}^{\lambda /2}$$

$$V^{{}^{\prime }}=\frac{e^{\lambda /2}W}{r^2}+\frac{2V}{r}$$

$$K^{{}^{\prime }}=H_0{\nu }^{{}^{\prime }}+H_0^{{}^{\prime }}$$

$$\begin{split}&H_{0}^{{}^{\prime\prime}}+H_{0}^{{}^{\prime}}\bigg{[}\frac{2}{r}+e^{\lambda}\Big{(}\frac{2m(r)}{r^{2}}+4\pi r(P-\rho)\Big{)}\bigg{]}\\ &+H_{0}\left[-\frac{6e^{\lambda}}{r^{2}}+4\pi e^{\lambda}\left(5\rho+9P+\frac{\rho+P}{{c_{s}}^{2}}\right)-{\nu^{\prime}}^{2}\right]=0,\end{split}$$

$$\delta {\Pi }_r^r=A{Y}_{lm}=\frac{2\mu }{3}\left[K-H_2+2r{e}^{-\lambda /2}\left(\frac{W^{{}^{\prime }}}{r^3}-\frac{3W}{r^4}\right)-\frac{l(l+1)V}{r^2}\right]Y_{lm}$$

$$\delta {\Pi }_{\mathcal{A}}^r=B{Y}_{lm,\mathcal{A}}=\mu \left[\frac{e^{-\frac{\lambda }{2}}}{r^2}W-r^2{e}^{-\lambda }\left(\frac{V^{{}^{\prime }}}{r^2}-\frac{2V}{r^3}\right)\right]Y_{lm,\mathcal{A}}$$

$$\delta {\Pi }_{\mathcal{B}}^{\mathcal{A}}=2\mu V{Y}_{lm}{|}_{\mathcal{B}}^{\mathcal{A}}+\frac{\mu }{3}\left[H_2-K-2r{e}^{-\lambda /2}\left(\frac{W^{{}^{\prime }}}{r^3}-\frac{3W}{r^4}\right)-2\frac{l(l+1)V}{r^2}\right]{\delta }_{\mathcal{B}}^{\mathcal{A}}{Y}_{lm},$$

$$W^{{}^{\prime }}=\frac{r^2{e}^{\lambda /2}}{2}(\frac{3}{2\mu }A-K+H_0)+\frac{3W}{r}+(16\pi \mu r^2+3)V{e}^{\lambda /2}$$

$$V^{{}^{\prime }}=\frac{e^{\lambda /2}W}{r^2}+\frac{2V}{r}-\frac{B{e}^{\lambda }}{\mu }$$

$$\frac{\Delta n}{n}=\left[-\frac{e^{-\lambda /2}}{r^2}{W}^{{}^{\prime }}-\frac{6V}{r^2}+\frac{H_2}{2}+K\right],$$

$$\Delta \rho =(\rho +P)\frac{\Delta n}{n}$$

$$\Delta P=c_s^2\Delta \rho .$$

$$\Delta P=\delta P+{\xi }^r{P}^{{}^{\prime }}=\delta P-\frac{(\rho +P){\nu }^{{}^{\prime }}}{2r^2}{e}^{-\lambda /2}W.$$

$$\delta P=(P+\rho )c_s^2\left[-\frac{3}{4\mu }A+\frac{3}{2}K-\frac{9V}{r^2}+\frac{e^{-\lambda /2}}{r^3}\left(-3+\frac{r{\nu }^{{}^{\prime }}}{2c_s^2}\right)W\right].$$

$$16\pi r^2{e}^{\lambda }(\delta P-A)$$

$$-16\pi \mu V{r}^2({\nu }^{{}^{\prime }}{)}^2+16\pi e^{\lambda }rB(2+r{\nu }^{{}^{\prime }}).$$

$$K^{{}^{\prime }}=H_0{\nu }^{{}^{\prime }}+H_0^{{}^{\prime }}+\frac{16\pi \mu (r{\nu }^{{}^{\prime }}+2)V}{r}-16\pi B{e}^{\lambda }.$$

$$H_2=H_0+32\pi \mu V.$$

$$-\delta P=\frac{e^{-\lambda }}{16\pi r}({\nu }^{{}^{\prime }}+{\lambda }^{{}^{\prime }})H_0-\frac{4\mu V}{r^2}-\frac{B}{2r}(4+r{\lambda }^{{}^{\prime }}+r{\nu }^{{}^{\prime }})-B^{{}^{\prime }}+\frac{A}{2}.$$

$$B^{{}^{\prime }}=\frac{e^{-\lambda }}{16\pi r}({\nu }^{{}^{\prime }}+{\lambda }^{{}^{\prime }})H_0-\frac{B}{2r}(4+r{\lambda }^{{}^{\prime }}+r{\nu }^{{}^{\prime }})-\frac{4\mu V}{r^2}+\delta P+\frac{A}{2}.$$

$$\begin{gathered} -r^2{H}_0^{{}^{\prime \prime }}+\left(\frac{1}{2}r{\lambda }^{{}^{\prime }}-r{\nu }^{{}^{\prime }}-2\right)r{H}_0^{{}^{\prime }}+r^2{\nu }^{{}^{\prime }}{K}^{{}^{\prime }}-\frac{1}{2}{r}^2{\nu }^{{}^{\prime }}{H}_2^{{}^{\prime }}+6e^{\lambda }{H}_0 \\ +[2(e^{\lambda }-1)-r({\lambda }^{{}^{\prime }}+3{\nu }^{{}^{\prime }})]H_2=-8\pi r^2{e}^{\lambda }(3\delta P+\delta \rho ). \end{gathered}$$

$$\begin{gathered} -r^2{H}_0^{{}^{\prime \prime }}+\left(\frac{1}{2}r({\lambda }^{{}^{\prime }}-{\nu }^{{}^{\prime }})-2\right)r{H}_0^{{}^{\prime }}+[6e^{\lambda }+2(e^{\lambda }-1)-r({\lambda }^{{}^{\prime }}+3{\nu }^{{}^{\prime }})+r^2({\nu }^{{}^{\prime }}{)}^2]H_0 \\ =8\pi \{-r^2{e}^{\lambda }(3\delta P+\delta \rho )+8\mu \left[1-e^{\lambda }+r\left({\nu }^{{}^{\prime }}+\frac{1}{2}{\lambda }^{{}^{\prime }}\right)-\frac{1}{4}(r{\nu }^{{}^{\prime }}{)}^2\right]V \\ +2r^2{\nu }^{{}^{\prime }}(\mu V{)}^{{}^{\prime }}+2r^2{\nu }^{{}^{\prime }}B{e}^{\lambda }\}. \end{gathered}$$

Full text

††preprint: LIGO-P1900071

I Introduction

The detection of gravitational waves (GWs) from the binary neutron star (BNS) merger event GW170817 has ushered in a new probe for constraining the equation of state (EOS) of neutron stars Abbott2017 ; Abbott2019 ; AbbottGRB . Post-Newtonian theory predicts that the inspiral signal from BNSs will carry an imprint of the EOS of neutron stars  Flanagan2008 ; Hinderer2008 ; Binnington2009 ; Damour2009 . This result has been vetted by Numerical Relativity simulations as well. One can use this imprint to study the properties of dense matter far from the nuclear saturation density with an event having significantly high signal-to-noise ratio in the future generation of detectors Hinderer2010 . Neutron stars are also believed to have a solid crust as their outermost layers Chamel2008 . The effect of the crust is a crucial ingredient for probing nuclear physics through GWs. The pioneering work of Carter and collaborators that introduced the theory of elastic solids in general relativity (GR) Carter1972 ; Carter1973 has paved the way for studying the effect of NS crust in a consistent relativistic framework. Many of the studies concerning the NS perturbation incorporating crust elasticity have used the Cowling approximation Yoshida2002 ; Flores2017 . By contrast, there exist only a few studies that have accounted for full GR effects in the analysis Schumaker1983 ; Finn1990 .

In one of the first attempts in this direction, Penner et al. tried to extract tidal information by employing an elastic crust Penner2011 and followed it up in another work to study the crustal failure during BNS inspiral  Penner2012 . However, they used mostly modest details of dense matter and rudimentary crust models. Therefore, the tidal behavior of different crust models inspired from various realistic nuclear interactions has not been very clear from their results. At this point, it should be mentioned that there have been other studies that investigate the tidal deformability of solid quark star using similar GR perturbative framework Lau2017 ; Lau2018 . The point of interest in them is mostly the phenomenology of a solid core of a star forming due to a deconfinement phase-transition at the center. In our case, we solely focus on the effect of a solid crust encapsulating a fluid core.

The elastic properties of NS crust strongly depend on the composition of matter across a range of sub-nuclear densities. In the outer crust, the nuclei are arranged in the form of a body-centered-cubic (BCC) lattice that is embedded in a noninteracting and degenerate electron gas Baym:1971pw ; Nandi:2010fp . As the density increases with depth, the neutron-drip point ($\sim 3\times 10^{-4}$ fm*-3*) is reached, which signals the beginning of the inner crust. In this region, the neutron-rich nuclei are arranged as a lattice immersed in inter-penetrating gas of free neutrons and electrons Baym:1971ax ; Negele1973 ; Haensel2001 ; Nandi2011 . The shear modulus is also higher in the inner crust. Thus, the inner crust contributes the most towards the tidal response due to shear. This region extends till the crust-core transition density ($\sim 8\times 10^{-4}$ fm*-3*). Complex structures (e.g., rod, slab, bubble, etc. – collectively known as nuclear pasta) are expected to occur in the inner crust as the matter gradually changes from crystalline to homogeneous phase, with increasing density 1983PhRvL..50.2066R ; 1984PThPh..71..320H ; 2016PhRvC..94b5806N ; 2018ApJ…852..135N . Beneath the inner crust, the outer core starts, with uniform nuclear matter. As the density grows even higher, one reaches the inner core that can have superfluid neutrons and, perhaps, even exotic matter like strange baryons, meson condensates, quark matter, etc. Glendenning2001 ; Page2006 . The characteristics of the matter in the core is highly speculative. Additionally, the EOS of the crust is qualitatively different and represents different physical conditions than that in the core. Therefore, one has to rely on proper matching of both EOSs at the crust-core interface. This is very crucial as it has been found that a proper thermodynamically consistent matching is required to avoid large uncertainties on the macroscopic properties of the star Fortin2016 . Even then there might be some ambiguity due to different choices of crust-core transition density. The best way is to use unified EOS models where the EOSs of the crust and the core are calculated within the same underlying theory. Hence, we employ unified EOSs in our work.

The main aim of this paper is to provide a comprehensive picture of the interplay between the perturbative response of the elastic crust of a NS and the nuclear physics of the constituents of the crustal matter using several unified EOS models. We have reworked the formalism of Penner et al. Penner2011 using the analysis of the perturbed quantities from Finn Finn1990 . We find that realistic EOSs of crusts, with a non-zero shear modulus, cause a small correction, of $\sim 1\%$, in the second Love number. This correction will likely be subdominant to the statistical error expected in LIGO-Virgo observations at their respective advanced design sensitivities, but rival that error in third generation detectors.

The paper is organized as follows. In Sec. II, we discuss the formalism for tidal deformation and derivation of the Love numbers in elastic relativistic stars. Thereafter, we present an overview of the EOS used in our calculation in Sec. III. We discuss our results in Sec. IV and summarize in Sec. V. Throughout our analysis we have assumed $c=G=1$.

II Analytical set up

In this section, we present the analytical formulation of our work. Our focus is to calculate the tidal deformation of neutron stars with a solid crust. A solid crust supports shear stress and as a result two different types of pulsation mode arise. The odd parity modes are called torsional modes which creates twist in the star. These are first discussed in 1983 by Schumaker and Thorne Schumaker1983 . Then in 1990 Finn Finn1990 presented a new set of even parity type modes for a solid star. Since our focus is to calculate the tidal deformation of a neutron star we only need static perturbation equations which are basically zero frequency modes of pulsation problem. The set of static polar perturbation equations for a solid star are first given by Penner et al Penner2011 . However, we found some inconsistencies in their equations, in particular, they don’t match with the zero frequency limit of pulsation equations given by Finn Finn1990 . Recently, Lau et al. Lau2018 has also pointed out same sort of inconsistencies. Therefore, we re-derive all the equations in this work. Throughout the paper, we follow the notation of Thorne and Campolattaro Thorne1967 as was also adopted by Finn Finn1990 , so that we can easily verify our equations with that of Finn’s at zero frequency limit.

Our aim is to quantify the effect of elasticity of crust in the tidal deformabilty of neutron stars. For this we first calculate the background of the star by solving the standard Tolman-Oppenheimer-Volkoff (TOV) equation. In the next step we consider a static linear perturbation of this background model which takes into account elastic crust. To study the linear perturbation we expand each components of fluid displacement vector and perturbed metric in terms of spherical harmonics. Each spherical harmonics is characterized by $l$, $m$ and parity, which can be either even $(-1)^{l}$ or odd $(-1)^{l+1}$. Under small amplitude motions these two parity decouple from each other and, hence, can be treated separately. Here we consider both cases individually and compute the deviation in Love number.

II.1 Background problem

The equilibrium of a static, spherically symmetric relativistic star is given by,

[TABLE]

where $\nu$ and $\lambda$ are two metric functions, and $\lambda$ can be expressed in terms of mass $m(r)$ inside a radius of r,

[TABLE]

Here we are interested in neutron stars which have a fluid core and a solid crust. We assume that in equilibrium configuration the contribution of shear stress due to the presence of solid crust vanishes Finn1990 ; Penner2011 . In reality, this assumption is not necessarily correct. But since we are interested in small amplitude perturbation we can think that the background shear is almost negligible and its contribution is important only in the perturbed configuration. Therefore, the contribution of shear stress only enters through the perturbed stress energy tensor. The advantage of this assumption is that it makes our background problem very simple as we can now use the perfect fluid stress-energy tensor to model the background star:

[TABLE]

where $u_{\alpha}$, $\rho(r)$ and $P(r)$ denote fluid 4-velocity, energy density and pressure, respectively, inside the star. Solving Einstein equation for this equilibrium configuration we arrive at TOV equation,

[TABLE]

[TABLE]

[TABLE]

For a cold neutron star it is reasonable to assume that this fluid does not exchange heat with the surroundings. Therefore, one can take the EOS to be a zero-temperature barotrope: $P=P(\rho)$. Given the EOS of neutron stars, Eqs. (4) and (6) can be solved to obtain their mass-radius relationship.

II.2 Even parity perturbation

II.2.1 Fluid perturbation equations

First we consider $l=2$, static, even-parity perturbations in the Regge-Wheeler gauge. Since we are only interested in quadrupole deformation, $l=2$ case is considered from the beginning. Also we further simplify equation of motion by choosing spherical harmonics with $m=0$. Under these assumptions perturbed metric becomes,

[TABLE]

The contravariant component of fluid displacement field takes the form

[TABLE]

[TABLE]

For the case of perfect fluid all the off-diagonal components of perturbed stress-energy tensor vanish. The non-vanishing components of perturbed stress-energy tensor are

$\delta T_{0}^{0}=\delta\rho$ ,

$\delta T_{i}^{i}=-\delta P$ .

For a barotrope we can assume the following form of perturbed pressure:

$\delta P=\frac{dP}{d\rho}\delta\rho=c_{s}^{2}\delta\rho$  ,

where $c_{s}^{2}$ is the speed of sound inside the star.

The set of equations which describes the fluid problem is given by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where the prime (${}^{{}^{\prime}}$) denotes derivative respect to $r$. Basically, only a single differential equation, namely of $H_{0}$ is sufficient to determine the tidal Love number of a fluid star. Rest of the three coupled differential equations are needed just to join the fluid core of the star with the solid crust.

II.2.2 Elastic perturbation equations

We assume our background star to be relaxed and unstrained. An elastic crust does not affect the equilibrium model. The contribution of elasticity comes only through the perturbed stress-energy tensor. Therefore, the perturbed stress energy tensor in solid medium becomes,

$\delta T_{\alpha\beta}=\delta T_{\alpha\beta}^{fluid}+\delta\Pi_{\alpha\beta}$  ,

where $\delta\Pi_{\alpha\beta}$ is anisotropic stress energy tensor. Detailed derivation of calculating this anisotropic stress energy tensor is given by Finn and Penner et al. Finn1990 ; Penner2011 . The non-vanishing components of $\delta\Pi_{\alpha\beta}$ are,

[TABLE]

[TABLE]

[TABLE]

where $\mathcal{A}$ and $\mathcal{B}$ run over the coordinates $\theta$ and $\phi$. The vertical ($|$) denotes covariant derivatives on two sphere. Also we have defined two terms $A$ and $B$ which are related to radial and tangential component of the anisotropic shear stress tensor. Now equation (14) and (15) can be used to integrate $W$ and $V$ (only $l=2$, $m=0$ case has been considered here):

[TABLE]

[TABLE]

In Regge-Wheeler gauge the perturbed number density takes the form Thorne1967 :

[TABLE]

and the corresponding Lagrangian changes in density and pressure are:

[TABLE]

[TABLE]

The Lagrangian change in pressure is related to the Eulerian change as:

[TABLE]

Combining (20), (21) and (22) we get an expression for the perturbed Euler pressure:

[TABLE]

The $[rr]$ component of perturbed Einstein tensor gives another expression for $\delta P$:

[TABLE]

By solving these two algebraic equations, (23) and (24), we calculate $\delta P$ and $A$. $[r\theta]$ component leads to equation of motion for $K$:

[TABLE]

Subtraction of $[\phi\phi]$ component from $[\theta\theta]$ component leads to:

[TABLE]

We write the sum of $[\theta\theta]$ and $[\phi\phi]$ component in terms of $A$ and $B$

[TABLE]

We can use this equation to integrate $B$:

[TABLE]

If we take the trace of perturbed Einstein equation we arrive at a second order equation for $H_{0}$:

[TABLE]

After plugging (25) and (26) in the above equation we get

[TABLE]

The five differential equations (17), (18), (25), (28) and (30) given above together with two algebraic equations (23) and (24) form a complete set of equations which describe the evolution of perturbed quantities in the elastic medium of the star.

II.2.3 Boundary condition at center and stellar surface

At the center of the star, all the perturbed quantities must be regular. For our study we take core of the star to be fluid, for which the boundary conditions were analyzed by Thorne and Campolattaro Thorne1967 . Here we just summarize their result. All the perturbed quantities are expanded in Taylor series about $r=0$ as:

[TABLE]

Using these expansions in equation (19) we get $W^{(0)}=-lV^{(0)}$. Next, by combining equations (24) and (28) for $\mu=0$ we obtain :

[TABLE]

It is straightforward to show that expansion of this equation about $r=0$ leads to $K^{(0)}=H^{(0)}_{0}$. Therefore, out of four constants appearing in equations (31)-(34), only two are independent. These two are fixed by the demand that the Lagrangian perturbation of pressure vanishes at the surface of the star.

II.2.4 Interface condition

We have derived the perturbation equations in the solid crust region. Now we need to find proper interface conditions to join them with the fluid perturbation equations in the core. The interface conditions are obtained from the equations of motion of fluid variables and Einstein Field equation (please see Finn1990 for the derivation). The continuity of intrinsic curvature demands that $H_{0}$, $K$, $W$ must be continuous at the interface:

[TABLE]

Again, continuity of extrinsic curvature imposes two additional boundary conditions:

[TABLE]

[TABLE]

Since, $W$ is continuous across the interface equation (39) reduces to :

[TABLE]

By noting that $A=0$ in the fluid core we obtain the value of radial stress at the interface as :

[TABLE]

where $\delta P_{f}$ and $\delta P_{i}$ are the Eulerian perturbations of pressure at the base of fluid core and at the interface, respectively. Using equations (23) and (28) we get the expression of $\delta P_{f}$ and $\delta P_{i}$ :

[TABLE]

[TABLE]

II.2.5 Calculation of electric tidal Love number

Our focus here is to calculate the electric love number of neutron stars consisting of a fluid core and an elastic crust. We first integrate the fluid perturbation equations starting from the center of the star to the core-crust junction, using the specified boundary conditions at the center. Next, we integrate the elastic perturbation equations from this junction to the surface. The starting point of the later integration is obtained by imposing the interface conditions at the core-crust junction. Now, in order to calculate the tidal Love number we have to match this internal solution with the external solution at the surface of the star. We suggest the reader to see Refs. Hinderer2008 ; Binnington2009 ; Damour2009 , where extensive details about the calculation of tidal Love number can be found. The value of tidal Love number can be computed in terms of $y$ and compactness parameter $C=\frac{M}{R}$ as :

[TABLE]

where $y$ depends on the value of $H_{0}$ and its derivative at the surface:

[TABLE]

II.3 Odd parity perturbation:

II.3.1 Fluid perturbation equations:

Magnetic tidal Love number were computed together by Binnington and Poisson Binnington2009 (BP) and Damour and Nagar Damour2009 (DN) back in 2009. In their calculation, BP assumed that tidal field varies slowly over the time, therefore it never throws the body out of hydrostatic equilibrium. Based on this assumption they derived all the perturbation equations using a static-fluid ansatz and from there they calculated magnetic tidal Love number. On the other hand, instead of re-deriving the perturbation equations, DN took the Cunningham, Price and Moncrief master function Cunningham1978 governing odd parity perturbation of Schwarszchild space-time and used a stationary perfect-fluid ansatz for stress-energy tensor. In 2015 Landry and Poisson Landry2015 (LP) revisited BP’s calculation by taking an irrotational state of the fluid which permits internal motions of fluid inside the body. They find that magnetic tidal Love number for this irrotational state are different from the magnetic Love number associated with static fluid. LP also found that there results agree with DN’s result since irrotational condition is automatically imposed by the stationary master function chosen by DN. This state of affair is recently re-examined by Pani et al. Pani2018 . In our work, we allow internal motion of fluid since it is the more realistic configuration to describe the fluid than the hydrostatic equilibrium scenario.

We consider here magnetic type perturbation for $l=2,m=0$ in the Regge-Wheeler gauge by a time dependent tidal field. However we assume tidal field varies very slowly over the time, therefore, we neglect all the time derivative appears in our field equations. But this slowly varying tidal field does have impact on the internal motion of the fluid which establishes irrotational state of it. Under the above mentioned assumption the perturbed metric becomes:

[TABLE]

The contravariant fluid displacement vector has the following form:

[TABLE]

where $U(r,t)$ is the fluid displacement function for odd parity perturbation. Perturbed four-velocities corresponds to this first order in displacement, are

[TABLE]

Since density and pressure are scalar they do not change under odd parity perturbation. However fluid four velocity will be shifted to $u^{\mu}$ to $u^{\mu}+\delta u^{\mu}$. In first order perturbation we note the following relation, $\delta u_{\mu}=g_{\mu\nu}\delta u^{\nu}+h_{\mu\nu}u^{\nu}$ to compute the components of $\delta u_{\mu}$,

[TABLE]

where, $v_{\mu}=g_{\mu\nu}\delta u^{\nu}$. Now, irrotational state of fluid implies $\delta u_{r}=0=\delta u_{\mathcal{A}}$ Landry2015 , where $\mathcal{A}$ runs over the coordinate $\theta$ and $\phi$. Consequently, the form of perturbed stress-energy tensor can be written as,

[TABLE]

Therefore, for the irrotational case, the $[t\phi]$ component of Einstein equation give us,

[TABLE]

For the static case, $v_{\phi}=0$ gives $\delta u_{\phi}=h_{t\phi}u^{t}$. In that case, the $[t\phi]$ component of Einstein equation give us,

[TABLE]

Notice that assumption of irrotational state of fluid changes the sign of $\frac{\lambda^{{}^{\prime}}+\nu^{{}^{\prime}}}{r}$ in the term which is proprtinal to $h_{0}$.

Since, we will be working with irrotational state of fluid, We need to solve (52) with the regular boundary condition at the center:

[TABLE]

where $h_{0}^{(0)}$ is an arbitrary constant.

II.3.2 Elastic perturbation equations:

In this case the non-vanishing components of the anisotropic stress energy tensor are:

[TABLE]

After including these anisotropic stress energy-tensor in the perturbed Einstein equations the $[t\phi]$ gives us:

[TABLE]

We solve equation (52) from the center to core-crust junction and equation (58) from there to the stellar surface where at the core-crust junction $h_{0}$ is continuous.

II.3.3 Magnetic Love number calculation:

Asymptotic behavior of equation (52) at large distances is given by,

[TABLE]

By matching its asymptotic solution with $g_{t\phi}$ component of metric in asymptotically mass-centered Cartesian coordinate (ACMC)Thorne1980 ; Cardoso2017 , we obtain the expression of magnetic Love number

[TABLE]

where $y=\frac{rh_{0}^{{}^{\prime}}}{h_{0}}$ is evaluated at the surface of the star and $C$ is the compactness.

III Equation of state

It has been shown Fortin2016 that for an unambiguous calculation of NS properties (especially radius and crust thickness) it is necessary to adopt unified equations of state (EOS), where the EOS of crust and core are obtained within the same many-body theory. As the crust thickness plays the key role in the present study we employ only unified EOSs here. We consider six unified EOSs: SLy4 Douchin:2001sv , KDE0V1 Agrawal:2005ix , SkI4 Reinhard:1995zz , NL3 Lalazissis:1996rd ; Grill:2014aea , NL3$\omega\rho$ Horowitz:2002mb ; Grill:2014aea and DDME2 Lalazissis:2005de ; Grill:2014aea . The first three are based on non-relativistic Skyrme interactions and are obtained from the CompOSE database compose ; Gulminelli:2015csa . The other three are derived from the relativistic mean-field (RMF) model. The RMF EOSs are not fully unified as the outer crust is not calculated within the same theory but is taken from Ref. Baym:1971pw . Since the most part of the outer crust is determined from the experimentally measured nuclear masses, the choice of it does not significantly affect the observables. The important properties of these EOSs are shown in Table 1. All of them are consistent with the observed maximum mass ($2.01\pm 0.04M_{\odot}$) of neutron stars Antoniadis:2013pzd . The shear moduli for all the EOSs are plotted in Fig. 1. They are calculated using the following expression 1991ApJ…375..679S ; 2016PhRvC..94b5801N :

[TABLE]

where $a=[3/(4\pi n_{i})]^{1/3}$, $n_{i}$ is the density of ions and $Z$ is the atomic number of the nucleus present.

IV Results

In this section we present our numerical findings for a set of neutron star EOSs that were discussed in Sec. III. First, for each of these EOSs we generated a set of equilibrium stellar configurations within the mass range of $1M_{\odot}$ to $2M_{\odot}$ by solving the TOV equations. In Fig. 2, crustal thickness is plotted w.r.t mass for each considered EOSs. Then for each star we integrate all the perturbed variables and calculate the corresponding Love number. For our numerical calculations, we use dimensionless variables in all the necessary differential equations, which are presented in the Appendix Appendix: Dimensionless form. All the differential equations are solved using fourth order Runge-Kutta method.

IV.1 Even parity perturbations:

IV.1.1 Non relativistic EOS

In Fig. 3, we plot the change in $k_{2}$ (in $\%$) due to the inclusion of the crust as a function of mass for three non-relativistic EOSs. The fractional change in $k_{2}$ is defined as $\Delta k_{2}/k_{2}^{\rm fluid}$, where $\Delta k_{2}=k_{2}^{\rm fluid}-k_{2}^{\rm crust}$; $k_{2}^{\rm fluid}$ and $k_{2}^{\rm crust}$ are, respectively, the tidal Love numbers of a purely fluid star and a star with an elastic crust. Since, the elastic crust would resist deformation, it is expected that $k_{2}^{\rm crust}<k_{2}^{\rm core}$, resulting in $\Delta k_{2}>0$. This is indeed the case as can be seen from Fig. 3. It is also observed that as the thickness of the crust increases, the change in Love number increases (for comparison, see Fig. 2). The change in $k_{2}$ is about $~{}0.1-0.4\%$ for KDE0V1 EOS, $~{}0.3-0.9\%$ for SLy4 EOS and $~{}0.1-0.2\%$ for SkI4 EOS. For a given mass, the increasing order of crustal thickness among these three non-relativistic unified EOSs is: SLy$4>$ KDE0V1 $>$ SkI4. A similar trend is seen for the change in $k_{2}$ in Fig. 3, which points to the fact that stars with bigger crusts would have lesser deformation, as expected.

IV.1.2 Relativistic EOS

Change in $k_{2}$ as a function of mass is plotted in Fig. 4 for three unified RMF EOSs. To investigate the importance of unified EOS, we also include EOSs obtained by matching a crust EOS with the core EOS in a thermodynamically consistent way Fortin2016 . All three RMF EOSs of core are matched to the BPS+BBP Baym:1971pw ; Baym:1971ax EOS of the crust at $n_{t}=n_{0}/2$, where $n_{0}$ is the saturation density of the core EOS (see table 1 for values). From Fig. 2 we see that the crust is bigger for the matched EOS than that of unified EOS for NL3. Whereas, for DDME2 and NL3$\omega\rho$ the scenario is opposite. The reason is for unified NL3, $n_{t}=0.0548$ fm*-3* and for the matched case $n_{t}=n_{0}/2=0.074$ fm*-3*. As the transition from crust to core happens later in the matched EOS, the crust is bigger in size for matched NL3. In contrast, the transition happens earlier for matched EOS in case of DDME2 and NL3$\omega\rho$, which leads to a reduction in the crustal thickness.

For non-relativistic EOSs we have already seen that larger crustal thickness gives larger deviation from the fluid value of $k_{2}$. A similar behavior seen to be in play here. However, the magnitude of the shear modulus also plays a role in the change of $k_{2}$. Models with higher magnitude of shear modulus will be less deformed, which will correspond to a larger change in $k_{2}$. A careful inspection of Figs. 1 and 4 supports this finding. Depending on the magnitude of crustal thickness and shear modulus, one of the effects dominates over the other. For example, unified DDME2 has larger crustal thickness than matched DDME2, but the shear modulus of the latter is higher than that of the first. However, we find that the change in $k_{2}$ is higher for matched DDME2. Since the difference in crustal thickness between these two EOSs is very small, the change in $k_{2}$ is mainly caused by the magnitude of the shear modulus in this case. For reference, we have included the values of the shear modulus at the bottom of the crust in the last column of Table 1. Overall, the change in $k_{2}$ due to the presence of a solid crust is between $\sim 0.4-1.3\%$ for all the RMF EOSs considered here. It is also noted that $\Delta k_{2}$ for matched EOS can considerably differ from that of unified EOS. Among the three RMF EOSs studied here we found that the difference is highest for the NL3 EOS and can be as large as $\sim 90\%$. This emphasizes the necessity of the use of unified EOSs in such calculations.

The above results show that even with realistic EOSs and realistic crustal models, shear modulus of solid crust has a small effect on the electric tidal Love number. The reason why the effect of shear modulus on $k_{2}$ is negligible can be understood from Fig. 5, where the profile of $y$ is plotted as a function of radius in the presence and absence of shear using unified DDME2 EOS. The dotted (red) and solid (black) curves represent the cases of with and without shear, respectively. It is seen that the values of $y$ mainly differ from the fluid case in the inner crust region, while in the outer crust region difference from the fluid case is negligible. This happens because the magnitude of shear modulus is much higher in the inner crust and as a result the inner crust has greater response towards the tidal field than the outer crust. However, as the value of $y$ at the surface only enters into the calculation of $k_{2}$ (see eq. 46), we do not observe any significant change in it.

IV.2 Numerical results for odd parity perturbations:

IV.2.1 Non-Relativistic EOS

In Fig. 6, we plot the percentage change in magnetic Love number w.r.t mass for three unified non-relativistic EOSs. The change in magnetic Love number is denoted as $\Delta j_{2}/j_{2}^{\rm fluid}$, where $\Delta j_{2}=\mid j_{2}^{\rm fluid}-j_{2}^{\rm crust}\mid$; $j_{2}^{\rm fluid}$ and $j_{2}^{\rm crust}$ are respectively magnetic Love number of a perfectly fluid star and a star with elastic crust. The value of magnetic Love number is itself negative and that’s why the absolute values of difference have been taken. Similar to the electric Love number, the magnetic love number is also found to be affected by both the values of crustal thickness and the shear modulus. The change in $j_{2}$ varies between $\sim 0.00005-0.0005\%$ for considered non-relativistic EOSs. This suggests deviation in magnetic Love number due to solid crust is negligible.

IV.2.2 Relativistic EOS

For RMF EOSs we observe similar type of changes in $j_{2}$. From Fig. 1 we see that the value of shear modulus is similar for all the matched EOSs but the crustal thickness is different (see, Fig 2) with values in the order : NL3$>$NL$3\omega\rho>$DDME2. We see similar order of change in $j_{2}$ for matched EOSs. On the other hand, the change in $j_{2}$ is dominated by the magnitude of shear modulus for the unified EOSs as both the shear modulus (see, Fig. 1) and the change in $j_{2}$ have same order: NL$3\omega\rho>$DDME2$>$NL3.  However, similar to non-relativistic EOSs, here also the changes in magnetic Love number are in between $\sim 0.0001-0.0009\%$ and hence, are practically negligible.

IV.3 Summary of the results:

Here we briefly summarize the key findings of our numerical studies. We have analyzed the effect of elastic crust on both electric and magnetic tidal Love numbers for a set of realistic equations of state and realistic models of shear modulus. In our study, we used 3 non-relativistic and 3 relativistic EOSs. Non-relativistic EOSs are based on Skyrme interactions and relativistic EOSs are constructed using relativistic mean field theory. We also used matched RMF EOS where core RMF EOSs are matched to BPS EOS of the crust at half of saturation density of the core EOS.

• •

Effect on electric Love number: Percentage change in $k_{2}$ is higher for relativistic EOSs than the non-relativistic ones. The reason is the RMF EOSs have larger crustal thickness. We observed that the EOSs with larger crustal thickness have a larger change in $k_{2}$ w.r.t the fluid case. Also we observe larger shear modulus corresponds to larger change in $k_{2}$. So we conclude that the crustal thickness and the magnitude of the shear modulus both have an effect on the electric tidal Love number to varying degrees depending on the EOS.

• •

Effect on magnetic Love number: Similar to electric Love number both the crustal thickness and magnitude of shear modulus have effect on the magnetic Love number. Note, however, that the magnetic Love number is much smaller than the electric one for all of these EOSs and therefore it is highly unlikely that their imprints will be observed in BNS waveforms, let alone these corrections, in the era of Advanced LIGO and Virgo.

• •

Comparison between Penner et al. and our analysis:

Finally, we compare the results of our analysis with that of Penner et al. Penner2011 . In their analysis they used a polytropic model for EOS and a simple profile of shear modulus that varies linearly with pressure. We implement the same EOS and shear modulus profile in our analysis and compare the obtained results with that of Pennal et al. in Fig. 8. We see that the changes in $k_{2}$ for our calculation is up to two orders of magnitude larger than their analysis. However, the change in $k_{2}$ is still small and the effect of solid crust in $k_{2}$ is unlikely to be observed by LIGO and Virgo detectors in the near future. This may, however, change in the detectors of the subsequent generation.

V Conclusion

In this paper, we have investigated the effect of elastic crust on the tidal deformation of neutron stars. We presented a complete set of static perturbation equations both for fluid core and solid crust of a neutron star. We verified here that our static perturbation equations are consistent with the zero frequency limit of Finn Finn1990 but at variance with some of the expressions in Penner et al. Penner2011 . Recent independent calculation by Lau et al. Lau2018 have also pointed out that the set of static perturbation equations given by Penner et al. Penner2011 are inconsistent with zero frequency limit of pulsation equations of Finn Finn1990 . Our paper should be seen as an extension of the work done by Penner et al. Penner2011 who used a simple model of neutron stars based on a polytropic EOS and a simple linear profile for crustal shear modulus. In this paper, we have investigated the effect of realistic EOSs and realistic model of shear modulus on the tidal deformability of neutron stars. We find that realistic EOSs and shear modulus can cause a change of $\sim 1\%$ in electric tidal Love number, much larger than found by Penner et al. While this change may not be of much consequence for LIGO-Virgo observations in the near future, they may be important for subsequent generations of detectors.

Acknowledgements

We would like to thank Sayak Datta and Kabir Chakravarti for helpful discussions. We also thank Wolfgang Kastaun for carefully reading the manuscript and making useful comments. This work was supported in part by the Navajbai Ratan Tata Trust. Bhaskar Biswas acknowledges the University Grant Commission (UGC) India, for the financial support as a senior research fellow.

Appendix: Dimensionless form

V.1 Background equations

We introduce following dimensionless variables:

[TABLE]

where $\rho_{c}$ and $p_{c}$ are central density and pressure, respectively. So, the TOV equations take the form

[TABLE]

with

[TABLE]

V.2 Perturbation equations for even parity

Here we introduce additional dimensionless variables as:

[TABLE]

Now we get following equations (${}^{\prime}=d/dx$):

[TABLE]

[TABLE]

Boundary conditions at the center (for $l=2$) become :

[TABLE]

Where $a$ and $c$ are constants. Interface conditions are :

[TABLE]

V.3 Perturbation equations for odd parity

Dimensionless form of odd parity perturbed equation is:

[TABLE]