Tidal deformability and other global parameters of compact stars with strong phase transitions
Magdalena Sieniawska Władysław Turczański Michał Bejger Julian Leszek Zdunik
 







Nicolaus Copernicus Astronomical Center, PAS, ul. Bartycka 18, 00-716, Warsaw, Poland
Astronomical Observatory, Adam Mickiewicz University, Poznań, Poland
APC, AstroParticule et Cosmologie, Université Paris Diderot, CNRS/IN2P3, CEA/Irfu, Observatoire de Paris, Sorbonne Paris Cité, F-75205 Paris Cedex 13, France




Global parameters of compact stars with phase transitions
Sieniawska et al.




 

Using parametric equations of state (relativistic polytropes and a simple quark bag model) 
to model dense-matter phase transitions, we study global, measurable astrophysical 
parameters of compact stars such as their allowed radii and tidal deformabilities. 
We also investigate the influence of stiffness of matter before the onset of the 
phase transitions on the parameters of the possible exotic dense phase.

The aim of our study is to compare the parameter space of the dense 
matter equation of state permitting phase transitions to a sub-space compatible 
with current observational constraints such as the maximum observable mass, 
tidal deformabilities of neutron star mergers, radii of configurations before 
the onset of the phase transition, and to give predictions for future observations.

We studied solutions of the Tolman-Oppenheimer-Volkoff equations for a flexible
set of parametric equations of state, constructed using a realistic description
of neutron-star crust (up to the nuclear saturation density), and relativistic
polytropes connected by a density-jump phase transition to a simple bag model
description of deconfined quark matter.  

In order to be consistent with recent observations of massive neutron 
stars, a compact star with a strong high-mass phase transition cannot have a radius 
smaller than 12 km in the range of masses 1.2 - 1.6 M_⊙. We also compare 
tidal deformabilities of stars with weak and strong phase transitions with 
the results of the GW170817 neutron star merger. Specifically, we study 
characteristic phase transition features in the Λ_1-Λ_2 relation, 
and estimate the deviations of our results from the approximate formulæ for 
Λ̃-R(M_1) and Λ-compactness proposed in the literature.  
 We find constraints on the hybrid equations of state to produce stable neutron stars on the twin branch. For the exemplary equations of state most of the high-mass twins occur for the minimum values of the density jump λ = 1.33-1.54; corresponding values of the square of the speed of sound are α = 0.7-0.37. We compare results with observations of gravitational waves and with the theoretical causal limit and find  that the minimum radius of a twin branch is between 9.5 and 10.5 km, and depends on the phase transition baryon density. For these solutions the phase transition occurs below 0.56 [fm^-3]. 


Tidal deformability and other global parameters of compact stars 
with strong phase transitions
    M. Sieniawska1Corresponding author: msieniawska@camk.edu.pl
       W. Turczański
       2
       M. Bejger3,1 
       J. L. Zdunik1
       
    Received ...; accepted ...
===============================================================================================================================================






§ INTRODUCTION




The interior composition of  neutron stars (NSs) is still not fully known today. Neutron stars
are so compact that the density in their central regions exceeds by far the
density of the atomic nuclei. At such extreme conditions matter may exist in a form
impossible to obtain and study in  terrestrial laboratories.  One of the
possibilities of  very dense matter is the deconfinement of  quarks. 

Astrophysically, we study the dense matter in the interiors of NSs by measuring
their parameters (masses, radii, etc.) and compare them to the theoretical
models of a structure and equation of state (EOS). The multiple EOS models, which include a deconfined quark-matter segment, were recently proposed in <cit.>, <cit.>, <cit.>, <cit.>, <cit.> and <cit.> <cit.>. A defining characteristic of these EOS  is  the phase transition between  normal matter and  quark matter. To be consistent with the observational constraints, all the EOS have to satisfy the maximum mass M_max criterion, which is now established to be approximately 2 M_⊙ (see  for PSR J0348+0432 and  for PSR J1614-2230, as well as  for the re-evaluation of the mass of the latter).

Creation of a new quark phase—the appearance of a new phase core in
the centre of a NS—may lead to a destabilisation 
of a part of the sequence of
NS configurations (where the sequence is labelled by a central EOS parameter,
such as central pressure, P_c). If the stability is regained further up the
sequence, it is considered  a detached (twin) branch of the compact stars (see right panel of Fig. <ref> for a schematic depiction of this
situation in the M(R) plane). In principle, both mass maxima (twin solutions) may exist around the maximum mass M_max≈2 M_⊙, and  is called the high-mass twins solutions. 

In this article we study the properties related to the twin branch of the NS using 
a simplified parametric EOS that allows us to analyse the space of NS solutions. 
Specifically, we study the influence of a strong (destabilising) phase transition 
on the tidal deformability of the NS, and apply these results to the recent measurements 
obtained by the LIGO-Virgo Collaborations from the GW170817 event <cit.>. For the recent studies concerning the tidal deformability see 
<cit.>, <cit.>, <cit.>, <cit.> and <cit.>.
We also aim to constrain how the part of the EOS before the phase
transition influences the parameters of the EOS after the phase
transition. In other words, we aim to put constraints on the observable
parameters of the NS (such as the radius that is potentially a measurable
parameter) for the masses below that at which the phase transition occurs, and
investigate how the shape of the M(R) is related to the properties of the
twin branch. 

