Systematic study on the quark-hadron mixed phase in compact stars

We investigate systematically the quark-hadron mixed phase in dense stellar matter, and its influence on compact star structures. The properties of quark matter and hadronic matter are fixed based on various model predictions. Beside adopting constant values, the surface tension $\Sigma$ for the quark-hadron interface is estimated with the multiple reflection expansion method and equivparticle model. To fix the structures of quark-hadron pasta phases, a continuous dimensionality of the structure is adopted as proposed by Ravenhall, Pethick, and Wilson. The corresponding properties of hybrid stars are then obtained and confronted with pulsar observations. It is found that the correlation between radius and tidal deformability in traditional neutron stars preserves in hybrid stars. For those permitted by pulsar observations, in almost all cases the quark phase persists inside the most massive compact stars. The quark-hadron interface plays an important role on hybrid star structures once quark matter emerges. The surface tension $\Sigma$ estimated with various methods increases with density, which predicts stiffer EOSs for the quark-hadron mixed phase and increases the maximum mass of hybrid stars. The EOSs of hybrid star matter are well constrained at densities $n\lesssim 0.8$ fm${}^{-3}$, while larger uncertainty is expected at higher densities.

We investigate systematically the quark-hadron mixed phase in dense stellar matter, and its influence on compact star structures. The properties of quark matter and hadronic matter are fixed based on various model predictions. Beside adopting constant values, the surface tension Σ for the quark-hadron interface is estimated with the multiple reflection expansion method and equivparticle model. To fix the structures of quark-hadron pasta phases, a continuous dimensionality of the structure is adopted as proposed by Ravenhall, Pethick, and Wilson. The corresponding properties of hybrid stars are then obtained and confronted with pulsar observations. It is found that the correlation between radius and tidal deformability in traditional neutron stars preserves in hybrid stars. For those permitted by pulsar observations, in almost all cases the quark phase persists inside the most massive compact stars. The quark-hadron interface plays an important role on hybrid star structures once quark matter emerges. The surface tension Σ estimated with various methods increases with density, which predicts stiffer EOSs for the quark-hadron mixed phase and increases the maximum mass of hybrid stars. The EOSs of hybrid star matter are well constrained at densities n < ∼ 0.8 fm −3 , while larger uncertainty is expected at higher densities.

I. INTRODUCTION
Due to the asymptotic freedom of strong interaction, the deconfinement phase transition is expected as one increases the density of hadronic matter. However, it is still unclear how such a transition takes place. Traditionally, for zero temperature cases, a first-order phase transition between hadronic matter (HM) and quark matter (QM) was predicted by various quark models, which indicates a quark-hadron mixed phase (MP) [1][2][3][4][5]. Adopting different surface tension values for the quark-hadron interface, it was found that the MP exhibits various structures [6]. For vanishing surface tensions and Coulomb interactions, the MP is comprised of HM and QM that satisfy the Gibbs condition [1]. If a moderate surface tension is employed, with the charged particles relocate themselves via charge screening effects, geometrical structures appear [6][7][8][9][10][11][12][13][14]. Those structures become unstable for enough large surface tensions, which leads to a bulk separation of quark and hadron phases, i.e., the Maxwell construction scenarios.
In our previous study [13], we have considered the possibility of constraining the surface tension from pulsar observations, where a first-order deconfinement phase transition was assumed. By adopting the covariant density functional TW99 [43] for nuclear matter and perturbation model [44] for quark matter, it was found that varying the surface tension value will have sizable effects on the radii and tidal deformabilities of 1.36-solarmass hybrid stars.
Nevertheless, due to the important roles played by many-body interactions as well as the emergence of hadrons other than nucleons, the properties of hadronic matter at densities larger than twice the nuclear saturation density are not very well constrained, where the differences between various predictions grow dramatically [45]. Meanwhile, even though the perturbation model gives reliable predictions at ultra-high densities [46], the properties of quark matter inside hybrid stars are poorly constrained. Under such circumstances, in the present work, we further extend our study by investigating systematically the hadron-quark deconfinement phase transition in dense stellar matter, where various combinations of models that describe QM and HM are adopted along with different values of surface tension.
To fix the structures of quark-hadron pasta phases, a continuous dimensionality of the structure is introduced as proposed by Ravenhall et al. [57]. The energy contribution due to the quark-hadron interface is treated with a surface tension Σ, for which we employ constant values as well as those estimated by the multiple reflection expansion method [58][59][60][61] and equivparticle model including both linear confinement and leadingorder perturbative interactions [62,63].
The EOSs of hybrid star matter are obtained, while the corresponding compact star structures are determined by solving the Tolman-Oppenheimer-Volkov (TOV) equation. For the EOSs of hybrid star matter consistent with pulsar observations, it is found that in almost all cases the quark phase takes place inside the most massive compact stars. Once quark matter emerges, we find that the quark-hadron interface plays an important role on hybrid star structures.
The paper is organized as follows. We present our theoretical framework in Sec. II, Sec. III, and Sec. IV. Two formalisms are adopted for the HM, i.e., the RMF model in Sec. II A and the variational method in Sec. II B. The equivparticle model, perturbation model, and NJL model for QM are introduced in Sec. III. The formalism in obtaining the structures of quark-hadron mixed phase is introduced in Sec. IV A, while the surface tension of quark-hadron interface is obtained in Sec. IV B. The numerical results are presented and discussed in Sec. V. Our conclusion is given in Sec. VI.

