1 Introduction↩︎

Neutron stars (NSs) are among the most compact astrophysical objects known and, to date, are the only objects that provide empirical access to the behavior of cold, extremely dense matter. Observational advances may make it possible to address the composition of matter at the densities found in the cores of these stars [1]–[3]. In order to do so, the observations must be analyzed in the context of microphysical models. The microphysical assumptions underlying these models affect the predicted equation of state (EoS) and, ultimately, the macroscopic and measurable properties of NSs thereby allowing microphysical questions to be addressed through observations of astronomical objects.

As the EoS cannot currently be computed directly from the fundamental theory of quantum chromodynamics (QCD) at NS densities [4], [5], a large number of microphysical models have been developed to facilitate such studies, some of which are available through publicly accessible databases [6]–[11]. These phenomenological models encode expectations based on specific physical assumptions rather than first-principles calculations with well-defined uncertainties, and are therefore typically valid only within a limited density range where they provide a meaningful description of matter. Ideally, this range should include all densities realised in NSs (typically a few times saturation density \(n_\mathrm{sat}\approx 0.16~\mathrm{fm}^{-3}\)).

As we show here, microphysical models also constrain the EoS beyond the density domains in which they can be directly applied. This somewhat unintuitive feature arises from the requirement that the EoS at NS densities be consistent with first-principles perturbative QCD (pQCD) calculations at asymptotically high densities.

For most physical quantities, pQCD does not provide robust information at intermediate or lower densities. The EoS, however, is a special case: mechanical stability, causality, and thermodynamic consistency restrict how it can interpolate between the densities described by the models and the high-density regime constrained by pQCD [12]. As a result, these requirements impose constraints even in regions where neither the phenomenological models nor pQCD calculations can be directly applied.

We explore the behaviour that different NS-matter models predict at densities above those reached in NSs by employing a recent non-parametric, model-agnostic prior that samples the space of EoSs connecting the model predictions at NS densities with the pQCD EoS at high densities [13]. By constructing a large number of interpolating EoSs, we investigate the implications of the underlying model assumptions at higher densities. In particular, we find that many models, when extended to the maximal densities realised in NSs, imply a relatively specific continuation of the EoS at higher densities. This continuation typically involves a strong phase-transition-like change in the behaviour of matter just above the highest densities reached in NSs. We argue that this information can be used to further constrain EoS models.

The paper is organized as follows. In Section 2.1, we discuss our method for sampling extensions of microphysical models up to the regime where pQCD calculations are converged. In Section 2.2, we discuss our selection of representative microphysical models to use in this work. In Section 3, we examine the EoS extensions and discuss what different microphysical models imply for thermodynamic behaviour on the unstable branch.

2 Methods↩︎

The NS-matter EoS, \(p(\mu)\) with \(p\) the pressure and \(\mu\) the baryon chemical potential, can be computed in pQCD at high (baryon number) densities of around \(n_\mathrm{pQCD}\sim 20-40\, n_\mathrm{sat}\) [14]–[16]. While no NS reaches these densities, the requirement to reach these results at high densities constrains how the EoS can behave at lower densities.

In particular, given a NS-matter model extending up to some termination density, \(n_\mathrm{term}\), the requirement to reach the pQCD constraint in a causal, stable and thermodynamically consistent manner limits how the EoS can behave at intermediate densities \({n_\mathrm{term}< n < n_\mathrm{pQCD}}\). The origin of these pQCD constraints was demonstrated in [12]. By sampling possible interpolations between the endpoint of the NS-matter model and the pQCD limit, we can reveal what kind of behaviour is required to connect these two regimes.