The article is composed as follows. In Sect. <ref> we describe 
the parametric EOS and the methods of obtaining the solutions of the hydrostatic 
equilibrium equations for the NSs. In Sect. <ref> we present 
the properties of NS sequences compatible with existing observational constraints, 
such as the M_max requirement and the tidal constraints of GW170817. 
Section <ref> contains discussion and summary.   




§ EQUATIONS OF STATE AND METHODS




In order to survey the space of solutions corresponding to M(R) sequences with 
high-density phase transitions, we employ a conservative approach and use the following simplified, parametric
EOS. We assume the knowledge of the low-density part of the EOS and adopt the SLy4 EOS 
description of <cit.> and <cit.> up to a
baryon density n_0, comparable to and typically larger than the nuclear
saturation density (n_s≡ 0.16 fm^-3). At n_0 a relativistic
polytrope <cit.> replaces the tabulated SLy4 EOS. The definition of the
pressure P and the energy-density ρ c^2 are standard,  

    P=κ n^γ,  ρ c^2 = P/γ -1 + nm_bc^2,
where κ is the pressure coefficient, γ is the index of the
polytrope, and m_b is the mass of the baryon
in this phase. The index γ is a parameter of choice; consequently, by
demanding the chemical and mechanical equilibrium at n_0, κ and
m_b are fixed. The polytrope ends at a density n_1 > n_0, and is
connected to a simple bag EOS <cit.>, characterising 
the quark matter. We use a linear pressure-density relation of <cit.>, 

    P(ρ)=α(ρ - ρ_*)c^2,    n(P) = n_* (1+α+1/αP/ρ_*c^2)^1/(α+1),
  
with α denoting square of the speed of sound in a quark matter, and ρ_* 
and n_*
the
energy density and baryon density of this matter at zero pressure, respectively. 
We assume that at the polytrope/bag boundary matter is softened by the
first-order phase transition defined by a density jump 
λ=n_2/n_1. Maxwell construction at this point results in a corresponding
mass-energy density jump ρ_2/ρ_1=λ+(λ-1) P_1/ρ_1 c^2. The values of n_2, ρ_2 determined by the definition of the transition point (n_1,ρ_1,P_1), together with a given λ yield, from Eq.<ref>, the values of ρ_*=ρ_2-P_1/α c^2 and n_*. 
The schematic pressure–density relation for such an EOS is
presented in the left panel of Fig. <ref>. The initial parameter
ranges are shown in Table  <ref>. 
 


   


Given the EOS, we solve the equations of hydrostatic equilibrium for a spherically 
symmetric distribution of mass <cit.>, 

    dP(r)/dr   =   -G/r^2(ρ(r)+P(r)/c^2)
    (M(r)+4π r^3 P(r)/c^2) 
       ×   (1-2GM(r)/c^2r)^-1,
 
supplied with the equation for one of the metric functions, 

    dν(r)/dr= -(2/P(r)+ρ(r)c^2) dP(r)/dr,
 
for a spherically symmetric metric of the form

    ds^2=e^ν(r) c^2 dt^2 - dr^2/1-2GM(r)/rc^2 - r^2(dθ^2 + sin^2θ dϕ^2),
 
and the equation for the total gravitational mass inside the radius r:

    dM(r)/dr=4 πρ(r) r^2.

We solve an additional equation for the tidal deformability of the
star, defined as 

    λ_td = 2/3R^5k_2.
 
It represents the reaction of the star on the external tidal field (such as
that in a tight binary system; e.g. ). Influence of the tidal field is obtained in the
lowest order approximation, by calculating the second (quadrupole) tidal Love
number k_2 <cit.>  

    k_2    = 8/5x^5 (1-2x)^2 ( 2-y + 2x(y-1))   
    ( 2x ( 6 -3y +3x(5y-8)) . 
       + 4x^3 ( 13 -11y + x(3y -2) + 2x^2(1+y) ) 
       .+ 3(1-2x)^2 ( 2 - y + 2x(y-1)) ln(1-2x) )^-1,
 
with the star's compactness x=GM/Rc^2, and y the solution of 

    dy/dr   = -y^2/r-1+4π G r^2/c^2(P/c^2-ρ)/(r-2GM(r)/c^2)y  
       + (2G/c^2(M(r) + 4π r^3 P/c^2)/√(r)(r-2GM(r)/c^2))^2 + 6/r-2GM(r)/c^2
       - 4π G r^2/c^2/r-2GM(r)/c^2(5ρ+ 9P/ c^2
    +( ρ   + P/c^2)^2 c^2/ρ dP/dρ),

evaluated at the stellar surface <cit.>. In the
following we  use the normalised value of the λ_td parameter, 

    Λ = λ_td( GM/c^2)^-5.
 
Concerning the current sensitivity of the detectors <cit.>, what
is actually measured is the effective tidal deformability Λ̃,
defined as 

    Λ̃ = 16/13(M_1 + 12 M_2) M_1^4 Λ_1 + 
            (M_2 + 12 M_1) M_2^4 Λ_2/(M_1 + M_2)^5,
 
