Astrophysical aspects of general relativistic mass twin stars

TL;DR

This paper explores the astrophysical implications of mass twin stars, which are compact stars with similar masses but different radii due to phase transitions, using an effective polytropic EoS model.

Contribution

It introduces a polytropic EoS model to study mass twins and analyzes their properties in both slow rotation and full general relativity frameworks.

Findings

Mass twins can exist with a strong first-order phase transition.

A third stable branch appears in the mass-radius diagram.

Constraints on maximum mass from GW170817 are discussed.

Abstract

In this chapter we will introduce an effective equation of state (EoS) model based on polytropes that serves to study the so called "mass twins" scenario, where two compact stars have approximately the same mass but (significant for observation) quite different radii. Stellar mass twin configurations are obtained if a strong first-order phase transition occurs in the interior of a compact star. In the mass-radius diagram of compact stars, this will lead to a third branch of gravitationally stable stars with features that are very distinctive from those of white dwarfs and neutron stars. We discuss rotating hybrid star sequences in the slow rotation approximation and in full general relativity and draw conclusions for an upper limit on the maximum mass of nonrotating compact stars that has recently be deduced from the observation of the merger event GW170817.

Equations

$$\frac{\Delta \epsilon }{{\epsilon }_c}\ge \frac{1}{2}+\frac{3}{2}\frac{P_c}{{\epsilon }_c}.$$

$$d{s}^2$$

$$T_{\mu }^{\nu }=(\epsilon +P)u_{\mu }{u}^{\nu }-P{\delta }_{\mu }^{\nu },$$

$$G_r^r-G_t^t$$

$$G_{\theta }^{\theta }+G_{\varphi }^{\varphi }$$

$$G_{\theta }^r$$

$$G_{\varphi }^t=8\pi G{T}_{\varphi }^t.$$

$$H(r,\theta )\equiv \int \frac{d{P}^{\prime }}{P^{\prime }+{\epsilon }^{\prime }}=\frac{1}{2}\ln [u^t(r,\theta )]+\mathrm{const},$$

$$d{s}^2$$

$$M$$

$$M_0$$

$$J$$

$$T$$

$$R_e={\kappa }^{1/2}{\bar{r}}_e{e}^{({\gamma }_e-{\rho }_e)/2},$$

$$Z_p=e^{-({\gamma }_p+{\rho }_p)/2}-1,$$

$$Z_f={\left(\frac{1-v_e}{1+v_e}\right)}^{1/2}\frac{e^{-({\gamma }_e+{\rho }_e)/2}}{1+{\hat{\omega }}_e{e}^{-{\rho }_e}}-1,$$

$$Z_b={\left(\frac{1+v_e}{1-v_e}\right)}^{1/2}\frac{e^{-({\gamma }_e+{\rho }_e)/2}}{1-{\hat{\omega }}_e{e}^{-{\rho }_e}}-1.$$

$$g_{\mu \nu }(r,\theta )=\sum _{j=0}^{\infty }{\left(\sqrt{\beta }\right)}^j{g}_{\mu \nu }^{(j)}(r,\theta ).$$

$$\displaystyle\begin{array}[]{ccc}e^{-\lambda(r,\theta)}&=&e^{-\lambda^{(0)}(r)}[1+\beta f(r,\theta)]+O(\beta^{2})\,,\\ e^{\nu(r,\theta)}&=&e^{\nu^{(0)}(r)}[1+\beta\Phi(r,\theta)]+O(\beta^{2})\,,\\ e^{\mu(r,\theta)}&=&r^{2}[1+\beta U(r,\theta)]+O(\beta^{2})\,,\end{array}$$

$$\displaystyle\begin{array}[]{ccc}\omega(r,\theta)&=&\sqrt{\beta}q(r,\theta)+O((\sqrt{\beta})^{3})\,.\end{array}$$

$$P(r,\theta )$$

$$\epsilon (r,\theta )$$

$$g_{\mu \mu }^{(j)}(r,\theta )$$

$${\lambda }^{(0)}(r)$$

$${\nu }^{(0)}(r)$$

$$\frac{d{P}^{(0)}(r)}{dr}=-G[P^{(0)}(r)+{\epsilon }^{(0)}(r)]\frac{m(r)+4\pi P^{(0)}(r)r^3}{r[r-2Gm(r)]},$$

$$m(r)=4\pi {\int }_0^r{\epsilon }^{(0)}(r^{\prime })r^{\prime 2}d{r}^{\prime },$$

$$M=m(R)=4\pi {\int }_0^R{\epsilon }^{(0)}(r^{\prime })r^{\prime 2}d{r}^{\prime },$$

$${\lambda }^{(0)}(r)$$

$${\nu }^{(0)}(r)$$

$$\frac{d{N}_B(r)}{dr}$$

$$d{s}^2$$

$$\frac{dA(r)}{dr}$$

Full text

Topics on Strong Gravity

Chapter 1 Astrophysical aspects of general relativistic mass twin stars

David Blaschke1,2,3, David Edwin Alvarez-Castillo1,

Alexander Ayriyan4,5, Hovik Grigorian4,5,6, Noshad Khosravi Lagarni1,7,

Fridolin Weber8,9

1* Bogoliubov Laboratory for Theoretical Physics, Joint Institute for Nuclear Research, Joliot-Curie street 6, 141980 Dubna, Russia

2 Institute of Theoretical Physics, University of Wroclaw, Max Born place 9, 50-204 Wroclaw, Poland

3 National Research Nuclear University (MEPhI), Kashirskoe Shosse 31, 115409 Moscow, Russia

4 Laboratory for Information Technologies, Joint Institute for Nuclear Research, Joliot-Curie street 6, 141980 Dubna, Russia

5 Computational Physics and IT Division, A.I. Alikhanyan National Science Laboratory, Alikhanyan Brothers street 2, 0036 Yerevan, Armenia

6 Department of Physics, Yerevan State University, Alek Manukyan street 1, 0025 Yerevan, Armenia

7 Department of Physics, Alzahra University, Tehran, 1993893973, Iran

8 Department of Physics, San Diego State University, 5500 Campanile Drive, San Diego, CA 92182, USA

9 Center for Astrophysics and Space Sciences, University of California, San Diego, La Jolla, CA 92093, USA *

Abstract. In this chapter we will introduce an effective equation of state (EoS) model based on polytropes that serves to study the so called ”mass twins” scenario, where two compact stars have approximately the same mass but (significant for observation) quite different radii. Stellar mass twin configurations are obtained if a strong first-order phase transition occurs in the interior of a compact star. In the mass-radius diagram of compact stars, this will lead to a third branch of gravitationally stable stars with features that are very distinctive from those of white dwarfs and neutron stars. We discuss rotating hybrid star sequences in the slow rotation approximation and in full general relativity and draw conclusions for an upper limit on the maximum mass of nonrotating compact stars that has recently be deduced from the observation of the merger event GW170817.

1 Introduction

Compact stars, the stellar remnants following the death of main sequence stars, have been the subject of investigation since the beginning of the last century. In particular, the determination of the internal composition of neutron stars is an open problem. Researching it involves many areas of physics, like nuclear, plasma, particle physics and relativistic astrophysics. Moreover, due to the enormous compactness (as expressed in the mass-radius ratio) of compact stars, these objects are extremely relativistic. Therefore, one can neither exclusively apply non-relativistic quantum mechanics nor classical Newtonian gravity to describe the observational properties of compact stars.

During the last decade important astronomical observations have shed light onto the nature of the dense, cold matter in the stellar interiors of compact stars. The detection of massive neutron stars, of about $2M_{\odot}$, has constrained the maximum density values in their cores and also revealed the stiff nature of the nuclear equation of state (EoS) at ultra-high densities. Strongly related to this issue, and one of the most interesting aspects of modern dense-matter physics, concerns the possible onset of quark deconfinement in the cores of compact stars.