To characterize how extreme the connection between the two limits is, a useful quantity is the pQCD tension index, defined in [17] as \[\label{eq:iqcd} \mathcal{I}_\mathrm{pQCD}\equiv \frac{p_\mathrm{pQCD}-p_\mathrm{term}-\Delta p_\text{min}}{\Delta p_\text{max}-\Delta p_\text{min}}.\tag{1}\] \(\mathcal{I}_\mathrm{pQCD}\) quantifies how close the bounds for all valid EoSs connecting the low- and high-density points are to the minimal (\(\mathcal{I}_\mathrm{pQCD}= 0\), with \(\Delta p_\text{min}\)) or maximal (\(\mathcal{I}_\mathrm{pQCD}= 1\), with \(\Delta p_\text{max}\)) pressure-difference constructions between \(\mu_\mathrm{term}\) and \(\mu_\mathrm{pQCD}\). These minimal and maximal constructions combine maximally causal segments with large first-order phase transitions. Explicit formulas for \(\Delta p_\text{min}\) and \(\Delta p_\text{max}\) are provided in [12].

2.1 Sampling of the EoS extensions↩︎

For any allowed model with \(\mathcal{I}_\mathrm{pQCD}\in (0,1)\), one can define a class of EoS extensions beyond the termination density.

Here, we construct such extensions using a Gaussian-process bridge (GPB) introduced in [13]. The method first samples the allowed functional space of EoS extensions in the interval \({n_\mathrm{term}< n < n_\mathrm{pQCD}}\) via a hierarchical self-similar refinement procedure, successively adding EoS points that remain consistent with all previously sampled points. The resulting EoSs span the allowed functional space between the low- and high-density limits and contain structures on all scales. The EoSs are then processed by diffusive filtering to impose a chosen correlation length for \(\mu(n)\). By construction, this procedure efficiently samples the full space of EoS extensions of a given low-density NS EoS that are thermodynamically consistent with the high-density pQCD constraints.

The uncertainty of the pQCD results are quantified by a choice of dimensionless renormalization scale \({X \equiv 3\bar\Lambda/ (2\mu)}\) with \(\bar\Lambda\) the renormalization scale in the modified minimal subtraction scheme. The value of \(X\) is usually chosen to minimize large logarithms in the perturbative results. We use for the pQCD EoS the scale-averaged result of [18], taking a log-uniform distribution of \(X\) values in the range \(X \in [1/2 , 2]\) [15]. We note that the apparent convergence of the pQCD result is slowest for the smallest values of \(X \approx 1/2\), rendering the small values less trustworthy. Extended discussions of the convergence of the pQCD results can be found in [15], [19], [20].

In practice, we generate the self-similar EoSs extensions between the chosen termination density of the model \(n_\mathrm{term}\), and the high-density limit at \({n_\mathrm{pQCD}= 30\, n_\mathrm{sat}}\), where we match to the pQCD results. As prescribed in [13], we diffuse the self-similar extensions between the termination density and \(40\,n_\mathrm{sat}> n_\mathrm{pQCD}\). Extending the upper diffusion region into the pQCD EoS range (from 30 to \(40\,n_\mathrm{sat}\)) enables a smooth variation of the sound speed when transitioning to the high-density regime. However, unlike in [13], we do not extend the diffusion into the model region at lower densities, allowing for sharp features in the EoS when transitioning from the model to the extensions. We consider this approach more conservative, as discussed later. For the final diffused extensions, we employ a hierarchical model of different logarithmically constant correlation lengths \(\sigma/n\), selecting them from within a uniform range \(\sigma/n \in [0.2, 0.4]\).

2.2 Model selection and representative models↩︎

We choose representative EoSs from the three classes of nucleonic models for the cold NS EoSs listed in the CompOSE database [6], namely microscopic calculations, non-relativistic density functional models, and relativistic density functional models. We focus on EoSs based on nuclear models fitted directly to experimental observables, i.e., nucleon-nucleon scattering data and properties of nuclei, binding energies, radii, and surface thickness. This excludes nuclear model parameterizations based on nuclear matter properties. As additional criteria, we check the EoSs for consistency with the well established EoS up to saturation density and that the EoS gives a maximum NS mass of at least \(2M_\odot\). In some cases, the parameter fits include the EoS of pure neutron matter. For each model, the corresponding reference and CompOSE entry are listed in 1.

