# 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}

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 region^{2}

n3lo414 | -15.79 | 0.171 | 223 | 17.4 | 0.066 | 0.33 |

n3lo450 | -15.50 | 0.161 | 244 | 17.2 | 0.064 | 0.32 |

## 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.

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}

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)).

## 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.

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}

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.

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}

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}

(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}

n3lo414 | |||
---|---|---|---|

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}

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 expanded^{9}

(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}

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 ).

### 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 by^{11}

(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.

### 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).

### 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.

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.