From the microscopic to the macroscopic world: from nucleons to neutron stars

TL;DR

This paper reviews how recent neutron star observations, especially gravitational waves, constrain nuclear interactions by connecting microscopic nucleon interactions to macroscopic neutron star properties.

Contribution

It provides a comprehensive review of first-principles calculations of the equation of state of nuclear matter and their implications for neutron star observations.

Findings

Constraints on nuclear interactions from gravitational wave data

Correlation between microscopic models and neutron star properties

Advances in first-principles calculations of nuclear matter

Abstract

Recent observations of neutron-star properties, in particular the recent detection of gravitational waves emitted from binary neutron stars, GW 170817, open the way to put strong constraints on nuclear interactions. In this paper, we review the state of the art in calculating the equation of state of strongly interacting matter from first principle calculations starting from microscopic interactions among nucleons. We then review selected properties of neutron stars that can be directly compared with present and future observations.

Figures

Tables

Table 1. Table 1: Fitting parameters for the neutron matter EOS defined above. The values of E s ​ y ​ m subscript 𝐸 𝑠 𝑦 𝑚 E_{sym} and L 𝐿 L are obtained by fitting the AFDMC results. The N 2 LO parameter represent the higher, middle and lower part of the read band of Fig.  4 .

Hamiltonian	$E_{\mathrm{sym}}$	$L$	$a$	$\alpha$	$b$	$\beta$
	(MeV)	(MeV)	(MeV)		(MeV)
AV8’	30.5	31.3	12.7	0.49	1.78	2.26
AV8’+UIX	35.1	63.6	13.4	0.514	5.62	2.436
N2LOup (1.0 fm)	34.8	56.2	13.19	0.51	5.66	2.12
N2LOlow (1.0 fm)	30.2	24.4	13.87	0.59	0.36	-0.05
N2LOmid (1.0 fm)	32.5	40.5	13.13	0.51	3.41	1.99

Equations

$$\frac{dP}{dr}$$

$$\frac{dm(r)}{dr}$$

$$\Psi(\tau)=\exp\big{[}-H\tau\big{]}\,\Psi_{v}\,,$$

$$H_{\mathrm{nuc}}=\sum _i\frac{p_i^2}{2m_N}+\sum _{i<j}{v}_{ij}+\sum _{i<j<k}{v}_{ijk}.$$

$$E_{PNM}(\rho )=a{\left(\frac{\rho }{{\rho }_0}\right)}^{\alpha }+b{\left(\frac{\rho }{{\rho }_0}\right)}^{\beta },$$

$$E_{\mathrm{sym}}(\rho )=E_{\mathrm{PNM}}(\rho )-E_{\mathrm{SNM}}(\rho ),$$

$$E_{sym}$$

$$L$$

$$E_{\mathrm{SNM}}(\rho )=E_{\mathrm{PNM}}(\rho )-E_{\mathrm{sym}}(\rho ).$$

$$E(\rho ,x)$$

$$E_{sym}(\rho )=C{\left(\frac{\rho }{{\rho }_0}\right)}^{\gamma },$$

$$L=3{\rho }_0\frac{\partial E_{\mathrm{sym}}(\rho )}{\partial \rho }=3{\rho }_0\frac{C\gamma }{\rho }{\left(\frac{\rho }{{\rho }_0}\right)}^{\gamma }\to L\sim C.$$

$$C=a+b+16.$$

$$P={\rho }^2\frac{\partial E_{\mathrm{SNM}}}{\partial \rho }{|}_{\rho ={\rho }_0}=0.$$

$$\gamma =\frac{a\alpha +b\beta }{C}=\frac{a\alpha +b\beta }{16+a+b}.$$

$$L=3(a\alpha +b\beta ).$$

$$E(\rho ,x)=E_{\mathrm{PNM}}(\rho )+(16+a+b){\left(\frac{\rho }{{\rho }_0}\right)}^{\frac{a\alpha +b\beta }{16+a+b}}\left[(1-2x{)}^2-1\right].$$

$$ϵ$$

$${\mu }_y$$

$${\mu }_n-{\mu }_p=4(1-2x)E_{sym}(\rho ).$$

$${\mu }_n-{\mu }_p={\mu }_e={\mu }_{\mu }$$

$${\mu }_e=(m_e^2+{\hslash }^2{k}_F^2{)}^{1/2}\approx \hslash (3{\pi }^2\rho x_e{)}^{1/3},$$

$${\mu }_{\mu }=[m_{\mu }^2+{\hslash }^2(3{\pi }^2\rho x_{\mu }{)}^{2/3}{]}^{1/2}.$$