Microscopic models that take into account the nuclear interactions either at the nucleon or quark level aim at providing a realistic hadronic or quark matter EoS, respectively. Neutron star matter must be thermodynamically consistent. Interestingly, due to the fast cooling of neutron stars after their birth in a supernova collapse the thermal contributions to the EoS do not contribute substantially and can safely be neglected Yakovlev et al. (2001). The thermodynamic system can therefore be described by three macroscopic variables: energy density $\varepsilon$, baryonic density $n$, and pressure $P$. A fourth quantity of great interest, the chemical potential, can then be obtained as $\mu=(P+\varepsilon)/n$.

In addition, the most basic conditions that the system must fulfill include global charge neutrality and $\beta$-equilibrium. The latter is derived from the reaction balance of beta decay and its inverse, the electron capture, due to the weak interactions and fixes the relation between the chemical potentials of different species in the system.

With the above conditions satisfied, the neutron star equation of state becomes an expression of the form $P(\varepsilon)$, where $\varepsilon=\varepsilon(n)$ and $P=P(n)$ acquire parametric forms. Furthermore, in order to compute the internal properties of compact stars, it is necessary to obtain internal pressure profiles. This will result in mass-radius relations that characterize an EoS. Neutron stars are extremely relativistic objects which require to be treated within Einstein’s general theory of relativity rather than simply Newtonian gravity, which may still be applicable to white dwarf stars. In this sense, our contribution addresses ”strong gravity” in a unique fashion. To give an example, for a pulsar of mass $M=2~{}M_{\odot}$ with a typical radius of around 12 km, the general relativistic correction factor amounts to $1/(1-2GM/R)=2$ Tolman (1939); Oppenheimer and Volkoff (1939), which is a 100% correction relative to Newtonian gravity!

In the following sections we will introduce an effective EoS model based on polytropes (Alvarez-Castillo and Blaschke, 2017) that serves to study the so called ”mass twins” scenario, where two compact stars have approximately the same mass but (significant for observation) quite different radii Glendenning and Kettner (2000). Stellar mass twin configurations are obtained if a strong first-order phase transition occurs in the interior of a compact star. In the mass-radius diagram of compact stars, this will lead to a third branch of gravitationally stable stars with features that are very distinctive from those of white dwarfs and neutron stars. The condition on the EoS that will lead to mass twins was first derived by Seidov (Seidov, 1971), see also (Schaeffer et al., 1983; Lindblom, 1998), namely that the central energy density $\varepsilon_{v}$, central pressure $P_{c}$, and the jump in energy associated with the phase transition $\Delta\varepsilon$ obey the relation

[TABLE]

When fulfilled, the corresponding compact star will suffer an instability of the same type as the maximum mass star of a stellar sequence. Most interestingly, stars with central densities higher than the density of the maximum-mass star become stable again if their gravitational masses obey $\partial M/\partial\varepsilon_{c}(0)>0$, thereby populating a third branch with stable mass twins.

2 Self-consistent set of field equations for stationary rotating and

tidally deformed stars

The geometrical description of the space-time structure curved due to the mass -energy of a compact star is given by the general metric form defining the interval between the infinitesimally close events,

[TABLE]

The curvature of the space-time is satisfying the Einstein field equations $G_{\nu}^{\mu}=8\pi GT_{\nu}^{\mu},$where $G_{\nu}^{\mu}=R_{\nu}^{\mu}-1/2\delta_{\nu}^{\mu}R$. Here $R_{\nu}^{\mu}$ is the Ricci curvature tensor and $R$ the scalar curvature. On the right hand side of the Einstein equation we have the energy -momentum tensor $T_{\nu}^{\mu}$ of the stellar matter and $G$ is the gravitational constant ($\hbar=c=1)$.

The metric tensor $g_{\mu\nu}(x)$ has the same symmetry as the matter distribution. Therefore, if one assumes that the star is static and not deformed the metric tensor is diagonal and depends only on the distance from the center of the star. In the case of stationary rotating stars the symmetry of the matter will be axial symmetric. In this case due to the rotational motion of the star the non-diagonal element $g_{t\phi}$ ($t$-time coordinate, $\phi$ - azimutal angle of the spherical coordinate system) will not be zero in the inertial frames connected with the star. The existence of such a term leads to the Lense-Thirring effect of frame dragging for the motion of bodies in the gravitational field of rotating compact relativistic stellar objects. However, in the case of small but static deformations the metric will be non-spherical but diagonal.

2.1 Einstein equations for axial symmetry

The general form of the metric for an axially symmetric space-time manifold in the inertial frame where the star center is at rest is

[TABLE]

where a spherically symmetric coordinate system has been used in order to obtain the Schwarzschild solution as a limiting case. This line element is time - translation and axial-rotational invariant; all metric functions are dependent on the coordinate distance from the coordinate center $r$ and altitude angle $\theta$ between the radius vector and the axis of symmetry.

The energy momentum tensor of stellar matter can be approximated by the expression of the energy momentum tensor of an ideal fluid

[TABLE]

where $u^{\mu}$ is the $4$-velocity of matter, $P$ the pressure and $\varepsilon$ the energy density.

Once the energy-momentum tensor (4) is fixed by the choice of the equation of state for stellar matter, the unknown metric functions $\nu$,$\lambda$, $\mu$, $\bar{\omega}$ can be determined by the set of Einstein field equations for which we use the following four combinations.

There are three Einstein equations for the determination of the diagonal elements of the metric tensor,

[TABLE]

and one for the determination of the non-diagonal element

[TABLE]

We use also one equation for the hydrodynamic equilibrium (Euler equation)

[TABLE]

where the gravitational enthalpy ${H}$ thus introduced is a function of the energy and/or pressure distribution.

2.2 Full solution for uniform rotational bodies

In this section we describe the method of solution employed by the RNS code written by Stergioulas and Friedman (1995), based on the method developed by Komatsu et al. (1989) that also includes modifications by Cook et al. (1994). In addition, the inclusion of quadrupole moments is due to Morsink based on the method by Laarakkers and Poisson (1999).

In order to study the full solutions for rotating compact stars the following metric is considered Cook et al. (1994):

[TABLE]

which just like Eq. (2.1) properly describes a stationary, axisymmetric spacetime. In addition, the matter source is chosen to be a perfect fluid described by the stress energy tensor $T^{\mu\nu}=(\epsilon+P)u^{\mu}u^{\nu}+Pg^{\mu\nu}$, where $u^{\mu}$ is the four-velocity of matter. Three of the solutions to the gravitational field equations are found by a Green function approach therefore leading to the determination of the metric potentials $\rho$, $\gamma$ and $\omega$ in term of integrals, whereas the $\alpha$ potential is found by solving a linear differential equation. Therefore, we find the corresponding numerical solutions to our compact star models by employing the RNS code. In the formulation of the problem, all the physical variables are written in dimensionless form by means of a fundamental length scale $\sqrt{\kappa}$, where $\kappa\equiv\frac{c^{2}}{G\epsilon_{0}}$ with $\epsilon_{0}\equiv 10^{15}$g cm*-3*.

The global parameters of the star are computed by means of the following expressions:

[TABLE]

where $\hat{\omega}\equiv\bar{r}_{e}\bar{\omega}$ and $\hat{\Omega}\equiv\bar{r}_{e}\bar{\Omega}$ with $\bar{r}_{e}$ as the coordinate radius of the equator. Moreover, all the resulting quantities can be written in terms of the auxiliary variables $\mu$ and s, defined as $\mu\equiv\theta$ and $\bar{r}\equiv\bar{r}_{e}\left(\frac{s}{1-s}\right)$, respectively. Consequently, the four metric functions acquire a dependence on the above variables: $\rho(s,\mu)$, $\gamma(s,\mu)$, $\omega(s,\mu)$, $\alpha(s,\mu)$. The remaining quantities are:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where the subscripts $e$ and $p$ denote evaluation at the equation and at the pole, respectively.

For the solutions of maximally rotating compact stars in numerical general relativity the version of RNS code has been employed which was available for download by the time of the writing of this contribution from the website http://www.gravity.phys.uwm.edu/rns/

