Constraints on Phase Transitions in Neutron Star Matter

1. Introduction

Within recent decades the observational data base for neutron stars has expanded substantially. Its detailed analysis sharpened the empirical constraints on the equation of state (EoS) of matter deep inside the cores of these extreme objects. In particular, the existence of heavy neutron stars with masses around and above two solar masses implies that the EoS, i.e. pressure $P\left(\epsilon \right)$ as function of energy density $\epsilon$, must be sufficiently stiff in order to support such massive compact stars against gravitational collapse. Some previously discussed simple forms of exotic matter could therefore be excluded when their corresponding EoS’s turned out to be too soft.

The quest for a possible phase transition in cold and dense baryonic matter, from hadronic to quark degrees of freedom, has been a topic of prime interest for a long time. In this context a summary will be given of the present state of knowledge about the neutron star EoS as inferred directly from the existing observational data. Of particular interest is the speed of sound in the cores of neutron stars. Its behaviour as a function of energy density is a sensitive indicator for phase transitions or continuous changes (crossovers) of degrees of freedom in the star’s composition.

Primary sources of information for the inference of the EoS are neutron star masses deduced from Shapiro delay measurements [1,2,3,4,5,6], combined determinations of masses and radii inferred from X-ray data detected by the NICER telescope [7,8,9,10,11], and gravitational wave signals from binary neutron star mergers observed by the LIGO and Virgo Collaborations [12,13,14]. These data have served as basic input for a large variety of Bayesian inference analyses [15,16,17,18,19,20,21,22,23,24,25,26,27,28,29] in search for a range of data-compatible equations of state with controlled uncertainties. Low-density constraints from nuclear physics around the equilibrium density of nuclear matter, ${\rho }_0=0.16$fm⁻³, are commonly introduced by referring to results of chiral effective field theory (ChEFT) [30]. At asymptotically high densities, perturbative QCD (pQCD) is applicable and provides a constraint to be matched in extrapolations far beyond the conditions realized in neutron star cores [31,32]. The Bayesian inference studies have recently been further extended [27], given the supplementary observational data reported in [33] (the black widow pulsar PSR J0952-0607 with a mass $M=2.35\pm 0.17M_{\odot }$ in solar mass units).

This report presents a survey of results obtained in the detailed inference analysis of Ref.[27] and discusses possible interpretations. A central topic is the evaluation of constraints and possible evidence for a first-order phase transition within the posterior uncertainties of the EoS deduced from the existing empirical neutron star data. Later parts of this manuscript refer to ongoing considerations about the structure and composition of dense matter in neutron star cores. In view of the constraints emerging from the analysis of the observational data, especially of the heaviest known neutron stars, we discuss the quest for a chiral phase transition, quark-hadron continuity and its realisation as a soft crossover, a Fermi liquid approach to neutron star matter, and other related issues.

2. Neutron star matter equation of state

2.1. Observational constraints

2.1.1. Neutron star masses and radii

Information about the masses of the heaviest neutron stars derives primarily from precise Shapiro time delay measurements of pulsars orbiting in binary systems with white dwarf companions. Three such massive objects, PSR J1614–2230 [1,3,4], PSR J0348+0432 [2] and PSR J0740+6620 [5,6], have been established in the past:

$$\begin{array}{ccc}& \mathrm{PSR}\mathrm{J}1614\text{–}2230& M=1.908\pm 0.016M_{\odot },\hfill \end{array}$$

(1)

$$\begin{array}{ccc}& \mathrm{PSR}\mathrm{J}0348+0432& M=2.01\pm 0.04M_{\odot },\hfill \end{array}$$

(2)

$$\begin{array}{ccc}& \mathrm{PSR}\mathrm{J}0740+6620& M=2.08\pm 0.07M_{\odot }.\hfill \end{array}$$

(3)

Recently the heaviest neutron star observed so far has been reported [33], a black widow pulsar,

$$\begin{array}{cc}\hfill \mathrm{PSR}\mathrm{J}0952-0607& M=2.35\pm 0.17M_{\odot }.\end{array}$$

(4)

This result was obtained based on the detection of the black widow’s companion star. PSR J0952-0607 is also one of the fastest rotating pulsars. With a spin period of 1.4 ms it requires radius-dependent corrections for rotational effects as described in [27]. For example, the rotating mass of a $2.35M_{\odot }$ star with an assumed radius of $R=12$ km decreases by 3% to an equivalent non-rotating mass of about 2.28 $M_{\odot }$.

Radii of neutron stars together with their masses can be inferred from X-ray profiles of rotating hot spot patterns measured with the NICER telescope. Two neutron stars have been investigated in this way [7,9]:

$$\begin{array}{ccc}\hfill \mathrm{PSR}\mathrm{J}0030+0451& M=1.{34}_{-0.16}^{+0.15}{M}_{\odot },& R=12.{71}_{-1.19}^{+1.14}\mathrm{km},\hfill \end{array}$$

(5)

$$\begin{array}{ccc}\hfill \mathrm{PSR}\mathrm{J}0740+6620& M=2.{072}_{-0.066}^{+0.067}{M}_{\odot },& R=12.{39}_{-0.98}^{+1.30}\mathrm{km}.\hfill \end{array}$$