$$x=x_e+x_{\mu }.$$

$$c_S=\sqrt{\frac{\partial p(ϵ)}{\partial ϵ}}.$$

$$F_{\mathrm{np}}(\rho ,x,T)=F_{\mathrm{virial}}(\rho ,x,T)g+F_{\mathrm{deg}}(\rho ,x,T)(1-g),$$

$$g=1/(1+3z_n^2+3z_p^2),$$

$$F_{\mathrm{deg}}(\rho ,x,T)$$

$$E_{\mathrm{sym}}=h(\rho )E_{\mathrm{PNM}}+[1-h(\rho )]E_{\mathrm{NS}}(\rho )-F_{\mathrm{Skyrme}}(\rho ,x=1/2,T=0),$$

$$h=\frac{1}{1+\exp [\gamma (\rho -3/2{\rho }_0)]},$$

$$\Lambda =\frac{2}{3}{k}_2{\left(\frac{c^2}{G}\frac{R}{M}\right)}^5,$$

$$\tilde{\Lambda }=\frac{16}{13}\left(\frac{(m_1+12m_2)m_1^4{\Lambda }_1}{m_{\mathrm{tot}}^5}+\frac{(m_2+12m_1)m_2^4{\Lambda }_3}{m_{\mathrm{tot}}^5}\right).$$

Full text

1 Introduction:

The idea that astronomical observations may provide insight into nuclear interactions originated in the 1950s [1, 2]. There have been a few landmark neutron-star observations in the past decade which can be directly connected to the interaction between nucleons. Two neutron-star mass measurements obtain results $\sim 2$ solar masses [3, 4] (see also [5]). Simultaneous information on neutron-star masses and radii is becoming available [6, 7] (see updates in Refs. [8, 9, 10]). Finally, the first detection of gravitational waves from a binary neutron-star merger, GW 170817, has also provided mass and tidal deformability constraints [11], and its electromagnetic counterpart has demonstrated that neutron-star mergers are an important source for r-process nuclei. These advances in neutron-star observations provide a unique opportunity to improve our knowledge of nuclear physics. Although the corresponding length scales are separated by many orders of magnitude, properties of neutron stars and nuclei are strongly connected. In particular, the equation of state (EOS) of the crust and the outer core is one of the main ingredients for neutron-star structure, determining radii, tidal deformabilities, and other properties of neutron stars. The EOS is obviously related to nuclear forces and properties of nuclei.

At the femtometer scale, it is of fundamental importance to understand how nucleons interact, how neutrinos interact in nuclear matter, how nuclei are formed and how their properties emerge. All these are among the main ingredients for astrophysical simulations, including simulations of supernova explosions and the evolution of neutron stars. Those simulations predict or explain observations of phenomena which can have exameter (100 light-year) length scales, including electromagnetic and gravitational-wave emission, nuclear abundances, ejecta from binary neutron-star mergers, and others; see Fig. 1 for a diagram summarizing the connections between microphysical models, simulations, and observations. There is a strong connection between the microscopic and macroscopic world, separated by many orders of magnitude. Information in Fig. 1 flows both ways: microphysical models can be used to describe observations, and observations constrain our knowledge of nuclear interactions, nuclear structure, and nuclear reactions.

In this paper, we will start by briefly discussing how nuclear forces are constructed and tested, how they are used to calculate the EOS, and then how the EOS is used to predict selected properties of neutron stars. We will then provide details regarding how Bayesian inference can be utilized to attack the “inverse problem”: the problem of constraining our model parameters from observational data.

2 The EOS

Most of the static and dynamical properties of neutron stars can be calculated once an Equation of State (EOS) describing the matter inside the stars is specified. At very low densities, in the outer crust of neutron stars, the matter is mainly composed of a lattice of ordinary nuclei in the iron region. With increasing density, the neutron chemical potential also increases and nuclei become extremely neutron-rich. At the interface with the inner crust, the neutron chemical potential is sufficiently high so that neutrons start to drip out of the nuclei, and the nuclei start to be surrounded by a sea of neutrons. Eventually, at the bottom of the inner crust, the geometry of nuclei begins to be deformed, forming the so-called “pasta” phase. At about half the saturation density, $\rho_{0}=0.16$ fm${}^{-3}\approx 2.7\cdot 10^{14}$ g/cm3, the nuclei completely melt, and the neutron-star core begins. Here, matter is composed of a uniform liquid of a large fraction of neutrons with a few protons, electrons, and eventually muons. At even higher densities, above $\approx 2\rho_{0}$, the composition of matter is basically unknown. In this inner core of the neutron star, many scenarios for the state of matter have been suggested, for example the formation of hyperons [12], quark matter [13], or other more exotic condensates. Fig. 2 shows the nature of these layers and gives a summary of properties of a typical neutron star.

