Two novel approaches to the hadron-quark mixed phase in compact stars

First-order phase transitions, like the liquid-gas transition, proceed via formation of structures such as bubbles and droplets. In strongly interacting compact star matter, at the crust-core transition, but also at the hadron-quark transition in the core, these structures form different shapes dubbed"pasta phases". We describe two methods to obtain one-parameter families of hybrid equations of state (EoS) which mimic the thermodynamic behavior of pasta phases in between a low-density hadron and a high-density quark matter phase, thus generalizing the Maxwell construction. The first method replaces the behavior of pressure vs. chemical potential in a finite region around the critical %chemical potential pressure of the Maxwell construction by a polynomial interpolation. The second method uses extrapolations of the hadronic and quark matter EoS beyond the Maxwell point to define a mixing of both with weight functions bounded by finite limits around the Maxwell point. We apply both methods to the case of a hybrid EoS with a strong first order transition that entails the formation of a third family of compact stars and the corresponding mass twin phenomenon. We investigate for both models the robustness of this phenomenon against variation of the single parameter, the pressure increment at the critical chemical potential which quantifies the deviation from the Maxwell construction. We also show sets of results for other compact star observables than mass and radius, namely the moment of inertia and the baryon mass.


Introduction
The understanding of the properties of dense matter in compact star interiors is a subject of current research. Recently, great progress in this direction has been achieved by the detection of the gravitational radiation which emerged from the inspiral phase of two coalescing compact stars, an event named GW170817 [1]. Since it was observed also in all other bands of the electromagnetic spectrum, it marked the birth of multi-messenger astronomy. Among the various obtained results, GW170817 has shed light on the properties of the equation of state (EoS) of compact star matter, namely on its stiffness, since through the constraints on the tidal deformability parameter Λ [2] from the LIGO-Virgo Collaboration (LVC) results one could estimate the maximum radius of a 1.4 M compact star to R 1.4,max = 13.6 km [3] and maximum mass of nonrotating compacts stars M TOV,max = 2.16 M [4]. Of great scientific interest is the phenomenon of a phase transition from hadronic matter to a deconfined quark phase in hybrid compact stars. Those stars are comprised of a deconfined quark matter core surrounded by a hadronic mantle. The nature of the deconfinement transition is a matter of debate. Whether it exhibits a jump in the thermodynamic variables or represents a crossover 1 is a question that is addressed to both, laboratory experiments as well as compact star observations. The possible mixed phase in a neutron star presents geometrical structures of different shapes, in the literature denoted as "pasta phases". The adequate description of the physics involved is a complicated problem where the geometrical properties of the structures as well as transitions between them must be taken into account (different methodologies can be found in [5][6][7][8][9][10]12]). In the case of the hadron-quark interface, the procedure is well explained in [11] (see also [12] for a recent work): one models several geometrical structures and finds the energetically most favorable ones in different density regions inside compact stars.
In this work we take a different route and introduce two types of phenomenological interpolations which aim at mimicking the effect of those geometrical structures on the thermodynmical behavior and simultaneously explore the whole range of densities in a unified way. A first realization the idea to describe the transition from the hadronic to the quark matter phase of matter in neutron stars by an interpolation in order to model a crossover-like behaviour was carried out in [13] and followed up in Refs. [14,15], where the jump of the EoS ε(P) was replaced by a smooth behavior using as an ansatz a tangens hyperbolicus function.