2.3 Perturbation approach to the solution

The problem of the rotation can be solved iteratively by using a perturbation expansion of the metric tensor and the physical quantities in a Taylor series with respect to a small positive parameter $\beta$. As such a parameter for the perturbation expansion we use a dimensionless quantity. One of possible physically motivated way is to take the ratio of the rotational or deformation energy to the gravitational one. The gravitational energy could be estimated for a homogeneous Newtonian star as $\beta=E_{\rm def}/E_{\rm grav}$. In case of rotating stars the deformability is connected with the induced centrifugal force and the expansion parameter is $\beta=E_{{\rm rot}}/E_{{\rm grav}}=(\Omega/\Omega_{0})^{2}$, where $\Omega_{0}^{2}=4\pi G\rho(0)$ with the mass density $\rho(0)$ at the center of the star. The choice of this parameter could be also motivated with the conditions when the problem is discussed. For example since for the stationary rotating stars can not have too high value of the angular velocity, because of mass shedding on Keplerian angular velocity $\Omega_{K}=\sqrt{GM/R_{e}^{3}}$ for the star with total mass $M$ and $R_{e}$ equatorial radius, the expansion gives sufficiently correct solutions already at $O(\Omega^{2})$. So the expansion parameter is naturally limited to values $\Omega/{\Omega}_{0}\ll 1$ by this condition of mechanical stability of the rigid rotation, because always $\Omega<\Omega_{K}={\Omega}_{0}/\sqrt{3}$. This condition is fulfilled not only for homogeneous Newtonian spherical stars but also for the relativistic configurations even with a possible hadron-quark (deconfinement) transition, which we are going to discuss later in this chapter, see Fig. 6 below. The perturbation approach to slowly rotating stars has been developed first by Hartle (1967); Hartle and Thorne (1968), and independently by Sedrakyan and Chubaryan (1968b, a) and is described in detail in Refs. Weber, F. (1999); Glendenning, N. K. (2000); Chubarian et al. (2000). Our notation and derivation in this section will follow Chubarian et al. (2000) while numerical solutions for the slow rotation ($\Omega^{2}$-) approximation are obtained with a code based on the improved Hartle scheme Weber and Glendenning (1992).

The expansion of the metric tensor in a perturbation series with respect to the slow rotation parameter $\beta$

can be expressed as

[TABLE]

According to the metric form for the axial symmetry in the linear approximation via $\beta$ parameter we introduce the notations describing explicitly the non perturbed (spherically symmetric case noted with upper index $(0)$) and perturbed terms (corresponding to $j=1,2$) in the metric,

[TABLE]

and for the frame dragging frequency $\omega$ the odd orders

[TABLE]

In the same way one performs a velocity expansion of the energy-momentum tensor, the pressure and energy density distributions, and of the kinetic energy,

[TABLE]

where $P^{(0)}$ and $\varepsilon^{(0)}$ denote the zero-order coefficients which correspond to non-deformed spherically symmetric stars.

Because of rotational symmetry the diagonal elements $g_{\mu\mu}^{(j)}(r,\theta)$ (no summation over $\mu$) of the metric coefficients can be written as ($j$ and $l$ are even values only)

[TABLE]

The same is true also for the non diagonal elements, but the angular dependence is different ($j$ and $l$ values are now odd only). The case for the non diagonal term will be investigated in section 3.1, where the moment of inertia will be discussed.

In the case of a deformed distribution of the matter, the external metric has the form of the Kerr metric. However this solution is for black holes and does not correspond to a realistic stellar models. Therefore the external as well the internal solutions of the field and matter distributions can only be obtained either via a perturbation approximation or via a completely numerical treatment.

2.4 Static spherically symmetric star models

As a first step in the perturbation approach of a slowly rotating star one needs to find the internal gravitational field, the mass and matter distributions, the total gravitational mass, the radius and all other characteristic properties (including the metric functions) of spherically symmetric stars. The solution for the metric functions in empty space (i.e., the external solutions) are given by the Schwarzschild solution.

[TABLE]

where $M$ is a constant of integration, which asymptotically is the Newtonian gravitational mass of the object.

These nonlinear equations, however, could be written in an elegant form suggested by Tolman, Oppenheimer and Volkoff, which are known as the TOV equations (Tolman, 1939; Oppenheimer and Volkoff, 1939), and solved in a way such that the internal solution matches the analytic external Schwarzschild solution. The TOV equation is given by (for a derivation, see, e.g., the textbook by Misner et al. (1973))

[TABLE]

where $P^{(0)}$ and $\varepsilon^{(0)}$ denote the equation of state (EoS) describing the stellar matter. The quantity $m(r)$, defined as

[TABLE]

stands for the amount of gravitational mass contained inside a sphere of radius $r$, with $r$ denoting the distance from the center of the star. The star’s total gravitational mass, $M$, is then given

[TABLE]

where $R$ denotes the radius of the star defined by $P(r=R)=0$. Physically, the TOV equation describes the balance of gravitational and internal pressure forces at each radial distance inside the star. Both forces exactly cancel each other inside a static stellar configuration, as described by the TOV equation.

The TOV equation is solved numerically, for a given model for the EoS, by choosing a value for the star’s central density and then integrating Eq. (30) out to a radial location where the pressure becomes zero. So for any fixed choice of the EoS, the stars form a one-parameter sequence (parameter $\varepsilon^{0}_{c}$). An entire family of compact stars is obtained by solving the TOV equation for a range of central densities which result in the mass-radius relationship of compact stars. It is characterized by the existence of a maximum mass star (several maximum mass stars if permitted by the model chosen for the EoS). The stars are stable against gravitational collapse if they are on the stellar branch for which $\partial M/\partial\varepsilon^{(0)}_{c}>0$. Stars on the stellar branch where $\partial M/\partial\varepsilon^{(0)}_{c}<0$ are unstable agains radial oscillations and will therefore not exist stably in the universe.

Each stellar model has unique solutions for $m(r)$, $P^{(0)}(r)$ and $\varepsilon^{(0)}(r)$ in terms of which the internal gravitational field (the metric coefficients) is defined as

[TABLE]

The internal field solutions are smoothly connected to the external field solutions at the stellar surface, $r=R.$ Once the internal pressure profiles are derived from the solution of the TOV equations, it is possible to compute other astrophysical quantities like baryonic mass, and also make the next step in the perturbation approach to define the moment of inertia and tidal deformabilities, which are of very great observational interest.

Of particular interest for astrophysical scenarios (stellar evolution) is the expression for the total baryon mass

[TABLE]

where $n(r)$ is the baryon number density and $N_{B}(R)$ is the total baryon number of the star. This number is a characteristic conserved quantity and is very important in discussions of evolutionary scenarios of compact stars, see, e.g., Bejger et al. (2017); Ayvazyan et al. (2013); Chubarian et al. (2000); Poghosyan et al. (2001). The functions of the spherically symmetric solution in Eqs. (23) and (26) can be found from Eq. (5) and Eq. (9) in zeroth order of the $\Omega$-expansion.

3 Tidal deformability of compact stars

The tidal deformability (TD) is a measure of the shape deformation property of the astrophysical object under the gravitational influence of another nearby object. To determine it in the first order we need to consider a modification of the space-time metric when the distribution of matter of the star becomes elliptic. According to the symmetries of the metric coefficients introduced in Eq. (LABEL:metric-1) we have even orders $j=0,2,\ldots$ for the diagonal elements111Notation corresponds to the work of Chubarian et al. (2000).. The the first correction corresponding to small deformations (small values of $\beta$) one can consider terms linear in $\beta$ or equivalently the $j=2$ perturbation approximation to the spherically symmetric star. We introduce some new notation and work under the assumption that $f_{2}(r)=-\Phi_{2}(r)=A(r)$ like in the expansion of the external solution, since $\nu^{(0)}(r)=-\lambda^{(0)}(r)$.

