Tidal Deformations of Hybrid Stars with Sharp Phase Transitions and Elastic Crusts

The recent detection of gravitational waves (GWs) from binary compact systems (Abbott et al. 2016a, 2016b, 2017b, 2017c, 2018, 2019a, 2020) has inaugurated the era of GW astronomy, and we now have the real possibility of addressing many issues for compact astrophysical objects. After the first detection of GWs from neutron stars (NSs; Abbott et al. 2017c), independent constraints from X-ray observations—which rely on pulse profile modeling (see, e.g., Özel et al. 2016; Miller et al. 2019; Riley et al. 2019; de Lima et al. 2020 and references therein)—have already been obtained, for instance on their equation of state (EOS) and radii (Annala et al. 2018; De et al. 2018; Most et al. 2018; Nandi & Char 2018; Paschalidis et al. 2018). Further constraints are possible due to inferences of tidal deformations (Hinderer 2008; Damour & Nagar 2009; Hinderer et al. 2010), which already have known upper limits (Abbott et al. 2017c, 2018, 2019b). Future GW detections hold the promise of stronger constraints on the NS parameters (Abbott et al. 2017a; Chirenti et al. 2017; Baiotti 2019). We focus here on whether GW observations may be able to provide evidence for elastic and hybrid phases in NS interiors.

Hybrid NSs are compact systems presenting both quark and hadronic phases. They might be possible in nature due to the existence of phase transitions in quantum chromodynamics (QCD; see Alford et al. 2008; Paschalidis et al. 2018; Bauswein et al. 2019; Alford et al. 2019). However, the order of this phase transition is unknown in the parameter range relevant for NSs (low temperatures and large chemical potentials; Alford et al. 2008). In addition, since the surface tension for dense matter is poorly known, the possibility of a mixed phase remains an open issue (Bombaci et al. 2016; Lugones & Grunfeld 2017). Astrophysical observations cannot currently exclude hybrid NSs (Bejger et al. 2017). However, this might be possible with future, more precise GW measurements from the late inspiral and also the postmerger phase (Most et al. 2018; Nandi & Char 2018; Bauswein et al. 2019). Another promising way of probing hybrid stars is through pulse profile modeling, timing, and polarimetry measurements, i.e., with the NICER (Özel et al. 2016) and eXTP (Zhang et al. 2016) missions. The former has already provided the first measurements of an NS radius with an uncertainty of around 10% at the 68% confidence level (see Bilous et al. 2019; Bogdanov et al. 2019a, 2019b; Guillot et al. 2019; Miller et al. 2019; Raaijmakers et al. 2019; Riley et al. 2019). The eXTP mission promises even more constraining measurements (∼5% uncertainties or less; Zhang et al. 2016). NICER's precision is already in an interesting range to probe the interiors of NS (see, e.g., Sieniawska et al. 2018), especially hybrid stars that allow for additional stable branches ("third family") with smaller radii than their hadronic counterparts (Alford & Sedrakian 2017; Alvarez-Castillo & Blaschke 2017; Christian & Schaffner-Bielich 2019).

It is already known that tidal deformations of hadronic stars with solid/elastic crusts and liquid cores hardly differ at all from perfect-fluid stars in ordinary cases (Penner et al. 2011; Gittins et al. 2020), which assume a continuous energy density at the interface between a (sizable) core and a solid crust. However, the situation may be different for hybrid stars. As is well known, the presence of density jumps in hybrid stars with sharp phase transitions strongly influences the tidal deformation in the case of perfect fluids (see, e.g., Damour & Nagar 2009; Postnikov et al. 2010; Han & Steiner 2019). To the best of our knowledge, studies have not been carried out that probe the consequences of combined density jumps and shear stresses in the hadronic phase of a hybrid NS. Therefore, we investigate how this issue impacts GW observations. This question is complementary to the existing studies of hybrid stars in the context of GW170817 (Paschalidis et al. 2018; Alford et al. 2019; Bauswein et al. 2019; Han & Steiner 2019; Li et al. 2020; Montaña et al. 2019; Sieniawska et al. 2019; Essick et al. 2020), which assume perfect-fluid systems. Natural motivations for our study are that a discontinuous quark-hadronic density profile might non-negligibly change the tidal deformability of solid/elastic hybrid stars, and a nonzero shear modulus right from the bottom of the hadronic phase might be seen as a simplistic model for an elastic mixed phase.

It is also known that solid aspects of a quark phase could significantly change the tidal deformation of NSs (Lau et al. 2017, 2019). QCD allows for such a possibility through the so-called LOFF phase, where shear stresses could be up to a thousand times larger than crustal ones (Mannarelli et al. 2007). It has been reported that for purely crystalline quark stars, tidal deformations could decrease up to 60% with respect to their perfect-fluid counterparts (Lau et al. 2017). On the other hand, if the quark phase is located inside a hadronic phase, tidal deformations are much subtler: (i) if the (elastic) quark phase extends beyond around 70% of the whole star, it still follows that the tidal deformations could change considerably; (ii) in other cases, changes with respect to a perfect fluid are expected to be small (Lau et al. 2019).

For the above case in particular, the lower limits on the tidal deformations coming from GW170817 (Radice et al. 2018; Coughlin et al. 2019; Kiuchi et al. 2019; Radice & Dai 2019) are important for probing hybrid star models and the solidity of the quark phase. In addition, maximum deformations of stars with elastic quark phases might lead to observable effects in current GW detectors (Haskell et al. 2007) and may allow for the estimation of breaking strains of solid stars (Ushomirsky et al. 2000). Although the issue of the maximum deformation of solid stars is very important for determining their maximum ellipticities in general, in this work, we only focus on the cases of small deformations of elastic hadronic regions of stars (solid/elastic crusts). The reason for doing so is that the elastic properties of such matter are currently better understood.

The article is arranged as follows: in Section 2, we review the main ingredients for calculating the tidal deformations of perturbed stars with static tidal deformations. Section 3 is devoted to presenting simple quark and hadronic models for a hybrid star, which can be used to estimate tidal deformations when part of the hadronic phase is elastic. In Section 4, we obtain the appropriate boundary conditions for hybrid stars presenting sharp phase transitions with no mixed phases. Our results regarding the tidal deformations of hybrid stars with solid crusts are presented in Section 5. We discuss the main issues raised by our work and provide the relevant conclusions in Section 6.

We work with geometric units. Our metric signature convention is +2. Unless otherwise stated, ΔX is defined as the Lagrangian perturbation (Shapiro & Teukolsky 1983; Andersson & Comer 2007; Penner et al. 2011; Andersson et al. 2019) of a physical quantity X, while δX is the corresponding Eulerian perturbation.

The formalism we follow has been laid down by Penner et al. (2011). Here, we briefly review it, focusing on the equations we will use (for further details, see Gittins et al. 2020). For the background hybrid stars, fluids are assumed to be perfect and, hence, described by the Tolman–Oppenheimer–Volkoff (TOV) system of equations, namely,

$$\begin{eqnarray}&&\displaystyle \frac{{dp}}{{dr}}=-\displaystyle \frac{\epsilon m}{{r}^{2}}\left(1+\displaystyle \frac{p}{\epsilon }\right)\left(1+\displaystyle \frac{4\pi {{pr}}^{3}}{m}\right){\left(1-\displaystyle \frac{2m}{r}\right)}^{-1},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{dm}}{{dr}}=4\pi {r}^{2}\epsilon \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d\nu }{{dr}}=-\displaystyle \frac{2}{\epsilon +p}\displaystyle \frac{{dp}}{{dr}},\,\,\,\,\,{e}^{\lambda (r)}:=\,{\left[1-\displaystyle \frac{2m(r)}{r}\right]}^{-1},\end{eqnarray} \tag{ 3 }$$

where p is the pressure, is the energy density, and m is the gravitational mass at radial distance r. The λ and ν functions are related to the background metric, assumed to be given by

$$\begin{eqnarray}&&{{ds}}^{2}=-{e}^{\nu (r)}{{dt}}^{2}+{e}^{\lambda (r)}{{dr}}^{2}+{r}^{2}(d{\theta }^{2}+{\sin }^{2}\theta d{\phi }^{2}).\end{eqnarray} \tag{ 4 }$$

The tidal deformations of stars are taken to be related to non-radial perturbations with zero frequency. We work within the Regge-Wheeler gauge (Regge & Wheeler 1957), in which case (even parity) metric perturbations (${h}_{{ab}}=\delta {g}_{{ab}}$) are diagonal and defined as

