nEoS: Neutron Star Equation of State from hadron physics alone

TL;DR

This paper provides a set of neutron star equations of state derived solely from hadron physics principles, enabling tests of gravity theories without relying on observational data or General Relativity.

Contribution

It offers publicly available EoS tables based on hadron physics, free from observational feedback, to test gravity theories with neutron star data.

Findings

EoS constrained by hadron physics and fundamental principles

Applicable for testing gravity theories independently of observations

Combines chiral potentials and pQCD results for a broad density range

Abstract

We contribute a publicly available set of tables and code to provide Equations of State (EoS) for matter at neutron star densities. Our EoSes are constrained only by input from hadron physics and fundamental principles, without feedback from neutron star observations, and so without relying on General Relativity. They can therefore be used to test General Relativity itself, as well as modified gravity theories, with neutron star observables, without logical circularity. We have adapted state of the art results from NN chiral potentials for the low--density limit, pQCD results for the asymptotically high-density EoS, and use monotonicity and causality as the only restrictions for intermediate densities, for the EoS sets to remain as model-independent as is feasible today.

Figures

Equations

$$\epsilon =n\left(M_N+\frac{E}{A}\right)$$

$$P=n^2\frac{d(E/A)}{dn}.$$

$$P_{\mathrm{QCD}}=\left[\frac{3}{4{\pi }^2}{\left(\frac{{\mu }_B}{3}\right)}^4\right]\left(0.9008-\frac{0.5034X^{-0.3553}}{{\mu }_B-1.452X^{-0.9101}}\right).$$

$$0\le \frac{dP}{d\epsilon }\le 1,$$

Full text

††: *J. Phys. G: Nucl. Phys.

*Focus Issue on Hadrons and Gravitational Waves after GW170817

Keywords: Neutron Matter, Equation of State, Tests of General Relativity, Neutron star EoS.

1 Introduction

Neutron stars have long offered possible tests of General Relativity at the energy density of nuclear matter [1, 2, 3], but the discovery of gravitational waves has accelerated the pace of activity [4, 5], and the multimessenger observation of neutron star inspirals [6] has opened an entire new field [7].

In relating neutron star measurements such as masses, radii, vibrational modes or gravitational wave emission to General Relativity (GR), a key component is the Equation of State (EoS) of the material in the neutron star.

A broadly held view is that hadron theory is hopelessly difficult and that the Equation of State of neutron stars will eventually be determined or at least strongly constrained from neutron star observations [7]. This would be akin to other instances, such as the determination of parton distribution functions, where theory computations, with strong model assumptions that thwart control of systematic uncertainties, are substituted by direct use of high–energy data–based extractions in predicting cross sections.

However, our point here [8] is that an EoS determined this way cannot be used to constrain theories beyond General Relativity because it already assumes GR (in its use of Tolman-Oppenheimer-Volkoff equations to calculate an observational mass–radius diagram, or in its use of the GR gravitational wave solutions, for example). To avoid circular reasoning, the community needs an effort to provide EoS purely from hadron theory. The purpose of this brief article is to report an endeavor in this direction that we have carried out and now put at the disposal of other researchers.

In section 2 we collect and analyze references to what, in our opinion, currently forms the most reliable source of knowledge on the EoS from first-principles hadron theory. This is in the first place perturbative Quantum Chromodynamics (pQCD) for asymptotically high density (presumably valid outside the range of densities where neutron stars are sensitive to the EoS, but constraining the EoS nonetheless) and Chiral Perturbation Theory (ChPT) for low and moderate densities where neutrons in medio behave not too differently from neutrons in vacuo.

In section 3 we briefly show how the basic principles of causality and monotonicity of the EoS help constrain it through the intermediate density regime most relevant for neutron star properties. Section 4 then describes our wrapper to put this information within reach of Numerical Relativity practitioners and others who may want to test General Relativity or modifications thereof with neutron stars (some tests have been run and reported in Section 5). Our conclusions and suggestions for further progress are then spelled out in section 6.

To summarize the outcome of the work: we provide a set of Equations of State that sample the uncertainty bands of Chiral Perturbation Theory (at low density), likewise sample the band of perturbative QCD (at high density) and interpolate between both at intermediate density with the only restrictions of causality and monotonicity.

