# Book Title

# Topics on Strong Gravity

## Chapter 1 Astrophysical aspects of general relativistic mass twin stars

David Blaschke,
David Edwin Alvarez-Castillo,

Alexander Ayriyan,
Hovik Grigorian,
Noshad Khosravi Lagarni,

Fridolin Weber

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

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

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

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

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

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

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

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

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.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 , 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 , baryonic density , and pressure . A fourth quantity of great interest, the chemical potential, can then be obtained as .

In addition, the most basic conditions that the system must fulfill include global charge neutrality and -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 , where and 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 with a typical radius of around 12 km, the general relativistic correction factor amounts to [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 , central pressure , and the jump in energy associated with the phase transition obey the relation

(1.0) |

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 , thereby populating a third branch with stable mass twins.

### 1.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,

(1.0) |

The curvature of the space-time is satisfying the Einstein field equations where . Here is the Ricci curvature tensor and the scalar curvature. On the right hand side of the Einstein equation we have the energy -momentum tensor of the stellar matter and is the gravitational constant (.

The metric tensor 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 (-time coordinate, - 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.

#### 1.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

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 and altitude angle 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

(1.0) |

where is the -velocity of matter, the pressure and the energy density.

Once the energy-momentum tensor (1.0) is fixed by the choice of the equation of state for stellar matter, the unknown metric functions ,, , 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,

(1.0) | |||||

(1.0) | |||||

(1.0) |

and one for the determination of the non-diagonal element

(1.0) |

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

(1.0) |

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

#### 1.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)]:

(1.0) | |||||

which just like Eq. (1.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 , where 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 , and in term of integrals, whereas the 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 , where with g cm.

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

(1.0) | |||||

(1.0) | |||||

(1.0) | |||||

where and with as the coordinate radius of the equator. Moreover, all the resulting quantities can be written in terms of the auxiliary variables and s, defined as and , respectively. Consequently, the four metric functions acquire a dependence on the above variables: , , , . The remaining quantities are:

(1.0) |

(1.0) |

(1.0) |

(1.0) |

where the subscripts and 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/

#### 1.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 . 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 . In case of rotating stars the deformability is connected with the induced centrifugal force and the expansion parameter is , where with the mass density 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 for the star with total mass and equatorial radius, the expansion gives sufficiently correct solutions already at . So the expansion parameter is naturally limited to values by this condition of mechanical stability of the rigid rotation, because always . 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. 1.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 (-) 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 can be expressed as

(1.0) |

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

(1.0) |

and for the frame dragging frequency the odd orders

(1.0) |

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,

(1.0) | |||||

(1.0) |

where and denote the zero-order coefficients which correspond to non-deformed spherically symmetric stars.

Because of rotational symmetry the diagonal elements (no summation over ) of the metric coefficients can be written as ( and are even values only)

(1.0) |

The same is true also for the non diagonal elements, but the angular dependence is different ( and values are now odd only). The case for the non diagonal term will be investigated in section 1.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.

#### 1.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.

(1.0) | |||||

where 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))

(1.0) |

where and denote the equation of state (EoS) describing the stellar matter. The quantity , defined as

(1.0) |

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

(1.0) |

where denotes the radius of the star defined by . 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. (1.0) 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 ). 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 . Stars on the stellar branch where are unstable agains radial oscillations and will therefore not exist stably in the universe.

Each stellar model has unique solutions for , and in terms of which the internal gravitational field (the metric coefficients) is defined as

(1.0) | |||||

(1.0) |

The internal field solutions are smoothly connected to the external field solutions at the stellar surface, 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

(1.0) |

where is the baryon number density and 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. (1.0) and (1.0) can be found from Eq. (1.0) and Eq. (1.0) in zeroth order of the -expansion.

### 1.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 for the
diagonal elements^{1}^{1}1Notation corresponds to the work of
Chubarian et al. (2000).. The the first correction corresponding
to small deformations (small values of ) one can consider terms
linear in or equivalently the perturbation approximation
to the spherically symmetric star. We introduce some new notation and
work under the assumption that like in
the expansion of the external solution, since
.

The non diagonal term could be taken to be zero, because we consider only the static case (the parameter defining the static deformation does not change the sign under time reversal ) and , so that we have

(1.0) | |||||

In this approximation, without loss of generality, one can set the values of , because we neglect the contribution of the deformation energy to the gravitational mass. The equations show that the prime symbol denoting the derivative of those quantities with respect to . The functions and obey the differential equations

(1.0) | |||||

(1.0) | |||||

Here, 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 and as , The is a constant that quantifies the deformation of the star which can be taken arbitrary. This constant corresponds to the choice of 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 and the following combination

(1.0) |

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

where in the expression is the compactness of the star ( is the ratio of gravitational radius to spherical radius).

The dimensionless tidal deformability parameter is defined as , a quantity defined for small tidal deformabilities. Here is the TD of the star with a gravitational mass , just as defined above. In addition, the love number is related to TD and defined as

(1.0) |

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

#### 1.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)

(1.0) |

where denotes the total, conserved angular moment. The quantity is given by

(1.0) |

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

(1.0) |

where is the invariant volume and . 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 accumulated in the sphere with radius , we obtain from Eq. (1.0)

Here the difference of the frame dragging frequency and the angular velocity . In general relativity, due to the Lense-Thirring law, rotational effects are described by

(1.0) |

This expression is approximated from the exact expression of the energy momentum tensor coefficient and the metric tensor in the axial symmetric case.

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

(1.0) |

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 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 of the star configuration the angular velocity of the rotation 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. (1.0), where the unknown function is which is defined by Eq. (1.0) and scaled such that it is independent of the angular velocity. Using the static solutions Eqs. (1.0) and (1.0), and the representation of by the series of the Legendre polynomials,

(1.0) |

one can see that this series is truncated and only the coefficient is nontrivial, i.e., for , see (Hartle, 1967; Chubarian et al., 2000). Therefore one can write down the equations for , which is more suitable for the solution of the resulting equation in first order

(1.0) |

which corresponds to Ref. (Hartle, 1967), where it was obtained using a different representation of the metric. In this equation we use the notation , where for , i.e., outside of stellar configuration.

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

(1.0) |

and solve (1.0) as a coupled set of first order differential equations, one for the moment of inertia (1.0) and the other (1.0) for the frame dragging frequency .

This system of equations is valid inside and outside the matter distribution. In the center of the configuration and . The finite value has to be defined such that the dragging frequency smoothly joins the outer solution

(1.0) |

at , and approaches in the limit . In the external solution (1.0) the constant is the total moment of inertia of the slowly rotating star and is the corresponding angular momentum. In this order of approximation, 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.

#### 1.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 .

For a more detailed description of the solutions of the field equations in the 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 -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 .

The deformation of the isobaric surfaces can be parameterised by the deformation shifts from the spherical shape. It describes the deviation from the spherical distribution as a function of radius for a fixed polar angle and is completely determined by

(1.0) |

since the expansion coefficients of the deformation are connected with the pressure corrections

(1.0) |

is the polynomial index in the angular expansion in Legendre polynomials, analogous to Eq. (1.0). The function is the radius where To avoid confusion, we denote the spherical radius as , which not anymore the actual radius, but rather which is the distance of the stellar surface from the center of the configuration at a polar angle . In particular, we define the equatorial radius as , the polar radius as , and the eccentricity as , 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 and represent it as a sum of several different contributions,

(1.0) |

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

(1.0) |

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

(1.0) | |||||

(1.0) | |||||

(1.0) |

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

(1.0) |

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.

### 1.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 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 and a relatively low critical energy density at the transition point to fulfil the Seidov criterion (1.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 (). With the CSS parametrization, these constraints could be fulfilled relatively straightforwardly by dialling and adjusting a sufficiently large value of by hand.

The question arose whether a hybrid star EoS describing a third family of compact stars with a maximum mass above 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 . 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.

#### 1.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 ():

(1.0) |

where is the polytropic index in each of the density regions labeled by . We fix 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 (). 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),

(1.0) | |||||

(1.0) | |||||

(1.0) |

where we fix the integration constant by the condition that . Now we can invert the above expressions to obtain

(1.0) |

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

(1.0) |

The above form of the pressure (1.0) 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 1.3.

The ACB4 EoS features a first order phase transition at a rather high nucleon number density value, fm that produces an instability in a 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 . In that case the instability occurs for stars with , 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.1 shows the resulting mass-relation curves for these two EoS featuring mass twins together with measurements and constraints regions.

#### 1.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 (). 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:

(1.0) |

with as free parameter that adds pressure to the mixed phase at . Generally, all parametrizations of the type shown in Eq. (1.0) for the mixed phase pressure are even order (, k=1, 2, …) polynomials which we refer to as . In order to smoothly match the EoS at and up to the -th derivative of the pressure the following conditions shall be fulfilled

(1.0) | |||||

(1.0) | |||||