$$\begin{eqnarray}\begin{array}{rcl}{h}_{{ab}} & = & {\rm{diag}}[{H}_{0}(r){e}^{\nu (r)},{H}_{2}(r){e}^{\lambda (r)},{r}^{2}k(r),\\ & & {r}^{2}{\sin }^{2}\theta k(r)]{Y}_{l}^{m}(\theta ,\phi ),\end{array}\end{eqnarray} \tag{ 5 }$$

where $({H}_{0},{H}_{2},k)$ are functions to be determined by the Einstein field equations and Y^m_l are the spherical harmonics. Without any loss of generality, and due to the spherically symmetric background, we work with m = 0 and constrain our analysis to quadrupole deformations (l = 2).

The fluid perturbations are assumed to be given by (Penner et al. 2011)

$$\begin{eqnarray}&&{\xi }^{r}=\displaystyle \frac{W(r)}{r}{P}_{2};\,\,{\xi }^{\theta }=\displaystyle \frac{V(r)}{{r}^{2}}\displaystyle \frac{{{dP}}_{2}}{d\theta };\,\,{\xi }^{\phi }=0,\end{eqnarray} \tag{ 6 }$$

where ${P}_{2}(\theta )\propto {Y}_{2}^{0}$ is the Legendre polynomial of second order and W(r) and V(r) are functions to be determined by the field equations and conservation laws. Covariant components of the above perturbations are obtained readily through ${\xi }_{a}={g}_{{ab}}{\xi }^{b}$. Note that there is no time dependence in the terms of Equation (6), meaning that only u_t is different from zero. From the usual condition ${u}_{a}{u}^{a}=-1$ and Equation (5), we then get

$$\begin{eqnarray}&&{u}_{t}={u}_{t}^{0}+\delta {u}_{t}={e}^{\tfrac{\nu }{2}}\left(1-\displaystyle \frac{1}{2}{H}_{0}{P}_{2}\right),\end{eqnarray} \tag{ 7 }$$

where δ denotes the Eulerian perturbation.

We also take the baryon number conservation into account. It geometrically implies that the Lagrangian change of the baryon number density is given by (Penner et al. 2011; Andersson et al. 2019)

$$\begin{eqnarray}&&{\rm{\Delta }}n=\delta n+{\xi }^{a}{n}_{;a}=-\displaystyle \frac{n}{2}{{ \mathcal P }}^{{ab}}{\rm{\Delta }}{g}_{{ab}},\end{eqnarray} \tag{ 8 }$$

where the semicolon $(;)$ represents the covariant derivative (Landau & Lifshitz 1975), and ${{ \mathcal P }}_{{ab}}:={g}_{{ab}}+{u}_{a}{u}_{b}$ is the projector onto the orthogonal directions of the four-velocity u^a. In addition,

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}{g}_{{cd}} & = & {h}_{{cd}}+{\xi }_{c;d}+{\xi }_{d;c}\\ & = & {h}_{{cd}}+{\partial }_{c}{\xi }_{d}+{\partial }_{d}{\xi }_{c}-2{{\rm{\Gamma }}}_{{cd}}^{a}{\xi }_{a},\end{array}\end{eqnarray} \tag{ 9 }$$

where ${{\rm{\Gamma }}}_{{cd}}^{a}$ represents the usual Christoffel symbols (connection coefficients; Landau & Lifshitz 1975). The Eulerian change of the energy-momentum tensor of the fluid is assumed to be given by (Penner et al. 2011; Andersson et al. 2019)

$$\begin{eqnarray}&&\delta {T}_{a}^{b}={\left(\delta {T}_{a}^{b}\right)}_{{\rm{perf}}}+\delta {{\rm{\Pi }}}_{a}^{b},\end{eqnarray} \tag{ 10 }$$

where the first term is the perturbation of the energy-momentum tensor of the perfect-fluid $[{(\delta {T}_{a}^{b})}_{{\rm{perf}}}={\rm{diag}}[-\delta \epsilon (r),\delta p(r),\delta p(r),\delta p(r)]{P}_{2}(\theta )]$ and the second term represents the contribution from shear stresses (Penner et al. 2011; Andersson et al. 2019) such that

$$\begin{eqnarray}&&\delta {{\rm{\Pi }}}_{a}^{b}=-\tilde{\mu }\left({{ \mathcal P }}_{a}^{c}{{ \mathcal P }}^{{db}}-\displaystyle \frac{1}{3}{{ \mathcal P }}_{a}^{b}{{ \mathcal P }}^{{cd}}\right){\rm{\Delta }}{g}_{{cd}},\end{eqnarray} \tag{ 11 }$$

where $\tilde{\mu }$ is the shear modulus. We note that this geometric formalism assumes perturbations with respect to an unstrained background configuration. For further details, see Andersson et al. (2019) and references therein.

As stressed by Penner et al. (2011), one can check that for time-independent perturbations $\delta {{\rm{\Pi }}}_{t}^{b}=0$. From Equations (11), (9), and (8), it also follows that

$$\begin{eqnarray}\begin{array}{rcl}\delta {{\rm{\Pi }}}_{r}^{r} & := & \delta {\tilde{{\rm{\Pi }}}}_{r}^{r}{P}_{2}=\displaystyle \frac{2\tilde{\mu }{P}_{2}}{3{r}^{2}}[{r}^{2}(k-{H}_{2})\\ & & +\,(4-r\lambda ^{\prime} )W-2{rW}^{\prime} -6V],\end{array}\end{eqnarray} \tag{ 12 }$$

where primes indicate radial derivatives. Also,

$$\begin{eqnarray}\begin{array}{rcl}\delta {{\rm{\Pi }}}_{\theta }^{\theta } & = & \displaystyle \frac{\tilde{\mu }}{3{r}^{2}}[\{-{r}^{2}(k-{H}_{2})-(4-r\lambda ^{\prime} )W+2{rW}^{\prime} \}{P}_{2}\\ & & +\,6V(2{P}_{2}-1)],\end{array}\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&\delta {{\rm{\Pi }}}_{\phi }^{\phi }=\displaystyle \frac{\tilde{\mu }}{3{r}^{2}}[\{-{r}^{2}(k-{H}_{2})-(4-r{\lambda }^{{\rm{{\prime} }}})W+2{rW}^{{\rm{{\prime} }}}\}\,{P}_{2}+6V],\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\delta {{\rm{\Pi }}}_{r}^{\theta }:=\delta {\tilde{{\rm{\Pi }}}}_{r}^{\theta }\displaystyle \frac{{{dP}}_{2}}{d\theta }=-\displaystyle \frac{\tilde{\mu }}{{r}^{3}}\displaystyle \frac{{{dP}}_{2}}{d\theta }[{e}^{\lambda }W-2V+{rV}^{\prime} ].\end{eqnarray} \tag{ 15 }$$

As expected, from Equations (12)–(14), the trace of $\delta {{\rm{\Pi }}}_{b}^{a}$ vanishes. Following Penner et al. (2011), for future reference, we define

$$\begin{eqnarray}&&{T}_{\theta }(r)\,:=\,-2\tilde{\mu }[{e}^{\lambda }W-2V+{rV}^{\prime} ]\end{eqnarray} \tag{ 16 }$$

and

$$\begin{eqnarray}&&{T}_{r}(r)\,:=\,\displaystyle \frac{2}{3}\tilde{\mu }[{r}^{2}(k-{H}_{2})+(4-r\lambda ^{\prime} )W-2{rW}^{\prime} -6V],\end{eqnarray} \tag{ 17 }$$

related to the polar and radial components of the traction in the star (${T}_{{ab}}{n}^{b}$, where n^b is the normal four-vector to a given hypersurface), respectively, which helps us to apply the relevant boundary conditions.

Now we write down the Einstein equations relevant to our analysis. When one subtracts the $[\theta \theta ]$ component of the Einstein equations from the $[\phi \phi ]$ component, with the help of Equations (13) and (14), it follows that

$$\begin{eqnarray}&&{H}_{2}(r)={H}_{0}(r)+32\pi \tilde{\mu }V(r),\end{eqnarray} \tag{ 18 }$$

confirming that solid stars spoil the equality between H₂ and H₀ (Penner et al. 2011; Krüger et al. 2015), which holds for perfect-fluid stars (Hinderer 2008; Hinderer et al. 2010). The $[r\theta ]$ component of the Einstein equations tells us that