2 State of the art computations

Historically, the $\beta$–equilibrated (mostly neutron) matter EoS for neutron stars has been obtained by feeding “realistic” nuclear potentials into the many body methods appropriate for continuous matter [9]. In the last decades, however, a new family of potentials has appeared, based on modern chiral interactions, that promise to improve the precision of low–energy QCD. (The LENPIC collaboration [10] has demonstrated that the spectroscopy of light nuclei can be successfully treated, so that a first–principles foundation for nuclear physics is ever closer). Such improvement is naturally applicable to the low–density EoS, and because the interactions are derived from QCD through effective theory, this makes it possible for the first time to obtain systematic uncertainties in the obtained EoS.

Among the chiral potentials, we have consulted several contemporary works [11, 12, 13, 14, 15] and extracted the EoS (pressure $P$ as function of the energy density $\varepsilon$) from their data via

[TABLE]

Here, $E/A$ is the energy per nucleon in nuclear matter, $n$ is the nucleon–number density, and $M_{N}$ the nucleon mass. Chiral perturbation theory is used only up to an energy density of order 200 MeV/fm3 ($0.2$ GeV4), so that the nucleons are not very relativistic (an exception is the work of [12] that extends to 375 Mev/fm3).

It is worth noting that chiral perturbation theory follows the systematic counting of nuclear–force terms, described for example in [16]. In these implementations, the precise counting in vacuo has been employed: assessments based on the counting modified for nuclear matter [17] are still needed [18] and it would be interesting to have systematic studies thereof. With this caveat in mind, we believe that the current status of the calculations already allows us to set reasonably systematic uncertainty bands to low–energy predictions.

The EoS in our collection use Chiral Perturbation Theory with baryons. We use ChPT calculations that go at least to NLO and reach in some cases up to N4LO. The higher order sets typically treat the two–body nucleon–nucleon interaction properly but use an approximation to the three-body force (or include it only at the lowest order). We acknowledge that this causes a systematic bias, lowering the repulsion of the nuclear interaction and thus the stiffness of the EoS (since it is known that three–body forces lower the attraction [19]). This is an aspect where systematic improvement is possible in the next few years. Further work by specialists should reaffirm our trust in the low-density uncertainty band, and address the systematic error introduced by simplifying three and many-body forces, which is making current EoS a tad too soft.

Further, the chiral interactions are cut off at a momentum transfer where the chiral expansion deteriorates and higher resonances or quark-level effects have to be taken into account. The EoS sets available in the literature use cutoffs of 450 MeV and either 500 or 600 MeV depending on the reference. In Ref. [12] the long–range part of the interaction is cut off in coordinate space instead of momentum space, and the sensitivity to the cutoff is exposed by providing two sets with $R=0.8$ m and 0.9 fm respectively (the short-range contact terms have a standard momentum cutoff $\Lambda=2/R$).

At high density, on the other hand, eventually one enters the regime of pQCD. This has been studied in [20, 21] and a new, partial N3LO calculation has been presented in [22] 111In fact, this group [23] has, in a study close to our philosophy, interpolated the EoS between their high density computation and a selected low-density one, but incorporating the constraints from neutron star observables. This makes their interpolated EoS less uncertain for conventional astrophysics applications, but not usable to constrain theories modifying General Relativity.. For the equation of state at baryon chemical potential above $\mu_{\rm match}\sim 2.6$ to 2.8 GeV we have made use of the earlier, completed, computations to provide constraints at high density on our sampling of the Equation of State.

The authors of [21] provide a handy parametrization of their field theoretical calculation. The pressure (in GeV4) is given in pQCD, as a function of the baryon chemical potential $\mu_{B}$ (in GeV) and a parameter $X$ proportional to the renormalization scale ($X\equiv 3\bar{\Lambda}/\mu_{B}$) by the simple expression

[TABLE]

The first factor between brackets corresponds to the free gas and the second factor stems from the $O(\alpha_{s}^{2})$ corrections. As recommended by those authors, we have taken as limits for the $X$ parameter the values 1 (low) and 4 (high) to define a band of allowed EoS in the high density regime. The renormalization scale for the two lines that delimit that band (see figure 2) is then, in the $\bar{\rm MS}$, given by $\bar{\Lambda}=\frac{\mu_{B}}{3}$ and $\bar{\Lambda}=\frac{4\mu_{B}}{3}$ respectively.