with M_1 and M_2 denoting the component masses. 
The equations are solved using a Runge-Kutta fourth-order numerical scheme with a
variable integration step (see e.g.  for details) for a
range of central parameters of the EOS (central pressure P_c) to obtain
global parameters of NSs: their masses and radii in the form of the mass-radius
M(R) sequence. 




§ FEATURES OF THE TWIN BRANCH




From all the possible M(R) sequences resulting from the parameter set
presented in Table <ref>, we select those consistent with
the following criteria:  

  
  * the maximum mass M_max is equal to or larger than the largest currently measured 
  NS mass: M_max≥ 2 M_⊙ <cit.>; 
  
  * strong (destabilising) phase transition in the EOS, i.e. an instability near the phase transition point that results in the local minimum 
  in the M(R) curve and the presence of the twin branch; 
  
  * unless stated otherwise, the maximum mass is located on the twin branch, i.e. M_max≥ M_ph (e.g. Sect. <ref> for twin sequences with M_max≤ M_ph); 
  
  * the EOS is causal (the speed of sound √(∂ P/∂ρ) is smaller 
  than the speed of light in vacuum).   




 §.§ Mass-radius diagram: minimum radius for the second branch 
and the tidal deformability
 














We are interested in establishing a lower limit on the radius R at densities
below the phase transition point for the EOS compatible with the twin branch
scenario. Since it is likely that a large fraction of the NS population has 
masses below M_max=2 M_⊙, the future measurements of their radii
will be indicative of whether the EOS is able to support a phase transition
at higher densities (high-mass twins scenario) in a way compatible with the observations of the massive NSs.
In order to obtain the lower bound on the radius, we adopt the α
parameter in the quark phase equal to its extreme value, α=1.
Characteristic shapes of the M(R) relations are presented in
Fig. <ref>, where  selected 
regions of the Λ parameters in the M(R) plane are also shown. The values 
of Λ decrease with M and increase with R thanks to the strong scaling with 
these values. A limit of 900 corresponds to the lower edge of the orange region. 
In principle, we can construct a specific EOS corresponding to limiting  
R_ph or M_ph still consistent with this bound. We conclude that the minimum radius R on the polytropic branch of the M(R) sequence that still allows for the twin branch scenario cannot be smaller than 12 km with the adopted SLy4 crust, for the astrophysically interesting range of masses (≃1.0 - 2.0 M_⊙). Softer pre-phase-transition polytropes which give R<12 km 
result either in M_max<2 M_⊙, or in M_ph > M_max.  In
addition to the polytropic index γ, factors that determine the radius
for stars below the high-density twin branch segment are the densities n_0 at
which the polytrope connects to the SLy4 EOS (low-density SLy4-polytrope
connection results in larger radii) and its ending point n_1.   This simple
exercise shows that the measurements of radii in this range of masses disclose
information pertinent to our understanding of the EOS at densities above the
nuclear saturation density. 

The thick dashed black line in Fig. <ref> 
approximates the minimum radius of the configurations before the phase transition, 

    R = R_SLy4(1M_⊙) + (M/M_⊙-1) km,
 
with R_SLy4(1M_⊙)=11.85 km. For comparison, we also plot the SLy4 EOS
curve, which resides entirely in the quark-core part of the M(R) diagram.
Figure <ref> shows the relation between the central pressure
P_c for configurations with M_ph and R(M_ph) at the end of the polytropic
branch; the left-side boundaries of the scatter plots correspond to the minimum radius 
line in Fig. <ref>; for example, stars at the polytropic branch 
with radii smaller than 12 km and corresponding masses smaller than 1.2 M_⊙ 
cannot exist with central pressures larger than approximately 7.5× 10^34 
dyne/cm^2. The relation between M_max and R(M_ph) is shown in Fig. <ref>. Small radii at the transition point support lower maximum masses on the twin branch.
    








 §.§ Tidal deformability of stars with weak and strong phase transitions




We now study the tidal deformability Λ parameters
(Eq. <ref>) using a subset of the parametric EOS. The ranges of 
astrophysically interesting Λ parameters are indicated in
Fig. <ref>. We note that it is virtually impossible 
to define a clear boundary between different Λ ranges on the M(R) plane. 
As the value of Λ depends sensitively on the EOS of dense matter, it cannot be 
treated as a simple function of M and R only. This is
particularly visible for the comparison of the SLy4 EOS in
Fig. <ref>. Here, depending on the EOS, configurations 
with the same compactness GM/Rc^2 but visibly different Λ occupy the 
same position on the M(R) diagram, e.g. for M=1.32 M_⊙ and R=11.75 km, 
the end of the magenta segment yielding Λ=450 is located in the middle of 
the Λ=300-450 range for the parametric EOS. 

Figure <ref> shows a M(R) diagram with the EOS
  colour-coded   according to their values of Λ(M_1.4). Here
the situation is even less clear in terms of defining the Λ on the M(R) 
plane, also because of the non-monotonic behaviour of M due to the strong
phase transition: a given EOS may yield several values of Λ(M_1.4)
corresponding to different R. We conclude that representing
Λ(M_1.4) on the M(R) diagram (as proposed by e.g. ) 
is not the best way of uncovering the relation between the Λ and other 
EOS functionals. Our results are based on a relatively simple model of the EOS; adding 
more degrees of freedom will only further complicate this picture.  