For microscopic calculations we choose the EoSs denoted as APR [27] and BL(chiral) [26]. The APR EoS is constructed from the Argonne nucleon-nucleon potential fitted to nucleon-nucleon scattering data and three-body forces plus further corrections for the calculation of the EoS. As a representative for an EoS built from chiral nuclear forces, again constrained by two-body and three-body nuclear data, we take the EOS BL(chiral).

The non-relativistic density functionals, or Skyrme-type models, adopted below are the parameter set BSk24 [22] and Sly4 [21], [33]. The EoS BSk24 derives from Hartree–Fock–Bogoliubov nuclear mass model fitted to the binding energy of nuclei over nearly the entire chart of nuclides and the EoS of pure neutron matter. The Skyrme model SLy4 has been fitted to properties of selected nuclei and the EoS of pure neutron matter.

Similar to the non-relativistic density functionals, our chosen relativistic density functional models DD2 [23], [34] and FSU2R [25] have been fitted to properties of nuclei. In addition, the parameters of the set FSU2R are chosen to give a reasonable description of pure neutron matter by fixing the slope parameter \(L\) and have been modified from the original FSU2 model [35] to arrive at a smaller radius for a \(1.4M_\odot\) NS, in accord with the constraint on the tidal deformability from NS merger gravitational wave event GW170817. The nuclear model DD2 includes density-dependent coupling constants which are tuned to describe the nonlinear density dependence of the nucleon self-energy extracted from relativistic Brueckner–Hartree–Fock calculations. By this additional input, the low density pure neutron matter EoS up to saturation density can be described.

In addition to these nuclear models we investigate the effect of possible exotic phases and modern, nonconventional approaches for the cold NS EoS. In particular, for the EoS with exotic matter we take QHC18 [28] as a representative of an EoS with quark matter in the core, DDHdeltaY4 [29] as an EoS with hyperonic matter appearing at high density, and CMF-7 [30], [36], [37] as a chiral mean-field model which includes hyperons and Delta-baryons as quasi-particle degrees of freedom. A quarkyonic model was also implemented using the publicly available code from [38], based on the model described in [32], with the following parameters: symmetry-energy slope \(L = 50\,\mathrm{MeV}\), shell parameter \(\Lambda = 1400\,\mathrm{MeV}\), and quark drip density \(n_t = 0.3\,\mathrm{fm}^{-3}\). Finally, we take an EoS using a holographic approach to describe the NS-matter EoS, namely the V-QCD EoS [31].

All EoSs used in this study correspond to zero-temperature EoSs (or the lowest-temperature slice available) in \(\beta-\)equilibrium. For general-purpose EoSs, the \(\beta-\)equilibrium composition was recalculated according to the CompOSE prescription.

3 Results↩︎

1 shows the possible extensions of the fully hadronic EoSs, assuming that the models remain valid up to the densities reached in a canonical \(1.4M_\odot\) NS. The EoSs are displayed in terms of the squared speed of sound, \(c_s^2\), and the (normalized) trace anomaly, \({\Delta \equiv (\epsilon - 3p)/(3\epsilon)}\), as functions of density. Both quantities characterize the properties of matter. The trace anomaly, \(\Delta\), has the advantage of being obtained by integrating \(n(\mu)\) and is therefore insensitive to local variations in \(n(\mu)\). By contrast, there is no thermodynamic constraint on how rapidly \(c_s^2(n)\) may vary.