3 Interpolation at intermediate densities

and construction of the sample EoS in the entire range

The intermediate density range is the most important to determine key properties of neutron stars such as the maximum allowed mass (known to be about two solar masses [24], with the highest mass claim [25] currently being $2.27^{+0.17}_{-0.15}M_{\odot}$). Yet it is also the region where a consistent expansion of QCD is not workable: we can only resort to basic theoretical properties, causality ($c_{s}<c$) and monotony ($c_{s}\geq 0$).

As we do not wish to use polytropic equations of state (assuming such form introduces a systematic bias, and we want to have different systematic uncertainties from other groups working in the field) we have adopted a simple numerical interpolation scheme on a square (logarithmic) grid in the $(\varepsilon,P)$ energy–density/pressure two–dimensional space.

We fix the low–density and high density limits of the intermediate band as specified shortly. At each limit we establish a window of acceptable pressures. At the low density end this is read off from the data $(\rho,E/A)$ in the references mentioned in section 2. At the high-density end it is computed from the parametrization of [21] given in Eq. (3). We generate random functions on the grid that interpolate between the small window of maximum pressures for the low-density EoS and the small window of pressures for the high-density EoS. The conditions of causality and monotony translate into

[TABLE]

which we simply impose by Von Neumann’s rejection method, discarding all generated EoS that do not satisfy this criterion for every $\varepsilon$. We impose no condition of smoothness nor monotony to the squared speed of sound (the derivative of the $P(\varepsilon)$ EoS function). This means that the speed of sound can vary dramatically, changing quickly from large values near 1 (strong repulsive interactions) to low values near zero (which would signal phase transitions, broadly expected in cold nuclear matter [26, 27, 28]). Since we are not constraining the number and nature of phase transitions in the intermediate density region, where several alternatives have been proposed (hadron phases with strange hyperons, 2SC condensates, various chiral phases, etc; see for example [26] for hybrid stars) we allow our family of EoSes to have widely disparate sound speeds as long as they are consistent with the minimum requirements in Eq. (4)  222Within General Relativity, we know that the EoS will have to be steep, with $c^{2}_{s}$ near 1 and definitely above $1/3$ in at least part of this intermediate region [29, 30, 31, 32], in order to support $2M_{\odot}$ neutron stars. .

We can easily determine the boundary of the region for which the constraint of Eq. (4) can hold true. The construction, that we very briefly describe, is shown in figure 1. Starting at the low–end of the low–density pressure interval, $P_{1}$, we take a line $L^{(0)}_{\rm low}$ of vanishing slope as the first segment of the boundary. The length of this straight boundary segment parallel to the energy–density axis is determined by its crossing with the steepest possible line $L^{(1)}_{\rm low}$ dropped from the low–end of the high–density pressure window, $Q_{1}$ in the figure. Algorithmically, we calculate the ($Q_{1}$,$Q_{2}$) interval from Eq. (3) as shown in figure 2. (Physically, the line $L^{(0)}_{\rm low}$ would correspond to a phase transition immediately following chiral perturbation theory, finishing in $L^{(1)}_{\rm low}$ a phase of very repulsive interactions.) Analogously, the upper density end of the band is determined by a maximum–slope line $L^{(1)}_{\rm high}$ from $P_{2}$ at low density, followed by a flat line $L^{(0)}_{\rm high}$ to $Q_{2}$ at high density.

Then, a grid is laid out in that band and for each energy density we choose a random value of the pressure from the values on the grid. The inequalities (4) are tested at each point in the grid and if the function fails the test at any point, it is rejected. The remaining EoS form a set that satisfies all constraints. We have executed the program numerous times and stored the resultant set of EoS at the website described in section 4.