(6)

All of the mass and radius data are included as input in setting constraints for the neutron star EoS. In the case of PSR J0740+6620 the analysis of the NICER measurement includes the Shapiro delay result. Only the data from NICER have therefore been used in the Bayesian inference data set. An alternative analysis of the NICER data has also been performed by a second independent group [8,10]. Both results are mutually compatible within their uncertainties.

2.1.2. Binary neutron star mergers and tidal deformabillities

Gravitational wave signals produced by the merger of two neutron stars in a binary have been detected by the LIGO and Virgo collaborations [12,13,14]. These signals are interpreted using theoretical waveform models that depend on the mass ratio of the two neutron stars, $M_2/M_1$, and a mass-weighted combination of their tidal deformabilities. Based on information from the GW170817 event the (dimensionless) tidal deformability for a 1.4 solar mass neutron star has been deduced [13]:

$$\begin{array}{c}\hfill {\Lambda }_{1.4}={190}_{-120}^{+390}.\end{array}$$

(7)

The GW170817 event was further investigated together with electromagnetic signals. In one such analysis [35] the following masses and tidal deformabilities were reported for the individual neutron stars in the binary :

$$\begin{array}{ccc}& M_1=1.{46}_{-0.09}^{+0.13}{M}_{\odot },& {\Lambda }_1={255}_{-171}^{+416},\hfill \\ & M_2=1.{26}_{-0.12}^{+0.09}{M}_{\odot },& {\Lambda }_2={661}_{-375}^{+858}.\hfill \end{array}$$

(8)

The event GW190425 which could have also been a binary neutron star merger, was also included in the inference analysis [14].

2.2. Inference of sound velocity and EoS in neutron stars

Based on the data in the previous subsection the squared sound velocity and the equation of state,

$$\begin{array}{c}\hfill c_s^2\left(\epsilon \right)=\frac{\partial P\left(\epsilon \right)}{\partial \epsilon }\mathrm{and}P\left(\epsilon \right)={\int }_0^{\epsilon }d{\epsilon }^{\prime }{c}_s^2\left({\epsilon }^{\prime }\right),\end{array}$$

(9)

of neutron star matter can be inferred using Bayesian methods. In our recent work [25,27] the ChEFT constraint at low baryon density is used in the Bayes inference procedure as a likelihood (not as a prior) within a conservative window of baryon densities [34], $\rho \lesssim 1.3{\rho }_0$. The equations of state are interpolated to pQCD results at asymptotically high densities [32].

The result of the inferred sound velocity in Figure 1 shows, within the $68\%$ credible bands, a rapid increase of $c_s^2$ beyond the conformal limit of 1/3 at energy densities relevant for a wide range of neutron stars with masses $M\sim$ 1.4 - 2.3 $M_{\odot }$. A quantification in terms of a corresponding Bayes factor [17,25,27] demonstrates indeed extreme evidence for $c_s^2$ exceeding the conformal bound inside neutron stars.

Figure 2 presents the inferred posterior bands for the EoS. Notably, the additional inclusion of the heaviest pulsar PSR J0952-0607 in the data set leads to a more rapidly rising pressure as compared to previous inference results which do not incorporate this data. Within the range of energy densities 0.5 - 0.8 GeV/fm³ reached in the core of a 2.3 solar mass star and corresponding to baryon densities between 3.0${\rho }_0$ and 4.6${\rho }_0$, the pressure even exceeds that of the time-honored APR EoS [36].

A further quantity of interest is the chemical potential associated with the conserved baryon number,

$$\begin{array}{c}\hfill \mu \left(\rho \right)=\frac{\partial \epsilon \left(\rho \right)}{\partial \rho }=\frac{P+\epsilon }{\rho },\end{array}$$

(10)

where the latter equality follows from the Gibbs-Duhem relation at zero temperature, $P+\epsilon =\mu \rho$. This baryon chemical potential can be computed for any given equation of state within the uncertainty bands in Figure 2. The resulting $\mu \left(\rho \right)$ bands in Figure 3 show a steep rise at baryon densities $\rho >2{\rho }_0$, indicating strongly repulsive correlations at work. The agnostic approach for the inference of the sound speed on which these results are based does not permit to distinguish between different species of constituents inside neutron stars, so $\mu$ represents the total chemical potential from all active degrees of freedom carrying baryon number: $\mu ={\sum }_i{x}_i{\mu }_i$, where $x_i={\rho }_i/\rho$ is the fraction corresponding to each baryonic species i. But irrespective of whether these species are nucleons, other baryonic composites or quarks, the empirically inferred behaviour of $\mu$ demonstrates that these degrees of freedom must be correlated by strongly repulsive forces as the density increases, in order to build up the necessary high pressure in the cores of neutron stars to support two solar masses and beyond. We return to this discussion again at a later stage.

The approach to conformality in strongly interacting matter at high baryon density has been an important theme in several recent investigations [38,39,40,41]. A key quantity is the trace anomaly measure,

$$\begin{array}{c}\hfill \Delta =\frac{g_{\mu \nu }{T}^{\mu \nu }}{3\epsilon }=\frac{1}{3}-\frac{P\left(\epsilon \right)}{\epsilon }.\end{array}$$