II. EFFECTIVE MODELS FOR HADRONIC MATTER
A. RMF model The Lagrangian density for infinite nuclear matter obtained with RMF model [47] reads Here the Dirac spinor Ψ i represents nucleons with the effective mass m * = m+g σ σ and isospin τ i . Three types of mesons are included to describe the interactions between nucleons, i.e., σ-, ω-, and ρ-mesons with their masses being m σ , m ω and m ρ , respectively. The baryon number density is given by n = n n + n p = i=n,p Ψ i γ 0 Ψ i . In this work, we adopt two different schemes for the density dependence of effective interaction strengths, i.e., the nonlinear self-couplings of σ and ω mesons in U (σ, ω) and the Typel-Wolter ansatz with density dependent coupling constants [43].
The nonlinear self-couplings for σ and ω mesons are where we have adopted the effective interaction TM1 [64], i.e., Shen EOS2 [65]. Meanwhile, it was shown that the slope of symmetry energy L = 110.8 MeV predicted by TM1 was too large according to various constrains from nuclear physics and pulsar observations, which can be reduced to L = 40 MeV by adding the cross coupling term This gives Shen EOS4 by adopting the effective interaction TM1e [66]. To further include the contribution of Λ hyperons, in Eq. (1) we add the following Lagrangian density [67][68][69][70], where m * Λ = m Λ + α σΛ g σ σ is the effective mass of the Λ hyperon. The ratio of coupling constants α ωΛ ≡ g ωΛ /g ω = 2/3 is predicted by the naive quark model [71], then α σΛ ≡ g σΛ /g σ = 0.621 is obtained by reproducing the binding energies of Λ-hyperon in Λ-hypernuclei, i.e., Shen EOS3 [65]. However, the obtained hyperonic EOS is too soft to support massive neutron stars, i.e., the hyperon puzzle. To resolve this, larger values of g ωΛ were adopted to provide more repulsive interaction from ω meson, where in this work (denoted as TM1Λ) we take α ωΛ = 1 and α σΛ = 0.887 [70].
Despite the great successes in describing finite nuclei with nonlinear self-couplings of mesons, a direct extension of the density functional to higher densities may cause problems of stability. Alternatively, we can adopt couplings that depend explicitly on densities, which can be derived from self-energies of Dirac-Brueckner calculations of nuclear matter [43,72]. We thus adopt the effective nucleon-nucleon interactions PKDD [73], TW99 [43], DDME2 [74], and DD2 [75], where U (σ, ω) = 0 and the density dependence of coupling constants g σ,ω,ρ [43] are obtained with Here n 0 represents the nuclear saturation density. Carrying out standard mean-field and no-sea approximations, one obtains the energy density E, chemical potentials µ B , and pressure P at given baryon density n. Then the EOSs for nuclear matter and hypronic matter can be obtained.

