Thermodynamics of isospin-asymmetric nuclear matter from chiral effective field theory

Thermodynamics of isospin-asymmetric nuclear matter from chiral effective field theory

Abstract

The density and temperature dependence of the nuclear symmetry free energy is investigated using microscopic two- and three-body nuclear potentials constructed from chiral effective field theory. The nuclear force models and many-body methods are benchmarked to properties of isospin-symmetric nuclear matter in the vicinity of the saturation density as well as the virial expansion of the neutron matter equation of state at low fugacities. The free energy per particle of isospin-asymmetric nuclear matter is calculated assuming a quadratic dependence of the interaction contributions on the isospin asymmetry. The spinodal instability at subnuclear densities is examined in detail.

I Introduction

Determining the thermodynamic equation of state (EoS) of nuclear matter is a central objective in modern nuclear theory. The isospin-asymmetry dependence of the EoS is essential for many phenomena in nuclear physics and astrophysics Steiner et al. (2005). Next-generation radioactive beam facilities Balantekin et al. (2014); Li et al. (2008) studying the reactions and structure of exotic neutron-rich isotopes in particular provide motivation to improve our microscopic description of highly isospin-asymmetric nuclear matter. To a certain extent, the isospin-asymmetry dependence of the EoS is described by the so-called symmetry free energy. In this work, starting from microscopic calculations of the EoS of isospin-symmetric nuclear matter and pure neutron matter using two- and three-body chiral nuclear interactions, we examine in detail the density and temperature dependence of the symmetry free energy. Furthermore, we construct the EoS of isospin-asymmetric nuclear matter, and study the behavior of the nuclear liquid-gas instability as the proton fraction is decreased. The present work is a first step toward the development of a chiral effective field theory thermodynamic equation of state across the temperatures, densities and isospin asymmetries relevant for describing astrophysical phenomena and the matter produced experimentally in heavy-ion collisions at moderate energies.

Chiral effective field theory (EFT) provides the basis for the study of strongly-interacting matter at the energy scales characteristic of normal nuclei Machleidt and Entem (2011); Holt et al. (2013); Epelbaum et al. (2009). In EFT, microscopic nuclear interactions are organized in a systematic expansion, with many-nucleon forces naturally included. The low-energy constants parametrizing the interactions are generally fixed by high-precision fits to nucleon-nucleon scattering phase shifts and properties of light nuclei. Employing chiral interactions in calculations of nuclear many-body systems then gives pure predictions without additional fine tuning, and theoretical uncertainties can be estimated Epelbaum et al. (2009); Bogner et al. (2010); Tews et al. (2013); Furnstahl et al. (2015); Sammarruca et al. (2015) by varying the resolution scale, the fitting procedures applied to fix the low-energy constants, and the chiral order of the nuclear potentials. In the case of isospin-symmetric nuclear matter, semiempirical constraints from the zero-temperature saturation energy, density and incompressibility as well as the critical point of the nuclear liquid-gas phase transition have been reproduced with low-momentum chiral nuclear forces in many-body perturbation theory Coraggio et al. (2014); Wellenhofer et al. (2014). This motivates a study of isospin-asymmetric nuclear matter, which is by comparison much less constrained by experimental data.

The symmetry free energy is defined as the difference between the free energy per particle in homogeneous isospin-symmetric matter (SNM) and pure neutron matter (PNM):

(1)

where is the temperature, is the total nucleon density, and is the isospin-asymmetry parameter (with the neutron/proton density). The free energy per particle of homogeneous nuclear matter with proton fraction can be written as

(2)

where . If isospin-symmetry breaking effects are neglected is an even function of . It has been validated in various microscopic many-body calculations (see e.g., Refs. Bombaci and Lombardo (1991); Drischler et al. (2014)) that the isospin-asymmetry dependence of at zero temperature is approximately quadratic to high accuracy over the entire range , i.e., . At finite temperatures however the free Fermi gas contribution to contains large terms with quartic and higher powers of , which we quantify explicitly in this work. By comparison, in initial calculations we have found that at the densities and temperatures relevant for the nuclear liquid-gas phase transition the interaction contributions give rise to weaker nonquadratic terms and will therefore be assumed to have a quadratic dependence on in this work. Future research will address the accuracy of this approximation in greater detail.