For each model, the black solid line extends up to the central density of a \(1.4M_\odot\) star. The continuation of each model to higher densities is shown by the red dashed line, terminating at a black point that represents the core of the maximally massive star. At high densities, \({n = 40\, n_\mathrm{sat}}\), the EoSs are given by the pQCD calculation, corresponding to \(c_s^2 \approx 1/3\) and a small \(\Delta \approx 0\). The lines are colored based on the value of the dimensionless pQCD renormalization scale \(X\) with three ranges equally spaced in \(\log(X)\), i.e., \({X_{\rm low} \in [0.5, 0.79]}\) (light purple), \({X_{\rm mid} \in [0.79, 1.26]}\) (magenta), and \({X_{\rm high} \in [1.26,2]}\) (dark violet). We note that the \(X_{\rm mid}\) range approximately corresponds to the ‘most-consistent’ prediction for the next-to-next-to-next-to-leading order pQCD EoS as presented in [16] and may be thought of as a proxy for the potential full pQCD result at this order.

The range between \(n(1.4 M_\odot)\) and \(40\,n_\mathrm{sat}\) is interpolated using a number of draws from the GPB described above. The transition to the interpolating GPB is—or at least can be—smooth, and a broad range of different possible intermediate behaviors are sampled. We see that there are consistent interpolations for all values of \(X\). Similarly, the resulting mass-radius curves show broad range of behaviors.

Extending these purely nucleonic models to densities reached in maximally massive stars (the TOV density), the situation is very different (see 2). The Sly4, PCP(BSK24), DD2, SFHO, and BL models are incompatible with the higher \(X_{\rm high}\)-range. For these models, the sampled interpolations for the less restrictive \(X_{\rm mid}\) and \(X_{\rm low}\) ranges are significantly constrained.

In the case of APR, the model EoS is terminated at a density below the maximum density reached in the TOV solution, specifically where the EoS first becomes acausal. This density corresponds to the central density of a NS with mass \(M\approx2.07M_\odot\). Since the termination occurs at a lower density, the resulting extrapolations span a broader range.

In all of these cases (besides FSU2R) the behavior of the matter must have an abrupt change immediately the TOV density; the matter must undergo a dramatic softening resembling a strong first-order phase transition at densities immediately above those reached in NSs, which extends for a density interval \(\Delta n \approx 15\, n_\mathrm{sat}\). This is most clearly seen as the significant change of the slope of the trace anomaly. Note that since this transition takes place above the TOV density, it does not cause the star to collapse; rather it is coincidental that the collapse and phase-transition-like behavior coincide, which seems highly unlikely. The change in thermodynamic behaviour is also reflected in the mass-radius plot, leading to a kink feature entering the unstable branch.

We note that diffusing into the model region would lead to even more restrictive behaviour, as it would require a smooth connection to the stiff EoS, extending the high-\(c_s^2\) behaviour of the model over a larger density range. In most cases, this would leave no viable extensions at all. Therefore, as mentioned before, allowing for abrupt behaviour at the TOV point is a more conservative way to treat the low-density limit.

In contrast, FSU2R permits extensions for the full range of \(X\). This is a consequence of its relatively soft behaviour and positive trace anomaly all the way up to the TOV density. Unlike the other hadronic models considered here, however, the parameters of FSU2R were already informed by astrophysical observations. It is important to note, that while this behaviour is compatible with pQCD constraints for larger values of \(X\), the model predicts large radii above 12 km, which is in tension with the NICER radius measurements of PSR J0614\(-\)​3329 [39] and PSR J0437\(-\)​4715 [40].

The representative set of non-hadronic EoS models is shown in 3, where each model is taken up to the TOV density. As evident from the figure, the behaviour required by the models at higher densities is generally less restricted than for purely hadronic models. However, BFH(QHC19-B), OPGR(DDHdeltaY4), and DS(CMF) remain sufficiently stiff that they require an abrupt change in the speed of sound at the TOV density, leading to effectively first-order-phase-transition-like extensions for the \(X_{\rm high}\) range. In contrast, for the V-QCD and Quarkyonic models, the softening already occurs before the TOV point. As a result, the required behaviour above the TOV density is considerably less restrictive and does not force first-order-phase-transition-like behaviour.

4 Conclusions↩︎

