NPLQCD collaboration

The determination of the internal structure of neutron stars presents a long-standing and important challenge for nuclear theory. Since neutron stars were first predicted in the 1930s and observed 30 years later, many models for the structure of their interiors have been proposed, including various phases of nuclear matter, mesonic condensates and hyperonic matter, and deconfined quark cores Ozel:2016oaf ; Oertel:2016bki ; Baym:2017whm ; Mannarelli:2019hgn ; Burgio:2021vgk ; Lattimer:2021emm . As observational data and terrestrial probes of the relevant nuclear densities are not sufficiently constraining, most of these possibilities for the neutron star equation of state (EoS) remain viable. From a theoretical perspective, it is expected that neutron star interiors can be described by the Standard Model of particle physics, however in a regime where the strong interactions are non-perturbative. The numerical technique of lattice quantum chromodynamics (LQCD) is applicable at such large couplings, but is beset by a notorious sign problem at nonzero baryon chemical potential deForcrand:2009zkb , prohibiting its direct application. Consequently, theoretical approaches to the nuclear EoS are based on models and interpolations between phenomenological constraints from nuclear structure and perturbative QCD (pQCD) calculations that are valid at asymptotically large chemical potentials (see for example, Refs. Annala:2023cwx ; Semposki:2024vnp ; Han:2022rug ; Jiang:2022tps ; Raaijmakers:2019dks ; Huth:2020ozf ; vanDalen:2005sk ; Typel:2005ba ; Kolomeitsev:2004ff ; Liu:2001iz ; Leonhardt:2019fua ). In light of this, any rigorous information that can impact such analyses is of paramount importance.

In this work, the first non-perturbative QCD constraint on the nuclear EoS with complete quantification of systematic uncertainties is presented. These calculations build upon the proof-of-principle, single lattice spacing results of Ref. Abbott:2023coj with improved methodology, increased statistical precision, an extrapolation to the continuum limit, and an interpolation to the physical quark masses, enabling a systematically-controlled result to be achieved for the first time. The pressure and other thermodynamic properties of low-temperature isospin-dense matter are determined over a wide range of densities and chemical potentials, spanning all scales from hadronic to perturbatively-coupled. At small values of the isospin chemical potential, $\mu_{I}$, the results agree with chiral perturbation theory ($\chi$PT) son_qcd_2001 ; Kogut:2001id at next-to-leading order (NLO) Adhikari:2019zaj ; Andersen:2023ofv . At large $\mu_{I}$, the results are seen to agree with pQCD with pairing contributions Fujimoto:2023mvc . The comparison of next-to-next-to-leading-order (NNLO) pQCD predictions Freedman:1976xs ; Freedman:1976dm ; Freedman:1976ub ; Baluni:1977ms ; Kurkela:2009gj ; Fujimoto:2023mvc for the pressure (partial next-to-next-to-next-to-leading-order results are also available Gorda:2021znl ) with the continuum limit of the LQCD calculations provides a determination of the superconducting gap as a function of $\mu_{I}$. This is seen to agree with the leading-order perturbative calculation of the pairing gap Fujimoto:2023mvc , but is more precise. The speed of sound in isospin-dense matter is also seen to significantly exceed the conformal limit of $c_{s}^{2}/c^{2}\leq 1/3$ over a wide range of $\mu_{I}$. A Bayesian model mixing approach that combines $\chi$PT, LQCD and pQCD provides a determination of the zero-temperature EoS for isospin-dense QCD matter valid at all values of the isospin chemical potential.

These results provide constraints for phenomenological frameworks seeking to describe the QCD phase diagram MUSES:2023hyz , and from simple path-integral relations Cohen:2003ut ; Cohen:2004qp ; Fujimoto:2023unl ; Moore:2023glb , the determination of the pressure in isospin-dense matter provides a non-perturbative, model-independent bound on the pressure of isospin-symmetric QCD matter at nonzero baryon chemical potential and hence on the nuclear EoS. The current results therefore provide a systematically-controlled QCD bound at all densities, and the impact on nuclear phenomenology is briefly discussed.

Thermodynamic relations: Thermodynamic quantities are accessed in this work by building an approximation to the grand canonical partition function, valid at low-temperature. The grand canonical partition function is defined at a temperature, $T=1/\beta$, and isospin chemical potential, $\mu_{I}$, by