The non diagonal term could be taken to be zero, because we consider only the static case $q(r,\theta)=0$ (the parameter defining the static deformation $\sqrt{\beta}$ does not change the sign under time reversal $t\rightarrow-t$) and $U_{2}(r)=K(r)$, so that we have

[TABLE]

In this approximation, without loss of generality, one can set the values of $f_{0}=\Phi_{0}=U_{0}=0$, because we neglect the contribution of the deformation energy to the gravitational mass. The equations show that $K^{\prime}(r)=A^{\prime}(r)+{\nu^{(0)}}^{\prime}(r)A(r)$ the prime symbol denoting the derivative of those quantities with respect to $r$. The functions $A(r)$ and $B(r)=dA/dr$ obey the differential equations

[TABLE]

Here, $c_{s}^{2}=dP/d\varepsilon$ is the square of the speed of sound, which is equivalent to the knowledge of the equation of state. The pressure profile provided by solving the TOV equations will complement the above equations.

The system is to be integrated with the asymptotic behavior of metric functions $A(r)=a_{0}r^{2}$ and $B(r)=2a_{0}r$ as $r\to 0$, The $a_{0}$ is a constant that quantifies the deformation of the star which can be taken arbitrary. This constant corresponds to the choice of $\beta$ as an external parameter. Since it cancels in the expression for the Love number and in all other quantities in consideration its value is not important. Using the solution on the surface at $r=R$ and the following combination

[TABLE]

it is possible to compute the $l=2$ Love number (Hinderer, 2008; Damour and Nagar, 2009; Binnington and Poisson, 2009; Yagi and Yunes, 2013; Hinderer et al., 2010):

[TABLE]

where $M/R$ in the expression $C=GM/R$ is the compactness of the star ($2C$ is the ratio of gravitational radius to spherical radius).

The dimensionless tidal deformability parameter is defined as $\Lambda=\lambda/M^{5}$, a quantity defined for small tidal deformabilities. Here $\lambda$ is the TD of the star with a gravitational mass $M$, just as defined above. In addition, the love number is related to TD and defined as

[TABLE]

In the investigations and observations of the process of neutron star merging the TD $\Lambda$ is a key parameter characterising the stiffness of equation of state of the stellar matter.

3.1 Moment of inertia

The moment of inertia is one of the main characteristics of the mechanical properties of the rotating body, therefore one needs to define it also for the relativistic objects as neutron star obeying the gravitational field contribution to the rotational motion. The baryonic mass is an important quantity often associated with explosive events, where it can be conserved while the gravitational mass of the star suffers modifications (Alvarez-Castillo et al., 2015; Bejger et al., 2017). The moment of inertia is related to the glitch phenomenon, which is a sudden spin-up in the general spin-down evolution of rotation frequencies, observed for some pulsars, see (Haskell and Melatos, 2015) and references therein. Moreover, it is expected to be measured in pulsar binaries, providing a strong constraint on the compact star EoS (Lattimer and Prakash, 2007).

In a very simplified way one can estimate the impact of relativistic effects on the moment of inertia from (Ravenhall and Pethick, 1994)

[TABLE]

where $J$ denotes the total, conserved angular moment. The quantity $C_{J}$ is given by

[TABLE]

which can be readily computed since only the knowledge of spherically symmetric quantities is required, which are easy to compute.

However, because of the deformation of the star due to rotation and the impact of the gravitational field on the rotational inertia of the star, the moment of inertia becomes a function of the rotational state, i.e., a function of the angular velocity or spin frequency of the neutron star. To take all these effects into account in the defining expression for the moment of inertia we will follow the steps of our perturbation approach.

By definition, the angular momentum of the star in the case of stationary rotation is a conserved quantity and can be expressed in invariant form

[TABLE]

where $\sqrt{-g}dV$ is the invariant volume and $g=\det||g_{\mu\nu}||$. For the case of slow rotation where the shape deformation of the rotating star can be neglected and using the definition of the moment of inertia $I_{0}(r)=J_{0}(r)/\Omega$ accumulated in the sphere with radius $r$, we obtain from Eq. (43)

[TABLE]

Here $\bar{\omega}$ the difference of the frame dragging frequency $-\omega$ and the angular velocity $\Omega$. In general relativity, due to the Lense-Thirring law, rotational effects are described by

[TABLE]

This expression is approximated from the exact expression of the energy momentum tensor coefficient $T_{\phi}^{t}$ and the metric tensor in the axial symmetric case.

We keep only two non-vanishing components of the 4-velocity

[TABLE]

because we assume that the star due to high viscosity (ignoring the super-fluid component of the matter) rotates stationarly as a solid body with an angular velocity $\Omega$ that is independent of the spatial coordinates. The time scales for changes in the angular velocity which we will consider in our applications are well separated from the relaxation times at which hydrodynamical equilibrium is established, so that the assumption of a rigid rotator model is justified.

Now besides of central energy density $\varepsilon(0)$ of the star configuration the angular velocity of the rotation $\Omega$ is an additional parameter of the theory.

As a next step, going beyond the spherically symmetric case that corresponds to the first-order approximation, we solve Eq. (8), where the unknown function is $q(r,\theta)$ which is defined by Eq. (25) and scaled such that it is independent of the angular velocity. Using the static solutions Eqs. (30) and (33), and the representation of $q(r,\theta)$ by the series of the Legendre polynomials,

[TABLE]

one can see that this series is truncated and only the coefficient $q_{0}(r)$ is nontrivial, i.e., $q_{m}(r)=0$ for $m>0$ , see (Hartle, 1967; Chubarian et al., 2000). Therefore one can write down the equations for $\bar{\omega}(r)=\Omega(1+q_{0}(r)/{\Omega}_{0})$, which is more suitable for the solution of the resulting equation in first order

[TABLE]

which corresponds to Ref. (Hartle, 1967), where it was obtained using a different representation of the metric. In this equation we use the notation $j(r)\equiv e^{-(\nu^{(0)}(r)+\lambda^{(0)}(r))/2}$, where $j(r)=1$ for $r>R_{0}$, i.e., outside of stellar configuration.

Using this equation one can reduce the second order differential equation (47) to the first order one

[TABLE]

and solve (47) as a coupled set of first order differential equations, one for the moment of inertia (43) and the other (48) for the frame dragging frequency $\bar{\omega}(r)$.

This system of equations is valid inside and outside the matter distribution. In the center of the configuration $I_{0}(0)=0$ and $\bar{\omega}(0)=\bar{\omega}_{0}$. The finite value $\bar{\omega}_{0}$ has to be defined such that the dragging frequency $\bar{\omega}(r)$ smoothly joins the outer solution

[TABLE]

at $r=R_{0}$, and approaches $\Omega$ in the limit $r\to\infty$. In the external solution (49) the constant $I_{0}=I_{0}(R_{0})$ is the total moment of inertia of the slowly rotating star and $J_{0}=I_{0}\Omega$ is the corresponding angular momentum. In this order of approximation, $I_{0}$ is a function of the central energy density or the total baryon number only. This solution remains connected to the spherically distributed matter and therefore does not differ too much from our previous expression, which uses this solution to incorporate the relativistic corrections.

However, to find the explicit dependence of the moment of inertia on the angular velocity, one needs to take the second step and account for the deformation of the stellar configuration, which, in the framework of our scheme, is a second-order correction.

3.2 Rotational deformation and moment of inertia

To calculate these contributions and the internal structure of the rotating star which is deformed due to centrifugal force one needs to return to our perturbation description and take into the corrections in the diagonal elements of the metric and the energy-momentum tensor, as in the equations above with the parameter $\beta=(\Omega/\Omega_{0})^{2}$.

For a more detailed description of the solutions of the field equations in the $\sim O(\Omega^{2})$ approximation we refer to the works of Hartle (1967) as well as Chubarian et al. (2000). Since these equations have a complicated form, here we will discuss only the qualitative meaning of the physical quantities concerning the star’s deformation and its moment of inertia.