The liquid-gas phase transition in isospin-asymmetric nuclear matter (ANM) involves isospin distillation: in the transition region the system separates into two phases whose proton concentrations deviate from the global value , with and for the case . These distillation effects are a generic property of first-order phase transitions in binary thermodynamic systems. If isospin-symmetry breaking effects are neglected, neutrons and protons are thermodynamically indistinguishable in SNM. Hence, thermodynamically SNM is a pure substance, and there is no isospin distillation for . The new features of the phase transition in ANM were discussed in detail in Refs. Ducoin et al. (2006); Hempel et al. (2013); Müller and Serot (1995); Chomaz et al. (2004); Margueron and Chomaz (2003) using different phenomenological models of the nuclear force. The transition region is comprised of regions of metastable and unstable single-phase equilibrium, corresponding to different dynamical phase separation mechanisms: nucleation and spinodal decomposition Binder (1987); Abraham (1979). In this work we focus mostly on the spinodal which delineates the inner region of thermodynamic instability where no metastable state can exist. The evolution of the unstable spinodal region with increasing isospin asymmetry is analyzed in terms of the trajectory of the critical temperature . Moreover, we determine the neutron drip point in cold nuclear matter and the fragmentation temperature above which no self-bound drop of liquid nuclear matter can exist. It should be emphasized that in the present paper we discuss the liquid-gas instability of infinite nuclear matter, i.e., bulk nucleonic matter without Coulomb interactions. The inclusion of surface energies and the Coulomb repulsion of protons is required for an accurate description of the matter produced in intermediate-energy heavy-ion collisions. In neutron stars and core-collapse supernovae the realization of the nuclear liquid-gas instability is strongly affected by the presence of a (highly incompressible) charge neutralizing background of electrons (and myons) Ducoin et al. (2007); Napolitani et al. (2007). In particular, the competition between nuclear and Coulomb interactions (frustration) entails the fomration of mesoscopic inhomogeneities with nontrivial spatial structures. These so-called pasta phases have been studied extensively in the literature Ravenhall et al. (1983); Hashimoto et al. (1984); Pethick and Potekhin (1998); Horowitz et al. (2004); Maruyama et al. (2005); Watanabe et al. (2005); Avancini et al. (2012); Schneider et al. (2013). Including these effects as well as the presence of few-nucleon bound-states Horowitz and Schwenk (2006a); Shen et al. (2010); Typel et al. (2010) at very low densities represents a future challenge.

The paper is organized as follows. In Sec. II we recall the main results for the EoS of SNM obtained in Ref. Wellenhofer et al. (2014), and show results for additional derived thermodynamic quantities, i.e., the entropy per nucleon and the internal energy per nucleon. In Sec. III we extend the calculations to pure neutron matter. The zero-temperature results are compared to those from recent quantum Monte Carlo simulations while the finite-temperature EoS at low densities is compared to the virial expansion. In Sec. IV we investigate the temperature and density dependence of the symmetry free energy, entropy, and internal energy. The thermodynamics of isospin-asymmetric nuclear matter is studied in Sec. V. In particular, we examine in detail the dependence on isospin asymmetry of the EoS of a free nucleon gas. Finally, Sec. VI provides a short summary.

Ii Isospin-symmetric nuclear matter

In Ref. Wellenhofer et al. (2014) we calculated the free energy per nucleon in infinite homogeneous SNM using the Kohn-Luttinger-Ward Kohn and Luttinger (1960); Luttinger and Ward (1960) many-body perturbation series including contributions up to second order, i.e.,

(3)