The low–density band below the $P_{1}$ and $P_{2}$ points (these furnishing the starting pressure interval for the interpolation to intermediate densities) in figure 1, is itself also randomly sampled from the lowest density by means of a random walk, imposing only the same very basic conditions. Here pressure grows very slowly between the two lines limiting the band, so exceeding the speed of light is not a concern ($c_{s}$ is typically 0.01-0.1). But in each step, in addition to not exceeding the side lines of the band, monotonicity of the pressure is also imposed. Additionally, if a point falls beyond the limiting lines $L^{(i)}_{\rm low}$ or $L^{(i)}_{\rm high}$ of the band it is also rejected. Those limiting lines are alternatively taken from each of the low–density sets available in the literature, at fixed order in perturbation theory. However they cover both cutoffs: we simply take a band broad enough that the uncertainty introduced by varying $\Lambda$ (or $R$) between the two quoted values is completely accounted for. If the computed line violates either causality or monotony, it is truncated.

For the high–density band, the procedure is similar. We stop the computation of the intermediate energy band upon reaching the line joining the pressures $Q_{1}$ and $Q_{2}$ corresponding to an approximate matching baryon chemical potential $\mu_{\rm match}=2.6$ GeV or $\mu_{\rm match}=2.8$ GeV as computed from $Q_{2}$ (to show sensitivity to the choice) beyond which the energy density is computed by perturbative QCD. By construction we end the intermediate band at a pressure that can be matched to the pQCD one. Finally, this last band is also sampled with a random walk from lower to higher densities, as the low–density one.

In figures 3 and 4 we show some examples of equations of state obtained following this procedure. These example EoSes span a band in $(\varepsilon,P)$ space that is broader than in [21] precisely because we have imposed no observational constraint from astrophysics.

4 Software package and sampled EoS

We have set up a website to publicly distribute this nEoS project at

http://teorica.fis.ucm.es/nEoS/ ; the website provides a “quick–start” menu with a few typical EoS constrained by alternative chiral potentials (viz. those from figures 3 and 4, classified by the order of the chiral potentials from NLO up to partially N4LO), with the interaction cutoff at two different momenta when available, and with the matching point to pQCD slightly varied at 2.6 or 2.8 GeV alternatively, to test the sensitivity to this parameter. This serves the purpose of code testing and can be useful for illustration purposes in various contexts. In this section we will use a sampler of these “quick–start” tables of the pair $(\varepsilon,P)$ to informally show what the user should expect of the nEoS sets.

We have also provide tables in .tar.gz files containing a larger sample of EoS for use in production work. They are again classified according to the chiral potential employed at low density (all of them match the essentially unique high–density computation of pQCD), and by the momentum cutoff up to which the ChPT computation is used.

A delicate choice is how to set the density of points in the grid: an excessive number of points will require much time to interpolate the EoS for possible future massive computational applications in general relativity and modifications thereof. (Practitioners there are often used to very fast polytrope EoS parametrizations and it is not clear how well the interpolated EoS will perform with a variety of computer codes.)

On the other hand, a low-density grid limits the sample to a small number of distinguishable equations of state (as there is a limited number of $(\varepsilon_{i},P_{j})$ pairs that can be sampled). We have opted for a grid with modest size, about 50 points along the $\varepsilon$ axis (logarithmically distributed between $2.6\times 10^{2}$MeV/fm3 and $4.3\times 10^{5}$MeV/fm3). This is enough to allow for flat plateaux that represent first-order phase transitions, as seen in the samples shown in figures 3 and 4, and to reach high enough energy density, above (3 GeV)4, to allow credible matching to pQCD.

On the negative side, this choice of a relatively sparse grid constrains the number of different EoS that can be constructed to $O(10^{5})$. To reduce to a minimum the probability that any two of the EoS that we generate and provide can be equal, we have limited the number to $10^{3}$ for each of the sets. The probability of repeating an EoS in the sample is then very small (a short script can be executed to check for repetitions).

5 Tests of the generated Equations of State