Using the selected EOS from Fig. <ref> as an example, we
discuss the tidal deformability Λ for the strong-phase transition EOS
by applying a measurement of the chirp mass, 

    ℳ= (M_1M_2)^3/5/(M_1 + M_2)^1/5,

and the component mass ranges from the one binary NS merger observed so far in gravitational waves, 
the GW170817 event <cit.>. We adopt the central value of the 
chirp mass, ℳ=1.188 M_⊙, and component masses in the low-spin 
prior estimation case: M_1 ∈ (1.36, 1.60), M_2 ∈ (1.17,
1.36) <cit.>. An improved analysis that gives consistent results is described in <cit.>, where revised values for the low-spin prior of the chirp mass and component masses are ℳ=1.186 M_⊙, M_1 ∈ (1.36, 1.60), and M_2 ∈ (1.16,1.36). A recent correction of the misprint reveal that Λ̃≲ 900 (Table IV in ), instead of Λ̃≲ 800 (). Nevertheless, due to the very small differences between initial and re-evaluated estimations of the mentioned values, and the fact that their estimation errors overlap (in the case of tidal deformability values, errors of the estimations are considerable and depend on the waveform model, e.g. for the symmetric TaylorF2 waveform model Λ̃ = 340^+580_-240 and for HPD PhenomDNRT Λ̃ = 300^+430_-230, according to ), we decided to use the results from the detection paper. 

To visualise the influence of a phase transition softening on Λ, we
manufacture a selection of the EOS that differ in the size of the density jump
λ=n_2/n_1, and otherwise share their parameters (M(R) diagrams 
shown in Fig. <ref>). The phase transition
from the polytropic segment of γ=4.5 to the linear EOS of α=1 occurs
at n_1=0.335 fm^-3 (crust-core transition from the SLy4 EOS to a
polytropic segment occurs at n_0=0.21 fm^-3). The parameters are selected
such that the onset of the softening occurs within the range of masses
estimated for the GW170817 event <cit.>. 





In Fig. <ref> we plot the Λ_1-Λ_2 relations for 
the component masses of the GW170817 event by assuming the estimated central value of 
the chirp mass, ℳ=1.188 M_⊙. Depending on the strength of the phase 
transition (measured by the size 
of the density jump λ), the Λ_1-Λ_2 relation may 
exhibit a non-monotonic behaviour, as shown in the left panel of 
Fig. <ref> for the blue curve. This feature is a direct 
consequence of a non-monotonic behaviour of the gravitational mass M as a function 
of the tidal deformability Λ (right panel of Fig. <ref>), 
which is connected to the instability caused by a sufficiently strong phase transition. 
Detecting this feature in the incoming observations of NS mergers will clearly signal 
a dense-matter softening in the range of central parameters corresponding to NS 
merger component masses.  





Additionally, we produced a set of four strong phase-transition M(R) 
curves, displayed in Fig. <ref>: their common parameters are γ=4.5,
n_1=0.335 fm^-3, λ=1.7, α=1. The EOS differ only in the
values of the SLy4-polytrope density matching point, n_0=0.16, 0.185, 0.21,0.235 fm^-3, respectively. These n_0 values were selected so that the phase
transition mass M_ph is placed below, above, or within the ranges of the
component masses for the low-spin prior estimation of the GW170817 event 
<cit.>. The Λ_1-Λ_2 relations are displayed 
in Fig. <ref>. Specifically, the green and blue curves are mirror 
rotated versions of themselves, reflecting the fact that the phase transition 
occurs in the range of M_1 or M_2, respectively. 

When both configurations (twin case) are present in the relevant range of masses (e.g. blue and green curves in Fig. <ref> and in the right panel of Fig. <ref>), the twin branch gives more compact M(R) configurations in comparison to
those on the polytropic branch, their Λ values are in general smaller,
and thus  preferable from the point of view of the GW170817 observations
<cit.>. Likewise, larger radii before the phase transition onset
(e.g. red curve in Fig. <ref>) are generally less favoured  by
the observation of GW170817. This may put limits on the proposed high-mass twin
branch solutions in which the onset of a phase transition occurs near 2M_⊙ <cit.>.
Figure <ref> visualises the requirement for a high-mass twin solution and the radius:  the higher the M_ph, the larger
the radius at M_ph. For M(R) dependencies presented in Fig. <ref> in the very large range of masses radius increases
with mass;  this property is typical for a very stiff EOS. As a consequence, high-mass twin configurations in the
range of GW170817 masses (before the phase transition) have larger radii and 
larger Λ than lower mass twins. 