In $\Omega^{2}$-approximation the shape of the star is an ellipsoid, and each of the equal-pressure (isobar) surfaces in the star is an ellipsoid as well. All diagonal elements of the metric and energy-momentum tensors could be represented as a series expansion in Legendre polynomials, as we have already discussed in the previous section where is has been noted that the only non vanishing solutions obeying the continuity conditions on the surface with the external solution of fields are those with $l=0,2$.

The deformation of the isobaric surfaces can be parameterised by the deformation shifts $R(r,\theta)-r=\Delta(r,\theta)$ from the spherical shape. It describes the deviation from the spherical distribution as a function of radius $r$ for a fixed polar angle $\theta$ and is completely determined by

[TABLE]

since the expansion coefficients of the deformation $\Delta_{l}(r)$ are connected with the pressure corrections

[TABLE]

$l\in\{0,2\}$ is the polynomial index in the angular expansion in Legendre polynomials, analogous to Eq. (28).

The function $R(R_{0},\theta)$ is the radius where $p(R(R_{0},\theta))=0.$ To avoid confusion, we denote the spherical radius as $R_{0}$, which not anymore the actual radius, but rather $R(R_{0},\theta)$ which is the distance of the stellar surface from the center of the configuration at a polar angle $\theta$. In particular, we define the equatorial radius as $R_{e}=R(R_{0},\theta=\pi/2)$, the polar radius as $R_{p}=R(R_{0},\theta=0)$, and the eccentricity as $\epsilon=\sqrt{1-(R_{p}/R_{e})^{2}}$, all three quantities characterizing the deformed shape of the star.

Using this the same approach we write the correction to the moment of inertia as $\Delta I(r)=I(r)-I_{0}(r)$ and represent it as a sum of several different contributions,

[TABLE]

Since these contributions are obtained from the exact expression of angular momentum in the integral form the first three contributions can also be expressed by integrals of the form

[TABLE]

where integration is taken from the angular averaged modifications of the matter distribution, the shape of the configuration and the gravitational fields,

[TABLE]

respectively. All quantities appearing have been determined from the Eq. (5) in second order approximation. The contribution of the change of the rotational energy to the moment of inertia is given by

[TABLE]

and includes the frame dragging contribution.

In the next sections of this chapter we will discuss results for the moment of inertia along with stability conditions. We note that a consistent discussion of the stability of rotating stars requires one to take into account the contribution of the rotational energy to the mass energy, as well as the corresponding corrections to the moment of inertia.

4 Models for the EoS with a strong phase transition

In this contribution, we focus on EoS models which describe a strong phase transition in the sense that upon solving the TOV equations with them compact star sequences are obtained which exhibit a third family branch in the mass-radius or mass-central (energy) density diagram which is separated from the second family of neutron stars by a sequence of unstable configurations. The possibility of the very existence of a third family of compact stars as a consequence of a strong phase transition in dense nuclear matter, together with a sufficient stiffening of the high-density matter that can be expressed by a strong increase in the speed of sound (but not violating the causality bound) has been discusses by Gerlach as early as 1968 Gerlach (1968). Let us note here that the very existence of such a third family of compact stars is an effect of strong, general relativistic gravity! Namely, that the compactification which accompanies the strong phase transition of the star leads to a reduction of the gravitational mass of the hybrid star configuration from which it only recovers (and thus escapes gravitational collapse) when after the transition the hybrid star core consists of sufficiently stiff high-density matter.

Such EoS lead to mass-radius relationships for hybrid stars that have been classified (D)isconnected or (B)oth in Ref. Alford et al. (2013). The ”D” topology consists of a hadronic and a hybrid star branch, both of which being gravitationally disconnected from each other. In contrast to this, the ”B” topology consists of a branch of stable hadronic stars followed by stable hybrid stars, which are gravitationally disconnected from a second branch of stable hybrid stars. For the introduction of this classification scheme, the constant-speed-of-sound (CSS) EoS was used in Alford et al. (2013) for describing the high-density matter, see also Ref. Zdunik and Haensel (2013) for the justification of its validity in the case of color superconducting quark matter EoS. The first demonstration that high-mass twin stars and thus a corresponding high-mass third family sequence with $M_{\rm max}>2~{}M_{\odot}$ were possible, has been given in Alvarez-Castillo and Blaschke (2013) The intricacy of an equation of state describing a third family of stars (with twin stars at high or low masses as a consequence) consists in the fact that one needs, on the one hand, a sufficiently large jump in energy density $\Delta\varepsilon$ and a relatively low critical energy density $\varepsilon_{c}$ at the transition point to fulfil the Seidov criterion (1) for gravitational instability (i.e., a stiff nuclear matter EoS has to be followed by a soft high-density one), while on the other the high-density EoS needs to become sufficiently stiff directly after the phase transition, without violating the causality condition ($c_{s}^{2}<1$). With the CSS parametrization, these constraints could be fulfilled relatively straightforwardly by dialling $c_{s}^{2}=1$ and adjusting a sufficiently large value of $\Delta\varepsilon$ by hand.

The question arose whether a hybrid star EoS describing a third family of compact stars with a maximum mass above $2~{}M_{\odot}$ could also be obtained when applying the standard scheme of a two-phase approach based on a realistic nuclear matter EoS and a microscopically well-founded quark matter EoS, both joined, e.g., by a Maxwell construction. A positive answer was given already in 2013, when two examples of this kind were presented in Blaschke et al. (2013a), where the excluded-volume corrected nuclear EoS APR and DD2 were joined with a quark matter EoS based on the nonlocal NJL model approach Blaschke et al. (2007); Benic et al. (2014), augmented with a density dependent repulsive vector meanfield that was constructed by employing a thermodynamically consistent interpolation scheme introduced in Ref. Blaschke et al. (2013b). Such an interpolation scheme, based on the nonlocal, color superconducting NJL model of Blaschke et al. (2007), but extended to address also a density-dependent bag-pressure that facilitates a softening of the quark matter EoS in the vicinity of the deconfinement transition, has recently been developed in Alvarez-Castillo et al. (2019), guided by a relativistic density functional (RDF) approach to quark matter Kaltenborn et al. (2017) with an effective confinement mechanism according to the string-flip model Horowitz et al. (1985); Ropke et al. (1986). This approach has been applied very successfully to describe a whole class of hybrid EoS with a third family branch fulfilling the modern compact star constraints Ayriyan et al. (2018); Alvarez-Castillo et al. (2019). In the RDF approach to the string-flip model of quark matter an essential element is the ansatz for the nonlinear density functional resembling confinement and embodying the aspect of in-medium screening of the string-type confining interaction within an excluded volume mechanism. Another nonlinearity term in the density functional embodies the stiffening of quark matter at high densities in a similar way as it was obtained from 8-quark interactions in the NJL model Benic (2014) that were an essential part of the description of high-mass twins in Benic et al. (2015). Both nonlinearities, due to (de)confinement and high-density stiffening, are mimicked in a rather flexible way by the twofold interpolation scheme suggested in Alvarez-Castillo et al. (2019) which can be reinterpreted as a generalization of the nonlocal NJL model with chemical-potential-dependent parameters.

Having discussed the successful approaches to construct EoS with a strong phase transition which account for third family branches of compact stars and can be recognized observationally by the mass twin phenomenon, we would like to mention which ingredients are indispensable for obtaining this feature and which approaches have failed to obtain it. An excellent illustration of the various possibilities to join by interpolation hadronic and quark matter phases which themselves have different characteristics of stiffness, can be found in the recent review of Ref. Baym et al. (2018). However, despite being quite general, the case of the third family branch and mass-twin compact stars could not be captured! The reason can be found elucidated in Ref. Alvarez-Castillo et al. (2019), where it is demonstrated for a representative set of hadronic as well as quark matter EoS of varying degree of stiffness, that either a phase transition in the relevant domain of densities is entirely absent or results in a hybrid star branch that is directly connected to the hadronic brach which therefore does not form a third family of stars. The generated patterns are very similar to the results of Ref. Orsaria et al. (2014) which also uses the nonlocal NJL model for describing the quark phase of matter. The clue to obtaining third family sequences within microscopically motivated studies is a subtle softening followed by a stiffening of quark matter that can be realized employing the thermodynamically consistent interpolation procedure between different parametrizations of the same quark EoS (e.g., varying the vector meson coupling strength) before applying a Maxwell-, Gibbs-, or pasta phase transition construction.