(11)

Conformal matter is characterized by an equation of state $P=\epsilon /3$, so the signature for its occurance is $\Delta \to 0$ over an extended energy density range. The posterior credible bands for the trace anomaly measure deduced from the data-based EoS are shown in Figure 4. Starting from $\Delta =1/3$ at zero energy density, the median of $\Delta$ turns negative at $\epsilon \sim 0.7{\mathrm{GeV}/\mathrm{fm}}^3$, entering a high-pressure domain with $P>\epsilon /3$. This crossing appears within the range of energy densities possibly realized in the cores of the heaviest neutron stars. At much higher energy densities beyond those displayed in Figure 3, pQCD then implies a switch back to positive $\Delta$ before the asymptotic limit $\Delta \to 0$ is approached.

2.3. Selected neutron star properties

Numerically solving the Tolman-Oppenheimer-Volkov (TOV) equations leads to posterior credible bands for the mass-radius relation of neutron stars as displayed in Figure 5. Notably, the median of $M\left(R\right)$ suggests an almost constant radius R, independent of M.

The density profiles $\rho \left(r\right)$ computed for neutron stars in the mass range 1.4 - 2.3 $M_{\odot }$ reach central baryon densities, ${\rho }_c=\rho (r=0)$, that are systematically below $5{\rho }_0$ in the cores of even the heaviest stars. For example,

$$\begin{array}{c}\hfill {\rho }_c\left(1.4M_{\odot }\right)=(2.6\pm 0.4){\rho }_0\mathrm{and}{\rho }_c\left(2.3M_{\odot }\right)=(3.8\pm 0.8){\rho }_0,\end{array}$$

(12)

at the 68% level [27]. In a baryonic picture of neutron stars, this implies that the average distance between baryons even in the highly compressed star center exceeds 1 fm, more than twice the characteristic short-range hard-core distance of 1/2 fm in the nucleon-nucleon interaction. We return to this point at a later stage.

For each given EoS, the TOV equations can be solved in combination with equations for the tidal deformability $\Lambda$. This leads to posterior probability bands, $\Lambda \left(M\right)$, as a function of neutron star mass as shown in Figure 6. The large uncertainties of the empirical values deduced from the analysis of the gravitational wave events do so far not allow to constrain the tidal deformability bands to a stronger degree.

3. Constraints on phase transitions in neutron stars

The nature and location of the transition from dense nuclear or neutron matter to quark and gluon degrees of freedom is still largely unknown. Lattice computations have been successful in exploring the QCD phase structure at high temperatures and vanishing baryon chemical potential, establishing that the transition is a continuous crossover at a temperature $T_c\simeq 155$ MeV. However, at nonzero baryon densities the notorious sign problem in the Euclidean action prevents a similarly systematic approach. As a consequence, studies of the phase diagram at low temperatures and high baryon densities are primarily based on models with varying assumptions about the relevant active degrees of freedom [42].

A key element of this discussion is the question about a possible phase transition from spontaneously broken to restored chiral symmetry in QCD. At high baryon densities beyond $\rho >2{\rho }_0$ various hypotheses have been explored, ranging from a chiral first-order phase transition to continuous hadron-quark crossovers [43,44,45,46,47,48]. The behaviour of the sound speed $c_s$ is a prime indicator for phase transitions or crossovers. This is schematically illustrated in Figure 7 which shows typical patterns of such transitions as they would show up in the sound speed. For example, a first-order phase transition with Maxwell construction is characterized by a region of constant pressure over an interval of density (or energy density) in which two phases coexist. Between the lower and upper end points of this phase coexistence interval the squared sound velocity $c_s^2=\partial P/\partial \epsilon$ jumps to zero and back. In a first-order phase transition with Gibbs construction the pressure in the mixed phase is not constant but $\partial P/\partial \epsilon$ still changes discontinuously [50]. A crossover has no such discontinuities but features a pronounced maximum in $c_s^2$. Purely nucleonic scenarios lead instead to a continuously rising sound speed.

3.1. Evidence against a very low squared sound speed in neutron stars

One indication for the possible appearance of a first-order phase transition would be a very low squared sound velocity $(c_s^2\le 0.1)$ within the range of energy densities relevant for neutron stars. Let us introduce $c_{s,min}$ as the minimum speed of sound following a maximum at lower densities. Given the inferred bands of probability distributions for $c_s^2$ in Section 2.2 and Figure 1, one can quantify the evidence for two competing scenarios: minimum $c_{s,min}^2>0.1$ versus $c_{s,min}^2\le 0.1$, in terms of corresponding Bayes factors, ${\mathcal{B}}_{c_{s,min}^2\le 0.1}^{c_{s,min}^2>0.1}$. A detailed analysis of such Bayes factors has been carried out in [27] for the minimum speed of sound up to different neutron star maximum masses, $M_{max}$. Bayes factors exceeding 10 or 100 indicate strong or extreme evidence, respectively, for the scenario where the minimum $c_s^2$ remains larger than 0.1. This is the case for all masses $M_{max}\le 2.1M_{\odot }$ as demonstrated in Table 1. Only for even heavier neutron stars does this evidence turn out to be lower. It is worth noting that the inclusion of the heavy black widow pulsar PSR J0952-0607 plays an important role strengthening the evidence against $c_s^2\le 0.1$ further in comparison with previous analyses which do not inlude this object [25]. Squared sound speeds $c_s^2>0.1$ in the cores of neutron stars with $M_{max}\le 2.0M_{\odot }$ were also found in [20,22] without inclusion of the black widow pulsar.