Baryonic chemical potential
hadronic phase quark phase mixed phase Figure 1. Schematic representation of the interpolation function P M (µ) obtained from the mixed phase constructions discussed in this work. For both interpolation methods discussed in the text it has to go though three points: P H (µ H ), P c + ∆P and P Q (µ Q ).
A systematic and thermodynamically consistent formulation was recently given in [16,17] where a parabolic interpolation function was introduced to replace the behavior of the hybrid EoS for a Maxwell transition. We shall denote this procedure as replacement interpolation method (RIM). The resulting hybrid EoS was then used to study the effect of the mixed phase on the properties of compact stars. A second realization of this concept has been worked out recently in [18], where instead of replacing the hadronic and quark matter branches of the hybrid EoS in the limits µ H < µ < µ Q (see Fig. 1) a mixing of these branches is defined using switch functions and a bell-shaped function for the pressure 1 Here the term "crossover" is used generically for a transition that does not proceed like in a Maxwell construction at a strictly constant pressure with a jump in (energy) density, but rather by a varying pressure in the transition region. It can thus be a generic crossover transition like in ferromagnetic systens under external magnetic field, but also a first order transition for several globally conserved charges which proceeds via formation of structures of different shapes (pasta phases). increment with an amplitude ∆P = ∆ P P c , where P c = P(µ c ) is the critical pressure of the Maxwell construction. This procedure is denoted as mixing interpolation method (MIM) in [18]. The free parameter ∆ P occurs in both methods with an equivalent influence on the behaviour of the EoS in the mixed phase region, in particular on its extension, see Fig. 1. We would like to note that in both methods a negative value of ∆ P would signal that a Maxwell construction using both input EoS P H (µ) for hadronic matter and P Q (µ) for quark matter would not make sense because it would describe a transition from quark matter at low densities (where P Q (µ) is not trustworthy) to hadronic matter at high densities (where P H (µ) is not trustworthy). For a discussion of this situation, see Ref. [19].
In this work we present a comparative study of the RIM and MIM approaches to construct mixed phases of the quark-hadron phase transition which mimic the thermodynamic behaviour of pasta phases. We discuss the similarities and differences of these two approaches and apply them to obtain a hybrid EoS under neutron star constraints for which we discuss the resulting hybrid star sequences and their properties. While the first approach (RIM) is rather intuitive and simple to realise as its properties just depend on the order of interpolating polynomial, the second approach (MIM) is based on a procedure of "mixing" the EoS of the two phases in the coexistence region and reminds in its properties on the physics of substitutional compounds as in the crust of compact stars, resulting in an intermediate stiffening effect.
The paper is structured as follows. In section 2 we start with the reference EOS for the present study, for which a four-polytrope ansatz is employed which features a hadronic phase (first polytrope), a constant pressure polytrope resembling a strong first order phase transition as described by a Maxwell construction (second polytrope) and two polytropes for quark matter phases at high densities. Next, in section 3 we introduce the RIM and MIM approaches to construct mixed phases when two reference EoS for the low-density (hadronic) and high-density (quark matter) phases are given. We discuss the speed of sound c s as the key characterizing property of the family of obtained hybrid EoS. Subsequently, in section 4 we discuss the similarities and differences between the hybrid star EoS of both approaches and show results for the macroscopic properties of compact stars. We motivate these results by the feasibility of detection by multi-messenger astronomy. Consequently, future detections of gravitational wave radiation emitted by of NS-NS or NS-BH mergers shall provide new constraints on both the star mass and radius. Moreover, the determination of the fate of the merger, whether it evolves via a prompt or delayed collapse into a black hole, can be used as an independent estimate on the mass and radius, as proposed in [20]. Up to now, tests for the current compact star models with the at present still single compact star merger event have been performed, e.g., in [17,18,21].

