Third family of compact stars within a nonlocal chiral quark model equation of state

A class of hybrid compact star equations of state is investigated that joins by a Maxwell construction a low-density phase of hadronic matter, modeled by a relativistic meanfield approach with excluded nucleon volume, with a high-density phase of color superconducting two-flavor quark matter, described within a nonlocal covariant chiral quark model. We find the conditions on the vector meson coupling in the quark model under which a stable branch of hybrid compact stars occurs in the cases with and without diquark condensation. We show that these hybrid stars do not form a third family disconnected from the second family of ordinary neutron stars unless additional (de)confining effects are introduced with a density-dependent bag pressure. A suitably chosen density dependence of the vector meson coupling assures that at the same time the $2~$M$_\odot$ maximum mass constraint is fulfilled on the hybrid star branch. A twofold interpolation method is realized which implements both, the density dependence of a confining bag pressure at the onset of the hadron-to-quark matter transition as well as the stiffening of quark matter at higher densities by a density-dependent vector meson coupling. For three parametrizations of this class of hybrid equation of state the properties of corresponding compact star sequences are presented, including mass twins of neutron and hybrid stars at 2.00, 1.39 and 1.20 $M_\odot$, respectively. The sensitivity of the hybrid equation of state and the corresponding compact star sequences to variations of the interpolation parameters at the 10% level is investigated and it is found that the feature of third family solutions for compact stars is robust against such a variation. This advanced description of hybrid star matter allows to interpret GW170817 as a merger not only of two neutron stars but also of a neutron star with a hybrid star or of two hybrid stars.