These limiting bounds on very low speeds of sound suggest that the appearance of a strong first-order phase transition with Maxwell construction is unlikely for most of the known neutron stars. An additional feature of a Maxwell-constructed first-order transition is an extended phase coexistence region. This will be examined more closely in the following subsection.

3.2. Evidence against a strong first-order phase transition in the cores of neutron stars

The relatively moderate baryon densities inferred in the cores of neutron stars together with the evidence against very small sound speeds can be complemented by an additional study that is specific to first-order phase transitions with Maxwell construction. Such phase transitions are characterized by a domain of phase coexistence that extends over a certain range of baryon densities, $\Delta \rho$. The width of this domain, $\Delta \rho /\rho$ (with $\rho$ the density at which the coexistence interval starts), is a measure of the strength of the phase transition, i.e. the magnitude of surface tension between the two coexisting phases. We refer to a `strong’ first-order phase transition if the width of the mixed phase is $\Delta \rho /\rho >1$. One example is the liquid-gas phase transition in symmetric nuclear matter at temperatures $T<15$ MeV [51,71]. In contrast a `weak’ first-order phase transition has $\Delta \rho /\rho$ small compared to unity.

The posterior credible bands inferred from neutron star data, as displayed in Figure 2, permit a systematic study of the maximum possible phase coexistence widths over the range of relevant energy densities. Starting from $P\left(\epsilon \right)$ the Gibbs-Duhem relation is used to reexpress pressure as a function of baryon density, $P\left(\rho \right)$. In this way posterior credible bands for $P\left(\rho \right)$ can be derived. These constrain the maximum possible widths of mixed phase domains, ${(\Delta \rho /\rho )}_{\max }$. An example is shown in Figure 8. In the analysis [27] it turns out that these maximum possible phase coexistence regions are narrow: ${(\Delta \rho /\rho )}_{\max }\simeq 0.2$ at the 68% level. Even within the 95% credible bands this width does not exceed 0.3. In fact one finds that ${(\Delta \rho /\rho )}_{\max }$ stays nearly constant as a function of baryon density $\rho$ (taken at the starting point of the mixed-phase interval) over the whole region of densities $\rho \simeq 2-5{\rho }_0$ relevant for neutron stars.

3.3. Intermediate summary

The conclusion drawn in the reported analysis is therefore that only a weak first-order phase transition with a mixed-phase width $\Delta \rho /\rho \le 0.2-0.3$ in a Maxwell construction can be realized inside neutron stars within the posterior credible bands of the inferred equation of state. A strong first-order phase transition is thus very unlikely to occur [27]. It is worth pointing out that this conclusion is entirely based on the analysis of observational data, independent of any specific EoS model. As will be discussed later, this rules out some relativistic mean-field models which suggested the occurance of a strong (chiral) first-order phase transition already at densities below those encounterd in neutron star cores. On the other hand the empirical uncertainty bands still leave room for a continuous transition such as a hadrons-to-quarks crossover [45,46,47,48]. The possiblity of a purely baryonic equation of state with monotonically rising speed of sound is also not excluded [25].

4. Phenomenology and models

The previous sections have shown how present and forthcoming neutron star data (masses, radii, tidal deformabilities) provide constraints for the equation of state and for possible phase transitions in cold, dense baryonic matter. The enhanced stiffness of the EoS at high baryon densities is a necessary feature: sufficiently high pressures are required to support neutron stars with masses around and beyond $2M_{\odot }$ at relatively large radii ($R\simeq 11-13$ km). The central core densities are not as extreme as was previously imagined. Even in the core of a 2.3 $M_{\odot }$ neutron star the baryon density at its center stays below five times the density of equilibrium nuclear matter (at 68% credible level, see Eq. (12)). The upper border of the 95% credible interval of the inferred central density ${\rho }_c\left(2.3M_{\odot }\right)/{\rho }_0=3.8_{-1.3}^{+1.6}$[27] lies slightly above $5{\rho }_0$, but with very low probability. This observation has consequences for the interpretation of the possible structure and composition of neutron star matter that we discuss in the following.

4.1. Reminder of low-energy nucleon structure and a two-scales scenario