In this study, we reveal what NS-matter models imply for the EoS behaviour beyond different termination densities by exploring the allowed range of possible valid extensions from the termination point to the perturbative-QCD limit. We find that for purely hadronic models, the behaviour remains largely unconstrained if the model is only trusted up to the central density of a \(1.4M_\odot\) NS. However, if one assumes that hadronic models remain valid up to the maximal density reached in stable NSs, the situation changes dramatically. In this case, the EoS becomes incompatible with higher values of the pQCD renormalization-scale parameter \(X\) and is forced into highly specific behaviour for allowed values of \(X\).

This behaviour is characterized by an abrupt change at the TOV point, with an effectively first-order-phase-transition-like density jump of \(\Delta n \approx 15\,n_\mathrm{sat}\). Importantly, this behaviour does not itself destabilize the star; rather, it emerges as a consequence of enforcing a causal, stable, and thermodynamically consistent connection to the pQCD limit. Such a coincidental behaviour across all viable extensions seems highly unlikely. In contrast, this behaviour is not enforced if the EoS softens before the TOV density, as demonstrated by models exhibiting crossover or first-order phase transitions, such as V-QCD.

While we investigated this behaviour only for a representative set of EoSs, the extremeness of the required extension to pQCD can be quantified for all zero-temperature, \(\beta\)-equilibrium EoSs available in the CompOSE database. This is shown in 4, where the measure of extremeness is the pQCD tension index, defined in 1 , shown as a function of central density. Values of \(\mathcal{I}_\mathrm{pQCD}>1\) indicate that the model is inconsistent with pQCD for a given value of \(X = 1\). Values of \(\mathcal{I}_\mathrm{pQCD}=1\), or slightly below, correspond to the highly specific behaviour associated with a forced first-order-like transition.

The purple crosses in 4 indicate the central density of a \(1.4M_\odot\) NS, where all models predict \({\mathcal{I}_\mathrm{pQCD}\approx0.5}\), corresponding to maximal versatility in the allowed behaviour. The green dots correspond to the TOV point for EoSs that include tabulated quark or hyperonic degrees of freedom, including the EoSs shown in 3, which are additionally highlighted in the figure.

The purple dots correspond to the TOV point of the purely hadronic EoSs. Their clustering close to the \(\mathcal{I}_\mathrm{pQCD}=1\) line suggests that the set of EoSs chosen for 2 is representative of the broader trend: purely hadronic EoSs trusted up to the TOV point generically lead to large first-order-phase-transition-like behaviour above it. The exceptions to this trend, such as the highlighted FSU2R, are typically softer EoSs at intermediate densities.

The EoSs showcased in this work have been widely used in binary NS-merger or core-collapse supernova simulations, as well as many other studies. Our results therefore imply that such simulations may carry an implicit assumption about the behaviour of QCD matter at higher densities: namely, that it undergoes a strong first-order-phase-transition-like change that happens to coincide with the TOV point. Since such behaviour appears highly artificial, avoiding it while remaining consistent with astrophysical constraints requires a softening of the EoS before the TOV point, for example through a crossover to quark matter [3], [41] or a first-order phase transition [42], [43] at lower densities.

In conclusion, we have demonstrated that pQCD constraints provide useful guidance for NS-EoS model building. Our results indicate that models with additional degrees of freedom, leading to a softening of matter below \(M_{\rm TOV}\), are preferred. Consequently, our work disfavors purely nucleonic matter and supports the appearance of additional degrees of freedom—without indicating a preference for any particular type—in stable massive NS.

Data Availability↩︎

An interactive web application based on the framework developed in this work is available at [44]. A Jupyter notebook implementation is also available on Zenodo [45]. The application generates a prior ensemble of allowed EoS extensions connecting a low-density neutron-star matter model to the high-density pQCD regime. Given a user-specified termination point, \((\mu_\mathrm{term}, n_\mathrm{term}, p_\mathrm{term})\), it constructs and visualizes the ensemble of allowed extensions.

Author Contributions↩︎

Authors are listed in alphabetical order.

References↩︎