The twin branch segments lie in the region compatible with the SLy4 EOS.
Even though the magenta EOS (n_0=0.235 fm^-3) and the blue EOS
(n_0=0.21 fm^-3) mass-radius relations cross the SLy4 EOS M(R)
sequence, their Λ values differ in these points, as shown in the right
panel of Fig. <ref>. For the SLy4-magenta EOS crossing at M≈
1.17 M_⊙, the ΔΛ≈ 120 (to be compared with the value
of Λ_Sly4≈ 970 at this point) and ΔΛ≈ 35 for
the SLy4-blue EOS crossing at M≈ 1.51 M_⊙ (Λ_Sly4≈
190 at this point). Generally, configurations with the same M and R (same
compactness M/R) have notably different tidal deformabilities depending on
their EOS. We plot the set of configurations from
Fig. <ref> on the Λ-compactness x=GM/Rc^2 plane
(Fig. <ref>), and compare the results with the fitting formula
proposed in <cit.> and recently updated in <cit.>,
x_fit=∑_i a_i(log(Λ))^i (see Eq. 78 and Fig. 15 in
). The deviation from the  fit is, for
the EOS used here, typically larger than 5% and increases above 10% for
Λ>100. The best fit to our set of EOS 
has the following parameters: a_0 = 0.353, a_1 =
-0.0359, and a_2 = 0.000790 (indicated by  a green line in
Fig. <ref>).