Hybrid star EoS with a third family and high-mass twins
Compact stars are traditionally divided into white dwarf (first family) and neutron star (second family) branches. Hybrid stars whose equation of state undergoes a sufficiently strong first order phase transition (large jump in energy density ∆ε) can populate a third family branch in the mass-radius diagram, separated from the second one by a sequence of unstable configurations. As a consequence, there appear so called mass twin configurations: the second and third family solutions overlap within a certain range of masses while the radii of any two stars with the same mass (mass twins) are very different. If the mass-twin phenomenon occurs at high masses ∼ 2 M then one speaks of high-mass twin (HMT) stars [22]. Depending on the critical pressure of the phase transition, the mass-twin phenomenon can occur also at lower masses such as the typical binary radio pulsar mass of ∼ 1.35 M , see [17,18,21], so that the corresponding twin star configuration become of relevance for the interpretation of GW170817. In the latter case, a mass ratio q = m 1 /m 2 = 1 of the merger would not entail that the merging stars have the same radii and internal structure! Would the mass-twin phenomenon (at whatever mass) be observed, this would entail that the QCD phase diagram has to possess at least one critical endpoint since for the study of the cold region of the QCD phase diagram the existence of a first order phase transition between hadron to quark matter had to be concluded. Since the high temperature region of the QCD diagram is kown to feature a crossover transition, compact stars can serve as a probe of the existence of a critical end point [23] and provide insight into the properties of matter in heavy ion collision conditions [24].
In order to study the effects of pasta phases at the hadron-quark matter interface in hybrid star interiors, we consider a piecewise polytropic EoS as previously used in various works [3,[25][26][27][28]. The polytropic representation used in the present work consists of four segments of matter at densities higher than saturation density n 0 = 0.15 fm −3 (n 0 n 1 < n < n 5 ).
Each density region is labeled by i = 1 . . . 4 with prefactor κ i and polytropic index Γ i . HMT stars require a rather stiff nucleonic EoS which here is represented by the first polytrope. The hadron-quark matter first-order phase transition is described by the second polytrope with constant pressure P tr = κ 2 and vanishing polytropic index (Γ 2 = 0). At higher densities the polytropes 3 and 4 represent a rather stiff quark matter EoS. The parameters for this HMT realisation are given in table 1. For the present Table 1. Parameters for the four-polytrope EoS of Ref. [28], called "ACB4" in Ref. [21]. The corresponding description is presented in Eq. 1 of the main text. The last column displays the maximum masses M max on the hadronic (hybrid) branch corresponding to region i = 1 (i = 4). In addition, the minimal mass M min in region i = 3 of the hybrid branch is displayed in that column.
valid for the respective regions (phases) i = 1 . . . 4, where for the constant pressure region i = 2 this formula collapses to P(µ = µ c ) = P c = κ 2 because of Γ 2 = 0. For applying the MIM below, it will be important that the pressure of the hadronic phase (i = 1) valid for µ < µ c can be extrapolated to the neighboring quark matter phase (i = 3) where µ > µ c and vice-versa. HMT star EoS fulfill the Seidov conditions over quantity values at the phase transition [29] ∆ε ε c for the third family of compact stars to exist. The choice of parameters for this EoS corresponds to a sufficiently stiff high-density region in order to prevent gravitational collapse while at the same time not violating the causality condition for the speed of sound c s < c. See [28] for details.

Mixed phase constructions
In this section we present the details of the interpolation descriptions for the mixed phase between the hadronic and quark matter phases. For this purpose we consider the chemical potential dependent pressures of both the hadronic (i = 1) and the neighboring quark matter (i = 3) phases: P H (µ), P Q (µ), respectively. As mentioned above, our polytropic HMTs EoS features a first order phase transition implemented in the form of a Maxwell construction at a critical chemical potential value µ c where pressures for both phases are equal: thus both phases are in thermodynamic equilibrium.

The replacement interpolation method (RIM)
In this mixed phase approach the relevant regions of both, the hadronic and quark matter EoS around the Maxwell critical point (µ c , P c ) are replaced by a polynomial function defined as where ∆ P is a free parameter representing additional pressure of the mixed phase at µ c . Generally, the ansatz (5) for the mixed phase pressure is an even order (N = 2k, k=1,2, ...) polynomial and it smoothly matches the EoS at µ H and µ Q up to the k-th derivative of the pressure, where N + 2 parameter values (µ H , µ Q and α q , for q = 1, . . . , N) can be found by solving the above system of equations, leaving one parameter (∆P) as a free parameter of this method. The simplest case of the RIM is the parabolic model for N = 2 which has been first introduced in [16,17], As usual, the parameters α 1 , α 2 , µ H and µ Q are found from the following system of equations involving quantities at the borders of the mixed phase, It is evident that the order of the interpolating function (5) will determine whether or not there are discontinuities for the derivatives of the function P M (µ). For instance, the square of the speed of sound, involves the second derivative of the pressure with respect to µ since n = ∂P/∂µ, see figure 2. The result is that for k = 1 the function (5) exhibits a clear discontinuity in the speed of sound at ε c and ε c + ∆ε, whereas in between these borders, the speed of sound slightly increases relative to the case of the Maxwell construction for which c 2 s = 0 in the mixed phase region. For k = 2, the mixed phase pressure (5) allows for a continuous speed of sound. However, it is connected at ε c and ε c + ∆ε to the speed of sound outside these borders with a jump in its derivative. At the order k = 3 and higher the speed of sound behaves smoothly without a jump in its derivative.

The mixing interpolation method (MIM)
This approach has recently been defined in Ref. [18], where the interpolation ansatz was based on trigonometric functions. Here we will use instead a polynomial ansatz for the interpolation that  ∆ P = 0% ∆ P = 1% ∆ P = 2% ∆ P = 3% ∆ P = 4% ∆ P = 5% ∆ P = 6% ∆ P = 7% ∆ P = 8% consists of a pair of functions f off and f on that will switch off and on the hadronic and quark parts of the equation of state, as well as an additional compensating function ∆ in order to eliminate thermodynamic instabilities, see figure 3. This interpolation is applied in the pµ plane within the range µ H ≤ µ ≤ µ Q . The pressure that interpolates between the hadron and quark phase at the phase transition reads Even though f off and f on might be any switching functions, our choice of definition consists of the following pair of left and right side polynomials: that together with the complementary functions f >,R (µ) = 1 − f <,R (µ) and f <,L (µ) = 1 − f >,L (µ) will complete the switch functions. The above coefficients α L , α R , β L and β R can be determined by the following conditions where the value of 1/2 is chosen for symmetric convenience. Consequently, the switching functions are defined as and furthermore obey f off/on (µ) = 1 − f on/off (µ).
In order to construct a proper dimensionless function ∆(µ) we introduce consisting of the functions whose coefficients are determined by the conditions Regarding ∆P as the only free external parameter, up to this moment we have 10 unknowns and 8 independent equations which leave us with the possibility to fix the second order derivative of P at µ H and µ Q in the following way:

Hybrid star EoS with mixed phases
The two interpolation methods presented above result in a thermodynamically consistent EoS. Knowing that n = ∂P/∂µ, the thermodynamic identity used to derived all the needed variables at zero temperature reads ε = −p + µ n.
(17) Figure 4 shows the resulting mixed phase interpolations for both approaches charaterised by the dimensionless pressure increment ∆ P = ∆P/P c that ranges from 1% to 8%, where ∆ P = 0 reproduces the Maxwell construction. Figure 5 shows pressure values depending on energy density. The first order phase transition via a Maxwell construction corresponds to the ∆P = 0 case with the pressure being constant in the mixed phase region. Furthermore, figure 6 shows the squared speed of sound for both approaches where the difference between them becomes obvious: while the MIM shows a peak in the mixed phase region the RIM shows a rather structureless behaviour in this region. This feature is a direct consequence of the functional form of the interpolation implemented by the two methods.

Compact star sequences
In order to compute the compact star internal pressure (energy density) profiles leading to mass-radius relations, we solve the Tolman-Oppenheimer-Volkoff (TOV) equations [30,31] derived in the framework of General Relativity for a static, spherically-symmetric compact star dP(r) dr = − G (ε(r) + P(r)) M(r) + 4πr 3 P(r) r (r − 2GM(r)) , with the boundary conditions P(r = R) = 0, M(0)=0 and M(R) = M that serve to determine the total stellar mass M and total stellar radius R once a central pressure P(0) = P(r = 0) (and with it the central energy density because P(ε) is known) is given as input. By increasing the central energy density values, a whole sequence of star configurations up to the one with the maximal mass can be obtained, thus populating the mass-radius diagram. Figure 7 shows compact star sequences for all models characterised by the ∆P value for both, the MIM and RIM approaches together with up-to-date constraints from astrophysical measurements. We can notice that for the lower values of ∆P < 6% the HMT phenomenon persists regardless which mixed phase interpolation method has been applied. In Fig. 8 we show the mass versus central energy density and the radius versus central pressure for both interpolation methods. For the MIM one observes a trace of the intermediate stiffening effect in the mass versus central energy density which is absent for the RIM. In addition, two other quantities of astrophysical interest are the total baryonic mass of the star that results from integrating the following equation and similarly, its moment of intertia [32] which are related to observational phenomena as well, like energetic emissions that might conserve baryon mass or moment of inertia dependent pulsar glitches. For a detailed discussion of the moment of inertia in the slow-rotation approximation, and for the hybrid star case see, e.g., [33][34][35], and references therein. In figure 9 we show the baryon mass versus radius and and the moment of inertia versus gravitational mass for the compact star sequences obtained in this work with both interpolation methods. When increasing the pressure increment from ∆ P = 0 to 8%, the sharp edges which are obtained for the Maxwell construction case get washed out. One observes no qualitative difference between the MIM and the RIM in the patterns of these families of sequences. For ∆ P > 5%, the second and third family branches in the M B vs. R diagrams get joined so that neutron star and hybrid star configurations form a connected sequence and the HMT phenomenon get lost. This effect is reflected in the I vs. M diagrams by the loss of multi-valuedness (the lowest branch up to the maximum mass of 2.11 M shall be ignored because it is unstable). From the M B vs. R diagrams one can read off which configuration on the hybrid star branch would be reached when the maximum mass neutron star configuration would collapse under conservation of baryon number. Comparing the gravitational masses of these two star configurations one can estimate the release of binding energy in this process, see Ref. [35].   [36] and PSR J0348+432 [37], respectively. The gray and orange bands denoted by M1 and M2 are the compact star mass windows for the binary merger GW170817. The green band corresponds to the 1.44 ± 0.07 M mass of PSR J0437-4715 whose radius is expected to be measured by NICER [38]. The hatched regions are excluded by GW170817: the star radius at 1.6 M cannot be smaller than 10.68 km [20] and for a 1.4 M the star has to have a radius smaller than 13.6km [3]. The maximum mass of compact stars shall be below 2.16 M [4].

Conclusions
In this work we have introduced two interpolation approaches to a mixed phase at the hadron-quark phase transition. An advantage of these two interpolation methods presented here over the construction employing hyperbolic tanges functions [14,15]   potentials of the mixed phase between the hadronic and the quark EoS, whereas the latter strictly converges only at infinity.
While each approach uses a different functional form both of them fulfill the same conditions at the border of the mixed phase. We have found that both methods can be distinguished by the behaviour of the speed of sound that they predict. The MIM approach motivated by the analogy with sequential phase transitions occurring for substitutional compounds in the neutron star crust finds an intermediate stiffening of the mixed phase EoS. The RIM approach does not exhibit this feature. In the case of the RIM approach, we have studied both, a fourth and sixth order polynomial interpolation. We found that the latter connects the hadron and quark EoS smoothly up to second derivatives which is being visible in smooth behaviour of the speed of sound. However, the differences in the neutron star properties for both polynomial orders are safely negligible.
The macroscopic properties of compact stars show for both mixed phase constructions a very similar systematic behaviour as the pressure increment ∆P is increased: the mass-radius relation smoothens out eliminating the gap between second and third branches, but only for the highest values we have considered. Up to ∆ P ∼ 5% the HMT phenomenon is robust against the mixed phase construction, regardless whether the MIM or RIM approach is used. For the mass versus central energy density, one observes a trace of the intermediate stiffening effect for the MIM which is absent for the RIM. For the other compact star quantities evaluated here, the baryonic mass and the moment of inertia, both interpolation methods display a similar type of behaviour when the pressure increment is varied.
The methods presented here can potentially be applied to the compact star crust-core transition as well. Just like at the hadron-quark boundary, the transition at the bottom of the crust may proceed via pasta phases dominated by Coulomb forces and surface tension effects [39]. Further astrophysical aspects of mixed phases inside neutron stars include potentially observable effects such as the rotational evolution, pulsar glitches, gravitational wave emission and cooling. They could be sufficiently sensible to the nature of the phase transition, proceeding via pasta phases or not, and thus provide potential signatures of the presence and extension of a mixed phase in compact stars.