$$\begin{eqnarray}&&k^{\prime} (r)=\displaystyle \frac{8\pi }{r}(4\tilde{\mu }V-{T}_{\theta })+{H}_{0}^{{\prime} }+{H}_{0}\nu ^{\prime} +16\pi \tilde{\mu }V\nu ^{\prime} .\end{eqnarray} \tag{ 19 }$$

The sum of $[\theta \theta ]$ and $[\phi \phi ]$ components, with the use of Equations (13), (14), (18), and (19) and the background equations, leads to

$$\begin{eqnarray}\begin{array}{rcl}\delta p & = & \displaystyle \frac{{T}_{r}}{2{r}^{2}}+\displaystyle \frac{{e}^{-\lambda }}{4{r}^{2}}\{{T}_{\theta }[r(\lambda ^{\prime} -\nu ^{\prime} )-2]-16{e}^{\lambda }\tilde{\mu }V-2{{rT}}_{\theta }^{\prime} \}\\ & & +\,\displaystyle \frac{{e}^{-\lambda }}{16\pi r}{H}_{0}(\lambda ^{\prime} +\nu ^{\prime} ),\end{array}\end{eqnarray} \tag{ 20 }$$

which determines the radial part of the Eulerian change of the pressure. The $[{rr}]$ component of the Einstein equations can be simplified with Equations (18) and (19), and implies another relation for k(r),

$$\begin{eqnarray}\begin{array}{rcl}4{e}^{\lambda }k & = & {H}_{0}\left[6{e}^{\lambda }-2-r(\lambda ^{\prime} +\nu ^{\prime} )+{\left(r\nu ^{\prime} \right)}^{2}\right]++{r}^{2}{H}_{0}^{{\prime} }\nu ^{\prime} \,+\\ & & +\,8\pi {e}^{\lambda }\left[8\tilde{\mu }V-3{T}_{r}+2{e}^{-\lambda }\tilde{\mu }{r}^{2}V{\left(\nu ^{\prime} \right)}^{2}\right]\,-\\ & & -\,4\pi {T}_{\theta }[2+r(\lambda ^{\prime} +\nu ^{\prime} )]+8\pi {{rT}}_{\theta }^{\prime} .\end{array}\end{eqnarray} \tag{ 21 }$$

Another relevant equation is the trace of the Einstein equations. By subtracting the $[{tt}]$ component from $[{rr}]+[\theta \theta ]+[\phi \phi ]$, one gets

$$\begin{eqnarray}&&\begin{array}{l}16\pi {r}^{2}{e}^{\lambda }(\delta \epsilon +3\delta p)=-2{H}_{0}\left[2-8{e}^{\lambda }+r(\lambda ^{\prime} +3\nu ^{\prime} )\right.\\ \quad \left.-\,{r}^{2}{\left(\nu ^{\prime} \right)}^{2}\right]-{H}_{0}^{{\prime} }r[4-r(\lambda ^{\prime} -\nu ^{\prime} )]-2{r}^{2}{H}_{0}^{\prime\prime} -16\pi {{rT}}_{\theta }\nu ^{\prime} +\\ \quad +32\pi \tilde{\mu }V[{r}^{2}{\left(\nu ^{\prime} \right)}^{2}-2r(\lambda ^{\prime} +2\nu ^{\prime} )+4{e}^{\lambda }-4]\\ \quad -\,32\pi {r}^{2}\nu ^{\prime} (\tilde{\mu }V)^{\prime} ,\end{array}\end{eqnarray} \tag{ 22 }$$

where we recall that $\delta \epsilon =\delta \epsilon (r)$ is the radial part of the Eulerian change of the energy density of the star. For completeness, the $[{tt}]$ component of the Einstein equations implies that

$$\begin{eqnarray}\begin{array}{rcl}8\pi {r}^{2}\delta \epsilon & = & -{e}^{-\lambda }{r}^{2}k^{\prime\prime} -{e}^{-\lambda }\left(3-\displaystyle \frac{1}{2}r\lambda ^{\prime} \right){rk}^{\prime} +2k\\ & & +\,{e}^{-\lambda }{{rH}}_{2}^{\prime} +\left[3+{e}^{-\lambda }(1-r\lambda ^{\prime} )\right]{H}_{2}.\end{array}\end{eqnarray} \tag{ 23 }$$

Consistency with the case where perturbations have zero frequency is only obtained when the NS EOS is barotropic (Penner et al. 2011), that is, when p = p(), which results in

$$\begin{eqnarray}&&\delta \epsilon =\left(\displaystyle \frac{\partial \epsilon }{\partial p}\right)\delta p\equiv \displaystyle \frac{1}{{c}_{s}^{2}}\delta p,\end{eqnarray} \tag{ 24 }$$

where c_s is the adiabatic speed of sound.

From thermodynamical considerations, one can easily show that (Penner et al. 2011)

$$\begin{eqnarray}&&{\rm{\Delta }}\epsilon \equiv \left(\delta \epsilon +\displaystyle \frac{W}{{{rc}}_{s}^{2}}p^{\prime} \right){P}_{2}=(p+\epsilon )\displaystyle \frac{{\rm{\Delta }}n}{n}.\end{eqnarray} \tag{ 25 }$$

From Equations (5), (6), (8), and (9), it follows that

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}n & = & -\displaystyle \frac{n}{2{r}^{2}}\left[{r}^{2}({H}_{2}+2k)+W(r\lambda ^{\prime} +2)\right.\\ & & +\,\left.2W^{\prime} r-12V\right]{P}_{2}.\end{array}\end{eqnarray} \tag{ 26 }$$

When the above equation is combined with Equations (24) and (25), one gets an additional equation for $\delta p$, namely,

$$\begin{eqnarray}\begin{array}{rcl}\delta p & = & -{c}_{s}^{2}\displaystyle \frac{p+\epsilon }{2{r}^{2}}\left[{r}^{2}({H}_{2}+2k)+W\left(r\lambda ^{\prime} +2-\displaystyle \frac{r\nu ^{\prime} }{{c}_{s}^{2}}\right)\right.\\ & & +\,\left.2W^{\prime} r-12V\right].\end{array}\end{eqnarray} \tag{ 27 }$$

Other equations, e.g., associated with the components of ${({T}_{b}^{a})}_{;a}=0$, may also be deduced. They are the, usually not obvious, consequence of the field equations; hence, it is enlightening to write them down. From ${({T}_{r}^{a})}_{;a}=0$, it follows that

$$\begin{eqnarray}&&\begin{array}{l}\delta {{\rm{\Pi }}}_{r}^{r}\left(\displaystyle \frac{6}{r}+\nu ^{\prime} \right)-(\epsilon +p){H}_{0}^{{\prime} }{P}_{2}+2(\delta {{pP}}_{2}+\delta {{\rm{\Pi }}}_{r}^{r})^{\prime} +\\ \ \ \times \,(\delta p+\delta \epsilon ){P}_{2}\nu ^{\prime} -12\delta {\tilde{{\rm{\Pi }}}}_{r}^{\theta }{P}_{2}+2p^{\prime} +(p+\epsilon )\nu ^{\prime} =0.\end{array}\end{eqnarray} \tag{ 28 }$$

Consider now ${({T}_{\theta }^{a})}_{;a}=0$, which implies that

$$\begin{eqnarray}&&\begin{array}{l}2\tilde{\mu }[{r}^{2}({H}_{2}-k)+18V+(r\lambda ^{\prime} -4)W+2{rW}^{\prime} ]\\ \ \ +3{r}^{2}\left\{\Space{0ex}{3.35ex}{0ex}2\delta p-{H}_{0}(p+\epsilon )+{e}^{-\lambda }{r}^{2}\left[\Space{0ex}{3.15ex}{0ex}2(\delta {\tilde{{\rm{\Pi }}}}_{r}^{\theta })^{\prime} \right.\right.\\ \quad \left.\left.+\,\delta {\tilde{{\rm{\Pi }}}}_{r}^{\theta }\left(\displaystyle \frac{8}{r}-\lambda ^{\prime} +\nu ^{\prime} \right)\right]\right\}=0.\end{array}\end{eqnarray} \tag{ 29 }$$

This equation reduces to $\delta p={H}_{0}(\epsilon +p)/2$ when $\tilde{\mu }=0$ (perfect-fluid case). This is in full agreement with Equation (20) because $8\pi {{re}}^{\lambda }(p+\epsilon )=(\lambda ^{\prime} +\nu ^{\prime} )$. These equations will be useful when we deduce the appropriate boundary conditions to the perturbations (Section 4).