I. INTRODUCTION
Recently, in the context of the observation of pulsars with high masses of about 2 M by Demorest et al. [1][2][3] and Antoniadis et al. [4], the question of the possible existence of a third family of compact stars was revived, because of its relation to strong phase transitions in dense matter [5].The question was asked: could it be that stars at this high mass appear as mass twins [6] or almost mass twins like in the case above, where stars have about the same mass but may have very different radii, pointing to a very different structure of their interiors.In such a case, an explanation would be that the larger star would be an ordinary neutron star while the more compact one would exhibit a core composed of high-density matter, e.g., quark matter, with a large density jump at the border between inner core and outer core.The condition on the magnitude of the jump in energy density for the instability of neutron star configurations to occur is given by the Seidov relation [7].In addition, the high-density matter has to be stiff enough to allow for a stable hybrid star branch, the so-called "third family" of compact stars.The observational verification of the existence of mass twins would therefore imply very important lessons for the properties of compact star matter: the existence of a strong phase transition to a new form of high-density matter and the necessity that this matter obeys a sufficiently stiff equation of state (EoS).It is not only the case of high-mass twins which are interesting in this context.Also the possibility of mass twins in the range of typical neutron star masses around 1.4 M is very interesting!In such a case, the onset of the phase transition shall occur sufficiently early, at densities reached in the center of a typical-mass neutron star.It shall be soft enough at those densities to allow for a large enough density jump and stiffen quickly at increasing density to prevent gravitational collapse on the third family branch at least until the mass constraint of about 2 M is reached without violating the causality constraint.Namely that the speed of sound should not exceed the speed of light.If these constraints could be fulfilled, the corresponding EoS could provide a viable scenario for the recently observed [8] binary neutron star merger event: at least one of the two stars could be a hybrid star from the third family branch which is sufficiently compact to make the binary system fulfill the condition on the tidal deformabilities derived from the observation of the inspiral in the LIGO gravitational wave detector, see [9,10] for a recent discussion of this scenario.
Up to now, the third family case was investigated with very schematic EoS for the high-density phase, like the bag model [6], the constant-speed-of-sound (CSS) model [9,[11][12][13][14][15][16][17], the multi-polytrope model [9,18,19], but also dynamical models for interacting quark matter like the Nambu-Jona-Lasinio (NJL) model with higher order quark interactions in the Dirac vector channel [20,21], or a relativistic density-functional model [10,22].The question arises whether the third family phenomenon could be obtained also for EoS derived from dynamical quark models which are closer to QCD in the sense that they take into account the running of the dynamical quark masses with 4-momentum and embody a (dynamical) confinement mechanism.Such models are provided by QCD Dyson-Schwinger equations in their generalization to finite temperature and chemical potential [23] and have recently been applied to study hybrid compact stars, see, e.g., [24,25] and references therein.The dynamical model for the gluon propagator mediating the nonperturbative quark interactions is still quite schematic and more realistic ones including a matching with perturbative QCD behaviour at large momentum transfer make applications at high densities for compact star matter rather involved.Therefore, a separable ansatz for the nonlocal dynamical interactions has been suggested [26,27] and was used for the description of hadron properties in the QCD vacuum, see [28][29][30][31] for early works.This ansatz allows to develop covariant nonlocal chiral quark models for low-energy QCD at finite temperatures and densities that share the running of the quark mass function [26,32] and the wave function renormalization of the quark propagator [33][34][35][36][37] with full QCD as probed in lattice QCD simulations [38].This approach allowed a description of the QCD phase diagram [39][40][41][42] and has been also applied to the description of hybrid compact stars with quark matter cores [43] For a recent discussion of the issues occuring in the description of phases of dense matter in compact stars under modern constraints see, e.g., Ref. [44].Models with the capability to address the occurrence of a third family of compact stars shall fulfill these conditions [20] • stiff hadronic EoS (to obtain a phase transition in the observed mass range of compact stars) • stiff high-density EoS (to generate a stable hybrid star branch which should reach 2 M ), • sufficient density jump at the phase transition (to produce the instability as necessary condition for the existence of a third family branch).
Difficulties arise when one attempts to describe systematically a possible medium dependence of the interaction model that would go beyond the covariant formfactor ansatz and address Lorentz symmetry breaking effects [45].Of particular importance for the description of cold degenerate QCD matter as in compact stars, a medium dependence of the vector meson interaction is a crucial effect as it directly related to the density dependence of the stiffness of QCD matter between the deconfinement transition and the perturbative QCD regime at asymptotically high densities.In order to capture such nonperturbative medium effects on the basis of the nonlocal, covariant formfactor model it has been suggested in Ref. [46] to apply an interpolation technique similar in spirit to the one introduced for a flexible description of the deconfinement transition region in cold degenerate matter by Masuda et al. [47,48] who followed earlier concepts of interpolation for the description of lattice QCD thermodynamics in the high-temperature region [49].A recent discussion of this interpolation technique can be found in [50] where also the necessity to capture softening effects from quark confinement has been pointed out.
In the present work we develop the interpolation technique for nonlocal covariant formfactor approaches further in order to fulfill the above requirements on models of compact star matter which would produce third family solutions for hybrid stars with quark matter cores that at the same time obey the constraints on their maximum mass to exceed the present observational bound of 2.01 ± 0.04 M [4] and on their compactness to be in accord with the bounds derived from the gravitational waves detected for the inspiral phase of the binary compact star merger GW170817 [8].To this end we suggest here a twofold interpolation approach that captures the density dependence of both, the confining effects (the vanishing of a bag-like negative pressure of the nonperturbative QCD medium) and the stiffness (the increase of the vector meson coupling strength).

II. THEORETICAL FORMALISM
In the present work we use a two-phase description in order to account for the transition from nuclear to quark matter.In the following two subsections we explain the theoretical approaches used to obtain the EoS for each of these phases.Throughout this work we use natural units with = c = k B = G = 1.

A. Nuclear matter equation of state
For the nuclear matter phase, we use the wellestablished relativistic density-functional model by Typel et al. [51] based on meson-exchange interactions and with the so-called "DD2" parameterization given in Ref. [52].It describes the properties of nuclear matter at saturation density and below very well, also in accordance with the chiral effective field theory approach [53]; see also Ref. [54].To improve the higher-density behavior, a generalized excluded volume effect is included according to Ref. [55].The DD2 p40 EoS features this correction by considering the available volume fraction Φ N for the motion of nucleons as density dependent in a Gaussian form and Φ N = 1 if n ≤ n 0 .Here, v = 16πr 3 N /3 is the van der Waals excluded volume for a nucleon with a hard-core radius r N and n 0 = 0.15 fm 3 is the saturation density of infinite, symmetric nuclear matter.The index "p40" with the DD2 parametrization denotes a positive excluded volume parameter of v = 4 fm 3 .This type of nuclear EoS has recently been extensively used in systematic studies of hybrid star models, see for instance [20,22,56].

B. Quark matter equation of state
For the description of the quark matter phase we consider a nonlocal chiral quark model of the Nambu-Jona-Lasinio type (nlNJL) which includes scalar quarkantiquark interaction, anti-triplet scalar diquark interactions and vector quark-antiquark interactions that was presented in detail in Ref. [43].Let us start writing the corresponding effective Euclidean action, that in the case of two light flavors reads [43] Here m c is the current quark mass, that we consider to be equal for u and d quarks.The nonlocal currents j S,D,V (x) are given in terms of operators based on a separable approximation of the effective one gluon exchange model (OGE) of QCD.These currents read In the above equation we have used ψ C (x) = γ 2 γ 4 ψT (x), Γ f = (1 1, iγ 5 τ ) and Γ D = iγ 5 τ 2 λ a , while τ and λ a , with a = 2, 5, 7, stand for Pauli and Gell-Mann matrices acting on flavor and color spaces, respectively.The function g(z) in Eqs. ( 5) is a covariant formfactor characterizing the nonlocality of the effective quark interaction [39].
The effective action in Eq. ( 2) might arise via Fierz rearrangement from some underlying more fundamental interactions, and is understood to be used -at the mean field level-in the Hartree approximation.In general, the ratios of coupling constants H/G S , G V /G S depend on such microscopic action.For example, for OGE interactions in the vacuum Fierz transformation leads to H/G S = 0.75 and G V /G S = 0.5.However, since the precise form of the microscopic interaction cannot be derived directly from QCD, this value is subject to rather large theoretical uncertainties.In fact, thus far there is no strong phenomenological constraint on the ratio H/G S , except for the fact that values larger than one, are quite unlikely to be realized in QCD since they might lead to color symmetry breaking in the vacuum.We introduce η = G V /G S for the dimensionless vector coupling strength and use it as a free parameter of the model rsponsible for the stiffness of the quark matter EoS at nonzero densities.Details of the values used in the present work will be given below.
To proceed it is convenient to perform a standard bosonization of the theory.Thus, we introduce scalar, vector and diquark bosonic fields and integrate out the quark fields.In what follows we will work within the mean field approximation (MFA), in which these bosonic fields are expanded around their vacuum expectation values and the corresponding fluctuations are neglected.The only nonvanishing mean field values in the scalar and vector sectors correspond to isospin zero fields, σ and ω respectively.Concerning the diquark mean field values, we will assume that in the density region of interest only the 2SC phase might be relevant, thus, we adopt the ansatz ∆5 = ∆7 = 0, ∆2 = ∆.
Next, we consider the Euclidean action at zero temperature and finite baryon chemical potential µ B .For such purpose we introduce different chemical potentials µ f c for each flavor and color, then the corresponding mean field grand canonical thermodynamic potential per unit volume can be written as where we have introduced different chemical potentials µ f c for each flavor and color.The inverse propagator S −1 is a 48 × 48 matrix in Dirac, flavor, color and Nambu-Gorkov spaces.A detailed description of the model and the explicit expression for the thermodynamic potential after calculating the determinant of the inverse of the propagator can be found in Ref. [43].Then, the mean field values σ, ∆ and ω can be obtained by solving the coupled equations In principle one has six different quark chemical potentials, corresponding to quark flavors u and d and quark colors r, g and b.However, there is a residual color symmetry (between red and green colors due to the the ansatz we considered) arising from the direction of ∆ in color space.Moreover, if we require the system to be in chemical equilibrium, it can be seen that chemical potentials are not independent from each other.In general, it is shown that all µ f c can be written in terms of three independent quantities: the baryonic chemical potential µ B , a quark electric chemical potential µ Qq and a color chem-ical potential µ 8 .The corresponding relations read The chemical potential µ Qq , distinguishes between up and down quarks, and the color chemical potential µ 8 has to be introduced to ensure color neutrality.
As we are interested in describing the behaviour of quark matter in the core of neutron stars we have to take into account the presence of electrons and muons, in addition to quark matter.We treat leptons as a free relativistic Fermi gas and the corresponding thermodynamic potential expression can be found in [43].In addition, it is necessary to take into account that quark matter has to be in beta equilibrium with electrons and muons through the beta decay reactions for l = e, µ.Thus, assuming that (anti)neutrinos escape from the stellar core, we have an additional relation between fermion chemical potentials, namely for c = r, g, b, µ e = µ µ = µ l .Finally, in the core of neutron stars we also require the system to be electric and color charge neutral, hence the number of independent chemical potentials reduces further.Indeed, µ l and µ 8 get fixed by the condition that charge and color densities vanish, where the expressions for the lepton densities ρ l and the quark densities ρ f c can be found in the Appendix of [43].In summary, in the case of neutron star quark matter, for each value of µ B one can find the values of ∆, σ, µ l and µ 8 by solving the gap equations (7), supplemented by the constraints of Eq. ( 13) for β− equilibrium and Eqs. ( 14) and ( 15) for electric and color charge neutrality, resp.This allows to obtain the quark matter EoS in the thermodynamic region we are interested in.

C. Twofold interpolation method
When asymptotic forms of the EoS are known, e.g., in the low-density regime of nuclear matter and in the high-density regime of quark matter, but they cannot readily be trusted in the intermediate range of densities where a quark-hadron phase transition is expected, it has been suggested in Ref. [47] to use an interpolation technique (a corrected version has been given in [48]).These works were inspired by the earlier developed interpolation technique [49] for the EoS in the vicinity of the crossover transition at finite temperatures and vanishing baryon densities.
Such an interpolation technique has been used subsequently for the nlNJL EoS (P (µ; η)) which depends on an a priori unknown value of the vector coupling strength parameter η.Since this parameter may actually vary as a function of the chemical potential, it has been suggested in [46] to use the interpolation technique in order to model the EoS with a varying vector coupling strength parameter within certain limits and in a certain range of chemical potentials.
In the present work we want to generalize this approach by using it for two purposes: 1. to model the unknown density dependence of the confining mechanism by interpolating the contribution of a bag pressure B between zero and a finite value at low densities in the vicinity of the hadronto-quark matter transition, and 2. to model a density dependent stiffening of the quark matter EoS at high density by interpolating between EoS for two values of the vector coupling strength, η < and η > .
The resulting doubly interpolated quark matter EoS reads where we have introduced two smooth switch-off functions, one that changes from one to zero at a lower chemical potential µ < with a width Γ < , and one that switches off at a higher chemical potential µ with a width Γ , whereby the corresponding switch-on functions are the complementary ones, The input EoS differ in their η-values, we have taken for those values at low (<) and high (>) chemical potentials η < and η > .Moreover, the bag constant B is introduced to enforce confinement effects in the quark EoS at low chemical potential.
More detailed discussions about the adoption of a bag constant (or bag function) in addition to the chiral quark model pressure see, e.g., Refs.[57][58][59][60].Values of B = 10 − 50 MeV/fm 3 are used in these references and will be taken as a guidance for adjusting this parameter in the present work.For a density-dependent bag constant in hybrid compact star physics applications, see also Ref. [61].
Note that for B = 0 one would not obtain third family solutions within the present approach.

D. Phase transition
In the present work, we shall use a simple Maxwell construction for the phase transition between the EoS for nuclear matter and quark matter phases described above.They shall separately fulfill charge-neutrality and beta equilibrium with electrons and muons.The two distinct phases are then matched according to the Gibbs conditions for phase equilibrium by requiring that temperatures, chemical potentials and pressures of the two phases coincide at the phase transition Technically, we plot the T = 0 isotherms of both phases in pressure over baryon chemical potential and merge them at the crossing point P c = P (µ c ). Outside the phase transition, the phase with higher pressure (lower grand canonical potential) is to be chosen as the physical one.
A more sophisticated construction of the phase transition considers the occurrence of structures (so-called "pasta phases") with an interplay of surface tension, Coulomb energy and charge-screening effects.For a recent work, see e.g.Ref. [62] and references therein.It has been found that for typical values of the surface tension the pasta phase construction yields an EoS very similar to that of the Maxwell construction [63].
The thermodynamics of pasta phases can very well be mimicked by applying interpolating constructions [10,64] for the transition between pure phases.These interpolations differ from the one employed in subsection II C since they are bound to a finite region of chemical potentials starting from the endpoint of the hadronic phase at µ H and ending at the onset of the quark matter phase at µ Q and enclosing the critical chemical potential µ c of the Maxwell construction.Moreover, these mixed phase constructions include a pressure increment ∆P at µ c relative to P c , which is absent for the interpolation in II C which is between EoS of the same phase for different values of one or several free parameters.

E. TOV equations, moment of inertia and tidal deformability
In this subsection, we present the set of equations that are to be solved for obtaining the numerical results on compact star structure and global properties shown in the next section, once the EoS is given.

TOV equations
In order to compute the internal energy density distribution of compact stars and thus derive the massradius relation we utilize the Tolman-Oppenheimer-Volkoff (TOV) equations for a static and spherical star in the framework of general relativity [65,66]: By considering that P (r = R) = 0, M = m(r = R) and P c = P (r = 0) we have the necessary boundary and initial conditions for a relativistic star with mass M and radius R, respectively.Once a central density c is chosen, the TOV solutions provide internal density profiles of a star.By increasing c (or equivalently P c ) for each star, the whole sequence up to the maximum mass, the mass-radius relation can be computed.A similar relation can be derived for the baryon number N B (r), enclosed in a distance r from the center of the mass distribution, Integrating this equation one ontains the baryon number N B = N B (R) of the star, from which its baryonic mass is obtained as M B = m N N B , where m N is the nucleon mass.

Moment of inertia
Another important stellar quantity is the moment of inertia.We compute the relativistic moment of inertia based on the approach presented in [67] For a detailed discussion of the moment of inertia in the slow-rotation approximation, and for the hybrid star case see, e.g., [21,68,69], and references therein.

Tidal Love number and deformability parameter Λ
In this section we briefly describe how to compute the tidal deformability (TD) of a compact star, based on the results of [70][71][72][73][74].We start by considering the dimensionless tidal deformability parameter Λ = λ/M 5 which is computed for small tidal deformabilities.Here λ is the stellar TD and M is the stellar gravitational mass as introduced before.λ is related to the so called love number The TD can be thought of a modification of the spacetime metric by a linear l = 2 perturbation onto the spherically symmetric star, where K (r) = H (r) + 2H(r)Φ (r), primes denoting derivatives with respect to r.The functions H(r), β(r) = dH/dr obey where f = d /dp is the equation of state.The above equations must be solved simultaneously with the TOV equations.The system is to be integrated outward starting near the center using the expansions H(r) = a 0 r 2 and β(r) = 2a 0 r as r → 0. a 0 is a constant that determines how much the star is deformed and turns out to be arbitrary since it cancels in the expression for the Love number.With the definition of the l = 2 Love number is found as where C = M/R is the compactness of the star.

III. NUMERICAL RESULTS AND DISCUSSION
In the present work we consider the case of a Gaussian form factor function g(z) which translates to a Gaussian regulator function g(p) = exp(−p 2 /p 2 0 ) in Euclidean 4-momentum space.The fixed parameters of the quark model are m c = 5.4869 MeV, p 0 = 782.16 MeV, G S p 2 0 = 19.804and H/G S = 0.75.The dimensionless vector coupling η and the bag pressure B are the two parameters that are used for calculating the quark matter EoS that enter the twofold interpolation which is determined by the four parameters µ < , Γ < and µ , Γ for the switch functiuons f < (µ) and f (µ), respectively.In table I we give the values for these three sets of these parameters which will be used in the numerical calculations of this work.I.
In Fig. 1 we show for the example of parameter set 2 of Tab.I the three quark matter EoS (in red thin lines) generated by the nonlocal chiral quark model that are used in the twofold interpolation procedure to obtain the resulting quark matter EoS shown by the thick red dashdotted line.
This procedure is an effective method to reconstruct a chemical potential and thus density dependence of the vector coupling strength η(µ) and the bag pressure B(µ) which determine the stiffening at high densities and the softening due to onset of confinement at low densities, respectively.A Maxwell construction is performed with the hadronic EoS shown by the blue dashed line which results in the black solid line for the quark-hadron hybrid EoS.The three parameter sets are adjusted such that the onset masses M c for the deconfinement phase transition in a compact star lie at 2.0, 1.39 and 1.20 M for set 1, set 2 and set 3, respectively, while the maximum mass on the hybrid star branch exceeds the value of 2.01 M measured for PSR J0348+432 [4].This is achieved by minimal variations in the low-density value η < of the vector coupling and the parameters of the switch functions f < (µ) and f (µ) while keeping the bag constant B = 35 MeV/fm 3 and η > = 0.12 constant.In Fig. 2 we show a comparison of the hybrid EoS resulting from the Maxwell construction with sets 1, 2 and 3.
Fig. 3 shows the squared speed of sound c 2 s = dP/dε as a function of the energy density for sets 1, 2 and 3 which exhibit the regions of the first order phase transition where c 2 s = 0 and fulfills the condition of causality c 2 s < 1 (in units of the speed of light squared).
In Fig. 4 we show a key result of this paper, the massradius relationships for the three hybrid EoS of Fig. 2 as obtained from the solution of the TOV equations.The dotted lines denote the unstable configurations that should not be realized in nature but guide the eye to the corresponding stable hybrid star sequence (third family) disconnected from the neutron star one (second family).The blue and red horizontal bands denote the mass measurement for PSR J0348+432 [4] and PSR J1614-2230 [3], resp.The grey and orange bands labelled M1 and M2 are the mass ranges for the compact stars in the binary merger GW170817 for which Ref. [76] has excluded radii smaller than 10.68 km of 1.6 M stars and Ref. [77] excludes radii exceeding 13.4 km at 1.4 M .The green band denotes the mass range 1.44 ± 0.07 M of PSR J0437-4715, the primary target of the radius measurement by NICER [75].
Fig. 5 shows baryon mass vs. radius for the compact star sequences corresponding to set 1, set 2 and set 3 In Fig. 6 the mass vs. central energy density as solutions of the TOV equations for the three parameter sets given in Tab.I is shown, line styles as in the previous figures.The rising sections of the lines denote stable sequences.
Fig. 7 shows the moment of inertia vs. mass for the three hybrid EoS parametrizations introduced in this work.Unstable configurations are shown by dotted lines.Note that the relationship is multivalued in the mass twin ranges.Hybrid stars on the third family branch have generally a smaller moment of inertia.Measurements of the moment of inertia of compact stars have the power to constrain the existing EoS models.Consider for instance the system PSR J0737-3039 [79] where both compact stars A and B have been observed as pulsars at the time of their discovery.Due to relativistic effects there exists a spin-orbit coupling that could eventually lead to the determination of I A , the moment of inertia of star A in the binary.Since the masses of both A and B stars are already accurately determined to M A = 1.3381 (7) M and M B = 1.2489(7)M [80]and the moment of inertia can be expressed in terms of only mass and radius The dotted lines denote the unstable configurations that should not be realized in nature but guide the eye to the corresponding stable hybrid star sequence (third family) disconnected from the neutron star one (second family).For comparison, the mass-radius sequence for the APR EoS is shown (red solid line) which is a standard EoS for nuclear astrophysics applications, concerning only the second family, without a third family branch.The blue and red horizontal bands denote the mass measurement for PSR J0348+432 [4] and PSR J1614-2230 [3], resp.The green band denoted the mass 1.44 ± 0.07 M of PSR J0437-4715, the primary target of the radius measurement by NICER [75].The grey and orange bands labelled M1 and M2 are the mass ranges for the compact stars in the binary merger GW170817.The hatched regions labelled "GW170817" are excluded by analyses of the merger event: Ref. [76] has excluded radii of 1.6 M stars smaller than 10.68 km, Ref. [77] excludes radii exceeding 13.6 km at 1.4 M , and Ref. [78] gives an upper limit for the maximum mass of nonrotating compact stars of 2.16 M .
and no other EoS parameters, the measurement of I A is of great importance as it allows for the simultaneous determination of mass and radius for the same object thus providing a datum of high discriminating power among all EoS models for compact stars [81,82].Therefore, we show in Fig. 7 the mass of PSR J0737-3039 (A) for which a measurement of I A at the prognosed 10% level could well discriminate between a neutron star or hybrid star case in our example of the set 3 EoS parametrization.
In Fig. 8 the dimensionless tidal deformability parameter Λ vs. compact star mass for the three sets of hybrid EoS parametrizations introduced in this work.The constraint Λ < 800 derived from GW170817 [8] for the mass range 1.16 -1.60 M is fulfilled for the third family solu-      8 for the EoS given by sets 1 -3 compared to the probability contours for the low-spin prior of the LVC analysis of GW170817 [8].
tions of set 2 and set 3 while the second family sequence of set 1 cannot fulfill this constraint.
The second key result obtained for the new class of hybrid star EoS introduced in this work is shown in Fig. 9.It is the plot of the tidal deformability parameters Λ 1 and Λ 2 of the high-and low-mass components of the binary merger, constructed from the Λ(M ) relation of Fig. 8 for the EoS given by sets 1 -3 compared to the probability contours for the low-spin prior of the LVC analysis of GW170817 [8].It shows that a stiff hadronic EoS like the one used in this work would be excluded by the tidal deformability constraint at the 90% confidence level when it is to be used in a binary neutron star merger scenario.The same EoS being part of the three sets of hybrid EoS suggested in this work lead to acceptable scenarios when at least one of the stars in the coalescing binary belongs to the third family branch of hybrid stars.

IV. CONCLUSIONS
We have suggested in this work a twofold interpolation procedure to obtain a class of hybrid compact star equations of state which support the existence of a third family of compact stars composed of color superconducting two-flavor quark matter in the core and a shell of asymmetric nuclear matter described within a relativis-tic meanfield model with nucleonic excluded volume effects.The quark matter equation of state is based on a nonlocal chiral quark model with a covariant separable interaction in the scalar meson and diquark channels to which a vector meson mean field with coupling constant η responsible for the stiffness and a bag constant B for quark confinement are superimposed.As the density dependence of η(µ) and B(µ) is a priori unknown, we use the twofold interpolation method to model the dependence of these parameters as a smooth interpolation between two values: η < ≤ η(µ) ≤ η > and B ≥ B(µ) ≥ 0 for increasing baryochemical potential.The phase transition to the hadronic DD2 p40 EoS is obtained by a Maxwell construction.
We have presented results for this class of hybrid equation of state and the corresponding properties of compact star sequences for three sets of parametrizations which result in a maximum mass of the second family (neutron star sequence) at 2.00, 1.39 and 1.20 M , respectively, and a third family of stable compact stars (hybrid star sequence) separated from the second one by a sequence of unstable configurations so that the phenomenon of mass twin stars is obtained.The hybrid star sequences include the observed maximum pulsar mass of 2.01 M as a necessary constraint for compact star EoS.
It is demonstrated that this advanced description of hybrid star matter allows to interpret GW170817 as a merger not only of two neutron stars but also of a neutron star with a hybrid star or of two hybrid stars.The latter two scenarios can be in accordance with the constraints on compactness from GW170817 when a binary neutron star merger scenario would be ruled out because of a too stiff hadronic equation of state.The NICER experiment on board of the International Space Station [75,83] has the potential to rule out too soft hadronic equations of state and thus support a merger scenario involving hybrid stars from a third family which can be described with the new class of EoS presented in this work.
phase transition to hadronic matter is obtained by a Maxwell construction with the hadronic matter EoS DD2 p40[55], the critical chemical potential µc and critical pressure Pc are obtained.The values of energy density and baryon number density corresponding to the onset of the first-order phase transition are εc and nc, resp.Upon solving the TOV equations with the hybrid EoS, the mass-radius relation for compact stars is obtained which reveals a maximum mass Mmax and a mass Mc at the onset of the phase transition in the center of the compact star.

Figure 1 :
Figure 1: Hybrid star EoS (black solid line) obtained by a Maxwell construction between the quark matter EoS for set 2 (red dash-dotted line) and the density-dependent relativistic meanfield EoS DD2 p40 for nuclear matter in β−equilibrium with electrons and muons.The doubly interpolated quark matter EoS is based on three parametrizations of the nonlocal NJL model: a soft (low vector coupling η) one with confinement (B = 0) at low densities (red dashed line), a soft one without confinement at intermediate densities (red dash-double-dotted line) and a stiff one (high η) at high densities (red double-dash-dotted line).The parameters of the switching functions are given in table I.

Figure 2 :
Figure 2: The equation of state for pressure vs. energy density resulting from the twofold interpolation method and Maxwell construction of the deconfinement phase transition illustrated in Fig. 1 for the parameters of set 2. Results for set 1, set 2 and set 3 are shown as solid, dashed and dash-dotted lines, respectively.For comparison, the APR EoS is shown (orange dotted line) which is a standard EoS for nuclear astrophysics applications

Figure 3 :
Figure 3: The squared speed of sound c 2 s as a function of the energy density for set 1, set 2 and set 3 exhibits the regions of the first order phase transition where c 2 s = 0 and fulfills the condition of causality c 2 s < 1 (in units of the speed of light squared).Line styles as in Fig. 2. For comparison, the APR EoS is shown (magenta dotted line) without deconfinement phase transition

Figure 4 :
Figure 4: Mass vs. radius for sequences of compact stars with the hybrid EoS of this work parametrized by set 1, set 2 and set 3, corresponding to different onset masses for the deconfinement transition.The dotted lines denote the unstable configurations that should not be realized in nature but guide the eye to the corresponding stable hybrid star sequence (third family) disconnected from the neutron star one (second family).For comparison, the mass-radius sequence for the APR EoS is shown (red solid line) which is a standard EoS for nuclear astrophysics applications, concerning only the second family, without a third family branch.The blue and red horizontal bands denote the mass measurement for PSR J0348+432[4] and PSR J1614-2230[3], resp.The green band denoted the mass 1.44 ± 0.07 M of PSR J0437-4715, the primary target of the radius measurement by NICER[75].The grey and orange bands labelled M1 and M2 are the mass ranges for the compact stars in the binary merger GW170817.The hatched regions labelled "GW170817" are excluded by analyses of the merger event: Ref.[76] has excluded radii of 1.6 M stars smaller than 10.68 km, Ref.[77] excludes radii exceeding 13.6 km at 1.4 M , and Ref.[78] gives an upper limit for the maximum mass of nonrotating compact stars of 2.16 M .

Figure 5 :
Figure 5: Baryon mass vs. radius for the compact star sequences corresponding to set 1, set 2 and set 3 parametrizations of the present EoS model.From this figure one can read off what drop in radius can be expected when a transition from the maximum mass of the second family branch to third family branch takes place under baryon number conservation, triggered for instance by mass accretion from a companion star (red arrows).

Figure 6 :
Figure 6: Mass vs. central energy density as solutions of the TOV equations for the three parameter sets given in Tab.I, line styles as in the previous figures.The rising sections of the lines denote stable sequences.

Figure 7 :
Figure 7: Moment of inertia vs. mass for the three hybrid EoS parametrizations introduced in this work.Unstable configurations are shown by dotted lines.Note that the relationship is multivalued in the mass twin ranges.Hybrid stars on the third family branch have generally a smaller moment of inertia.The mass of PSR J0737-3039 (A) is shown for which a measurement of I is expected at the 10% level.

Figure 8 :
Figure 8: Dimensionless tidal deformability parameter Λ vs. compact star mass for the three sets of hybrid EoS parametrizations introduced in this work.The constraint Λ < 800 derived from GW170817 [8] for the mass range 1.16 -1.60 M is fulfilled for the third family solutions of set 2 and set 3 while the second family sequence of set 1 cannot fulfill this constraint within a binary neutron star scenario.

Figure 9 :
Figure 9: Tidal deformability parameters Λ1 and Λ2 of the high-and low-mass components of the binary merger, costructed from the Λ(M ) relation of Fig.8for the EoS given by sets 1 -3 compared to the probability contours for the low-spin prior of the LVC analysis of GW170817[8].