QCD phase transition and equation of state of stellar strong interaction matter via Dyson-Schwinger equation approach

We study the phase structure and phase transition of cold dense QCD matter via the Dyson-Schwinger equation approach. We take the rainbow approximation and the Gaussian-type gluon model. In order to guarantee that the quark number density begins to appear at the nuclear liquid-gas phase transition chemical potential, we propose a chemical potential dependent modification factor for the gluon model. We find that for the iso-symmetric quark matter, the modification reduces the chemical potential of the phase coexistence region of the first--order phase transition. We also implement the relativistic mean field theory to describe the hadron matter, and make use of the Maxwell and Gibbs construction method to study the phase transition of beta--equilibrium and charge neutral matter in compact stars. The results show that the phase transition will not happen in case of the Gaussian--type gluon model without any modification. The results also indicate that the upper boundary of the coexistence region should be larger than the current Nambu solution existing region. We also calculate the mass-radius relation of the compact stars, and find that the hadron-quark phase transition happens at too high chemical potential so that the maximum mass of the compact star is hardly affected by the hadron-quark phase transition.

It is well known that lattice QCD is principally only valid at zero chemical potential, since it endures the notorious "sign problem" at finite chemical potential. Although there are extrapolation methods such as the Taylor expansion [69][70][71][72][73], the imaginary chemical potential [74][75][76][77][78][79][80], and other approach (e.g., Ref [81]), the lattice QCD can still only deal with small chemical potential µ/T < 2. Therefore, the theoretical study of the phase transition at zero temperature and finite chemical potential is not as clear as that in zero chemical potential.
The problem also exists on the experimental side. There are colliders such as RHIC, FAIR, NICA and HIAF which aim at studying the structure of hot or warm dense QCD matter, but terrestrial experiments are now not able to create cold dense QCD matter. The way to study the possible phase transition at large chemical po-tential has then to resort the astronomical observations. There have been several very heavy compact stars with a mass over ∼ 2M ⊙ discovered [82][83][84][85][86][87], and the detection of gravitational wave from binary neutron star merger also provide constraints on the equation of state (EOS) and hence the possible phase transition inside compact stars [88][89][90][91][92][93][94][95][96][97]. Still, the detail of the phase transition is hidden inside compact stars.
However, in spite of the fact that the DS equation approach recovers successfully many properties of the DCSB phase, it is still far from being able to describe the real-world hadron matter. For example, one of the most important properties of hadron matter is the emergence of matter from vacuum, which corresponds to the firstorder nuclear liquid-gas phase transition. In DS equation calculations, it has been shown that the quark number density and quark condensate remains the same as in vacuum up to a critical chemical potential [21,[125][126][127], but this critical chemical potential is not in accordance with the liquid-gas phase transition chemical potential. Therefore, in this paper, we propose a chemical potential dependent modification on the coupling strength of the model in the approach to describe the phase transition in the cold dense matter, and construct the EOS of the hybrid star matter. The mass-radius relation of the highly massive compact star can be described well.
The remainder of this paper is organized as follow. After this introduction, in Sec. II, we reiterate briefly the DS equation approach at zero temperature and finite chemical potential, and propose the chemical potential dependent modification on the model. In Sec. III, we present the numerical results with the DS equation being solved, including the phase transition of iso-symmetric matter as well as the β-equilibrium and charge neutral cold dense matter. We also calculate the EOS and the mass-radius relation of the compact star. In Sec. IV, we give a summary and some remarks. In this section, we describe briefly the DS equation approach at zero temperature and finite chemical potential.
The starting point of the DS equation approach is the gap equation: where S(p; µ) is the quark propagator, Σ(p; µ) is the renormalized self-energy of the quark: where Λ is the translationally regularized integral, Λ is the regularization mass-scale. g(µ) is the strength of the coupling, D ρσ is the dressed gluon propagator, Γ a σ is the dressed quark-gluon vertex, λ a is the Gell-Mann matrix, and m q is the current mass of the quark. Z 1,2 are the renormalization constants. In this paper, we take a model that ultraviolet integration is finite, so we can move the renormalization point to infinity and take Z 1,2 = 1.
At finite chemical potential, the quark propagator can be decomposed according to the Lorentz structure as: with u = (0, iµ). A complete decomposition should include another term proportional to σ µν , but this term contributes little, and is usually omitted [21,67]. At zero chemical potential, a commonly used ansatz for the dressed gluon propagator and the dressed quarkgluon interaction vertex is: G(k 2 ) is the effective interaction to be introduced in the model, and Γ σ is the quark-gluon vertex. In this paper, we adopt the rainbow approximation which is the leading-order term in a symmetry preserving approximation scheme [128,129], For the interaction part, we adopt the Gaussian type effective interaction : where D and ω are the parameters of the model. In this paper we take ω = 0.5 GeV and D = 1.0 GeV 2 which are determined with the pion mass m π = 0.14 GeV and decay constant f π = 0.091 GeV with m q = 5 MeV for u and d quark [130]. And the current mass of the strange quark is m q = 115 MeV by fitting the kaon mass m K = 0.492 GeV [131].
The DS equation has multiple solutions. In chiral limit (zero current quark mass) and zero chemical potential, there is a solution with B(p) ≡ 0, and a solution with B(p) > 0, which are called Wigner solution, Nambu solution, respectively. The Wigner solution corresponds to the dynamical chiral symmetry (DCS) phase, where the quarks are massless. The Nambu solution corresponds to dynamical chiral symmetry breaking (DCSB) phase, since the mass function M (p) = B(p)/A(p) acquires a non-zero value. Although there are hints that there might exist a quarkyonic phase where the quark is confined but the chiral symmetry is preserving, it is usually believed that the DCSB (Nambu) solution corresponds to the confined hadron phase, and the DCS (Wigner) solution corresponds to the deconfined quark phase.
At zero temperature and small quark chemical potential, the DCSB phase should stay in the vacuum ground state until the chemical potential is larger than some critical value. This property is known as "Silver-Blaze" property of QCD. At baryon chemical potential µ B = 923 MeV, which is the proton mass minus the binding energy, there will be a liquid-gas phase transition where nucleons begin to emerge from vacuum, and the Silver-Blaze property will be broken beyond this chemical potential.

B. Quark Number Density and Silver-Blaze Property
After solving the DS equation at some quark chemical potential, the number density of the quarks can be obtained through [119] : where f q is the distribution function and reads where the trace is for the spinor index.
The integration in Eq. (9) can be converted to a con-tour integral on the complex plane ofp 4 : where In the last line of Eq. (10), the integral −∞ ∞ is zero since the kernel is an odd function, and the last two integrals are also zero because the denominator M diverges at infinity. Therefore, the value of distribution function f q is determined by the poles inside the contour. The pole corresponds to the zero in the denominator M:p A schematic figure of the pole and the contour is shown in the upper panel of Fig. 1. For Nambu solution, since B is non-zero, we can see from Eq. (12) that the imaginary part of the first pole, µ a , is finite for arbitrary value of p 2 , and for the chemical potential µ < µ a , the integral in Eq. (10) is zero since there is no singularities inside the contour, and the quark number density remains zero. When µ > µ a , the contour begins to include the pole, and the number density becomes nonzero. This is a manifestation of the Silver-Blaze property.

C. Modification for the Gluon Model
In DS equation calculations, the Silver-Blaze property is maintained. The condensation and number density of the Nambu solution remains the same as in vacuum for small quark chemical potential (see, e.g. Ref. [119,127,132]). However, the critical chemical potential where quark number density becomes non-zero relies on the vertex and the gluon model, and is usually not in accordance with the nuclear liquid-gas phase transition chemical potential.
We notice that, by altering the coupling constant D in Eq. (7), the critical chemical potential can be changed. The schematic characteristic of the procedure is shown in the lower panel of Fig. 1. For a fixed chemical potential, the pole may lie outside the contour of Eq. (10), and the number density is zero. However, by decreasing the coupling constant D, the imaginary part of the pole should also decrease. When the pole moves into the contour, the number density will becomes nonzero.
In order to correctly recover the liquid-gas phase transition chemical potential, we are going to add a chemical potential dependence to the coupling D by multiplying a modification function h(µ) to the coupling constant D defined in Eq. (7), so as to have where D(µ = 0) = 1.0 GeV 2 . There are three constraints on the modification function h(µ): First, we should have h(µ = 0) = 1 since the vacuum value is fixed by meson properties. Second, h(µ = ∞) = 0 to take into account the asymptotic freedom at large chemical potential. And we also requires that h(µ) guarantees that the critical chemical potential is µ B,c = 923 MeV which is in accordance with the emergence of nuclear matter.
Previously, in the study of compact star matter, an exponentially damping function is proposed [120-122, 124, 133] as: where α is a parameter and ω is the same as in the original gluon model. This function is a monotonically decreasing function. In general, a larger coupling constant D corresponds to a larger quark mass function M (p), and a larger critical chemical potential. Meanwhile it reduces the value µ B,c for the phase transition to happen [133]. However, for some vertex and gluon propagator, the µ B,c might be less than 923 MeV (see, e.g.Ref. [132]), which indicates that a non-monotonic modification function is required. In this work, we introduce another type of modification function: The value of the parameters α and β are fixed by requiring µ B,c = 923 MeV (i.e., µ q,c = 307.6 MeV). In the following, we will denote the modification function in Eq. (14), in Eq.(15) as α type, β type, respectively.

A. Determination of parameter
The Nambu solution of the quark's DS equation has a relatively large mass function. Therefore at small chemical potential, the number density of quark in Nambu solution should be zero. The mass function is related to the coupling constant D, so for a fixed chemical potential, the quark number density relies definitely on the coupling strength D.
For D/D 0 = 1, where D 0 = 1.0 GeV 2 is the coupling strength in vacuum and fixed with the meson properties, the number density should be zero [119], which satisfies the Silver-Blaze property. However, because of numerical error, the calculated number density is non-zero and numerically unstable. Therefore, it is hard to determine the critical coupling strength via the number density directly.
Another way to determine the critical coupling strength is by searching for the mass pole. As we have stated in the last section, the emergence of the quark number density corresponds to a mass pole in the propagator. The denominator M, as defined in Eq. (11), is a function of | p| and p 4 , for a fixed value of µ = µ q,c . When the pole enters the contour of the integral in Eq. (10), as illustrated in the lower panel of Fig. 1, there should be a zero M for a certain value of | p| and p 4 , or equivalently, the maximum value of 1/|M| should be divergent. The obtained maximum value of the 1/|M| as a function of coupling constant is shown in Fig. 2. From the figure, it is apparent that 1/|M| has a pole at D c = 0.716D 0 , where the quark number density begins to appear.
The parameter α and β in Eq. (14) and Eq. (15) is fixed by requiring that h(µ q,c ) × D 0 = D c , and we have α = 0.883 and β = 2.714. The calculated variation behavior of the modification coefficient h with respect to the quark chemical potential µ q (i.e., µ B /3) is shown in Fig. 3. The black solid line corresponds to the result with the α-type modification function defined in Eq. ( 14) and the red dashed line corresponds to the β-type modification function defined in Eq. ( 15). Both of the two presently obtained modification functions are monotonically decreasing. However, at small chemical potential, the modification h α (µ) is smaller than the h β (µ), while at large chemical potential h β (µ) decreases to almost zero more rapidly. The two lines cross at µ q = 307.6 MeV, which is required by our assumption. After fixing the parameters, we can solve the DS equations with the Gaussian-type dressed gluon model defined in Eq. (7) without any modification, with the α-type modification and the β-type modification, respectively.

B. Phase Transition of the Iso-symmetric Matter and the Coexistence Region
Because of the asymptotic nature of QCD, there must exist a deconfined phase transition at high density (chemical potential). Regardless of the possible existence of quarkyonic phase, this phase transition is also the chiral phase transition from DCSB phase to DCS phase.
The order parameter of the chiral phase transition reads usually the condensate of quark, which is the trace of the quark propagator. The chiral susceptibility, which is the derivative of the quark condensate with respect to the current quark mass, is taken as the criteria of the phase transition.
The chiral susceptibility is related to the stability of the system [15,26]. If the chiral susceptibility is positive, the system is stable or meta-stable, and the system is unstable if the chiral susceptibility is negative.
In vacuum, the chiral susceptibility of the Nambu solution, which corresponds to the DCSB phase, is positive, while the chiral susceptibility of the Wigner solution, which corresponds to the DCS phase, is negative(see, e.g. Ref. [15]). This means that the DCS quark matter is unstable in vacuum, and the system should be described by the Nambu solution. When the chemical potential gets large enough, the Nambu solution will disappear, and the chiral susceptibility of the Wigner solution is positive, this means that the quarks gets deconfined from hadrons and the matter should be in quark phase.
In chiral limit, the quark condensate and the chiral susceptibility is well defined. In case of nonzero current quark mass, however, both the quark condensate and the chiral susceptibility with the original definition in tracing the propagator diverges drastically. There have been several subtract schemes to eliminate the divergence (see, e.g., Refs. [2,3,12,15,21,26,66], but it is quite convenient to simply make use of the Lorentz scalar part of the quark propagator at zero momentum, i.e., taking B(p = 0) as a representative of the order parameter, and χ m ≡ dB(p=0) dm0 as chiral susceptibility, as has been done in many literatures. Specifically, it has been shown that, for the phase transition at zero chemical potential and high temperature, the value of the pseudo-critical temperature would be slightly different by using this criterion [26]. At zero temperature and high chemical potential, however, the chiral susceptibility will diverge at the phase boundary, and the phase transition chemical potential determined with this definition is not different from those fixed via other forms of the defination.
In Fig. 4, we show the calculated chiral susceptibilities of both the Nambu and the Wigner solutions with the different gluon models, and we have neglected the solution where the chiral susceptibility is negative. As we can see from the figure, for all the three models, there exists a (quark) chemical potential region where both the Nambu solution and the Wigner solution have positive chiral susceptibility, and this region is usually referred to as the "phase coexistence region", where the DCSB hadron matter and the DCS quark matter appear simultaneously.
In details, in case of the Gaussian-type gluon model without chemical potential dependence, the coexistence region is µ q ∈ [0.322, 0.500] GeV. In case of the Gaussiantype gluon model with the α-type modification, the coexistence region is µ q ∈ [0.271, 0.309] GeV, and for the gluon model with the β-type modification, the coexistence region is µ q ∈ [0.276, 0.308] GeV. Both the upper and lower boundaries of the coexistence region are smaller than the corresponding one in case of without the modification, and the coexistence region gets narrower.  14), and the blue lines represent the results using the β-type modification function defined in Eq. (15). The solid lines correspond to that of the Nambu solution, and the dashed lines correspond to that of the Wigner solution.
This is because our modification coefficient generally reduces the interaction strength, and it will be easier for the quark to be released from hadrons.

C. β-equilibrium Quark Matter
In the last subsection, we have shown the calculated phase boundary and phase transition for the isosymmetric matter, i.e., the chemical potentials of the u and d quarks are the same. However, it is now impossible for us to get high density matter in terrestrial experiment. The most common way to study the phase structure and the phase transition in the high density region resorts then studying the matter in compact stars. The compact star matter is usually in β-equilibrium and charge neutral. Therefore, we are now going to study the β-equilibrium and charge neutral quark matter.
To study of the property of compact star matter composed of both hadrons and quarks (hybrid matter), it is presently common to implementing the theoretical approach of hadron matter and that of quark model separately to study hadron, quark phase, respectively, and combine them with construction schemes.
In this section, we calculate the property of quark matter by solving the DS equation of quark, and only keep the Wigner solution, since the Nambu solution should corresponds to the hadron phase. For the hadron phase, however, we take the relativistic mean field (RMF) theory with TW-99 parameterization [134], which is described in the appendix.
The calculated quark chemical potential dependence of the number density of u and d quarks in the Wigner solution is shown in Fig. 5. It is clear that the result with the β-type modification has the largest number density. This is physically reasonable, because the β-type modification has the smallest coupling strength in the relevant chemical potential region, and the dressed mass is generally smaller corresponding to a smaller coupling strength, and therefore the number density of the quarks with smaller dressed mass is larger.
In the study of compact star matter, it is also necessary to take into consideration the effect of strange quarks. The study of the solutions of DS equation has revealed that, for a fixed coupling strength, there exists a critical current quark mass, above which the Wigner solution will disappear [2,3,6]. For a coupling strength fixed in vacuum, the strange quark mass is well above the critical mass, which means that the strange quark does not have the Wigner solution for the Gaussian-type gluon model without modification, whose coupling strength is the same as that in vacuum at any chemical potential. In case of the gluon model with the modifications, however, since the Wigner solution must exist at large enough chemical potential where the coupling strength is quite small, we can then get the strange quark number density above some critical chemical potential.
The calculated number density of the the strange quarks as a function of quark chemical potential is shown in Fig. 6. One can notice from Fig. 6 that, in case of the gluon model with the α-type modification, the critical chemical potential for the strange quark to appear is µ s,c = 0.461 GeV, while for the gluon model with the β-type modification, the critical chemical potential is µ s,c = 0.420 GeV. The critical chemical potential in case of the β type modification is smaller, because it has smaller coupling strength in the relevant chemical potential region. Also, the strange quark density for the β-modification is larger, with the same reason for that the u and d quark number density is larger for the β-type modification.
After calculating the number density of the quarks, the pressure of each flavor of the quark at zero temperature can be obtained by integrating the number density: Mathematically, the starting point of the integration µ q,0 can be any value, but in our calculation, it is more convenient to choose µ u/d,0 to be the value of the left boundary of the coexistence region, and µ s,0 to be the critical chemical potential where the Wigner solution for the s quark to appear. In Ref. [119], the pressure difference between the Nambu and Wigner solutions is calculated by using the "steepest-descent" approximation, and since the Nambu solution corresponds to the vacuum at µ q,0 , the initial value of the integration can be chosen as P q (µ q,0 ) = P W (µ q,0 )−P N (µ q,0 ), where P W and P N is the pressure of the Wigner solution, Nambu solution, respectively. By adopting the result from Ref. [119], we have P q (µ q,0 ) = 2.58×10 −4 , 4.32×10 −4 and 4.19×10 −4 GeV −4 in case of the Gaussian-type gluon model without modification, with the α-type modification and with the β-type modification, respectively. For the s quark, we simply choose P s (µ s,0 ) = 0 as in Ref. [120]. The total pressure of the quark matter is the sum of the pressure of each flavor of the quarks: whereP and In compact star matter, the contributions from leptons is also important. Usually, the electron and muon are considered. The number density for the leptons is: where k 2 F l = µ 2 l − m 2 l for l = e − , µ − . In this work, we take m e = 0.511 MeV and m µ − = 105 MeV.
The quark matter in a compact star should be in βequilibrium and electric charge neutral, so we have: And we have the baryon density and chemical potential as: With these relations, we can calculate the properties of the β-equilibrium and charge neutral quark matter with a given baryon chemical potential (baryon density).
Except for solving the DS equation for the quark matter and the RMF for the hadron matter, we also need to make use of some construction scheme to describe the hadron-quark phase transition. The condition for the phase transition to occur is the chemical potential and the pressure in the two phase are equal, i.e., where the footnote H and Q denote the hadron and quark sector, and the pressure as a function of baryon chemical potential is determined by the hadron and quark model we implemented. By solving the above equation, we can get the chemical potential corresponding to the hadronquark phase transition. The above described construction scheme is called the "Maxwell construction", which provides a straight forward way to describe the phase transition, and is widely used. However, it only considers one chemical potential, the baryon chemical potential. In the compact star matter, there are two conservation numbers, the baryon number and the charge number. In turn, there are two chemical potentials, the baryon chemical potential and the charge chemical potential. We should then implement the Gibbs construction to take into account both of the chemical potentials.
In Gibbs construction, it assumes that there is a mix phase region (phase coexistence region) where both quark and hadron exist simultaneously, and the baryon chemical potential (µ B ) and the charge chemical potential(µ e ) are the same in both the quark and the hadron phases.
In the phase coexistence region, the pressure of the two phases are the same. And though the two phases are no longer charge neutral separately, there will be a global charge neutral condition. If we define the quark fraction χ ∈ [0, 1], the phase transition conditions can be expressed as: where p H , p Q is the pressure of the hadron, the quark phase, respectively, which is a function of both the µ B and the µ e . n c H , n c Q is the charge density of each of the two phases, which is determined by the corresponding hadron and quark model.
Then, combining Eqs. (26) and (27), together with the field equations of the two phases, we can obtain the µ B and µ e with any given quark fraction χ. By taking χ = 0, 1, we can calculate the phase boundary of the coexistence region under the charge neutral and βequilibrium condition.
The calculated baryon chemical potential dependence of the pressure of each kind matter is shown in Fig. 7. From the figure, we can see that, among the three quark models, the one with the β-type modification has the largest pressure, which is in accordance with that the βtype modification has the largest quark number density.
In Fig. 7, the black solid line and gray dotted line represent the results of the calculated hadron matter pressure via the TW-99 RMF model described in the appendix. The black solid line corresponds to the nuclear matter that only proton, neutron and leptons are considered, while the gray dotted line is the one with the contribution of the hyperons having been taken into consideration. As we have mentioned, for the hadron-quark phase transition to take place, the pressure and the chemical potential of hadron matter and those of the quark matter must be simultaneously the same, respectively. This means that the P -µ B curve of the hadron matter and that of the quark matter must have a cross point. However, Fig. 7 manifests clearly that, in case of the Gaussian-type gluon model without any modification, the pressure of the quark matter is too small and there is no cross point with the line of hadron matter. This means that, for the Gaussian-type gluon model without modifications, the coupling strength in the relevant chemical potential region is too strong and should be excluded, since the hadron-quark phase transition is forbidden in the model, which is physically unlikely. This further indicates the necessity of taking into consideration of the modification factor. For the Gaussian-type gluon model with the α-and β-type modification, the phase transition chemical potential for the β-equilibrium and charge neutral matter is µ B = 1.94 GeV, 1.71 GeV, respectively. And the gray dotted line does not have cross point with all the quark lines, which means that the phase transition will not happen after the inclusion of hyperons. Therefore, we will not include the hyperon effect in the following. By implementing the Gibbs construction, we can obtain that the phase coexistence region of the βequilibrium and global charge neutral matter is µ B ∈ [1.69, 2.02] GeV, µ B ∈ [1.48, 1.79] GeV in case of the model with the α-type, β-type modification, respectively, as shown in Fig. 7.
The obtained phase coexistence region and the phase boundary here is different from what we have shown in Sec. III B, but with some simple relations: where µ G,1 , µ G,2 is the baryon chemical potential cor-responding to the left, the right boundary of the coexistence region from the Gibbs construction, µ W,c is the quark chemical potential corresponding to the divergence of the chiral susceptibility of the Wigner solution (DCS phase), and µ N,c corresponds to the divergence of the chiral susceptibility of the Nambu solution (DCSB phase). The first inequality in Eq.(28) is satisfied well in our calculation, while the last inequality is not satisfied. This deficiency comes from both the insufficiency of the hadron model and the quark model.
By solving the DS equation of the quark(s), we are calculating the property of uniform quark matter. However, the quarks in the hadron matter is not uniformly distributed since they are clustered as hadrons. So the Nambu solution of DS equation is different from the real world hadron matter after the appearance of quark number density.
The RMF model also has its deficiency. The parameters of the hadron matter model are fixed by fitting the property of the hadron matter at saturation density or lower densities where terrestrial experiments are able to create. Beyond the saturation density, different models give distinct results. Still, the RMF model is more reliable than the Nambu solution of the DS equation in describing the hadron matter, since the solution of the DS equation has not taken into account the effects of the non-uniform distribution of the quark, the surface of the cluster, and so forth. We then conclude that the right boundary of the phase coexistence region of the matter is at µ N,c ≥ µ G,2 /3 = 0.673 GeV, µ N,c ≥ 0.597 GeV for the α-type modified, β-type modified gluon model, respectively.

D. Equation of State and Mass-radius Relation of Compact Star
In order to study the hadron-quark phase transition at zero temperature and high chemical potential (density), one has to make use of the compact star observations to check our theoretical results since it is not able to generate such dense matter on earth.
The most important observable of compact stars is the maximum mass, which is related directly to the EOS P = P (ε) of the dense matter.
The energy density ε of the dense matter under Maxwell construction is: where ε H and ε Q is the energy density of the charge neutral hadron matter and quark matter, respectively. µ B,c is the phase transition baryon chemical potential. The pressure of the dense matter under Maxwell construction is: where P H is the pressure of the hadron matter and P Q is the pressure of the quark matter. For Gibbs construction, the EOS should be decomposed into three parts: µ B ≤ µ G,1 region, µ G,1 ≤ µ B ≤ µ G,2 region and µ B ≥ µ G,2 region, where µ G,1 and µ G,2 is the left and right boundary of the phase coexistence region (mixed phase).
The energy density of the mixed phase consists of the contribution of the two phases.
where the footnote M, Q and H corresponds to the mixed, the quark, the hadron phase, respectively. For a fixed µ B , µ e is calculated by solving Eq.(26) and Eq. (27). And the pressure for the mixed phase reads: The energy density and the pressure under the Gibbs construction are: The calculated EOSs of the pure hadron, pure quark and hybrid (mixed phase) matter are illustrated in Fig. 8. As we can recognize from Fig. 8, the pure quark matter has a softer EOS than the pure hadron matter, and both the two different construction schemes connect the EOSs of the hadron matter and the quark matter. In case of the Maxwell construction, there is a density region where the pressure is constant, which corresponds to a phase transition with a constant chemical potential. However, in compact stars, the pressure must have a gradient in order to resist its own gravity, so that such a constant pressure region will not appear inside compact star. The phase transition under the Maxwell construction corresponds to a sudden change in the energy density between the quark core and the hadron shell.
The mass-radius relation of a compact star can be calculated by solving the Tolman-Oppenheimer-Volkov (TOV) equation: where G is gravitational constant and m(r) is the mass inside a radius r: Then given the EOS as input, and with a given center density, we can integrate the TOV equation from inside out to get the mass and radius of the compact star. For the pure hadron star and the hybrid star, at small density region, we make use of the Baym-Pethick-Sutherland (BPS) EOS [135]. For the pure quark star, we integrate to the surface where the pressure is zero.
The calculated mass-radius relation of the compact star in the cases of our consideration is shown in Fig. 9. From Fig. 9, one can notice apparently that the pure quark star has a much smaller maximum mass than the pure hadron star and hybrid star, and the maximum mass is 1.07M ⊙ , 0.95M ⊙ for the α-, β-type modification, respectively.
For pure hadron star, the maximum mass is 2.06M ⊙ . From Fig. 9, we can see that at large radius and small mass, the mass-radius curves for hybrid star and pure hadron star are exactly the same, because for hybrid star with relatively small mass, the inner density is not large enough for the quark matter to appear.
When the mass of the hybrid star is large enough, the quark matter begins to appear in the composing matter of the star, and the mass-radius curves begin to deviate from the pure hadron one. As can be seen from Fig. 9, for both α and β modification, and for both Gibbs and Maxwell construction, the deviation begin to appear at over ∼ 2M ⊙ , which means that the hadron-quark phase transition happens at very high density and the massradius relation is determined mainly by the EOS of the hadron sector. If the quark sector is described by DS equation with the α-type modification, the maximum mass will be 2.06M ⊙ for both Gibbs and Maxwell construction. If the quark sector is described by DS equation with the β-type modification, the maximum will be 2.00M ⊙ for Gibbs construction and 2.05M ⊙ for Maxwell construction.

IV. SUMMARY
The Dyson-Schwinger equation has been implemented to study the phase transition and the phase structure of QCD matter. However, the possible dependence of coupling strength on the quark chemical potential is not well studied. In this paper, we take into account several types of the chemical potential dependence, and make use of the nuclear liquid-gas phase transition to fix the parameter(s).
For the modification factor we propose, the boundary of the phase coexistence region of the iso-symmetric matter is reduced than that without modification. The phase coexistence region is µ q ∈ [0.322, 0.500] GeV for the Gaussian-type gluon model without any modification, and becomes µ q ∈ [0.271, 0.309] GeV, [0.276, 0.308] GeV in case of the α-, β-type modification, respectively.
However, by taking into account the β-equilibrium and charge neutral condition for the matter inside compact stars, and comparing the results from the DS equation calculation for the quark sector and the RMF hadron model, we find that the quark pressure is too low for the Gaussian-type gluon model without any modification, and the hadron-quark phase transition is prohibited, which is physically unlikely. This means that the adoption of modification is necessary.
By implementing the RMF model for hadron matter and the DS equation calculation for the quark sector, we find that the phase transition chemical potential under the Maxwell construction is µ B = 1.94, 1.71 GeV for the α-, β-type modification, respectively. The phase coexistence region under the Gibbs construction is µ B ∈ [1.69, 2.02] GeV, [1.48, 1.79] GeV in case of the α-, βtype modification, respectively. The upper boundary of this mixed phase region is not compatible with the iso-symmetric result, and discrepancy comes from the uncertainty of both the DS equation approach and the RMF model in the phase coexistence region. Since the RMF model is more accurate in describing the hadron phase, we conclude that this boundary indicates that µ N,c ≥ 0.673 GeV, 0.579 GeV, for the α-β-type modification, respectively, where µ N,c is the quark chemical potential at which the Nambu solution (DCSB phase) disappears.
We also obtain the mass-radius relation for compact stars in the cases we considered. The maximum mass of the pure hadron star is 2.06M ⊙ . If the quark sector is described using α-type modification, the maximum mass of the hybrid star is the same as that of pure hadron star. If the quark sector is described by implementing βtype modification, the maximum mass of the hybrid star is 2.00M ⊙ and 2.05M ⊙ for the Gibbs construction, the Maxwell construction, respectively. Therefore, for hybrid stars, the maximum mass is almost the same since the phase transition happens at very high density. However, the maximum mass of the pure quark star is only about 1M ⊙ . In this work, the possible appearance of hyperons in the dense star matter is not considered, because it prevents from the happening of the hadron-quark phase transition. This problem can be fixed by implementing the 3-window construction method instead of the Gibbs and Maxwell constructions, as has been done in our previous work [133].
Although we have made use of the liquid-gas phase transition chemical potential to fix the parameters of the modification function in our DS equation claculation, it is still impossible to recover all the properties of the hadron matter from first-principle calculation. For example, there should be a density gap at liquid-gas phase transition. Also, the saturation density of nuclear matter cannot be recovered by the Nambu solution of the DS equation. To fix these problems, taking into account the dressing effect of the quark-gluon interaction vertex or/and the interface effect between quarks and hadrons (clustered quarks) in the DS equation approach should be an efficient way. The related work is under progress. The relativistic mean field (RMF) theory [136][137][138] is a successful theory which describes the property of the hadron matter in the density region relevant in compact stars [139,140]. There are many different RMF models, with different parameterizations which are calibrated by the properties of dense nuclear matter. However, most of the models satisfy only a small number of experimental data. In Ref. [141], hundreds of RMF models are taken to fit the experimental data, and the TW-99 model [134] is found to be one of the best model that can reproduce most of the nuclear matter properties. Also, TW-99 model generates an EOS that is stiff enough to support a 2-solar-mass neutron star [142], and we have taken such a model to construct a massive hybrid star in our previous work [133].
The Lagrangian of the TW-99 model for the hadron matter is where L B is the Lagrangian of free baryons, which reads where i = p, n stands for proton and neutron. If we want to take into consideration the effect of hyperons, we can take i = p, n, Λ, Σ ±,0 , Ξ −,0 for the baryon octect. L M in Eq.(A1) is the Lagrangian of mesons, where σ, ω µ , and ρ µ are the isoscalar-scalar, isoscalarvector and isovector-vector meson field, respectively.
The L int in Eq. (A1) is the Lagrangian describing the interactions between baryons which are realized by exchanging the mesons: where g iB for i = σ, ω, ρ are the coupling strength parameters between baryons and mesons, which depend on the baryon density.
In some other literatures, the self-interaction of σmeson, the cross interaction between different kind mesons and the effect of the isovector-scalar δ-meson are included explicitly (see, e.g., the review in Refs. [140,141]). In TW-99 parameterization, however, all these terms are taken zero, and their effects are represented in the density dependence of the coupling constants. For nucleons, the coupling constants are where n B is the baryon density, n s is the saturation nuclear matter density and x = n B /n s . The density function can be chosen as [134]: where the parameters a i , b i , c i , d i and g iN (ρ sat ) are fixed by fitting the properties of the nuclear matter at the saturation density, and their values are listed in Table I. For hyperons, we represent them with the relation between the hyperon coupling and the nucleon coupling as: . On the basis of hypernuclei experimental data, we choose them as those in Refs. [139,142]: χ σ = 0.7, χ ω = χ ρ = 0.783.
The L lep is the Lagrangian for leptons, which are treated as free fermions: and we include only the electron and muon in this paper. The field equations can be derived by differentiating the Lagrangian. Under RMF approximation, the system is assumed to be in the static, uniform ground state. The partial derivatives of the mesons all vanish, only the 0component of the vector meson and the 3rd-component of the isovector meson survive and can be replaced with the corresponding expectation values. The field equations of the mesons are then: where τ 3B is the 3rd-component of the isospin of baryon B.
The Σ R µ is called the "rearranging" term, which appears because of the density-dependence of the coupling constant, and reads where j µ =Ψ B γ µ Ψ B is the baryon current. Under the EoM of Eq. (A12), the baryons behave as quasi-particles with effective mass and effective chemical potential: One can then get the baryon (number) density: where k F B = µ * 2 B − m * 2 B is the Fermi momentum of the particle, γ B = 2 is the spin degeneracy. And the scalar density is: The density of a kind of leptons is the same as that for baryons, except that the effective mass and the effective chemical potential should be replaced with the corresponding mass and chemical potential of the lepton.
The matter in the star composed of hadrons should be in β-equilibrium. Since there are two conservative charge numbers: the baryon number and the electric charge number, all the chemical potential can be expressed with the baryon chemical potential and the electron chemical potential: where B and Q is the baryon number, electric charge number for the particle i, respectively. Since neutron has one baryon number and zero charge number, we can take µ B = µ n where µ n is the neutron chemical potential. Then, combining Eqs. (A9), (A10), (A11), (A14), (A15), (A16), (A17), (A18), (20) and (A19), together with the charge neutral condition: one can determine the ingredients and the properties of the hadron matter with any given baryon density n B . The EOS of the hadron matter can be calculated from the energy-momentum tensor: The energy density ε is: where the contribution of the baryon B to the energy density is: The contribution of the leptons to the energy density can be written in the similar form as baryons with a spin degeneracy parameter γ l = 2, except that the effective mass and effective chemical potential should be replaced with those of the leptons, respectively.
As for the pressure of the system, we can determine that with the general formula: (A24)