Neutron star equation of state: QMF modeling and applications
A. Li, Z.-Y. Zhu, E.-P. Zhou, J.-M. Dong, J.-N. Hu, C.-J. Xia
NNeutron star equation of state: QMF modeling andapplications
A. Li a , Z.-Y. Zhu a,b , E.-P. Zhou c , J.-M. Dong d , J.-N. Hu e,f , C.-J. Xia g,h a Department of Astronomy, Xiamen University, Xiamen 361005, China; [email protected] b Institute for Theoretical Physics, Frankfurt am Main 60438, Germany c Max Planck Institute for Gravitational Physics (Albert Einstein Institute), Am M¨uhlenberg1, Potsdam-Golm, 14476, Germany d Institute of Modern Physics, Chinese Academy of Sciences, Lanzhou 730000, China e School of Physics, Nankai University, Tianjin 300071, China f Strangeness Nuclear Physics Laboratory, RIKEN Nishina Center, Wako, 351-0198, Japan g School of Information Science and Engineering, Zhejiang University Ningbo Institute ofTechnology, Ningbo 315100, China h Advanced Science Research Center, Japan Atomic Energy Research Institute, Tokai,Ibaraki 319-1195, Japan
Abstract
Because of the development of many-body theories of nuclear matter, the long-standing, open problem of the equation of state (EOS) of dense matter maybe understood in the near future through the confrontation of theoretical cal-culations with laboratory measurements of nuclear properties & reactions andincreasingly accurate observations in astronomy. In this review, we focus onthe following six aspects: 1) providing a survey of the quark mean-field (QMF)model, which consistently describes a nucleon and many-body nucleonic systemfrom a quark potential; 2) applying QMF to both nuclear matter and neu-tron stars; 3) extending QMF formalism to the description of hypernuclei andhyperon matter, as well as hyperon stars; 4) exploring the hadron-quark phasetransition and hybrid stars by combining the QMF model with the quark mattermodel characterized by the sound speed; 5) constraining interquark interactionsthrough both the gravitational wave signals and electromagnetic signals of bi-nary merger event GW170817; and 6) discussing further opportunities to studydense matter EOS from compact objects, such as neutron star cooling and pulsarglitches.
Preprint submitted to Elsevier July 13, 2020 a r X i v : . [ nu c l - t h ] J u l . Introduction The equation of state (EOS) of dense stellar matter is a problem for bothnuclear physics and relativistic astrophysics and has been greatly promoted bythe detection of gravitational waves from the GW170817 binary neutron star(NS) merger event [1] . Multimessenger observations of NS mergers [2] canprovide information for determining the EOS of supranuclear matter [3, 4, 5] andthat can possibly constrain the phase diagram of the quantum chromodynamics(QCD) [6, 7, 8, 9].In NSs, nuclear matter is present in beta equilibrium from very low den-sity to several times the saturation density ( ρ ≈ .
16 fm − ) and is extremelyneutron-rich [10, 11, 12, 13]. One assumes that there is one theoretical modelthat can correctly explain the nuclear matter data of different physical situa-tions obtained in both laboratory nuclear experiments [e.g., 15, 16, 17, 18, 19]and astronomical observations [e.g., 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30].However, this is a demanding task. It not only requires the theoretical models toextrapolate from lower density/temperature/isospin to unknown regions at highdensity/temperature/isospin [31, 32] but also depends on the relevant degreesof freedom of the problem, from nucleons to exotic particles [33, 34, 35, 71, 36],even dark matter particles [e.g., 37, 38].In this paper, we follow a widely used relativistic mean-field (RMF) ap-proach [39] based on an effective Lagrangian with meson fields mediating stronginteractions between quarks, which we call the quark mean-field (QMF) model [40,41]. It self-consistently relates the internal quark structure of a nucleon and ahyperon to the RMFs arising in nuclear and hyperonic matter, respectively, andhas been employed extensively in the calculations of finite (hyperon-)nuclei andinfinite dense matter [40, 42, 45, 46, 47, 48, 49, 50, 43, 44]. We focus on theEOS that have been developed so far, testing the QMF predictions concerningthe constraints from experiments. We also illustrate the developments of this arXiv page: http://blogs.cornell.edu/arxiv/2017/10/16/gw170817/ ab initio approach, such as the Brueckner theory [e.g.,51, 52], the chiral effective field theory [e.g., 53, 54], the quantum Monte Carlomethod [e.g., 55, 56], and the variational method [e.g., 57], which starts frommicroscopic nucleon-nucleon potentials explicitly including many-body forces.As a comparison, we include some results based on these ab initio many-bodyapproaches.The paper is organized as follows. In Sec. 2, we introduce QMF models byintroducing the confinement potential of the constituent quarks for a nucleon.Sec. 3 is then devoted to the NS properties based on the QMF EOSs. In Sec.4, we demonstrate how strange baryons, e.g., hyperons, are incorporated in theQMF model and discuss the hyperon puzzle with the obtained hyperon starmaximum mass. We also discuss hybrid stars and strange quark stars (QSs) byintroducing quark matter models. This is followed by the discussions of the NSbinary in Sec. 5. Other opportunities for studying EOS are given in Sec. 6,including NS cooling and pulsar glitches. Sec. 7 contains the main conclusionsand future perspectives of this review.
2. EOS models from the quark level within QMF
In 1988, Guichon [58] developed a novel model for nuclear matter to treat thechanges in the nucleon properties of nuclear matter, i.e., the European MuonCollaboration (EMC) effects. This model is similar to the RMF model, but thescalar and the vector meson fields couple not with the nucleons but directly withthe quarks. Then, the nucleon properties change according to the strengths ofthe mean fields acting on the quarks, and the nucleon is dealt with in termsof the MIT bag model [59]. The Guichon model was extended by Thomas andhis collaborators under the name of the quark-meson coupling (QMC) model.Excellent reviews on the QMC model can be found in the literature [60, 61];see also, e.g., [62, 63, 64, 65] for some of the latest improvements. Taking analternative model for the nucleon, the quark potential model [66], Toki and his3ollaborators constructed the QMF model [41]. For a more detailed comparisonof these two models, we refer to [40, 42]. Briefly, the bag model assumes thenucleon is constituted by bare quarks in the perturbative vacuum, i.e., currentquarks, with a bag constant to account for the energy difference between theperturbative vacuum and the nonperturbative vacuum, while in the potentialmodel, the nucleon is described in terms of the constituent quarks, which couplewith the mesons and gluons. We shall first introduce the potential model andthen introduce the QMF formalism.
In the MIT bag model, the quarks inside the nucleon are confined by a bag,which ensures that the quarks can only move freely and independently inside thenucleon through an infinite potential well. In the potential model, quarks areconfined by a phenomenological confinement potential, where the polynomialforms are widely used. A harmonic oscillator potential is usually adopted, withwhich the Dirac equation can be solved analytically, U ( r ) = 12 (1 + γ )( ar + V ) , (1)where the scalar-vector form of the Dirac structure is chosen for the quark con-finement potential and the parameters a and V are determined from the vacuumnucleon properties. When the effect of the nuclear medium is considered, thequark field ψ q ( (cid:126)r ) satisfies the following Dirac equation:[ γ ( (cid:15) q − g ωq ω − τ q g ρq ρ − (cid:126)γ · (cid:126)p − ( m q − g σq σ ) − U ( r )] ψ q ( (cid:126)r ) = 0 , (2)where σ , ω , and ρ are the classic meson fields. g σq , g ωq , and g ρq are the cou-pling constants of σ, ω and ρ mesons with quarks, respectively. τ q is the thirdcomponent of the isospin matrix, and m q is the constitute quark mass at approx-imately 300 MeV. The nucleon mass in the nuclear medium can be expressedas the binding energy of three quarks, defined by the zeroth-order term aftersolving the Dirac equation E N = (cid:80) q (cid:15) ∗ q . The quarks are simply confined in a4wo-body confinement potential. Three corrections are taken into account in thezeroth-order nucleon mass in the nuclear medium, including the contribution ofthe center-of-mass (c.m.) correction (cid:15) c . m . , pionic correction δM πN and gluoniccorrection (∆ E N ) g . The pion correction is generated by the chiral symmetry ofQCD theory and the gluon correction by the short-range exchange interactionof quarks. Finally, the mass of the nucleon in the nuclear medium becomes M ∗ N = E N − (cid:15) c . m . + δM πN + (∆ E N ) g . (3)The nucleon radius is written as (cid:104) r N (cid:105) = 11 (cid:15) (cid:48) q + m (cid:48) q (3 (cid:15) (cid:48) q + m (cid:48) q )( (cid:15) (cid:48) q − m (cid:48) q ) . (4)where (cid:15) (cid:48) q = (cid:15) ∗ q − V / , m (cid:48) q = m ∗ q + V /
2. The effective single quark energy isgiven by (cid:15) ∗ q = (cid:15) q − g qω ω − τ q g qρ ρ , and the effective quark mass is given by m ∗ q = m q − g σq σ . By reproducing the nucleon mass and radius ( M N , r N ) in freespace, we determine the potential parameters ( a and V ) in Eq. (1). We obtain V = − . a = 0 . − with m q = 300 MeV by fitting M N = 939 MeV and r N = 0 .
87 fm [15].
In the above section, we construct the nucleon at the quark level with theconfinement potential and the pion and gluon corrections. Next, we would liketo connect such nucleons in a nuclear medium with nuclear objects, such asnuclear matter and systems of finite nuclei. A good bridge is the RMF model atthe hadron level, which is developed based on the one-boson exchange potentialbetween two nucleons. The effective nucleon mass from the quark model isinserted into the RMF Lagrangian. The nucleon and meson fields are solvedself-consistently, and then, the properties of the nuclear many-body system areobtained. We mention here that the nucleons are treated as point-like particleseven though a quark model is used to describe the structure of the nucleon. Inmany-body calculations, the structure of the nucleon only modifies the effectivemass of a nucleon, i.e., Eq. (3). 5e consider the σ, ω and ρ mesons exchanging in the Lagrangian [42, 43, 44],and the cross-coupling from the ω meson and ρ meson is introduced to achievea reasonable slope of symmetry energy (see Sec. 2.3) [67], L = ψ (cid:0) iγ µ ∂ µ − M ∗ N − g ωN ωγ − g ρN ρτ γ (cid:1) ψ −
12 ( ∇ σ ) − m σ σ − g σ − g σ + 12 ( ∇ ρ ) + 12 m ρ ρ + 12 g ρN ρ Λ v g ωN ω + 12 ( ∇ ω ) + 12 m ω ω , (5)where g ωN and g ρN are the nucleon coupling constants for ω and ρ mesons. Fromthe simple quark counting rule, we obtain g ωN = 3 g ωq and g ρN = g ρq . The calcu-lation of the confined quarks gives the relation of the effective nucleon mass M ∗ N as a function of the σ field, g σN = − ∂M ∗ N /∂σ , which defines the σ coupling withnucleons (depending on the parameter g σq ). m σ = 510 MeV, m ω = 783 MeV,and m ρ = 770 MeV are the meson masses. In this Lagrangian, we alreadyconsider the static approximation on the mesons so that their time componentsare neglected. The spatial part of the ω meson disappears for the time rever-sal symmetry. The infinite nuclear matter has translational invariance, whichfurther removes the partial part of the coordinate space.The equations of motion of nucleons and mesons can be generated by theEuler-Lagrangian equation from the Lagrangian,( iγ µ ∂ µ − M ∗ N − g ω ωγ − g ρ ρτ γ ) ψ = 0 , (6) m σ σ + g σ + g σ = − ∂M ∗ N ∂σ (cid:104) ¯ ψψ (cid:105) , (7) m ω ω + Λ v g ωN g ρN ωρ = g ωN (cid:104) ¯ ψγ ψ (cid:105) , (8) m ρ ρ + Λ v g ρN g ωN ρω = g ρN (cid:104) ¯ ψτ γ ψ (cid:105) . (9)where ρ S = (cid:104) ¯ ψψ (cid:105) = 1 π (cid:88) i = n,p (cid:90) p iF dpp i M ∗ N (cid:112) M ∗ N + p i , (10) E iF = (cid:113) M ∗ N + ( p iF ) , (11)6 ∗ ω = m ω + Λ v g ωN g ρN ρ , m ∗ ρ = m ρ + Λ v g ρN g ωN ω . (12) p nF ( p pF ) is the Fermi momentum for a neutron (proton), ρ = (cid:104) ¯ ψγ ψ (cid:105) = ρ p + ρ n ,and ρ = (cid:104) ¯ ψτ γ ψ (cid:105) = ρ p − ρ n , which equals 0 in symmetric nuclear matter.Then, the energy density and pressure, with arbitrary isospin asymmetry β =( ρ n − ρ p ) /ρ , can be generated by the energy-momentum tensor, ε = 1 π (cid:88) i = n,p (cid:90) k iF (cid:113) k + M ∗ N k dk + 12 m σ σ + 13 g σ + 14 g σ + 12 m ω ω + 12 m ρ ρ + 32 Λ v g ρN g ωN ρ ω , (13) P = 13 π (cid:88) i = n,p (cid:90) k iF k (cid:112) k + M ∗ N dk − m σ σ − g σ − g σ + 12 m ω ω + 12 m ρ ρ + 12 Λ v g ρN g ωN ρ ω . (14)where we have written the meson field with their mean-field values denoted by σ , ω , and ρ . We subtract the nucleon mass from the energy density (Eq. (13)) to studythe binding energy per nucleon,
E/A = ε/ρ − M N . The parabolic approximationis usually applicable, and the energy per nucleon can be written as E/A ( ρ, β ) = E/A ( ρ, β = 0) + E sym ( n ) β + ... (15)and it is sufficient for performing the calculations only for symmetric nuclearmatter and pure neutron matter. E/A ( ρ, β = 0) can be expanded around thesaturation density, E/A ( ρ,
0) =
E/A ( ρ ) + 118 K ρ − ρ ρ + ... (16)where K is the incompressibility at the saturation point. The symmetry energy E sym ( ρ ) can be expressed in terms of the difference between the energies per7article of pure neutrons ( β = 1) and symmetric ( β = 0) matter, E sym ( ρ ) ≈ E/A ( ρ, − E/A ( ρ, E sym ( ρ ) can beexpanded around the saturation density ρ as follows: E sym ( ρ ) = E sym ( ρ )+ dE sym dρ ( ρ − ρ ) + 12 d E sym dρ ( ρ − ρ ) + ... (17)and the following parameters can be defined, where all have an energy dimension(MeV), E sym = E sym ( ρ ) , (18) L = 3 ρ ( dE sym dρ ) ρ , (19) K sym = 9 ρ ( d E sym dρ ) ρ . (20) E sym ( ρ ) can also be written as E sym ( ρ ) = E sym + 13 L ρ − ρ ρ + 118 K sym ( ρ − ρ ρ ) + ... (21)In laboratory experiments, the symmetry energy E sym ( ρ ) can be studiedby analyzing the neutron skin [e.g., 68], the different isovector nuclear excita-tions [e.g., 69], and the data on heavy-ion collisions such as isospin diffusion andthe isotopic distribution in multifragmentation processes [e.g., 70]. The largeamount of novel exotic nuclei produced in the laboratory and the developmentof radioactive ion beams have greatly stimulated new research projects on sym-metry energy [71, 72, 73, 74]. We mention here that in the following discussion,we only discuss up to the second expansion terms in both the binding energy(Eq. 15) and the symmetry energy (Eq. 17); see, e.g., [75, 76] for detailed dis-cussions on the higher order terms and the suitability of a nuclear EOS for upto the high density matter possible in NSs. Some of the latest constraints onhigher order terms are also discussed in, e.g., [77, 78, 79, 80]. There are six parameters ( g σq , g ωq , g ρq , g , c , Λ v ) in this Lagrangian (Eq. (1))to be determined by fitting the saturation density ρ and the corresponding val-ues at the saturation point of the binding energy E/A , the incompressibility8 able 1: Properties of nuclear matter at saturation predicted by the EOSs employed in thisstudy, in a comparison with the empirical ranges. The BCPM EoS, named after the Barcelona-Catania-Paris-Madrid energy density functional [52], is based on the microscopic Brueckner-Hartree-Fock (BHF) theory [51]. The BSk20 and BSk21 EoS belong to the family of Skyrmenuclear effective forces derived by the Brussels-Montreal group [84]. The high-density part ofthe BSk20 EoS is adjusted to fit the result of the neutron matter APR EOS [57], whereas thehigh-density part of the BSk21 EOS is adjusted to the result of the BHF calculations usingthe Argonne v18 potential plus a microscopic nucleonic three-body force. The TM1 EOS isbased on a phenomenological nuclear RMF model with the TM1 parameter set [85], as well asthe GM1 EOS, which uses a different parameter set [86]. The number density n is in fm − .The energy per baryon E/A and the compressibility K , as well as the symmetry energy E sym and its slope L at saturation, are in MeV. The empirical values are taken from [17, 18, 71, 81]. ρ E/A K E sym L EoS (fm − ) (MeV) (MeV) (MeV) (MeV)QMF 0.16 -16.00 240.00 31.00 40.0BCPM 0.16 -16.00 213.75 31.92 53.0TM1 0.145 -16.26 281.14 36.89 110.8BSk20 0.159 -16.08 241.4 30.0 37.4BSk21 0.158 -16.05 245.8 30.0 46.6APR 0.16 -16.00 247.3 33.9 53.8GM1 0.153 -16.32 299.2 32.4 93.9Empirical 0 . ± . − . ± . ±
20 31 . ± . . ± . K , the symmetry energy E sym , the symmetry energy slope L and the effective(Landau) mass M ∗ N ( ≈ . M N ). In particular, we use the most preferred valuesfor ( K, E sym , L ) as recently suggested by [71, 81], namely, K = 240 ±
20 MeV, E sym = 31 . ± . L = 58 . ± . K = 328 MeV [49], and we choose thiscase as well for a comparison. To study the effect of r N , we varied this parame-ter from the intermediate value 0 .
87 fm [15] by approximately 10% according toour model capability: r N = 0.80 fm, 0.87 fm, and 1.00 fm. This covers both ofthe most recent experimental analyses of the rms − radius of the proton chargedistribution: 0.879 ± ± a and V ) by repro-ducing ( m N , r N ) and then determine QMF many-body parameters by reproduc-ing the saturation properties of nuclear matter ( ρ , E/A, E sym , K, L, M ∗ N /M N ),which is shown in the first line of Table 1. Six EOS models from other theo-9 B . E . [ M e V ] (a) SNM
328 260240 .
80 fm0 .
87 fm1 .
00 fm (b) PNM
K [MeV] = 328 260240 ρ N /ρ P [ M e V / f m ] (c) SNM HIC 1 2 3 4 5 6 ρ N /ρ (d) PNM HIC + sym _ softHIC + sym _ stiff Figure 1: (Color online) Binding energy (B.E.) and pressure as a function of the numberdensity for symmetric nuclear matter (SNM) and pure neutron matter (PNM). The calcula-tions are performed for fixed symmetry parameters E sym = 31 and L = 60 MeV and differentcases of incompressibility K at saturation: K = 240 , ,
328 MeV. The results with differentnucleon radii of 0.80, 0.87, and 1.00 fm, chosen from the CODATA values and two recentexperiments [15, 82, 83], are shown by the solid, dashed, and dotted curves, respectively.Heavy-ion collisions (HIC) are expected to go through a quarkCgluon plasma (QGP) phase,where matter is strongly interacting, resulting in the development of collective motion. TheEOS results for SNM and PNM lie inside the boundaries obtained from the analysis of thecollective flow in HIC [16], which are shown with two density-dependent cases of symmetryenergy (light blue for the stiff case and dark blue for the soft case). The radius of the nucleonis shown to have limited effects on the nuclear matter EOSs even at high density. Takenfrom Zhu & Li [43]. retical frameworks are also listed, together with the empirical ranges in the lastrow.The binding energy and pressure from the QMF are displayed in Fig. 1 forsymmetric nuclear matter and pure neutron matter with different nucleon radii.The EOS results within the QMF fulfill the flow constraints from heavy-ioncollisions for both symmetric nuclear matter and pure neutron matter. Thenucleon radius has a weak effect on the nuclear matter even at high density.We address other important aspects before closing this section: • T emperature : The above discussions are only for the zero-temperaturecase, below ∼ ∼
50 MeV. Althoughthe matter is expected to cool down on timescales of 10 − − − secondsand 1 −
10 seconds, respectively, the thermal effects cannot be ignored,especially in the study of dynamic processes [e.g., 87, 88]. However, for theequilibrium configurations of cold NSs, the EOSs are not affected muchby finite temperature. For example, the temperature influence on themaximum mass is very limited, and there is an increase in the NS radiusfor a fixed amount of gravitational mass [e.g., 89]. • M eson - coupling parameters : The present calculations are structured tobe renormalizable to fix the coupling constants and the mass parametersby the empirical properties of nuclear matter at saturation. They canalso be determined by fitting the ground-state properties of closed-shellnuclei. In the latter case, a substantial stiff EOS with an extremely highincompressibility is usually obtained, ∼
328 MeV, which is not consistentwith recent experimental results [81] (as seen in Table 1). An alternativelylow compressibility usually cannot describe the finite nuclei with a properspin-orbit coupling. • Beyond mean - f ield : As a starting point, we choose the mean-field ap-proximation, which should be reasonably good at very high densities (a fewtimes the nuclear matter density). There have been studies that demon-strate that the isoscalar Fock terms could be important for the predictionof NS properties (see, e.g., Zhu et al. [90] for a study based on relativisticHartree-Fock theory). In such models, the Lorentz covariant structure iskept in full rigor, which guarantees all well-conserved relativistic symme-tries. Additionally, the attractive Fock term introduced in the frameworkof QMC could effectively decrease the incompressibility at the saturationpoint [91]. 11 . Neutron star NSs with typical masses M ≈ − M (cid:12) (where M (cid:12) is the mass of the sun, M (cid:12) = 1 . × g) and radii on the order of R ≈
10 km have many extremefeatures that are unique in the universe [10, 11, 12, 13] and lie outside the realmof terrestrial laboratories, such as rapid rotation, extremely strong magneticfields, superstrong gravitation, interior superfluidity and superconductivity, andsuperprecise spin period. These intriguing features have aroused much interestfrom researchers of many branches of contemporary physics as well as astron-omy because of their importance to fundamental physics. However, informationregarding the NS interior has not yet been sufficiently revealed through the cur-rent observations due to the complexity of the NS system and many uncertainfactors [92]. It is time to combine the efforts from different communities anddiscuss mutual interests and problems. In this section, we introduce the basicinsights into NSs, in particular the global properties such as the mass, radius,and tidal deformability of the star, which have a one-to-one correspondence toits underlying EOS and are usually used as a tool to connect nuclear physics toastrophysics for the study of dense matter above the nuclear saturation density.[e.g., 44, 80, 64, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107,108, 109, 110, 111].A wide range of matter density from ∼ . − in the star atmosphere tovalues larger than ∼ g cm − in the star core is encountered in these objects.Theoretically, the global properties are studied by using the overall EOSs asbasic input and ignoring their thin atmosphere ( ∼ . −
10) cm, where hot X-rays originate. The observations of massive NSs [21, 22, 23, 24, 25] have alreadyruled out soft EOSs that cannot reach 2 M (cid:12) . Here, this serves as a criterion forthe selection of the NS (core) EOSs. The saturation properties of the employedcore EOSs are collected in Table 1, with the empirical ranges listed in the lastrow. The determination of the EOS above the saturation density representsone of the main problems in NS study because first principle QCD calculationsare difficult to perform in such a many-body system. In most of the model12alculations available in the literature, a central density as high as (2 − ρ isfound for the maximum mass, and one or more types of strangeness-driven phasetransitions (hyperons, kaons, Delta isobars or quarks) may take place in the NSs’innermost parts, e.g., [47, 89, 90, 112, 113, 114, 115, 116, 158, 117, 118]. NSswith exotic phases are discussed in Sec. 4. In this section, we restrict ourselvesto normal nuclear matter. In the outer crust, at densities below ∼ g cm − , nuclei arrange themselvesin a Coulomb lattice mainly populated by Fe nuclei. At higher densities (10 g cm − − × g cm − ), the nuclei are stabilized against beta decay bythe filled Fermi sea of electrons and become increasingly neutron-rich. Thecomposition of the outer crust is mainly determined by the nuclear masses,which are experimentally measured close to stability, whereas the masses of thevery neutron-rich nuclei are not known, and they have to be calculated usingnuclear models.The inner crust is a nonuniform system of more exotic neutron-rich nuclei,degenerate electrons, and superfluid neutrons. The density range extends from ∼ × g cm − to the nuclear saturation density 2 . × g cm − , atwhich point the nuclei begin to dissolve and merge together. Nonsphericalnuclear structures, generically known as nuclear ”pasta”, may appear at thebottom layers of the inner crust. In fact, one of NSs’ irregular behaviors, the glitch , is closely related to the inner crust EOS and the crust-core transitionproperties [e.g., 119, 120, 121, 122]. The crust is also crucial for NS cooling[119].It may be necessary to calculate all EOS segments (outer crust, inner crust,and liquid core) using the same nuclear interaction, the so-called ”unified”EOS [e.g., 123, 84, 52, 124], since matching problems in nonunified EOS couldcause nontrivial conflicts in the predictions of the stars’ properties [125]. Fig. 2shows the crust EOS for the different theoretical approaches in Table 1. Weobserve that all outer crust EOSs display a similar pattern, with some differ-13 igure 2: Various EOSs for the outer crust (left) and inner crust (right). Among them, BCPM,TM1, BSk20, and BSk21 are unified NS EOSs, namely, all EOS segments (outer crust, innercrust, liquid core) are calculated using the same nuclear interaction. The BPS (NV) EOSfor the outer (inner) crust part is indicated by the black dotted line. The BPS outer crustEOS is based on a semi-empirical mass formula for matter from 10 g cm − to 3.4 × gcm − [126], whereas the NV inner crust EOS is based on quantal Hartree-Fock calculationsfor spherical Wigner-Seitz cells [127]. ences around the densities where the composition changes from one nucleus tothe next one. Only the TM1 EOS, based on an RMF model, shows a slightlydifferent trend due to the semiclassic-type mass calculations, in which A and Z vary in a continuous way, without jumps at the densities associated with achange in the nucleus in the crust. A is the number of nucleons in the nucleus,and Z is the atomic number. On the other hand, the energy in the inner crustis largely determined by the properties of the neutron gas; hence, the neutronmatter EOS plays an important role. Moreover, the treatment of complicatednuclear shapes, in a range of average baryon densities between the crust and thecore, produces some uncertainties in the EOS of the inner crust, where somedifferences are visible.In Fig. 3, we show the above discussed EOSs, with the full symbols indicatingthe transition point from the inner crust to the core for each chosen EOS. TheAPR and GM1 EOSs have to be matched with an inner crust EOS, which isat variance with the unified EOSs (BCPM, TM1, BSk20, and BSk21), and weachieve this by imposing that the pressure is an increasing function of the energydensity. It is evident that the matching of the GM1 core (dotted black line) to14 igure 3: Various EOSs for the NS core from the low-density inner crust indicated withsymbols. In addition to the models in Fig. 2, we include another model within the RMF,the GM1 EOS [86], as well as the present QMF model. The inner crust EOS of NV is alsoincluded. the TM1 crust (solid gray line) shows nonsmooth behavior in the dP/dρ slopecompared to the matching to the BCPM and NV crust. Since the crust effectswere shown to be more important for distorted fast-rotating stars than for staticstars [128], later in Sec. 3.4 on rotating NSs, we discuss three widely used crustEOSs (TM1, BCPM, NV + BPS) that are matched with one core EOS (GM1).Note that the above crust is based on the ground state approximation forzero-temperature matter, which can only be applied to an isolated NS born ina core-collapse supernova explosion. It is assumed that during the process ofcooling and crystallization, the plasma maintains nuclear equilibrium. Conse-quently, when the matter becomes strongly degenerate, the structure and EOSof the crust can be approximated well via cold-catalyzed matter. For an NScrust formed by accreted plasma from the companion star in a low-mass X-raybinary, the outermost layer of the accreted plasma undergoes thermonuclearflashes, observed as X-ray bursts, during the active stages. The layers deeperthan a few meters are at T < × K, becoming increasingly neutron-richdue to electron capture and neutrino emissions and finally dissolving in the15iquid core. After the fully accreted crust is formed, the layered structure ofthe crust ceases to evolve and becomes quasistationary, with matter elementsmoving inwards due to compression and undergoing exothermic nuclear trans-formations [129]. There is a microscopic model for a fully accreted crust [130]that calculates the EOS and distribution of deep crustal heating sources by fol-lowing the nuclear evolution of an element of matter consisting initially of X-rayashes under quasistatic compression from 10 g cm − to 10 g cm − (crust-coreinterface). To study the structure of NSs, we have to calculate the composition andEOS of cold, neutrino-free, catalyzed matter. We require that the NS containscharge-neutral matter consisting of neutrons, protons, and leptons ( e − , µ − )in beta equilibrium. Additionally, since we are looking at NSs after neutrinoshave escaped, we set the neutrino chemical potentials equal to zero. The en-ergy density of NS matter can be written as a function of the different partialdensities, ε ( ρ n , ρ p , ρ e , ρ µ ) = ρM N + ρE ( ρ n , ρ p ) /A + ρ µ m µ + 12 m µ (3 π ρ µ ) / π + (3 π ρ e ) / π (22)where we use ultrarelativistic and nonrelativistic approximations for the elec-trons and muons, respectively, from textbooks [10]. Then, the various chemicalpotentials µ i of the species ( i = n, p, e, µ ) can be computed, µ i = ∂ε/∂ρ i , (23)which fulfills beta-equilibrium, µ i = b i µ n − q i µ e (24)( b i and q i denote the baryon number and charge of species i ). Supplementedwith the charge neutrality condition, (cid:88) i ρ i q i = 0 (25)16he equilibrium composition ρ i ( ρ ) can be determined at the given baryon density ρ , and finally, the EOS is P ( ρ ) = ρ d ( ε/ρ ) dρ = ρ dεdρ − ε = ρµ n − ε (26)for the interior of NSs.The NS stable configuration in hydrostatic equilibrium can be obtained bysolving the Tolman-Oppenheimer-Volkoff (TOV) equation [131, 132] for thepressure P and the enclosed mass m , dP ( r ) dr = − Gm ( r ) ε ( r ) r (cid:104) P ( r ) ε ( r ) (cid:105)(cid:104) πr P ( r ) m ( r ) (cid:105) − Gm ( r ) r , (27) dm ( r ) dr = 4 πr ε ( r ) . (28) G is the gravitational constant. Starting with a central mass density ε ( r = 0) = ε c , we integrate out until the pressure on the surface equals that correspondingto the density of iron. This gives the stellar radius R , and the gravitationalmass is then m ( R ) = 4 π (cid:90) R drr ε ( r ) (29)For the description of the NS crust, we usually join the EOS P ( ε ) with theNV EOSs of Negele and Vautherin in the medium-density regime [127] andthose of Baym-Pethick-Sutherland for the outer crust [126]. After solving theTOV equations, we can obtain the maximum mass M TOV and the mass-radiusrelation for comparison with astrophysical observations.
Currently, the EOS of SNM ( β = 0) is constrained relatively well. Mat-ter with nonzero isospin asymmetry remains unknown, largely due to the un-certainty in the symmetry energy. Conflicts remain for the symmetry energy(especially its slope) despite significant progress in constraining the symmetryenergy around and below the nuclear matter saturation density [72, 73]. Thesymmetry energy slope characterizes the density dependence of the symmetry17
00 1000 E [MeV / fm ] P [ M e V / f m ] L =20 MeV L =25 MeV L =30 MeV L =35 MeV L =40 MeV L =60 MeV L =80 MeV R [km] . . . . . . M [ M (cid:12) ] PSR J0348+0432PSR J1614 − M max < . M (cid:12) Λ <
400 Λ <
800 Λ > B H P < ∞ C a u s a li t y L = 40 MeV L = 60 MeV L = 80 MeV Figure 4: (Left) NS EOSs and (right) mass-radius relation within QMF with different valuesof symmetry energy slope L , with more L cases shown in the left panel than in the right panel.The shaded region is the favored region from ab initio calculations at the subsaturation densityin chiral effective field theory [53] and from [95]. They are causal and fulfill the two-solar-mass constraint of heavy pulsars ( M TOV > M (cid:12) ) and the tidal deformability constraintof binary merger event GW170817 (Λ . ≤ . M (cid:12) star. Also shown are thelatest NICER measurements from the pulse-profile modeling of the accretion hot spots of theisolated millisecond pulsar PSR J0030+0451 [27, 28]. The general constraints from the blackhole limit, the Buchdahl limit and the causality limit are also included. The figure showsthat the radius sensitively depends on the symmetry energy slope with the maximum massonly slightly modified. A smaller L (softer symmetry energy) leads to a smaller radius. Allcases of L = 20 −
80 MeV lie within the Λ ≤
800 boundary [1] and fulfill the updated limitΛ . = 190 +390 − [180] using the PhenomPNRT waveform model at a 90% confidence level. L in the range of 20 −
80 MeV. The QMFparameters are fitted to reproduce the saturation properties in Table 1, withthe other five parameters ( ρ , E/A, K, E sym ) unchanged [44]. The TOV massof the star hardly changes with changing L and fulfills the recent observationalconstraints of three massive pulsars for which the masses are precisely mea-sured [21, 22, 23, 24, 25]. There is a strong positive correlation between theslope parameter and the radius of a 1 . M (cid:12) star [for more discussion, see, e.g.,96, 133]. However, a small dependence is found in Hornick et al. [133]. Thecases of L ≈ −
60 MeV in our QMF model may be more compatible with theneutron matter constraint [53]. Capano et al. [134] found that the radius of a1 . M (cid:12) NS is R . = 11 . +0 . − . km (90% credible interval), assuming a descrip-tion in terms of nuclear degrees of freedom remains valid up to 2 ρ . The recentNICER measurements of PSR J0030+0451 [27, 28] might indicate L (cid:38) . ≤
800 for a 1 . M (cid:12) star [1].The tidal deformability describes the magnitude of the induced mass quadrupolemoment when reacting to a certain external tidal field. It is zero in the blackhole case. The dimensionless tidal deformability Λ is related to the compactness M/R and the Love number k through Λ = k ( M/R ) − (see more discussionlater in Sec. 5.1).To study the effects of the symmetry energy slope L in more detail, wepresent the resulting Love numbers (tidal deformabilities) as a function of themass and the compactness for different L in Fig. 5 (Fig. 6). In Fig. 5, k first increases and then decreases with mass and compactness. In Fig. 6, Λmonotonously decreases with the mass and compactness. The increase in k . . . . . . M [ M (cid:12) ] . . . . . . . . . k L = 20MeV L = 40MeV L = 60MeV L = 80MeV 0 .
05 0 .
10 0 .
15 0 .
20 0 . M/R . . . . . . . k L = 20MeV L = 40MeV L = 60MeV L = 80MeV Figure 5: Love numbers as a function of the mass (left) and compactness (right) for four EOSswith different values for the symmetry energy slope L (20, 40, 60, 80 MeV). k first increasesand then decreases with mass and compactness. The increase in k (below ∼ M (cid:12) ) is due tolarge radii and a large portion of soft crust matter. The vertical line and colorful dots indicate M = 1 . M (cid:12) . Taken from Zhu et al. [44]. . . . . M [ M (cid:12) ] Λ L = 20MeV L = 40MeV L = 60MeV L = 80MeV 0 .
05 0 .
10 0 .
15 0 .
20 0 . M/R Λ L = 20MeV L = 40MeV L = 60MeV L = 80MeV Figure 6: Same as Fig. 5 but for the tidal deformabilities. Λ monotonously decreases withthe mass and compactness. Similar to Fig .5, the large values of Λ for small masses (below ∼ M (cid:12) ) are due to large radii and a large portion of soft crust matter. Taken from Zhuet al. [44]. ∼ M (cid:12) ) are due to large radiiand a large portion of soft crust matter. If no crust is considered (e.g., an EOSdescribed by a pure polytropic function), k still decreases monotonously withmass and compactness. Further loud gravitational-wave measurements frommerging binary NSs would provide data with good precision for learning moreabout the slope parameter as well as the NS structure. Moreover, the final fateof the merger, i.e., prompt or delayed collapse to a black hole or a single NS star,depends on the EOS, as well as the amount of ejected matter that undergoesnucleosynthesis of heavy elements. These discussions are presented in Sec. 5. NSs are usually rotating, and the rotational periods P of rapidly rotatingNSs (pulsars) could provide restrictions on the EOSs and their evolution pro-cesses when combined with the mass constraint. When rapidly rotating, an NSis flattened by the centrifugal force, and the TOV equation, suitable for a staticand spherically symmetric situation, cannot correctly describe the rotating stel-lar structure. We assume NSs are steadily rotating and have an axisymmetricstructure. Based on the axial symmetry, the space-time metric used to model arotating star can be expressed as ds = − e ν dt + e α dr + e α r dθ + e β r sin θ ( dφ − ωdt ) , (30)where ν, α, β and ω is the function of r, θ . The matter inside the star is approx-imated by a perfect fluid, and the energy-momentum tensor is given by T µν = ( ε + p ) u µ u ν − pg µν , (31)where ε, p and u µ are the energy density, pressure and four-velocity, respectively.To solve Einsteins field equation for potentials ν, α, β and ω , Komatsu et al.[137] transformed the Einstein equation from differential equations to integralsby using the Green function method. In this form, the asymptotic flatnesscondition, which is the boundary condition of the Einstein equation, can be21 igure 7: NSs’ masses as a function of central energy density (left) and radius (right) forthree cases of crust EOSs (TM1, BCPM, NV + BPS) matching one GM1 core EOS, with thedetailed EOS matching data shown in Table 1. The calculations are performed for both thestatic case and Keplerian rotating case. The maximum masses and radii, as well as the centraldensities, hardly depend on how the inner crusts are described for NSs heavier than 1 . M (cid:12) . satisfied automatically. This method for solving the Einstein equation is writtenas a standard code. This is the well-tested RNS code . Using tabulated EOSs,the stationary and equilibrium sequences of rapidly rotating, relativistic starscan then be computed in general relativity [see more detail about the code in,e.g., 137, 138, 139].The Keplerian (mass-shedding) frequency f K is one of the most studiedphysical quantities for rotating stars. An EOS that predicts Kepler frequenciesthat are smaller than the observed rotational frequencies is to be rejected, asit is not compatible with observation. An empirical formula was proposed inLattimer & Prakash [11], f K = f (cid:18) MM (cid:12) (cid:19) (cid:18) R
10 km (cid:19) − , (32)where M is the gravitational mass of the Keplerian configuration, R is the M , and f is a constant thatdoes not depend on the EOS. An optimal prefactor f = 1080 Hz was foundin [140, 141] for NSs as well as hybrid stars. See more discussion in [140, 141]regarding the justification of the functional form of Eq. (32) and its valid range.The calculated highest spin frequencies f K are all higher than 1000 Hz, whilethe current observed maximum is f = 716 Hz [20] for PSR J1748-2446a in theglobular cluster Terzan 5. A possible reason for this discrepancy is that the starfluid is suffering from r -mode instability [142]. A simple estimation showed thatthis would lower the maximum frequency by ∼ . M (cid:12) .Generally, rotation increases both the gravitational mass and the radius.Based on the EOSs collected in Table 1, rotation can increase the star’s gravi-tational mass up to ∼ − ∼ . M (cid:12) inthe APR case. Additionally, the star becomes flattened, and the correspondingcircumferential radius is increased up to ∼ − ∼ − M (cid:12) , the radius increase is more pronounced, reaching ∼ − ∼ − ∼ − ρ to ∼ − ρ , which is due to the effect of the centrifugalforce, effectively stiffening the EOS. We show in Fig. 8 the gravitational massas a function of the central density at various fixed rotation frequencies basedon the QMF EOS. In static stars, QMF18 gives M TOV = 2 . M (cid:12) at a central23
10 15 ε [ × g cm − ] . . . . . . . . . M [ M (cid:12) ] Figure 8: NSs’ masses as a function of the central energy density with the QMF EOS atvarious fixed rotation frequencies 1 /P =0.4-1.8 kHz. The lower blue curve is the static case,and the upper red curve corresponds to the Keplerian frequencies at different rotating cases.The change in the maximum mass M crit with frequency is indicated with a dashed blackcurve. In static stars, QMF gives M TOV = 2 . M (cid:12) at a central density ρ c = 6 . ρ with acorresponding radius 10 . f K = 1699 Hz, the maximum mass andcorresponding radius with QMF are 2 . M (cid:12) and 14 . ρ c = 8 . ρ .Curves with a fixed baryonic mass of M b = 2 . , . , . , . M (cid:12) are also shown with nearlyhorizontal gray curves. density ρ c = 6 . ρ with a corresponding radius of 10 . f K = 1699 Hz and with QMF18, the maximum mass is 2 . M (cid:12) with acorresponding radius of 14 . ρ c = 8 . ρ .One of the most interesting rotating stars is the so-called “supramassive”star, which exists only by virtue of rotation. It is well known that the onset ofthe instability of the static sequence is determined by the condition dM/dρ c = 0,i.e., the curve should stop at the maximum value of gravitational mass M .In the rotating case, the above criteria have to be generalized, i.e., a stellarconfiguration is stable if its mass M increases with increasing central densityfor a fixed angular momentum J . Therefore, the onset of the instability, whichis called the secular axisymmetric instability, is expressed by (cid:18) ∂M∂ρ c (cid:19) J = 0 . (33)Since rotation increases the mass M that a star of a given central density cansupport, the static configuration with the baryon mass M b > M TOVb (the baryon24
10 15 ε [ × g cm − ] . . . . . . . . . M [ M (cid:12) ] J=1.0J=1.5J=2.0J=2.5J=3.0J=3.5J=4.0
10 11 12 13 14 15 16 17 R [km] . . . . . . . . . M [ M (cid:12) ] J=0.5J=1.0J=1.5J=2.0J=2.5J=3.0J=3.5
Figure 9: (Left) NS masses as a function of the central energy density and (right) mass-radius relations with the QMF EOS at various fixed angular momenta J . The lower bluecurve is the static case, and the upper red curve corresponds to the Keplerian frequencies indifferent rotating cases. In the left panel, the change in the maximum mass M crit with theangular momentum is indicated with a dashed black curve. A spin-down star, losing angularmomentum over its evolution, follows the lines with fixed baryonic mass M b , shown by thenearly horizontal gray curves for M b = 2 . , . , . , . M (cid:12) . mass of a TOV mass star) does not exist. Such sequences are supramassive starsthat are doomed to collapse as they lose energy and angular momentum duringtheir spin-down, following the nearly horizontal line of fixed baryonic mass M b .We show in Fig. 9 the NS mass as a function of the central energy densityas well as the mass-radius relations with the QMF EOS at various fixed angularmomenta J . The lower blue curve is the static case, and the upper red curvecorresponds to the Keplerian frequencies in different rotating cases. In the leftpanel, the change in the maximum mass M crit with the angular momentum isindicated with a dashed black curve. There may be a universal relation between M crit /M TOV and j/j K [e.g., 143, 144, 145, 146, 147, 148] that does not dependon the specific choice of EOS or the f K value, M crit M TOV = 1 + a (cid:18) jj K (cid:19) + a (cid:18) jj K (cid:19) (34)where j = J/M is the dimensionless angular momentum and the coefficientsare a = 1 . × − and a = 7 . × − [143].Note that the above discussions focus only on the case of rigid rotation,25hile differential rotation can be much more efficient in increasing the maximumallowed mass. In differentially rotating stars, the high-density inner part mayrotate faster than the low-density outer part, so the inner part can be supportedby rapid rotation without the equator having to exceed the Keplerian limit [e.g.,149]. While both rigid and differential rotation follow axisymmetry, there arecases when a rotating NS breaks its axial symmetry if the rotational kineticenergy to gravitational binding energy ratio, T / | W | , exceeds a critical value.The abovementioned r -mode instability could also trigger NSs’ motion with off-axis symmetry. It is presently unclear whether such configurations of NSs canactually be realized in practice [e.g., 150]. Overall, it is especially importantto calculate models of rotating stars to better understand the observations ofbinary merger events (see details in Sec. 5.2).
4. EOS with exotic particles
While around the saturation densities ρ = ρ , the matter inside an NS con-sists only of nucleons and leptons, at higher densities, several other species ofparticles may appear due to the fast increase in the baryon chemical poten-tials with density [33, 35, 36], just because their appearance is able to lowerthe ground state energy of the dense nuclear matter phase. Among these newparticles are strange baryons, namely, the Λ , Σ , ± , Ξ , − hyperons. Other species(such as kaons and Delta isobars) might also appear in stellar matter, which thepresent paper does not cover. Generally, the presence of one species of strangeparticle is found to push the onset of other species of strange particles to higherdensities, even out of the physically relevant density regime [e.g., 65, 112].It is necessary to generalize the QMF study of the nuclear EOS with theinclusion of hyperons [e.g., 45, 46, 47, 49, 50, 48]. The density thresholds ofhyperons are essentially determined by the masses and their interaction. Themass of Λ( uds ) is 1116 MeV. The masses of Σ + ( uus ), Σ ( uds ), Σ − ( dds ) are1189, 1193, and 1197 MeV, respectively. The masses of Ξ ( uss ) and Ξ − ( dss )26re 1315 and 1321 MeV, respectively. From hypernuclei experiments in the lab-oratory [19], we know that Λ-nucleus and ΛΛ interactions are attractive, whileΣ-nucleus interactions are repulsive. Additionally, the nature of the Ξ-nucleusinteraction has been suggested to be attractive [151]. Theoretically, any effectivemany-body theories should respect the available hypernuclei data before pro-ceeding with other sophisticated studies [35, 36]. The adopted hyperon-mesoncouplings need to at least reproduce unambiguous hypernuclear data.At the mean-field level, the single Λ, Σ, Ξ potential well depths in symmetricnuclear matter are U ( N )Λ , Σ , Ξ ∼ −
30, 30, −
14 MeV at the saturation density, re-spectively. In the extended QMF model [46, 47], we introduce different confiningstrengths for the s quarks and the u, d quarks in the corresponding Dirac equa-tions (under the influence of the meson mean fields). The confining strength ofthe u, d quarks is constrained by finite nuclei properties, and that of the s quarksis constrained by the well-established empirical value of U ( N )Λ ∼ −
30 MeV. Themass difference among baryons is generated by taking into account the spin cor-relations E ∗ B = (cid:80) i e ∗ i + E B spin , and the spin correlations of the baryons are fixedby fitting the baryon masses in free space. In addition, the spurious c.m. motionis removed through the usual square root method as M ∗ B = (cid:112) E ∗ B − (cid:104) p (cid:105) . Thecontribution of the σ meson is contained in the effective mass M ∗ Λ , Σ , Ξ , and the ω and ρ mesons couple to the baryons with the following coupling constants: g ωN = 3 g qω , g ω Λ = cg ω Σ = 2 g qω , g ω Ξ = g qω (35) g ρN = g qρ , g ρ Λ = 0 , g ρ Σ = 2 g qρ , g ρ Ξ = g qρ (36)where a factor c is introduced before g ω Σ for a large Σ ω coupling, in contrastwith the Λ − ω coupling, to simulate the additional repulsion on the Σ − nucleonchannel, and U ( N )Σ = 30 MeV at the nuclear saturation density. The basicparameters are the quark-meson couplings ( g qσ , g qω , and g qρ ), the nonlinear self-coupling constants ( g and c ), and the mass of the σ meson ( m σ ) [for moredetail regarding the formalism and the parameters, see 40, 47]. With such aparameter set, the saturation properties of nuclear matter can be described [47].The values of a single Ξ hyperon in nuclear matter are obtained as U ( N )Ξ = − igure 10: (a) Single hyperon potential, (b) fractions of leptons and baryons, (c) EOS, and(d) mass-radius relations for NSs with hyperons within QMF. The cases without hyperons inthe star’s core are also shown in the lower panels. The maximum mass of QMF EOS withouthyperons is slightly lower than 2 M (cid:12) due to the absence of a high-order vector coupling termfor effective nuclear interaction from an earlier work [47]. When hyperons are included, themass is largely reduced and well below the observational 2-solar-mass limit. The hyperonpuzzle is also present in many microscopic studies based on developed realistic baryon-baryoninteractions [e.g., 55, 113, 153, 154]. Adapted from Hu et al. [47]. MeV, consistent with the BNL-E885 experiments [151]. The density depen-dences of the single hyperon (Λ , Σ , Ξ) potentials are shown in Fig. 10(a).Regarding the EOS of hyperonic matter, the baryon contributions can beobtained through the mean-field ansatz from the Lagrangian (including hyper-ons) [40, 47]. Electrons are again treated as a free ultrarelativistic gas, whereasthe muons are relativistic, as in Eq. (22). The total EOS can be calculated fora given composition of baryon components. This allows the determination ofthe chemical potentials of all species, which are the fundamental input for the28quations of chemical equilibrium: µ n = µ Λ = µ Σ = µ Ξ (37) µ e = µ µ (38) µ n − µ e = µ p = µ Σ + (39) µ n + µ e = µ Σ − = µ Ξ − (40)The above equations must be supplemented with two other conditions, i.e.,charge neutrality and baryon number conservation. These are ρ p + ρ Σ + = ρ e + ρ µ + ρ Σ − + ρ Ξ − , (41) ρ = ρ n + ρ p + ρ Λ + ρ Σ + ρ Σ + + ρ Σ − + ρ Ξ − + ρ Ξ . (42)Finally, the actual detailed fraction Y i = ρ i ( ρ ) of the dense matter is determinedfor each fixed baryon density ρ , as shown in Fig. 10 (b). In the low-density region(until ρ < .
21 fm − ), the proton fraction ρ p /ρ is well below 1 /
9, which fulfillsthe astrophysical observation that direct URCA cooling might not occur at tooof low densities [152]. With the properly chosen Λ , Σ and Ξ hyperon potentials,Λ is the first hyperon appearing at ∼ − ρ . Then, Ξ − hyperons appear at ∼ ρ followed by Ξ hyperons at ∼ ρ . The fractions of hyperons increase withdensity. At densities above ∼ ρ , the fractions of Λ and Ξ − are almost thesame as the fractions of protons and neutrons. Σ − , however, does not appearuntil very high density of 2 . − .In Fig. 10 (c), we show the pressure of beta-equilibrated matter as a functionof the energy density. The solid curve represents the EOS including the hyper-ons, and the dot-dashed curve is the EOS without hyperons. The EOS becomessofter with the presence of strangeness freedom. The NS properties are calcu-lated by using the EOSs with/without hyperons obtained from the QMF model,and the NS mass-radius relations are plotted in Fig. 10 (d). It is found that themaximum mass of the NSs including hyperons is approximately 1 . M (cid:12) , muchlower than that without hyperons, which is below the observational limit. Thisis the so-called hyperon puzzle, which is also found in many microscopic studiesbased on developed realistic baryon-baryon interactions [e.g., 55, 113, 153, 154].29ince hyperons are not present in nuclear matter, they cannot be determinedfrom the nuclear matter properties. Moreover, the analysis of experimentaldata on hypernuclei shows that we cannot fix these parameters in a unique way.How can a sufficiently stiff high-density EOS generate a heavy hyperon starwith properly reproduced nuclear matter properties at the saturation density?There may be three approaches forward:1. Three-body hyperon interactions can be introduced in microscopic studiesor high-order meson fields in effective calculations. If they are repulsive, astiff enough hyperon star EOS can be obtained by increasing the repulsionas the density increases, and there is a possibility of massive hyperon starswith a central density > ∼ ρ , [e.g., 91, 155, 156, 157]. This is a naturalsolution based on the known importance of three-body nucleon forces innuclear physics;2. Larger maximum masses can be produced through a transition to anotherphase of dense (quark) matter in the stellar core at high densities [e.g.,117, 118, 158]. This approach is presented in Sec. 4.2;3. A separate branch of pulsar-like objects can be introduced to accountfor the heavy ones, for example, QSs made of free quarks [e.g., 159, 160,161]. Unlike NSs, which are bound by gravity, QSs are bound by stronginteractions; therefore, they have opposite M-R dependence. This is theso-called two-branch scenario [e.g., 162, 163], which is discussed in Sec. 4.3. The matter inside the NS core possesses densities ranging from a few times ρ to one order of magnitude higher. At such densities, the hadronic mattermight undergo a phase transition to quark matter, and a hybrid NS with aquark matter or mixed core can be formed. However, the exact value of thetransition density to quark matter is unknown and still a matter of recent debatenot only in astrophysics but also within the theory of high energy heavy-ioncollisions. Additionally, it is not obvious whether the information on the nuclearEOS from high energy heavy-ion collisions can be related to the physics of NS30nteriors. The possible quark-gluon plasma produced in heavy-ion collisionsis expected to be characterized by low baryon density and high temperature,while the possible quark phase in NSs appears at high baryon density and lowtemperature. Nevertheless, we must be careful that the transition cannot occurat too low of density (below the nuclear saturation density ρ ) to maintainconsistency with the current experimental data of heavy-ion collisions.The possibility of the existence of strange quark matter in NS high-densitycores is of special interest in the present era of GW astronomy [e.g., 9, 101,164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174]. Presently, we have nounified models to address the hadron phase and the quark phase, and it is stillnot clear whether the change in the hadron phase corresponding to that in thequark phase is a crossover or a first-order transition. Here, we analyze a specificexample in the context of a first-order transition (at transition pressure P trans )to express the experimental constraints in model-independent terms. For thehadronic sector, we use the above QMF model. For the high-density quarkphase, we utilize the CSS parametrization [175], exploiting the fact that for aconsiderable class of microscopic quark matter models, the speed of sound isweakly density-dependent, e.g., [176, 177, 178]. The present scheme can onlydiscuss the transition that occurs at a sharp interface (Maxwell construction)between bulk hadronic matter and quark matter, i.e., the quark-hadron surfacetension is high enough to disfavor mixed phases (Gibbs construction). It hasbeen shown that a strong first-order phase transition with a sharp interface isthe most promising scenario to be tested or distinguished from pure hadronicmatter by future observations [9, 101, 171]. We tend to find that the observationof a two-solar-mass star and the accurate measurement of the typical NS radiusconstrain the CSS parameters, including the squared speed of sound in the high-density phase c , the hadron-quark phase transition density ρ trans , and thediscontinuity in the energy density at the transition ∆ ε/ε trans , where ρ trans ≡ ρ NM ( P trans ) , ε trans ≡ ε NM ( P trans ). 31 igure 11: EOSs (left) and mass-radius relations (right) for hybrid stars (colorful curves) at afixed discontinuity in energy density at the transition ∆ ε/ε trans = 0 . ρ trans /ρ = 1 . , . , , , , ,
6, with the symmetry energy slope 60 MeV. Thesquared speed of sound is fixed at c = 1 in the quark matter. The corresponding NSresults within the QMF are shown for comparison (black curves). The shaded region is thefavored region from the maximal model [54], with the underlying EOSs constrained at lowdensities from EFT, facilitating the complete allowed parameter space for the speed of soundabove the saturation density n and enforcing the LIGO/Virgo constraint from GW170817(70 ≤ ˜Λ ≤ For a given nuclear matter EOS ε NM ( P ), the full CSS EOS is ε ( P ) = ε NM ( P ) P < P trans ε NM ( P trans ) + ∆ ε + c − ( P − P trans ) P > P trans
We perform the calculation by varying c from the causality limit ( c =1) to the conformal limit ( c = 1 /
3, the value for systems with conformalsymmetry that may be applicable to relativistic quarks). It is worth mentioningthat perturbative QCD calculations exhibit quark matter with c values ofapproximately 0.2 to 0.3 [179]. We use units where (cid:126) = c = 1. In Fig 11,we show representative EOSs P ( ε ) for dense matter with a sharp first-orderphase transition and the corresponding mass-radius relations. The two-solar-mass lower limit for maximum gravitational mass is explicitly indicated in themass vs. radius plot. We include the curves with increasing transition densityfrom 1 . ρ to 6 ρ at a fixed energy density discontinuity and speed of soundin quark matter, and the nuclear matter EOS is chosen to be the QMF modelresult with the symmetry energy slope L = 60 MeV. We mention here that L ∼ −
60 MeV is the preferred range within QMF as indicated by the ab32nitio calculations (shown in Fig. 4). We see that a lower transition density(pressure), therefore a stiffer EOS, leads to a heavier hybrid star. The smallesthybrid star is typically the heaviest.Systematically, we carry out calculations for the mass-radius of hybrid starsspanning the whole parameter space of the speed of sound in a domain witha transition density up to ρ trans = 6 ρ and an energy density discontinuity upto ∆ ε = 1 . ε trans . The calculations are performed using two values for thesymmetry energy slope parameter of L = 30 MeV and L = 60 MeV, and theresults are shown in Fig. 12.Figure 12 displays the correlation of the radius of a 1 . M (cid:12) hybrid star R . with the transition density ρ trans /ρ (upper panels) and the maximummass M TOV (lower panels). In general, there exists an anticorrelation between R . and ρ trans /ρ and a correlation between R . and M TOV . A conservativeupper limit of 13 . R . can be obtained with different analyses [e.g.,94, 95, 97, 125, 180, 181, 182, 183]. In the upper panel, the upper limit of 13 . R . corresponds to ≈ . ρ for L = 30 MeV and ≈ . ρ for L = 60MeV. The possible onset of a first-order phase transition below such densitiesmight be strongly disfavored. From the mass measurement of heavy pulsars, wecan put lower limits on R . by making use of the R . − M TOV correlation inthe lower panel. The two-solar-mass constraint leads to a lower limit of ≈ . R . < . M TOV < . M (cid:12) .We conclude this section by further discussing the following aspects: • Submillisecond rotation : It is commonly believed that only self-boundstable QSs may rotate rapidly with a submillisecond period [188]. How-ever, it is suggested that pulsars rotating with approximately half a mil-lisecond period could also be interpreted as hybrid stars [189], with NSscontaining a metastable deconfined quark phase at their centers. This33 . 0 1 . 5 2 . 0 2 . 5 3 . 01 01 11 21 31 41 5
S o l i d : L = 6 0 M e V D a s h e d : L = 3 0 M e V
R1.4 [km] r t r a n s / r D e / e t r a n s R1.4 [km] M T O V [ M (cid:1) ]c
2Q M
Figure 12: Upper: Radius of a 1 . M (cid:12) star vs. the transition density, with the energy densitydiscontinuity explicitly indicated; Lower: Radius of a 1 . M (cid:12) star vs. the maximum mass, withthe sound speed explicitly indicated. The results are shown for two kinds of symmetry energyslopes 60 MeV (solid lines) and L = 30 MeV (dashed lines). The horizontal line represents aconservative upper limit of 13 . . M (cid:12) star with or without a phase transition [e.g.,94, 95, 97, 125, 180, 181, 182, 183]. The vertical line in the lower panel represents the two-solar-mass constraint. There are cases when no 1 . M (cid:12) star is possible, shown by breaks inthe curves. Taken from Miao et al. [184]. • M ixed phase : In the present work, we adopt the simple Maxwell con-struction. The Gibbs construction provides a realistic model of the phasetransition between the hadronic and quark phases inside the star [190],yielding a range of baryon densities where both phases coexist, which pro-vides an EOS containing a pure hadronic phase, a mixed phase, and apure QM region [e.g., 118, 191]. The pressure is the same in the twophases to ensure mechanical stability, while the chemical potentials of thedifferent species are related to each other, satisfying chemical and betastability. Both the hadron and quark phases are separately charged whilepreserving total charge neutrality [190]. As a consequence, the pressure isa monotonically increasing function of the density. The realization of themixed phase depends on the nuclear surface tension, which is currently anunknown parameter [192]. The Gibbs treatment is the zero surface tensionlimit of the calculations, including finite-size effects. It was demonstratedthat the influence of different constructs on the maximum mass value israther small [191].
We now turn to the description of the bulk properties of uniform quark mat-ter. The strange quark matter is composed of up ( u ), down ( d ) and strange ( s )quarks with charge neutrality maintained by the inclusion of electrons (hereaftermuons as well if present): 23 ρ u − ρ d − ρ s − ρ e = 0 , (43)The baryon number conservation,13 ( ρ u + ρ d + ρ s ) = ρ, (44)35s also satisfied with n being the baryon number density. Due to the weakinteractions between quarks and leptons, d → u + e + ˜ ν e , u + e → d + ν e ; s → u + e + ˜ ν e , u + e → s + ν e ; s + u ↔ d + u , The β -stable conditions µ s = µ d = µ u + µ e should be fulfilled in neutrino-freematter. The energy density and pressure include both contributions from quarksand leptons, and those of leptons can be easily calculated by the model of anideal Fermi gas such as in the NS matter case.In the density regime achieved inside compact stars, the dense matter prop-erties cannot be calculated directly from the first principle lattice QCD or per-turbative QCD. The latter is only applicable at ultrahigh densities beyond therange of compact stars. The current theoretical description of quark matteris based on phenomenological models [e.g., 159, 160] and burdened with largeuncertainties. In the following, we consider the nonperturbative contributionsfrom perturbative QCD [193]. For simplicity, we use the pQCD thermodynamicpotential density to the order of α s [194],Ω pt = Ω + Ω α s , (45)withΩ = (cid:88) i = u,d,s g i m i π (cid:26)(cid:20) (cid:18) ¯Λ m i (cid:19) + 4 (cid:21) [ u i v i − ln( u i + v i )] + 3 [ u i v i − ln( u i + v i )] − v i (cid:27) , (46)where u i ≡ µ i /m i and v i ≡ (cid:112) u i −
1. The coupling constant α s and quarkmasses m i run with the energy scale and can be determined by [194], α s (¯Λ) = 1 β L (cid:18) − β ln Lβ L (cid:19) , (47) m i (¯Λ) = ˆ m i α γ β s (cid:20) (cid:18) γ β − β γ β (cid:19) α s (cid:21) . (48)Here, L ≡ ln (cid:16) ¯Λ Λ (cid:17) , and we take the MS renormalization point Λ MS = 376 . m u = 3 . m d = 8 MeV, andˆ m s = 158 MeV. The parameters for the β -function and γ -function are β = π (11 − N f ), β = π (102 − N f ), γ = 1 /π , and γ = π ( − N f ) [196](The formulas are for arbitrary N f , and in this study, N f = 3). It is not clearhow the renormalization scale evolves with the chemical potentials of quarks,and we use ¯Λ = C (cid:80) i µ i , with C = 1 − pt + B ≡ Ω pt + B QCD + ( B − B QCD ) exp (cid:34) − (cid:18) (cid:80) i µ i − µ (cid:19) (cid:35) (49)where we take B = 40 ,
50 MeV / fm [59] for the calculations and ∆ µ = ∞ indicates no medium effect for the bag parameter. If α S and m u,d,s run with theenergy scale as reported by the particle data group [195], the maximum mass ofQSs does not reach ∼ M (cid:12) . In such cases, the dynamic rescaling of the bag con-stant with finite ∆ µ is essential, which basically originates from nonperturbativeeffects such as chiral symmetry breaking and color superconductivity [6, 7, 199]. B QCD = 400 MeV / fm is demanded by the dynamic equilibrium condition atthe critical temperature of the deconfinement phase transition and is obtainedby equating the pressures of the QGP ( − B QCD + 37 π T /
90) and pion gas( π T /
30) at T = T c ( ∼
170 MeV).At given chemical potentials µ i , the pressure P , particle number density ρ i ,and energy density ε are determined by: P = − Ω , (50) ρ i = g i π (cid:0) µ i − m i (cid:1) − ∂ Ω ∂µ i α s + ρ , (51) ε = Ω + (cid:88) i µ i ρ i . (52)37 Q M F
P (MeV/fm3) e ( M e V / f m ) c m a xQ M / c p Q C D cQM/c e ( M e V / f m ) Figure 13: EOSs (left) and sound speed c QM (right) of SQM. The EOSs are generated in theperturbation model, fulfilling the available astrophysical constraints of mass [21, 22, 23, 24, 25],radius [27, 28] and tidal deformability [29] for QSs. They are compared with the results of theperturbative QCD (red curve) without nonperturbative corrections. The EOS of the nuclearmatter obtained with the QMF (blue curve) is shown in the left panel. The maximum soundspeeds c maxQM are explicitly indicated. The horizontal lines in the c QM plot show the conformallimit. To ensure a large mass for QSs above two solar mass, the sound speed is necessarilylarge, c QM /c > .
68. Taken from [169]. ρ = − C (cid:88) i (cid:18) ∂ Ω ∂m i + ∂ Ω ∂m i α s (cid:19) d m i d¯Λ+ C ∂ Ω ∂ ¯Λ α s + C d α s d¯Λ − ∂B∂µ i . (53)In Fig. 13, we show the EOSs generated in the perturbation model, whichfulfill the available astrophysical constraints of mass [21, 22, 23, 24, 25], ra-dius [27, 28] and tidal deformability [29] for QSs. They are compared with theresults of pQCD without nonperturbative corrections, taking C = 1 − B = 0. The EOSs for the nuclear matter obtained with the QMF models arealso shown. [200] pointed out that if the two-solar-mass constraint is combinedwith the hadronic matter EOS below and around the nuclear saturation density, c QM might first increase then decrease after reaching a maximum (maybe evenup to 0 . c ) and finally approach the conformal limit c/ √ c maxQM )ranges from 0.68 c to c , similar to previous results for NSs [193, 202, 203, 204].We show in Fig. 14 the various properties of strange quark matter based onthe perturbation model using the exemplary parameters of C = 3 . , B = 40MeV/fm , and ∆ µ = 800 MeV. The composition, binding energy, pressure,sound velocity, and scaled bag parameter are shown as functions of the baryondensity or chemical potential. Note that in the binding energy plot, the con-dition that the strange quark matter be the absolute stable strong-interactionsystem, requiring that at P = 0 , E/A ≤ M ( Fe /
56) =930 MeV, is fulfilled.Fig. 15 shows the predicted properties of the QSs, including the mass, Lovenumber, and tidal deformability.Previously, although the quark star EOS models could reach 2 solar mass,they required a too small surface density (not much larger than the nuclear satu-ration density) in some cases, and a larger maximum mass meant an even smallersurface density (because of the anticorrelation between the two [205]), for ex-ample, the CDDM2 model [206] and the PMQS3 model [205]. Then, the radius39 r i/ r r ( f m - 3 ) dus 0 . 0 0 . 4 0 . 8 1 . 29 0 09 5 01 0 0 01 0 5 01 1 0 01 1 5 01 2 0 0 E/A (MeV) r ( f m - 3 ) M (
F e ) / 5 6
P (MeV) m b ( M e V ) c o n f o r m a l l i m i t cQM/c r ( f m - 3 ) B (MeV/fm3) r ( f m - 3 ) P (MeV) r ( f m - 3 ) Figure 14: Various properties of betastable strange quark matter, including the quark frac-tions, the binding energy, the pressure, the sound velocity, and the scaled bag parameter,which are plotted as a function of the baryon density or chemical potential. The calculationsare performed based on the perturbation model using the parameters C = 3 . , B = 40MeV/fm , and ∆ µ = 800 MeV. The horizontal line in the E/A plot indicates that of themost stable Fe nucleus. The horizontal line in the c QM plot shows the conformal limit. . 0 0 . 5 1 . 0 1 . 50 . 00 . 51 . 01 . 52 . 02 . 5 M (M (cid:1) ) r ( f m - 3 ) M (M (cid:1) ) R ( k m )0 . 1 5 0 . 2 0 0 . 2 5 0 . 3 0 0 . 3 50 . 0 00 . 0 50 . 1 00 . 1 50 . 2 00 . 2 5 1 . 4 M (cid:2)(cid:1) s t a r k2 M / R1 . 4 M (cid:2)(cid:1) s t a r 0 . 1 5 0 . 2 0 0 . 2 5 0 . 3 0 0 . 3 511 01 0 01 0 0 0 L M / R
Figure 15: Various properties of QSs based on the EOS in Fig. 14, including the gravitationalmass, the Love number, and the tidal deformability, which are plotted as a function of centraldensity, radius, or compactness. The crosses in the upper two panels show where the maximummass is reached. The locations of a 1 . M (cid:12) star are explicitly indicated in the lower two panels. (and the tidal deformability) exceeded the observational values [1, 27, 28, 29].There models were not welcomed by particle physicists studying hadrons (forwhich experiments have been established studying the nonperturbative effects)because in such a density realm, the quarks are thought to be very dilute andare very possibly confined inside hadrons. In the real world, we do obtain nu-clear matter rather than quark matter around the nuclear saturation density.The present perturbative model with an in-medium bag can achieve both a rea-sonable surface density and a maximum mass as large as 2 . M (cid:12) . The predictedproperties of dense matter ( c QM , Γ) and quark stars (
R, M ), as well as the EOSof Fig. 14, are collected in Table 2.We conclude this section by further discussing the following aspects: • QS vs. hybrid stars : Although it is known that the degree of freedom ishadronic around the nuclear saturation density, the QCD phase state forcold, dense matter of intermediate densities is unfortunately unknown, anda great deal of effort is being applied in the communities of astrophysics,41 able 2: QS EOS with proper sound velocity behavior and the predicted properties of densematter ( c QM , Γ) and quark stars (
R, M ). The calculations are done based on the perturbationmodel using the parameters C = 3 . , B = 40 MeV/fm , and ∆ µ = 800 MeV. See Sec. 4.3for details. ρ (fm − ) ε (MeV/fm ) P (MeV/fm ) c QM /c Γ R (km) M/M (cid:12) . M (cid:12) . The similarity of the sound speed of the hadron-quarkmixed phase with that of the pure quark matter in the intermediate densityregion of ∼ − ρ , a particular shape with a peak, further complicatesthe distinguishing of QS from hybrid stars [169]. • T wo branch scenario : Because of the tension of a low tidal deformability(190 +390 − [180]) and a high maximum mass (2 . +0 . − . M (cid:12) for the presentlyheaviest pulsar [25] and ≤ . M (cid:12) based on the numerical simulationstudies on NS binary mergers [211, 212, 213]) for a certain EOS in theNS model, binary QSs have been proposed as the possible scenario for theGW170817 event [100, 214]. A binary QS merger for some binary configu-rations could eject amounts of matter (comparable to the binary NS case)to account for the electromagnetic observations in the optical/infrared/UVbands (namely, kilonova) [215]. A magnetar with QS EOS is preferred asthe post-merger remnant for explaining some groups of short gamma-rayburst (SGRB) observations [206, 205]. It has been suggested that becauseof this discrepancy (if confirmed), small and large stars of the same masscould coexist as hadronic and quark matter stars [162, 163]. • Comments on the maximum mass of N S/QS : Various microscopic cal-culations of NS matter (without strangeness) indicate a possible upperlimit of ∼ . − . M (cid:12) for the NS maximum mass, for example, Bruecknertheory calculations [153] and quantum Monte Carlo calculations [56]. Ex-43tic particles, if they are present, usually lower the limit as a result of theextra degrees of freedom during the phase transition (while the appearanceof quarks might increase the limit in the case of crossover [104, 201]). Thequark-hadron crossover EOS gives a maximum mass of 2 . M (cid:12) [201]. Thebound of M TOV (cid:46) . − . M (cid:12) may be applicable to QSs. For example,the maximum mass of QSs is 2 . M (cid:12) (2 . M (cid:12) with color-flavor-locked su-perfluity [216]) within the MIT bag model [100]. The present perturbationmodel yields 2 . M (cid:12) with the peak sound speed ( c maxQM ) approaching thespeed of sound (we expect an increase in the value including the uncertainsuperfluity of ∼ . M (cid:12) ). These high theoretical limits on the maximummass are higher than (but close to) the observational bound of pulsars ofapproximately 2 . M (cid:12) [25]. A looser upper limit based on extreme causalEOSs may be in the range of M TOV < . − . M (cid:12) [54, 185, 186, 187].The observations of accreting black holes, on the other hand, revealeda paucity of sources with masses below 5 M (cid:12) [e.g., 217, 218, 219, 220].Presently, binary mergers involving one or two companions have massesthat fall into the so-called mass gap range (3 − M (cid:12) ) that are hard todistinguish [e.g., 221, 222, 223].
5. Neutron star binary
The gravitational waves (GWs) detected from binary neutron star (BNS)merger event GW170817 [224], as well as its electromagnetic (EM) counter-parts [225], announced the beginning of the multimessenger astronomy era.In addition to hinting at the origin of SGRB [1, 226] and revealing the siteof r-process nucleosynthesis [2, 227], our knowledge of the EOS of cold densematter at supranuclear densities has been greatly enriched. In the past year,various studies have been performed to constrain the EOS of dense matter, ei-ther by putting constraints on observable characteristics of NSs [i.e., radius ortidal deformability; see e.g., 95, 97, 180] or by connecting the constraint withmodel parameters in nuclear physics [i.e., symmetry energy slope or neutronskin parameter; see e.g., 44, 228]), which could be tested by nuclear physics44xperiments. Some studies also go beyond the conventional NS scenario andput constraints on compact star models involving strong interaction phase tran-sitions [e.g., 97, 100, 171, 181].Those constraints mainly come from 3 aspects from a BNS merger event.First, during the inspiral stage, unlike binary black hole (BBH) mergers, defor-mation is induced for each NS due to the tidal field of the companion, providingadditional dissipation of the orbital energy and angular momentum and henceaccelerating coalescence [224]. This deformation therefore leaves a detectablesignature in the GW signal of the late inspiral stage, from which we can learnabout the tidal deformability of the NS EOSs [229]. Second, the detection ofSGRB hints at a delayed collapse to a BH for the merger remnant [230, 231].This interpretation of the SGRB observation provides information on the maxi-mum mass of a nonrotating configuration for the NS EOS (namely, M TOV ). Forinstance, the EOS should not be too stiff; otherwise, the remnant supramassiveNS lives much longer [232] in the magnetar central engine model [233, 234].However, if the EOS is too soft, the merger might result in a prompt collapse toa BH. In such occasions, the magnetic field might not be enhanced sufficiently(by a differentially rotating NS remnant) and thus not able to launch a jet [212].Third, the features of the transient optical/infrared/UV observations (namely,the kilonova) powered by the radio activity of the neutron-rich elements in theejecta depend directly on the mass, velocity and electron fraction in the ejecta,which is related to the properties of the EOS for the merging NS.In this section, we briefly review the information we have learned about theEOS of NSs from the BNS merger events GW170817 and GRB170817A as wellas AT2017gfo.
The finite size effects of NSs alter the late inspiral GW signal comparedwith that of the BBH case [229, 235]. Through the leading order, the GWobservations constrained a combination of the tidal deformability for each NSin the binary (Λ and Λ ) [224]. 45Λ = 1613 (12 q + 1)Λ + (12 + q ) q Λ (1 + q ) , (54)in which q = M /M is the mass ratio of the binary. The dimensionless tidaldeformability of each star is Λ = 23 k ( RM ) , (55)where k is the tidal Love number describing the fraction between the inducedquadrupole moment of a star and the external tidal field and R and M arethe radius and mass of the star, respectively. On the other hand, k can beobtained for a given EOS for any given mass and hence can be tested with theobservation of GW170817.Practically, the tidal deformability is fitted to the GW observation togetherwith other parameters [29]. For instance, in the Taylor F2 post-Newtonianaligned-spin model, 13 parameters need to be fitted, including 7 extrinsic pa-rameters (sky location, distance of the source, polarization angle, inclinationangle and coalescence phase and time) and 6 intrinsic parameters (mass of eachNS, tidal deformability of each NS and the aligned spin of each NS). Therefore,the uncertainties in the estimation of other parameters weakly correlate withthe determination of the tidal deformability. Hence, the constraint on the tidaldeformability is normally made with certain assumptions.For instance, as seen in [224], the assumption in the spin parameter of the NScould significantly affect the interpretation of the mass of each NS, thus furtheraffecting the constraint on the tidal deformability. Through the assumption thatthe Λ of each NS vary independently, the first constraint on Λ and Λ couldbe obtained under different spin priors. The result favors the softer EOS, i.e.,the EOS that predicts more compact stars. Another analysis assumes a uniformprior for ˜Λ, which sets an upper limit of 800 on ˜Λ in the low-spin case and 700in the high spin case. Alternatively, through the expansion of Λ( M ) arounda certain M , constraints can be directly placed on the tidal deformability of acertain mass star. This constraint is Λ(1 . ≤ . ≤
800 in the low-spin case. 46ollow-up analysis further improves these constraints under more assump-tions. For example, in [180], instead of assuming an independent and uniformprior for the tidal deformability of each star, the EOS of each star in the mergingbinary is assumed to be the same. Consequently, the area of the 90% confidenceregion in the Λ -Λ parameter space shrinks by a factor of 3. This also im-proves the determination of Λ(1 .
4) to 190 +390 − . In [29], a lower cut-off frequencyof 23 Hz is used instead of the 30 Hz in the original analysis. Although thetidal effects mainly affect the GW signals above several hundred Hz, a lowerfrequency cut off allows for the better determination of other parameters, henceimproving the measurement of the tidal deformability. In [183], it was pointedout that under the assumption that two stars in a binary system have a commonEOS, there is an approximate relation between the tidal deformability of eachstar, i.e., Λ / Λ = q , where q is the mass ratio. With the aid of this relation,once the assumption in the mass ratio of the binary is made, the tidal deforma-bility can be further constrained. In [183], the improved analysis shows that ˜Λis 222 +420 − for a uniform mass prior and 245 +453 − for a mass prior inferred fromobserved double neutron star systems and 233 +448 − for a mass prior informed byall galactic neutron star masses within the 90% credibility level.We can directly test existing EOS models by simply calculating the tidaldeformability and comparing it with the observational constraint. Nevertheless,more insight could be gained regarding the EOS of neutron-rich matter by moresystematically studying the impact of the EOS (i.e., p as a function of ε , where p and ε are the pressure and rest mass density of the matter) on the tidaldeformability. Such interpretations are presented in [e.g., 95, 97, 180]. Forinstance, in [180], the logarithm of the adiabatic index of the EOSs is treatedas a polynomial of the pressure for the high density EOS, namely, Γ = Γ( p ; γ i )and γ i = ( γ , γ , γ , γ ) are free parameters. For densities below half the nuclearsaturation density, the EOS is connected to the SLy EOS [236]. The samplingof the EOS then consists of uniformly sampling γ i in certain intervals. Foreach of the EOS samples, the mass radius relation and tidal deformability couldbe theoretically obtained and constrained by the observation of both the tidal47 . . . . . . . a B / e ff [ M e V ] fl a v o r li ne fl a v o r li ne . M (cid:12) . M (cid:12) Λ ( . ) = Λ ( . ) = Λ ( . ) =
10 20 30 40 50 60 70 80 90 100 ∆ [MeV] B / e ff [ M e V ] fl a v o r li n e . M (cid:12) . M (cid:12) Λ ( . ) = Λ ( . ) = Λ ( . ) = (2) a = 0 . Figure 16: Parameter space for QS EOS models within the MIT bag model ( B eff , a , ∆) fornormal QSs (left) and superfluid QSs (right) obtained by combining the GW170817 constrainton Λ(1.4), the two-solar-mass constraint on M TOV and the stability window for quark matter.With the constraint of the ”2 flavor” line, we ensure that normal atomic nuclei do not decayinto nonstrange quark matter. With the constraint of the ”3 flavor” line, we ensure thatstrange quark matter is more stable than normal nuclear matter, namely, Bodmer-Witten’sconjecture [207, 208]. The perturbative QCD correction parameter a characterizes the degreeof the quark interaction correction, with a = 1 corresponding to no QCD corrections (Fermigas approximation). a = 0 .
61 is chosen to be close to the calculated result with differentchoices of the renormalization scale [237]. The effective bag constant ( B eff ) also includes a phe-nomenological representation of nonperturbative QCD effects. Due to the strong correlationbetween M TOV and Λ(1 . .
4) from the two-solar-masslimit, namely, ∼
510 ( ∼ deformability and 2 solar mass pulsars [21, 22, 23, 24, 25]. The constraint on theneutron star radius is R = 11 . +1 . − . km for the merging binary of GW170817.Similar analysis can be found in, e.g., [95], which shows that the radius of a1.4 solar mass NS is in the range of [9 . , .
8] km. However, it is worth notingthat such analysis might be affected by the choice of EOS priors. In [97], itwas pointed out that when the prior for possible twin star (for which there is ahadron-quark phase transition inside the star) branch EOSs is considered, theradius becomes less constrained, i.e., R . ∈ [8 . , .
74] km.In addition to systematic studies on parameterized EOS priors, phenomeno-logical models can be applied to interpret tidal deformability constraints. Ac-cording to [228], a better upper limit for neutron skin effects is obtained com-pared with that of the experiment done by PREX, and better results could beachieved with future GW observations and terrestrial nuclear physics experi-ments. Both the symmetry energy parameter and the symmetry energy slope48 igure 17: Postmerger product fractions for (a) stable star, (b) supermassive star and (c) blackhole for the NS and QS EOS models, labelled with the star type plus the corresponding M TOV :Unified BSk21 (red line labelled NS2.28), nonunified GM1 (blue line labelled NS2.37), andMIT model (green and purple lines labelled QS2.08 and QS2.48). In panel (b), the observed22% constraint for supermassive stars from Gao et al. [240] is shown by the horizontal line forcomparison. The vertical dotted line in the same panel is for a typical initial period of 1 ms.Taken from Li et al. [205]. are better constrained with respect to previous nuclear physics studies [44]. Un-der the assumption that GW170817 originates from a binary quark star (BQS)merger, the quark interaction parameters are studied in [100]. Fig. 16 shows theresults of both normal and superfluid QSs. It is worth noting that due to thefinite surface density of QSs, a surface correction needs to be taken into accountwhen calculating Λ of a QS [238]. Therefore, a QS could have a larger TOVmaximum mass without violating the tidal deformability constraint comparedwith those of NSs [239].
BNS mergers have long been proposed as the central engine of SGRBs [226].This suggestion has been verified by GRB170817A detected by Fermi/GBMand INTEGRAL/SPI-ACS, which accompanies the detection of GW170817.According to the time of the merger implied from the chirp signal, there is a49 . ± .
05 s delay for the onset of the SGRB [2]. The detection of GRB170817Anot only helps the determination of the location of the source, which allows forabundant follow-up observations in other bands, but also provides useful in-formation about the post-merger evolution of the merger event, thus providingconstraints on the EOSs [e.g., 240, 241, 242]. Exemplary fractions of the outcomeof the binary are shown in Fig. 17 using both NS and QS EOSs. The dependenceon the EOS, as well as the initial period, is evident [205]. It is found that thefraction of stable star (panel (a)) is determined by the static maximum mass M TOV . The fractions of supermassive star (panel (b)) and black hole (panel (c))are further sensitive to the initial period since the fast-rotating configurationsof the star have to be taken into account for them [206].Depending on the TOV maximum mass of the NS EOS and the total massof the merging binary, there could be 4 different outcomes after the merger: • if the total mass of the binary system ( M tot ) is much larger than M TOV ,the direct formation of a black hole (BH) on a dynamic time scale, namely,prompt collapse, occurs. The total binary mass, above which promptcollapse can occur, is denoted as the threshold mass ( M thres ); • if M tot is smaller than M thres but larger than the maximum mass that canbe reached by uniformly rotating NSs (denoted M supra ), then a short-livedNS could exist as a remnant supported by differential rotation . Thedifferential rotation dissipates due to magnetorotational instabilities aswell as viscosity within a timescale of ∼
100 ms, and then the NS collapsesto a BH; • if M tot is smaller than M supra but larger than M TOV , the remnant is along-lived supramassive NS. The uniformly rotating NS could still loseangular momentum by magnetic dipole radiation, but it takes a muchlonger time to sufficiently reduce the angular momentum to induce the Such NSs are called hypermassive NSs (HMNS). NSs that can be supported by onlyuniform rotation are called supramassive NSs (SMNS) • if M tot is smaller than M TOV , a stable NS remnant exists.The GW observations do not provide any hint to which scenario applies tothe case of GW170817 due to the lack of post-merger GW observations [224,243]. The electromagnetic counterparts, on the other hand, can indicate whathappens after the merger of the two NSs.A very robust interpretation is that scenario 1) should be excluded due to theSGRB detected. According to Ruiz et al. [212], in the BH central engine modelfor SGRBs, enhancement in the magnetic field of the merger remnant due tothe differential rotation of the hypermassive NS is essential for jet formation. Aprompt collapse results in a magnetic field that is too weak to explain the SGRBobservations. Hence, the detection of GRB170817A directly implies M tot It has long been suggested that the BNS merger is an important site for theproduction of heavy elements in the Universe [227]. R-process nucleosynthesisis expected to occur in the neutron-rich matter ejected during the merger. Theradioactive decay of such neutron-rich isotopes could power a transient in opti-cal/UV/IR, i.e., a kilonova [248, 249]. Such a transient event (AT2017gfo) wasdetected several hours after the merger time of GW170817 [2], the luminosity,spectrum and light curve of which are consistent with the prediction of the kilo-nova model. Such a kilonova detection not only enriches our knowledge aboutthe abundance of heavy elements in the Universe but also greatly increases ourunderstanding of NS EOSs.The observational properties, for example, the peak luminosity and peaktime, of the kilonova are closely related to the ejecta properties (c.f. [250]). Theabundance of lanthanides (atomic numbers from 58 to 71) is strongly related tothe electron fraction ( Y e ) of the ejecta. On the other hand, the opacity of theejecta is mainly determined by the lanthanides in it, and is hence related to theelectron fraction of the ejecta. The overall luminosity is related to the amount ofradioactive heavy elements and thus the total mass of the ejecta. The ejecta isexpanding at a certain velocity and becomes translucent after a period of time.Therefore, the characteristic duration of a kilonova is related to the velocity ofthe ejecta. To summarize, the mass, velocity and electron fraction of the ejectaare key parameters for understanding the observations of AT2017gfo.The ejecta during a BNS merger mainly consists of two components. Thefirst component is the so-called dynamic ejecta, which is normally more neutron-rich (lower Y e ). Part of the dynamic ejecta is due to the tidal torque during theinspiral [251, 252]; hence, it has a lower temperature and very low Y e (smaller53han 0 . − . 2) and is spatially distributed in the equatorial plane of the binary.Another part of the dynamic ejecta results from shock during coalescence (alsocalled shock-driven ejecta) [253, 254]. Due to the higher temperature at coales-cence, this part of the dynamic ejecta normally has a slightly higher electronfraction ( Y e > . 25) [255, 256] and can expand in the polar direction. In additionto the dynamic ejecta, the neutrino emissions from the remnant before collaps-ing to BH as well as the viscosity could further drive more ejecta (wind-drivenejecta) from the disc surrounding the remnant [257, 258, 259]. Due to neutrinoirradiation, this part of the ejecta has a broader distribution of Y e , which couldbe as high as 0.5 [260, 261, 262]. Clearly, the amount of wind-driven ejectais dependent on the lifetime of the remnant NS. For instance, in the case of aprompt collapse, the wind-driven ejecta could be significantly suppressed.The kilonova observation following GW170817 has shown clear evidence oftwo distinct ejecta components [263, 250, 264] , an early rising ( ∼ . M (cid:12) with a relativelylarger electron fraction Y e > . 25 and velocity v blue ∼ . − . c . For the redcomponent observed at later times, in total, approximately 0 . M (cid:12) lanthanide-rich ( Y e < . 25) ejecta is needed, with a lower velocity of v r ed ∼ . − . c .The inferred properties can be used to constrain the EOS of the merged NSs,although this constraint is quite model dependent.One property we can use to constrain NS EOS models is the mass of theejecta, as it is related to the properties of the merging binary. As summarized in[265], stiffer EOS models (for which the tidal deformability is larger) typicallyhave a smaller amount of tidal-induced dynamic ejecta than softer EOS models.However, softer EOSs normally eject more dynamic ejecta overall because of amore violent coalescence and hence eject more shock-driven ejecta. The amountof wind-driven ejecta depends on the lifetime of the merger remnant before Note that there are studies arguing a model with 3 components (c.f. [341]). M TOV of the NS EOS.Ideally, the details of the ejecta property for BNS mergers with different EOSmodels could be determined through extensive numerical simulations. However,to fully resolve the wind-driven ejecta, a very long-term post-merger simulationwith the implementation of viscosity and neutrino cooling is required, which isnormally extremely time consuming and not affordable in practice. Neverthe-less, conservative estimations can still be made. In [266], it has been found thatthe total dynamic ejecta plus all the mass in the disc surrounding the remnanthas a positive correlation with the ˜Λ parameter of the binary. As not all thematter in the remnant disc is ejected, the dynamic ejecta plus the remnantdisc mass has to be larger than the inferred mass of the ejecta to explain theobservation of AT2017gfo, which is approximately 0.05 solar mass. This setsup a lower limit for the dynamic ejecta mass plus the remnant disc mass andhence a lower limit for the binary tidal deformability. In [266], this conservativeconstraint is ˜Λ > < 400 could still be consistent with the luminosity of AT2017gfoin terms of the ejecta mass. Therefore, this lower limit must be considered withcaution.Another implication of the observation of AT2017gfo is the fate of the mergerremnant (the 4 scenarios mentioned in the previous subsection). This observed”blue” component of the ejecta clearly rules out the possibility of a prompt col-lapse, in which case there is a negligible amount of high Y e shock-driven ejectaand wind-driven ejecta. In such a case, the kilonova observation should be dom-inated by the tidal ejecta and thus should be red. Distinguishing whether theremnant is long-lived is very uncertain and model dependent. In [244], it wassuggested that if a long-lived SMNS is produced, then a significant amount ofthe rotational kinetic energy of the SMNS is injected into either a collimatedrelativistic jet or the ejecta within tens of seconds. This extra energy injec-tion is considered to be inconsistent with the observations of GRB170817A and55T2017gfo. Therefore, the authors of [244] believe an HMNS or very short-livedSMNS is produced in the remnant and put a similar upper limit on M TOV of2.17 M (cid:12) . Nevertheless, in [246], it was shown that with a long-lived SMNS asthe merger remnant, both the early and late emission of AT2017gfo can be ex-plained without requiring an unrealistically low opacity and high ejecta mass.Similar arguments can also be found in [245]. A long-lived SMNS remnant isbelieved to be able to provide strong neutrino emissions to reduce the lanthanidecontamination in our line of sight as well as to produce enough ejecta with rel-atively high speed in the post-merger phase. Particularly, a temporal featureobserved 155 d after the merger in the X-ray afterglow provides a more directhint supporting a long-lived remnant. Considering the possibility of differentmerger outcomes, the allowed range for the TOV maximum mass for NSs couldactually be larger [268]. We have summarized some of the many studies on NS EOSs in light ofGW170817 and its EM counterparts. Nevertheless, it is worth noting that thephase diagram of strong interactions at supranuclear densities is not yet clearlyunderstood due to the nonperturbative nature of QCD at low energy scales.Apart from conventional NS models, other models involving strong interactionphase transitions are suggested, e.g., twin stars or strange stars [175, 269, 270].As is summarized below, a binary quark star (BQS) scenario could be totallyconsistent with the observation of GW170817 and its EM counterparts. Due tothe self-bound nature, QSs have a very large surface density. This leads to manysignificant differences between QSs and NSs. For example, when supported byuniform rotation, QSs can reach an even higher maximum mass with respect totheir TOV maximum mass (40% more) than NSs (20% more) [271]. The finitesurface density leads to a correction when calculating the tidal deformability[135]. QSs could reach a much higher T / | W | ratio when rotating, which couldlead to more significant GW radiation in the post-merger phase [100]. As aresult, the above analysis on NS models should not be directly applied to QSs.56t is interesting and useful to understand the constraints on QS models fromwhat we have learned from GW170817 and its EM counterparts and how todistinguish between NS/QS models in the multimessenger era.Qualitatively, the tidal deformability measurement constrains the stiffnessof QSs, similar to the case of NSs. Stiffer EOSs normally reach higher M TOV but also have larger tidal deformability than a softer EOS due to the larger sizeof the star described by a stiffer EOS. There is an overall positive correlationbetween M TOV and Λ in NS models [95]. Investigating QS properties based onthe MIT bag model reveals a similar relation between M TOV and Λ(1 . 4) [150].However, the quantitative results are quite different. In [95], creating NS EOSmodels stiff enough to reach M TOV ∼ M (cid:12) with Λ(1 . ≤ 800 was found to beimpossible. In [239], it was shown that a self-bound strange star model can bestiff enough to reach M TOV > M (cid:12) without violating the tidal deformabilityconstraint.It was suggested that a BQS merger should result in a clean environmentwith little or no hadronic dynamic ejecta [272]. It is not easy to verify thisargument with numerical simulations, as the density discontinuity on the QSsurface is difficult to handle numerically. Nevertheless, in [215], several BQSmerger simulations were performed with the SPH method, and it was shownthat with some binary configurations, a BQS merger could eject a comparableamount of (quark) matter in the case of BNS mergers. According to [273, 274],under certain conditions (i.e., if the baryon number of an ejected quark nuggetis smaller than a critical value), quark matter could efficiently evaporate intonormal nucleon matter and contribute to the kilonova observation [214, 275].In addition, as uniformly rotating QSs can reach a higher maximum mass, thepost-merger remnant is likely to be longer lived than the case of the BNS merger.It has been found that a magnetar with QS EOS as the post-merger remnantis better for understanding the internal X-ray plateau observations followingSGRBs [206]. In addition, both differentially rotating and uniformly rotatingtriaxial QSs are found to be sufficient GW emitters [100, 149], which could betested by future GW observations. 57 10 11 12 13 14 15 16 R [km] . . . . . . . . M [ M (cid:12) ] APRSFH0DD2ALF2 Figure 18: The combined constraints on NS EOSs from the multimessenger observationsof GW170817, GRB170817A and AT2017gfo. The gray shaded regions are excluded by thetidal deformability measurement of GW170817 [243, 95]. The red shaded regions are forbiddenbecause the GRB170817A and AT2017gfo observations exclude a prompt collapse after merger[94]. The solid horizontal line is the 2.01 solar mass lower limit for M TOV according to theobservation of the massive pulsar [22], and the dotted horizontal line is the 2.17 solar massconstraint. Note that if the SGRB is powered by a BH central engine (or magnetar centralengine), the dotted horizontal line is an upper limit (or lower limit). The M-R relation ofseveral commonly used NS EOSs is shown in the figure. Another interesting possibility is a BNS merger that leads to a conversionof the merger remnant to a QS. In such a case, the inspiral GW signal anddynamic ejecta properties should be exactly the same as the case of a normalBNS merger, whereas the post-merger behavior could be quite different. On theone hand, if the phase transition occurs partially inside the star (i.e., only thehigh-density core part of the remnant), a softening of the EOS occurs, hencereducing the lifetime of the merger remnant as well as shifting the f peak in thepost-merger GW signal to a higher frequency compared with those of the purelynucleonic merger remnant case [166, 168]. On the other hand, if the entire starcould be converted to a QS after merger, which results in a stiffening of the EOS,the lifetime of the remnant is longer, and the f frequency is smaller [276]. Inthis scenario, the time delay between the merger and the SGRB is believed tobe the time needed for converting the surface of the remnant to quark matter,which significantly reduces the baryon load in the environment, thus helping58he formation of a collimated jet [277].Overall, the possibilities of a BQS merger or BNS merger with a QS remnantare consistent with the multimessenger observations of GW170817. However,current knowledge of the details of mergers involving QSs is quite limited. Moresimulations need to be performed for a better and more complete understandingin the future. With the help of numerical results and more future observations,whether QS could be involved or formed in a merger event could be distin-guished, particularly according to the post-merger GW signals. To summarize, the multimessenger observation of BNS merger GW170817has greatly increased our knowledge about the EOS of dense matter. The mostrobust constraint is from the tidal deformability encoded in the GW signalduring the inspiral. Such a tidal deformability measurement could translateinto a constraint on the radius of the NS at a given mass. The EM counterpartscontain large amounts of information on the EOS models, the most reliable ofwhich is to exclude the prompt collapse scenario. This provides independentconstraints on the neutron star radius for a given mass. Other constraints onthe lifetime of the remnant NS, however, depend on the central engine modelof the SGRB and are not reliable. The constraints could be totally opposite inthe BH central engine model and magnetar model, and the current observationscould not reliably rule out either possibility for GW170817. We summarizedall the constraints mentioned above, which is shown in Fig. 18. Nevertheless,GW170817 is just the beginning of the multimessenger era. As an increasingnumber of GW signals and EM counterparts from BNS mergers are detected inthe future, our knowledge of NS and even QS EOS models will be enriched. Inparticular, if the time delay between the merger and the collapse to BH couldbe more robustly determined by future observations (either by post-merger GWsignals or by more extensive EM counterpart observations), the TOV maximummass can be crucially inferred, leading to a much better understanding of theEOS of dense neutron-rich matter. 59 . Other opportunities from compact objects With the ever-increasing accuracy of observational instruments, more detailsof the signals emitted by NSs can be quantitatively monitored. Apart fromthe measurements of NS masses, radii and tidal deformabilities, high-densityNS models can be confronted with the surface temperatures of isolated NSsof known or estimated ages and the thermal photon luminosity of the X-raytransients in quiescence with an estimated time-averaged accretion rate on theNS [92, 152, 279, 278, 280]. The NS EOS determines the stellar structure, aswell as the effective masses and superfluid gaps of baryons, and is thereforecrucial for the heat capacity and neutrino emission rate [e.g., 152, 281, 282].In the Newtonian framework, the energy balance equation for NS cooling iswritten as [278], dE th dt = C v dTdt = − L ν − L γ + H, (58)where T and C v are the stellar internal temperature and the total heat capac-ity, respectively. The loss of the thermal energy E th occurs through neutrinoemission (total luminosity L ν ) and photon emission (total luminosity L γ ). H represents all possible energy sources to heat the star, such as the decay of themagnetic field energy stored in stars. Current simulations of thermal evolutionare usually based on a general relativistic formulation, and some robust codeshave already been established , which comprises all the relevant cooling reac-tions: direct URCA (DU) ( n → p + l + ν l , p + l → n + ν l ), modified URCA(MU) ( n + N → p + N + l + ν l , p + N + l → n + N + ν l ), nucleon-nucleonbremsstrahlung (NNB) ( N + N → N + N + ν + ν ), Cooper pair breaking andformation (PBF) processes ( N + N → [ N N ] + ν + ν ). Exemplary cooling curvesof a 1 . M (cid:12) NS [283] are displayed in Fig. 19.In the quiescent state of X-ray transients, the accreted matter sinks graduallyin the interior of the NS and undergoes a series of nuclear reactions [129]. These igure 19: Cooling curves of a canonically isolated NS within the minimal coolingparadigm [338], without including fast neutrino emissions, charged meson condensate, hy-perons, or confinement quarks in canonical NSs. The stellar structure is built with the APREOS [57]. The calculations are carried out in three cases for a comparison: without any Z factors (the Fermi surface depletion due to the SRC), with Z factors only in superfluidity, with Z factors both in superfluidity and neutrino emission and with Z factors. Dong et al. [283]found that the SRC effect reduces the neutrino emissivity for the DU, MU, NNB and PBFprocesses, as well as the heat capacity of the stellar interior [283]. reactions release some heat, which propagates into the whole NS, inwardly heat-ing the core and outwardly emitted in the form of photons at the surface. This isthe so-called deep crustal heating. The heating curves of X-ray transients can bederived, relating the L γ in quiescence to the estimated time-averaged accretionrate ˙ M [284]. The relevance of the pasta phase, which is beyond the neutrondrip density, to explaining some X-ray transients (if confirmed) might be re-garded as smoking-gun evidence of the NS model for pulsar-like objects [e.g.,285] and disfavors the alternative QS model.Above all, a reliable theory for NS cooling, in combination with accurateobservations, is indispensable for gaining important information about the stel-lar interior. The complexity of NS systems, such as anisotropic magnetic fieldsand the compositions of the stellar core and envelope, is not controlled welltheoretically, currently rendering the task of providing reliable and quantitativepredictions extremely difficult, and considerable effort might be achieved in thefuture. 61 .2. Pulsar glitch and glitch crisis A glitch (sudden spin-up) is a common phenomenon in pulsar observationsand was discovered during pulsar timing studies in the Vela pulsar [287]. Sincethen, the number of known glitches has greatly increased, with more than 555glitches now known in more than 190 pulsars. The observed glitches are collectedby the Jodrell Bank Observatory and the ATNF Pulsar Catalogue . Duringglitches, the pulsar spin frequency Ω suddenly increases over a very short timeand then usually relaxes to its preglitch rate over a longer time. The glitchsize, often defined as the relative increase in the spin frequencies during glitches,∆Ω g / Ω, has a bimodal distribution ranging from ∼ − to 10 − with peaks at ∼ − and 10 − [288, 289]. The glitches in young pulsars, including magnetars,are generally large [290]. However, the youngest pulsars, e.g., the Crab pulsar,tend to have more frequent and smaller glitches. Various authors have usedobserved glitch properties as a probe to investigate the pulsar inner structure,i.e., the EOS of dense matter [291].The physical mechanism behind glitches is suggested to be a sudden de-crease in the NS moment of inertia, which could result from the coupling anddecoupling between the (observed) charged component (rotating slower; labelledas index p ) and the superfluid component (rotating faster; labelled as index n ) [292]. The fractional moment of intertia I n /I p is related to the glitch activ-ity A g = (1 /T )( (cid:80) ∆Ω p ) / Ω p as follows,2 τ c A g (cid:46) I n I p , (59)where T is the total data span for glitch monitoring and (cid:80) (∆Ω g ) is the sumof all observed glitch frequency jumps. τ c = − Ω p / p is the characteristicage of the pulsar. The glitch activity A g can be estimated for systems thathave exhibited at least two glitches of similar magnitude, like the Vela pul-sar. The glitch observations from the Vela pulsar place a constraint on the An “antiglitch”, i.e., an abrupt spin-down, has also been detected [286]. igure 20: Fractional moments of inertia as a function of the stellar mass for both the chargedcomponent I p and the (crustal) superfluid component I n , with four cases of NS EOSs (BCPM,BSk21, BSk20, Av18*). Av18* indicates the nonunifed EOS of ”Av18 + NV + BPS”. Takenfrom Li & Wang [331]. fractional moment of inertia, , which is I n /I p (cid:38) . 6% [293, 294]. It has beenargued that “entrainment” of the neutron superfluid by the crystalline struc-ture of the crust greatly reduces its mobility, increasing the lower limit from1 . 6% to ∼ 7% and making it very difficult for the nuclear EOSs to fulfill witha normal M > . M (cid:12) NS [122, 294, 295]. This is clearly shown in Fig. 20.Consequently, the unpinning of the crustal superfluid is insufficient to accountfor large glitches. This is the so-called glitch crisis problem. Other mechanisms,e.g., the unpinning of core superfluid neutrons, may be required. However, themobility of superfluid neturons are related to the effective neutron mass, whichhas been discussed actively in the literature, see e.g., [296]. According to thecalculation in Watanabe & Pethick [296], the constraint for the fractional mo-ment of inertia is I c /I ≥ . − . 3% Then, an NS of M (cid:46) . M (cid:12) NS would beacceptable, and there is no need to invoke the core superfluid. However, this isan open problem, and more detailed work has to be done.At present, we are still far from a thorough understanding of the general63icture of glitches; for example, is there a connection between the stellar interiorand the magnetosphere of a star? How can various types of post-glitch behaviorbe explained? Is there an alternative model besides the superfluid model? Theoriginal starquake model [297] suggested that the change in the moment ofinertia was due to relaxation in the NS (solid) crust to the current equilibriumoblateness but has difficulty explaining the large glitches observed from theVela pulsar. In a solid star model, the whole body of the star, rather thanonly the crust, is in a solid state. In such cases, the glitch amplitude could beexplained [298, 299]. It is a challenge to quantitatively describe glitch behaviorssince the related physical processes are complicated. There has been progressin the determination of NS properties in the literature [e.g., 300, 301, 302, 303]. 7. Summary and perspectives Although NSs were anticipated in the early thirties and discovered as pul-sars in the late sixties, the state of their liquid interiors remains unclear dueto a lack of a good understanding of QCD at low energy scales. The currentand upcoming multimessenger observatories [e.g., 304, 305, 306, 307, 308, 309]will continue improving the detection of pulsars together with the precise mea-surements of their masses and radii. Laboratory experiments will provide anemerging understanding of nuclear matter EOS and the transition to deconfinedquark matter. Hopefully, the high density behavior of the NS EOS can be de-termined soon, shedding light on many unsolved problems in nuclear physicsand high-energy astrophysics [e.g., 310, 311, 312, 313, 314, 315].Since the compact star EOS is such a demanding problem, it is necessaryto combine efforts from different communities and discuss mutual interests andproblems [e.g., 134, 316, 317]. Additionally, it is important to establish newquantitative results testable by experiments/observations. In this work, the mi-croscopic physics of dense matter are modelled within the QMF, which connectsthe structure of a nucleon to the EOS of infinite nuclear matter, with a widerange of experimental and observational data available for use. The parameterspace of the QMF has been constrained well for describing NSs, following the64resent robust measurements of mass, radius, and tidal deformability. The pureNS maximum mass is approximately 2 . M (cid:12) , with a satisfying reproduction ofthe nuclear matter properties around the saturation density. The results have amodest dependence on the model parameters. Based on the available hyperonu-clei data, the hyperon puzzle is present, and we need to understand better howhyperon three-body interaction plays a role to understand more clearly whetherhyperons are relevant in NSs (especially the heavy ones). The CSS parametriza-tion allows us to handle the high-density cores of NSs in a model-independentway. After demonstrating how the NSs’ mass and radius depend on the CSSparameters for the phase transition of deconfined quark matter, we find a safeupper limit for the hybrid star maximum mass at approximately 3 . M (cid:12) basedon the extreme causality EOS, similar to previous studies. In particular, theNS/QS maximum mass is predicted to be approximately 2 . − . M (cid:12) from vari-ous model calculations, as well as analysis on the merger remnant of GW170817.Therefore, if more massive pulsars above 2 . M (cid:12) or fewer massive black holes be-low 5 M (cid:12) are found, what is their nature? Is there a mass gap between NSs andblack holes and why? Such problems are mysteries to be solved. The detailedfeature of the sound speed in quark matter is explored in a perturbative model,and an enhancement in the sound speed is necessary to fulfill the two-solar-massconstraint of pulsars, located at intermediate densities, indicating that the pairof quarks starts to play a nontrivial role at such densities. Quark superfluidsshould not be ignored in the quark matter relevant to compact stars. This isalso consistent with a lower bound for tidal deformability that is more consistentwith the GW170817-like event in another MIT model than in the case withoutsuperfluid. In the future, even if we understand the stiffness of the EOS, afurther challenge is the particle degree of freedom in cold, dense matter. Thiscould help us understand the nonperturbative properties of low-energy QCD (orthe parameters of an effective model). This is work for the future.65 cknowledgements We are grateful to the members of the XMU neutron star group. We arethankful to the anonymous referee for his or her beneficial comments and HongShen, Jirina Stone, Anthony Thomas for insightful discussions. This work wassupported by the National Natural Science Foundation of China (Grant Nos.11873040, 11775276, 11705163, and 11775119), the Youth Innovation PromotionAssociation of the Chinese Academy of Sciences, the Natural Science Foundationof Tianjin, the China Scholarship Council (Grant No. 201906205013), and theNingbo Natural Science Foundation (Grant No. 2019A610066). References [1] Abbott, B. P. et al., Oct. 2017. GW170817: Observation of GravitationalWaves from a Binary Neutron Star Inspiral. PhRvL, 119 (16), 161101.[2] Abbott, B. P. et al., Oct. 2017. Multi-messenger Observations of a BinaryNeutron Star Merger. ApJ, 848 (2), L12.[3] Abbott, B. P. et al., Feb. 2017. Exploring the sensitivity of next generationgravitational wave detectors. CQGra, 34 (4), 044001.[4] Baiotti, L., Nov. 2019. Gravitational waves from neutron star mergers andtheir relation to the nuclear equation of state. PrPNP, 109, 103714.[5] Chirenti, C., Gold, R., Miller, M. C., Mar. 2017. Gravitational Wavesfrom F-modes Excited by the Inspiral of Highly Eccentric Neutron StarBinaries. ApJ, 837 (1), 67.[6] Alford, M., Schmitt, A., Rajagopal, K., Sch¨afer, T., Oct. 2008. Colorsuperconductivity in dense quark matter. RvMP, 80 (4), 1455–1515.[7] Baym, G., Hatsuda, T., Kojo, T., Powell, P. D., Song, Y., Takatsuka, T.,May 2018. From hadrons to quarks in neutron stars: a review. RPPh,81 (5), 056902. 668] Paschalidis, V., Yagi, K., Alvarez-Castillo, D., Blaschke, D. B., Sedrakian,A., Apr. 2018. Implications from GW170817 and I-Love-Q relations forrelativistic hybrid stars. PhRvD, 97 (8), 084038.[9] Han, S., Mamun, M. A. A., Lalit, S., Constantinou, C., Prakash, M., Nov.2019. Treating quarks within neutron stars. PhRvD, 100 (10), 103022.[10] Shapiro, S. L., Teukolsky, S. A., 1983. Black holes, white dwarfs, andneutron stars : the physics of compact objects. bhwd.book[11] Lattimer, J. M., Prakash, M., Apr. 2004. The Physics of Neutron Stars.Sci, 304 (5670), 536–542.[12] Baldo, M., Burgio, G. F., Feb. 2012. Properties of the nuclear medium.RPPh, 75 (2), 026301.[13] Graber, V., Andersson, N., Hogg, M., Jan. 2017. Neutron stars in thelaboratory. IJMPD, 26 (8), 1730015–347.[14] Guichon, P. A. M., Saito, K., Rodionov, E., Thomas, A. W., Feb. 1996.The role of nucleon structure in finite nuclei. NuPhA, 601 (3), 349–379.[15] Mohr, P. J., Newell, D. B., Taylor, B. N., Sep. 2016. CODATA rec-ommended values of the fundamental physical constants: 2014*. RvMP,88 (3), 035009.[16] Danielewicz, P., Lacey, R., Lynch, W. G., Nov. 2002. Determination ofthe Equation of State of Dense Matter. Sci, 298 (5598), 1592–1596.[17] Wang, M., Audi, G., A. H., W., F. G., K., MacCormick, M., Xu, X.,Pfeiffer, B., Dec. 2012. The Ame2012 atomic mass evaluation. ChPhC,36 (12), 003.[18] Angeli, I., Marinova, K. P., Jan. 2013. Table of experimental nuclearground state charge radii: An update. ADNDT, 99 (1), 69–95.6719] Feliciello, A., Nagae, T., Sep. 2015. Experimental review of hypernu-clear physics: recent achievements and future perspectives. RPPh, 78 (9),096301.[20] Hessels, J. W. T., Ransom, S. M., Stairs, I. H., Freire, P. C. C., Kaspi,V. M., Camilo, F., Mar. 2006. A Radio Pulsar Spinning at 716 Hz. Sci,311 (5769), 1901–1904.[21] Demorest, P. B., Pennucci, T., Ransom, S. M., Roberts, M. S. E., Hes-sels, J. W. T., Oct. 2010. A two-solar-mass neutron star measured usingShapiro delay. Natur, 467 (7319), 1081–1083.[22] Antoniadis, J., et al., Apr. 2013. A Massive Pulsar in a Compact Rela-tivistic Binary. Sci, 340 (6131), 448.[23] Fonseca, E., et al., Dec. 2016. The NANOGrav Nine-year Data Set: Massand Geometric Measurements of Binary Millisecond Pulsars. ApJ, 832 (2),167.[24] Arzoumanian, Z. et al., Apr. 2018. The NANOGrav 11-year Data Set:High-precision Timing of 45 Millisecond Pulsars. ApJS, 235 (2), 37.[25] Cromartie, H. T., et al., Jan. 2020. Relativistic Shapiro delay measure-ments of an extremely massive millisecond pulsar. NatAs, 4, 72–76.[26] ¨Ozel, F., Freire, P., Sep. 2016. Masses, Radii, and the Equation of Stateof Neutron Stars. ARA&A, 54, 401–440.[27] Miller, M. C., et al., Dec. 2019. PSR J0030+0451 Mass and Radius fromNICER Data and Implications for the Properties of Neutron Star Matter.ApJ, 887 (1), L24.[28] Riley, T. E., et al., Dec. 2019. A NICER View of PSR J0030+0451: Mil-lisecond Pulsar Parameter Estimation. ApJ, 887 (1), L21.[29] Abbott, B. P. et al., Jan. 2019. Properties of the Binary Neutron StarMerger GW170817. PhRvX, 9 (1), 011001.6830] Abbott, B. P. et al., Feb. 2020. Model a comparison from LIGO-Virgodata on GW170817’s binary components and consequences for the mergerremnant. CQGra, 37 (4), 045006.[31] Prakash, M., Lattimer, J. M., Pons, J. A., Steiner, A. W., Reddy, S.,2001. Evolution of a Neutron Star from Its Birth to Old Age. LNP, 578,364.[32] Burgio, G. F., Fantina, A. F., 2018. Nuclear Equation of State for CompactStars and Supernovae. ASSL, 457, 255.[33] Weber, F., Mar. 2005. Strange quark matter and compact stars. PrPNP,54 (1), 193–288.[34] Li, A., Burgio, G. F., Lombardo, U., Peng, G. X., Jan. 2008. Exotic Phasesin Neutron Stars. IJMPE, 17 (9), 1635–1647.[35] Gal, A., Hungerford, E. V., Millener, D. J., Aug. 2016. Strangeness innuclear physics. RvMP, 88 (3), 035004.[36] Tolos, L., Fabbietti, L., Feb. 2020. Strangeness in Nuclei and NeutronStars. preprint (arXiv:2002.09223).[37] Li, A., Huang, F., Xu, R.-X., Sep. 2012. Too massive neutron stars: Therole of dark matter? APh, 37, 70–74.[38] Das, H. C., Kumar, A., Kumar, B., Biswal, S. K., Nakatsukasa, T., Li, A.,Patra, S. K., Feb. 2020. Effects of Dark Matter on Nuclear and NeutronStar Matter. preprint (arXiv:2002.00594).[39] Nikˇsi´c, T., Vretenar, D., Ring, P., Jul. 2011. Relativistic nuclear energydensity functionals: Mean-field and beyond. PrPNP, 66 (3), 519–548.[40] Shen, H., Toki, H., Apr. 2000. Quark mean field model for nuclear matterand finite nuclei. PhRvC, 61 (4), 045205.6941] Toki, H., Meyer, U., Faessler, A., Brockmann, R., Dec. 1998. Quark meanfield model for nucleons in nuclei. PhRvC, 58 (6), 3749–3752.[42] Zhu, Z.-Y., Li, A., Hu, J.-N., Shen, H., Feb. 2019. Quark mean-field modelfor nuclear matter with or without bag. PhRvC, 99 (2), 025804.[43] Zhu, Z.-Y., Li, A., Mar. 2018. Effects of the nucleon radius on neutronstars in a quark mean field model. PhRvC, 97 (3), 035805.[44] Zhu, Z.-Y., Zhou, E.-P., Li, A., Aug. 2018. Neutron Star Equation of Statefrom the Quark Level in Light of GW170817. ApJ, 862 (2), 98.[45] Shen, H., Toki, H., Sep. 2002. Study of Λ hypernuclei in the quark mean-field model. NuPhA, 707 (3), 469–476.[46] Hu, J. N., Li, A., Shen, H., Toki, H., Jan. 2014. Quark mean-field modelfor single and double Λ and Ξ hypernuclei. PTEP, 2014 (1), 013D02.[47] Hu, J. N., Li, A., Toki, H., Zuo, W., Feb. 2014. Extended quark mean-fieldmodel for neutron stars. PhRvC, 89 (2), 025802.[48] Hu, J., Shen, H., Nov. 2017. Single-Ξ − hypernuclei within a quark mean-field model. PhRvC, 96 (5), 054304.[49] Xing, X., Hu, J., Shen, H., Oct. 2016. Quark mean field model with pionand gluon corrections. PhRvC, 94 (4), 044308.[50] Xing, X., Hu, J., Shen, H., May 2017. Quark mean field model with pionand gluon corrections for Λ and Ξ hypernuclei and neutron stars. PhRvC,95 (5), 054310.[51] Baldo, M., Nov. 1999. Nuclear Methods And The Nuclear Equation OfState. in: Baldo M. (Ed.), Series: International Review of Nuclear Physics.World Scientific Publishing Co. Pte. Ltd., vol. 8.[52] Sharma, B. K., Centelles, M., Vi˜nas, X., Baldo, M., Burgio, G. F., Dec.2015. Unified equation of state for neutron stars on a microscopic basis.A&A, 584, A103. 7053] Hebeler, K., Lattimer, J. M., Pethick, C. J., Schwenk, A., Aug. 2013.Equation of State and Neutron Star Properties Constrained by NuclearPhysics and Observation. ApJ, 773 (1), 11.[54] Tews, I., Margueron, J., Reddy, S., Jun. 2019. Confronting gravitational-wave observations with modern nuclear physics constraints. EPJA, 55 (6),97.[55] Lonardoni, D., Lovato, A., Gandolfi, S., Pederiva, F., Mar. 2015. HyperonPuzzle: Hints from Quantum Monte Carlo Calculations. PhRvL, 114 (9),092301.[56] Gandolfi, S., Lonardoni, D., Lovato, A., Piarulli, M., Jan. 2020. Atomicnuclei from quantum Monte Carlo calculations with chiral EFT interac-tions. preprint (arXiv:2001.01374).[57] Akmal, A., Pandharipande, V. R., Ravenhall, D. G., Sep. 1998. Equationof state of nucleon matter and neutron star structure. PhRvC, 58 (3),1804–1828.[58] Guichon, P. A. M., Jan. 1988. A possible quark mechanism for the satu-ration of nuclear matter. PhLB, 200 (3), 235–240.[59] Degrand, T., Jaffe, R. L., Johnson, K., Kiskis, J., Oct. 1975. Masses andother parameters of the light hadrons. PhRvD, 12 (7), 2060–2076.[60] Guichon, P. A. M., Stone, J. R., Thomas, A. W., May 2018. Quark-Meson-Coupling (QMC) model for finite nuclei, nuclear matter and be-yond. PrPNP, 100, 262–297.[61] Saito, K., Tsushima, K., Thomas, A. W., Jan. 2007. Nucleon and hadronstructure changes in the nuclear medium and the impact on observables.PrPNP, 58 (1), 1–167.[62] Bohr, H., Moszkowski, S. A., Panda, P. K., Providˆencia, C., daProvidˆencia, J., Feb. 2016. QMC approach based on the Bogoliubov in-71ependent quark model of the nucleon. International Journal of ModernPhysics E 25 (2), 1650007.[63] Stone, J. R., Guichon, P. A. M., Reinhard, P. G., Thomas, A. W., Mar.2016. Finite Nuclei in the Quark-Meson Coupling Model. PhRvL, 116 (9),092501.[64] Motta, T. F., Kalaitzis, A. M., Anti´c, S., Guichon, P. A. M., Stone, J. R.,Thomas, A. W., Jun. 2019. Isovector Effects in Neutron Stars, Radii, andthe GW170817 Constraint. ApJ, 878 (2), 159.[65] Motta, T. F., Thomas, A. W., Guichon, P. A. M., Mar. 2020. Do Deltabaryons play a role in neutron stars?. PhLB, 802, 135266.[66] Isgur, N., Karl, G., Dec. 1978. P-wave baryons in the quark model.PhRvD, 18 (11), 4187–4205.[67] Horowitz, C. J., Piekarewicz, J., Jun. 2001. Neutron Star Structure andthe Neutron Radius of Pb. PhRvL, 86 (25), 5647–5650.[68] Tagami, S., Yasutake, N., Fukuda, M., Yahiro, M., Mar. 2020. Determina-tion of symmetry energy from experimental and observational constraints;prediction on CREX. preprint (arXiv:2003.06168).[69] Danielewicz, P., Lee, J., Feb. 2014. Symmetry energy II: Isobaric analogstates. NuPhA, 922, 1–70.[70] Li, B.-A., Chen, L.-W., Ko, C. M., Aug. 2008. Recent progress and newchallenges in isospin physics with heavy-ion reactions. PhR, 464 (4-6),113–281.[71] Oertel, M., Hempel, M., Kl¨ahn, T., Typel, S., Jan. 2017. Equations ofstate for supernovae and compact stars. RvMP, 89 (1), 015007.[72] Li, B.-A., Ramos, `A., Verde, G., Vida˜na, I., Feb. 2014. Topical issue onnuclear symmetry energy. EPJA, 50,7273] Baldo, M., Burgio, G. F., Nov. 2016. The nuclear symmetry energy.Progress in Particle and Nuclear Physics 91, 203–258.[74] Zhang, Y., Liu, M., Xia, C.-J., Li, Z., Biswal, S. K., Feb. 2020. Constraintson the symmetry energy and its associated parameters from nuclei toneutron stars. preprint(arXiv:2002.10884).[75] Dutra, M., Louren¸co, O., S´a Martins, J. S., Delfino, A., Stone, J. R.,Stevenson, P. D., Mar. 2012. Skyrme interaction and nuclear matter con-straints. PhRvC, 85 (3), 035201.[76] Dutra, M., Louren¸co, O., Avancini, S. S., Carlson, B. V., Delfino, A.,Menezes, D. P., Providˆencia, C., Typel, S., Stone, J. R., Nov. 2014.Relativistic mean-field hadronic models under nuclear matter constraints.PhRvC, 90 (5), 055203.[77] Malik, T., Alam, N., Fortin, M., Providˆencia, C., Agrawal, B. K., Jha,T. K., Kumar, B., Patra, S. K., Sep. 2018. GW170817: Constraining thenuclear matter equation of state from the neutron star tidal deformability.PhRvC, 98 (3), 035804.[78] Zhang, N.-B., Li, B.-A., Jul. 2019. Implications of the Mass M= { } { -0.10 } ˆ { +0.11 } M (cid:12) of PSR J0740+6620 on the Equation of State of Super-dense Neutron-rich Nuclear Matter. ApJ, 879 (2), 99.[79] Zhang, N.-B., Li, B.-A., Mar. 2019. Extracting nuclear symmetry ener-gies at high densities from observations of neutron stars and gravitationalwaves. European Physical Journal A 55 (3), 39.[80] Zimmerman, J., Carson, Z., Schumacher, K., Steiner, A. W., Yagi,K., Feb. 2020. Measuring Nuclear Matter Parameters with NICER andLIGO/Virgo. preprint (arXiv:2002.03210).[81] Shlomo, S., Kolomietz, V. M., Col`o, G., Oct. 2006. Deducing the nuclear-matter incompressibility coefficient from data on isoscalar compressionmodes. EPJA, 30 (1), 23–30. 7382] Arrington, J., Sick, I., Sep. 2015. Evaluation of the Proton Charge Ra-dius from Electron-Proton Scattering. Journal of Physical and ChemicalReference Data 44 (3), 031204.[83] Pohl, R., et al., Jul. 2010. The size of the proton. Natur, 466 (7303),213–216.[84] Potekhin, A. Y., Fantina, A. F., Chamel, N., Pearson, J. M., Goriely,S., Dec. 2013. Analytical representations of unified equations of state forneutron-star matter. A&A, 560, A48.[85] Shen, H., Toki, H., Oyamatsu, K., Sumiyoshi, K., Jul. 1998. Relativis-tic equation of state of nuclear matter for supernova and neutron star.NuPhA, 637 (3), 435–450.[86] Glendenning, N. K., Moszkowski, S. A., Oct. 1991. Reconciliation ofneutron-star masses and binding of the Lambda in hypernuclei. PhRvL,67, 2414–1417.[87] Li, A., Liu, T., Jul. 2013. Revisiting the hot matter in the center of gamma-ray bursts and supernovae. A&A, 555, A129.[88] Li, A., Liu, T., Gubler, P., Xu, R.-X., Mar. 2015. Revisiting the boilingof primordial quark nuggets at nonzero chemical potential. APh, 62, 115–121.[89] Li, A., Zhou, X. R., Burgio, G. F., Schulze, H. J., Feb. 2010. Protoneutronstars in the Brueckner-Hartree-Fock approach and finite-temperature kaoncondensation. PhRvC, 81 (2), 025806.[90] Zhu, Z.-Y., Li, A., Hu, J.-N., Sagawa, H., Oct. 2016. ∆ (1232) effectsin density-dependent relativistic Hartree-Fock theory and neutron stars.PhRvC, 94 (4), 045803.[91] Stone, J., Guichon, P. A. M., Matevosyan, H. H., Thomas, A. W., Aug.2007. Cold uniform matter and neutron stars in the quark meson-couplingmodel. NuPhA, 792 (3-4), 341–369.7492] Lattimer, J. M., Prakash, M., Apr. 2007. Neutron star observations: Prog-nosis for equation of state constraints. PhR, 442 (1-6), 109–165.[93] Lattimer, J. M., Nov. 2012. The Nuclear Equation of State and NeutronStar Masses. ARNPS, 62 (1), 485–515.[94] Bauswein, A., Just, O., Janka, H.-T., Stergioulas, N., Dec. 2017. Neutron-star Radius Constraints from GW170817 and Future Detections. ApJ,850 (2), L34.[95] Annala, E., Gorda, T., Kurkela, A., Vuorinen, A., Apr. 2018.Gravitational-Wave Constraints on the Neutron-Star-Matter Equation ofState. PhRvL, 120 (17), 172703.[96] Lim, Y., Holt, J. W., Aug. 2018. Neutron Star Tidal Deformabilities Con-strained by Nuclear Theory and Experiment. PhRvL, 121 (6), 062701.[97] Most, E. R., Weih, L. R., Rezzolla, L., Schaffner-Bielich, J., Jun. 2018.New Constraints on Radii and Tidal Deformabilities of Neutron Stars fromGW170817. PhRvL, 120 (26), 261103.[98] Nandi, R., Char, P., Apr. 2018. Hybrid Stars in the Light of GW170817.ApJ, 857 (1), 12.[99] Zhao, T., Lattimer, J. M., Sep. 2018. Tidal deformabilities and neutronstar mergers. PhRvD, 98 (6), 063020.[100] Zhou, E.-P., Zhou, X., Li, A., Apr. 2018. Constraints on interquark in-teraction parameters with GW170817 in a binary strange star scenario.PhRvD, 97 (8), 083015.[101] Han, S., Steiner, A. W., Apr. 2019. Tidal deformability with sharp phasetransitions in binary neutron stars. PhRvD, 99 (8), 083014. 2[102] Raithel, C. A., ¨Ozel, F., Psaltis, D., Apr. 2018. Tidal Deformability fromGW170817 as a Direct Probe of the Neutron Star Radius. ApJ, 857 (2),L23. 75103] Tsang, C. Y., Tsang, M. B., Danielewicz, P., Fattoyev, F. J., Lynch,W. G., Sep. 2019. Insights on Skyrme parameters from GW170817. PhLB,796, 1–5.[104] Wang, Q., Shi, C., Yan, Y., Zong, H.-S., Dec. 2019. Crossover hadron-quark transition with a modified interpolation method and constraintsfrom tidal deformability of GW170817. preprint (arXiv:1912.02312).[105] Zhou, Y., Chen, L.-W., Nov. 2019. Ruling Out the Supersoft High-densitySymmetry Energy from the Discovery of PSR J0740+6620 with Mass2 . +0 . − . M (cid:12) . ApJ, 886 (1), 52.[106] Essick, R., Landry, P., Holz, D. E., Mar. 2020. Nonparametric inferenceof neutron star composition, equation of state, and maximum mass withGW170817. PhRvD, 101 (6), 063007.[107] Ferreira, M., Fortin, M., Malik, T., Agrawal, B. K., Providˆencia, C., Feb.2020. Empirical constraints on the high-density equation of state frommultimessenger observables. PhRvD, 101 (4), 043021.[108] G¨uven, H., Bozkurt, K., Khan, E., Margueron, J., Jan. 2020. Multi-messenger and multi-physics bayesian inference for GW170817 binary neu-tron star merger. preprint (arXiv:2001.10259).[109] Louren¸co, O., Bhuyan, M., Lenzi, C. H., Dutra, M., Gonzalez-Boquera,C., Centelles, M., Vi˜nas, X., Apr. 2020. GW170817 constraints analyzedwith Gogny forces and momentum-dependent interactions. PhLB, 803,135306.[110] Louren¸co, O., Dutra, M., Lenzi, C. H., Biswal, S. K., Bhuyan, M.,Menezes, D. P., Feb. 2020. Consistent Skyrme parametrizations con-strained by GW170817. EPJA, 56 (2), 32.[111] Traversi, S., Char, P., Pagliara, G., Feb. 2020. Bayesian Inference of DenseMatter Equation of State within Relativistic Mean Field Models usingAstrophysical Measurements. preprint (arXiv:2002.08951).76112] Li, A., Zuo, W., Mi, A.-J., G, B., Jul. 2007. Hyperon hyperon interactionin relativistic mean field model. ChPhy, 16 (7), 1934–1940.[113] Burgio, G. F., Schulze, H. J., Li, A., Feb. 2011. Hyperon stars at finitetemperature in the Brueckner theory. PhRvC, 83 (2), 025804.[114] Zuo, W., Li, A., Li, Z. H., Lombardo, U., Nov. 2004. Nuclear three-bodyforce effect on a kaon condensate in neutron star matter. PhRvC, 70 (5),055802.[115] Li, A., Burgio, G. F., Lombardo, U., Zuo, W., Nov. 2006. Microscopicthree-body forces and kaon condensation in cold neutrino-trapped matter.PhRvC, 74 (5), 055801.[116] Li, A., Zuo, W., Schulze J., H., Lombardo, U., Dec. 2008. Hot NuclearMatter Equation of State and Finite Temperature Kaon Condensation.ChPhL, 25 (12), 4233–4236.[117] Li, A., Peng, G.-X., Lombardo, U., Mar. 2009. Deconfinement phase tran-sition in neutron star matter. ChPhC, 33 (S1), 61–63.[118] Li, A., Zuo, W., Peng, G. X., Mar. 2015. Massive hybrid stars with afirst-order phase transition. PhRvC, 91 (3), 035803.[119] Chamel, N., Haensel, P., Dec. 2008. Physics of Neutron Star Crusts. LRR,11 (1), 10.[120] Piekarewicz, J., Fattoyev, F. J., Horowitz, C. J., Jul. 2014. Pulsar glitches:The crust may be enough. PhRvC, 90 (1), 015803.[121] Li, A., Jul. 2015. Glitch Crisis or Not: a Microscopic Study. ChinesePhysics Letters 32 (7), 079701.[122] Li, A., Dong, J. M., Wang, J. B., Xu, R. X., Mar. 2016. Structures ofthe Vela Pulsar and the Glitch Crisis from the Brueckner Theory. ApJS,223 (1), 16. 77123] Fantina, A. F., Chamel, N., Pearson, J. M., Goriely, S., Nov. 2013. Neu-tron star properties with unified equations of state of dense matter. A&A,559, A128.[124] Antic, S., Stone, J. R., Thomas, A. W., Dec. 2019. Neutron starsfrom crust to core within the Quark-meson coupling model. preprint(arXiv:1912.03815)[125] Fortin, M., Providˆencia, C., Raduta, A. R., Gulminelli, F., Zdunik, J. L.,Haensel, P., Bejger, M., Sep. 2016. Neutron star radii and crusts: Uncer-tainties and unified equations of state. PhRvC, 94 (3), 035804.[126] Baym, G., Pethick, C., Sutherland, P., Dec. 1971. The Ground State ofMatter at High Densities: Equation of State and Stellar Models. ApJ,170, 299.[127] Negele, J. W., Vautherin, D., Jun. 1973. Neutron star matter at sub-nuclear densities. NuPhA, 207 (2), 298–320.[128] Stergioulas, N., Jun. 2003. Rotating Stars in Relativity. LRR, 6 (1), 3.[129] Fortin, M., Taranto, G., Burgio, G. F., Haensel, P., Schulze, H. J., Zdunik,J. L., Apr. 2018. Thermal states of neutron stars with a consistent modelof interior. MNRAS, 475 (4), 5010–5022.[130] Haensel, P., Zdunik, J. L., Mar. 2008. Models of crustal heating in accret-ing neutron stars. A&A, 480 (2), 459–464.[131] Tolman, R. C., Feb. 1939. Static Solutions of Einstein’s Field Equationsfor Spheres of Fluid. PhRv, 55 (4), 364–373.[132] Oppenheimer, J. R., Volkoff, G. M., Feb. 1939. On Massive Neutron Cores.PhRv, 55 (4), 374–381.[133] Hornick, N., Tolos, L., Zacchi, A., Christian, J.-E., Schaffner-Bielich, J.,Dec. 2018. Relativistic parameterizations of neutron matter and implica-tions for neutron stars. PhRvC, 98 (6), 065804.78134] Capano, C. D., Tews, I., Brown, S. M., Margalit, B., De, S., Kumar, S.,Brown, D. A., Krishnan, B., Reddy, S., Mar. 2020. Stringent constraintson neutron-star radii from multimessenger observations and nuclear the-ory. NatAs, 1–8.[135] Postnikov, S., Prakash, M., Lattimer, J. M., Jul. 2010. Tidal Love numbersof neutron and self-bound quark stars. PhRvD, 82 (2), 024016.[136] Read, J. S., et al., Aug. 2013. Matter effects on binary neutron star wave-forms. PhRvD, 88 (4), 044042.[137] Komatsu, H., Eriguchi, Y., Hachisu, I., Mar. 1989. Rapidly rotating gen-eral relativistic stars. I - Numerical method and its application to uni-formly rotating polytropes. MNRAS, 237, 355–379.[138] Cook, G. B., Shapiro, S. L., Teukolsky, S. A., Feb. 1994. Rapidly RotatingPolytropes in General Relativity. ApJ, 422, 227.[139] Stergioulas, N., Friedman, J. L., May 1995. Comparing Models of RapidlyRotating Relativistic Stars Constructed by Two Numerical Methods. ApJ,444, 306.[140] Haensel, P., Zdunik, J. L., Bejger, M., Lattimer, J. M., Aug. 2009. Ke-plerian frequency of uniformly rotating neutron stars and strange stars.A&A, 502 (2), 605–610.[141] Wei, J. B., Chen, H., Burgio, G. F., Schulze, H. J., Aug. 2017. Rotat-ing hybrid stars with the Dyson-Schwinger quark model. PhRvD, 96 (4),043008.[142] Andersson, N., Kokkotas, K. D., Jan. 2001. The R-Mode Instability inRotating Neutron Stars. IJMPD, , 10 (4), 381–441.[143] Breu, C., Rezzolla, L., Jun. 2016. Maximum mass, moment of inertia andcompactness of relativistic stars. MNRAS, 459 (1), 646–656.79144] Urbanec, M., Miller, J. C., Stuchl´ık, Z., Aug. 2013. Quadrupole momentsof rotating neutron stars and strange stars. MNRAS, 433 (3), 1903–1909.[145] Yagi, K., Yunes, N., Jul. 2013. I-Love-Q: Unexpected Universal Relationsfor Neutron Stars and Quark Stars. Sci, 341 (6144), 365–368.[146] Yagi, K., Yunes, N., Jul. 2013. I-Love-Q relations in neutron stars andtheir applications to astrophysics, gravitational waves, and fundamentalphysics. PhRvD, 88 (2), 023009.[147] Haskell, B., Ciolfi, R., Pannarale, F., Rezzolla, L., Feb. 2014. On theuniversality of I-Love-Q relations in magnetized neutron stars. MNRAS,438 (1), L71–L75.[148] Pappas, G., Apostolatos, T. A., Mar. 2014. Effectively Universal Behav-ior of Rotating Neutron Stars in General Relativity Makes Them EvenSimpler than Their Newtonian Counterparts. PhRvL, 112 (12), 121101.[149] Zhou, E., Tsokaros, A., Ury¯u, K., Xu, R., Shibata, M., Aug. 2019. Differ-entially rotating strange star in general relativity. PhRvD, 100 (4), 043015.[150] Zhou, E., Tsokaros, A., Rezzolla, L., Xu, R., Ury¯u, K., Jan. 2018. Uni-formly rotating, axisymmetric, and triaxial quark stars in general relativ-ity. PhRvD, 97 (2), 023013.[151] Khaustov, P. et al., May 2000. Evidence of Ξ hypernuclear production inthe C(K − ,K + ) Be reaction. PhRvC, 61 (5), 054603.[152] Yakovlev, D. G., Kaminker, A. D., Gnedin, O. Y., Haensel, P., Nov. 2001.Neutrino emission from neutron stars. PhR, 354 (1-2), 1–155.[153] Rijken, T. A., Schulze, H. J., Feb. 2016. Hyperon-hyperon interactionswith the Nijmegen ESC08 model. EPJA, 52, 21.[154] Bombaci, I., Jan. 2017. The Hyperon Puzzle in Neutron Stars. hspp.conf,101002. 80155] Yamamoto, Y., Furumoto, T., Yasutake, N., Rijken, T. A., Oct. 2014. Hy-peron mixing and universal many-body repulsion in neutron stars. PhRvC,90 (4), 045805.[156] Haidenbauer, J., Meißner, U. G., Kaiser, N., Weise, W., Jun. 2017.Lambda-nuclear interactions and hyperon puzzle in neutron stars. EPJA,53 (6), 121.[157] Logoteta, D., Vida˜na, I., Bombaci, I., Nov. 2019. Impact of chiral hyper-onic three-body forces on neutron stars. EPJA, 55 (11), 207.[158] Peng, G. X., Li, A., Lombardo, U., Jun. 2008. Deconfinement phase tran-sition in hybrid neutron stars from the Brueckner theory with three-bodyforces and a quark model with chiral mass scaling. PhRvC, 77 (6), 065807.[159] Li, A., Xu, R.-X., Lu, J.-F., Mar. 2010. Strange stars with different quarkmass scalings. MNRAS, 402 (4), 2715–2719.[160] Li, A., Peng, G.-X., Lu, J.-F., Apr. 2011. Strange star candidates revisedwithin a quark model with chiral mass scaling. RAA, 11 (4), 482–490.[161] Haensel, P., Zdunik, J. L., Schaefer, R., May 1986. Strange quark stars.A&A, 160 (1), 121–128.[162] Wiktorowicz, G., Drago, A., Pagliara, G., Popov, S. B., Sep. 2017. StrangeQuark Stars in Binaries: Formation Rates, Mergers, and Explosive Phe-nomena. ApJ, 846 (2), 163.[163] Drago, A., Pagliara, G., Jan. 2018. Merger of Two Neutron Stars: Pre-dictions from the Two-families Scenario. ApJ, 852 (2), L32.[164] Alford, M., Han, S., Schwenzer, K., Nov. 2019. Signatures for quark matterfrom multi-messenger observations. JPhG, 46 (11), 114001.[165] Aloy, M. A., Ib´a˜nez, J. M., Sanchis-Gual, N., Obergaulinger, M., Font,J. A., Serna, S., Marquina, A., Apr. 2019. Neutron star collapse and81ravitational waves with a non-convex equation of state. MNRAS, 484 (4),4980–5008.[166] Bauswein, A., Bastian, N.-U. F., Blaschke, D. B., Chatziioannou, K.,Clark, J. A., Fischer, T., Oertel, M., Feb. 2019. Identifying a First-OrderPhase Transition in Neutron-Star Mergers through Gravitational Waves.PhRvL, 122 (6), 061102.[167] Most, E. R., Papenfort, L. J., Dexheimer, V., Hanauske, M., Schramm, S.,St¨ocker, H., Rezzolla, L., Feb. 2019. Signatures of Quark-Hadron PhaseTransitions in General-Relativistic Neutron-Star Mergers. PhRvL, 122 (6),061101.[168] Weih, L. R., Most, E. R., Rezzolla, L., Aug. 2019. Optimal Neutron-starMass Ranges to Constrain the Equation of State of Nuclear Matter withElectromagnetic and Gravitational-wave Observations. ApJ, 881 (1), 73.[169] Xia, C., Zhu, Z., Zhou, X., Li, A., Jun. 2019. Sound velocity indense stellar matter with strangeness and compact stars. preprint(arXiv:1906.00826).[170] Gomes, R. O., Char, P., Schramm, S., Jun. 2019. Constraining Strangenessin Dense Matter with GW170817. ApJ, 877 (2), 139.[171] Chatziioannou, K., Han, S., Feb. 2020. Studying strong phase transitionsin neutron stars with gravitational waves. PhRvD, 101 (4), 044019.[172] Fischer, T., Wu, M.-R., Wehmeyer, B., Bastian, N.-U. F., Mart´ınez-Pinedo, G., Thielemann, F.-K., Mar. 2020. Core-collapse supernova ex-plosions driven by the hadron-quark phase transition as rare r processsite. preprint (arXiv:2003.00972).[173] Tonetto, L., Lugones, G., Mar. 2020. Discontinuity gravity modes in hy-brid stars: assessing the role of rapid and slow phase conversions. preprint(arXiv:2003.01259). 82174] Nunna, K. P., Banik, S., Chatterjee, D., Feb. 2020. Signatures ofstrangeness in neutron star merger remnants. preprint (arXiv:2002.07538).[175] Alford, M., Han, S., Prakash, M., Oct. 2013. Generic conditions for stablehybrid stars. PhRvD, 88 (8), 083013.[176] Zdunik, J. L., Haensel, P., Mar. 2013. Maximum mass of neutron starsand strange neutron-star cores. A&A, 551, A61.[177] Alford, M., Burgio, G. F., Han, S., Taranto, G., Zappal`a, D., Oct.2015. Constraining and applying a generic high-density equation of state.PhRvD, 92 (8), 083002.[178] Ranea-Sandoval, I. F., Han, S., Orsaria, M. G., Contrera, G. A., We-ber, F., Alford, M. G., Apr. 2016. Constant-sound-speed parametrizationfor Nambu-Jona-Lasinio models of quark matter in hybrid stars. PhRvC,93 (4), 045812.[179] Kurkela, A., Romatschke, P., Vuorinen, A., May2010. Cold quark matter.PhRvD, 81 (10), 105021.[180] Abbott, B. P. et al., Oct. 2018. GW170817: Measurements of NeutronStar Radii and Equation of State. PhRvL, 121 (16), 161101.[181] Monta˜na, G., Tol´os, L., Hanauske, M., Rezzolla, L., May 2019. Constrain-ing twin stars with GW170817. PhRvD, 99 (10), 103009.[182] Burgio, G. F., Drago, A., Pagliara, G., Schulze, H. J., Wei,J. B., Jun. 2018. Are Small Radii of Compact Stars Ruled out byGW170817/AT2017gfo? ApJ, 860 (2), 139.[183] De, S., Finstad, D., Lattimer, J. M., Brown, D. A., Berger, E., Biwer,C. M., Aug. 2018. Tidal Deformabilities and Radii of Neutron Stars fromthe Observation of GW170817. PhRvL, 121 (9), 091102.[184] Miao, Z. Q., Li, A., Zhu Z. Y. and Han, S., 2020. Constraining hadron-quark phase transition parameters within the quark-mean-field model83sing multi-messenger observations of neutron stars. preprint (arXiv:2006.00839)[185] Kalogera, V., Baym, G., Oct. 1996. The Maximum Mass of a NeutronStar. ApJ, 470, L61.[186] Rhoades, C. E., Ruffini, R., Feb. 1974. Maximum Mass of a Neutron Star.PhRvL, 32 (6), 324–327.[187] Brecher, K., Caporaso, G., Feb. 1976. Obese ‘neutron’ stars. Natur,259 (5542), 377–378.[188] Glendenning, N. K., Jan. 1990. Fast Pulsars, Strange Stars:. AN Oppor-tunity in Radio Astronomy. MPLA, 5 (27), 2197–2207.[189] Burgio, G. F., Schulze, H. J., Weber, F., Sep. 2003. On the maximumrotational frequency of neutron and hybrid stars. A&A, 408, 675–680.[190] Glendenning, N. K., Aug. 1992. First-order phase transitions with morethan one conserved charge: Consequences for neutron stars. PhRvD,46 (4), 1274–1287.[191] Maruyama, T., Chiba, S., Schulze, H.-J., Tatsumi, T., Dec. 2007. Hadron-quark mixed phase in hyperon stars. PhRvD, 76 (12), 123015.[192] Xia, C.-J., Maruyama, T., Yasutake, N., Tatsumi, T., May 2019.Constraining quark-hadron interface tension in the multimessenger era.PhRvD, 99 (10), 103017.[193] Kurkela, A., Fraga, E. S., Schaffner-Bielich, J., Vuorinen, A., Jul. 2014.Constraining Neutron Star Matter with Quantum Chromodynamics. ApJ,789 (2), 127.[194] Fraga, E. S., Romatschke, P., May 2005. Role of quark mass in cold anddense perturbative QCD. PhRvD, 71 (10), 105014.84195] Olive, K. A., Particle Data Group, Aug. 2014. Review of Particle Physics.ChPhC, 38 (9), 090001.[196] Vermaseren, J. A. M., Larin, S. A., van Ritbergen, T., Feb. 1997. The4-loop quark mass anomalous dimension and the invariant quark mass.PhLB, 405, 327–333.[197] Burgio, G. F., Baldo, M., Sahu, P. K., Santra, A. B., Schulze, H. J., Jan.2002. Maximum mass of neutron stars with a quark core. PhLB, 526 (1-2),19–26.[198] Maieron, C., Baldo, M., Burgio, G. F., Schulze, H. J., Aug. 2004. Hybridstars with the color dielectric and the MIT bag models. PhRvD, 70 (4),043010.[199] Buballa, M., Feb. 2005. NJL-model analysis of dense quark matter [reviewarticle]. PhR, 407 (4-6), 205–376.[200] Bedaque, P., Steiner, A. W., Jan. 2015. Sound Velocity Bound and Neu-tron Stars. PhRvL, 114 (3), 031103.[201] Baym, G., Furusawa, S., Hatsuda, T., Kojo, T., Togashi, H., Nov. 2019.New Neutron Star Equation of State with Quark-Hadron Crossover. ApJ,885 (1), 42.[202] Alsing, J., Silva, H. O., Berti, E., Jul. 2018. Evidence for a maximummass cut-off in the neutron star mass distribution and constraints on theequation of state. MNRAS, 478 (1), 1377–1391.[203] Tews, I., Carlson, J., Gandolfi, S., Reddy, S., Jun. 2018. Constraining theSpeed of Sound inside Neutron Stars with Chiral Effective Field TheoryInteractions and Observations. ApJ, 860 (2), 149.[204] McLerran, L., Reddy, S., Mar. 2019. Quarkyonic Matter and NeutronStars. PhRvL, 122 (12), 122701. 85205] Li, A., Zhu, Z.-Y., Zhou, X., Jul. 2017. New Equations of State for Post-merger Supramassive Quark Stars. ApJ, 844 (1), 41.[206] Li, A., Zhang, B., Zhang, N.-B., Gao, H., Qi, B., Liu, T., Oct. 2016.Internal x-ray plateau in short GRBs: Signature of supramassive fast-rotating quark stars? PhRvD, 94 (8), 083010.[207] Bodmer, A. R., Sep. 1971. Collapsed Nuclei. PhRvD, 4 (6), 1601–1606.[208] Witten, E., Jul. 1984. Cosmic separation of phases. PhRvD, 30 (2), 272–285.[209] Alford, M., Braby, M., Paris, M., Reddy, S., Aug. 2005. Hybrid Stars thatMasquerade as Neutron Stars. ApJ, 629 (2), 969–978.[210] Wei, W., Irving, B., Salinas, M., Kl¨ahn, T., Jaikumar, P., Dec. 2019.Camouflage of the Phase Transition to Quark Matter in Neutron Stars.ApJ, 887 (2), 151.[211] Rezzolla, L., Most, E. R., Weih, L. R., Jan. 2018. Using Gravitational-wave Observations and Quasi-universal Relations to Constrain the Maxi-mum Mass of Neutron Stars. ApJ, 852 (2), L25.[212] Ruiz, M., Shapiro, S. L., Tsokaros, A., Jan. 2018. GW170817, general rel-ativistic magnetohydrodynamic simulations, and the neutron star maxi-mum mass. PhRvD, 97 (2), 021501.[213] Shibata, M., Zhou, E., Kiuchi, K., Fujibayashi, S., Jul. 2019. Constrainton the maximum mass of neutron stars using GW170817 event. PhRvD,100 (2), 023015.[214] Lai, X.-Y., Yu, Y.-W., Zhou, E.-P., Li, Y.-Y., Xu, R.-X., Feb. 2018. Merg-ing strangeon stars. Research in Astronomy and Astrophysics 18 (2), 024.[215] Bauswein, A., Janka, H. T., Oechslin, R., Pagliara, G., Sagert, I.,Schaffner-Bielich, J., Hohle, M. M., Neuh¨auser, R., Jul. 2009. Mass Ejec-86ion by Strange Star Mergers and Observational Implications. PhRvL,103 (1), 011101.[216] Alford, M., Rajagopal, K., Wilczek, F., Jan. 1999. Color-flavor locking andchiral symmetry breaking in high density QCD. NuPhB, 537 (1), 443–458.[217] Bailyn, C. D., Jain, R. K., Coppi, P., Orosz, J. A., May 1998. The MassDistribution of Stellar Black Holes. ApJ, 499 (1), 367–374.[218] ¨Ozel, F., Psaltis, D., Narayan, R., McClintock, J. E., Dec. 2010. The BlackHole Mass Distribution in the Galaxy. ApJ, 725 (2), 1918–1927.[219] Farr, W. M., Sravan, N., Cantrell, A., Kreidberg, L., Bailyn, C. D., Man-del, I., Kalogera, V., Nov. 2011. The Mass Distribution of Stellar-massBlack Holes. ApJ, 741 (2), 103.[220] Kreidberg, L., Bailyn, C. D., Farr, W. M., Kalogera, V., Sep. 2012. MassMeasurements of Black Holes in X-Ray Transients: Is There a Mass Gap?ApJ, 757 (1), 36.[221] Wyrzykowski, (cid:32)L., Mandel, I., Apr. 2019. Constraining the masses ofmicrolensing black holes and the mass gap with Gaia DR2. preprint(arXiv:1904.07789).[222] Tsokaros, A., Ruiz, M., Shapiro, S. L., Sun, L., Ury¯u, K., Feb. 2020.Great Impostors: Extremely Compact, Merging Binary Neutron Stars inthe Mass Gap Posing as Binary Black Holes. PhRvL, 124 (7), 071101.[223] Abbott, B. P. et al., Jan. 2020. GW190425: Observation of a Compact Bi-nary Coalescence with Total Mass ∼ . M (cid:12) . preprint (arXiv:2001.01761).[224] Abbott, B. P. et al., Dec. 2017. Estimating the Contribution of DynamicalEjecta in the Kilonova Associated with GW170817. ApJ, 850 (2), L39.[225] Abbott, B. P. et al., Dec. 2017. Search for Post-merger GravitationalWaves from the Remnant of the Binary Neutron Star Merger GW170817.ApJ, 851 (1), L16. 87226] Narayan, R., Paczynski, B., Piran, T., Aug. 1992. Gamma-Ray Bursts asthe Death Throes of Massive Binary Stars. ApJ, 395, L83.[227] Eichler, D., Livio, M., Piran, T., Schramm, D. N., Jul. 1989. Nucleosyn-thesis, neutrino bursts and γ -rays from coalescing neutron stars. Natur,340 (6229), 126–128.[228] Fattoyev, F. J., Piekarewicz, J., Horowitz, C. J., Apr. 2018. Neutron Skinsand Neutron Stars in the Multimessenger Era. PhRvL, 120 (17), 172702.[229] Flanagan, ´E. ´E., Hinderer, T., Jan. 2008. Constraining neutron-star tidalLove numbers with gravitational-wave detectors. PhRvD, 77 (2), 021502.[230] Lawrence, S., Tervala, J. G., Bedaque, P. F., Miller, M. C., Aug. 2015. AnUpper Bound on Neutron Star Masses from Models of Short Gamma-RayBursts. ApJ, 808 (2), 186.[231] Murguia-Berthier, A., Montes, G., Ramirez-Ruiz, E., De Colle, F., Lee,W. H., Jun. 2014. Necessary Conditions for Short Gamma-Ray Burst Pro-duction in Binary Neutron Star Mergers. ApJ, 788 (1), L8.[232] Ai, S., Gao, H., Dai, Z.-G., Wu, X.-F., Li, A., Zhang, B., Li, M.-Z., Jun.2018. The Allowed Parameter Space of a Long-lived Neutron Star as theMerger Remnant of GW170817. ApJ, 860 (1), 57.[233] Dai, Z. G., Wang, X. Y., Wu, X. F., Zhang, B., Feb. 2006. X-ray Flaresfrom Postmerger Millisecond Pulsars. Sci, 311 (5764), 1127–1129.[234] Metzger, B. D., Quataert, E., Thompson, T. A., Apr. 2008. Short-durationgamma-ray bursts with extended emission from protomagnetar spin-down.MNRAS, 385 (3), 1455–1460.[235] Hinderer, T., Apr. 2008. Tidal Love Numbers of Neutron Stars. ApJ,677 (2), 1216–1220.[236] Douchin, F., Haensel, P., Dec. 2001. A unified equation of state of densematter and neutron star structure. A&A, 380, 151–167.88237] Fraga, E. S., Pisarski, R. D., Schaffner-Bielich, J., Jun. 2001. Small, densequark stars from perturbative QCD. PhRvD, 63 (12), 121702.[238] Damour, T., Nagar, A., Oct. 2009. Relativistic tidal properties of neutronstars. PhRvD, 80 (8), 084035.[239] Lai, X., Zhou, E., Xu, R., Apr. 2019. Strangeons constitute bulk strongmatter: Test using GW 170817. EPJA, 55 (4), 60.[240] Gao, H., Zhang, B., L¨u, H.-J., Feb. 2016. Constraints on binary neutronstar merger product from short GRB observations. PhRvD, 93 (4), 044065.[241] Ma, P.-X., Jiang, J.-L., Wang, H., Jin, Z.-P., Fan, Y.-Z., Wei, D.-M., May2018. GW170817 and the Prospect of Forming Supramassive Remnantsin Neutron Star Mergers. ApJ, 858 (2), 74.[242] Gao, H., Ai, S.-K., Cao, Z.-J., Zhang, B., Zhu, Z.-Y., Li, A., Zhang,N.-B., Bauswein, A., Jan. 2020. Relation between gravitational mass andbaryonic mass for non-rotating and rapidly rotating neutron stars. FrPhy,15 (2), 24603.[243] Abbott, B. P. et al., Oct. 2017. Gravitational Waves and Gamma-Raysfrom a Binary Neutron Star Merger: GW170817 and GRB 170817A. ApJ,848 (2), L13.[244] Margalit, B., Metzger, B. D., Dec. 2017. Constraining the Maximum Massof Neutron Stars from Multi-messenger Observations of GW170817. ApJ,850 (2), L19.[245] Shibata, M., Fujibayashi, S., Hotokezaka, K., Kiuchi, K., Kyutoku, K.,Sekiguchi, Y., Tanaka, M., Dec. 2017. Modeling GW170817 based on nu-merical relativity and its implications. PhRvD, 96 (12), 123012.[246] Li, S.-Z., Liu, L.-D., Yu, Y.-W., Zhang, B., Jul. 2018. What Powered theOptical Transient AT2017gfo Associated with GW170817? ApJ, 861 (2),L12. 89247] Piro, L., et al., Feb. 2019. A long-lived neutron star merger remnantin GW170817: constraints and clues from X-ray observations. MNRAS,483 (2), 1912–1921.[248] Li, L.-X., Paczy´nski, B., Nov. 1998. Transient Events from Neutron StarMergers. ApJ, 507 (1), L59–L62.[249] Metzger, B. D., et al., Aug. 2010. Electromagnetic counterparts of com-pact object mergers powered by the radioactive decay of r-process nuclei.MNRAS, 406 (4), 2650–2662.[250] Kasen, D., Metzger, B., Barnes, J., Quataert, E., Ramirez-Ruiz, E., Nov.2017. Origin of the heavy elements in binary neutron-star mergers from agravitational-wave event. Natur, 551 (7678), 80–84.[251] Ruffert, M., Janka, H. T., Takahashi, K., Schaefer, G., Mar. 1997. Coalesc-ing neutron stars - a step towards physical models. II. Neutrino emission,neutron tori, and gamma-ray bursts. A&A, 319, 122–153.[252] Rosswog, S., Liebend¨orfer, M., Thielemann, F. K., Davies, M. B., Benz,W., Piran, T., Jan. 1999. Mass ejection in neutron star mergers. A&A,341, 499–526.[253] Oechslin, R., Janka, H. T., Marek, A., May 2007. Relativistic neutronstar merger simulations with non-zero temperature equations of state. I.Variation of binary parameters and equation of state. A&A, 467 (2), 395–409.[254] Hotokezaka, K., Kiuchi, K., Kyutoku, K., Okawa, H., Sekiguchi, Y.-i.,Shibata, M., Taniguchi, K., Jan. 2013. Mass ejection from the merger ofbinary neutron stars. PhRvD, 87 (2), 024001.[255] Wanajo, S., Sekiguchi, Y., Nishimura, N., Kiuchi, K., Kyutoku, K., Shi-bata, M., Jul. 2014. Production of All the r-process Nuclides in the Dy-namical Ejecta of Neutron Star Mergers. ApJ, 789 (2), L39.90256] Sekiguchi, Y., Kiuchi, K., Kyutoku, K., Shibata, M., Taniguchi, K., Jun.2016. Dynamical mass ejection from the merger of asymmetric binary neu-tron stars: Radiation-hydrodynamics study in general relativity. PhRvD,93 (12), 124046.[257] Metzger, B. D., Piro, A. L., Quataert, E., Oct. 2008. Time-dependentmodels of accretion discs formed from compact object mergers. MNRAS,390 (2), 781–797.[258] Dessart, L., Ott, C. D., Burrows, A., Rosswog, S., Livne, E., Jan. 2009.Neutrino Signatures and the Neutrino-Driven Wind in Binary NeutronStar Mergers. ApJ, 690 (2), 1681–1705.[259] Fern´andez, R., Metzger, B. D., Oct. 2013. Delayed outflows from blackhole accretion tori following neutron star binary coalescence. MNRAS,435 (1), 502–517.[260] Metzger, B. D., Fern´andez, R., Jul. 2014. Red or blue? A potential kilo-nova imprint of the delay until black hole formation following a neutronstar merger. MNRAS, 441 (4), 3444–3453.[261] Perego, A., Rosswog, S., Cabez´on, R. M., Korobkin, O., K¨appeli, R.,Arcones, A., Liebend¨orfer, M., Oct. 2014. Neutrino-driven winds fromneutron star merger remnants. MNRAS, 443 (4), 3134–3156.[262] Martin, D., Perego, A., Arcones, A., Thielemann, F. K., Korobkin, O.,Rosswog, S., Nov. 2015. Neutrino-driven Winds in the Aftermath of aNeutron Star Merger: Nucleosynthesis and Electromagnetic Transients.ApJ, 813 (1), 2.[263] Cowperthwaite, P. S. et al., Oct. 2017. The Electromagnetic Counterpartof the Binary Neutron Star Merger LIGO/Virgo GW170817. II. UV, Opti-cal, and Near-infrared Light Curves and a comparison to Kilonova Models.ApJ, 848 (2), L17. 91264] Tanaka, M. et al., Dec. 2017. Kilonova from post-merger ejecta as anoptical and near-Infrared counterpart of GW170817. PASJ, 69 (6), 102.[265] Radice, D., Perego, A., Zappa, F., Bernuzzi, S., Jan. 2018. GW170817:Joint Constraint on the Neutron Star Equation of State from Multimes-senger Observations. ApJ, 852 (2), L29.[266] Radice, D., Perego, A., Hotokezaka, K., Fromm, S. A., Bernuzzi, S.,Roberts, L. F., Dec. 2018. Binary Neutron Star Mergers: Mass Ejection,Electromagnetic Counterparts, and Nucleosynthesis. ApJ, 869 (2), 130.[267] Kiuchi, K., Kyutoku, K., Shibata, M., Taniguchi, K., May 2019. Revisitingthe Lower Bound on Tidal Deformability Derived by AT 2017gfo. ApJ,876 (2), L31.[268] Ai, S., Gao, H., Zhang, B., Dec. 2019. What constraints on the neutronstar maximum mass can one pose from GW170817 observations? preprint(arXiv:1912.06369).[269] Drago, A., Lavagno, A., Pagliara, G., Feb. 2014. Can very compact andvery massive neutron stars both exist? PhRvD, 89 (4), 043014.[270] Alcock, C., Farhi, E., Olinto, A., Nov. 1986. Strange Stars. ApJ, 310, 261.[271] Gondek-Rosi´nska, D., Bulik, T., Zdunik, L., Gourgoulhon, E., Ray, S.,Dey, J., Dey, M., Nov. 2000. Rapidly rotating compact strange stars.A&A, 363, 1005–1012.[272] Haensel, P., Paczynski, B., Amsterdamski, P., Jul. 1991. Gamma-RayBursts from Colliding Strange Stars. ApJ, 375, 209.[273] Alcock, C., Farhi, E., Sep. 1985. Evaporation of strange matter in theearly Universe. PhRvD, 32 (6), 1273–1279.[274] Madsen, J., Heiselberg, H., Riisager, K., Nov. 1986. Does strange matterevaporate in the early Universe \ ? PhRvD, 34 (10), 2947–2955.92275] Paulucci, L., Horvath, J. E., Benvenuto, O., Jan. 2017. Nucleosynthesisin Strange Star Mergers. IJMPS, 45, 1760042.[276] De Pietri, R., Drago, A., Feo, A., Pagliara, G., Pasquali, M., Traversi, S.,Wiktorowicz, G., Aug. 2019. Merger of Compact Stars in the Two-familiesScenario. ApJ, 881 (2), 122.[277] Drago, A., Lavagno, A., Metzger, B. D., Pagliara, G., May 2016. Quark de-confinement and the duration of short gamma-ray bursts. PhRvD, 93 (10),103001.[278] Page, D., Geppert, U., Weber, F., Oct. 2006. The cooling of compactstars. NuPhA, 777, 497–530.[279] Page, D., Reddy, S., Nov. 2006. Dense Matter in Compact Stars: The-oretical Developments and Observational Constraints. ARNPS, 56 (1),327–374.[280] Potekhin, A. Y., Pons, J. A., Page, D., Oct. 2015. Neutron Stars—Coolingand Transport. SSRv, 191 (1-4), 239–291.[281] Li, A., Hu, J. N., Shang, X. L., Zuo, W., Jan. 2016. Nonrelativistic nucleoneffective masses in nuclear matter: Brueckner-Hartree-Fock model versusrelativistic Hartree-Fock model. PhRvC, 93 (1), 015803.[282] Shang, X. L., Li, A., Miao, Z. Q., Burgio, G. F., Schulze, H. J., Jun. 2020.Nucleon effective mass in hot dense matter. PhRvC, 101 (6), 065801.[283] Dong, J. M., Lombardo, U., Zhang, H. F., Zuo, W., Jan. 2016. Role ofNucleonic Fermi Surface Depletion in Neutron Star Cooling. ApJ, 817 (1),6.[284] Yakovlev, D. G., Levenfish, K. P., Haensel, P., Aug. 2003. Thermal stateof transiently accreting neutron stars. A&A, 407, 265–271.93285] Deibel, A., Cumming, A., Brown, E. F., Reddy, S., Apr. 2017. Late-timeCooling of Neutron Star Transients and the Physics of the Inner Crust.ApJ, 839 (2), 95.[286] Archibald, R. F., Kaspi, V. M., Ng, C. Y., Gourgouliatos, K. N., Tsang,D., Scholz, P., Beardmore, A. P., Gehrels, N., Kennea, J. A., May 2013.An anti-glitch in a magnetar. Natur, 497 (7451), 591–593.[287] Radhakrishnan, V., Manchester, R. N., Apr. 1969. Detection of a Changeof State in the Pulsar PSR 0833-45. Natur, 222 (5190), 228–229.[288] Espinoza, C. M., Lyne, A. G., Stappers, B. W., Kramer, M., Jun. 2011.A study of 315 glitches in the rotation of 102 pulsars. MNRAS, 414 (2),1679–1704.[289] Yu, M., Manchester, R. N., Hobbs, G., Johnston, S., Kaspi, V. M., Keith,M., Lyne, A. G., Qiao, G. J., Ravi, V., Sarkissian, J. M., Shannon, R.,Xu, R. X., Feb. 2013. Detection of 107 glitches in 36 southern pulsars.MNRAS, 429 (1), 688–724.[290] Manchester, R. N., Jan. 2018. Pulsar glitches and their impact on neutron-star astrophysics. preprint (arXiv:1801.04332).[291] Haskell, B., Melatos, A., Jan. 2015. Models of pulsar glitches. IJMPD,24 (3), 1530008.[292] Anderson, P. W., Itoh, N., Jul. 1975. Pulsar glitches and restlessness as ahard superfluidity phenomenon. Natur, 256 (5512), 25–27.[293] Link, B., Epstein, R. I., Lattimer, J. M., Oct. 1999. Pulsar Constraintson Neutron Star Structure and Equation of State. PhRvL, 83 (17), 3362–3365.[294] Andersson, N., Glampedakis, K., Ho, W. C. G., Espinoza, C. M., Dec.2012. Pulsar Glitches: The Crust is not Enough. PhRvL, 109 (24), 241103.94295] Chamel, N., Jan. 2013. Crustal Entrainment and Pulsar Glitches. PhRvL,110 (1), 011101.[296] Watanabe, G., Pethick, C. J., Aug. 2017. Superfluid Density of Neutronsin the Inner Crust of Neutron Stars: New Life for Pulsar Glitch Models.PhRvL, 119 (6), 062701.[297] Baym, G., Pethick, C., Pines, D., Ruderman, M., Nov. 1969. Spin Upin Neutron Stars : The Future of the Vela Pulsar. Natur, 224 (5222),872–874.[298] Zhou, E. P., Lu, J. G., Tong, H., Xu, R. X., Sep. 2014. Two types ofglitches in a solid quark star model. MNRAS, 443 (3), 2705–2710.[299] Lai, X. Y., Yun, C. A., Lu, J. G., L¨u, G. L., Wang, Z. J., Xu, R. X.,May 2018. Pulsar glitches in a strangeon star model. MNRAS, 476 (3),3303–3309.[300] Ho, W. C. G., Espinoza, C. M., Antonopoulou, D., Andersson, N., Oct.2015. Pinning down the superfluid and measuring masses using pulsarglitches. SciA, 1 (9), e1500578–e1500578.[301] Pizzochero, P. M., Antonelli, M., Haskell, B., Seveso, S., Jul. 2017. Con-straints on pulsar masses from the maximum observed glitch. NatAs, 1,0134.[302] Ashton, G., Lasky, P. D., Graber, V., Palfreyman, J., Aug. 2019. Ro-tational evolution of the Vela pulsar during the 2016 glitch. NatAs, 3,1143–1148.[303] Ge, M. Y., Lu, F. J., Yan, L. L., Weng, S. S., Zhang, S. N., Wang, Q. D.,Wang, L. J., Li, Z. J., Zhang, W., Aug. 2019. The brightening of thepulsar wind nebula of PSR B0540-69 after its spin-down-rate transition.NatAs, 3, 1122–1127. 95304] Wei, J. Y., et al., Oct. 2016. The Deep and Transient Universe in theSVOM Era: New Challenges and Opportunities - Scientific prospects ofthe SVOM mission. preprint (arXiv:1610.06892).[305] ¨Ozel, F., Psaltis, D., Arzoumanian, Z., Morsink, S., Baub¨ock, M., Nov.2016. Measuring Neutron Star Radii via Pulse Profile Modeling withNICER. ApJ, 832 (1), 92.[306] Li, T. et al., Mar. 2018. Insight-HXMT observations of the first binaryneutron star merger GW170817. SCPMA, 61 (3), 31011.[307] Li, D., et al., Apr. 2018. FAST in Space: Considerations for a Multibeam,Multipurpose Survey Using China’s 500-m Aperture Spherical Radio Tele-scope (FAST). IMMag, 19 (3), 112–119.[308] Ray, P. S. a., Mar. 2019. STROBE-X: X-ray Timing and Spec-troscopy on Dynamical Timescales from Microseconds to Years. preprint(arXiv:1903.03035).[309] Watts, A. L. et al., Feb. 2019. Dense matter with eXTP. SCPMA, 62 (2),29503.[310] Liu, T., Lin, C.-Y., Song, C.-Y., Li, A., Nov. 2017. a comparison of Grav-itational Waves from Central Engines of Gamma-Ray Bursts: Neutrino-dominated Accretion Flows, Blandford-Znajek Mechanisms, and Millisec-ond Magnetars. ApJ, 850 (1), 30.[311] Zhou, X., Tong, H., Zhu, C., Wang, N., Dec. 2017. Dependence of pulsardeath line on the equation of state. MNRAS, 472 (2), 2403–2409.[312] Yasin, H., Sch¨afer, S., Arcones, A., Schwenk, A., Dec. 2018. Equation ofstate effects in core-collapse supernovae. preprint (arXiv:1812.02002).[313] Lau, S. Y., Leung, P. T., Lin, L. M., Jan. 2019. Two-layer compact starswith crystalline quark matter: Screening effect on the tidal deformability.PhRvD, 99 (2), 023018. 96314] Deng, Z.-L., Gao, Z.-F., Li, X.-D., Shao, Y., Mar. 2020. On the Formationof PSR J1640+2224: A Neutron Star Born Massive? ApJ, 892 (1), 4.[315] Shen, H., Ji, F., Hu, J., Sumiyoshi, K., Mar. 2020. Effects of SymmetryEnergy on the Equation of State for Simulations of Core-collapse Super-novae and Neutron-star Mergers. ApJ, 891 (2), 148.[316] Zhang, N.-B., Li, B.-A., Xu, J., Jun. 2018. Combined Constraints on theEquation of State of Dense Neutron-rich Matter from Terrestrial NuclearExperiments and Observations of Neutron Stars. ApJ, 859 (2), 90.[317] Dietrich, T., Coughlin, M. W., Pang, P. T. H., Bulla, M., Heinzel, J., Issa,L., Tews, I., Antier, S., Feb. 2020. New Constraints on the SupranuclearEquation of State and the Hubble Constant from Nuclear Physics – Multi-Messenger Astronomy. preprint (arXiv:2002.11355).[318] Baldo, M., Burgio, G. F., Centelles, M., Sharma, B. K., Vi˜nas, X., Sep.2014. From the crust to the core of neutron stars on a microscopic basis.PAN, 77 (9), 1157–1165.[319] Centelles, M., Roca-Maza, X., Vi˜nas, X., Warda, M., Mar. 2009. NuclearSymmetry Energy Probed by Neutron Skin Thickness of Nuclei. PhRvL,102 (12), 122502.[320] Chrien, R. E., Dover, C. B., Jan. 1989. Nuclear systems with strangeness.ARNPS, 39 (39), 113–150.[321] Fraga, E. S., Kurkela, A., Vuorinen, A., Feb. 2014. Interacting QuarkMatter Equation of State for Compact Stars. ApJ, 781 (2), L25.[322] Frederico, T., Carlson, B. V., Rego, R. A., Hussein, M. S., Mar. 1989.Quark structure of the nucleon and quantum hadrodynamics. JPhG,15 (3), 297–302.[323] Gamba, R., Read, J. S., Wade, L. E., Jan. 2020. The impact of the crustequation of state on the analysis of GW170817. CQGra, 37 (2), 025008.97324] Gittins, F., Andersson, N., Pereira, J. P., Mar. 2020. Tidal deformationsof neutron stars with elastic crusts. preprint (arXiv:2003.05449).[325] Kl¨ahn, T., Fischer, T., Sep. 2015. Vector Interaction Enhanced Bag Modelfor Astrophysical Applications. ApJ, 810 (2), 134.[326] Kojo, T., Powell, P. D., Song, Y., Baym, G., Feb. 2015. PhenomenologicalQCD equation of state for massive neutron stars. PhRvD, 91 (4), 045003.[327] Lattimer, J. M., Lim, Y., Jul. 2013. Constraining the Symmetry Parame-ters of the Nuclear Interaction. ApJ, 771 (1), 51.[328] Lattimer, J. M., Prakash, M., Mar. 2001. Neutron Star Structure and theEquation of State. ApJ, 550 (1), 426–442.[329] Lattimer, J. M., Steiner, A. W., Feb. 2014. Constraints on the symmetryenergy using the mass-radius relation of neutron stars. EPJA, 50, 40.[330] Lattimer, J. M., Steiner, A. W., Apr. 2014. Neutron Star Masses andRadii from Quiescent Low-mass X-Ray Binaries. ApJ, 784 (2), 123.[331] Li, A., Wang, R., Aug. 2018. Pulsar glitch and nuclear EoS: Applicabilityof superfluid model. In: Weltevrede, P., Perera, B. B. P., Preston, L. L.,Sanidas, S. (Eds.), IAUS, 337, 360.[332] Li, B.-A., Han, X., Nov. 2013. Constraining the neutron-proton effectivemass splitting using empirical constraints on the density dependence ofnuclear symmetry energy around normal density. PhLB, 727 (1-3), 276–281. 9.[333] Li, B.-A., Steiner, A. W., Nov. 2006. Constraining the radii of neutronstars with terrestrial nuclear laboratory data. PhLB, 642 (5-6), 436–440.[334] Li, Z. H., Schulze, H. J., Aug. 2008. Neutron star structure with modernnucleonic three-body forces. PhRvC, 78 (2), 028801.98335] Link, B., Jul. 2014. Thermally Activated Post-glitch Response of the Neu-tron Star Inner Crust and Core. I. Theory. ApJ, 789 (2), 141.[336] Margaritis, C., Koliogiannis, P. S., Moustakidis, C. C., Feb. 2020. Speed ofsound constraints on maximally rotating neutron stars. PhRvD, 101 (4),043023.[337] M¨uller, H., Serot, B. D., Feb. 1996. Relativistic mean-field theory and thehigh-density nuclear equation of state. NuPhA, 606, 508–537.[338] Page, D., Lattimer, J. M., Prakash, M., Steiner, A. W., Dec. 2004. Mini-mal Cooling of Neutron Stars: A New Paradigm. ApJS, 155 (2), 623–650.[339] Patrignani, C. et al., Oct. 2016. Review of Particle Physics. ChPhC,40 (10), 100001.[340] Pearson, J. M., Chamel, N., Goriely, S., Ducoin, C., Jun. 2012. Innercrust of neutron stars with mass-fitted Skyrme functionals. PhRvC, 85 (6),065803.[341] Perego, A., Radice, D., Bernuzzi, S., Dec. 2017. AT 2017gfo: AnAnisotropic and Three-component Kilonova Counterpart of GW170817.ApJ, 850 (2), L37.[342] Raaijmakers, G., et al., Dec. 2019. A Nicer View of PSR J0030+0451:Implications for the Dense Matter Equation of State. ApJ, 887 (1), L22.[343] Samajdar, A., Dietrich, T., Feb. 2020. Constructing Love-Q-Relationswith Gravitational Wave Detections. preprint (arXiv:2002.07918).[344] Schaeffer, R., Zdunik, L., Haensel, P., Sep. 1983. Phase transitions instellar cores. I - Equilibrium configurations. A&A, 126 (1), 121–145.[345] Schaffner-Bielich, J., Gal, A., Sep. 2000. Properties of strange hadronicmatter in bulk and in finite systems. PhRvC, 62 (3), 034311.99346] Schaffner-Bielich, J., Hanauske, M., St¨ocker, H., Greiner, W., Jan. 2002.Phase Transition to Hyperon Matter in Neutron Stars. PhRvL, 89 (17),171101.[347] Soares-Santos, M. et al., Oct. 2017. The Electromagnetic Counterpart ofthe Binary Neutron Star Merger LIGO/Virgo GW170817. I. Discovery ofthe Optical Counterpart Using the Dark Energy Camera. ApJ, 848 (2),L16.[348] Vida˜na, I., Polls, A., Ramos, A., Engvik, L., Hjorth-Jensen, M., Sep.2000. Hyperon-hyperon interactions and properties of neutron star matter.PhRvC, 62 (3), 035801.[349] Watts, A. L., et al., Apr. 2016. Colloquium: Measuring the neutron starequation of state using x-ray timing. RvMP, 88 (2), 021001.[350] Wei, J. B., Burgio, G. F., Schulze, H. J., Apr. 2019. Neutron star coolingwith microscopic equations of state. MNRAS, 484 (4), 5162–5169.[351] Weissenborn, S., Chatterjee, D., Schaffner-Bielich, J., Jun. 2012. Hyper-ons and massive neutron stars: Vector repulsion and SU(3) symmetry.PhRvC, 85 (6), 065802.[352] Weissenborn, S., Sagert, I., Pagliara, G., Hempel, M., Schaffner-Bielich,J., Oct. 2011. Quark Matter in Massive Compact Stars. ApJ, 740 (1), L14.[353] Whittenbury, D. L., Carroll, J. D., Thomas, A. W., Tsushima, K., Stone,J. R., Jun. 2014. Quark-meson coupling model, nuclear matter constraints,and neutron star properties. PhRvC, 89 (6), 065801.[354] Yakovlev, D. G., Pethick, C. J., Sep. 2004. Neutron Star Cooling. ARA&A,42 (1), 169–210.[355] Zhang, K., Hirayama, T., Luo, L.-W., Lin, F.-L., Feb. 2020. Compact starof holographic nuclear matter and GW170817. PhLB., 801, 135176.100356] Zhang, K., Huang, G.-Z., Lin, F.-L., Feb. 2020. GW170817 and GW190425as Hybrid Stars of Dark and Nuclear Matters. preprint (arXiv:2002.10961).[357] Zhang, Z., Chen, L.-W., Sep. 2015. Electric dipole polarizability in Pbas a probe of the symmetry energy and neutron matter around ρ0