B. Variational methods
The variational method for uniform nuclear matter was developed in Refs. [76][77][78], where the nuclear Hamiltonian composed of a two-body potential V ij and threebody potentials V ijk are given by Adopting the Argonne v18 (AV18) two-body nuclear potential [79] and the Urbana IX (UIX) three-body nuclear force [80,81], the free energy per nucleon of uniform nuclear matter is predicted by the cluster variational method using the Jastrow wave function [49]. Then the equation of states for nuclear matter (denoted as VM) can be obtained, which was discussed in detail in Ref. [49]. For hyperonic EOS (VMΛ), we adopt the results presented in Ref. [82] with three body forces of hyperons. A more sophisticated variational method with the Fermi Hypernetted Chain calculations was performed for symmetric nuclear matter (SNM) and pure neutron matter (PNM) by Akmal, Pandharipande, and Ravenhall (APR) [48], where the aforementioned realistic nuclear Hamiltonian and Jastrow wave function were adopted. The energy density of nuclear matter obtained in Ref. [48] are fixed by the fitted formula Here δ = (n n −n p )/n is the isospin asymmetry and ν p,n = (3π 2 n p,n ) 1/3 the Fermi momentum of nucleons. A more detailed description on the parameters p i and functional form g(n, δ) can be found in the original publication [48].
In this work, we adopt the most comprehensive case employing the AV18 two-body nuclear potential and UIX three-body interaction [80,81] with relativistic corrections.
C. The EOSs of nuclear/hyperonic matter Finally, for the hadronic phase, we adopt in total 10 different EOSs, i.e., 8 nuclear EOSs (TM1e, TM1, PKDD, TW99, DDME2, DD2, VM, APR) and 2 hyperonic EOSs (TM1Λ and VMΛ). These EOSs are predicted by both RMF model with various effective interactions and variational methods started from realistic baryon interactions. The corresponding saturation properties are indicated in Table I and compared with the constraints from terrestrial experiments and nuclear theories [83], which give the binding energy ε ≈ 16 MeV, the incompressibility K = 240 ± 20 MeV [84], the symmetry energy S = 31.7 ± 3.2 MeV and its slope L = 58.7±28.1 MeV [85,86] around n 0 ≈ 0.15-0.16 fm −3 and δ = 0. The uncertainties may be further reduced if the constraints from the GW170817 binary neutron star merger event [36,42] are included [87], e.g., a recent estimation suggests K = 250.23 ± 20.16 MeV, S = 31.35 ± 2.08 MeV and L = 59.57 ± 10.06 MeV [88]. In general, TM1 and PKDD slightly overestimate K, S, and L, while TM1e predicts reasonable symmetry energy properties. The VM EOS has the smallest S and L but still lie within the permitted ranges.
At larger densities, in Fig. 1 we present the pressures of SNM and PNM as functions of baryon number density, which are compared with the constraints from the flow data of heavy ion collisions [89]. It is found that the EOSs of nuclear matter predicted by TM1e, TM1, PKDD, DDME2 and DD2 are slightly stiffer than those constrained from the flow data of heavy ion collisions [89]. Nevertheless, the emergence of the quark phase may ease the tension and reduce the stiffness of EOSs effectively.
The EOS of neutron star matter can be obtained by further including the contributions of electrons and muons, where their energy densities take the form of free Fermi gas with  The pressures of symmetric nuclear matter (SNM) and pure neutron matter (PNM) predicted by various nuclear theories, which are compared with the experimental constraints from the flow data [89].
Here g e,µ = 2 is the degeneracy factor and x e,µ ≡ ν e,µ /m e,µ with ν e,µ being the Fermi momentum of leptons, which predicts their number densities n e,µ = ν 3 e,µ /3π 2 . The total energy density of neutron star matter is obtained with E = E HM + E e + E µ . Then the pressure is determined by P = i µ i n i − E with the chemical potential µ i = ∂E ∂ni . In Fig. 2 we present the EOSs of neutron star matter, which are obtained by simultaneously fulfilling the β-stability condition and local charge neutrality condition.
Based on the EOSs indicated in Fig. 2, the correspond- ing structures of compact stars are obtained by solving the TOV equation while the tidal deformability is estimated with Here the gravity constant is taken as G = 6.707 × 10 −45 MeV −2 , while k 2 is the second Love number and is obtained from the response of the induced quadrupole moment Q ij in a static external quadrupolar tidal field [90][91][92]. Note that at n < 0.08 fm −3 we have adopted the EOSs presented in Refs. [93][94][95], which account for the crusts of neutron stars. For the cases of TM1, TM1Λ, TM1e, VM, and VMΛ, the crust EOSs were previously obtained, i.e., Shen EOSs [65,66] and VM EOSs [49,82]. However, instead of using those EOSs, we still adopt the crust EOSs presented in Refs. [93][94][95] since the variations on neutron star structures are relatively small. The obtained mass, radius, and tidal deformability are presented in Fig. 3 and compared with astrophysical observations, where the maximum masses and radii of 1.4M neutron stars are indicated in Table I. All the maximum masses of compact stars predicted by various EOSs in Fig. 2 are consistent with the observational mass (2.14 +0. 10 −0.09 M ) of PSR J0740+6620 [28]. Nevertheless, we should mention that the velocity of sound will exceed c for APR, VM, and VMΛ at n ≥ 0.87, 0.90, and 1.08 fm −3 , which are reached in center regions of the massive compact stars indicated in Fig. 3. The tidal deformabilities predicted by TM1e,  Table I. PKDD, DDME2, and DD2 slightly exceed the constraint 302 ≤Λ ≤ 720 from the GW170817 binary neutron star merger event [36][37][38][39][40][41], which coincide with the experimental constraints on SNM from the flow data [89] in Fig. 1. Note that the recent radius measurements of PSR J0030+0451 with the equatorial radius R eq = 11.52-14.26 km and mass M = 1.18-1.59 M obtained via pulse-profile modeling in the NICER mission [34,35] do not constrain the EOSs adopted here.

A. Equivparticle model
In the equivparticle model, the quarks are treated as quasi-free particles with density dependent equivalent masses. Taking into account both the linear confinement and leading order perturbative interactions, the quark mass scaling is given by [52] Here m i0 is the current mass of quark flavor i (i = u, d, s) [96] and n ≡ i=u,d,s n i /3 the baryon number density. The confinement parameter D is connected to the string tension σ 0 , the chiral restoration density ρ * , and the sum of the vacuum chiral condensates q qq 0 . Meanwhile, the perturbative strength parameter C is linked to the strong coupling constant α s . Due to the uncertainties in relevant quantities, we do not know the exact values of D and C. Nevertheless, it has been estimated that √ D approximately lies in the range of (147, 270) MeV [51] and C < ∼ 1.2 [52]. In this work, we At zero temperature, the energy density E QM = i=u,d,s E 0 i (ν i , m i ) and particle number density n i = ν 3 i /π 2 are identical to the cases of free Fermi gas with E 0 i given by Eq. (9) and g u,d,s = 6. Note that in Eq. (9) we have adopted the mass scaling of Eq. (13) for quarks, i.e., m i ≡ m i (n). The pressure is determined by

