Nuclear and neutron-star matter from local chiral interactions

We report a quantum Monte Carlo calculation of the equation of state of symmetric nuclear matter using local interactions derived from chiral effective field theory up to next-to-next-to-leading order fit to few-body observables only. The empirical saturation density and energy are well reproduced within statistical and systematic uncertainties. We have also derived the symmetry energy as a function of the density, finding good agreement with available experimentally derived constraints at saturation and twice saturation density. We find that the corresponding pressure is also in excellent agreement with recent constraints extracted from gravitational waves of the neutron-star merger GW170817 by the LIGO-Virgo detection.

Introduction. The nuclear equation of state (EOS) is of great interest for nuclear physics and nuclear astrophysics. At proton fractions x ∼ 0.5, for the so-called symmetric nuclear matter (SNM), the EOS sets the bulk properties of atomic nuclei and determines where atomic nuclei saturate. At lower proton fractions, x 0.1, the nuclear EOS determines the properties of neutron stars. The energy difference of nuclear matter at different proton fractions, e.g., between SNM and pure neutron matter (PNM) with x = 0, is governed by the nuclear symmetry energy S(n), with n the baryon density. The nuclear symmetry energy is a fundamental physical quantity that affects a range of neutron-star properties, such as cooling rates, the thickness of the neutron-star crust, the massradius relation [1][2][3], and the moment of inertia [4,5], and is deeply connected to properties of atomic nuclei, e.g., the dipole polarizability, the giant dipole resonance, and the neutron skin of neutron-rich nuclei [6]. Hence, it is possible to infer properties of the nuclear EOS at low densities and larger proton fractions from laboratory experiments with atomic nuclei, e.g., in the future Facility for Rare Isotope Beams (FRIB) at Michigan State University. The very neutron-rich nuclear matter that forms neutron stars, on the other hand, cannot be directly probed in any terrestrial experiment. Hence, neutron stars provide a unique and complementary laboratory for dense nuclear matter at extreme conditions.
Understanding the properties of the nuclear EOS has recently become even more critical with the advent of gravitational-wave (GW) astronomy and the first-ever detection of gravitational waves from a binary neutronstar merger by the LIGO-Virgo collaboration in 2017 [7,8]. This GW detection, together with the identification of such mergers as a site of r-process nucleosynthesis, has attracted tremendous interest. Many expected future observations during the current and future LIGO observing runs, like the first detection of a neutron-star black-hole merger event (S190814bv), will be of similar high impact. Reliable and precise calculations of the nuclear EOS in neutron-rich regimes and for sufficiently high densities are needed to accurately extract information from such astrophysical observations. In addition, reliable calcula-tions of the nuclear EOS at various proton fractions are needed in order to compare to experiments extracting the symmetry energy. Such calculations require a systematic theory for strong interactions that enables estimations of theoretical uncertainties, in combination with advanced many-body methods to solve the complicated nuclear many-body problem.
In this work, we address this timely, high-impact topic, and present quantum Monte Carlo (QMC) calculations of both neutron and symmetric nuclear matter using local chiral effective field theory (EFT) interactions [9,10], in order to constrain the nuclear EOS and the symmetry energy with robust uncertainty estimates. Our goal is to compare QMC results directly to observational/experimental constraints, rather than to other many-body calculations of nuclear matter [11][12][13][14][15][16][17], which will require a more detailed discussion on methods and interactions, outside the scope of this work.
Hamiltonian and method. We describe the nuclear many-body system as a collection of pointlike interacting nucleons with a nonrelativistic Hamiltonian H that includes two-body (v ij ) and three-body (V ijk ) potentials [18]. In this work we employ local interactions derived from chiral EFT up to next-to-next-to-leading order (N 2 LO) [9,10,19]. These give a reasonable description of nucleon-nucleon phase shifts within uncertainty estimates up to laboratory energies of 500 MeV, except in the triplet D waves. They include consistent two-and threebody potentials, and have been used to study the groundstate properties of nuclei [10,[20][21][22], few-neutron systems [23,24], and neutron-star matter [9,10,[25][26][27][28][29]. In particular, we consider the interactions with coordinatespace cutoff R 0 = 1.0 fm and two different parametrizations of the three-body contact term V E , namely E1 and Eτ . E1 employs the identity operator 1 between particle i and j, while Eτ employs the isospin operator structure τ i · τ j . The difference of the two parametrizations is due to regulator and cutoff artifacts. A comparison of both interactions allows one to quantify the uncertainty due to the regularization scheme and scale. We do not consider the cutoff R 0 = 1.2 fm because it leads to severe regulator artifacts in atomic nuclei [22]. For lo-arXiv:1912.09411v2 [nucl-th] 7 May 2020 cal chiral interactions, the short-range contact term V E , together with the three-body one-pion-exchange-contact term V D , are characterized by two low-energy couplings c D and c E which are fit to few-body observables, namely the α-particle binding energy and the spin-orbit splitting in the neutron-α P -wave phase shifts (see Refs. [10,22] for details). However, the interactions are capable of describing ground-state properties of nuclei up to (at least) A = 16 [21], and the E1 potential can simultaneously predict properties of neutron matter compatible with astrophysical observations [27].
The starting point of all QMC methods is the choice of a wave function representing the system, typically expressed as the trial state |Ψ T = F |Φ , where F is the correlation operator acting between pairs and triplets of particles. F incorporates strong spin and isospin dependence into the trial state, as induced by the employed nuclear Hamiltonian. For infinite matter, the term |Φ is built from a Slater determinant of plane waves φ k (i) = e ik·ri with momenta discretized in a finite box whose volume is determined by the chosen baryon density n and number of particles A [30]. The infinite system is then realized by applying periodic boundary conditions [31]. Finite-size corrections to the energy results are included as described in Refs. [30,32]. All the parameters of Ψ T are chosen by minimizing the variational energy as described in Ref. [33].
In this work, we make use of the auxiliary field diffusion Monte Carlo (AFDMC) method [18,22,34] that allows one to access the true ground state of a nuclear system by evolving the initial trial state in imaginary time according to the projection operator exp[−(H − E T )τ ]|Ψ T , where E T is a normalization parameter. In the limit of infinite imaginary time, higher-energy components in the trial state are filtered out, and the system is projected onto the ground state. Such imaginary-time evolution is performed by sampling both spatial coordinates and spin/isosospin configurations, the latter via a Hubbard-Stratonovich transformation. For local chiral interactions, the propagation of three-body terms is carried out via an effective Hamiltonian as described in Refs. [21,22]. The sign problem is initially suppressed by evolving the trial wave function using the constrained-path approximation [35]. An unconstrained evolution is then performed until the sign problem dominates and the variance of the results becomes severely large. Finally, expectation values are evaluated over the sampled configurations to compute the relevant observables [22]. Note that, similarly to the case of atomic nuclei, the approximate propagation of three-body terms and the unconstrained evolution for infinite matter are robust and under control. The associated uncertainties are included in the final Monte Carlo uncertainty estimate.
The AFDMC method has been used in the past to calculate the EOS of PNM using both phenomenological and local chiral interactions [2, 10, 25-29, 31, 36, 37], and attempts at calculating the EOS of asymmetric nuclear matter have been carried out for a simplified phenomenological model [30]. Here, we perform a study of the EOS of PNM and SNM using local chiral interactions up to N 2 LO. We consider, respectively, 66 neutrons and 28 nucleons in a periodic box described by the trial state |Ψ T , where the spin/isospin-dependent two-and three-body correlations are expressed as a sum of linear and quadratic spin/isospin operators as in Ref. [22]. Note that finite size effects in SNM are smaller than in PNM [30], hence, a smaller particle number is sufficient.
Due to the high computational cost of performing derivatives of the trial wave function, no spin-orbit correlations are typically included in the AFDMC wave function. However, as reported in Ref. [38], the largest contribution to the total energy given by spin-orbit terms can be obtained by using a simplified spin-orbit correlation that can be implemented in the AFDMC wave function by substituting the plane wave e ik·ri with where β is a variational parameter, and f LS ij is a spin-orbit radial function obtained from the solution of Schrödingerlike equations in the relative distance r ij as described in Ref. [18]. Equation (1) defines the so-called spinbackflow correlations, and is analogous to the implementation of standard backflow correlations in fermionic systems [39]. Spin-backflow correlations only imply spin rotations among the components of the Slater determinant |Φ [38], which makes them computationally cheap. As shown in a simplified case for PNM in Ref. [38], such correlations greatly improve the quality of the wave function, and their contribution to the total energy is not negligible.
Results. The AFDMC results for the EOS of PNM and SNM at N 2 LO are shown in the left panel of Fig. 1 for the E1 (red) and Eτ (blue) parametrizations. A table with the energy results as a function of the density is available in the Supplemental Material [43].
We have found that spin/isospin-dependent correlations yield a negligible improvement of the total energy in PNM, while, similarly to the case of atomic nuclei [22] Table I). Colored bands represent the uncertainties of the many-body calculations, which include both statistical Monte Carlo errors and the uncertainties coming from the truncation of the chiral expansion. In panel (a), the green box indicates the empirical saturation point [16]. In panel (b), we show experimental constraints on the symmetry energy below saturation density from Ref. [40], at saturation [41] and twice saturation density [for S(nsat) = 31 MeV] [42]. The dashed black curve is the Fermi gas result.
tions, of the order of < 2% at saturation density, but requires large computational resources.
The solid curves in Fig. 1 are fit to the AFDMC results according to the relations [44,45] E PNM (n) = a n n sat where n sat = 0.16 fm −3 is the empirical saturation density, n 0 and E 0 are saturation density and saturation energy for the given Hamiltonian, and K 0 , Q 0 , and Z 0 are the incompressibility, skewness, and kurtosis parameters. For SNM we fit the AFDMC energies above n = 0.12 fm −3 , since clustering is expected to appear at lower densities. Assuming that the system behaves as uniform matter over the whole density range, which is not a realistic picture for low density nuclear matter (hence the dashed curves in Fig. 1), we enforce E SNM (0) = 0 by adjusting Z 0 accordingly. All the fitting parameters, together with the empirical values, where available, are reported in Table I. In Fig. 1, colored bands represent the uncertainties of the many-body calculation, which include both statistical Monte Carlo errors and the uncertainties coming from the truncation of the chiral expansion. The latter is evaluated according to the prescription by Epelbaum et al. [46]. In this work we consider the average momentum  [16,44] are shown for comparison.
Par.   [20], k F being the Fermi momentum of PNM or SNM, and Λ b = 500 MeV [10]. The difference of the two bands indicates an additional source of uncertainty due to the regularization scheme. The EOS of SNM for the E1 potential saturates at a slightly higher density n = 0.22(1) fm −3 and higher energy E = −13.96(8) MeV compared to the empirical point, while the incompressibility K 0 = 223(16) MeV lies within the expected range [44]. The skewness parameter is poorly constrained, but is consistent, for instance, with the analysis carried out in Refs. [44,47], where the authors considered terms up to n 4 to fit the SNM EOS. Note that, by considering in Eq. (3) only terms up to n 3 , and by constraining the skewness parameter to impose the passage to zero, n 0 , E 0 , and K 0 stay within the previously identified ranges, while Q 0 changes to −145(108) MeV. This value is similar to that extracted from Skyrme parametrizations, where a similar EOS up to n 3 terms is used (see, for instance, Ref. [48]).
The new PNM EOS is consistent with earlier AFDMC results using local chiral interactions [10], where simplified wave functions were used and the unconstrained propagation was not performed. Hence, the PNM EOS for the E1 interaction remains stiff enough to be compatible with astrophysical observations, while the Eτ potential is too attractive at high density, resulting in negative pressure above 0.20 fm −3 . This behavior is similar to the SNM case, where saturation is reached at high density [n = 0.36(1) fm −3 ], with an energy below the empirical saturation value [−17.29(9) MeV]. The incompressibility [K 0 = 184(64) MeV] is within the expected range, but the uncertainty is large. The skewness parameter is also very poorly constrained.
The nuclear symmetry energy S and its slope L as a function of the density are defined as: Please note, that the symmetry energy as defined Eq. (4) is similar to the quadratic term in the isospin asymme- S vs L (red) and SPNM vs LPNM (purple) at the empirical saturation density nsat for the E1 potential compared to experimental constraints from nuclear masses [51], the neutron-skin thicknesses of Sn isotopes [52], the dipole polarizability of 208 Pb [53], giant dipole resonances (GDR) [54], isotope diffusion in heavy ion collisions (HIC) [55], and from isobaric analog states and isovector skin (IAS+∆R) [56]. The areas denoted by red-dashed lines are theoretical constraints from Ref. [1] (HS) and Ref. [2] (GCR), while the area denoted by black-dashed lines is the inference of Ref. [3] (SG). The thick black line shows the unitary-gas constraint from Ref. [57]. Figure adapted from Ref. [57]. try expansion if quartic contributions are small as expected [49,50]. The AFDMC results for the symmetry energy are reported in the right panel of Fig. 1 for both the local chiral interactions considered in this work. The E1 potential predicts values of the symmetry energy compatible with the available constraints at both saturation [41] and twice saturation density [42] (see also Table II). For the Eτ model, compatibility is only marginal at saturation density, and deteriorates at higher density due to the behavior of the corresponding PNM EOS above 0.20 fm −3 .
In Table II we report the values of the symmetry energy and its slope at n sat and 2n sat as obtained from Eqs. (4) and (5). For the E1 interaction, L is compatible with the empirical values, within the estimated uncertainties. The Eτ potential, instead, predicts a too low value for L, as a consequence of the corresponding too soft PNM EOS. If the empirical saturation was achieved for the employed interactions, n 0 ≡ n sat , the symmetry energy and its slope at n sat would be completely determined by the PNM EOS: S PNM = a + b − E sat , with E sat = −15.86 MeV, and L PNM = 3(aα + bβ), . Pressure of PNM (red band) and matter in βequilibrium (red hatched band) as a function of density for the E1 interaction. The green area is the pressure extracted from the GW signal GW170817 for EOSs reproducing a 2M neutron star [58]. see values in Table II. Differences between (S, L) and (S PNM , L PNM ), due to extra energy and pressure contributions from the SNM EOS, are clearly shown in Fig. 2, where L is plotted versus S at the empirical saturation density. Shaded areas are calculated by sampling thousands of curves within the uncertainty bands of PNM and SNM, calculating S and L from these samples, and plotting the resulting densities. Colored bands with labels are experimental constraints as in Ref. [57].
Finally, in Fig. 3, for the E1 interaction we show the pressure as a function of the density for PNM (red solid curve) and β-equilibrated matter (red dashed curve), where the latter is obtained consistently from our results for PNM and SNM, including the uncertainty bands. The green area is the pressure extracted by the LIGO-Virgo collaboration from the GW signal GW170817 [58]. At low densities, the LIGO extraction is stitched to the SLy EOS, which is why the uncertainty band decreases. We find that our calculations lead to pressures that are compatible with but at the lower bound of the LIGO extraction.
Summary. We have performed QMC calculations of symmetric nuclear matter and the symmetry energy using realistic nuclear interactions from chiral EFT. For our local chiral interactions at N 2 LO, we find saturation at ≈ 1.4n sat and ≈ −14 MeV, but our results overlap with the empirical saturation point within uncertainties. Our results for the symmetry energy and its density behavior agree with previous inferences from experimental data up to 2n sat , and lead to a pressure for β-equilibrated matter in agreement with inferences from the neutron-star merger GW170817. We stress that the interactions employed here have been fit only to few-body systems, i.e., nucleon-nucleon scattering, and A = 4 and 5 nuclei, and it is quite remarkable that the EOS we have calculated is reproducing well the constraints from GWs. Table III. Equation of state of PNM and SNM at N 2 LO for coordinate-space cutoff R0 = 1.0 fm (see Fig. 1 of the main text). Results for the two operator structures E1 and Eτ are shown. First error is the statistical Monte Carlo uncertainty. Second error is the uncertainty coming from the truncation of the chiral expansion according to the prescription of Epelbaum et al. [46] (see main text for details). Densities are in fm −3 , energies in MeV/A.