The full description of the EOS in the whole range of densities encountered inside a neutron star, i.e. up to several times nuclear saturation density, is a formidable task. Especially for the inner core, most EOS used in astrophysical simulations necessarily make some model-dependent assumptions. However, once an EOS $\epsilon(p)$ is specified, the mass-radius relation of a non-rotating neutron star can be easily calculated by solving the Tolman-Oppenheimer-Volkoff (TOV) equations:

[TABLE]

where $P$ and $\epsilon$ are the pressure and the energy density (including the rest mass energy density contribution), $m(r)$ is the gravitational mass enclosed within a radius $r$, and $G$ is the gravitational constant. The solution of the TOV equations for a given central density gives the profiles of $\rho(r)$, $\epsilon(r)$ and $P(r)$ as functions of radius $r$, and also the total radius $R$ and mass $M=m(R)$, specified by the condition $P(R)=0$. By varying the input central density, the mass-radius (MR) relation is mapped out. The resulting MR relation is in one-to-one correspondence with the input EOS (as summarized in Fig. 3), which offers the possibility to measure the EOS, and thus the properties of strongly interacting matter up to a few times saturation density, by mapping out the MR relation observationally. In addition to properties of isolated neutron stars, also other dynamical properties of neutron stars and binary neutron-star mergers depend on the EOS, including the gravitational-wave spectrum, the amount of ejected material, and others.

An ideal starting point for the EOS of neutron-star matter is the EOS of pure neutron matter (PNM), which is a homogeneous system that contains only neutrons. The EOS of PNM can be calculated using sophisticated many-body methods once a nuclear Hamiltonian is specified. Among other methods, these include variational methods based on the cluster expansion [21], many-body perturbation theory [22], the coupled-cluster method [23], and Quantum Monte Carlo methods [24, 25]. In this paper, we will focus on recent results obtained with the Auxiliary Field Diffusion Monte Carlo (AFDMC) method, which was originally introduced by Schmidt and Fantoni [26], and is ideally suited to study neutron matter [27, 28].

The main idea of QMC methods is to evolve a many-body wave function in imaginary-time:

[TABLE]

where $\Psi_{v}$ is a variational ansatz of the many-body wave function and $H$ is the Hamiltonian describing the system. In the limit of $\tau\rightarrow\infty$, $\Psi$ approaches the ground-state of $H$. The evolution in imaginary-time is performed by sampling configurations of the system using Monte Carlo techniques, and expectation values are evaluated over the sampled configurations. For more details see for example Refs. [28, 29, 30, 31, 32, 33].

In addition to the many-body method, one needs the nuclear Hamiltonian as input, which describes the interactions among nucleons. In AFDMC, nuclei and neutron matter are described by non-relativistic point-like particles interacting via two- and three-body forces:

[TABLE]

The first two-body potential that has been extensively used with the AFDMC method is the phenomenological Argonne AV8’ potential [34], that is a simplified form of the Argonne AV18 potential [35]. Although simpler to use in QMC calculations, the AV8’ potential provides almost the same accuracy as AV18 in fitting NN scattering data. In addition to such two-body potentials, it was shown that one has to include also a three-body interaction to being able to describe nuclear systems accurately, e.g., to correctly describe the binding energy of light nuclei [28]. However, the three-body force is not as well constrained as the NN interaction. The Urbana IX (UIX) three-body force has been originally proposed to be used in combination with the Argonne AV18 and AV8’ [36] potentials. Although it slightly underbinds the energy of light nuclei, it has been extensively used to study the equation of state of nuclear and neutron matter [21, 30, 37]. The AV8’+UIX Hamiltonian has been used in many works, and provided a fairly good description of neutron star properties [37].

However, such phenomenological interactions suffer from certain shortcomings. Most importantly, they do not enable to estimate reliable theoretical uncertainties and cannot be systematically improved. These shortcomings can be addressed with the advent of chiral effective field theory (EFT), which offers a systematic expansion of nuclear forces that allows for theoretical uncertainties [38, 39].