We have solved the Tolman-Oppenheimer-Volkoff equations of hydrostatic equilibrium to provide a feeling for these bands of EoS, choosing one of the sets of example EoS (that with a low–density limit taken from [14].

The resulting tidal deformability [34] and mass–radius diagrams are plotted in figures 5 and 6. While the absolute deformability $\lambda$ is the proportionality constant between the quadrupole tensor deformation of the star $Q$, and the tidal field causing it, $\mathcal{E}$, i.e. $Q_{ij}=\lambda\mathcal{E}_{ij}$, it has become customary to normalize it differently. Following the aLIGO-VIRGO discovery of GW170817, the literature has been skewing towards the dimensionless quantity $\Lambda=\frac{\lambda}{M^{5}}$, and this is what we plot in figure 5. (the gravitational wave data constrains the mass-weighted tidal deformability $\tilde{\Lambda}$ that is a property of the binary system, not of each individual star; when $M_{1}=M_{2}$ in the binary system, $\tilde{\Lambda}\equiv\Lambda$).

The dimensionless tidal deformability varies across a wide range. It is quite uncorrelated from the neutron star mass and radius, so that it can be thought of as a third, independent observable: in GR this can be used to further constrain the EoS. It would be interesting to repeat the computation of the tidal deformability in modified gravity along the lines of [35] with the nEoS sets to see how broad the swath of lines is and whether a constraint on the allowed separation from Einstein’s theory is possible with the full nEoS sets or, on the contrary, the maximum mass of a neutron star remains the more constraining observable.

For comparison we also added the computation of [1] with the in-medio counting of [17] (but for which no uncertainty band has been yet calculated).

The two sets shown (with dashed and dotted lines respectively; they correspond to sets 7a and 7b in our collection at http://teorica.fis.ucm.es/nEoS) have slightly different matching baryon chemical potentials $\mu_{\rm match}=2.6$ and 2.8 GeV respectively, for the transition to pQCD matter. The sensitivity to this choice is small unless one lowers $\mu_{\rm match}$ to significantly lower values, say 2 GeV; then the pQCD EoS becomes quite steep (unlike the slow logarithmic running for higher densities) and this affects the interpolated intermediate energy density.

The plotted lines look somewhat more ragged than usually reported computations; this is because of the nonsmooth EoS in the set (as we are not employing soft polytrope–like interpolation at all), so the various curves occasionally change direction.

The basic feature of a maximum neutron star mass below three solar masses is there, which can help constrain variations of General Relativity. The differences among computed radii (and consequently among tidal deformabilities) are larger than thought but we hope upcoming theory efforts will reduce this. The radius-independent quantities such as the Love number $k_{2}$ are much less-sensitive. It is straightforward to extend these examples with the EoS listed in http://teorica.fis.ucm.es/nEoS/ or to proceed to systematical investigations with the larger sets with 1000 equations each, also stored there.

6 Conclusions

6.1 Results

In this concise work we have reported an extensive set of $\beta$-equilibrated neutron–star matter Equations of State based on first principles, on Chiral Perturbation Theory and on perturbation theory in QCD. In addition to avoiding systematic model biases as much as currently possible, we have not imposed any astrophysical constraints, so that the family of EoS truly samples the current state of the art hadron theory prediction.

This means that uncertainty bands for observables calculated from our nEoS family of equations will be broad, as illustrated with the examples in figures 5 and 6, and we look forward to theory improvements to reduce these.

We provide in all 18000 equations of state that satisfy all theory constraints, parceled in subsets of 1000 EoS each. These are organised as follows.

• •

Sets 1 through 3: low density band constrained by [11] at NLO, N2LO and N3LO respectively, spanning cutoffs of 450 and 600 MeV. Provided as two subsets each, one with $\mu_{B}=2.6$ GeV (subsets 1a, 2a and 3a) and another with 2.8 GeV (subsets 1b, 2b, 3b).

• •

Sets 4 through 6: low density band constrained by [12] at NLO, N3LO and N4LO, with cutoff parameters $R=0.9$ and 1 fm. Also provided as two subsets each, likewise denoted by (a) and (b) lettering.

• •

Set 7: low density band constrained by [14], subsets 7a and 7b as above depending on the $\mu_{B}$ matching to pQCD.

• •

Sets 8 through 10: low density band constrained by [13], at NLO, N2LO, N3LO respectively and also divided in (a) and (b) subsets each.

A systematic comparison of all the sets with respect to several simple observables is envisioned for the near future, but it may be that a different collaboration will want to undertake the endeavor, perhaps more efficiently than us, since the automation of neutron star calculations is a subfield in itself [36], so we have opted for releasing the sets for public examination.

For an example of such work, let us plot in figure 7 a comparison between two same–order subsets, namely (1a) and (4a) in the nEoS site, both using NLO input from chiral computations from different collaborations [11, 12].

We can appreciate in the figure’s left plot, the rendering of the nEoS sets based on (1a) (Sammarruca et al.) and (4a) (Hu et al.) how the first starts at higher pressure, so that the low–density band only reaches up to around 200 MeV/fm3 in energy density, while the second, that is initially less steep, has a low–energy band extending up to 375 MeV/fm3. (The low–density band becomes doubtful when the typical neutron momentum approaches the cutoff). At intermediate densities, both subsets start looking very similar, but subset (4a) never regains the stiffness lost at initial energy densities.

In consequence, subset (1a) reaches higher neutron–star masses and larger radii for a typical $1.4M_{\odot}$ neutron star: we find $R_{1a}\simeq(12.8\pm 2.1)$km, $R_{4a}\simeq(9\pm 1.8)$ km respectively. A lot of predictive power can be gained from hadron physics alone if the low–density investigations of the EoS in Chiral Perturbation Theory converge to a solidified result.

6.2 Discussion and outlook

There are several ways in which this project could be enhanced.

One aspect would be to incorporate results from lattice gauge theory. We have not done this yet because currently lattice gauge calculations have to be extrapolated from zero baryon density by using a Taylor expansion around $\mu_{B}=0$ and assuming analyticity in the baryon chemical potential [37]. This may be reliable for small $\mu_{B}$ at finite temperature, but still relies on a further model–dependent Hadron Resonance Gas extrapolation to reach zero temperature from the actual lattice calculations at $T\sim 100$ MeV. The method can well fail at higher densities, where singularities might appear, particularly near a phase transition. Additionally, the region of parameter space covered, after all extrapolations, would coincide at best with that where we have effective (chiral) theory working. The only real gain by using lattice gauge theory would be to have the low-density EoS directly connected to $\alpha_{s}$, the strong coupling constant at the quark level, without resource to the coefficients of the effective Lagrangian; even this much looks like a far goal given the present status of work, that is making strides only for heavy quarks [38]. In the practice of ChPT, the coefficients are fitted to $NN$, $\pi N$ and other low–energy data in laboratory experiments. For our purposes, since both approaches, QCD computation of the chiral coefficients, or accelerator data fitting, are independent of General Relativity, it does not make a difference as long as the Equation of State becomes known. We will probably revisit this decision at a later date when/if progress in lattice gauge theory brings significant advantage.

A second caveat worth mentioning is that we have not incorporated the outer crustal EoS, described by (significantly messy) nuclear phases with drops of actual nuclei, perhaps “pasta” phases, as opposed to a simpler, fixed-$r$ homogeneous, $\beta$-equilibrated hadron continuum as described by chiral perturbation theory. It is well-known that this omission is of no concern for the mass of the star, as the crust is relatively much lighter, but it does bring–in an underestimate of the star radius by about 0.5 km [39]; but given that the uncertainty in experimental data is still about 2 km for the radius, and the spread of our EoS band gives easily a theory uncertainty of 4 km, which is eight times bigger, it is not urgent to incorporate it. We may do so at a future date if a way of providing a relatively model–independent estimate is pointed out to us. Meanwhile, we point the reader out to [40] for an appraisal of the nuclear–end EoS.

At the highest densities, it is expected [41] that quark matter is in the Color-Flavor-Locked (CFL) phase, with Cooper pairing involving all three quark flavors. The pressure then contains a term proportional to the new scale, the pairing gap $\Delta$, of the form $P_{2}\propto\Delta^{2}\mu_{B}^{2}$. This is subleading to the $\mu_{B}^{4}$ scaling of the ungapped pQCD computation in [20] and is presently not included.

Finally, a compelling avenue for future work is to try to incorporate the effect of a finite temperature into nEoS, which is necessary for neutron star mergers, perhaps using data from heavy ion collisions [42].

We are looking forward to addressing some of these and other interesting topics in future revisions of the EoS sets if they become of use for the community of researchers trying to test General Relativity with neutron stars.

Acknowledgments

FJLE acknowledges the assistance of Iván Sayago Galván in the early stages of this investigation. Work supported by grants MINECO:FPA2016-75654-C2-1-P (Spain); the EU COST Action PHAROS (CA16214); and by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics under Award Number #DE-FG02-05ER41375. AW acknowledges support by the Schrödinger Fellowship J 3800-N27 of the Austrian Science Fund (FWF).