QCD equations of state and speed of sound in neutron stars

Neutron stars are cosmic laboratories to study dense matter in Quantum Chromodynamics (QCD). The observable mass-radius relations of neutron stars are determined by QCD equations of state, and can reflect the properties of QCD phase transitions. In the last decade there have been historical discoveries in neutron stars, the discoveries of two-solar mass neutron stars and neutron star merger events, which have imposed tight constraints on equations of state. While a number of equations of state are constructed to satisfy these constraints, a theoretical challenge is how to reconcile those constructions with the microphysics expected from the hadron physics and in-medium calculations. In this short article we briefly go over recent observations and discuss their implications for dense QCD matter, referring to QCD constraints in the low and high density limits, QCD-like theories, and lattice QCD results for baryon-baryon interactions.


I. INTRODUCTION
Quantum Chromodynamics (QCD) is the theory of strong interactions with quarks and gluons being elementary degrees of freedom (d.o.f.) [1,2]. In a normal circumstance they are confined into hadrons, e.g., protons, neutrons, pions, and so on. These composite particles are effective d.o.f. to describe the physics of nuclei or the vacuum. The situation changes when we heat a gas of hadrons or compress matter made of nuclei; when hadrons overlap, quarks and gluons become natural d.o.f. to describe the physics in extreme environments [3]. These extreme matters appear in the early universe and in astrophysical objects such as neutron stars [4].
Our understandings seem matured for some extreme environments. In hot QCD, a heated hadron resonance gas (HRG) continuously transforms into a quark-gluonplasma (QGP), in spite of apparent differences of these two phases. Experimentally such hot matter has been created by heavy ion collisions at high energy in which the colliding energy is converted into heat. Combined analyses of the experiments, ab-initio lattice simulations of QCD, and model calculations, together have formed a plausible picture of hot QCD matter [5].
Another extreme environment is cold, dense matter of QCD (for reviews, see, e.g. [6][7][8]). Compressing matter but keeping temperature low should lead to quark matter, as proposed in 70's [9,10]. How nuclear matter changes into quark matter is still unknown, as theoretical framework for such transitions is not established; in particular the sign problem prevents us from performing the lattice Monte-Carlo simulations. Only exceptions are (i) the domain of nuclei around baryon density, n B = n 0 0.16 fm −3 , where we have nuclear experiments, and (ii) very high density, n B ∼ 40 n 0 , where the asymptotic freedom in QCD allows weak coupling calculations [11,12]. Cold matter denser than nuclear matter cannot be created in terrestrial experiments such as heavy ion collisions which inevitably produce heat.
In this respect neutron stars (NSs) are unique cosmic laboratories to study cold, dense matter of QCD beyond the nuclear regime. A NS is an extremely dense object which contains a solar mass (M ) within the radius of ∼ 12 km [13]. Having a very large energy within a small domain, the gravitational is very strong. A star would collapse into a blackhole (BH) unless the gravitational force is balanced with the matter pressure. The structure of NSs is determined by the Einstein equation coupled to the QCD (plus electroweak) energy-momentum tensor, which is reduced to the Tolman-Oppenheimer-Volkoff (TOV) equation for a spherical star. The electroweak interactions affect the matter composition through the charge neutrality and β-equilibrium conditions on equations of state (EoS). The mass-radius (M -R) relations of NSs are the most basic observables which have one-toone correspondence with neutron star EoS where QCD plays dominant roles. The observational determination of M -R curves is not straightforward, as we have to look for good signals from the universe. Nevertheless, the last decade had historical discoveries: the existence of two-solar mass (2M ) pulsars [14][15][16][17][18] and a detection of gravitational waves (GWs) from a NS merger (GW170817) [19] with the associated electromagnetic (EM) counterparts [20][21][22][23][24][25][26]. These findings allow us to delineate the properties of dense QCD matter. The on-going and up-coming observational programs will be listed in Sec.II.
We are getting better constraints on M -R relations or EoS, but they do not directly tell the matter composition. We need to look into the microphysics behind EoS, taking into account the constraints from the nuclear physics near the saturation density and the finite size effects of hadrons. Taking typical estimates on these quantities (reviewed in Sec.III), we can infer the density domain for a given point of M and R. For a canonical NS with the mass ∼ 1.4M , the density is 2-3n 0 , and hence may be characterized within the hadronic regime, while for 2M NSs the core density may go beyond the density, n B ∼ 5n 0 , at which baryons of the radii 0.5-0.8 fm overlap, and hence the core may accommodate quark matter. If a massive NS indeed contains quark matter, there should be hadron-to-quark matter phase transitions near the core. The nature of the phase transition has not been determined from the current observations; modeling based on different orders of phase transitions can be arranged to pass the observational constraints (see Sec. IV). But the observations have certainly constrained the strength of the phase transition and our model building. The first order phase transition, if exists at low temperatures, would imply the existence of the QCD critical end point phase at finite temperature; this is one of important targets in the beam energy scan program of heavy ion experiments [27,28].
In this article we begin with reviewing observational constraints (Sec.II), and combine them with theoretical considerations on matter properties (Sec.III). Particular attention is payed for the speed of sound which differentiates dense QCD matter from conventional matters (Sec.IV). We take the natural unit c = = k B = 1.