where the sum is over all states, $s$, and $E_{s}$ and $I_{z}(s)$ correspond to the energy and $z$-component of isospin of a given state, respectively. Since states of different $I_{z}$ but the same $I$ are approximately degenerate, contributions from states with $I_{z}<I$ are suppressed by ${\cal O}(e^{-\beta\mu_{I}})$ relative to those with $I_{z}=I$ and can therefore be neglected. Additionally, only the ground state for each isospin contributes at low temperature. The summation can therefore be approximated in terms of these $I_{z}=I$ ground states, which can be labeled by their isospin charge $n=I=I_{z}$, and truncated at some $n_{\rm max}$ giving

$E_{0}=0$ is chosen, and this approximation is valid for values of $\mu_{I}$ such that the truncation at $n_{\rm max}$ does not affect the result significantly.

For an observable $\mathcal{O}$ that only depends on the energy and isospin charge of the system, the thermodynamic expectation value of $\mathcal{O}$ can be computed as

The energy density can be computed using the quantity $E_{n}/V$, while the isospin-charge density can be computed from the expectation value of $n/V$, where $V$ is the volume of the system. Derivatives of observables can also be computed using

which results from directly differentiating Eq. (3). This leads to the following expressions for the pressure,

and speed of sound defined by the isentropic derivative of the pressure with respect to the energy density, $\epsilon$,

Previous work Detmold:2008fn ; Detmold:2010au ; Detmold:2011kw ; Detmold:2012pi2 ; Abbott:2023coj studied isospin-dense matter through a canonical partition function approach by using the thermodynamic relation

to determine the isospin chemical potential from the extracted energies. Other studies have added the isospin chemical potential directly to the QCD action Kogut:2002tm ; Kogut:2002zg ; Brandt:2017oyy ; Bornyakov:2021mfj ; Brandt:2022hwy . These studies probe isospin dense QCD at $\mu_{I}\lesssim 2$ and focus primarily on nonzero temperature. The results in this work are consistent with the low-temperature results in those works, but span a larger range of chemical potentials.

The primary advantage of the method used here in comparison to the method of Ref. Abbott:2023coj is that $\mu_{I}$ enters as an input to the calculation of thermodynamic quantities rather than being derived from the isospin charge of the LQCD data and therefore is not subject to statistical and systematic uncertainties. A more detailed comparison to the approach of Ref. Abbott:2023coj is discussed in the Supplementary Material

Color-superconducting gap: At large isospin chemical potential, asymptotic freedom guarantees the validity of pQCD and the resulting prediction of a color-singlet superconducting state at zero temperature son_qcd_2001 . In this state, Cooper-pairs of quark–anti-quark fields condense, leading to a superconducting gap with order parameter $\langle\overline{d}_{a}\gamma_{5}u_{b}\rangle=\delta_{ab}\Delta$, where $a$ and $b$ are color indices. The gap $\Delta$ can be computed perturbatively son_qcd_2001 ; Fujimoto:2023mvc , with the next-to-leading order result given by

where $\tilde{b}=512\pi^{4}g^{-5}$ and $g=\sqrt{4\pi\alpha_{s}(\mu_{I})}$ is the strong coupling at the scale $\mu_{I}$. Notably, the prefactor of $1/g$ in the exponent of Eq. (8) is smaller than the analogous coefficient in the baryon-density case by a factor of $1/\sqrt{2}$ son_qcd_2001 , leading to an exponential enhancement of the gap and its effects in isospin-dense QCD. If pQCD is reliable for a given $\mu_{I}$ and $\mu_{B}=3\mu_{I}/2$, then the isospin-dense gap bounds the baryonic gap, in which there is significant phenomenological interest Kurkela:2024xfh .

The nontrivial background in the presence of the gap induces a change in the pressure of the system. This change can be computed perturbatively, as has been done at NLO in Ref. Fujimoto:2023mvc with the result

This difference allows for an indirect extraction of the gap by comparing the lattice QCD pressure with the pressure derived in pQCD without pairing.

LQCD calculations: Following the methods and analysis techniques developed in Ref. Abbott:2023coj , the energies of systems of isospin charge $n\in\{1,\ldots,6144\}$ are determined from two-point correlation functions

calculated on four ensembles of gauge field configurations whose parameters are shown in Table 1. Here, $\pi^{-}({\bf x},t)=\pi^{+}({\bf x},t)^{\dagger}=-\overline{d}({\bf x},t)\gamma_{5}u({\bf x},t)$ is an interpolating operator built from $u$ and $d$ quark fields that creates states with the quantum numbers of the $\pi^{-}$. The correlation functions are computed using the symmetric polynomial method of Ref. Abbott:2023coj from sparsened Detmold:2019fbk quark propagators computed using a grid of 512 source locations on a single timeslice of each configuration.