We have described here the state-of-the-art modeling of EoS with a strong phase transition that are based on microscopic models of high-density (quark) matter. Besides these, there is a large number of simple EoS parametrizations that are also in use for discussing third family sequences fulfilling the constraint of a high maximum mass of the order of $2~{}M_{\odot}$. These are basically the classes of CSS based models and multi-polytrope approaches. Without attempting completeness, we would like to mention Refs. Alford et al. (2015); Alford and Han (2016); Christian et al. (2018); Alford and Sedrakian (2017); Paschalidis et al. (2018); Christian et al. (2019); Montana et al. (2019); Han and Steiner (2019) from the class of CSS models. A particularly interesting work is the extension by Alford and Sedrakian Alford and Sedrakian (2017), who demonstrated that also a fourth family of hybrid stars can be obtained and besides mass-twins there are also mass-triples possible. The multi-polytrope EoS have been a workhorse for numerical relativity studies of astrophysical scenarios since many years. The approach to constrain the multi-polytrope parameters from observations of masses and radii as introduced by Read et al. Read et al. (2009) has then been developed further with great resonance in the community by Hebeler et al. (2010, 2013). While in Read et al. (2009) one already finds a one parameter set describing high-mass twin stars that have not yet become a matter of interest in the community, the third family branch in the multi-polytrope aproaches has been mainly overlooked (see, e.g., Raithel et al. (2016); Miller et al. (2019)), but was digged up in Ref. Alvarez-Castillo and Blaschke (2017) where it was found that one should have at least a four-polytrope ansatz and suitably chosen densities for the matching of the polytrope pieces of the EoS, see also Annala et al. (2018); Paschalidis et al. (2018); Hanauske et al. (2018).

An important issue when discussing strong first-order phase transitions is the appearance of structures of finite size, like bubbles and droplets in the boiling/condensation transitions of the water-vapour transformations. In general, different shapes in the new phase are possible like spherical, cylindrical and planar structures, which have been dubbed ”pasta structures”. Their size and thermodynamical favorability depends on the surface tension between the subphases and the effects of the Coulomb interaction, including screening. The resulting pressure in the mixed phase is then no longer constant as in the Maxwell construction case, but also not as dramatically changing when the surface tension is neglected Glendenning (1992). For details concerning the quark-hadron transition pasta phases under neutron star constraints see, e.g., Refs. Na et al. (2012); Yasutake et al. (2014); M. Spinella et al. (2016).

It is interesting to note that a simple one-parametric parabolic approximation of the pressure versus chemical potential dependence can give a satisfactory agreement with a full pasta phase calculation and that the single parameter can be directly related to the surface tension Maslov et al. (2018). It could be demonstrated that the third family feature of an EoS with strong phase transition is rather robust against pasta phase effects, see Ayriyan et al. (2018).

In the present work, we will use the multi-polytrope approach to the EoS describing a third family of compact stars and also discuss the effect of mimicking the pasta structures in the mixed phase by a polynomial interpolation. The results shall not depend qualitatively on these simplifying assumptions but be of rather general nature and also applicable to more realistic types of EoS as discussed above.

4.1 Multi-polytrope approach to the EoS

In this section we present an EoS model that features a strong first-order phase transition from hadron to quark matter. They are labeled “ACB” following (Paschalidis et al., 2018) and consist of a piecewise polytropic representation (Read et al., 2009; Hebeler et al., 2013; Raithel et al., 2016) of the EoS at supersaturation densities ($n_{1}<n<n_{5}\gg n_{0}$):

[TABLE]

where $\Gamma_{i}$ is the polytropic index in each of the density regions labeled by $i=1\dots 4$. We fix $\Gamma_{1}$ such that a stiff nucleonic EoS provided in (Hebeler et al., 2013) can be described. The second polytrope shall correspond to a first-order phase transition therefore in this region the pressure must be constant, given by $P_{c}=\kappa_{2}$ ($\Gamma_{2}=0$). The remaining polytropes, in regions 3 and 4, that lie above the phase transition shall correspond to high-density matter, like stiff quark matter.

In order to compute the remaining thermodynamic variables of the EoS, we utilize the formulae given in the Appendix of Ref. (Zdunik et al., 2006),

[TABLE]

where we fix the integration constant $C$ by the condition that $\varepsilon(n\to 0)=m_{0}\,n$. Now we can invert the above expressions to obtain

[TABLE]

so that the chemical potential dependent pressure for the polytrope EoS (58) can be written as

[TABLE]

The above form of the pressure (63) is suitable to perform a Maxwell construction of a first-order phase transition between the hadron and quark phases. The model parameters for the mass twin cases that we consider in this work are given in table 4.1.

The ACB4 EoS features a first order phase transition at a rather high nucleon number density value, $n_{2}=0.3174$ fm*-3* that produces an instability in a $2.0M_{\odot}$ neutron star, providing an example of the high-mass twins phenomenon. On the contrary, the ACB5 EoS presents a phase transition that occurs at the lower value of $n_{2}=0.284~{}\mathrm{fm}^{-3}$. In that case the instability occurs for stars with $1.4M_{\odot}$, providing a scenario for low-mass twin stars. This low density value for the phase transition of about two times saturation density, is particularly feasible in neutron star matter, where the effect of the isospin asymmetry manifests in the so called asymmetry energy which stiffens the EoS with respect to the symmetric case, equal number of protons and neutrons in hadronic matter. Figure 1 shows the resulting mass-relation curves for these two EoS featuring mass twins together with measurements and constraints regions.

4.2 EoS including mixed phase effects (pasta phases)

In this section we introduce a mixed phase approach to mimic pasta structures in regions of both the hadronic and quark EoS around the Maxwell critical point ($\mu_{c},P_{c}$). The method used is the replacement interpolation method (RIM) (Abgaryan et al., 2018) that consists of replacing the EoS in the aforementioned domain by a polynomial function:

[TABLE]

with $\Delta_{P}$ as free parameter that adds pressure to the mixed phase at $\mu_{c}$. Generally, all parametrizations of the type shown in Eq. (64) for the mixed phase pressure are even order ($N=2k$, k=1, 2, …) polynomials which we refer to as $G_{k}$. In order to smoothly match the EoS at $\mu_{H}$ and $\mu_{Q}$ up to the $k$-th derivative of the pressure the following conditions shall be fulfilled

[TABLE]

with $\alpha_{q}$, $\mu_{H}$ and $\mu_{Q}$ being determined by the above system of equations.

For the sake of simplicity we employ the parabolic model $G_{1}$ of the RIA as introduced in (Ayriyan and Grigorian, 2018; Ayriyan et al., 2018):

[TABLE]

where the parameters $\alpha_{1}$, $\alpha_{2}$, $\mu_{H}$ and $\mu_{Q}$ are to be determined as described above, from the continuity conditions at the Maxwell construction critical point:

[TABLE]

Figure 3 shows a schematic representation of the RIM method based on the Maxwell construction between the hadronic and quark EoS.

In addition, figure 4 shows the mixed phase equations of state for both low and high mass twins. The effect of the mimicked geometrical structures is quantified by the $\Delta_{P}$ parameter. It is evident that the order of the $G$ function will result in whether or not there are discontinuities for the derivatives of the $P_{M}$ function. For instance, the square of the speed of sound, $c_{s}^{2}$, is proportional to the second derivative of $G$ with respect to $\mu$, see figure 5.