C. NJL model with vector interactions
In the mean-field approximation, the Lagrangian density of a SU(3) NJL model is given by where the constituent quark mass reads At T = 0, the thermodynamic potential density of quark matter predicted by the NJL model is determined by with ω 0 i given by Eq. (15) and E 0 i by Eq. (9). Here a constant E 0 is introduced to ensure Ω QM = 0 in the vacuum. Λ is the three dimensional momentum cutoff to regularize the vacuum part, and µ * i the effective chemical potential which is connected with the true chemical potential via Based on the thermodynamic potential density in Eq. (22), the chiral condensate is given by σ i = ∂ΩNJL ∂Mi and quark number density n i = ν 3 i /π 2 . The energy density and pressure are then obtained with E QM = Ω QM + i µ i n i and P QM = −Ω QM . In this work, two different sets of parameters are adopted, i.e., the sets HK (Λ = 631 For the vector coupling G V , the Fierz-transition predicts G V = 0.5G S , while in this work we take it as a free parameter with G V = 0, 0.5G S , G S , and 1.5G S .

D. General discussion on the quark EOSs
In contrast to nuclear matter cases, we have little constraints on the properties of quark matter at intermediate densities. At ultra-high densities (n > ∼ 40n 0 ), however, Quantum Chromodynamics (QCD) can be solved with perturbative approaches [46]. The corresponding EOS at highest densities is then expected to be reproduced by the perturbation model (pQCD) explained in Sec. III B. The NJL model, on the other hand, is a low-energy model for QCD, where the gluons are integrated out while retaining only local quark interactions.
The corresponding coupling constants of NJL model are then fixed by reproducing the masses of π, K, η and the π decay constant [55,56]. The equivparticle model carries similar traits of quasiparticle model [100,101], where the results of pQCD at highest densities can be reproduced with the parameter C in Eq. (13) depending explicitly on α s [52]. Meanwhile, the linear confinement of quarks are well treated with an inversely cubic mass scaling in the equivparticle model [102]. Note that in this work we have neglected the effects of color superconductivity [103], which shall be considered in our future studies.
With the energy contributions of leptons determined by Eq. (9), the EOSs of quark matter in compact stars can be obtained by simultaneously fulfilling the βstability condition µ u + µ e = µ d = µ s , µ e = µ µ and local charge neutrality condition i q i n i = 0 with q i (q u = 2/3, q d = q s = −1/3 and q e = q µ = −1) being the charge of particle type i. The corresponding energy per baryon of quark matter predicted by various quark models are then presented in Fig. 4, which include 46 EOSs of quark matter, i.e., 6 of them obtained with equivparticle model, 32 with perturbation model, and 8 with NJL model. Note that stiffer EOSs are obtained with larger C, C 1 , G V and smaller ∆µ in those models.

A. Inhomogeneous structures
If the surface tension of quark-hadron interface Σ is smaller than the critical value Σ c , inhomogeneous structures of the mixed phase will emerge, i.e., the pasta phases. Adopting linearization for the charge densities of quark and hadron phases, the critical surface tension Σ c can be estimated with [10] where µ H,Q e0 are the electron chemical potential and λ H,Q D the Debye screening length of hadronic matter (H) and quark matter (Q) fulfilling both the β-stability condition and local charge neutrality condition. To obtain the properties of quark-hadron pasta phases, we adopt the formalism with a continuous dimensionality proposed by Ravenhall et al. [57], where the energy density is determined by with Here χ, r, and n I ch are the volume fraction, radius, and charge density of phase I, and E I,II the corresponding energy densities. The continuous dimensionality d lies in the range 1 ≤ d ≤ 3, where d = 1, 2, 3 represent the slab, rod/tube, droplet/bubble phases, respectively. For the case d = 2, Eq. (27) could yield the correct expression containing a logarithmic term [104]. The global charge neutrality condition χn I ch + (1 − χ)n II ch = 0 is fulfilled for the two phases in the cell. In Eq. (25), the term E s represents the energy contribution of the quarkhadron interface, while E C corresponds to the Coulomb energy per unit volume. Note that due to the charge screening effects, the local densities n I,II , n I,II ch , and E I,II should in principal vary with space coordinates. For simplicity, we neglect such effects and the densities in each phase are assumed to be constants, which may affect our estimations on the sizes of the inhomogeneous structures at large Σ (close to Σ c ). However, for smaller Σ, the nonuniform distributions of charged particles in each phase become insignificant and assuming constant densities gives a fairly well description. In any cases, the negligence of charge screening effects has little impact on the obtained EOSs of MP.
The structures of MP can be fixed by minimizing the energy density in Eq. (25) at a given total baryon number density n = χn I + (1 − χ)n II . By taking derivatives of E t with respect to each independent parameters (r, χ, d, n I , n I ch ) and equate them to zero, one obtains the following equations: Then the structures of MP are obtained by simultaneously fulfilling those equations, while the exact phase state (I = H or Q) is determined for the case that gives a smaller E t . The quark fraction χ Q is then fixed by In practice, to further simplify our calculation, we expand the thermodynamic quantities with respect to µ e at a given baryon chemical potential µ B as was done in Ref. [13], where the chemical potential of each particle species is given by Here B i (B p = B n = 1, B u = B d = B s = 1/3, and B e = B µ = 0) is the baryon number and q i (q p = 1, q n = 0, q u = 2/3, q d = q s = −1/3 and q e = q µ = −1) the charge of particle type i. At given µ B and µ e , the pressure and energy densities are obtained with Here P 0 , E 0 , and µ e0 are the pressure, energy density, and electron chemical potential corresponding to those in Figs. 2 and 4, while the derivatives n ch = ∂n ch ∂µe , E = ∂E ∂µe , E = ∂ 2 E ∂µ 2 e are taken at µ e = µ e0 . The Debye screening length is related to n ch with λ D ≡ (−4παn ch ) −1/2 . According to the basic thermodynamic relations, the charge density and baryon number density are obtained with n ch = n ch (µ e − µ e0 ) and n = (E + µ e n ch + P )/µ B .
By counting the number of depleted quarks, the average effects due to quark depletion are treated with a modification to the density of states, i.e., the multiple reflection expansion (MRE) method [58][59][60][61]. Consider only the surface term, the modification for each quark flavor i (i = u, d, s) reads [58] dN MRE where N MRE i is the negative number of depleted quarks, p i the momentum of quarks, and S the surface area. The corresponding contribution to the surface tension for each quark species i is then obtained by equating the pressure Here x i ≡ ν i /m i with ν i being the Fermi momentum of quarks. The surface tension predicted by the MRE method is then obtained with Σ = Σ MRE = i=u,d,s Σ MRE i . Nevertheless, the MRE method tends to overestimate the surface tension by twice the value obtained with equivparticle model [62,63]. This is mainly due to the different confinement potential adopted in those models, where the bag mechanism of the MRE method introduces an infinite wall that results in a sharp density discontinuity. Since the potential between quarks is proportional to the distance instead of a wall [119], a more realistic scenario was obtained with equivparticle model where confinement can be reached with density dependent quark masses in Eq. (13). A smoothly varying quark density is then obtained on the interface, where the surface tension was found to be connected with the density of quark matter by Σ ≈ 14.3n Q + 1.3 (in MeV/fm 2 ) [63]. Meanwhile, it was shown that the surface tension predicted by the MRE method coincides with equivparticle model if we introduce a dampening factor, i.e., Σ = 0.3 i=u,d,s Σ MRE i . Note that the surface tension obtained by the equivparticle model [62,63] was for the quark-vacuum interface. The contributions from the hadron phase can be roughly included by replacing the density of quark matter n Q by the density difference ∆n ≡ |n Q − n H | between the two phases. To avoid complications in minimizing the energy density in Eq. (25), we fix Σ at a given total baryon number density n for all cases considered here. Nevertheless, the surface tension estimated with the MRE method or equivparticle model varies with density, which will alter the baryon chemical potential and pressure with The corresponding EOSs of MP will thus become stiffer if dΣ dn > 0, which is the case in our current study. Due to the uncertainties in Σ, in this work we adopt 9 different values, i.e., • Σ = 0 with Gibbs construction; • Σ = 5, 20, 50 MeV/fm 2 ; • Σ = 0.5Σ c with Σ c predicted by Eq. (24); • Σ = Σ MRE and Σ = 0.3Σ MRE ; • Σ = 14.3∆n + 1.3; • Σ > Σ c with Maxwell construction.
We have adopted both the Gibbs and Maxwell constructions in the two extreme scenarios with Σ = 0 and Σ > Σ c , a detailed description on those phase construction schemes can be found in our previous publication [13].

V. RESULTS AND DISCUSSIONS
By equating the pressures of hadronic matter in Fig. 2 and quark matter in Fig. 4, we obtain the critical chemical potential µ T B at which deconfinement phase transition occurs. In Fig. 5 7, 180). In our previous study [13], we have found a linear correlation Σ c = 0.23(µ T B − 930) + 19 with Σ c in MeV/fm 2 and µ T B in MeV, which is indicated in Fig. 5 with a black solid line. However, such a linear correlation fails to reproduce most of the current results in Fig. 5. In particular, we notice that the slope and intercept of the line vary with the adopted EOSs for both HM and QM. The critical surface tension Σc estimated with Eq. (24) as a function of the chemical potential µ T B on the occurrence of deconfinement phase transition. The symbol color indicates the hadronic EOSs adopted, while the size and shape represent the adopted quark model. The permitted maximum masses of hybrid stars predicted by various combinations of hadronic matter EOSs and quark matter EOSs. For NJL model, the parameter set (HK [55] or RKH [56], GV /GS) is indicated explicitly. Similarly, the parameter set (C, √ D in MeV) for equivparticle model is presented as well.
where the corresponding structures of hybrid stars are determined by solving the TOV equation (10). Meanwhile, the tidal deformabilities of those stars are estimated with Eq. (12). In Fig. 6 we present the obtained tidal deformability (Λ 1.4 ) as a function of radius (R 1.4 ) for 1.4M compact stars, which shows a correlation between those two observable. In general, the traditional neutron stars indicated with open squares have the largest radius and tidal deformability, which will decrease as we include the quark phase. In most cases, for a given hadronic EOS, there are stronger correlations between Λ 1.4 and R 1.4 in hybrid stars, while the inclusion of hyperons has little impact on the relation. For traditional neutron stars, in Fig. 6 the black curve indicates the relation Λ 1.4 = 5.9 × 10 −5 R 6. 26 1.4 obtained in Ref. [120], while the correlation between the maximum mass M max and Λ 1.4 found in Ref. [88] is not observed due to the first-order deconfinement phase transition.
By comparing with the observational mass (2.14 +0.10 −0.09 M ) of PSR J0740+6620 [28] as well as the tidal deformability constraint 70 ≤ Λ 1.4 ≤ 580 from the GW170817 binary neutron star merger event [42], we obtain the permitted combinations of hadronic and quark EOSs along with different values of surface tensions. In addition, we require the hadron-quark transition density n T at µ T B is larger than 0.2 fm −3 according to heavy-ion collision phenomenology, while the sound speed of hybrid star matter dose not exceed c.  Fig. 7 we present the maximum masses of hybrid stars that are consistent with the aforementioned constraints. For the choices of hadronic EOSs, only APR, VM, VMΛ, and TW99 persist due to the constraint of the tidal deformability as indicated in Fig. 3. For APR, VM, and VMΛ, the maximum mass is reduced since we require the the sound speed v < c. In most of the cases, the obtained M max increases with Σ, while a slight deviation is observed for a few cases if Σ is not constant, e.g., Σ = Σ MRE or Σ = 14.3∆n + 1.3. As indicated in Fig. 11, the obtained Σ increases with density, so that the corresponding EOSs of hybrid star matter are stiffer. In such cases, the maximum mass of hybrid stars may be even larger than those obtained with Maxwell construction at Σ > Σ c , e.g., the combinations of VM and equivparticle model. A comparison between the values of M max obtained with VM and VMΛ shows that the structures of the most massive compact stars are altered by the emergence of hyperons, despite the deconfinement phase transition occurred in the center of a compact star. This can be easily identified according to the hyperon and quark fractions indicated in Fig. 9, where the Λhyperon appears before QM at n ≈ 0. much larger than HM at n < ∼ 0.8 fm −3 , so that χ Q is reduced and the onset densities of QM are larger than hyperons.
For those where the quark matter properties are obtained with perturbation model, the corresponding hadron-quark transition density n T is small for large enough C 1 . In such cases, beside the most massive ones, the structures of 1.4M hybrid stars are also altered by the quark phase. This can be observed in Fig. 8, where the maximum masses and tidal deformabilities of 1.4M hybrid stars are presented. At C 1 > ∼ 2.5, the obtained maximum mass M max and tidal deformability Λ 1.4 usually increase with Σ, while there are few exceptions, e.g., PKDD, DDME2, and DD2. Note that in our previous study using TW99, sizable variation on the tidal deformability was observed with respect to Σ [13], which is not evident in Fig. 8 since we have ruled out the cases with M max < 2.05 M . If the surface tension Σ increases with density, the tidal deformability and maximum mass of hybrid stars may become larger. The obtained EOSs of hybrid star matter are softer for larger ∆µ, which predicts smaller M max . Nevertheless, it is worth mentioning that in most cases increasing Σ will increase M max more evidently in comparison with decreasing ∆µ. For the cases of VM and VMΛ, adopting small C 1 will give similar conclusion as in Fig. 7 since the corresponding quark matter is too unstable to completely exclude hyperons, where the maximum mass is reduced with the emergence of hyperons. However, if we take C 1 = 3, as indicated in Fig. 9, quark matter becomes more stable so that hyperons are suppressed with χ Λ < ∼ 0.06 due to a deconfinement phase transition. The structures of hybrid stars are thus hardly affected by hyperons. In any cases, hyperons do not appear in 1.4M compact stars with the central density n central < ∼ 0.57 fm −3 so that Λ 1.4 are the same for those obtained with VM and VMΛ.
With the permitted combinations of hadronic and quark EOSs along with different values of surface tensions indicated in Figs. 7 and 8, in Fig. 10 we present the corresponding radius, dimensionality, and quark fraction of MP in hybrid star matter as functions of baryon number density. After the quark matter (phase I) appears and forms the droplet phase, the obtained dimensionality will later decrease from d = 3 to d = 1 as we increase the density. Then the phases I and II switch and the dimensionality increases from d = 1 to 3. For the cases with dimensionality lies in between (1 < d < 3), the radius r varies smoothly, while sudden variations are observed during the transition to d = 1 or 3. It is interesting to note that the structures of quark-hadron pasta phases vary smoothly by treating the dimensionality as a continuous variable with d = 1-3, which may resemble the evolution of intermediate structures of droplet and rod, slab and tube found in the quantum molecular dynamics (QMD) simulations [121] as well as in the fully three-dimensional calculations adopting RMF model and Thomas-Fermi approximation [122]. An early emergence of quark matter is observed if we adopt DD2, DDME2, PKDD, and TM1e for HM and perturbation model for QM, which is necessary in order to reduce the tidal deformabilities of 1.4M hybrid stars that was otherwise too large for traditional neutron stars. In general, the quark fraction χ Q increases with density, while there are several cases where χ Q varies non-monotonically if we adopt the perturbation model for QM. At certain choices of parameters, a mixed phase may even appear after the formation of quark phase with χ Q = 1, i.e., a retrograde transition with QM→MP→QM, which is mainly caused by adopting Eq. (19) for the quark phase. The transition with HM→MP→HM is also observed for few cases such as those obtained with NJL model, which is an artifact since the color confinement is not accounted for and we thus take χ Q = 0, i.e., assuming a single phase of HM. The variation of χ Q becomes more drastic if larger surface tension values were adopted and the density range of MP shrinks. Meanwhile, the obtained radii of phase I are usually on the order of fm and increase with Σ, which was discussed extensively in previous studies [6][7][8][9][10][11][12][13].
In Fig. 11 we present the surface tension values estimated with various methods, which are much smaller than the critical surface tension Σ c indicated in Fig. 5. This is a strong indication that the quark-hadron mixed phase prefers to form inhomogeneous structures inside hybrid stars, where the energy reduction Σ c S arises from the relocation of charged particles is larger than the surface energy ΣS. Note that we have fixed the surface tension value at a given total baryon number density n, so that one does not need to worry about the variations of Σ in minimizing the energy density in Eq. (25). In general, the surface tension values predicted by various methods are increasing with µ B , where the formula Σ = 14.3∆n + 1.3 gives the smallest Σ. The corresponding EOSs for MP are thus stiffer than those obtained with constant surface tension values, which increase the maximum masses of hybrid stars as indicated in Figs. 7 and 8. Note that at µ B ≈ 1400 MeV, for few cases the perturbation model predicts smaller baryon number density of quark matter than that of hadronic matter, which causes fluctuations if the surface tension is estimated with Σ = 14.3∆n+1.3. At larger µ B , the relation 0.3Σ MRE ≈ 14.3∆n + 1.3 can be fulfilled approximately, which coincide with the predictions of equivparticle model [62,63]. With the surface tension values indicated in Fig. 11, the obtained structures of quark-hadron mixed phase follow the same trend as indicated in Fig. 10, where the radius r increases with Σ.
The EOSs of hybrid star matter are presented in Fig. 12, where the variation of pressure is insignificant at n < ∼ 0.8 fm −3 despite the emergence of quark matter in many cases. At larger densities, the uncertainty grows and the pressure obtained with NJL model as well as few cases of equivparticle model are much larger than those of perturbation model. Nevertheless, we mention that the pressure obtained with perturbation model is expected to be reasonable at highest densities since a perturbation expansion with respect to α s becomes more reliable [44]. The EOSs of quark-hadron mixed phase is sensitive to the surface tension value Σ, which become softer with smaller density range if Σ is larger. The corresponding structures of hybrid stars are thus varying with Σ as well, where the maximum mass and tidal deformation are altered as indicated in Figs. 7 and 8.
Based on the EOSs presented in Fig. 12, the structures of compact stars can be obtained by solving the TOV equation (10). In the left panel of Fig. 13 we present the corresponding M -R relation of hybrid stars. The obtained radii of 1.4-solar-mass compact stars range from 11.1 km to 12.9 km, which are consistent with the radius measurements of J0030+0451 [34,35] and binary neutron star merger event GW170817 [42]. Meanwhile, the total mass of quark matter and quark-hadron mixed phase is obtained with M QM + M MP = Rc 0 4πE t r 2 dr, where R c is the critical radius that at r > R c the quark fraction χ Q reduces to 0. The corresponding fraction (M QM + M MP )/M is then indicated in the right panel of Fig. 13. The obtained M -R relation is identical to those in Fig. 3 before a deconfinement phase transition in the center. Once the quark phase emerges, the radius of a compact star becomes smaller and the compactness increases. For those with a smaller onset densities of quark phase (e.g., TM1e/PKDD/DDME2/DD2 & perturbation model), the fraction (M QM + M MP )/M increases quickly starting at M ≈ 0.2-1 M and approaches to almost 1 at M = M max . A third family of compact stars [123,124] is observed for the case with a combination of PKDD & perturbation model (C 1 = 3, ∆µ = 770 MeV) & Σ > Σ c , where a jump of 152 MeV fm −3 in the energy density from nuclear matter to quark matter is predicted. Note that if we take a larger ∆µ for Eq. (19), the third family of compact stars indicated in Fig. 13 will not be permitted by astrophysical observations. In almost all combinations of hadronic and quark EOSs along with different values of surface tensions, the quark phase persists in the most massive compact stars, where the fraction (M QM + M MP )/M ranges from 0 to almost 1.
Aside from the observational constraints on the mass, radius, and tidal deformability, the thermal evolution of compact stars also provides important information on their internal composition [125]. Based on the thermal emission, kinematic measurements, spin period and its derivative, both the surface temperatures and ages of compact stars can be estimated [126,127]. According to various observational data, the theoretical cooling models suggest that rapid cooling due to the direct Urca (DU) processes should not occur in typical neutron stars within the mass range 1-1.5 M [128][129][130]. The Hybrid star masses as functions of the central density corresponding to those indcated in Fig. 13. The open and full circles represent the critical mass and central density beyond which the direct Urca processes of rapid cooling may take place.
DU process in nuclear matter involves the β-decay and electron capture processes of nucleons, i.e., n → p + e − +ν e and p + e − → n + ν e . The quark analogs of the nucleon DU processes are d → u + e − +ν e and u + e − → d + ν e . Those processes will occur inside compact stars once the momentum conservation is fulfilled, i.e., the triangle inequalities ν n ≤ ν p + ν e and ν d ≤ ν u + ν e with ν i being the Fermi momentum [131]. If the strangeness is involved, the DU processes such as Λ → p + e − +ν e and s → u + e − +ν e should also take effects. However, we neglect those processes here since their neutrino emissivities are expected to be less than that of nucleon/quark DU processes [125] while hyperons appear at rather large densities as indicated in Fig. 9. In Fig. 14 the hybrid star masses as functions of the central density are presented, where the open and full circles mark the critical mass (M DU ) and central density (n DU ) fulfilling the triangle inequalities. For stars with masses larger than M DU , it was shown that the neutrino emissivity is enhanced significantly by the DU processes [132], which cool the stars too rapidly within just a few years [128].
The obtained critical mass and central density corresponding to NJL model and equivparticle model are usually large with M DU > ∼ 2M and n DU > ∼ 0.8 fm −3 , which are consistent the observational thermal evolution of compact stars [128][129][130]. If a large C 1 is adopted for the perturbation model, QM will emerge at small densities and lead to the quark DU processes. Meanwhile, as indicated in Fig. 14, the nucleon DU processes will take place if we adopt PKDD for HM, where the triangle inequality is fulfilled due to the reduction of electron chemical potentials with the appearance of QM. A more detailed investigation on different combinations of hadronic EOSs and surface tensions is indicated in Fig. 15, where the critical masses M DU for both nucleon and quark DU processes are presented. For these where DU processes never take place, we take M DU = M max . By applying the constraint M DU > 1.5M , it is found that most combinations of hadronic EOSs with quark EOSs determined by perturbation model at C 1 = 3 are not permitted due to an early emergence of QM at small hadron-quark transition densities n T < ∼ 0.3 fm −3 , which lead to the quark DU processes. Meanwhile, larger values of n T are obtained with the hadronic EOSs VM, VMΛ and APR, and consequently the quark DU processes do not occur if large surface tension values are adopted. In such cases, an early emergence of QM at n T < ∼ 0.3 fm −3 is prohibited by the DU criterion. At C 1 = 2.5, the small surface tension values Σ = 0 and 14.3∆n + 1.3 are excluded for the hadronic EOSs VM and VMΛ due to the occurrence of quark DU processes. Similarly, Σ = 20 MeV/fm 2 , Σ MRE , and 0.3Σ MRE are not permitted for the hadronic EOS APR. As indicated in the upper panel of Fig. 15, the hadronic EOS PKDD is excluded since nucleon DU processes always occur in typical compact stars, which rules out the third family of compact stars in Fig. 13. It is worth mentioning that the color superconductivity of quark matter will effectively hinder the quark DU processes [125], so that the cases with M DU < 1.5M in the lower panel of Fig. 15 may not necessarily lead to a fast cooling and the tension with the observational data can be eased. For example, if QM forms a two-flavor superconducting phase, the cooling history of a hybrid star with a large quark core may be consistent with the X-ray data [133]. Note that in the extreme scenario where hybrid stars are comprised almost entirely of QM in the color-flavor-locked phase, heat capacity would be too low to be consistent with observations [134,135].