The relevant ground-state energies, $E_{n}$, are determined from an analysis of the $t$-dependence of the effective energy functions

where $a$ is the lattice spacing, $\vartheta_{n}(t)$ and $\sigma_{n}(t)$ are the mean and standard deviation of $\log C_{n}(t)$, and the second equality is under the assumption of log-normality Abbott:2023coj . $N_{b}=2000$ bootstrap resamplings are used on each ensemble to assess the statistical uncertainties and address correlations. As in Ref. Abbott:2023coj , on each bootstrap the energy is given by the value of the effective mass on a randomly chosen time inside the effective mass plateau region.

Given the energies determined on each ensemble for systems of isospin charge $n\in\{1,\ldots,6144\}$, Eqs. (5) and (6) are used to determine the pressure, the energy density, and speed of sound. The action used in these calculations is perturbatively improved, so discretization effects are ${\cal O}(a^{2},g^{2}a)$. The mass dependence of quantities evaluated over the range of quark masses used in the calculations is expected to be described linearly in $m_{ud}\sim m_{\pi}^{2}$. Each quantity ${X\in\{P/P_{\rm SB},\epsilon/\epsilon_{\rm SB},c_{s}^{2}/c^{2}\}}$ (where $P_{\rm SB}=\epsilon_{SB}/3=\mu_{I}^{4}/32\pi^{2}$ is the pressure of a Stefan-Boltzmann gas) is fit with forms including arbitrary $\mu_{I}$ dependence and terms ${\cal O}(a^{2},a^{2}\mu_{I},a^{2}\mu_{I}^{2},(m_{\pi}^{2}-\overline{m}_{\pi}^{2}))$, where $\overline{m}_{\pi}=m_{\pi^{+}}=139.57039$ MeV ParticleDataGroup:2022pth , with coefficients independent of $\mu_{I}$. The systematic uncertainty from the extrapolation is assessed by combining fits with all possible subsets of the terms above through model averaging Rinaldi:2019thf ; NPLQCD:2020ozd ; Jay:2020jkz . The systematic uncertainty and the statistical uncertainty are combined under the bootstrap procedure. The lattice determinations of the pressure, energy density and speed of sound are shown in the Supplementary Material, along with further details.

The calculated pressure is shown in Fig. 1 for each lattice ensemble as well as for the continuum-limit, physical-quark-mass interpolation. The LQCD results for the different volumes, lattice spacings and quark masses used in the calculations are seen to agree with each other within uncertainties and also with the physical mass, continuum limit extraction. The LQCD pressure agrees with NLO $\chi$PT at small values of the chemical potential and can be compared with NNLO pQCD Freedman:1976dm ; Freedman:1976ub ; Freedman:1976xs ; Baluni:1977ms ; Kurkela:2009gj ; Fujimoto:2023mvc at large values of the chemical potential. The mild tension seen between LQCD and pQCD for $\mu_{I}\in[1500,2250]$ MeV potentially indicates the presence of a superconducting gap Fujimoto:2023mvc (it could alternatively be an indication of the breakdown of pQCD, although NLO and NNLO pQCD results are in agreement over this range).

The corresponding speed of sound is seen to exceed the conformal limit of $c_{s}^{2}/c^{2}=1/3$ over a wide range of the isospin chemical potential (see Fig. 3 below). While this behavior was seen in Refs. Brandt:2022hwy ; Brandt:2022fij ; Abbott:2023coj (and also in two-color QCD Iida:2022hyy ; Iida:2024irv ), the present results confirm that this is not a lattice artifact and that such behavior is possible in strongly interacting QCD matter. This suggests that the assumption that the speed of sound remains below this value in baryonic matter is questionable.

Given the LQCD calculation of the pressure, a determination of the superconducting gap can be made by subtracting the pQCD calculation of the pressure in the absence of the gap. In the range of chemical potentials where pQCD is a controlled expansion, this determines the gap, accurate to the same order as the perturbative subtraction. Figure 2 shows the extracted gap found using the NNLO pQCD pressure subtraction as well as a comparison to the pQCD gap in Eq. (8) evaluated at scales $\bar{\Lambda}=\mu_{I}\times\{0.5,1.0,2.0\}$ as a guide to uncertainty. As can be seen, the gap extracted from the LQCD calculations is in agreement with the perturbative gap for $\mu_{I}\in[1500,3250]$ MeV but is considerably more precisely determined than the uncertainty from perturbative scale variation. Since there is agreement with the perturbative estimate, the gap is also most likely larger than the corresponding gap for baryonic matter.