Here, counts the number of interaction insertions, corresponds to a nonrelativistic free nucleon gas, and is a correction term which together with reproduces the properties of a relativistic free nucleon gas over a wide range of densities and temperatures Fritsch et al. (2002). The first- and second-order terms and receive contributions from both the two-body and the three-body nuclear force. The second-order term includes (temperature and density dependent) self-energy corrections, and has been evaluated by approximating the three-nucleon interaction with a temperature and density dependent effective two-body potential (for details see Refs. Wellenhofer et al. (2014); Holt et al. (2010); Hebeler and Schwenk (2010); Hebeler et al. (2011); Holt et al. (2009). Explicit formulas for the different contributions in Eq. (3) are given in Ref. Wellenhofer et al. (2014). The effective one-body chemical potential is in one-to-one correspondence with the nucleon density via

(4)

where is the isospin projection quantum number and is the average nucleon mass. For Eq. (3) to be sufficiently converged at second order in , low-momentum interactions have to be used, i.e., interactions with restricted resolution in coordinate space (corresponding to an ultraviolet cutoff in momentum space). The various sets of N3LO (i.e., fourth order in the chiral expansion) two-body and N2LO three-body chiral low-momentum interactions used in Ref. Wellenhofer et al. (2014) correspond to different regularization methods, resolution scales , and low-energy constants. For interactions constructed at resolution scales appropriate perturbative behavior was found. The SNM equation of state obtained from the sets of two- and three-body potentials denoted by n3lo414 () and n3lo450 (), respectively, (see Refs. Coraggio et al. (2007, 2014, 2013) for details) agree with empirical constraints from the zero-temperature saturation energy, density and incompressibility Brockmann and Machleidt (1990); Blaizot (1980); Khoa et al. (1997); Youngblood et al. (1999), and with estimates for the critical point of the nuclear liquid-gas phase transition obtained through the analysis of data from multifragmentation, fission and compound nuclear decay experiments Karnaukhov et al. (2008); Natowitz et al. (2002); Ogul and Botvina (2002); Elliott et al. (2013). The values of these quantities obtained from n3lo414 and n3lo450 in Ref. Wellenhofer et al. (2014) are displayed in Table 1. 1

Figure 1: (Color online) Results for the free energy per nucleon , the pressure , the entropy per nucleon and the internal energy per nucleon in isospin-symmetric nuclear matter. The uncertainty bands correspond to calculations using two different sets of chiral low-momentum two- and three-body interactions, n3lo414 (solid lines) and n3lo450 (dash-dot lines). The unstable spinodal region is marked out explicitly. The critical point is shown as a circle (full circle for n3lo414, open circle for n3lo450). The zero-temperature endpoint of the low-density part of the spinodal is located at .

From the free energy per nucleon the pressure and the entropy per nucleon follow via standard thermodynamic relations:

(5)

The internal energy per nucleon is given by . The results for these quantities are shown in Fig. 1 for temperatures in the range . The spinodal region2 where the homogeneous (i.e., single-phase constrained) system is unstable with respect to infinitesimal density fluctuations is shown explicitly. Note that for low temperatures the region of negative pressure extends into the metastable region (cf. also Fig. 14), which is a generic property of liquids that are self-bound at low temperatures and a well-known feature of superheated molecular liquids Debenedetti (1996).

                                                     
n3lo414             -15.79         0.171         223          17.4          0.066         0.33
n3lo450             -15.50         0.161         244          17.2          0.064         0.32
Table 1: Zero-temperature saturation energy , density and incompressibility , as well as the critical temperature , density and pressure of the liquid-gas phase transition in isospin-symmetric nuclear matter from the sets of chiral two- and three-body nuclear interactions n3lo414 and n3lo450.

Iii Pure neutron matter

The free energy per particle in PNM, , is obtained by restricting the isospin sum(s) in Eq. (4) and in the different contributions in Eq. (3). Moreover, in PNM the three-body contributions proportional to the low-energy constants and are absent Holt et al. (2010). The results for the free energy per particle, the pressure, the entropy per particle and the internal energy per particle in PNM are shown in Fig. 2. Note that the uncertainty bars obtained by varying the resolution scale (n3lo414 vs. n3lo450) increase with temperature and are significantly reduced as compared to the SNM results in Ref. Wellenhofer et al. (2014), which is due to the decreased magnitude of nuclear interactions in PNM. At very low temperatures the internal energy per particle increases monotonically with increasing density (as required by the absence of a liquid-gas instability in PNM), but otherwise there is a local minimum at finite density.

Figure 2: (Color online) Results for the free energy per particle , the pressure , the entropy per particle and the internal energy per particle in pure neutron matter. The solid lines show the results from n3lo414, the dash-dot lines the n3lo450 results. The thick dashed lines correspond to the model-independent virial equation of state (VEoS) determined from neutron-neutron scattering phase shifts. The VEoS lines end where the fugacity is .

At very low energies, where higher partial waves are unimportant, the interaction between neutrons is characterized by the neutron-neutron scattering length . In the regime where is large compared to the interparticle separation, (with the Fermi momentum), a perturbative approach to neutron matter is not reliable. The model-independent virial equation of state (VEoS) computed by Horowitz and Schwenk in Ref. Horowitz and Schwenk (2006b) from neutron-neutron scattering phase-shifts provides a benchmark for perturbative calculations of low-density neutron matter at nonzero temperature. In the virial expansion, the grand canonical expressions for the pressure and the density are expanded in powers of the fugacity , leading to

(6)

where is the chemical potential, and is the neutron thermal wavelength. The second virial coefficient is given by

(7)

where is the sum of the isospin-triplet elastic scattering phase shifts at laboratory energy . From the pressure and density as functions of the fugacity the free energy per particle , entropy per particle , and internal energy per particle follow again from standard thermodynamic relations (see Ref. Horowitz and Schwenk (2006b) for details). The results for these quantities are shown as green dashed lines in Fig. 2.3 One sees that in the case of , and there are almost no visible deviations between the VEoS and the perturbative results. This seemingly perfect agreement is however misleading, because the discrepancies corresponding to the different treatment of the interactions in the virial and the perturbative approach are overpowered by the large size of the (nonrelativistic) free Fermi gas contribution. The deviations are more transparent in the case of the internal energy per particle due to cancellations of the free Fermi gas terms in the free energy and entropy. The virial and perturbative results are closer at larger temperatures, since the EoS is less sensitive to the physics of large scattering lengths at higher momentum scales.

Figure 3: (Color online) Interaction contribution to the internal energy per particle, , in pure neutron matter at and low densities. The different lines correspond to the results from microscopic chiral nuclear interactions at first order (labeled “HF”) and second order in Eq. (3) as well as the virial expansion truncated at second order (VEoS) and with uncertainty bands obtained by estimating the third-order term. The n3lo414 and n3lo450 results are almost identical.

The differences between the virial and the perturbative results for the internal energy per particle are examined more closely in Fig. 3 for . To depict the deviations more clearly we have subtracted the noninteracting contributions, i.e., the quantity shown is . The virial results include uncertainty bands obtained from estimating the neglected third virial coefficient as . We also show the perturbative results at the Hartree-Fock level [first order in Eq. (3)]. One sees that compared to the Hartree-Fock results the inclusion of second-order contributions leads to much closer agreement with the virial expansion. The second-order calculation still slightly underpredicts the attractive interaction contributions, in contrast to the pseudopotential approach based on nucleon-nucleon scattering phase shift data that was explored in Ref. Rrapaj et al. (2015). We conclude that while the perturbative approach cannot fully capture the large scattering length physics of low-density neutron matter, the resulting errors are reasonably small when second-order contributions are included.

In recent years, the zero-temperature EoS of PNM from chiral nuclear interactions has been studied by numerous authors within various many-body frameworks Tews et al. (2013); Krüger et al. (2013); Wlazłowski et al. (2014); Hebeler and Schwenk (2014); Roggero et al. (2014); Gezerlis et al. (2013); Hagen et al. (2014); Carbone et al. (2014); Coraggio et al. (2013); Baardsen et al. (2013); Kohno (2013); Hebeler and Schwenk (2010); Gezerlis et al. (2014). We compare our results to results obtained from perturbative calculations with various chiral interactions by the Darmstadt group (red band in Fig. 8 in Ref. Krüger et al. (2013)) in Fig. 4. In addition to the N2LO chiral three-neutron forces, their calculations also include all N3LO three- and four-neutron interactions. The uncertainty bands in their results were obtained by allowing large variations of the low-energy constants parameterizing the many-neutron forces. One sees that the (almost overlapping) results from n3lo414 and n3lo450 lie within these bands. In Fig. 4 we also show results obtained from auxiliary-field quantum Monte Carlo simulations with chiral N3LO two-body (AFQMC [NN]) and N3LO two-body plus N2LO three-body forces (AFQMC [NN+3N]) by Wlazłowski et al. Wlazłowski et al. (2014). The perturbative and the AFQMC results are very similar at densities , where both are in close agreement with the (fixed-node) quantum Monte Carlo calculations (based on the AV18 potential) of Gezerlis and Carlson Gezerlis and Carlson (2008). However, at higher densities the EoS predicted by the AFQMC calculations (with three-body forces included) is significantly more repulsive. This discrepancy may be (partly) related to systematic errors in the AFQMC treatment (cf. also Ref. Gezerlis et al. (2014)).

Figure 4: (Color online) Energy per particle in pure neutron matter at zero temperature, , obtained from various many-body methods (see text for details). The inset magnifies the behavior at very low densities where quantum Monte Carlo simulations, labeled “QMC [AV18]”, are expected to be most accurate.

Iv Symmetry free energy, entropy and internal energy

From the results for the free energy per particle in homogeneous SNM and PNM the symmetry free energy is obtained via Eq. (1). The symmetry entropy and internal energy are related to the symmetry free energy via and . The results for , and are shown as functions of density at different temperatures in the left column of Fig. 5. In the insets we show the noninteracting contribution to these quantities, i.e.,

(8)

in the case of the symmetry free energy. In the right column of Fig. 5 we show , , and as functions of temperature at different densities.

Figure 5: (Color online) Left column: results for the symmetry free energy , the symmetry entropy times temperature , and the symmetry internal energy , plotted as functions of density. The insets show the noninteracting (free Fermi gas) contribution to the different symmetry quantities. Right column: , , and as functions of temperature at different densities. The lines are interpolated, with calculated data points at .

One sees that in the considered range of densities and temperatures, is a monotonic increasing function of density and temperature. The density and temperature dependence of is more involved. At low densities decreases monotonically with , but for densities a local minimum is found at . 4 The results from n3lo414 and n3lo450 are very similar for densities well below nuclear saturation density, but at higher densities the dependence on the resolution scale becomes significant. In particular, the decrease in the slope of with increasing density is more pronounced in the n3lo450 results. The dependence of approximately balances that of , and as a result their sum, the symmetry internal energy , increases with density but varies only very little with temperature. At densities near nuclear saturation density the deviations of from its value at zero temperature are below .

In Fig. 6 we show the symmetry quantities with the noninteracting contributions subtracted, i.e.,

(9)

as functions of temperature at different densities. In both cases the interaction contributions tend to counteract the temperature dependence of noninteracting contributions, cf. the insets in Fig. 5. In the case of (and also ) the noninteracing contributions dominate, but in the case of the size of the noninteracting contributions and the ones from chiral nuclear interactions is more balanced, and the dependence of both contributions approximately cancels each other, leading to the observed approximate temperature independence at densities near nuclear saturation density.

Figure 6: (Color online) Temperature dependence of the interaction contributions to the symmetry free energy, , and the symmetry internal energy, at different densities. The lines are interpolated, with calculated data points at .

In Fig. 7 we compare our results for the symmetry (free) energy at zero temperature to the results obtained by Drischler et al. Drischler et al. (2014) from calculations of the EoS of neutron-rich matter using several renormalization group–evolved chiral nuclear interactions. For comparison we also show the results from microscopic calculations within a variational approach by Akmal et al. Akmal et al. (1998) based on the AV18 two-body and the Urbana UIX three-body potential. 5 While the results of Drischler et al. are compatible with our results, the calculations by Akmal et al. predict a symmetry energy that deviates visibly from the n3lo414 and n3lo450 results. In Fig. 7 we also show recent empirical constraints obtained from the analysis of isobaric analog states and neutron skins (IAS+NS) Danielewicz and Lee (2014). One sees that the n3lo414 and n3lo450 results lie in the IAS+NS bands in the entire constrained density region .

Figure 7: (Color online) Symmetry (free) energy as a function of density at zero temperature, . The results from the chiral nuclear interactions n3lo414 and n3lo450 are compared to those of Drischler et al. Drischler et al. (2014) and Akmal et al. Akmal et al. (1998). Also shown are empirical constraints from the analysis of isobaric analog states and neutron skins (IAS+NS).

For densities close to nuclear saturation density the symmetry (free) energy at zero temperature is usually expanded around in terms of :

(10)

where is called the slope parameter, and the symmetry incompressibility. The density where the ground state energy per particle in isospin-asymmetric nuclear matter has a local minimum is related to the parameters in the above expansion via (cf. Ref. Chen et al. (2009)).6 The corresponding incompressibility obeys the approximate relation

(11)

where is usually called the isobaric incompressiblity. In recent years, much effort has been invested in determining the parameters in Eqs. (10). The empirical values of and to a lesser degree also are relatively well constrained (values from Lattimer and Lim (2013), see also Newton et al. (2013, 2014); Lattimer and Steiner (2014); Hebeler and Schwenk (2014); Tews et al. (2013)), whereas experimental determinations of suffer from large uncertainties. For instance, from measurements of neutron skin thicknesses Centelles et al. (2009) the value was obtained, which is compatible with the giant monopole resonance measured in Sn isotopes Li et al. (2007) giving . Theoretical studies using a selection of Skyrme interactions however led to an estimate of Chen et al. (2009).7 Our results for , and are given in Table 2. They are in agreement with the mentioned constraints.

n3lo414
n3lo450
Table 2: Value of the (zero-temperature) symmetry energy at saturation density , the slope parameter , and the isobaric incompressibility , extracted from the results obtained from the sets of chiral nuclear two- and three-body interactions n3lo414 and n3lo450.

V Thermodynamics of isospin-asymmetric nuclear matter

In this section we examine the thermodynamic equation of state of isospin-asymmetric nuclear matter (ANM). The free energy per particle in homogeneous ANM is calculated as follows. The dependence of the nonrelativistic free Fermi gas contributions on the isospin asymmetry is treated exactly, while the relativistic correction term and the interaction contributions are assumed to have a quadratic dependence on the isospin asymmetry:

(12)

This approach is motivated in Sec. V.1, where we investigate in detail the dependence on isospin asymmetry of the noninteracting contributions and .

In Sec. V.2 we then discuss the construction of the spinodal in ANM and present our results for the trajectory of the critical temperature. The dependence of the neutron and proton chemical potentials on isospin asymmetry is examined in Sec. V.3, and we determine the neutron drip point in cold nuclear matter. Finally, in Sec. V.4 we show results for the free energy per particle and pressure in isospin-asymmetric nuclear matter and determine the values of and where an isolated drop of liquid nuclear matter becomes unstable.

v.1 Isospin dependence of free nucleon gas

Here we examine the dependence of the noninteracting contributions to the free energy per particle in homogeneous ANM.8 We compute the four leading terms in an expansion of and in powers of . The results show that the accuracy of the quadratic approximation for decreases significantly with increasing temperature, which necessitates the exact calculation of this contribution. The relativistic correction term on the other hand can be safely approximated via .

The noninteracting contributions to the free energy density, and , can be expressed in terms of polylogarithms , i.e.,

(13)
(14)

where are the neutron and proton effective one-body chemical potentials, and . For given values of , and the effective one-body chemical potentials are uniquely determined by .

The noninteracting contribution to the free energy per particle, , as a function of temperature , nucleon density and isospin asymmetry can be expanded9 in powers of around its value in SNM ():

(15)

The various expansion coefficients are given by

(16)

Setting in Eq. (15) we obtain the noninteracting symmetry free energy as the sum of the above coefficients, i.e.,

(17)

Here, we have introduced the weight factors as a means to specify the relative size of the different expansion coefficients.

The rules of multivariable calculus lead to the following expression for the derivative of order of , , at fixed density and temperature:

(18)

where the functions are defined recursively as

(19)

with the Kronecker delta. The expressions to start the recursion are

(20)

One then obtains for the expression

(21)

where , with the nucleon effective one-body chemical potential, which is uniquely determined by .

The results for the first weight factor and the ratios for are displayed in Fig. 8. Note that the limits and do not commute. The limit of the symmetry coefficients at finite temperature is given by

(22)

which comes entirely from the logarithmic terms, , in the expression for . The asymptotic behavior of the weight factors for is determined by the fact that the and limits of coincide.10 As apparent from Fig. 8, the asymptotic behavior sets in at relatively low values of , causing the convergence rate of the quadratic expansion of to decrease significantly with increasing temperature. This behavior is entirely caused by the presence of the logarithmic terms in the nonrelativistic contribution, i.e., by the first term proportional to the sum of the effective one-body chemical potentials, . The limit of the second term in proportional to equals , which is independent of ; the convergence rate of the expansion in powers of of this term alone increases very strongly with increasing temperature. Similarly, the limit of equals , and the quadratic approximation of the relativistic correction term becomes increasingly accurate with increasing temperature.

For comparison, in Fig. 9 we show the results for the weight factors in the analogous expansion in powers of of the noninteracting contributions to the internal energy per particle, and . The noninteracting contributions to the neutron and proton internal energy densities are given by

(23)
(24)

Both and become independent of as , and the convergence rate of the expansion of both terms increases strongly with increasing temperature (the increase is significantly more pronounced in the case of , which is reflected in the increase of the deviations between relativistically improved and the nonrelativistic results with increasing values of ).

Figure 8: (Color online) Temperature dependence of the first weight factor and the ratios for at different densities, corresponding to the expansion of in powers of . The full lines show the results with the relativistic correction term included, the dotted lines the nonrelativistic results. Note that the deviations between the relativistically improved and the nonrelativistic results decrease with increasing values of , indicating the opposite convergence behavior of the expansion in powers of of and , respectively.
Figure 9: (Color online) Temperature dependence of the first weight factor and the ratio , corresponding to the expansion in powers of of the noninteracting contributions to the internal energy per particle, and . Note the logarithmic scale in the second plot.

v.2 Spinodal and critical temperature

The onset of spinodal instability is associated with the violation of a number of equivalent stability criteria Beegle et al. (1974, 1975), which are derived from the fundamental principle of maximum entropy. In the canonical representation the corresponding stability requirement is that the free energy density at fixed temperature is a convex function of the component densities and , implying that the Hessian matrix has no negative eigenvalues. In the unstable region delineated by the spinodal the (analytical) free energy density is concave; in the metastable region between the spinodal and binodal (coexistence boundary) the free energy density is locally convex and the system is protected against phase separation by a nucleation barrier. The Hessian matrix is given by

(25)

Its eigenvalues are given by

(26)

The signs of the eigenvalues are invariant under (linear) basis transformations. From the data they are readily evaluated using as independent density parameters (nucleon density) and (isospin asymmetry density). In this basis the Hessian matrix becomes diagonal at with eigenvalues and (this result depends on the neglect of isospin-symmetry breaking effects). Hence, in SNM the region inside the spinodal corresponds to a negative isothermal compressibility , where is a stability criterion for a pure substance.

The exact expressions for the nonrelativistic free Fermi gas contribution to the Hessian matrix components are given by11

(27)
(28)

where again and . In the limit the proton density vanishes, , and the proton effective one-body chemical potential diverges (at finite ), , thus . Hence, the exact calculation of the nonrelativistic free Fermi gas contribution leads to divergent behavior of the Hessian components in the limit of vanishing proton concentration. The same divergent behavior is obtained in the exact calculation of at zero temperature. The unstable region then vanishes at a value for all values of . This constraint is lost if the free Fermi gas contribution is approximated by truncating the expansion in powers of of at a finite order, e.g., at first order as in the usual quadratic isospin approximation.

The evolution of the critical temperature where (for a given value of ) the unstable concave region vanishes is depicted in Fig. 10. The results from n3lo414 and n3lo450 are very similar; in both cases the critical lines end approximately at an isospin asymmetry or a proton concentration . The value exceeds the critical line endpoints obtained in Refs. Ducoin et al. (2006); Fedoseew and Lenske (2015); Wu and Ren (2011) using different phenomenological models. We note that in nuclear matter with the coexistence region does not vanish at the critical temperature but at a higher temperature , the so-called maximum temperature Hempel et al. (2013); Ducoin et al. (2006). The existence of a neutron drip point (see Sec. V.3) entails that at zero temperature the binodal extends to over a finite region of densities or pressures. The trajectory of the maximum temperatures therefore reaches its zero-temperature endpoint at vanishing proton fraction .

For comparison, in Fig. 10 we also show the trajectories of the temperature where the region with negative isothermal compressibility vanishes at fixed . The exact expression for the nonrelativistic free Fermi gas contribution to is given by

(29)

For both n3lo414 and n3lo450 the trajectories end at approximately or , which exceeds the values from Ref. Fiorilla et al. (2012) obtained in an in-medium chiral perturbation approach and from Ref. Drews and Weise (2015) obtained by applying the functional renormalization group to a chiral nucleon-meson model.

Figure 10: Trajectories of the critical temperature determined from n3lo450 and n3lo414. The trajectories end at . Also shown are the trajectories of the temperatures where the region with negative isothermal compressibility vanishes. The calculated data points are shown explicitly.

v.3 Stable self-bound liquid

From the data the neutron and proton chemical potentials are obtained via

(30)

The results for and are displayed in Fig. 11 for temperatures . One sees that increases and decreases with . The chemical potentials at finite diverge as , but at zero temperature for . The origin of this feature is the collapse of Fermi-Dirac distribution functions into Heaviside step functions at , which eliminates the logarithmic divergence of the free energy per particle at vanishing density Wellenhofer et al. (2014); Ducoin et al. (2006).

The Gibbs conditions for the coexistence of two bulk phases \@slowromancapi@ (liquid) and \@slowromancapii@ (gas) in mutual thermodynamic equilibrium are

(31)

Whereas at finite liquid-gas equilibrium corresponds to finite values of density and proton fraction in both phases, the vanishing of at vanishing density entails that at the Gibbs conditions for the neutron and proton chemical potentials can in most cases not be satisfied, leading to a gas phase that is either empty (vacuum) or contains only neutrons Ducoin et al. (2006); Lattimer and Ravenhall (1978). The neutron drip point, , is given by the value of where the neutron chemical potential at vanishing temperature and pressure becomes positive. For isospin asymmetries an isolated drop of cold liquid nuclear matter is stable (in equilibrium with the vacuum), defining a stable self-bound state. As seen from Fig. 11, in our equation of state neutron drip occurs at an isospin asymmetry or a proton concentration , which is similar to results obtained with effective Skyrme interactions Lattimer and Ravenhall (1978); Lamb et al. (1981).

Figure 11: (Color online) Neutron and proton chemical potentials, and , in (homogeneous) isospin-asymmetric nuclear matter at temperatures , calculated using n3lo414 (solid lines) and n3lo450 (dash-dot lines). Stable self-bound states are shown as thick bright red lines with circles (full circles for n3lo414, open circles for n3lo450); the lines end at the neutron drip point .

v.4 Metastable self-bound liquid

The zero-temperature results for the (ground state) energy per particle and the pressure are displayed in Fig. 12 as functions of the nucleon density for different values of . The trajectory of the points where the energy per particle has a local minimum and the pressure is zero is shown explicitly. These points correspond to the properties of a drop of cold liquid nuclear matter surrounded by vacuum. For the cold drop is stable (the local energy minimum lies on the binodal, cf. Sec. V.3), and for the local energy minimum lies in the metastable region between the binodal and the spinodal. We refer to the point where the trajectory of the local energy minima encounters the spinodal as the fragmentation point (FP). The energy per particle at neutron drip and at the fragmentation point is and , respectively. For comparison we follow the local energy minima also into the unstable spinodal region; i.e., we show also the point where both derivatives of the (analytical) energy per particle vanish (saddle point, SP) at as well as the local energy minimum at . For the energy per particle is positive at all (finite) densities, and for the pressure is a semipositive definite function of density.

Figure 12: (Color online) Energy per particle and pressure in homogeneous isospin-asymmetric nuclear matter at zero temperature, calculated using n3lo414 (solid lines) and n3lo450 (dash-dot lines). The trajectories of the local energy minima are shown as thick dark gray (bright red below neutron drip) lines with circles (full circles for n3lo414, open circles for n3lo450). The trajectories end at the fragmentation point . The inset magnifies the behavior of the pressure at low densities.

The finite-temperature results for the free energy per particle and the pressure are shown in Fig. 13 for temperatures . At finite the trajectories of the local free energy minima lie entirely in the metastable region; an isolated drop of hot liquid nuclear matter has to be stabilized by a surrounding nucleon gas. At the trajectory ends at , defining the fragmentation temperature for nuclear matter with proton concentration . The saddle point is located at for n3lo414 and for n3lo450. For no local free energy minimum exists; for the analytical free energy per particle is a monotonic increasing function of density for all values of (cf. Fig. 15).

The relation between the spinodal, the binodal, and the trajectories of the local free energy minima is illustrated in Fig. 14. The two plots in Fig. 14 represent isoplethal () and isothermal cross sections of the respective surfaces (spinodal, binodal, surface of local free energy minima) in space (cf. also Refs. Müller and Serot (1995); Ducoin et al. (2006)). In the second plot we also show the surface with divergent isothermal compressibility , which corresponds to the violation of the stability criterion for a one-component system. Hence, the saddle point (SP) where both derivatives of the free energy per particle with respect to the nucleon density vanish coincides with the fragmentation point (FP) where a liquid drop becomes unstable only for where nuclear matter behaves like a pure substance. For the more restrictive two-component stability criteria are needed Beegle et al. (1974) ( is not a relevant stability criterion in that case) and the SP is located in the interior of the spinodal.

Finally, in Fig.