II. OBSERVATIONAL CONSTRAINTS
M -R relations -In order to obtain a M -R point, we assign a core density n c , and then integrate matter from the core to the surface until the pressure reaches zero, P (r = R) = 0 [13]. The M -R points at different core densities form a M -R curve for a given EoS. M -R curves also depend on rotations, but in most cases they can be treated in perturbative treatments [29,30]. But the rotation goes beyond the perturbative regime for objects after the NS mergers (see below). For later convenience we call the maximum mass of a non-rotating star M TOV , and the radius of 1.4M NS R 1. 4 .
In NSs, EoS at various densities contribute. But the shape of M -R curves can be largely characterized by EoS at fiducial densities [31]. Shown in Fig. 1 is a typical M -R curve. A low mass star has a core plus a loosely bound crust. The rapid reduction of R within a small increase in M is due to the compression of dilute crust matter. For heavier stars the crust is very thin, and the size of the dense core, which is hard to compress, is observed. This turning point is seen around n c ∼ n 0 . The curve passes the domain of canonical 1.4M NSs at n c ∼ 2-3n 0 , and reaches M > 2M at n c 5n 0 . This illustration suggests that the overall radii of NSs are determined by EoS at low density, n B = 1-3n 0 , while M TOV by EoS at high density, n B 5n 0 . Now we list up observational constraints.
2M Pulsars -One of the historical measurements in the last decade was the Shapiro delay measurement of the pulsar PSR J1614-2230, announced in 2010. The mass was initially estimated as 1.97 +0.04 −0.04 M [14], and after long term observations the estimate has been updated to 1.908 +0.016 −0.016 M (68%CL) [16]. The second precisely measured 2M NS is PSR J0348+0432 whose mass is 2.01 +0.04 −0.04 M [17]. A new Shapiro delay measurement was done for the PSR J0740+6620 where the mass is estimated to be M = 2.14 +0. 10 −0.09 M [18]. These results have established the 2M constraints that require high density EoS to be stiff.
NS mergers -Another historical event is a NS merger, GW170817, from which GWs were detected by the aLIGO and Virgo [19]. Furthermore, this event happened at a rather close distance, 43.8 +2.9 −6.9 Mpc, allowing the measurements of EM signals [20][21][22][23][24][25][26], with which the merger event was analyzed from various aspects. The observations have been compared to general relativistic numerical simulations for various EoS [32,33].
The merger event experiences various stages, and each of which has distinct signals (for a review, e.g. Ref. [34]). In an early spiral phase two NSs are widely separated, and each can be treated as a point particle. This stage informs us the total mass. With long term emissions of GWs, two NSs come close enough to deform the shape of NSs. This is a tidally deformed phase. The deformation of the NSs adds an extra gravity whose net effect is attractive and hence accelerates the merging process. This affects the GW form patterns. The degree of the deformation is characterized by the tidal deformability which has strong correlations with a NS radius; a NS with a larger radius is more easily deformed, leading to a larger tidal deformability. In the end NSs collide. The stage after the collision is called a post merger phase, a highly dynamical stage. There are various signals; the gammaray burst, blue-and red-ejecta, neutrinos, and GWs at high frequency, 1-4 kHz, which reflects rapid oscillations of a compact object [35][36][37]. At beginning, the merged object is differentially rotating, and this rotation prevents the object from immediate gravitational collapse. After short time this differential rotation is braked by viscous effects and magnetic field effects. After this loss of centrifugal effects, the object exceeding the maximum NS mass collapses to a BH. The ejecta in a post-merger phase is sensitive to the compactness of merged objects and the time scale of the gravitational collapse; for example the amount of ejecta would be too little if the merger immediately collapses or is too compact.
In the event GW170817, the total mass before the merger was estimated in good accuracy, 2.73-2.78M , and the GWs in early spiral and tidally deformed phases were measured [19]. Referring to GW templates from numerical simulations, the tidal deformability (or radii of neutron stars) was constrained. The aLIGO-Virgo collaboration estimated radius to be 11.9 +1.4 −1.4 km for each star [38], which is largely consistent with estimates by other studies [39,40]. The current detector did not have enough sensitivity to detect the GWs in the post merger phase. But the EM detection provides information of the post merger phase; the estimated amount of ejecta of ∼ 10 −2 M , which is correlated with the lifetime and compactness of the merger, is used to put the upperand lower-bound on the M max , and the lower-bound on R 1.4 [41], although the bounds depend on our interpretation on the fate of the merger. The first and standard one is the hypermassive NS (HMNS) scenario, in which the merger becomes a metastable star with differential rotations and collapses. The second one is the supramassive NS (SMNS) scenario, in which a HMNS changes into a long-lived, rigidly rotating NS whose maximal mass is greater than the corresponding non-rotating star by a factor 1.2. The threshold mass M th = (2.73-2.78M )/1.2 2.28-2.32M appears as M TOV < M th in the HMNS scenario [42][43][44][45], while as M TOV > M th in the SMNS scenario [46,47]. Based on the HMNS picture, the lower-bound of the radius was also estimated [41,48], although later it was found that the condition can be significantly relaxed by allowing a sufficiently large M TOV for the long life time [49].
After the discovery of GW170817, several (candidates of) NS mergers have been found. The second detection is the GW190425 event with a large total mass 3.4 +0. 3 −0.1 M [50]. This system is exotic because at least one of NSs should be different from canonical NSs with 1.4M . The possibility of BH-NS binary is also not excluded [51,52]. The EM counterparts have not been detected and the tidal deformability is not well constrained. This is in part because this event happens at a far distance, 159 +69 −71 Mpc (about four times larger than the GW170817). Another possible reason is that the total mass is so large that the merger promptly collapsed without producing much ejecta. If so, the absence of the HMNS stage (which happens if the mass is less than ∼ 1.5M TOV ) implies M TOV 2. 27M 3.4M /1.5. This estimate is consistent with the HMNS interpretation of GW170817.
There is also a candidate of BH-NS merger, GW190814 [53], in which the mass for the heavier object is 22.2-24.3 M and should be a BH, while for the lighter one 2.50-2.67M . The second one can be either a BH or (rapidly rotating) NS. If the object is a BH, this is a low mass BH in addition to the previously found candidate with the mass 3.3 +2. 8 −0.7 M [54]. If the object is a rotating NS, the EoS must be stiff enough to make the 2.6M NS stable; this implies M TOV 2.08-2.23M 2.50-2.67M /1.2 [55].
These findings are already quite remarkable, but this is just the beginning. Now we have three powerful GW detectors LIGO-Virgo-KAGRA (KAGRA has just started in 2019), and the IndIGO will start in 2023.
There is also the project to build the third generation of detectors, Einstein-Telescope, for 2030s [56]. Finding more events will eventually lead to the statistical analyses, and some events may allow clear-cut interpretations. Also, if a supernova (SN) event happens in our galaxy as in SN1987A, it is possible to obtain GWs [57] together with significant amount of neutrinos; in Hyper-K [58], the expected count of neutrinos is ∼ 10 4 times greater than for SN1987A. The galactic SN happens a few times in a century, so it might happen within next 10 years.
Neutron Star Interior Composition ExploreR (NICER) -The NICER, launched in 2017, have taken data of X-rays from various pulsars to measure the NS radii and masses simultaneously. The pre-NICER analyses (for a review, see [59]) contained rather difficult estimates of the distance to NS, the assumption of black-body radiation, and atmospheric compositions, which introduced the systematic errors. The NICER got rid of the black-body assumption and distance estimate, by following the time evolution of hotspots on NS surface and performing phase-resolved spectroscopy [59]. Two teams in the NICER individually estimated the radius of the pulsar PSR J0030+0451, reporting R 1 = 13.02 +1.24 . Both estimates tend to be a bit larger than those from GW170817.
In addition, the NICER plans to measure the 2M pulsars, PSR J1614-2230 with the mass 1.91M [16] and PSR J0740+6620 with 2.14M [18]. These measurements will constrain the radii in the high mass region, and hence constrain high density part of EoS.