Spontaneouly broken chiral symmetry as the long-wavelength manifestation of QCD governs the low-energy structure and dynamics of hadrons including nucleons and pions. As chiral Nambu-Goldstone bosons, pions play a distinguished role in this context. Models based on chiral symmetry often view the nucleon as a complex system of two scales [52]: a compact hard core that hosts the three valence quarks and thus incloses the baryon number, and a surrounding quark-antiquark cloud in which pions figure prominently as the `soft’ degrees of freedom.

Such a two-scales scenario is manifest in empirical form factors and sizes of the nucleon. Consider for example the proton and neutron electromagnetic form factors and their slopes at zero momentum transfer which determine the corresponding mean squared radii. The empirical r.m.s. proton charge radius, ${\langle r_p^2\rangle }^{1/2}=0.840\pm 0.003$ fm, combined with the slope of the neutron electric form factor, $\langle r_n^2\rangle =-0.105\pm 0.006$ fm², gives the isoscalar and isovector mean-square radii of the nucleon, $\langle r_{S,V}^2\rangle =\langle r_p^2\rangle \pm \langle r_n^2\rangle$ with resulting values [53]:

$$\begin{array}{c}\hfill \sqrt{\langle r_S^2\rangle }\simeq 0.78\mathrm{fm},\sqrt{\langle r_V^2\rangle }\simeq 0.90\mathrm{fm}.\end{array}$$

(13)

Each of the nucleon form factors related to a current with index i has a representation in terms of a dispersion relation,

$$\begin{array}{c}\hfill G_i\left(q^2\right)=G_i\left(0\right)+\frac{q^2}{\pi }{\int }_{t_0}^{\infty }dt\frac{\mathrm{Im}{G}_i\left(t\right)}{t(t-q^2-iϵ)}.\end{array}$$

(14)

Mean-square radii are given as:

$$\begin{array}{c}\hfill \langle r_i^2\rangle =6\frac{d{G}_i\left(q^2\right)}{d{q}^2}{|}_{q^2=0}=\frac{6}{\pi }{\int }_{t_0}^{\infty }\frac{dt}{t^2}{S}_i\left(t\right),\end{array}$$

(15)

where the distribution $S_i\left(t\right)=\mathrm{Im}{G}_i\left(t\right)$ represents the spectrum of intermediate hadronic states through which the external probing field couples to the respective nucleon current. For example, the isovector charge radius reflects the interacting two-pion cloud of the nucleon governed by the $\rho$ meson and a low-mass tail extending down to the $\pi \pi$ threshold, $t_0=4m_{\pi }^2$. The isoscalar charge radius is related to the three-pion spectrum that is strongly dominated by the narrow $\omega$ meson [54] and starts at $t_0=9m_{\pi }^2$. The isoscalar charge form factor of the nucleon, $G_S^E\left(q^2\right)$ (with $G_S^E\left(0\right)=1$), is particularly suitable to discuss a delineation between `core’ and `cloud’ parts of the nucleon [55,56]. The vector meson dominance principle implies, in its simplest version, a representation of the form:

$$\begin{array}{c}\hfill G_S^E\left(q^2\right)=\frac{F_B\left(q^2\right)}{1+|q^2|/m_{\omega }^2}.\end{array}$$

(16)

The form factor $F_B\left(q^2\right)$ of the baryon number distribution in the nucleon core acts as a source for the $\omega$ field that propagates with its mass $m_{\omega }$. Introducing the mean-square radius of the baryon core, $\langle r_B^2\rangle =6\frac{d{F}_B\left(q^2\right)}{d{q}^2}{|}_{q^2=0}$, the mean-square isoscalar charge radius becomes

$$\begin{array}{c}\hfill \langle r_S^2\rangle =\langle r_B^2\rangle +\frac{6}{m_{\omega }^2}.\end{array}$$

(17)

Using $m_{\omega }=783$ MeV the estimated baryonic core radius is:

$$\begin{array}{c}\hfill \sqrt{\langle r_B^2\rangle }\simeq 0.47\mathrm{fm}.\end{array}$$

(18)

A nucleon core size of about 1/2 fm is characteristic of chiral `core + cloud’ models. It also keeps up in more detailed treatments of the spectral distributions governing the nucleon form factors [53]. Inclusion of additional $\varphi$ meson and $\pi \rho$ continuum contributions in the spectral function of $G_S^E\left(q^2\right)$ move the core radius to just slightly larger values.

Consider as another example the form factor associated with the axial vector current of the nucleon. The corresponding mean-square axial radius deduced from neutrino-deuteron scattering data is reported as [57]:

$$\begin{array}{c}\hfill \langle r_A^2\rangle =(0.46\pm 0.22){\mathrm{fm}}^2.\end{array}$$

(19)

A schematic axial vector dominance picture would assign a dominant part of the `cloud’ contribution to this form factor through the spectrum of the $a_1$ meson with its large width. If an approximate scale of this `cloud’ part is identified with $\delta \langle r_A^2\rangle \sim 6/m_a^2$ using the physical $a_1$ mass, $m_a\simeq 1.23$ GeV, one finds

$$\begin{array}{c}\hfill {\sqrt{\langle r_A^2\rangle }}_{\mathrm{core}}={[\langle r_A^2\rangle -\delta \langle r_A^2\rangle ]}^{1/2}\simeq 0.55\mathrm{fm}\end{array}$$

(20)

with an estimated uncertainty of about 25%.

Yet another interesting piece of information is the mass radius of the proton deduced from $J/\Psi$ photoproduction data [58]. It involves the form factor of the trace, $T_{\mu }^{\mu }$, of the nucleon’s energy-momentum tensor and is supposed to be dominated by gluon dynamics at the center of the nucleon. Low-mass (e.g. two-pion) components are suppressed. The result quoted in [58],