In this first approach to hybrid stars, we assume a simplified model that contains some aspects of more realistic systems. For the quark phase, we mostly assume the MIT bag-like EOS (Alford et al. 2005), such that the grand thermodynamic potential is given by

$$\begin{eqnarray}&&{\rm{\Omega }}=-p=-\displaystyle \frac{3}{4{\pi }^{2}}{a}_{4}{\mu }^{4}+\displaystyle \frac{3}{4{\pi }^{2}}{a}_{2}{\mu }^{2}+B,\end{eqnarray} \tag{ 30 }$$

where $\mu \equiv ({\mu }_{u}+{\mu }_{s}+{\mu }_{d})/3$ is the averaged quark chemical potential, a₄ is a constant accounting for the strong interactions between quarks, a₂ is another constant that takes into account quark finite masses, quark pairing (color superconductivity), etc., and B is a third constant encompassing the nontrivial vacuum of QCD (Alford et al. 2005). In the above expression, we have ignored the grand thermodynamic potential of electrons, ${{\rm{\Omega }}}_{e}$, because it is usually much smaller than that of the quarks (Pereira et al. 2018). From Equation (30), one can easily obtain other relevant thermodynamic quantities, such as the baryon number density ${n}_{b}\doteq -1/3(\partial {\rm{\Omega }}/\partial \mu )$ and $\epsilon =3\mu {n}_{b}-p$. From the above relations, one could isolate μ and get the following EOS (Pereira et al. 2018):

$$\begin{eqnarray}&&p(\epsilon )=\tfrac{1}{3}(\epsilon -4B)-\displaystyle \frac{{a}_{2}^{2}}{12{\pi }^{2}{a}_{4}}\left[1+\sqrt{1+\displaystyle \frac{16{\pi }^{2}{a}_{4}}{{a}_{2}^{2}}(\epsilon -B)}\right],\end{eqnarray} \tag{ 31 }$$

which generalizes the usual bag model. For the hadronic phase, we consider either an effective polytropic EOS of the form $p=K{\epsilon }^{2}$ (n = 1 of $p=K{\epsilon }^{1+\tfrac{1}{n}}$, K = 100 km² (Penner et al. 2011)) or work with hybrid star models involving the SLy4 EOS (Douchin & Haensel 2001; Sieniawska et al. 2019) and a relativistic mean field theory using the NL3 parameterization (Pereira et al. 2018), connected to the BPS EOS (Baym et al. 1971) for lower densities. For the model containing the SLy4 EOS, following Sieniawska et al. (2019), we assume that it is connected to a relativistic polytropic EOS (Tooper 1965) at a baryon density n₀ (freely chosen) with pressure $p={\kappa }_{{ef}}{n}_{b}^{\gamma }$ and mass-energy density $\epsilon ({n}_{b})=p({n}_{b})/(\gamma -1)+{n}_{b}{{mc}}^{2}$, where ${\kappa }_{{ef}}$, γ, and m_b are free parameters (m_b denoting the mass of the baryon in that phase, and n_b denoting the baryon density), and extends up to n₁ (transition baryon density of choice), where there is a sharp phase transition with a given energy density jump (free parameter). For baryon densities larger than n₁, we assume a quark phase described by a simple approximation to the bag model (Zdunik 2000), in the form of $p=\alpha (\epsilon -{\epsilon }_{* })$, where α is a constant representing the square of the speed of sound, and _* is a free parameter denoting the density of the quark phase at zero pressure. Continuity of the pressure and specification of the energy density at the hadronic-quark interface uniquely fixes _*. For the particular models we will use, we set γ = 4.5, n₁ = 0.335 fm⁻³ (above twice the nuclear saturation density, ${\epsilon }_{{\rm{sat}}}=2.7\times {10}^{14}$ g cm⁻³) and α = 1 to also investigate properties of stiff quark matter (see Figure 9 of Sieniawska et al. 2019 and the text therein for more details).

Table 1 summarizes the main aspects of the hybrid star models we use in this work. We assume the Maxwell construction (constancy of the pressure and the chemical potential at a sharp phase-splitting surface; Bauswein et al. 2019) to determine the aspects of background hybrid stars. For the HS1 and HS3 models, we take the density jump at the quark-hadronic phase interface as a free parameter. In the case of the HS1 model in particular, given an η: = _q/_h−1 (possible values for η are in the range 0–2.0), where _q is the density at the top of the quark phase and _h is the density at the bottom of the hadronic phase, from p_q = p_h at the quark-hadronic phase-transition surface, one can easily find either _q or _h, which then allows for the determination of the phase-transition pressure and the chemical potential. For the HS2 model, η is uniquely found by the Maxwell construction.

Table 1.  Main Aspects of the Hybrid Star Models Used in Our Work

Hybrid	Quark Model	η	Hadronic
Model	$\left(\displaystyle \frac{{B}^{\tfrac{1}{4}}}{{\rm{MeV}}},{a}_{4},\tfrac{{a}_{2}^{\tfrac{1}{2}}}{{\rm{MeV}}}\right)$	$\left(\tfrac{{\epsilon }_{q}}{{\epsilon }_{h}}-1\right)$	Model
HS1	(137, 0.40, 100)	free	Polyt. (n = 1)
HS2	(140, 0.55, 100)	≈0.45	BPS+NL3
HS3	Bag with ${c}_{s}^{2}=1$	free	SLy4+Polyt. (γ = 4.5)

Download table as:  ASCIITypeset image

We take the quark phase to be a perfect fluid, even though it is also theoretically possible that the quark phase is elastic (Mannarelli et al. 2007). Unless otherwise stated, we assume that the hadronic phase is elastic only for densities below ≈ 2 × 10¹⁴ g cm⁻³ (≈2/3 of the nuclear saturation density; Ducoin et al. 2011). For the shear modulus, we take the simple linear model (Chamel & Haensel 2008; Penner et al. 2011)

$$\begin{eqnarray}&&\tilde{\mu }(r)=\kappa p(r)+{\tilde{\mu }}_{0},\end{eqnarray} \tag{ 32 }$$

where κ and ${\tilde{\mu }}_{0}$ are given constants.

Naturally, our models (for the EOSs and the shear modulus) are simplistic but they are expected to give us the main aspects that hybrid stars with elastic hadronic regions (solid crusts) should present, which could be used as input for more accurate descriptions. We assume that the solid crust ends at the density 10⁷ g cm⁻³, where the envelope (or the ocean) of an NS is supposed to start (Haensel et al. 2007; Penner et al. 2011).³

In order to integrate the Einstein equations and obtain the tidal deformations, one needs to impose appropriate boundary conditions at the center, the liquid-elastic interfaces, and the surface of the star. We describe them here.

At the center of a hybrid star, we impose that all functions are regular. In other words, we expand a given function F as (Krüger et al. 2015)

$$\begin{eqnarray}&&F(r)={F}_{0}+\displaystyle \frac{1}{2}{F}_{2}{r}^{2}+{ \mathcal O }({r}^{4}),\end{eqnarray} \tag{ 33 }$$

where the coefficients F₀ and F₂ are to be determined from the field equations (the symmetry of the field equations render the first-order terms null; Krüger et al. 2015). One also needs to expand the background quantities , p, λ, and ν in the same fashion. One can always express second-order corrections to the background (F₂ in Equation (33)) in terms of the solutions to the TOV equations; for the explicit expressions, see Krüger et al. (2015). Naturally, when working numerically, one could directly find second-order corrections to the background in terms of their second derivatives. When zeroth and second-order coefficients to the perturbations are found, one has the initial conditions to start the core integration.

For a liquid-elastic interface, boundary conditions could be easily found from aspects of the extrinsic curvature. We restrict ourselves to the case of no surface degrees of freedom (continuity of the extrinsic curvature, see Poisson 2004). We assume that a liquid-elastic interface is described by ${\rm{\Psi }}=r-{R}_{{le}}-{\xi }^{r}=0$, where R_le is the background liquid-elastic phase-transition radius.⁴ Standard calculations for the extrinsic curvature (K^a_b, $\{a,b\}=t,\theta ,\varphi$; Finn 1990; Gittins et al. 2020) lead to the following nontrivial components at $r={R}_{{le}}$:

$$\begin{eqnarray}&&\delta {K}_{t}^{t}=-4{e}^{-\tfrac{1}{2}\lambda ({R}_{{le}})}{P}_{2}{\left({H}_{0}^{{\prime} }+\displaystyle \frac{1}{2}{H}_{2}\nu ^{\prime} -\displaystyle \frac{W}{r}\nu ^{\prime\prime} +\displaystyle \frac{1}{2}\displaystyle \frac{W}{r}\nu ^{\prime} \lambda ^{\prime} \right)}_{r={R}_{{le}}}\end{eqnarray} \tag{ 34 }$$

and

$$\begin{eqnarray}&&\delta {K}_{\theta }^{\theta }=\delta {K}_{\phi }^{\phi }=-\displaystyle \frac{4}{{R}_{{le}}^{2}}{e}^{-\tfrac{1}{2}\lambda ({R}_{{le}})}{P}_{2}{\left({{rH}}_{2}-{r}^{2}k^{\prime} +W\lambda ^{\prime} \right)}_{r={R}_{{le}}}.\end{eqnarray} \tag{ 35 }$$

The continuity of the extrinsic curvature implies that the terms in parentheses in the above equations have null jumps. (Note that the background components of the extrinsic curvature (Gittins et al. 2020) are automatically continuous when the Maxwell construction is assumed.) From Equation (34), one readily finds that

$$\begin{eqnarray}&&[{H}_{0}^{{\prime} }{]}_{-}^{+}=8\pi {e}^{\lambda ({R}_{{le}})}\displaystyle \frac{W}{{R}_{{le}}}[\epsilon {]}_{-}^{+}-16\pi {\mu }^{+}{V}^{+}\nu ^{\prime} ,\end{eqnarray} \tag{ 36 }$$

where we have made use of $[\nu ^{\prime\prime} -\nu ^{\prime} \lambda ^{\prime} /2{]}_{-}^{+}=8\pi {e}^{\lambda }[\epsilon {]}_{-}^{+}$ (see, e.g., Landau & Lifshitz (1975) for the precise expression for $\nu ^{\prime\prime} -\nu ^{\prime} \lambda ^{\prime} /2$) and Equation (18). In addition, "+" ("−") relates to the bottom of an elastic phase (top of the liquid phase), which implies that ${\tilde{\mu }}^{-}=0$. The continuity of the induced metric at $r={R}_{{le}}$ also leads to $[W{]}_{-}^{+}=0$ (see Gittins et al. 2020 for further details). From Equation (6), it follows that $[{\xi }^{r}{]}_{-}^{+}=0$. It turns out that the null discontinuity of ${K}_{\theta }^{\theta }$ does not add a new jump condition. This could be readily seen from Equations (19) and (36) and the fact that $[\lambda ^{\prime} {]}_{-}^{+}=8\pi {{re}}^{\lambda }[\epsilon {]}_{-}^{+}$.

The interface condition of other functions could be more easily found when the perturbation equations are promoted to distributions.⁵ For alternative and equivalent results based on the continuity of the induced metric and extrinsic curvature, see Gittins et al. (2020). Given that we have second derivatives of H₀ (see Equation (22)), the promotion of the Einstein equations to distributions would only make sense if the jump of H₀ at R_le is null, i.e., $[{H}_{0}{]}_{-}^{+}=0$. Equation (28) then shows us that $[\delta p\,+p^{\prime} W/r+\delta {\tilde{{\rm{\Pi }}}}_{r}^{r}{]}_{-}^{+}=0$, which, from Equations (17) and (12), implies that $[{T}_{r}+{r}^{2}\delta p+{rWp}^{\prime} {]}_{-}^{+}=0$. From Equation (29), we have that $[\delta {\tilde{{\rm{\Pi }}}}_{\theta }^{r}{]}_{-}^{+}=0$, which, due to Equation (16), leads to $[{T}_{\theta }{]}_{-}^{+}=0$. From the above jump conditions, one has $[k{]}_{-}^{+}=0$, as can be easily seen from Equation (19) when promoted to a distribution.

At the surface of the star (r = R) and outside of it, the reasoning of Hinderer (2008), Damour & Nagar (2009), and Hinderer et al. (2010) ensue. In other words, if ${H}_{0}^{{\prime} }$ is continuous at the star's surface, one can obtain the Love number k₂ in terms of $y\equiv {{RH}}_{0}^{{\prime} }(R)/{H}_{0}(R)$ as (Hinderer 2008; Damour & Nagar 2009; Hinderer et al. 2010; Penner et al. 2011)

$$\begin{eqnarray}&&{k}_{2}=\displaystyle \frac{8{C}^{5}{\left(1-2C\right)}^{2}[2+2C(y-1)-y]}{5\{2C[6-3y+3C(5y-8)]+4{C}^{3}[13-11y+C(3y-2)+2{C}^{2}(1+y)]+3{\left(1-2C\right)}^{2}[2-y+2C(y-1)]\mathrm{ln}(1-2C)\}},\end{eqnarray} \tag{ 37 }$$

where $C\equiv M/R$ is the compactness of the background hybrid star. If ${H}_{0}^{{\prime} }$ is not continuous at the stellar surface, for instance due to a density discontinuity or a nonzero shear modulus, as in Equation (36), then k₂ also changes due to the discontinuity of y (Damour & Nagar 2009; Hinderer et al. 2010). From k₂, one can construct the physically relevant tidal deformation ${\rm{\Lambda }}\doteq 2/3{(M/R)}^{-5}{k}_{2}$. We will focus mainly on this quantity for the ease of comparison with the GW constraints.

We turn now to a simplification of boundary conditions at a liquid-elastic interface. We recall that in this case, "+" ("−" ) relates to the elastic (liquid) phase. It is easy to show that the condition $[{T}_{r}+{r}^{2}\delta p+{rWp}^{\prime} {]}_{-}^{+}=0$ implies

$$\begin{eqnarray}&&\begin{array}{c}\begin{array}{l}{V}^{+}=\displaystyle \frac{1}{16\pi {\tilde{\mu }}^{+}{R}_{le}^{2}{\nu }^{{\rm{{\prime} }}2}}[4{e}^{\lambda }\{{k}_{le}+4\pi {R}_{le}{W}_{le}({p}_{-}^{{\rm{{\prime} }}}-{p}_{+}^{{\rm{{\prime} }}})\}\\ \,\,-\,{R}_{le}^{2}{H}_{0}^{{\rm{{\prime} }}+}{\nu }^{{\rm{{\prime} }}}+{H}_{0}(2+2{e}^{\lambda }\{4\pi {R}_{le}^{2}({p}^{-}+{\epsilon }^{-})-3\}-{R}_{le}^{2}{\nu }^{{\rm{{\prime} }}2})],\end{array}\end{array}\end{eqnarray} \tag{ 38 }$$

with ${k}_{{le}}:={k}^{-}$ (Equation (21) evaluated at the top of the perfect-fluid phase where $\tilde{\mu }=0$) and ${W}^{+}={W}^{-}:={W}_{{le}}$. In this equation, it is tacit that $\nu ^{\prime} \equiv \nu ^{\prime} ({R}_{{le}})$ and $\lambda \equiv \lambda ({R}_{{le}})$. The condition $[{T}_{\theta }{]}_{-}^{+}=0$ and Equations (38) and (36) actually lead to an algebraic system of equations to be solved for ${H}_{0}^{{\prime} +}$ and ${V}^{+}$, for instance, and it is not difficult to show that its solution is

$$\begin{eqnarray}\begin{array}{rcl}{W}_{{le}} & = & {\rm{unconstrained}},\\ {V}^{+} & = & \displaystyle \frac{1}{2}\left({R}_{{le}}{V}_{+}^{{\prime} }+{e}^{\lambda ({R}_{{le}})}{W}_{{le}}\right),\\ {H}_{0}^{{\prime} +} & = & {H}_{0}^{{\prime} -}-16\pi {V}^{+}{R}_{{le}}{\tilde{\mu }}^{+}\nu ^{\prime} ({R}_{{le}})+\displaystyle \frac{8\pi }{{R}_{{le}}}{e}^{\lambda ({R}_{{le}})}[\epsilon {]}_{-}^{+}{W}_{{le}}.\end{array}\end{eqnarray} \tag{ 39 }$$

For the solid crust-liquid envelope interface ($r={R}_{{ce}}$), similarly to Equation (39), the relevant boundary conditions are as follows:

$$\begin{eqnarray}\begin{array}{rcl}{T}_{r}^{-} & = & {R}_{{ce}}^{2}[\delta p{]}_{-}^{+}+{W}_{{ce}}{R}_{{ce}}[p^{\prime} {]}_{-}^{+},\\ {T}_{\theta }^{-} & = & 0\to {V}^{-}=\displaystyle \frac{1}{2}\left({R}_{{ce}}{V}_{-}^{{\prime} }+{e}^{\lambda ({R}_{{ce}})}{W}_{{ce}}\right),\\ {H}_{0}^{{\prime} +} & = & {H}_{0}^{{\prime} -}+16\pi {V}^{-}{R}_{{ce}}{\tilde{\mu }}^{-}\nu ^{\prime} ({R}_{{ce}})+\displaystyle \frac{8\pi }{{R}_{{ce}}}{e}^{\lambda ({R}_{{ce}})}[\epsilon {]}_{-}^{+}{W}_{{ce}},\\ {k}^{+} & = & {k}^{-},\end{array}\end{eqnarray} \tag{ 40 }$$

where now the "−" ("+") stands for the top of the solid crust (base of the ocean) and ${W}^{+}={W}^{-}:={W}_{{ce}}={\rm{unconstrained}}$. Given that ${k}^{+}$, $\delta {p}^{+}$ are related to a perfect fluid, it is not difficult to show that the first, third, and fourth equations from Equation (40) result in only one nontrivial condition to be fulfilled.

In summary, we have two free parameters (e.g., ${V}_{+}^{\prime}$ and W_le) at the base of the solid crust fulfilling Equations (39) and two constraints at the top of the solid crust (e.g., ${V}^{-}$ and ${H}_{0}^{{\prime} -}$) compatible with Equations (40). Thus, the problem is well posed, and its solution is unique. For the ocean, the boundary conditions for the perfect-fluid tidal deformation integrations are the continuity of H₀ at the crust-envelope interface and ${H}_{0}^{{\prime} +}$ from Equation (40). We assume in our calculations that the density jump at the interface between the top of the crust and the base of the envelope is negligible.⁶ The density jump at the interface separating the perfect-fluid hadronic phase from the solid crust is also taken as insignificant (Ushomirsky et al. 2000; Chamel & Haensel 2008).

The constraints we have derived cover all sets of equations that one might want to integrate in the solid crust. In our calculations, we choose to work with $(k,{H}_{0},V,W)$. Due to the re-scaling invariance of the equations, solutions for the perturbations should depend on an arbitrary amplitude that naturally does not affect physical quantities like the Love numbers. Given that we have two arbitrary quantities at the base of the solid crust (e.g., W and $V^{\prime}$) to be fixed by two boundary conditions at the top of the crust, a shooting method should be used to integrate the equations in the elastic region. After the integration fulfilling all boundary conditions, one is just left to evaluate Equation (37).

For completeness, we note that the constraints in Equations (39), or Equations (40), are only meaningful if ${\tilde{\mu }}^{+}\ne 0$ (${\tilde{\mu }}^{-}\ne 0$). If this is not the case, then boundary conditions valid for perfect fluids should be used. More specifically, the continuity of the radial traction and the extrinsic curvature at a perfect fluid–perfect fluid (liquid–liquid) interface (at $r={R}_{{ll}}$) imply that

$$\begin{eqnarray}\begin{array}{rcl}{W}_{{ll}} & = & -\displaystyle \frac{{\left[r{\left(\delta p\right)}_{{ll}}\right]}_{-}^{+}}{{\left[{p}_{0}^{{\prime} }\right]}_{-}^{+}}={\left[\displaystyle \frac{{H}_{0}{r}^{3}}{2\left(m+4\pi {r}^{3}p\right)}{e}^{-\lambda }\right]}_{r={R}_{{ll}}},\\ {\left({H}_{0}^{{\prime} +}\right)}_{{ll}} & = & {\left({H}_{0}^{{\prime} -}\right)}_{{ll}}+8\pi {e}^{\lambda ({R}_{{ll}})}\displaystyle \frac{{W}_{{ll}}}{{R}_{{ll}}}[\epsilon {]}_{-}^{+}.\end{array}\end{eqnarray} \tag{ 41 }$$

We use this constraint on ${({H}_{0}^{{\prime} })}_{{ll}}$ for any tidal deformation calculation involving liquid–liquid interfaces in hybrid stars. In these cases, naturally, we solve the the perfect-fluid equation for H₀ as in Hinderer (2008).

In our analysis, we have taken $\kappa =1.5\times {10}^{-2}$ and ${\tilde{\mu }}_{0}=0$ (physically, ${\tilde{\mu }}_{0}\lesssim {10}^{-25}$ cm⁻²). These parameters are reasonable, given that the shear moduli should be around 1% of the crust pressure (Chamel & Haensel 2008; Cutler et al. 2003). As a first cross-check of our equations, we have applied them to the case investigated in Penner et al. (2011). We have worked with the coupled system of equations for $({H}_{0},k,W,V)$ and applied the standard fourth-order Runge–Kutta method. The number of steps used in the solid crust integration is 10⁶, and we have kept an accuracy of 12 significant digits (and a much higher precision). The values obtained for k₂ in the case of elastic crusts is smaller than that in the perfect-fluid case, and relative changes were ${10}^{-6}\mbox{--}{10}^{-3}$, one order of magnitude smaller than the results quoted in Penner et al. (2011). This is in full agreement with the approach of Gittins et al. (2020). We have varied the value of $\tilde{\mu }$ in the solid crust and have found that indeed the solution to H₀ converges to the perfect-fluid case when shear stresses go to zero. This shows that the system of equations we have chosen does not present numerical instabilities.

We have solved the same system of equations for the models presented in Table 1, taking some representative values for the free parameters. Figure 1 shows some of the corresponding M–R relations. We stress that the HS1 model is simply an example of a hybrid EOS and has been used because it is convenient for numerical checks and is expected to give us the main results of the problem of tidal deformations of hybrid stars with elastic crusts. We have chosen the HS1 parameters so that the largest masses of the hybrid stars are above two solar masses, and the tidal deformations of 1.4 M_⊙ stars are smaller than the upper limits found by LIGO/VIRGO (Abbott et al. 2017c, 2019b). The same has been done for the parameters of the somewhat more realistic HS2 and HS3 models (see Table 1).⁷ The masses marking the appearance of the quark phase (the cusps of the mass–radius (M(R)) relations) also vary, reaching larger values for the HS3 EOSs (Sieniawska et al. 2019). In our study, we focus on the cases where $\partial M/\partial {\epsilon }_{c}\gt 0$, where _c is the central density of a family of solutions for an EOS. We do not enter into the details of the stability of hybrid stars with elastic phases but rather assume that the conditions valid for perfect fluids are also valid here.

Zoom In Zoom Out Reset image size