Equation of state for isospin-dense mater: The continuum-limit lattice QCD calculations presented above span isospin chemical potentials from just above the pion mass to values where pQCD appears to converge. Consequently, by combining the LQCD results with $\chi$PT and pQCD, the zero temperature EoS of isospin-dense matter can be described for all $\mu_{I}$ with uncertainties quantified using Bayesian inference. The functional dependence of each overlapping theoretical constraint on $\mu_{I}$ is modelled by a correlated Gaussian distribution. The ensemble of constraints is combined via a Gaussian Process (GP), following similar work for the nuclear EoS Semposki:2024vnp ; Melendez:2019izc ; Essick:2020flb ; Mroczek:2023zxo . Theoretical uncertainties of $\chi$PT are estimated from the difference between the NLO and LO results, and uncertainties in pQCD are assessed from scale variation over $\bar{\Lambda}\in\mu_{I}\times[0.5,2.0]$. Figure 3 shows the GP-model results for the speed of sound in comparison to the three theoretical inputs. A data file for evaluating this model accompanies this article. With a complete quantification of the isospin-dense equation of state, phenomenological implications such as the existence of pion stars Brandt:2018bwq ; Andersen:2022aig ; Stashko:2023gnn and the isospin effects that distinguishes pure neutron matter from symmetric nuclear matter STEINER2005325 ; LI2008113 can be further investigated.

Constraining the nuclear equation of state: The partition function of two-flavor QCD with an isospin chemical potential, $\mu_{I}$, can be written in terms of the path integral

where $A$ is the gluon field, $\mathcal{D}(\mu)\equiv D\!\!\!\!/+m-\mu_{\mathrm{q}}\gamma_{0}$ is the Dirac operator with quark chemical potential $\mu_{q}$, $S_{G}$ is the gauge action, and $\int_{\beta}[dA]$ indicates integration over gauge fields with period $\beta$ in the temporal direction. As first shown in Refs. Cohen:2003ut ; Cohen:2004qp , this partition function bounds the partition function of two-flavor QCD with equal chemical potentials for $u$ and $d$ quarks

as

By the monotonicity of the logarithm, the above inequality directly translates into an inequality between the pressures of the two media as a function of the energy density. Consequently, the isospin-dense EoS bounds the EoS for symmetric nuclear matter. At large values of the quark chemical potentials, where pQCD is valid, this bound becomes tight as differences between the partition functions enter only at ${\cal O}(\alpha_{s}^{k})$ for $k\geq 3$ Moore:2023glb . This bound was explored in Ref. Fujimoto:2023unl based on the previous lattice QCD results Abbott:2023coj at a single lattice-spacing and unphysical quark masses. Here, Fig. 4 presents updated bounds based on the continuum limit lattice QCD results at the physical quark masses, $\chi$PT, and perturbative QCD through the GP-model. While the bounds from isospin-dense matter do not significantly constrain phenomenological nuclear equations of state within the uncertainties that are typically presented Annala:2023cwx ; Semposki:2024vnp ; Han:2022rug ; Jiang:2022tps ; Raaijmakers:2019dks ; Huth:2020ozf ; vanDalen:2005sk ; Typel:2005ba ; Kolomeitsev:2004ff ; Liu:2001iz ; Leonhardt:2019fua , the bounds are independent of modeling uncertainties that enter the nuclear EoS in the regions that are unconstrained by nuclear structure or pQCD calculations. GP models without lattice QCD constraints result in significantly larger uncertainties in the position of the bound, in particular in the lower right corner of the red band of Fig. 4.

Summary: In this letter, a determination of the equation of state of isospin-dense matter for the complete range of isospin chemical potential at zero temperature is presented for the first time. To achieve this, continuum limit LQCD calculations are combined with pQCD calculations and $\chi$PT through a model-mixing approach in overlapping regions of isospin chemical potential. Comparison to pQCD enables a determination of the superconducting gap, and QCD inequalities translate the isospin-dense EoS into rigorous bounds on the nuclear EoS relevant for astrophysical environments.