III. DELINEATING DENSE MATTER
In order to delineate the properties of NS matter, we need studies based on physical pictures on the effective d.o.f. This will remain true even after the M -R relations are precisely determined or the sign problem for lattice simulations is solved. For instance, in case of hot QCD, the transition from a hadronic to a QGP was shown to be crossover by lattice simulations [62], but the detailed physical picture has emerged after the HRG [63] and pQCD [64], with clear-cut effective d.o.f, were used as baselines to diagnose the lattice data [5]. It is also desirable to have such baselines to characterize the NS data. In this section we divide the density domain into four categories, based on plausible effective d.o.f. [65].
For n B 2n 0 -At low density quarks form baryons, whose interactions are mediated by quark-exchanges in the color-singlet channel (meson exchange). Around ∼ n 0 the strange baryons are absent, and we consider only nucleons (the possibility of very early onset of deconfined quark matter will be mentioned in Sec.IV). The matter is dilute enough for nucleons to exchange only few mesons. There are several methods to compute EoS. Two popular approaches seem quite popular: (i) relativistic mean-field (RMF) calculations with in-medium effective interactions; and (ii) non-relativistic many-body calculations based on microscopic bare nuclear forces.
The RMF approach uses relativistic Lagrangian written in terms of baryon and meson fields [66]. The (inmedium) coupling constants are arranged to reproduce the nuclear physics near the saturation and finite nuclei. The clear advantage is its simplicity and flexibility; the method allows us flexible analyses of new experimental data. (For EoS reproducing wide range of phenomenology, see e.g. [67][68][69]). Another utility is that the theory includes relativistic effects by construction, and can be used to describe changes in d.o.f. through mean-field variations. Meanwhile the disadvantage lies in its systematics; the form of Lagrangian is not well-constrained, as there is no good reason to truncate higher orders of fields due to the lack of counting schemes. The EoS can be sensitive to the choice of higher order couplings.
The approach (ii) first prepares bare two-and threebody forces to reproduce experiments for few-body systems, and then use them for many-body calculations which should be non-perturbative to properly handle soft nucleonic excitations from the nucleon Fermi surface (however see also Ref. [70] for perturbative approaches). The advantage is its systematics; there are counting schemes based on either ranges (in potential models [71,72]) or momenta of interactions (Chiral Effective Field Theory [73,74]). Uncertainties in determining the forces can be directly converted into the EoS, and this sharpens questions on many-body effects. Twobody nuclear forces below the energy 350 MeV are well-constrained by two nucleon scattering experiments, while there are more uncertainties in three nucleon interactions, especially at short distance. In this respect the pure neutron matter computations have less uncertainties because the Pauli blocking does not allow three neutrons to overlap. It is challenging to include other d.o.f. toward higher density, the possible changes in nucleon properties, and relativistic effects which are not small already at ∼ 2n 0 [71].
For 2n 0 n B 5n 0 -Around 2n 0 , many-body forces become increasingly important, casting a doubt on the validity of nuclear matter descriptions. Since nuclear forces are mediated by quark exchanges, the importance of many-body forces may lead to the structural changes in nucleons as well [75]. Meanwhile, the density 5n 0 is presumably not high enough for establishing quark matter. The domain 2n 0 n B 5n 0 is not well-understood theoretically even at conceptual level, as the effective d.o.f. is not clear-cut. There are only a few works [75][76][77] that explicitly address the physics in this subtle domain. A recent work [77] discusses how quark wavefunctions can delocalize through quark exchanges among baryons. Assuming that typical quark exchanges take place when meson clouds (valence quark-antiquark pairs) of the thickness ∼ 0.7 fm overlap, some modes begin to percolate already around n B ∼ 1.8n 0 (called soft-deconfinement [77]), and at higher density more modes are gradually delocalized.
For practical construction of EoS, there are at least three descriptions. The first is to extrapolate hadronic EoS to high density. Considering the finite size of hadrons this description would be problematic. The other two descriptions allow the appearance of quarks. The first is to use nucleonic EoS to 2-3n 0 , and switch to quark matter EoS through the first order phase transition, or through the mixed phases in which the transition accompanies various clusters [78][79][80]. The second possibility is the quark-hadron continuity in which hadronic matter is smoothly connected to quark matter. More on this picture will be given in Sec.IV.
For 5n 0 n B 40n 0 -In this regime baryons overlap and the effective d.o.f. should be quarks and gluons, but the density is not high enough for pQCD descriptions. This domain has not been studied in detail. The density 5n 0 0.8 fm −3 is close to those inside of a single hadron, and hence it is reasonable to expect the validity of quasi-particle descriptions [81,82] as in constituent quark models for hadrons. Such constituent quark based discussions [83,84] can also explain the baryon-baryon interactions at short distance measured in the lattice simulations [85]. Furthermore, the energy density vs mechanical pressure inside of a proton [77,86] are found reasonably consistent with neutron star EoS. The mechanical pressure in a proton can be extracted from the gravitational form factor, and can be studied through deeply virtual Compton scatterings [87]. The precise determination of the gravitational form factor is one of important targets in the electron ion collider.
Another hint may be obtained from studies in two-color QCD; for this system lattice MonteCarlo simulations are doable [88][89][90]. There have been seminal works on the phase diagram, EoS, and order parameters. The obvious difference from the three-color case is that baryons are diquarks with which the baryonic matter starts with the Bose Einstein condensation. But such difference may not be essential when density is high; eventually quarks inside of diquarks manifestly establish the Fermi sea, and the condensation remains only around the edge of the Fermi sea as in the BCS theory. Indeed this behavior has been observed on the lattice [88]. The critical temperature of the diquark condensation is high, T s 100 MeV, even at a quark chemical potential of ∼ 1 GeV or n B ∼ 50n 0 . If we accept the BCS relation T s 0.57∆ with ∆ being the size of the gap, then ∆ 180 MeV. This would suggest that gluons remain non-perturbative, and the studies of gluon propagators seem consistent with this point of view [91,92].
For n B 40n 0 -In this domain the pQCD should be valid. The great advantage of this framework is that it contains the error-estimator; by varying the order of α s and/or renormalization scales, one can infer the importance of higher order corrections. The pioneering calculations were done already in 70's [11], and more systematic error estimates have been carried out in Ref. [12]. The result shows that α s expansion does not converge well at 40n 0 , indicating that the matter should be strongly correlated. This pQCD EoS have been used as boundary conditions to construct general curves connecting pQCD and nuclear EoS.