VI. CONCLUSION
In this work we investigate systematically the possible hadron-quark deconfinement phase transition in dense stellar matter, and its influence on compact star structures. For the hadronic phase, we adopt in total 10 different EOSs, i.e., 8 nuclear EOSs (TM1e [66], TM1 [64], PKDD [73], TW99 [43], DDME2 [74], DD2 [75], VM [49], APR [48]) and 2 hyperonic EOSs (TM1Λ [70] and VMΛ [82]), which are predicted by relativistic-meanfield model [47] and variational method with realistic baryon interactions [48,49]. For the quark phase, we adopt 46 EOSs predicted by equivparticle model [50][51][52], perturbation model [44,53,54], and NJL model with vector interactions [55,56]. With the properties of both hadronic matter and quark matter fixed, the structures of quark-hadron mixed phase are obtained assuming a continuous dimensionality as proposed by Ravenhall et al. [57]. The energy contribution due to the quarkhadron interface is treated with a surface tension Σ, where we have taken constant values for Σ as well as those estimated by the multiple reflection expansion method [58][59][60][61] and equivparticle model including both linear confinement and leading-order perturbative interactions [62,63]. The critical surface tension Σ c that accounts for the energy reduction due to the relocation of charged particles is estimated for various combinations of quark and hadronic EOSs. It is found that in most cases we have Σ < Σ c , where inhomogeneous structures for the quark-hadron mixed phase are favored.
As we increase the density of hadronic matter, quark matter will emerge and forms a quark-hadron mixed phase. By minimizing the energy density at given baryon number density, we have obtained the radius, dimensionality, and quark fraction of MP. It is found that the obtained radius normally ranges from ∼1 fm to ∼10 fm, and is increasing with Σ. The radius evolves more smoothly with density if the dimensionality changes continuously. Adopting various combinations of hadronic and quark EOSs along with different values of surface tensions, the quark fraction usually increases monotonically and turns into a pure quark phase. The corresponding EOSs for hybrid star matter are obtained, which predict the structures of compact stars by solving the TOV equation. It is found that the correlation between radius and tidal deformability in traditional neutron stars [88,120] preserves in hybrid stars. Once quark matter emerges inside compact stars, the quark-hadron interface plays an important role on their structures. The surface tension Σ estimated with the multiple reflection expansion method or equivparticle model increases with density, which predicts stiffer EOSs for the quark-hadron mixed phase and increases the maximum mass of hybrid stars. The hyperons are suppressed if we adopt a quark model that predicts relatively small energy per baryon of quark matter at small densities. Based on various constraints of nuclear physics, causality limit, and pulsar observations, we obtain the permitted parameter sets that are consistent with observation. It is found that the quark phase persists inside the most massive compact stars in almost all the permitted cases. Meanwhile, comparing with higher density regions, the variation of pressure is insignificant at n < ∼ 0.8 fm −3 irrespective of the emergence of quark matter. The current constraints can be further improved based on the thermal evolution of compact stars, which rules out an early emergence of quark matter at densities smaller than 0.3 fm −3 in the absence of color superconductivity.