Figure <ref> shows the resulting Λ̃
(Eq. <ref>) in a full range of the GW170817 event low-spin
component masses (1.36-1.60 M_⊙ for M_1) as a function of the radius
R(M_1) of the more massive component for the selected strong phase transition
EOS. These results are compared with the approximate formula for
Λ̃ proposed in <cit.>, based on a set of the
`standard' realistic EOS used in <cit.>. The twin branch spans
a considerably wider range in R for the allowed component masses of GW170817
than the standard EOS, e.g. the SLy4 EOS, for which the radii are
relatively constant. This feature results in a much broader range of allowed
Λ̃ values in comparison to the standard EOS: for comparison 
see Fig. 1 of <cit.>, where Λ̃ varies to a much
smaller extent than here. The behaviour of the red EOS (but also the violet EOS, in
which the phase transition occurs below the lower bound on the M_2 mass) is
qualitatively similar to the SLy4 EOS, but since the radii of the  components  change
more than in the case of the SLy4 EOS in the appropriate mass ranges, the
response of the star to the tidal deformation is larger, as the Λ
paramteter is a sensitive function of the radius. We expect a similar behaviour 
for any other EOS with a similar M(R) sequence. The deviation from the fit is of
the order of 100 in Λ̃ and of several hundred  metres in
R(M_1). The character of the curves depends on whether the phase transition
occurs within the range of M_2 or M_1 (blue and green lines). In the first
case, we note a discontinuity in Λ̃ related to the
phase-transition induced instability; in the second case, however, the R(M_1)
also behaves non-continuously because  a range of M_1 for
which two solutions with different values of  R is possible.  








 §.§ Features of the quark twin branch for given polytropic outer layer



In the following we present our study of the polytrope-quark matter EOS in order to connect 
the properties of the quark part of the M(R) diagram to its less-dense, outer part 
i.e. the polytrope. The last M(R)
configuration of the polytropic EOS (corresponding to the transition point) 
is denoted by C_ph. For the
sake of presentation we  use a specific C_ph configuration with 
M_ph=2  M_⊙. It was selected because of its relation
to the current maximum observed mass, which we  assume to be equal
M_max=2 M_⊙. 



Figure <ref> shows a complimentary set of features. Here we
want to study the influence of the C_ph configuration on the twin  branch.
We fix the α parameter to the maximum (α=1). The stiffness of the polytropic
EOS  (related to the value of γ) is directly connected with the derivatives dM/dρ_c or dM/dR at C_ph, and  has a direct impact on the presence of the twin branch.  Stiffer EOS at C_ph
(higher γ resulting in higher value of dR/dM) favours the presence and size of the twin branch. Additionally, we investigate the point before the onset of phase transition. Figure <ref> shows the α versus P relation at the end of the polytope (related to the n_1 parameter from Table <ref> for each EOS). The colours correspond to the M(R) relations in Fig. <ref>.
For illustrative purposes to investigate the conditions at which the stable twin
branch may be present, we select several mass-radius relations with M_max >
M_ph (configurations denoted dark red, DR, and dark blue, DB), M_max =
M_ph (red, R, and blue, B), M_max < M_ph (orange, O, and cyan, C), and one sequence without a twin branch (yellow, Y).  The parameters of these
configurations are given in Table <ref>.


   
      


Figure <ref> demonstrates the behaviour of the M(R) relations
for the selected values of the speed of sound parameter α (right panel)
and the density jump parameter λ (left panel). Only critical values of
α and λ parameters for which stable twin branches can be
formed are shown. As expected, the value of the sound velocity of a quark matter √(α) influences the stiffness of the
core and ability to stop gravitational collapse due to instability resulting from softening by the phase transition. On the other hand,  the stiffness λ influences the size of the instability region between polytropic one-phase stars and hybrid stars with a quark core. For the EOS with smaller radii, at the point of
the phase transition, it is harder to produce a stable twin branch. This means that all twin solutions for the R and Y EOS occupy a similar area on the
mass-radius plane, while twin branches for B and DB span a much wider space. From
the observational point of view, according to  Figs. <ref>
and <ref>, observations of the NS with smaller radii  would put stronger constraints on the parameters of the possible phase transition to quark matter and the quark matter EOS, e.g. values of α and λ, as well as
the macroscopic properties of the NSs, e.g. the  minimum  radius R_min(M_max), 
corresponding to the maximum mass on the twin branch. 


   


The onset and parameters of the polytrope have a crucial influence on the radius
for configurations below 2 M_⊙ and on the behaviour of the twin branch. The
smaller the radius (the higher the value of n_1, the smaller the value of
γ), the higher the value of α needed to yield a stable twin branch (compare the lowest red and blue sequences in the left panel of
Fig. <ref>, where lower limit for the R EOS equals α=0.78 for λ=1.60, and α=0.60 for λ=1.90 for the B EOS).  A similar situation, but for the λ parameter, is
presented in the right panel of the same figure: stable twin branches for the DB EOS can
be formed for values of λ between 1.58 and 2.5, while for the Y EOS the allowed 
λ range is much smaller, and lies between 1.34 and 1.48 (both for α=1).

Figure <ref> illustrates the conditions for a stable twin branch.
For this, local mass minimum and maximum, M_min and M_max, is required.
A marginal case corresponds to M_min=M_max. These configurations
correspond to the zeros of the approximately linear relation between M_max
and the difference M_max-M_min in Fig. <ref> for the values
of n_0, n_1, and γ as for the B EOS. We note that for smaller values of λ: i) a smaller value of α is required to produce a stable twin branch, ii) higher values of M_max are reached, iii) (M_max - M_min) obtains higher values.


   


For a given C_ph configuration we can plot a λ(α) relation, defining the onset of a twin branch, from
the crossing point of the curve plotted in Fig. <ref> with
a horizontal line (y=0). For a given λ this condition defines 
the critical value of α. For the presence of a new, stable twin branch,
matter in a dense phase should be sufficiently stiff: α>α_crit. 


Another astrophysically important condition corresponds to the presence of stellar configurations
on the twin branch with baryon masses higher than the maximum mass of one-phase star M_ph.
In this case there is a natural evolutionary track which allows for the
creation of twin stars. The condition M_B,max=M_B,ph defines for each λ the  limiting
value of α - α_Mmax.
This condition is stronger than just the presence of stable twin branch  and
α_Mmax>α_crit.
Strictly speaking this condition corresponds to the comparison of baryon masses. However, the 
relative difference in the gravitational mass of normal and twin sequence for the same M_B is of the order of ∼10^-3, 
and within this accuracy we can consider the condition M_max=M_ph=2M_⊙, which is equivalent to the crossing
points of curves presented in Fig. <ref> with vertical line M_max=2 M_⊙. Conditions for the α and λ parameters for which a stable twin branch is possible are shown in Fig. <ref>. 
A schematic explanation of this plot is presented in Fig. <ref>. 
Stable twin solutions are not possible to the left of the dotted vertical lines, nor below the solid curves. Between the solid curves and the dashed lines, a stable hybrid branch can be formed with M_max < 2M_⊙. Dashed lines correspond to the situation in which M_ph=M_max=2M_⊙, 
while the regions between the dashed and dotted curves correspond to M_max>M_ph. 

We also indicate a region of the  “Both” class proposed by <cit.>: it lies between the vertical dotted and vertical solid lines. These configurations are not immediately destabilised after the phase transition;  both maxima in the mass-radius relation occur for a finite size quark-phase core. For the Y EOS, all twins are of the “Both” class. Limiting λ for the DB, B, and C EOS at which “Both” configurations can still be formed ≈1.67, while for the Y, O, R, and DR EOS it equals ≈1.70. A clear separation between these two types of EOS is, in this context, the effect of the separation of the radii at the phase-transition point.









Parameter n_1 decides where the transition point C_ph occurs. Here we
investigate the properties of the twin branch with respect to the density at which
the polytrope ends.

For the set of EOS from Table <ref>, we changed
parameters n_1 and chose only the configurations with stable twin branches, as
shown in Fig. <ref>. The configurations with larger radii at the polytropic 
branch (the C, B, DB EOS) span a wider range of n_1 
for which our assumptions are fulfilled (see also Fig. <ref>). 
This kind of EOS support situation, where the phase transition point occurs
at low n_1 and where very massive, stable NSs (with M_max even higher than 3
M_⊙) can exist on the quark branch. We note that all possible
R_min(M_max) for the DR EOS are in the range ≈ 9.8-10.5 km, while for the C
EOS this range is much wider: ≈ 9.7-14.8 km (shown in
Fig. <ref> as thick lines and in Fig. <ref> as a
function of n_1). We compare these results with the causal limit (see e.g.
Appendix D in ):

    2GM/Rc^2≈ 0.71.

The maximum mass M_max is then related to the density ρ_* as

    M_max = 2.116 · (ρ_*,15)^-1/2 M_⊙,

where ρ_*,15 = ρ_*/10^15 g/cm^3 (see Eq. <ref>). 
The M(R) relations for all the EOS are limited by Eq. <ref> and asymptotically converge to the causal limit 
with decreasing n_1.





Minimum radii corresponding to M_max are indicated  by thick lines in Fig. <ref>, with R_min between 9.5 and 10.5 km, depending on the EOS. For low values of n_1 the radii decrease approximately along the causal limit, reach a value of R_min, and become bigger again for high values of n_1. 

We also investigate our results in the context of Eq. <ref>, as  shown in Fig. <ref>. For the configurations with dominant contribution of the quark phase (like the C EOS), the results converge to the causal limit. If the phase transition starts at higher densities (corresponding to ρ_*≈ 0.9× 10^15 g/cm^3), calculated curves diverge from the causal limit. The divergence is especially strong for the O, R, and DR EOS, for which twin branches can be formed for relatively small radii at the transition point. According to the recent analysis of the GW170817 event in <cit.>, these EOS are more favourable than C, B, and DB, due to the smaller radii for configurations without phase transition in the centre  and with masses relevant to GW170817 (M∼ 1.2 - 1.6 M_⊙).

On the twin branch the minimum radius of the hybrid star  depends on the density at phase transition n_1 (see Fig. <ref>). For the adopted EOS models the maximum n_1 for which a stable, twin branch exists is about 0.56 [fm^-3].







§ CONCLUSIONS




In our work we consider stable NSs configurations, described by the EOS constructed  using the SLy4 prescription of the crust, joined with the relativistic polytrope, which ends in the transition point where the quark matter phase occurs. This internal structure may reflect as a twin branch in the mass-radius relation, which means that objects with  different R might have the same M.

Recent detection and improved analysis of the gravitational wave event, originated as a NS-NS merger and named GW170817, delivered information about tidal deformabilities of the compact objects. By using the multiple hybrid EOS, we show that the behaviour of the Λ(M,R) is not straightforward, as regions with different Λ values occupy the same space in the M(R) plane. We point out that representing Λ(M_1.4) on the M(R) plane, as  was recently used in the literature, is not optimal for examining the M(R) dependency of the tidal deformability on the EOS. We also show that the non-monotonic behaviour of Λ_1 - Λ_2 relation depends on the size of the density jump λ, and directly reflects the non-monotonic behaviour of M(Λ) function. An even more complicated situation occurs when the dependency on the SLy4 crust-polytrope connection point is included.

We find that objects on the twin branch are more compact (in comparison to the polytropic branch) and give smaller Λ values. Such solutions are more favoured from the point of view of the GW170817 analysis. In the case of a high-mass twin branch (where the transition point occurs around 2 M_⊙), stars in the range of GW170817 masses have larger R and larger Λ than lower-mass twins. Because the NS with the same compactness can have different tidal deformability (due to the different EOS), we show that the deviation from the Λ - GM/Rc^2 fit proposed in the literature, for the broad range of the hybrid EOS, is around 5% for Λ < 100 and increases to above 10% for the higher values of Λ. We find a similar inaccuracy in the fit on the Λ̃ versus R(M_1) plane, where discontinuous relations for the hybrid EOS span a wider range in R than the standard ones (e.g. the SLy4 EOS).

Finally, we focus on the influence of the transition point on the twin branch. Approximately linear relations for M_max - M_ph versus M_max (for a fixed λ and changeable α) are used to establish conditions for a density jump and speed of sound, that allow for the creation of different types of twin configurations (i.e. “Both” NSs). For the exemplary EOS, most of the twins with M_max>M_ph can be formed for at least λ = 1.33-1.54 and corresponding values of α = 0.7 - 0.37, while “Both” stars occupy the λ≲ 1.7 range. Generally, the smaller the R on the polytropic branch, the more fine-tuned the conditions need to be  to produce a stable twin branch. Our whole set of solutions asymptotically converges to the limit of the causal EOS with  decreasing values of n_1, and it can be approximated by the causal limit when the central density at the beginning of the twin branch ρ_*≲ 0.9× 10^15 g/cm^3 (when contribution to the EOS of the quark phase is dominant). We also show that the minimum radius on the twin branch has a size between 9.5 and 10.5 km and depends on the n_1 value. We note that analyses were performed for exemplary EOS, and at most two varying parameters at once from the parameter set (n_0, n_1, γ, α, λ), when other parameters were constant. This strategy simplified investigation and allowed us to determine the constraints mentioned above. Nevertheless, a full analysis is needed and is dedicated to  future studies. 


We acknowledge support from the Polish National Science Centre via grant nos.
2014/13/B/ST9/02621, 2014/14/M/ST9/00707, and 2016/22/E/ST9/00037, and from the
European Union’s Horizon 2020 research and innovation programme under grant
agreement no. 653477. WT acknowledges the support of the Nicolaus Copernicus
Astronomical Center 2015 summer studies programme.





[Abbott et al.(2017a)]Abbott2017a
Abbott, B. P., Abbott, R., Abbott, T. D., et al., 2017a, Physical Review
Letters, 119, 161101 

[Abbott et al.(2018a)]Abbott2018a
Abbott, B. P., Abbott, R., Abbott, T. D., et al., 2018a, arXiv:1805.11581 

[Abbott et al.(2018b)]Abbott2018b
Abbott, B. P., Abbott, R., Abbott, T. D., et al., 2018b, arXiv:1805.11579

[Alford et al.(2013)]Alford2013 
Alford, M. G., Han, S. & Prakash, M. 2013, Phys. Rev. D 88, 083013

[Alford & Han (2016)]Alford2016 
Alford, M. & Han, S. Eur. Phys. J. A (2016) 52: 62.

[Alford & Sedrakian(2017)]Alford2017 
Alford, M. G. & Sedrakian, A. 2017, Phys. Rev. Lett. 119, 161104

[Annala et al.(2017)]Annala2017
Annala, E., Gorda, T., Kurkela, A., & Vuorinen, A. 2018, Phys. Rev. Lett., 120, 172703

[Antoniadis et al.(2013)]Antoniadis2013 
Antoniadis, J., Freire, P. C. C., Wex, N., et al. 2013, Science, 340, 448

[Benić et al.(2015)]Benic2015
Benić, S., Blaschke, D., Alvarez-Castillo, D. E., Fischer, T., & Typel,
S. 2015, , 577, A40 

[Buballa et al.(2014)]Buballa2014
Buballa, M., Dexheimer, V., Drago, A., Fraga, E., Haensel, P., Mishustin, I., Pagliara, G., Schaffner-Bielich, J., Schramm, S., Sedrakian, A. & Weber, F.  2014, J. Phys. G41, 123001 

[Chodos et al.(1974)]Chodos1974 
Chodos, A., Jaffe, R. L., Johnson, K., Thorn, C. B., & Weisskopf, V. F. 1974, 
, 9, 3471 

[Christian et al.(2018)]Christian2018 
Christian, J-E., Zacchi, A. & Schaffner-Bielich, J. 2018, Eur. Phys. J. A 54: 28

[Demorest et al.(2010)]Demorest2010 
Demorest P.B., et al., Pennucci T., Ransom S.M., Roberts M.S.E., Hessels
J.W.T., 2010, Nature, 467, 1081

[Douchin & Haensel(2001)]DouchinH2001
Douchin, F., & Haensel, P. 2001, , 380, 151

[Farhi & Jaffe(1984)]Farhi1984
Farhi, E., Jaffe, R.L., 1984, , 30, 2379

[Flanagan & Hinderer(2008)]FlanaganH2008 
Flanagan, É. É., & Hinderer, T. 2008, , 77, 021502 

[Fonseca et al.(2016)]Fonseca2016 
Fonseca, E., Pennucci, T. T., Ellis, J. A., et al. 2016, arXiv:1603.00545 

[Haensel et al.(1999)]Haensel1999 
Haensel, P., Lasota, J. P., & Zdunik, J. L. 1999, , 344, 151 

[Haensel & Pichon(1994)]HaenselP1994 
Haensel, P., & Pichon, B. 1994, , 283, 313 

[Haensel, Potekhin & Yakovlev(2007)]HPY2007book 
Haensel, P., Potekhin, A.Y., Yakovlev, D.G., Neutron Stars 1 Equation of
state and structure (Springer, New York 2007)

[Kaltenborn et al.(2017)]Kaltenborn2017 
Kaltenborn, M.A.R., Bastian, N.U.F., Blaschke, D.B. 2017, , 96, 056024


[Love(1911)]Love1911 
Love, A. E. H., 1911, "Some Problems of Geodynamics", (Cambridge University Press, Cambridge)

[Maselli et al.(2013)]Maselli2013 
Maselli, A., Cardoso, V., Ferrari, V., Gualtieri, L., & Pani, P. 2013, , 88, 023007 

[Mellinger et al.(2017)]Mellinger2017 
Mellinger, R. D., Weber, F., Spinella, W., Contrera, G. A. & Orsaria, M. G. 2017, Universe, 3, 5

[Most et al.(2018)]Most2018 
Most, E. R., Weih, L. R., Rezzolla, L., & Schaffner-Bielich, J. 2018, arXiv:1803.00549 

[Nandi & Char(2018)]NandiC2018 
Nandi, R., & Char, P. 2018, , 857, 12 

[Paschalidis et al.(2018)]Paschalidis2018 
Paschalidis, V., Yagi, K., Alvarez-Castillo, D., Blaschke, D. B., 
& Sedrakian, A. 2018, , 97, 084038 

[Press et al.(1992)]Press1992 
Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, 
Cambridge: University Press, 2nd ed.

[Raithel et al.(2018)]Raithel2018 
Raithel, C., Özel, F., & Psaltis, D. 2018, arXiv:1803.07687

[Ranea et al.(2016)]Ranea2016 
Ranea-Sandoval, I. F., Han, S., Orsaria, M. G., Contrera, G. A., Weber, F., & Alford, M. G.  2016, Phys. Rev. C 93, 045812

[Tolman(1939)]Tolman1939 
Tolman, R. C. 1939, Physical Review, 55, 364

[Tooper(1965)]Tooper1965
Tooper, R. F. 1965, , 142, 1541 

[Typel & Blaschke(2018)]Typel2018 
Typel, S., & Blaschke, D.,2018, Universe 4, 32


[Oppenheimer & Volkoff(1939)]OppenheimerV1939 
Oppenheimer, J. R., & Volkoff, G. M. 1939, Physical Review, 55, 374

[Van Oeveren & Friedman(2017)]VanOeverenF2017 
Van Oeveren, E. D., & Friedman, J. L. 2017, , 95, 083014

[Yagi & Yunes(2017)]YagiY2017 
Yagi, K., & Yunes, N. 2017, , 681, 1




[Zdunik(2000)]Zdunik2000 
Zdunik, J. L. 2000, , 359, 311 