In this approach, the relevant degrees of freedom are nucleons that can interact via explicit pion exchanges or via short-range contact interactions. The relevant diagrams entering in the nucleon-nucleon interaction are systematically organized in powers of $p/\Lambda_{b}$, where $p$ is the typical momenta of nucleons in the given nuclear system, i.e. similar to the pion mass $m_{\pi}\approx 140$ MeV, and $\Lambda_{b}\approx 500-600$ MeV [40] is the so-called breakdown scale, where the chiral EFT expansion is expected to loose its validity due to the increasing importance of shorter-range physics, i.e., new degrees of freedom. Chiral interactions include long-range pion-exchange physics explicitly, while short-range physics is described by a general operator basis consistent with all symmetries of the fundamental theory, Quantum Chromodynamics. The long-range pion-exchange terms entering into the chiral EFT potentials are fully determined by $\pi$-nucleon scattering data, while the parameters associated to the contact terms, called low-energy constants (LECs), are typically constrained by fitting nucleon-nucleon scattering data. For more details see Ref. [38]. Hamiltonians from chiral EFT have been recently combined with QMC methods [41, 42, 43, 44, 31], and give very reasonable predictions of properties of nuclei up to A=16, including energies, radii, and momentum distributions, and neutron-$\alpha$ scattering.

In addition, these interactions give a reasonable description of PNM [43, 45]. We show the PNM EOS calculated using AFDMC with local chiral interactions and the AV8’+UIX interaction in Fig. 4. The red and gray bands represent AFDMC results using three different parametrizations of the three-body force constrained in light nuclei, although they give different results in pure neutron matter due to regulator artifacts; see Ref. [43] for more details. In particular, one of them gives negative pressure in pure neutron matter between $1-2\rho_{0}$, and is represented by the gray band in Fig. 4. The other two parametrization give instead the red band in the figure, and can be employed to describe neutron stars [45].

The AFDMC results are compared to other models that are commonly used in astrophysical simulations: the Lattimer-Swesty EOS with incompressibility 220 [46], the TM1, SFHo, and SFHx EOSs [47], the FSU and NL3 EOSs [48], and the DD2 EOS [49]. We find that AFDMC calculations put strong constraints on the EOS of pure neutron matter.

The AFDMC results for the EOS of PNM can be conveniently parametrized using the functional

[TABLE]

that represents the energy per particle (without the rest mass energy) as a function of neutron density. We can define the symmetry energy as

[TABLE]

where $E_{\rm SNM}(\rho)$ is the EOS of symmetric nuclear matter, which is not known. However, by requiring that symmetric nuclear matter saturates at an energy $E_{\rm SNM}(\rho_{0})=-16$ MeV, and using the expression for $E_{\rm PNM}(\rho)$ above, the symmetry energy and slope can be obtained by

[TABLE]

Thus, for $E_{\rm SNM}(\rho)$ we find

[TABLE]

Assuming a quadratic expansion in the proton fraction $x=\rho_{p}/\rho$, the EOS at a finite proton fraction is given by

[TABLE]

The latter equation gives the expected results for pure neutron matter ($x=0$) and symmetric nuclear matter ($x=1/2$).

A very common parametrization for the symmetry energy is given by

[TABLE]

and using this definition, the slope parameter $L$ is linear with $C=E_{sym}(\rho_{0})$:

[TABLE]

From Eq. (6), we can obtain $C$:

[TABLE]

We can put an additional constraint on our simple ansatz by requiring that the pressure of SNM is zero at saturation, i.e.,

[TABLE]

This leads to the condition

[TABLE]

Note that by combining equations, we also find

[TABLE]

With these simple assumptions, the general form of the EOS as a function of density and proton fraction becomes:

[TABLE]

The parametrizations for selected EOSs are reported in Table 1, together with the corresponding symmetry energy and its slope.

While $L$ is now strongly constrained by chiral EFT and QMC, these constraints complement those from X-ray observations of neutron stars. Ref. [50] found $43<L<52$ MeV to within 68% confidence based on an analysis of neutron-star observations using Eq. 4.

3 Neutron star properties

As described above, once the EOS is specified it is easy to calculate the mass-radius relation of a neutron star. However, since the neutron star does not consist of pure neutron matter, one has to find ways of extending the microscopic PNM calculations to neutron-star conditions. In the following, we will discuss a few possibilities.

3.1 Effect of leptons

Typical neutron-star properties can be calculated directly from the PNM EOS, but such an approach misses the effects of the neutron-star crust and the remaining protons. Therefore, a more realistic EOS should contain these effects which can be estimated from the PNM EOS.

By starting from Eq. 15 we can solve for $x(\rho)$ by imposing $\beta$-equilibrium between neutrons, protons, electrons, and muons. With

[TABLE]

where $\epsilon$ is the total energy density, $m_{n}$ is the neutron mass, $m_{p}$ is the proton mass, and $\mu_{y}$ the chemical potential with $y=n,p$. We can easily obtain

[TABLE]

Charge neutrality requires that

[TABLE]