$$\begin{array}{c}\hfill {\langle r^2\rangle }_{\mathrm{mass}}^{1/2}=(0.55\pm 0.03)\mathrm{fm},\end{array}$$

(21)

is once again remarkably close to an assumed `core’ size scale of $\sim 1/2$ fm.

These empirical considerations motivate a picture of the nucleon as containing a compact `hard’ core in which the valence quarks and their baryon number are confined, and a `soft’ surface of quark-antiquark pairs forming a mesonic cloud. This structure has implications for the behaviour of nucleons in dense baryonic matter. With a core size $R_{\mathrm{core}}\sim 1/2$ fm and a cloud range given e.g. by the proton charge radius, $R_{\mathrm{cloud}}\sim 0.84$ fm, there is a significant separation of volume scales in vacuum: ${(R_{\mathrm{cloud}}/R_{\mathrm{core}})}^3\sim 5$.

This scale separation is expected to increase further in dense baryonic matter, for the following reasons. The properties of the soft multi-pion cloud are closely tied to spontaneously broken chiral symmetry and the approximate Nambu-Goldstone boson nature of the pion. The size of this cloud is expected to increase with density, along with the decreasing in-medium pion decay constant, $f_{\pi }^{*}\left(\rho \right)$, which acts as a chiral order parameter. The baryonic core, on the other hand, is governed by gluon dynamics, without a leading connection to chiral symmetry in QCD. This core is therefore expected to be quite stable against changes from increasing density up until the compact hard cores begin to touch and overlap.

What arises in this way is a two-scales scenario for dense baryonic matter as described in [47] and sketched in Figure 9. At $\rho \simeq {\rho }_0$ the tails of the meson clouds of nucleon pairs overlap, resulting in two-body exchange forces. As the average distance between nucleons decreases with increasing density, around $\rho \gtrsim 2-3{\rho }_0$, the soft clouds of $q\overline{q}$ pairs start to be delocalized. Their mobility expands over larger distances in a percolation process involving larger numbers of nucleons. In conventional nuclear physics terms, this corresponds to the emergence of many-body forces, the strength of which grows with increasing density. At that stage the baryonic cores are still separated but subject to increasingly repulsive Pauli principle effects. The cores begin to touch and overlap at average nucleon-nucleon distances $d\lesssim 1$ fm, corresponding to densities $\rho \gtrsim 5{\rho }_0$. Further compression of baryonic matter would still have to overcome the strongly repulsive hard core in the nucleon-nucleon interaction. Recalling the inferred credible intervals of the central densities in heavy neutron stars, Eq. (12), one concludes that a phase transition to valence quark deconfinement does not seem likely in a two-scales picture and under the conditions in the cores of neutron star, unless they are extremely heavy.

An interesting and closely related result emerges from a detailed analysis of $y-$scaling in electron-nucleus scattering at large momentum transfers ($|\overrightarrow{q}|\gtrsim 1$ GeV) and low energy transfers [59]. The persistently observed $y-$scaling under these kinematical conditions implies that the electrons still scatter from strongly correlated pairs of nucleons, rather than quarks, at short distances. The conclusion drawn in [59] is that at local densities as large as five times ${\rho }_0$, nuclear matter still appears to behave as a collection of nucleons.

4.2. Quark-hadron continuity and crossover

While a strong first-order phase transition in neutron star matter seems unlikely based on current empirical observations, a continuous crossover from hadrons to quarks is still possible within the present data-driven constraints. Such a scenario is realised for example in the QHC21 equation of state [48]. It features a smooth interpolation between low and high densities from nuclear to quark matter regimes. The quark matter phase is described by a three-flavor Nambu–Jona-Lasinio (NJL) model including pairing (diquark) degrees of freedom and a strongly repulsive vector coupling between quarks. In order to build up the necessary pressure at high densities, this vector coupling must be comparable to or larger than the strength of the standard scalar-pseudoscalar interaction in the NJL model. The density at which the interpolated turnover to quark matter takes place in the QHC21 EoS is chosen as $\rho \gtrsim 3.5{\rho }_0$, within the range of central densities that can be reached in $M\gtrsim 2M_{\odot }$ neutron stars. A crossover from hadrons to quarks could also involve an intermediate phase of quarkyonic matter [45,49].

4.3. Chiral symmetry restoration in dense matter: from first-order phase transition to crossover

The quest for chiral symmetry restoration at high baryon densities, a transition from the spontaneously broken Nambu-Goldstone realisation to the unbroken Wigner-Weyl mode, has been a persistent issue for a long time. A possible first-order chiral phase transition and the existence of a corresponding critical end point in the QCD phase diagram have always been topics of prime interest [60,61,62]. Early hypotheses concerning the occurance of a first-order phase transition were frequently based on Nambu–Jona-Lasinio type (NJL) models in mean-field approximation [63,64,65], later extended by incorporating confinement aspects through added Polyakov loop degrees of freedom [66,67,68].

The empirical constraints on strong first-order phase transitions in dense neutron star matter as described in Section 3.2 are strikingly at variance with previous mean-field (e.g. NJL model) predictions. These suggested that a first-order chiral phase transition should already appear at relatively low baryon densities, $\rho \sim 2-3{\rho }_0$. A possible explanation for this discrepancy can be found in [70,71] where a chiral nucleon-meson (ChNM) field theory was used to explore the effects of fluctuations beyond mean-field (MF) approximation. The starting point is a relativistic chiral Lagrangian, $\mathcal{L}(\Psi ;\pi ,\sigma ;v_{\mu })$, shaped around a linear sigma model with nucleons $\left(\Psi \right)$, pions $\left(\pi \right)$ and a scalar $\left(\sigma \right)$ field. Short-range dynamics is parametrized in terms of heavy isoscalar and isovector vector fields $\left(v_{\mu }\right)$. The expectation value $\langle \sigma \rangle$ of the scalar field acts as a chiral order parameter normalized in the vacuum to the pion decay constant, $f_{\pi }\simeq 92$ MeV. Two classes of fluctuations beyond MF were then systematically studied: first, vacuum fluctuations which introduce an additional term proportional to ${\sigma }^{4}ln(\sigma /f_{\pi })$ in the MF partition function; and secondly, fluctuations involving pion loops and nucleon particle-hole excitations. The vacuum fluctuations can be included in an extended mean-field (EMF) approximation [69]. Fluctuations involving pion and nucleon loops are computed using non-perturbative functional renormalization group (FRG) methods. The parameters of the ChNM model, in particular those related to short-distance dynamics, are fixed to reproduce empirical nuclear physics data [70,71].

Figure 10 demonstrates the important role of fluctuations beyond mean-field approximation for the chiral order parameter $\langle \sigma \rangle$. In symmetric nuclear matter the mean-field approximation of the ChNM model correctly reproduces the first-order liquid-gas phase transition at low density. But the chiral order parameter also displays a strong first-order chiral phase transition with a Maxwell-constructed phase coexistence region starting already below $2{\rho }_0$. For neutron matter, which has no liquid-gas phase transition, the MF approximation nevertheless predicts a first-order chiral phase transition at densities around $\rho \sim 3{\rho }_0$, well within the range of densities realised in neutron stars. However, in both nuclear and neutron matter, the inclusion of fermionic vacuum fluctuations (i.e. the effect of the ground state zero-point energy) in the extended mean-field (EMF) approximation converts the first-order chiral phase transition into a smooth crossover and shifts it to densities $\rho >5{\rho }_0$. This effect is further enhanced by the additional fluctuations included in the full FRG calculation as demonstrated in Figure 10. As a result, chiral symmetry restoration is relegated to very high baryon densities beyond the inferred core densities in even the heaviest neutron stars (see Eq. (12)). A comparable impact of fluctuations on the phase structure is seen in alternative chiral models [72,73].

Another approach based on chiral symmetry is the parity-doublet model. In this model the active coupled baryonic degrees of freedom are the nucleon with spin-parity $1/2^{+}$ and its chiral partner with spin-parity $1/2^{-}$, where the latter is identified with the $N^{*}\left(1535\right)$ resonance. Spontaneous chiral symmetry breaking in vacuum manifests itself in the mass splitting of these two states, while the in-medium restoration pattern of this symmetry is signalled by the N and N* masses becoming degenerate. A recent detailed analysis [74] of the chiral order parameter in this model using extended mean-field approximation finds a chiral phase transition in nuclear matter, but at extremely high densities $(\rho >10{\rho }_0)$, far beyond the density scales reached in neutron stars.

4.4. Dense baryonic matter: a Fermi liquid picture

The constraints on the equation of state inferred from the current empirical data still permit an interpretation of neutron star core matter in terms of baryonic degrees of freedom, such as a system dominated by neutrons [36] with small fractions of protons and perhaps hyperons [75,76,77]. The inferred baryon chemical potential $\mu \left(\rho \right)=\frac{\partial \epsilon \left(\rho \right)}{\partial \rho }$ shown in Figure 3 does not distinguish between different species of baryons. Its behaviour displays nonetheless increasingly strong correlations at high densities. It is instructive to analyze the gross properties of this state of matter using the Landau theory of relativistic Fermi liquids [78]. Here we perform a schematic study assuming "neutron-like" quasiparticles [79], with baryon number $B=1$ and a density-dependent mass $m\left(\rho \right)$, ignoring other small admixtures in the composition of the dense medium. These quasiparticles are characterized by their (relativistic) Landau effective mass $m_L^{*}$ at Fermi momentum $p_F={\left(3{\pi }^2\rho \right)}^{1/3}$,

$$\begin{array}{c}\hfill m_L^{*}\left(\rho \right)=\sqrt{p_F^2+m^2\left(\rho \right)},\end{array}$$

(22)

together with an effective potential, $U\left(\rho \right)$, so that the baryon chemical potential can be written [79]:

$$\begin{array}{c}\hfill \mu \left(\rho \right)=m_L^{*}\left(\rho \right)+U\left(\rho \right).\end{array}$$

(23)

The median of the $\mu \left(\rho \right)$ posterior credible bands in Figure 3 is now taken as a guiding starting point to extract the baryonic quasiparticle properties. With an educated ansatz for $m\left(\rho \right)$, the density dependence of the potential $U\left(\rho \right)$ can then be deduced and further discussed. One possible choice is to take $m\left(\rho \right)$ (with $m\left(0\right)=939$ MeV) from the non-perturbative FRG calculation employing the chiral nucleon-meson field-theoretical model [70] discussed in Section 4.1. The resulting Landau effective mass $m_L^{*}\left(\rho \right)$ is shown in Figure 11 together with the potential $U\left(\rho \right)=\mu \left(\rho \right)-m_L^{*}\left(\rho \right)$. The 95% credible band of $\mu \left(\rho \right)$ in Figure 3 leads to an uncertainty of about 15% for U at high densities. It is instructive to fit the resulting quasiparticle potential by a series in powers of baryon density for $\rho \lesssim 5{\rho }_0$ (with ${\rho }_0=0.16{\mathrm{fm}}^{-3}$ as before):

$$\begin{array}{c}\hfill U\left(\rho \right)=\sum _n{u}_n{\left(\frac{\rho }{{\rho }_0}\right)}^n.\end{array}$$

(24)

The coefficients fitted to the median of $\mu \left(\rho \right)$ are:

$$u_1=90.9\mathrm{MeV},u_2=15.3\mathrm{MeV},u_3=3.2\mathrm{MeV},u_4=-0.4\mathrm{MeV}.$$

(25)

This pattern reflects a hierarchy of many-body correlations, recalling that the term linear in density represents two-body interactions, the term of order ${\rho }^2$ corresponds to short-range three-body forces, and so forth. The role of the repulsive N-body terms with $N>2$ is quite significant: at $\rho \simeq 4{\rho }_0$, corresponding to an average distance of about 1 fm between the baryonic quasiparticles, these terms contribute as much as the two-body forces to the potential U and generate the strong pressure to support heavy neutron stars. Of course these statements rely on the ansatz for the density dependent mass $m\left(\rho \right)$ which is guided by FRG calculations of pure neutron matter. A small fraction of protons in beta-equilibrated matter will not change this picture substantially. However, neutron star core compositions that deviate qualitatively from this simplified picture may lead to different conclusions.

Finally, consider the dimensionless Fermi liquid parameters, $F_0$ and $F_1$, of the spin-averaged quasiparticle interaction. In terms of the quasiparticle mass and potential, they are given as [79]:

$$F_0\left(\rho \right)=\frac{p_F}{{\pi }^2}\left[m\left(\rho \right)\frac{\partial m}{\partial \rho }+m_L^{*}\left(\rho \right)\frac{\partial U}{\partial \rho }\right],F_1\left(\rho \right)=-\frac{3U\left(\rho \right)}{\mu \left(\rho \right)}.$$

(26)

Further useful relations are: $1+F_0=N\left(0\right)(\partial \mu /\partial \rho )$ with the density of quasiparticle states at the Fermi surface, $N\left(0\right)=m_L^{*}{p}_F/{\pi }^2$, and $1+F_1/3=1-U/\mu =m_L^{*}/\mu$.

The result in Figure 12 shows a strongly increasing $F_0$ at high densities. This reflects once again the growing importance of many-body correlations as matter gets more and more compact. Such repulsive correlations are responsible for the increase of the sound velocity beyond its canonical conformal limit, $c_s>1/\sqrt{3}$, as seen in Figure 1. The Fermi liquid parameter $F_1$ is smaller in magnitude and has a negative slope indicating the decreasing effective mass at higher densities.

The Landau parameter $F_0$ displays the typical behavior of a strongly correlated Fermi system. It is interesting though to observe that in comparison with the leading Landau parameters in liquid ³He, the correlations in neutron star core matter are not extraordinarily strong. Values of $F_0\approx 9.3$ are reported for ³He at zero pressure, and $F_0\approx 68.2$ at a pressure of 27 bar [80]. Accordingly, the Fermi liquid parameters for matter in the density range, $\rho \lesssim 5{\rho }_0$, inside even the heaviest neutron stars, are significantly smaller. In this perspective the dense baryonic medium encountered in the center of neutron stars is perhaps not as extreme as often imagined.

5. Concluding remarks

The neutron star data base, recently extended by adding the heaviest so far observed pulsar (PSR J0952-0607), requires a much stiffer neutron star matter equation of state compared to previous analyses. As a consequence, baryon densities reached in the cores of even very heavy neutron stars do not exceed about five times the density of normal nuclear matter. In chiral model descriptions of the nucleon which suggest a scale separation between a valence quark core and a mesonic cloud, the compact cores with radii of about 1/2 fm would begin to touch and overlap only in the inner cores of extremely heavy neutron stars. The appearance of a strong first-order phase transition becomes unlikely under such conditions. This is consistent with the empirical results which suggest that extended phase coexistence regions and small minimum sound speeds are not favoured by the current data. This interpretation is further endorsed by the fact that fluctuations beyond mean-field approximation convert a possible first-order chiral phase transition into a smooth crossover at high densities. At the same time, a hadrons-to-quarks crossover scenario, or alternatively, baryonic matter as a strongly correlated relativistic Fermi liquid, remain both possible within the data-driven inferred constraints.