The result is that $G_{1}$ presents a clear discontinuity in the speed of sound at $\varepsilon_{c}$ and $\varepsilon_{c}+\Delta\varepsilon$, whereas in between i.e., in the latent heat region, the speed of sound slightly increases. On the contrary, the construction $G_{2}$ allows for a continuous speed of sound, however it is not smoothly connected at $\varepsilon_{c}$ and $\varepsilon_{c}+\Delta\varepsilon$. Only $G_{3}$ is capable of joining smoothly the speed of sound between the hadron and quark EoS at the critical points.

5 Results

5.1 TOV solutions for mixed phase models

In figure 7 we show the results of the mass-radius diagram as a solution of the TOV equations for the equations of state ACB4 (left panel) which exhibits high mass twin stars and ACB5 (right panel) which describes low mass twins, depending on the value of the mixed phase parameter $\Delta_{P}$. In the insets we give a magnified view on the region of the maximum mass of the hadronic branch of the sequence, where the dotted lines indicate the unstable solutions that qualify the corresponding EoS as one with a third family. We can read off to the accuracy of the given 1% steps what the critical value for $\Delta_{P}$ is when the disconnected second and third families would merge to a connected hybrid star sequence.

While for the case of ACB4 the variation of $\Delta_{P}$ does not affect the mass-radius diagram in the mass region of the compact star merger GW170817, the corresponding variation for ACB5 leads to strong effects in that mass region. We therefore consider the tidal deformability in both cases in the next subsection.

5.2 Tidal deformability predictions

Together with the solution of the TOV equations, one can solve for the dimensionless tidal deformability $\Lambda(M)$ for the given EoS. After that, one can construct the corresponding lines in the $\Lambda_{1}-\Lambda_{2}$ diagram of the binary compact star merger GW170817 for which the individual masses $M_{1}$ and $M_{2}$ of the two compact stars fulfill the constraint derived from the detected gravitational wave signal of the inspiral phase of the merger Abbott et al. (2017). These lines can be overlaid to the constraint derived from the LVC observation, as shown in Fig. 7. As to be expected from Fig. 6, only in the case of ACB5 we can note an effect of the mixed phase construction while the results for ACB4 are inert against changes of the mixed phase parameter, because it influences the mass-radius diagram in a region of masses that is inaccessible to the gravitational wave signal of the inspiral phase and the effects of tidal deformation. Moreover, we notice that ACB4 is too stiff an EoS to fulfil the compactness constraint from GW170817. The EoS ACB5, however, with the early onset of the phase transition, becomes a soft EoS due to the mixed phase effects and for the largest values of the mixed phase parameter $\Delta_{P}$ is similar to a soft hadronic EoS despite the fact that the compact stars consist of extended regions of quark matter in pure or mixed phases.

In the following section, we will consider the effects of fast rotation on the sequences of hybrid star solutions and shall obtain a qualitative difference in the characteristics of pure phase (hadronic) and hybrid stars concerning their maximum masses which are relevant for the discussion of the phenomenology of binary compact star mergers and their implications for the state of superdense matter.

5.3 Rotating compact star solutions

In this subsection we present the numerical solutions for rotating hybrid star sequences in full GR equations for axial symmetry as obtained with the RNS code described in subsection 2.2 and in the perturbative expansion up to order $\Omega^{2}$ (slow rotation approximation) that was explained in subsection 2.3. We relate these solutions to those of the TOV equations for the static case of spherical symmetry discussed in the previous subsection 2.4.

In Fig. 8 we show the gravitational mass as a function of the central energy density for both multipolytrope EoS, ACB4 (with a deconfinement phase transition at high-mass, upper panels) and ACB5 (with a transition at the typical compact star mass of $1.4~{}M_{\odot}$, lower panels) for nonrotating (blue solid lines) and maximally rotating ($\Omega=\Omega_{K}$) stars in full GR (red dashed lines) and slow rotation approximation (orange dash-dotted line). In the left panels for the Maxwell construction case the jump in the central energy density by about a factor two at the onset of the transition is clearly seen and such an amount of latent heat is sufficient, according to the Seidov criterion Eq. (1), to trigger a gravitational instability which occurs in the region of densities where

[TABLE]

At the onset of this instability the gravitational mass of the star has reached the maximum attainable on the second family of compact stars, denoted by $M_{\max,2}$. For the Maxwell construction case, this mass is degenerate with that for the onset of the phase transition, $M_{\rm onset}$. Due to the absence of a pressure gradient in the interval of densities corresponding to the mixed phase, this phase is not realized in compact stars in this case. The EoS with mass twin compact star sequences are characterized by the fact that the instability criterion (74) is fulfilled in a finite interval of densities which is then followed by another stable, rising branch of sequences, the so-called third family of compact stars. This behavior defines two more characteristic masses: $M_{\rm min,3}$ at the lower and $M_{\rm max,3}$ at the upper turning point, see the left panels of Fig. 8.

For the mixed phase constructions depicted in the right panels of Fig. 8, the pressures at the onset and the end of the mixed phase are not identical and thus, due to the corresponding pressure gradient, a mixed phase can be realized in the star and the degeneracy between $M_{\rm onset}$ and $M_{\max,2}$ is lifted. On the other hand, by our choice of the value of the mixed phase parameter $\Delta_{P}$ close to the limiting value for which the second and the third families of compact stars would get connnected, the instability vanishes and thus $M_{\rm min,3}$ joins $M_{\max,2}$, so that for a slightly larger value of $\Delta_{P}$ both these masses can no longer be identified since only the second family survives which for $M>M_{\rm onset}$ consists of hybrid stars.

These four characteristic masses for a given mass twin compact star EoS are given in Tab. 5.3 for the static case obtained by solving the TOV equations. In the last column the absolute maximum of the mass-radius curve for the given EoS is listed. In tables 5.3 and 5.3 these five characteristic masses are listed for the sequences of stars rotating at the Kepler frequency $\Omega_{K}$ which are obtained from solutions of the axisymmetric Einstein equations in the $\Omega^{2}$ approximation and in full General Relativity, respectively.

Inspecting the rotating star sequences for the close-to-critical mixed phase parameter $\Delta_{P}$, we observe that due to the rotation the star branches with mixed phase and pure quark matter core can get disconnected so that the phenomenon of a third family reappears. Vice-versa, upon spin-down from a supramassive star configuration at maximal rotation frequency (created, e.g., in a binary neutron star merger) which is stable on the hadronic or mixed phase branch may end up either as a black hole or on the hybrid star branch for such mixed-phase EoS. The scenario of a delayed collapse to a black hole is of special importance for interpreting GW170817 and will therefore be discussed below in further detail. In this context appears the question whether between the above-introduced characteristic masses at maximal and at zero rotation frequency hold EoS-independent, so-called universal relations that have been investigated for the maximum mass of hadronic EoS Bozzola et al. (2018); Rezzolla et al. (2018) and recently also for hybrid EoS including the Maxwell construction cases of ACB4 and ACB5 Bozzola et al. (2019). We shall come back to this issue below.

Here we like to remark that the maximum masses on the second and third family branches correspond to stars with very different central (energy) densities. This may be the clue to understanding the fact that the increase in mass for stars on the more compact third family branch is smaller than for stars on the secnd family one because of their smaller radii and thus smaller moment of inertia (57) and rotational energy, see the gravitational mass vs. equatorial radius in Fig. 9. For a more quantitative discussion of the effects of rotation on the masses of the sequences, we extract from the tables 5.3-5.3 the ratios of the characteristic masses on the rotating sequences to those on the static ones in table 5.3 and table 5.3 for the $\Omega^{2}$ approximation and the full GR solution, respectively. For completeness, we give the ratio of the characteristic masses between the two rotation solutions in table 5.3.

6 Implications for the phenomenology of compact stars