We take electrons to be relativistic and degenerate:

[TABLE]

and for the muons

[TABLE]

We then calculate all the fractions by imposing the charge neutrality, chemical potential as above, and

[TABLE]

Homogeneous matter in $\beta$-equilibrium is a valid model for sufficiently high densities, where nuclei are not present. Therefore, the AFDMC EOS is used for $\rho\geq\rho_{\rm crust}=0.08$ fm*-3*. For the low-density EOS, describing the crust of a neutron star, results of earlier works can be used; see e.g. Ref. [51] and [52]. For the AV8’+UIX EOS, we show the MR relations for the PNM EOS and when we assume $\beta$ equilibrium in Fig 5, where also the proton and lepton fractions are presented. We can see that the effect of protons is giving a small correction to the neutron-star radius given the current uncertainties in neutron-star observations [50].

3.2 High-density extrapolations

The approach in the previous section assumes a description in terms of nucleons and leptons to be valid in the whole neutron star. However, while this assumption might be true, the EOS might also explore more extreme density behavior at higher densities, as produced by, e.g., strong phase transitions to exotic forms of matter.

To being able to analyze neutron-star properties, general extrapolations schemes have to be used which are constrained by nuclear-physics input at low densities as well as observational constraints. Such general extrapolation schemes can be based on piecewise polytropes, where polytropes are line segments in $(\log\epsilon,\log P)$ space; see e.g., Refs. [53, 7, 54]. However, Ref. [55] showed that line segments in $(\epsilon,P)$ space can more easily represent models which have a phase transition to exotic matter in the neutron-star core.

In contrast, in Refs. [56, 57] another general extrapolation scheme starting from the PNM EOS from chiral EFT was developed, that used the speed of sound, $c_{S}$ in neutron stars. This scheme was based on the initial work of Refs. [58, 59], but represents an extension of these models by exploring all allowed parameter space for the speed of sound $c_{S}$, defined as

[TABLE]

In particular, models are constrained by the PNM EOS up to a certain density $n_{\rm tr}$ which is varied between 1-2 $\rho_{0}$. This PNM EOS is extended to $\beta$-equilibrium and includes a crust as discussed in Ref. [60]. From the resulting neutron-star EOS, the speed of sound is computed up to $n_{\rm tr}$. Beyond this density, many possible paths in the $c_{S}-n$ plane are explored by randomly sampling several points $c_{S}^{2}(n)$ between $n_{\rm tr}$ and 12 $\rho_{0}$, and connecting them by linear segments. During this procedure, it is enforced that $0\leq c_{S}\leq c$ and that the resulting EOSs are sufficiently stiff to support a two-solar-mass neutron star [3, 4]. For more details on this extrapolation scheme, see Ref. [61].

We show the resulting MR regions in Fig. 6, where we compare with the AV8’+UIX result from the previous section as well as the same model EOS of Fig. 4. We find that chiral EFT input can place strong constraints on the MR relation, ruling out too stiff model EOS, e.g., the NL3 parametrization. In particular, the density range between 1-2 $\rho_{0}$ is very important to reduce the uncertainty in the MR plane and improved calculations in this density range with smaller uncertainties will be useful to pin down the MR relation of neutron stars.

4 Connecting the Microphysics with the Multi-messenger Observations

Multi-messenger astronomy requires a strong foundation of microphysics in order to fully interpret the observations. The first evidence of this fact was the observation of both photons and neutrinos from supernova 1987A. This observation confirmed the basic picture of stellar evolution, which required significant nuclear and neutrino physics input[e.g. 62, 63, 64, 65, 66, 67]. More recently, the first gravitational-wave observation of a binary neutron-star merger, GW170817, has further highlighted the tight connection between microphysics, including the nuclear EOS, and multi-messenger astronomy observations [e.g. 11, 68].

Neutron-star mergers have long been thought to be the progenitors of short gamma-ray bursts (GRBs) and a significant, if not dominant, site of r-process nucleosynthesis [e.g. 69, 70, 71, 72, 73, 74, 75, 76, 77]. This basic picture was confirmed by GW170817, which was observed in coincidence with a short GRB (GRB170817A) from the same location in the sky [e.g. 78, 79, 80]. Subsequent followup observations with instruments and telescopes spanning the entire electro-magnetic spectrum revealed a kilonova [e.g. 78, 81, 82, 83, 84, 85]. Kilonovae are rapidly fading (timescale of a few weeks) optical and infrared transients that are powered by the radioactive decay of the newly synthesized r-process elements [e.g. 86, 87, 88, 89, 90]. They are the smoking gun signatures of r-process nucleosynthesis events.

