Unbiased interpolated neutron-star EoS at finite T for modified gravity studies

Neutron stars and their mergers provide the highest-density regime in which Einstein’s equations in full (with a matter source) can be tested against modified theories of gravity. But doing so requires a priori knowledge of the Equation of State from nuclear and hadron physics, where no contamination from computations of astrophysics observables within General Relativity has been built in. We extend the nEoS uncertainty bands, useful for this very purpose, to finite (but small) temperatures up to T=30\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T=30$$\end{document} MeV, given that the necessary computations in ChPT and in pQCD are already available in the literature. The T-dependent band boundaries will be provided through the COMPOSE repository and our own website.


Eqn. of State of Neutron Matter at finite T
The Equation of State of neutron star matter is the basic microscopic input necessary for multiple computations of interest in astrophysics [1]. At low, nuclear densities it is well constrained from laboratory data [2], but there is no universally accepted theory pathway for intermediate energy densities of order the hadronic scale, typically (100 MeV) 4 and above.
Numerous models with various matter contents have been deployed for this intermediate density region, for example including hyperons [3], hybrid approaches [4] including quark matter [5], employing the holographic conjecture linking QCD to a gravity dual [6], and more.
A common approach is to employ astrophysics observables to constrain the equation of state, that therefore contains information from nuclear and particle physics, but also of General Relativity (GR), used to eliminate parts of the otherwise allowed parameter space [7][8][9][10]. The inconvenience of this method is that the general relativistic interpretation a e-mail: fllanes@fis.ucm.es (corresponding author) of the astrophysical observables employed to constrain the equation of state becomes entrained in the EoS; thereafter, it is not safe to employ such equation to constrain modifications of GR.
For example, there is great interest in constraining f (R) theories and variations thereof [11][12][13], or closely related scalarizations [14][15][16], or one could ask about whether the parameters of GR are sensitive to a large stress-energymomentum density [17], or pursue nonlocal gravity with compact objects [18] among many possibilities.
But the band of allowed equations of state compatible with these theories is presumably larger (since they have additional parameter freedom in the gravity sector) than in General Relativity. Thus, employing reduced uncertainty EoS-bands incorporating astrophysical observables leads to improper constraints on the beyond-GR theories.
A separate question, that we do not address, is how to avoid assuming General Relativity at short scales in the interpretation of the low-energy nuclear data themselves (for example, Fisher and Carlson [19] have looked at constraining Non-Minimally Coupled Gravity from nuclear properties). In this work we adopt the line of thought that microscopic nuclear physics can be read-off an inertial frame with negligible tidal forces, and thus only the large accumulation of mass at neutron-star scales can cause separations from General Relativity (in a large region of intense stress-energy T μν and intense curvature).
For such application it seems a necessity to provide EoS that are state of the art from the point of view of nuclear and particle physics, but that are free of astrophysical bias, if they are to be used to employ neutron star data for constraining beyond-GR theories. We have provided precisely such family of Equations of State, with controlled theoretical uncertainty, relying only on first principle approaches (causality, thermodynamic stability) and perturbation theory in the appropriate density domains (chiral perturbation theory, perturbative QCD) in the nEoS sets [20] (that can be downloaded from https://teorica.fis.ucm.es/nEoS/). Meanwhile, the detection of a neutron-star merger and the spur of improved fits to ejecta [21] and simulations of the matter-remains after the collision that followed [22,23] have made the need for modern EoS at finite temperature more poignant. We thus look at extending those sets to finite-T .

Low density ChPT input
Whereas there are numerous equations of state for neutron stars extracted from chiral perturbation theory (ChPT) at T = 0, full-temperature calculations do not abound.
In this work we have leaned on the computation of the Idaho group [24] for the low-density, chiral-Lagrangian tail of the computation. Pure neutron matter is a good approximation in this low-density regime as it is not too far from equilibrium. As it is close to the nonrelativistic limit, energy densities are dominated by the neutron mass contribution and it is essential to keep m n .
Their obtained Helmholtz potential depends on a N 4 L O treatment of the nucleon-nucleon potential; the three-body N N N potential is only partially treated. Their computation is for pure neutron matter alone, and this is not a bad approximation for the low to moderate densities where the approach is valid. The authors provide the Helmholtz free energy per nucleon, from which P follows directly (e.g. [28] for discussion) The convergence of the order by order chiral computations is very nice, and the effect of the temperature, interestingly, is seen to be opposite in E/A and F/A. In fact, Sammarruca, Machleidt and Millerson [24] find that E/A increases with the temperature at fixed density, probably reflecting the fact that T softens the step of the Fermi-Dirac occupation function, upgrading some neutrons from below to above the Fermi sea level of cold matter. Since neutrons of higher momentum have larger kinetic energy and see stronger (repulsive) chiral interactions, it is natural to expect higher energy. However, this is compensated by the second term of Eq. (1), so that the free energy is actually decreasing with temperature. Its derivative, however, which is the relevant quantity to compute the pressure, is not as sensitive to small T except at very low densities.  [24][25][26][27], with T increasing from bottom to top in 10 MeV steps. As the number density rises, the uncertainty grows and the bands do overlap, so that the effect of temperatures much smaller than the densities is less visible The resulting EoS 1 , that we adopt as our low-density bands, based on [24] are shown in Fig. 1.
There are many nuclear model-dependent computations, such as [29] based, for example, on Skyrme-type interactions, that corroborate the increase in E/A with temperature. Most computations are reasonably valid in the low-momentum regime, see Fig. 2, but the chiral perturbation theory bands that we adopt showcase the necessity of assigning systematic errors to these other computations. In the future, one might even be able to dispose of the need for a nuclear potential altogether and directly compute the equation of state from low-energy nucleon-nucleon scattering data, further reducing systematics [30].

High density pQCD input
While hadron interactions are strong and complicated at intermediate momentum, where the coupling is large, for asymptotically large density the smallness of α s (μ) allows the deployment of perturbative Quantum Chromodynamics (pQCD). Thus, all of our EoS terminate inside the highdensity band obtained from pQCD computations and shown in Fig. 3.
That band is obtained from perturbative computations at finite density following the Helsinki group work [31][32][33][34] and also [33,35] for the T → 0 limit. At such high chemical potentials of order 2.5 GeV and above, the quark masses are negligible: weak equilibrium is adopted consistent with flavour symmetry and charge neutrality.
Though we have checked our results with the simple approximate parameterization proposed in [32], our nEoS-T Fig. 2 Comparison of the N4LO ChPT computation [24] and the Skyrme-based computation of [29]. While the comparison is reasonable, if a bound on a beyond-GR (e.g. post-Newtonian) coefficient is obtained, what systematic error should be ascribed to the model computation? For this application it is better to have known uncertainty bands even if the EoS is less accurate for certain properties within General Relativity. Also, having a well-defined band can help model-makers understand the agreement of computed EoS with certain basic principles band follows the detailed formulae earlier reported in [33,34] that we have reprogrammed.
The span of the pQCD uncertainty band is an attempt at controlling the systematics of the computation by varying the renormalization scale ∈ (0.5, 2)¯ around the valuē The scale variation should be reabsorbed in the nonconformal terms of higher, non-considered orders of perturbation theory without affecting the energy density or pressure. However, a dependence arises which is an artifact of truncating the weak-coupling expansion. There is some arbitrariness in the choice of scale-uncertainty interval that becomes smaller with a higher order of perturbation theory (the exact series should be independent of the renormalization scale). What this band really does is to orient us in the validity of the perturbative computation: note in the figure how, for chemical potentials below 2 GeV, the uncertainty swells and the computation becomes unusable. In practice, we vertically slice the band at a baryon chemical potential μ B = 2.67 GeV (this is the nearest point of our grid to the characteristic quark chemical potential that the Helsinki group has been employing, about μ q = 0.87 GeV; for comparison, another line at 2.5 GeV is given). The intersection of that vertical line with the two solid lines (blue and red online), that provide the scale sensitivity in the interval (0.5 , 2 ) act as "goalposts" through which any EoS coming from lower density has to pass; those that do not, are automatically discarded. It marks the border between the high-density region (computed in pQCD) and the intermediate density region (interpolated). This particular point of our grid is chosen as a compromise between the necessity of the uncertainty being small enough but the density being as low as possible to reduce its difference to that of actual neutron stars.
These pQCD EoS entail a three-loop computation of the energy per quark. Recently, attempts have been made at narrowing the uncertainty band by means of the renormalization group [36] but we have been conservative and not incorporated these results yet to avoid introducing some bias or model dependence, at least until further research and independent confirmation reassures us that the procedure is safe. The same can be said of Dyson-Schwinger [37] related approaches, that may be insightful but not yet optimized to provide uncertainty bands for beyond-GR searches.
It is clear to us that physical neutron stars in General Relativity will not reach such high baryon chemical potentials as described in this section. But there are two reasons to anchor the computation of the EoS in the pQCD one. The first is that modified-gravity theories might actually reach such densities, so we need to provide them with input there. The second, broader and also applicable to GR, is that the EoS of hadron matter at attainable densities is constrained by its necessary behaviour at higher ones, due to causality and stability ( d P dε ∈ [0, 1)).

The nEoS sets at finite Temperature
We are now ready to mount the interpolation. We take the low-density ChPT band and start sampling it from lower to higher values of n; after its maximum density is reached, we continue sampling a broadening band sandwiched between the extreme lines allowed by the conditions of causality d P/dε ≤ 1 and thermodynamic stability d P/dε ≥ 0 between the ChPT and the pQCD bands. This is accomplished by Von Neumann's rejection, with the program described in [20]. If at any step in the numeric advance from (ε i , P i ) → (ε i+1 , P i+1 ) any of those two conditions is even locally violated, the EoS is rejected and we start anew.
Thousands of such equations, exploring all the systematics of the ChPT cutoff, the pQCD renormalization scale, and the different computations of the low-density limit, are provided for zero temperature in our website http://teorica.fis.ucm.es/ nEoS.
In this work, the extension to (small) finite temperature, we adopt a different strategy. Instead of providing multiple samples of the band (all workable EoS compatible with hadron constraints), we distribute the extreme boundary lines of the nEoS bands formed from (a) the N 4 L O chiral perturbation theory computation at 0, 10, 20, 30 MeV in section 2; (b) the pQCD computation described in Sect. 3; and (c), the interpolation between both limits for the corresponding temperatures, satisfying stability and causality.  Figure 4 shows the band up to the end of the intermediate interpolated zone (the higher density pQCD part is not show to avoid reducing the scale too much; it is the same as in Fig. 3).
In addition to our computer server just mentioned, we will be making these band extremes available through the COMPOSE website https://compose.obspm.fr/ where many other EoS sets for neutron stars can be found. We also provide the geometric mean of the band extremes at largest and lowest pressures, as an example of a typical EoS through the middle of the allowed band.
Additionally, in this article we show a few examples of allowed EoS that pass through the band in Fig. 5. These examples illustrate a feature of linearly interpolated and stepwise generated EoS: the possibility of having flat sections of zero derivative that can represent first order phase transitions. The longer these flat segments are in the graph, the larger the latent heat of the transition, with the nEoS band introducing an upper bound to the latent heat that is possible in hadron physics [39].

Comparison to the additive approach
Often used is the additive (or hybrid) approach [38] that approximates pressure and energy density at finite temperature by the T = 0 counterpart for cold nuclear matter and an additive thermal correction = cold + (T ) P = P cold + P(T ) in terms of a simple adiabatic proportionality factor , which is a constant of the approach, as is the following ratio: This constancy is an assumption that we here briefly examine; it has been recently employed in [23], together with a more sophisticated so-called M * parametrization of the relation between P(T ) and ε(T ). The authors find that the thermal pressure at T = 20 MeV is of the same size of the cold-matter pressure at saturation density, but the ratio falls rapidly so that by ε = 10ε sat , the thermal pressure is a 10% correction at most.
We have plotted Eq. (4) for T = 30 MeV in Fig. 6 following the top of the band (that is the most relevant for maximum mass neutron stars). The plot is quite bumpy, as one can identify the matching points between density regimes around 0.27 Fig. 6 Test of the additive approach to thermal effects on the neutron star EoS, evaluating Eq. (4) for T = 30 MeV and 15 GeV/fm 3 as well as the transition between causalitydominated (maximum slope of 1) and stability-dominated (minimum slope of 0) around 5 GeV/fm 3 as the pQCD phase is approached. Once these non-analytic features have been discounted (that can presumably also appear for the one true equation of state if it presents phase transitions between different states of nuclear matter), we can focuse on the "safe" regions.
They have very substantial values of the derivative of | − 1|. In the low-density regime, the variation of the function plotted, and thus the uncertainty of the additive approach, is of order 20%, which is reasonable for many applications.
The change is however of order 100% along the top line of Fig. 4 for energy densities in the few hundreds and thousands of MeV/fm 3 . This is the density region where the core of neutron stars should be in General Relativity and slightly modified theories thereof.
For the flat stretch of pressure (akin to a first order phase transition) around 10,000 MeV/fm 3 , the ratio | − 1| falls by about a factor 3, so assuming that it is constant is a poor approximation. This is a density regime, however, that will only be relevant for actual neutron star physics in theories that significantly separate from General Relativity (on the weaker side, as reaching such large densities requires to slow stellar collapse down).
So we conclude that the additive approach is reasonable for practical applications at relatively low-densities (still, large compared to T ) within GR and becomes more dubious in other circumstances.

Discussion
We have presented an extension of the nEoS sets [20] to (small) finite temperature. For this we use as input the finite-T chiral computation of [24] that determines the low-density EoS of pure neutron matter. At the high-density limit, in contrast, we have used the pQCD computation of [34] with μ u = μ s , i.e, employing SU (3) flavor symmetry; that is, we interpolate between two apparently different regimes as the weak interactions are concerned. This is a theoretical choice based on the information at hand at the low-density limit, where pure neutron matter is a good approximation, and trying to stay close to β-equilibrium at high densities, where the strange quark mass is a correction. In any case, through the central zone of intermediate densities, the range of pressures allowed at a given density due to the stabilitycausality band is very large, and this makes the choice of flavor-breaking chemical potential a correction. This can be seen, for example, from the expression valid for cold quark matter that shows the large μ expansion, with subleading terms carrying the strange quark mass, the CFL gap if such asymptotic phase is formed, the ground-state bag constant B or equivalent, etc. It is the first term that dominates the expansion, and it is flavor-blind. Future work might address βequilibrium through the entire range of densities. Perhaps useful to orient the reader as to what hadron physics currently allows at zero temperature is the traditional M(R) plot that we sketch in Fig. 7 for the usual General Relativity case.
The top, black line corresponds to the stiffest possible EoS allowed by hadronic and theoretical constraints at T = 0, the lower one (red online) to the softest (its bend is out of the plot due to the scale used). They are respectively above and below the astrophysical constraints on masses and radii for neutron stars in this (R, M) plane, within General Relativity (but perhaps not its extensions) but this is no place to describe such data and we refer to the extensive literature.
One could argue that this plot alone would exclude, within General Relativity, the two extreme lines of the nEoS band. The stiffest line reaches masses above 3 M , in contrast to even optimistic approaches that do not foresee neutron star masses above 2.6 M from statistical analysis of stellar populations [40]. The softest line on the other hand does not even reach one solar mass (!). Thus, hadron physics has much room to improve by itself. Of course, if GR is taken for granted, a large swath of the nEoS band must disappear to make the stiffest EoS softer and the softest one, stiffer, as marked in the figure.
Stronger constraints from the high-density pQCD side (for example, through extending perturbation theory with a functional approach as has recently been reported) would affect the former, while extending the ChPT band to higher densities would much affect the second. The reason is that pQCD is quite soft itself, with c 2 s ∼ 1/3, while ChPT is rather stiff, and there is no reason to expect a dramatic change of slope Basically, both extremes seem excluded in the framework of General Relativity, meaning that hadron physics computations are less competitive than astrophysical ones for now. However, modifications of GR could bring one or both of these curves back to agree with data, so they should be included in studying modified gravity with incremental improvements that do not seem far-fetched and can be expected in the next few years. Eventually there should be an intermediate-density ceiling that ChPT cannot break. But even small increases in its density reach, due to the control of higher order corrections (particularly their threeand more body parts) will have a noticeable impact.
Perhaps there will also be applicable constraints from the Relativistic Heavy Ion Collisions programme. For example, there is a constraint based on proton flow [41] that is directly relevant to symmetric matter 2 . However, extrapolation to neutron-star matter is required, channelled by a band between two parametrizations of the nuclear symmetry energy that are not obviously based on first-principles physics. We have not incorporated it following the principle of being very cautious about imposing constraints on the EoS that eventually gravity investigators can use to constrain beyond-GR theories, to avoid an unwarranted overconstraining. We may include it in the future if we become convinced that the procedure is safe.
In turn, the finite-T EoS is of importance to extract information from the neutron-star like object formed in mergers of two ordinary neutron-stars such as GW170817, that can be hot (due to friction) neutron-matter ephemeral systems supported by fast rotation. They can be used, for example [42], to bracket the maximum neutron star mass with gravitational wave data.
We have examined temperatures T = 0 through T = 30 MeV, following the existing [24] ChPT computations. The temperature dependence here is appreciable (Fig. 1): every 10 MeV, the central value of the N 4 L O chiral computation is outside the uncertainty band of the neighbour-ing temperatures. Therefore, T does have to be taken into account into computations involving densities at the nuclear scale. The high-density computations, on the other hand, are not too sensitive to T , as corresponds to the scale hierarchy T μ B (30 MeV versus 2.6 GeV and above, Fig. 3). The T -dependence becomes remarkable at lower chemical potentials where the perturbative expansion is less reliable anyway and should rather not be used. Finally, in the intermediate region (Fig. 5), the effect of the temperature is also small when compared with the large swath of parameter space allowed between the low-pressure limit of stability and the high-pressure limit of causality.
Should the user believe that that small temperature dependence can make a difference in her calculation, our work provides it through the entire density range. It is common among researchers working on modified gravity models to employ ad-hoc polytrope lines from estimates in the sixties; one can surely do better and we hope to at least provide a stepping stone. Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.

Data Availability Statement
The manuscript has associated data in a data repository. [Authors' comment: Data is available at the following web addresses https://teorica.fis.ucm.es/nEoS/nEoST.html (finite Temperature) https://teorica.fis.ucm.es/nEoS/ (zero Temperature) and will eventually be provided at the COMPOSE repository https://compose. obspm.fr/.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

Appendix
After this work was first released, a new preprint [43] by Komoltsev and Kurkela has shown how constraining causal- Fig. 8 Effect of bringing into the (ε, P) plane the constraint from causality in the (n, μ) plane, at T = 0. We plan to explore this constraint in upcoming work ity in the (n, μ) plane, together with the high-density pQCD band, further restricts (and significantly so) the band of allowed EoS. While we defer a detailed study of this welcome new information to a future investigation, we did not want to pass this opportunity to explore the effect, even if cursorily.
The new constraint, called "integral constraint" in that publication, stems from the inequality The difficulty to visualize it is that now a second thermodynamic variable, μ, in addition to P, needs to be taken into account when constructing the ε energy density. The nEoS band in the ( , P) plane, which is the useful one for T μν and the Einstein equations, becomes a more complicated tube with an additional dimension. In effect, the area under the EoS curve in the (n, μ) plane [43] has to equal the increase in the pressure, The constraint is then mapped to the (ε, P) plane by means of ε = −P + μn at fixed number density n.
We have made an attempt at a meaningful representation in Fig. 8.
The dotted lines in the figure represent the nEoS band with available information from the (ε, P) plane only. The tree shaded bands are further restricted, in a very significant way, especially at low pressures, by (n, μ)-plane causality, and differ from each other by the different μ H chosen upon matching at the pQCD band (blue, orange and red dots on the right upper corner, respectively). The effect is important and will be incorporated in future issuances of the sets.