We investigate the consequences of a strong phase transition in the EoS for dense compact star matter for sequences of configurations in the mass-radius as well as mass-central (energy) density plane, with and without rotation. While for isolated pulsars even the highest known spin frequencies are well below the Kepler frequency so that no strong modification of the TOV solution occurs, in the era of multi-messenger astronomy, with compact star mergers being accessible to detection by their GW signal from the inspiral phase and soon also from the postmerger state, the ”übermassive” Espino and Paschalidis (2019) as well as supramassive star solutions play a role in the interpretation of the observations. While the former are solutions for differentially rotating configurations at the mass shedding limit, the latter are uniformly rotating objects with a frequency close to the Kepler one Shibata et al. (2019). Recently, for this case EoS independent, so-called ”universal” relationships have been derived which relate the maximum mass of the supramassive star sequence to the static one from the solution of the TOV equation. Such relationships, once confirmed, are particularly useful in order to make general predictions or draw conclusions from merger phenomenology for constraints limiting the EoS properties. In this context we would like to mention the upper limit on the maximum mass of neutron stars that has been extracted from the GW signal and phenomenology of GW170817 Shibata et al. (2017); Margalit and Metzger (2017); Rezzolla et al. (2018).

It has been known since long Haensel et al. (2007) that there is a relationship between the maximum mass of uniformly rotating cold neutron stars at the maximum frequency and the maximum mass of the TOV equation solution for static stars as

[TABLE]

where recently in Ref. Breu and Rezzolla (2016) the universality of this relationship was confirmed for a very large set of hadronic EoS (without a deconfinement phase transition) with the coefficient $\alpha=1.20$. The hypothesis that the relation (75) can be extended to include hybrid stars with a strong phase transition and even with third family sequences has recently been investigated in Ref. Bozzola et al. (2019) and following the argumentation of Ruiz et al. (2018); Most et al. (2018); Rezzolla et al. (2018) it leads to a limitation for the maximum mass of static neutron stars as

[TABLE]

The value of $M_{\rm max}=2.591~{}M_{\odot}$ was extracted for the core mass of the compact star merger GW170817 in Ref. Rezzolla et al. (2018). We confirm the finding of Bozzola et al. (2019) that the coefficient spans a range of values for which we find $1.16<\alpha<1.33$, see the table 5.3. The lower limit in (76) comes from the new high-mass pulsar PSR J0740+6620 for which a mass $2.17^{+0.11}_{-0.10}~{}M_{\odot}$ has been determined by Cromartie et al. (2019) by measuring the Shapiro delay, and the upper limit in our case is $2.23~{}M_{\odot}$ for the lowest value of $\alpha=1.16$. In order to not come in conflict with the pulsar mass measurement of (Cromartie et al., 2019), there is an upper limit for the admissible value of $\alpha=1.25$, corresponding to the lower limit at the $1\sigma$ level of the PSR J0740+6620 mass, $M_{\rm max}^{\rm TOV}=2.07~{}M_{\odot}$. Taking the central value, $M_{\rm max}^{\rm TOV}=2.17~{}M_{\odot}$, would correspond to $\alpha=1.20$, see Fig. 10. This figure illustrates one of the main findings of this contribution, the dependence of the coefficient $\alpha$ in Eq. (75) on the central (energy) density of the stellar configuration that can be fitted by a linear regression to the values we determined at the positions corresponding to the characteristic masses to be

[TABLE]

where $a_{ACB4-4}=1.38\pm 0.07$, $b_{ACB4-4}=0.12\pm 0.01$ fm3/GeV and $a_{ACB5-2}=1.37\pm 0.07$, $b_{ACB5-2}=0.16\pm 0.02$ fm3/GeV.

We may conclude that only those states of matter are allowed for the inner core of a compact star at maximum mass of $2.07~{}M_{\odot}$ ($2.17~{}M_{\odot}$) which belong to a high-density region with $\varepsilon\geq 0.78$ GeV/fm3 ($\varepsilon\geq 1.12$ GeV/fm3).

Anyway, the main effect of the strong phase transition is a higher compactness of the high mass stars than in the purely hadronic case which moreover goes along with a smaller mass increase due to maximal rotation than in the purely hadronic case which entails the increase of the upper limit for the maximum mass relative to the purely hadronic case discussed in Shibata et al. (2017); Margalit and Metzger (2017); Rezzolla et al. (2018). Thus the high-density phase transition removes a certain tension from the discussion of the upper limit for the maximum mass and could be used as an argument in favour of the suggestion that a strong phase transition actually takes place in compact stars!

7 Summary and Conclusions

Stimulated by the unprecedented progress in observational astronomy, compact stars have become superb astrophysical laboratories for a broad range of physical studies. This is particularly the case for neutron stars, since their observables carry information about the fundamental building blocks of matter and even of the fabric of space itself. Against this background we did present in this book chapter a systematic investigation of the properties of compact stellar mass twins (i.e., the so-called third family of compact stars), which, according to theory, may exist in the mass-radius region between neutron stars and stellar-mass black holes. Particular emphasis is given to modeling the rotational properties of compact mass twins for multi-polytrope models for the equation of state of ultra-dense stellar matter that fulfill the constraint established for the maximum mass of a neutron star. The main results of our investigation can be summarized as follows:

1. The existence of mass twins invariably signals the existence of a strong phase transition in ultra-dense matter. The extreme softening of the equation of state caused by the strong phase transition increases the gravitational field so much that the star becomes gravitationally unstable over a certain range of densities, where just a certain fraction of the star’s core is in the new high-density phase. Eventually the stars becomes stable again when about half of the matter in its core is in the new phase of matter. The new stars have gravitational masses that are less than the maximal mass of the hadronic stellar sequence.

The observation of mass twins would indicate a strong phase transition so that from the existence of a strong first-order phase transition in one corner of the QCD phase diagram and a crossover behavior in another, one could conclude that at least one critical endpoint (CEP) must exist. This would be very reassuring for large scale experimental heavy-ion collision programs set up for the search for the CEP.

2. The mixed phase construction which mimics the pasta phase is in accordance with a full pasta calculation. The result is a ”smearing” of the phase transition over a certain pressure region, which is similar to the Gibbs construction in matter with more than one conserved charge and where charge conservation need not be fulfilled locally but rather globally Glendenning (1992). This construction makes not only the approach to the strong phase transition more realistic, but has also great advantages for the numerical realization of phenomenological scenarios of the phase transition in rapidly rotating stars studied in numerical relativity.

3. Lastly, we address the conjecture of an upper limit on the maximum masses of nonrotating compact stars from the phenomenology of GW170817 and its associated kilonova event. The conjecture is based on a quasi-universal relation between the maximum masses of uniformly rotating stars at the maximum frequency and that of static TOV solutions for the same EoS, which was demonstrated to hold for neutron star EoS without a phase transition Breu and Rezzolla (2016). The stellar mass at which the high-density phase transition (such as deconfinement of quarks) sets in is currently unknown. But if the transition would occur below the maximum mass of the TOV solution, the quasi-universal relation will have to be revisited and the conclusions for the upper limit on the maximum mass be revised. We have found a quantitative criterion for the minimal central energy density in the maximum-mass configuration of a compact star that would correspond to the core of GW170817 after dynamical mass ejection Rezzolla et al. (2018). Thus the EoS at high densities must be effectively soft, either in the form of a relatively soft hadronic EoS or as a hybrid EoS with a phase transition, since a too stiff a hadronic EoS would lead to heavy hadronic stars too dilute to fulfill the constraint derived by us.

With these prospects for strong gravity and strong phase transitions in compact stars, we are looking forward to the next series of exciting discoveries in the just opened era of multimessenger astronomy.

Acknowledgements

We would like to thank Andreas Bauswein for helpful discussions and Vasilis Paschalidis for pointing out caveats in the discussion of the maximum mass constraint from GW170817. We are grateful to Cesar Zen Vasconcellos for the invitation to contribute to this book and for his steady encouragement to complete the writeup despite interfering obligations. The work of D.B., A.A., H.G. has been supported by the Russian Science Foundation under grant 17-12-01427; F.W. acknowledges support by the U.S. National Science Foundation under Grant PHY-1714068. The authors are grateful to the European COST Actions CA15213 ”THOR” and CA16214 ”PHAROS” for supporting their networking activities.