Even though GW170817 confirmed some aspects of the basic picture of neutron-star mergers and the role they play in r-process nucleosynthesis, many details remain unresolved. Over the coming years, LIGO/VIRGO will observe more neutron-star mergers and we should be able to find the associated kilonovae for at least some of them [e.g. 93, 94, 95, 96, 97]. A detailed understanding of the microphysics will be crucial in accurately modelling the rich multi-messenger observational data from these events to fully understand them and extract system parameters that are not directly observable. Specifically, microphysics and the nuclear EOS directly impact the amount and possibly morphology of the mass ejecta and accretion disk [e.g. 98, 99, 100, 101, 102, 103, 104, 105, 68, 106, 107] that forms around the compact central remnant. They also determine how long a hypermassive neutron star (HMNS) could live, if one forms, and this lifetime affects the amount of ejecta blown away from the disk and its neutron-richness, which in turn determines what elements the r-process can synthesize [e.g. 108, 92]. An example of the effect of the lifetime of the HMNS is shown in Fig. 7. Neutrino interactions in general are very important to determine the composition of the ejecta and subsequent nucleosynthesis [e.g. 109, 110, 111, 112, 113, 114, 115, 116]. All of these properties affected by microphysics directly influence the kilonova lightcurve, nucleosynthetic yields, and possibly the gravitational wave signal. Therefore, we need to get the microphysics right in order to draw meaningful conclusions from neutron-star merger and kilonova observations about their intrinsic properties and how they enrich the galaxy with heavy elements [e.g. 117].

4.1 Constraining the Microphysics of Merger Simulations

with GW170817

Neutron-star merger simulations require a three-dimensional EOS table: a description of several thermodynamic quantities as a function of the baryon density, $\rho$, electron fraction, $Y_{e}$, and the temperature, $T$. Until recently, only about a dozen of these EOS tables was available. These EOS tables explore only a tiny fraction of the large space of EOSs which appear to be compatible with our current knowledge of the nucleon-nucleon interaction. For a review, see Ref. [118].

This state of affairs changed when Ref. [119] released an open-source code for EOS tables built upon the Skyrme interaction. This code allows one to fully explore a large space of EOSs. In particular, a merger simulation may systematically probe the sensitivity of the observables generated by the merger simulation to the parameters in the Skyrme interaction. This code, however, is limited by the applicability of the Skyrme interaction.

In Ref. [120] we presented a new class of phenomenological EOS for homogeneous nucleonic matter. The EOS is constructed to simultaneously match (i) second order virial expansion coefficients from nucleon-nucleon scattering phase shifts at high temperature and low density, (ii) experimental results of nuclear mass and radii, (iii) QMC calculations for neutron matter at saturation density, (iv) astrophysical observational measurements on neutrons star radii, (v) theoretical calculations with chiral perturbation field theory at finite temperature near the saturation density. It allows for computing the variation in the thermodynamic quantities based on the uncertainties of nuclear interaction.

Our free energy per particle can be written as

[TABLE]

where $x$ is the number of protons per baryon (assumed here to be equal to $Y_{e}$), $f_{\rm virial}$ is the virial free energy contribution [121] and $f_{\rm deg}$ is the free energy for degenerate matter. The function is defined by

[TABLE]

where $z_{i}\equiv\exp(\mu_{i}/T)$ are the fugacities. Since the virial expansion is valid when $z_{n}^{2},z_{p}^{2}\ll 1$, this functional form gives $g\approx 1$, thus the dominant contribution is from $F_{\rm virial}$. Otherwise, if $z_{n}$ or $z_{p}$ are large, $F_{\rm np}\approx F_{\rm deg}$. The free energy per particle of degenerate matter is further defined assuming quadratic expansion

[TABLE]

with $\delta=1-2x$. The 1000 parameter set of Skyrme model was chosen from UNEDF collaboration fitted to several nuclear mass, charge radii and pairing energies using Bayesian inference [122]. The $F_{\rm hot}$ is finite temperature results based on Kohn-Luttinger-Ward pertubation series fitted to Skyrme functional form [123, 124]. The symmetry energy is defined by

[TABLE]

where we interpolate between Eq. (4) near saturation density $\rho_{0}$ and a polynomial fit to neutron star observational data [7] above $2\rho_{0}$ using a function $h$, defined as

[TABLE]

where $\gamma=20.0\,{\rm fm}^{3}$. At zero temperature, the EOS can be compared to Eq. (15).