IV. SPEED OF SOUND
The current major constraints concern with M TOV and R 1.4 . While they constrain the EoS at different densities, these high and low density parts constrain each other through the causality constraint and thermodynamic stability; these conditions can be summarized into 0 ≤ c 2 s ≤ 1 where c 2 s = ∂P/∂ε is the squared speed of sound. The causality constraint becomes severer in the case with a smaller NS radius and a larger M TOV , see Refs. [93,94] for comprehensive analyses. Based on the estimates from GW170817 and NICER, below we assume that R 1.4 12-13 km, excluding very stiff EoS for 1-3n 0 . This puts the upperbound on the strength of possible first order phase transition from hadronic to quark matter; if the first order phase transition is too strong, the EoS just after the phase transition are soft, and with the presence of c 2 s ≤ 1 the pressure cannot grow rapidly enough to satisfy the required stiffness at high density.
This situation differentiates the hadronic-quark transition at high density from the hadron-QGP transition at high temperature. In the latter the c 2 s never exceeds 1/3 and has a dip in the crossover region; near the transition, many kinds of non-relativistic resonances, with the masses much larger than temperature, can appear because the entropic effects overcome the Boltzmann suppression. A non-relativistic gas naturally has small c 2 s than a pion gas, and this tendency continues until those resonances merge into a QGP, recovering the conformal behavior. In contrast, cold dense matter does not enjoy such entropic effects, and c 2 s can behave quite differently. Since c 2 s in dense QCD may be quite unique, EoS are often parameterized by c 2 s [95,96]. Shown in Fig.2 are (a bit exaggerated) sketches of P -ε, M -R, and c 2 s -n B relations for three characteristic types of EoS. (In c 2 s -n B relations we also indicate typical nucleonic curves by dashed lines.) The leftmost panel shows the hybrid EoS with the first order phase transition. This case has been most intensively studied. Provided that nuclear EoS are trustable to n B ∼ 1.5n 0 , the stiffness must grow so fast that EoS can remain stiff even after the first order phase transition [97]. Starting with c 2 s ∼ 0.1 in the nuclear domain, its value exceeds the conformal limit 1/3. The presence of the phase transition is reflected in a kink in the M -R curve. In a radical case it may lead to the third family of stars [79,[98][99][100].
The middle case of Fig.2 includes a very rapid growth in c 2 s , from c 2 s ∼ 0.1 to 1/3 at low density. Assuming quick stiffening to happen at 1.1n 0 n B 1.5n 0 , the resulting EoS begin to get stiffened at low density, achieving the required stiffness at high density without invoking c 2 s greater than 1/3. It has been known [101,102] that quark EoS with c 2 s = 1/3 can lead to stiff EoS, if we choose a suitable location for the onset (and neglect the interplay with nuclear matter). The behavior of c 2 s 1/3 is quite different from pure hadronic description, leading to a recent proposal that differentiates quark and hadronic matter by the index γ = c 2 s ε/P ; quark matter is characterized by γ 1.75 [103]. Substituting c 2 s = 1/3, the condition becomes P/ε 0.19; in the domain 1-2n 0 this ratio is substantially larger than the hadronic case with the pressure suppressed by the nucleon mass. The size of quark matter core is naturally large.
In the rightmost of Fig.2, nuclear descriptions are trusted to ∼ 2n 0 , so c 2 s does not increase as drastically as in the leftmost and middle panels; with soft EoS at low density, the 2M constraint disfavors first order phase transitions and favors c 2 s 1/3 in some domain; there must be a domain in which c 2 s has a peak that relaxes to the conformal limit, 1/3, at very large density [104][105][106]. With this picture and the finite size of baryons, the smooth connection of 2-5n 0 domain leads to the picture of quark-hadron continuity [105][106][107] (whose original proposal comes from the indistinguishability of the hadronic superfluid matter and color-superconducting (CSC) quark matter based on the order parameters [108]). Some crossover EoS have been constructed by interpolating hadronic EoS with gapless quark matter [105] or CSC quark matter [109].
The above three cases remain all compatible with the current observations. What is common for all these cases is the rapid growth in c 2 s . It is challenging to construct a theory which leads to such a growing behavior; we note that rapid reduction of c 2 s can be easily described by the second order-like phase transition, but the rapid enhancement cannot; this can be seen from the expression (1) With n B and baryon chemical potential µ B being continuous at a transition point, the phase with the bigger susceptibility (or smaller c 2 s ) is favored. Thus we need mechanisms other than phase transitions with discontinuity.
One way to cause rapid stiffening is to phenomenologically assume the rapidly growing strong repulsion among baryons. More microscopic arguments were proposed in [110] based on the quarkyonic matter hypothesis [111] in which high density matter has the quark Fermi sea but baryonic Fermi surface. Baryons on top of the Fermi sea may be regarded as relativistic baryons with threequarks collectively moving in the same direction, unlike baryons at low density where the moving directions of three quarks are opposite one another, leaving a small baryon momentum (and pressure) but a large mass den- sity [112]. Changes from the non-relativistic to relativistic regime may be a microscopic origin of phenomenological repulsion used in hadronic models. The discussion was initially given for two-flavor matter, but recent extension included the strangeness, with the concept of shortrange correlation to make hyperons relativistic [113].

V. SUMMARY
The progress in NS observations is rapidly improving the constraints on the QCD matter. But even more exciting discoveries would come in next 10 years, thanks to on-going and forthcoming observational programs. Systematic studies are being performed in numerical simulations for NS mergers and SNs.
The 2M constraint, together with low density constraints from nuclear physics and NS radii, suggests that the core density of NS should reach 5n 0 . Considering the radii of baryons of 0.5-0.8 fm, the dynamics at quark level should be discussed near the NS core.
The natures of hadron-quark phase transitions have not been determined. But the strength of possible first order phase transitions is being constrained. If we find a heavier and more compact NS, the constraint becomes tighter. On the other hand, if we establish kinks in M -R curves and associated phenomena, it would readily confirm the existence of the first order phase transition.
The physics of hadron-quark transitions in dense QCD, which are supposed to occur around 2-5n 0 , are difficult to analyze theoretically because of the fuzzy d.o.f., but a possible strategy in near future is to improve the constraints at 1-2n 0 and 5-40n 0 . The latter domain would not directly show up in NS phenomenology, but still it constrains EoS in the NS domain.
The speed of sound in dense QCD is likely to be very different from conventional matters. Whatever natures of hadron-quark phase transitions, the c 2 s must either exceed the conformal value 1/3 or increase very rapidly at density 1-2n 0 . This tendency is certainly different from the hadron-QGP crossover phase transition, and from usual non-relativistic condensed matter with c 2 s 1.
In this article we could not touch thermal effects and general lepton fraction around the core region. Their impacts on the mechanical properties of NS cores may be limited (see however, Ref. [114] for large latent heat), but they certainly affect the composition around the cores. They can be imprinted, e.g., in neutrino emissions [115][116][117]. The composition and temperature effects are sensitive to the phase structures, and provide us with very useful tools to diagnose the structure of matter. These physics also have relations to the physics being explored by the low-energy heavy ion collisions which are also expected to have dramatic progress in 2020s due to the activation of new experimental programs [27,28].

ACKNOWLEDGMENTS
This work was supported by NSFC grant 11650110435.