Figure 2 shows the expected relative tidal deformation changes $[({{\rm{\Lambda }}}_{{\rm{perf}}}-{{\rm{\Lambda }}}_{{\rm{sol}}})/{{\rm{\Lambda }}}_{{\rm{perf}}}]$ for 1.4M_⊙ hybrid stars (${{\rm{\Lambda }}}_{{\rm{perf}}}^{1.4}\approx 640$) with large quark phases (∼80% of the star's radius) and small elastic crusts (∼10% of R). We use the HS1 model for these calculations, given its flexibility with density jumps. Relative changes are very small, of the order of 10⁻³–10⁻²%. We find a weak dependence of the relative tidal deformation change on the density discontinuity for values smaller than the critical density jump, above which the whole hadronic phase is elastic. For the parameters chosen, η_crit ≈ 0.91. The sharp change at this density appears as the boundary conditions for the tidal perturbations are changed due to the absence of a perfect-fluid layer (of hadronic matter) separating the quark phase from the elastic crust. The above-mentioned weak dependence is mostly due to the fact that a density discontinuity is accompanied by a change of radius and mass of the star, which almost counteracts the tidal deformation changes.

Zoom In Zoom Out Reset image size

Table 2 shows some consequences (all satisfying the condition $\partial M/\partial {\epsilon }_{c}\gt 0$) of the models of Table 1, for instance, their relative tidal deformation changes with respect to perfect-fluid hybrid stars. We have assumed that the shear modulus of all EOSs is given by Equation (32), and we have taken κ and ${\tilde{\mu }}_{0}$ to be the same as in the beginning of this section. We find that for all models, the relative tidal deformation changes are of the order of 2%–4% for solid hadronic phases larger than ∼60% of the star's radius. For the other cases, the changes are much smaller. We conclude that, in addition to the EOS itself, relative tidal deformation changes depend mostly on the sizes of the elastic phases. As expected, $({\rm{\Delta }}{\rm{\Lambda }}/{{\rm{\Lambda }}}_{{\rm{perf}}})$ decreases with the decrease of the thickness of the elastic crust. We comment on aspects of the last column of Table 2 in the following.

Table 2.  Hybrid NS Properties for the Models of Table 1

Hybrid	η	c	M	Rquark	R	${\rm{\Delta }}{R}_{{\rm{had}}}^{{\rm{sol}}}/R$	${{\rm{\Lambda }}}_{{\rm{perf}}}$	${\rm{\Delta }}{\rm{\Lambda }}/{{\rm{\Lambda }}}_{{\rm{perf}}}$	${\left({\rm{\Delta }}{\rm{\Lambda }}/{{\rm{\Lambda }}}_{{\rm{perf}}}\right)}_{{\rm{had}}}$
Model		$({\epsilon }_{{\rm{sat}}})$	(M⊙)	(km)	(km)	%
HS1	0.0	1.530	0.43	1.06	11.71	38.1	2.558 × 105	2.32 × 10−3	3.56 × 10−2
''	0.16	1.485	0.37	1.19	11.90	45.0	5.984 × 105	5.25 × 10−3	3.56 × 10−2
''	0.30	1.46	0.33	1.02	11.86	51.5	1.148 × 106	9.92 × 10−3	3.56 × 10−2
''	0.45	1.9	0.95	9.28	11.84	16.6	4444.43	9.84 × 10−5	3.59 × 10−3
''	0.45	2.2	1.27	10.36	12.30	12.0	1041.92	4.29 × 10−5	1.89 × 10−3
''	0.45	2.358	1.40	10.69	12.43	10.7	629.5	3.25 × 10−5	1.53 × 10−3
''	0.45	3.0	1.73	11.29	12.62	8.0	178.61	1.65 × 10−5	9.26 × 10−4
''	0.45	5.0	2.02	11.33	12.27	5.8	42.24	7.57 × 10−6	5.74 × 10−4
''	0.60	1.442	0.27	2.25	11.90	64.6	3.314 × 106	2.34 × 10−2	3.47 × 10−2
''	0.70	''	0.25	2.69	11.80	67.7	4.392 × 106	2.63 × 10−2	3.39 × 10−2
''	0.80	''	0.23	3.00	11.67	70.1	5.695 × 106	2.83 × 10−2	3.33 × 10−2
''	0.90	''	0.22	3.03	11.61	73.9	7.223 × 106	3.07 × 10−2	3.56 × 10−2
''	0.91	''	0.22	3.26	11.52	71.7	7.388 × 106	3.25 × 10−2	3.56 × 10−2
''	1.0	''	0.20	3.43	11.36	69.8	8.96 × 106	3.19 × 10−2	3.19 × 10−2
''	1.5	''	0.16	3.99	11.47	61.9	1.924 × 107	2.81 × 10−2	2.81 × 10−2
HS2	0.45	1.625	0.62	1.68	14.41	43.2	8.019 × 104	3.81 × 10−3	3.57 × 10−2
''	0.45	2.1	1.09	8.49	13.29	22.0	2593.44	3.30 × 10−4	1.13 × 10−2
''	0.45	2.483	1.40	9.7	13.24	16.3	643.21	1.30 × 10−4	6.46 × 10−3
''	0.45	3.1	1.69	10.43	13.13	12.5	196.28	6.86 × 10−5	4.03 × 10−3
''	0.45	5.1	1.98	10.68	12.55	9.0	45.282	3.36 × 10−5	2.42 × 10−3
HS3	0.77	3.925	1.12	2.71	12.32	14.1	1607.36	2.28 × 10−5	3.50 × 10−2
''	0.77	4.3	1.30	6.05	11.40	10.6	351.49	8.25 × 10−5	2.25 × 10−2
''	0.77	4.456	1.40	6.64	11.22	9.3	195.47	6.42 × 10−6	1.86 × 10−2
''	0.77	4.6	1.49	7.03	11.09	8.4	123.99	5.88 × 10−6	1.60 × 10−2
''	0.77	5.0	1.68	7.71	10.88	6.7	48.033	1.20 × 10−7	1.14 × 10−2
''	0.79	4.01	1.32	3.26	12.64	11.9	745.94	1.51 × 10−5	3.46 × 10−2
''	0.81	4.11	1.55	3.65	12.99	10.0	356.99	7.57 × 10−6	3.45 × 10−2
''	0.83	4.359	1.88	4.62	13.43	7.8	122.55	5.55 × 10−6	3.37 × 10−2
''	0.83	5.2	2.01	7.09	12.25	6.0	32.229	3.18 × 10−6	2.42 × 10−2

Note. Solid (elastic) phases start at $2\times {10}^{14}\,{\rm{g}}$ cm⁻³ or at the base of the hadronic phase (related to the relative tidal deformation changes of the last column) and end at 10⁷ g cm⁻³ (for more explanation of the models, see the text). In our notation, ${\rm{\Delta }}{\rm{\Lambda }}:={{\rm{\Lambda }}}_{{\rm{perf}}}-{{\rm{\Lambda }}}_{{\rm{sol}}}$. The radius of the quark phase is denoted by R_quark whereas ${\rm{\Delta }}{R}_{{\rm{had}}}^{{\rm{sol}}}:=R-{r}_{\epsilon =2\times {10}^{14}}$ is the thickness of the elastic phase. For some η and the HS1, HS2, and HS3 models, we have calculated relative tidal deformation changes for different masses (central densities, _c) in order to check their dependencies on other quantities. For clarity, they are split by horizontal lines. All results for tidal deformations are related to $\tilde{\mu }$ as in Equation (32), with ${\tilde{\mu }}_{0}=0$ and $\kappa =1.5\times {10}^{-2}$. For the HS3 model, we have used the same parameters of Figure 1.

Download table as:  ASCIITypeset image

Figure 3 shows some relative tidal deformation changes for selected η of Table 2 as a function of the thickness of the elastic hadronic phase. We find that (different) cubic fits (of the form ${ax}+{{bx}}^{2}+{{cx}}^{3}$) could describe $({\rm{\Delta }}{\rm{\Lambda }}/{{\rm{\Lambda }}}_{{\rm{perf}}})$ as a function of $({\rm{\Delta }}{R}_{{\rm{had}}}^{{\rm{sol}}}/R)$ for all hybrid models. We note that the fits agree with the expected condition that ${\rm{\Delta }}{\rm{\Lambda }}\to 0$ when the thickness of the elastic crust goes to zero. We have also found similar dependencies of the relative tidal deformation change as a function of $(1/C){\rm{\Delta }}{R}_{{\rm{had}}}^{{\rm{sol}}}/R$ for each model, where we recall that C is the compactness.

Zoom In Zoom Out Reset image size

Table 2 also explores the cases where the elastic hadronic phase comprises a significant portion of a hybrid star. In order to consider this case, we have assumed a nonzero shear modulus (around 1%–2% of the local pressure) right from the bottom of the hadronic phase. The results for the relative tidal deformation changes are given in the last column of this table. At this level, these results are merely indicative, but they are motivated by the presence of an elastic mixed phase in a hybrid star (Sotani et al. 2013). In this case, we find that very different models (for example, the HS1, HS2, and HS3 models) roughly agree with the "largest" relative tidal deformation changes, which are also around 3%–4% when the solid hadronic region comprises over 75% of the star's radius.

When realistic shear moduli are considered, the results for the tidal deformation changes might change. For instance, the SLy4 shear modulus (Haskell et al. 2006; Sotani et al. 2007; Zdunik et al. 2008) is around one to five times smaller than that of Equation (32) in the range where the hadronic phase is elastic. This means, roughly speaking, that similar changes are expected for the tidal deformation quantities that we calculated before. Indeed, numerical calculations with the SLy4 shear modulus for all of the models of Table 1 lead to this level of decrease in $({\rm{\Delta }}{\rm{\Lambda }}/{{\rm{\Lambda }}}_{{\rm{perf}}})$. Therefore, relative tidal deformation changes of hybrid NSs are not expected to be larger than 3%–4% of their perfect-fluid values.

5.1. Phase-splitting Surface Quantities

An interesting issue in hybrid stars is the possibility of inducing, due to perturbations, surface quantities on the phase-splitting interfaces (see, e.g., Pereira & Lugones 2019). Surface quantities induced by perturbations could significantly change the boundary conditions, which, in turn, may change the tidal deformations of elastic stars. We show in this section, nonetheless, that the Einstein equations for quasi-static perturbations do not allow the induction of surface quantities at phase-splitting surfaces.

From the Darmois-Israel formalism (see, e.g., Pereira et al. 2014), we know that surface quantities are, in general, necessary for the consistent match of two spacetimes, and they are related to the presence of a surface energy-momentum tensor S^a_b of the form (Pereira et al. 2014)

$$\begin{eqnarray}&&{S}_{b}^{a}=-\displaystyle \frac{1}{8\pi }[{K}_{b}^{a}-{\delta }_{b}^{a}K{]}_{-}^{+},\end{eqnarray} \tag{ 42 }$$

where we recall that K^a_b is the extrinsic curvature and $K={K}_{a}^{a}={K}^{{ab}}{\bar{h}}_{{ab}}$ (${\bar{h}}_{{ab}}$ is the induced metric on the hypersurface splitting two given spacetimes) is its trace. From Equations (34) and (35), it is easy to see that S^a_b should be related to a perfect fluid, so that one could write

$$\begin{eqnarray}&&{S}_{{ab}}=(\sigma +{ \mathcal P }){u}_{a}{u}_{b}+{ \mathcal P }{\bar{h}}_{{ab}},\end{eqnarray} \tag{ 43 }$$

where

$$\begin{eqnarray}&&\sigma =-\displaystyle \frac{1}{4\pi }[{K}_{\theta }^{\theta }{]}_{-}^{+}\end{eqnarray} \tag{ 44 }$$

and

$$\begin{eqnarray}&&{ \mathcal P }=\displaystyle \frac{1}{8\pi }[{K}_{t}^{t}+{K}_{\theta }^{\theta }{]}_{-}^{+}=\displaystyle \frac{1}{8\pi }[{K}_{t}^{t}{]}_{-}^{+}-\displaystyle \frac{\sigma }{2}.\end{eqnarray} \tag{ 45 }$$

The physical reason for the perfect-fluid nature of the surface energy-momentum tensor is the axial symmetry of the spacetimes (perturbed core and crust).

Equation (19) allows us to easily find $[k^{\prime} {]}_{-}^{+}$, which, from Equation (44), relates $[{H}_{0}^{{\prime} }{]}_{-}^{+}$ to σ. The discontinuity of ${H}_{0}^{{\prime} }$, $[{H}_{0}^{{\prime} }{]}_{-}^{+}$, on the other hand, can be determined by Equations (45) and (34). When one solves the associated equation, it follows that ${ \mathcal P }=0$. Since ${ \mathcal P }$ is a monotonic function of σ in ordinary cases, we have σ = 0, too. Thus, the intrinsic curvature components should always be continuous at any phase-transition interface.

We have investigated tidal deformations of hybrid stars with elastic hadronic phases, building on previous works and extending them to the case of sharp phase transitions, which is a possibility for hybrid stars. In the case of a shear modulus around 1% of the hadronic pressure, stars with small elastic hadronic phases (large perfect-fluid quark phases) lead to negligible changes in Λ. However, there are cases where the tidal deformation changes with respect to perfect fluids could be as large as 2%–4%. These cases relate to elastic hadronic regions larger than ∼60% of the stellar radius. This would imply, for instance, small quark phases and configurations with large enough density jumps at the interface separating the quark phase from the hadronic phase or even small quark phases with lower phase-transition pressures. One would expect these cases to maximize the relative tidal deformation changes because they are exactly the ones where a solid phase would comprise the largest fraction of the star.

Although ≃5% changes in the tidal deformability will be difficult to register for current GW detectors, for signals emitted from typical distances of 100 Mpc, the issue of detecting the elastic hadronic regions of stars may become relevant for future detectors, such as the Einstein Telescope (Punturo et al. 2010). Indeed, in the recent review of the Einstein's Telescope science case (Maggiore et al. 2020), the expected increase of sensitivity translates into more than one order-of-magnitude larger values for the signal-to-noise. In the limit of high signal-to-noise, the Fisher-information matrix approximation suggests that the errors on the measured quantities scale inversely proportional to the signal-to-noise (Vallisneri 2008; Chatziioannou et al. 2017; Jiménez Forteza et al. 2018). For example, for a GW170817-like event observed by the Einstein Telescope, the errors of the measured component mass-weighted tidal deformability $\tilde{{\rm{\Lambda }}}$ (equal to ${300}_{-230}^{+420}$ in the low-spin prior case; Abbott et al. 2019b) would be sufficiently small to allow for studies of the elastic properties.

We have also found that there is a weak dependence of the tidal deformations of solid stars on density jumps. The reason seems to be a competition between larger-density discontinuities and smaller masses (radii) for stars with a given central density. However, a sudden change might appear whenever there is not a layer of perfect fluid separating the quark phase from the solid crust. This might happen when density jumps are large or the critical density favoring dense matter to arrange as a lattice is high enough. We recall that this critical density is not known (it is EOS dependent), but it is believed to be smaller than the nuclear saturation density (Ducoin et al. 2011). This means that one might expect the presence of a layer of hadronic perfect-fluid matter separating a quark phase from the elastic crust even in the case of sharp phase transitions. Thus, from Figure 2, one would expect that the relative tidal deformation changes are very small for systems with large perfect-fluid regions.

Hybrid stars where an elastic hadronic phase would start right after the quark phase (and would extend until the ocean) could also be seen as a simplistic model for a solid star with an elastic mixed phase. Although we have not made any direct calculations in this regard, our numerical integrations with nonzero shear stresses right after the quark phase indicate that elastic mixed phases could lead to non-negligible tidal deformation changes (around 2%–4% or even larger). One would expect this, given that the shear modulus of a mixed phase might be much larger than the crustal shear modulus, still roughly 1%–2% of the local pressure (Sotani et al. 2013). However, there is an important, though subtle, point in this case. In the presence of a mixed phase, one would expect physical quantities such as the energy density and shear modulus to be continuous anywhere in the star. If one takes the quark phase as a perfect fluid, then $\tilde{\mu }$ should go continuously from zero to a nonzero value inside the mixed phase. This would render equations ill-defined at the quark-mixed phase interface since there are terms of the form $1/\tilde{\mu }$ in the perturbation equations. Clearly, the main source of the problem is the perfect-fluid approximation. If one takes the shear modulus to be small but finite, the tidal deformation equations for the mixed phase would be well defined. We leave this issue to be investigated elsewhere. In this case, boundary conditions for the problem change, which could have a non-negligible effect on tidal deformation changes when compared to perfect fluids. It would also be of interest to investigate the case of dynamical tides (see, e.g., Andersson & Pnigouras 2019, 2020; Schmidt & Hinderer 2019 and references therein). Although they would change tidal deformations by some percent in the case of perfect fluids (Andersson & Pnigouras 2020), dynamical tides might have a larger impact on hybrid stars with solid phases due to aspects of the fundamental mode.

We have explicitly shown that surface degrees of freedom cannot be induced on phase-splitting surfaces. This is expected due to the assumption of no reactions around phase transitions and the static nature of the problem. It implies that there would be no way of significantly changing the tidal deformations of solid crusts from the results we have found when the Einstein equations are taken into account. However, surface degrees of freedom could play an important role in more exotic cases with background degrees of freedom, such as boson stars or gravstars (Johnson-McDaniel et al. 2018).

Our results suggest that when the accuracy of tidal deformations is around 5%, even in the most conservative cases, shear stresses should not be disregarded. In the worst case, they would be an important part of systematic uncertainties for some models. It would be of interest to advance the analysis proposed here for more realistic hybrid EOSs in order to single out particularities lost in our idealized approach. Most importantly, one should also consistently calculate the shear modulus for a given EOS.

We acknowledge the constructive discussions with Leszek Zdunik, Paweł Haensel, Morgane Fortin, and Lami Souleiman. This work was partially supported by the Polish National Science Centre grant No. 2016/22/E/ST9/00037. J.P.P. is also thankful for the partial support given by FAPESP (grant Nos. 2015/04174-9 and 2017/21384-2). N.A. acknowledges financial support from STFC via grant No. ST/R00045X/1.