The advantage of the formalism above is that it allows us not only to compute the EOS of homogeneous nucleonic matter over the full range of densities, electron fractions, and temperatures, but it allows us to describe the probability distribution of the EOS. Formulating the EOS in this way allows us to easily determine the impact that observations might have on the EOS. This is demonstrated in Fig. 8 where the correlation between the free energy per baryon at $\rho=3\rho_{0}$ and the tidal deformability of a 1.4 $\mathrm{M}_{\odot}$ neutron star is shown. The tidal polarizability describes how a neutron star deforms under an external gravitational field, as produced by a companion star. It is given by

[TABLE]

with the tidal Love number $k_{2}$ which has to be solved together with the TOV equations; see, e.g., Ref. [125]. As is expected, the free energy of neutron-rich matter is strongly corrlated with the tidal deformability, but the uncertainty in the free energy of isospin-symmetric matter is not impacted by constraints on the tidal deformability.

4.2 EOS constraints from tidal polarizabilities

The analysis of the gravitational-wave signal from the inspiral phase of two merging neutron stars allows to place constraints on properties of the binary system, e.g., the chirp mass or mass ratio, as well as on properties of individual neutron stars. For a neutron-star binary, the binary tidal polarizability can be defined as

[TABLE]

This quantity, $\tilde{\Lambda}$, can be constrained from the GW signal of a NS merger. For the first NS merger observed, GW170817, this tidal polarizability originally was constrained to be $\tilde{\Lambda}\leq 800$ at a 90% confidence level [126] but later modified after several reanalyses [127, 128] to values around $70\leq\tilde{\Lambda}\leq 700$.

Qualitatively, neutron stars with larger radii can exhibit stronger deformations than small NS and, therefore, also have larger tidal polarizabilities. Hence, the tidal polarizability is closely related to the structure of a NS and can help to place constraints to the EOS. Using the chiral EFT results discussed in Sec. 2 up to saturation density (left panel of Fig. 6), in Fig. 9 we show the resulting MR and EOS envelopes if $\tilde{\Lambda}$ would be constrained to the 4 specified ranges, stretching from $0\leq\tilde{\Lambda}\leq 800$. While the observation GW170817 basically allows for a total radius range of $8.4-13.6$ km for a typical $1.4M_{\odot}$ neutron star, we see that more precise observations in future might help to drastically constrain the radii of neutron stars and, therefore, the EOS. For example, constraining the tidal polarizability to lie in an interval with a width $\Delta\tilde{\Lambda}\approx 200$, reduces the radius uncertainty to $\approx 2$ km for a typical neutron star.

4.3 Bayesian Inference: Using Neutron Star Observations to

Determine the EOS

A method to directly determine the EOS from neutron-star mass and radius observations was first described in Ref. [129]. Shortly thereafter, Ref. [7] presented an alternate method using Bayesian inference. Ref. [7] also obtained the first quantitative results, obtaining uncertainties in the EOS around 30% just above the nuclear saturation density. Bayesian inference is still the tool of choice for extracting the EOS, although many works employ simple Monte Carlo methods which are identical to Bayesian results with trivial priors.

One critical question (which may be soon answered) is whether or not the EOS contains a strong phase transition: a region in the EOS where the derivative of the pressure with respect to the energy density (identical to the speed of sound) is nearly zero. Because the pressure must be continuous in the star, a strong phase transition creates a thin region in the neutron star where the energy density increases strongly with decreasing radius. Strong phase transitions of this type can result in mass-radius curves which have multiple branches, as described in Ref. [130] and recently explored in Ref. [131].

Below, we show that multiple branches in the mass-radius curve cannot be easily handled in the Bayesian methods which have been previously used in the literature and we show how this difficulty can be fixed. The traditional method for extracting the mass-radius curve from $N$ neutron star mass and radius determinations (each expressed in the form of a two-dimensional probability distribution ${\cal D}(R,M)$) is to use a likelihood function of the form

[TABLE]

where $\{p_{j}\}$ are the parameters of the EOS and the TOV equations are solved to obtain the function $R(M_{i},\{p_{j}\})$. (This likelihood is analogous to Eq. (15) of Ref. [129] and Eqs. (29) and (31) of Ref. [7], but the formalism in these two works is a bit different, see Ref. [132] for an interesting discussion comparing the two works.) The critical feature of these likelihoods is that they are both expressed as integrals over the gravitational mass. It is thus clear that these likelihood functions are identically zero for horizontal mass-radius curves.

In writing Eq. (30) we specifically made a choice to parametrize the M-R curve with the gravitational mass rather than the radius. Alternatively, we could have chosen the radius instead

[TABLE]

where the TOV equations now provide the function $M(R_{i},\{p_{j}\})$. This likelihood is not equivalent to Eq. (30); it results in different posterior distributions and, given an EOS parametrization, will be maximized at a different location in parameter space. To see this, note that Eq. (30) can be rewritten as

[TABLE]

which is not equal Eq. (31).

Choosing either Eq. (30) or Eq. (31) is further problematic because M-R curves like those from EOSs with strong phase transitions can result in two (or even three) configurations with the same gravitational mass yet different radii. The only way to use the likelihood in Eq. (30) is to break up the integral into different parts for each branch of the M-R curve

[TABLE]

where $M_{0}<M_{1}$ and $M_{2}<M_{3}$ but $M_{2}<M_{1}$. This form of the likelihood makes the numerical evaluation of the likelihood more awkward. One possible resolution is to use the central neutron star pressure, $P$ instead of the gravitational mass $M$ or radius $R$, i.e.

[TABLE]

This integral is simpler because the functions $R(P_{i},\{p_{j}\})$ and $M(P_{i},\{p_{j}\})$ are never multiply-valued and the pressure is continuous within the star. This form has the additional advantage that neither “horizontal” nor “vertical” mass-radius curves are arbitrarily given zero likelihood, though this choice is still part of the prior distribution. Note that choosing the central energy density, instead of the central pressure, as suggested by Eq. (14) in Ref. [8] also leads to multiple integrals as in Eq. (33) when strong phase transitions create regions where $dP/d\varepsilon=0$.

Ref. [133] recently clarified why the ambiguity between Eqs. (30) and (31) arises. It originates in the fact that we are attempting to match a one-dimensional model curve (the size of the EOS parameter space is not important for this discussion) to a two-dimensional data set ${\cal D}(R,M)$. In the language of differential geometry, we are embedding the curve into the data space and the ambiguity we found above is related to the choice of how we perform this embedding. The general form requires a metric ($g$) to specify how distances are defined

[TABLE]

where $j$ and $k$ range from 1 to 2 and $X_{1}\equiv M$ and $X_{2}\equiv R$ and for simplicity we can assume the metric is independent of $i$. The ambiguity in choosing the likelihood function is now explicit in the choice of the metric, $g_{jk}$. In the context of Bayesian inference, this metric choice is to be regarded as part of the prior probability. Choosing $\lambda_{i}=P_{i}$ and a trivial metric $g_{jk}=\delta_{jk}$ gives Eq. (34). It is clear, however, that other choices of the metric (the elements may depend on $\lambda$ or $X$) will give different results.

The implication, in the context of future gravitational wave observations of neutron star mergers, is that unambiguous determinations of the EOS will have to wait until sufficient data is obtained as to make the posterior distributions of quantities of interest (such as EOS parameters) are not strongly dependent on the choice of prior probability. This prior probability includes the EOS parametrization and the prior probabilities for the EOS parameters [18] as well as the prior choice of metric in Eq. (35).

5 Summary

In this paper, we have reviewed current calculations of the EOS of neutron-star matter starting from AFDMC calculations of pure neutron matter using interactions from chiral EFT. We have presented the EOS of PNM with theoretical uncertainties and explained, how to extend these PNM calculations to neutron-star conditions. We then used this to explore current theoretical uncertainties for the mass-radius relation of neutron stars.

We have then discussed how this microphysics can be connected to multimessenger observations of neutron-star mergers, using several EOS models or Bayesian inference.

Neutron-star mergers offer an ideal way to constrain the EOS of strongly interacting matter, and thus, nuclear interactions. To pin down the nuclear interaction, theorists working on a solid theoretical description of microphysics, computational astrophysicists who simulate neutron-star mergers and supernovae, and observers have to work hand-in-hand to reliably extract constraints from future merger observations.

Acknowledgments

The work of S.G. and I.T. was supported by the U.S. DOE under contract DE-AC52-06NA25396 and by the LANL LDRD program. The work of S.G. was also supported by the NUCLEI SciDAC program and by the DOE Early Career Research Program. the work of JL was supported by the Laboratory Directed Research and Development program of Los Alamos National Laboratory under project number 20190021DR. The work of XD and AWS was supported by DOE SciDAC grant DE-SC0018232. The work of MA and AWS was supported by NSF grant PHY 1554876. This work was supported by the US Department of Energy through the Los Alamos National Laboratory and used resources provided by the Los Alamos National Laboratory Institutional Computing Program. Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of U.S. Department of Energy (Contract No. 89233218CNA000001). We also used resources provided by NERSC, which is supported by the US DOE under Contract DE-AC02-05CH11231, and by the Jülich Supercomputing Center.

Figs. 2 and 3 are open source and available at https://github.com/awsteiner/nstar-plot.