Liquid-gas phase transition in nuclear matter from realistic many-body approaches

The existence of a liquid-gas phase transition for hot nuclear systems at subsaturation densities is a well established prediction of finite temperature nuclear many-body theory. In this paper, we discuss for the first time the properties of such phase transition for homogeneous nuclear matter within the Self-Consistent Green's Functions approach. We find a substantial decrease of the critical temperature with respect to the Brueckner-Hartree-Fock approximation. Even within the same approximation, the use of two different realistic nucleon-nucleon interactions gives rise to large differences in the properties of the critical point.


I. INTRODUCTION
The study of the thermal properties of nuclear matter in terms of realistic many-body approaches has received little attention in the literature, in spite of its potential applications in astrophysics and in the physics of heavy ion collisions. Dense hadronic matter appears in some astrophysical scenarios and it can be particularly hot in the very early stages after the creation of a proto-neutron star in a type-II supernova explosion. Typical temperatures for such systems are in the range of 20 − 60 MeV [1]. It is possible that some astrophysical observables are influenced by the presence of such high temperatures, because of the modifications induced to the bulk and microscopic properties of dense matter. As an example, the gravitational wave spectrum generated in a neutron star merger might be influenced by the temperature dependence of the equation of state [2]. Also, the cooling of the neutron star after its birth is dominated by neutrino emission. This involves a series of processes which are particularly sensitive to the formation of Cooper pairs in the nuclear medium [3]. Consequently, astrophysical observations could help us to constrain the temperature dependence of the microscopic and bulk properties of dense matter.
There is however another way to access the thermodynamical (TD) properties of nuclear and hadronic systems. "Hot" nuclei are created in experimental facilities, in the collisions of heavy ions at intermediate energies [4]. There has been an increasing effort to interpret experimental data concerning these heavy-ion collisions in terms of the equation of state of nuclear matter [5,6]. More specifically, multifragmentation reactions at intermediate energies are used to access the properties of thermally equilibrated "blobs" of nuclear matter [7]. The evidence gathered in these experiments points towards the existence of a liquidgas phase transition for nuclear systems at densities below the empirical saturation density, ρ 0 = 0.16 fm −3 , and temperatures around ∼ 6 − 9 MeV [8]. Empirically, the liquid-gas phase transition is usually discussed in terms of a plateau in the caloric curve for different types of reactions at different energies [9], although the presence of such plateau may also be explained in terms of a density and temperature dependent effective mass [10]. Statistical models, which assume an equilibrated thermal freeze-out density, have had certain success in describing experimental results [11]. In any case, a common underlying idea in all these discussions is the assumption that some sort of thermal equilibrium is reached at a given stage of the multifragmentation reaction.
A first step towards a full theoretical understanding of the thermal properties of nuclear systems can be achieved by studying symmetric nuclear matter (an ideal, infinite, highdensity system composed of the same amount of neutrons and protons interacting via the strong force) at finite temperatures, because the features of the liquid-gas phase transition are in general easier to study in the homogeneous system. This is obviously a very crude approximation and can only be taken as a guide for a better theoretical understanding of the thermal properties of dense matter. As an example, the low-density phase of dense nuclear matter is not a homogeneous gas of nucleons and should instead be described in terms of droplets of light nuclei (deuterons, alpha particles). In the case of real experiments, moreover, finite size effects are as important as the bulk properties in determining the critical behavior. Finite nuclei can only be excited up to some limiting temperature above which Coulomb effects together with the decrease in surface tension lead to their thermal dissociation [12,13]. All in all, it is not an easy task to link the properties of the phase transition in the homogeneous case to those of finite nuclei [14,15,16]. In this paper, however, we shall concentrate on studying the liquid-gas phase transition in the ideal case of nuclear matter. Our basic goal is to discuss to what extent the properties of this transition depend on the many-body description chosen as well as on the underlying NN interaction.
Traditionally, the studies of nuclear systems at finite temperature have been performed using effective interactions within mean-field theories, either relativistic [17,18] or nonrelativistic [19,20]. The temperature dependence in these approaches arise mainly from the modification of the zero temperature step-like momentum distributions, which become Fermi-Dirac distributions. As a consequence, the mean-field and the bulk properties have a rather simple dependence on temperature. Thermal effects on the correlations between the strongly interacting nucleons are completely ignored. The effective interactions, which are fitted to describe the bulk properties of nuclei at zero temperature and thereby account for correlation effects in a phenomenological way, have no temperature dependence. This is in striking contrast with the more microscopically founded many-body calculations [21,22,23], where the in-medium interaction is found by using some sort of Pauli blocking. Such blocking effects are weakened by temperature, therefore giving rise to non-trivial temperature dependences in the microscopic and macroscopic properties of dense matter.
There are only a few realistic many-body calculations for nuclear matter at finite temperature [24,25,26]. Some of these have consistency problems problems, since they are not truly based on finite temperature many-body theory, but on naive extensions of zero temperature approaches to the non-zero temperature domain. The Green's functions approach is however based on the perturbative expansion of the single-particle propagator at finite temperature [27] and therefore it is a well-grounded approach which allows for systematic improvement. In contrast, the variational approach is based on an explicit incorporation of the two-body nucleon-nucleon (NN) correlations in the nuclear wave function. Its extension to finite temperatures has traditionally relied on the "frozen" correlation approximation, i.e. the correlation functions have been assumed to be the same at zero and at finite temperature [24,28,29]. This might be adequate for the high-density low-temperature phase, but it remains to be seen if it offers an appropriate description of the high-temperature phase. Recently there has been a substantial effort to discuss the more formal aspects associated to the variational approach at finite temperatures [30]. Let us also note that the variational approach is restricted by nature to deal with local potentials, while the SCGF can in principle deal with any sort of realistic two-body force. The Brueckner-Hartree-Fock (BHF) approach has also been used in the study of hot nuclear systems [25], although the Brueckner-Bethe-Goldstone expansion, in which the BHF approach is based, is only valid at zero temperature. The standard finite temperature generalization of this approach is, in a way, phenomenological and relies on the replacement of all the step-function momentum distributions of the zero temperature case with Fermi-Dirac ones. This finite temperature extension is however not well-defined at a fundamental level, since it does not take into account the contributions of anomalous diagrams [31]. A consistent Brueckner-like approach at finite temperature is given by the Bloch-de Dominicis formalism [32,33,34], which has only recently been applied to the nuclear many-body problem at finite temperatures by Baldo and coworkers [13,22]. Relativistic effects have also been explored in the framework of an extension to finite temperature of the Dirac-Brueckner-Hartree-Fock approach [35,36].
A good alternative for a suitable microscopic many-body description of hot correlated systems is provided by the perturbation expansion of the one-body Green's functions [37], which relies on the generalization of the Wick theorem at finite temperature [38]. Due to the strong short-range repulsion of the NN interaction, the minimum meaningful scheme which can describe nuclear matter is provided by the ladder approximation. This arises from a decoupling of the three-body Green's function in terms of one-and two-body propagators [38] and can be cast as a set of self-consistent equations that describe the in-medium mod-ifications of the nucleon due to the presence of the surrounding nucleons. These equations lead to an approximation that goes beyond the mean-field and the quasi-particle pictures, i.e. the off-shell effects and the fragmentation of single-particle states are fully taken into account. This Self-Consistent Green's Functions (SCGF) approach is well established in nuclear physics [39] and has already been applied to study the microscopic properties of nuclear matter at finite temperatures [40,41,42,43]. A major motivation for these studies has been the fact that, at non-zero temperature, one can avoid the numerical and physical problems associated to the neutron-proton pairing instability [44,45,46]. The SCGF, however, can also be used to study the TD properties of the system in the normal phase, by making use of the Luttinger-Ward (LW) formalism. This approach leads to thermodynamically consistent results, once the effects of correlations in the entropy have been carefully taken into account [26,47,48]. In the following we shall use this formalism to study the properties of the liquid-gas phase transition in nuclear matter. Note however that the effect of three-body forces is not included in our approach and, as a consequence, the saturation properties of nuclear matter are not well reproduced. The lack of three-body forces will also have an impact on the liquid branch of the phase transition, thus modifying to a certain extent our predictions for the critical properties of nuclear matter. Our results should be considered as a theoretical study quantifying the importance of short-range correlations on the liquid-gas phase transition.
In the next section, we briefly summarize the SCGF approach at finite temperature and we discuss under which approximations the standard BHF at finite temperatures can be obtained from it. The third section will deal with the application of the LW formalism to the calculations of the TD properties of a correlated system of nucleons. The results for the liquid-gas phase transition will be discussed in the fourth section. Finally, a brief summary will be given in Section V.

II. SCGF AT FINITE TEMPERATURE
The key quantity in many-body Green's functions theory is the single-particle propagator which, in the grand-canonical ensemble, is defined according to: where we have introduced the density matrix operator: and the partition function: In the previous equations, β = 1/T denotes the inverse temperature and µ is the chemical potential. T stands for a time-ordering operator and the traces Tr{·} are to be taken over all energy and particle number eigenstates. One can express the single-particle propagator in Fourier-energy space by means of the spectral decomposition: where the finite-temperature equivalent of the T = 0 hole spectral function is given by the Lehmann representation: The main difference with respect to the zero temperature case comes from the average over the thermal bath in the initial states. A similar definition holds for A > (k, ω), with the replacement a k → a † k . In thermal equilibrium, both spectral functions are related by the Kubo-Martin-Schwinger relation: In contrast to the zero temperature case, the energy domains of A < (k, ω) and A > (k, ω) are not separated by the Fermi energy and both spectral functions are defined for all energies.
The total spectral function, A(k, ω), is given by the sum of the two functions, A < and A > , and therefore it can be expressed in terms of the values of G close to the real axis: (where we have introduced the notation ω + = ω +iη). Since the spectral function completely determines the one-body propagator via the previous equation, all the one-body properties of the system can be expressed in terms of it.
In the medium, the single-particle Green's function can be obtained from Dyson's equation: where Σ(k, ω) denotes a complex self-energy. The self-energy accounts for the interactions of a particle with the other particles in the medium. It fulfills the following spectral decomposition: where z is a complex variable and the term Σ HF (k) is a real energy-independent generalized Hartree-Fock contribution: with n(k) being the single-particle momentum distribution: and where f (ω) = e β(ω−µ) + 1 −1 stands for the Fermi-Dirac distribution. The imaginary part of the self-energy, necessary to compute the second term in Eq. (9), is obtained by letting z → ω + in Eq. (9) and it is related to the effective two-body NN interaction in the medium (the so-called scattering T -matrix): Note the presence of a Bose- of the symmetric treatment of particle-particle and hole-hole correlations.
The effective in-medium interaction is calculated in the ladder approximation. This accounts for the repeated scattering of particles in the medium and it is well-suited for the low-density strong-interaction regime of interest for nuclear matter [27]. The T -matrix is determined by the solution of the integral equation: where the intermediate propagator accounts for the propagation of two non-interacting but dressed nucleons: To reduce the dimensionality of Eq. (13), one usually relies in the standard partial wave decomposition. An extra simplification is achieved by using an angle average of the twobody propagator with respect to the center of mass and relative momentum of the two colliding particles [48,49].
Equations (9)-(14) form a closed set of equations that can be solved self-consistently. In terms of numerics, it is advantageous to work at constant density, and therefore we supplement the previous set of equations with the normalization of the momentum distribution: where ν denotes the degeneracy of the system (ν = 4 in the case of symmetric nuclear matter). Once convergence is reached in the self-consistent procedure for a given temperature and density, one has access to the spectral function A(k, ω), which, loosely speaking, describes the probability of finding a nucleon in the medium with momentum k and energy ω. At this point, one can calculate several micro-and macroscopic properties of the system.
The momentum distribution, for instance, can be computed using Eq. (11). The energy per particle is also accessible from the Galitskii-Migdal-Koltun sum-rule: which is valid for a Hamiltonian with only two-body interactions. The importance of selfconsistency in the calculations stems from the fact that it immediately leads to the conservation of both micro-and macroscopic properties [50]. In addition, it guarantees the fulfillment of the sum-rules for the one-body spectral function [51,52].
In this paper, we shall make comparisons between the finite temperature SCGF and BHF approaches. The latter can be formally derived from the first by performing some particular approximations. First, one has to assume that, for a given momentum, all the strength of the spectral function is accumulated in one energy, the BHF single-particle energy. This simplifies the calculation of the non-interacting twobody propagator of Eq. (14), which becomes a finite-temperature Pauli blocking factor involving both particle-particle and hole-hole propagation. Since in the BHF approach only intermediate particle-particle states are considered, the phase space factor needs to be . Finally, in the BHF self-energy, one does not take into account the contribution of the Bose function in Eq. (12). After the BHF equations are iterated and consistency is reached, one obtains a single-particle spectrum and an in-medium G-matrix interaction, the real part of which is used to obtain the energy per particle of the system: The SCGF and the BHF approaches with two-body NN interactions cannot reproduce the saturation properties of nuclear matter at zero temperature, due to the lack of repulsive contributions, most probably those coming from three-body forces [22,53]. The region of interest for liquid-gas phase transition studies is however in the low-density regime (critical densities are typically 1 2 to 1 3 of ρ 0 ) and this density regime should not be strongly modified by the presence of three-body forces. Estimates of the importance of the three-body forces in the liquid-gas regime have been performed within the BHF approximation, indicating small modifications to the critical properties [22,53].
To assess the model-dependence of the properties of this transition, we will compare the results obtained with two different NN interactions, the CDBONN [54] and the Argonne V18 [55] potentials. Although both of them reproduce the scattering phase-shifts up to about 300 MeV, they have very different short range cores, off-shell structure and tensor components.
The many-body calculations depend on these details and therefore the properties of dense matter, in particular the critical properties of the liquid-gas phase transition, will be different for the two interactions. To our knowledge, this is the first time that two different realistic interactions are used within the SCGF approach to study the liquid-gas coexistence and the critical properties. In all the calculations quoted in the following and for both the SCGF and the BHF approximations, partial waves up to J = 8 have been included. The in-medium effective interactions have been computed with J ≤ 4 and the Born approximation has been used for J > 4.

III. THERMODYNAMICAL PROPERTIES OF NUCLEAR MATTER
A complete TD description of the system requires the computation of the free energy, F = E −T S. As indicated in the previous section, the internal energy in the SCGF approach can be calculated from the one-body propagator via the Galitskii-Migdal-Koltun sum-rule.
Therefore, a suitable method for the calculation of the entropy is required to describe the thermodynamics of the system. The LW formalism can be used to find an expression of the grand-canonical potential, Ω, in terms of dressed propagators [31,56,57]. Because this expression of Ω is stationary with respect to variations of the one-body Green's functions, one can easily compute the entropy from the derivative S = − ∂Ω ∂T µ . This entropy can be split in two terms, S = S DQ + S ′ . The first one is the so-called dynamical quasi-particle (DQ) entropy density: defined as the convolution of a statistical factor, , and a spectral function B(k, ω), related to the single-particle spectral function, A(k, ω), and the self-energy by the following equation: This expression takes into account the correlations of the dressed particles in the medium, since it includes finite width effects. In this paper, we shall make the assumption that S ′ is negligible. As a matter of fact, it has been shown by Carneiro and Pethick that its effects are constrained by phase space restrictions [57]. This assumption is also confirmed by the fact that the results ignoring the contribution of S ′ are thermodynamically consistent [47], as we shall see in the following. The free energy will be computed from the difference of the Galitskii-Migdal-Koltun sum-rule energy and the DQ entropy, F = E GM K − T S DQ . In the BHF approach, the entropy should not include any effect due to the widening of the quasi-particle peak and therefore it will be computed from the mean-field expression: The free energy per particle, together with the energy per particle and the chemical potential, are shown in Figs. 1 and 2 as a function of the density at constant temperature.
In Fig. 1 we display the results at T = 8 MeV for the Argonne V18 interaction, while in Fig. 2 we consider the CDBONN interaction at T = 10 MeV. Both temperatures are sufficiently below the corresponding critical temperatures, as we will see. Note that in Let us first discuss the differences between the two many-body approaches. It is already well-established that hole-hole propagation, which is included in the SCGF but not in the BHF approach, yields a repulsive contribution to the energy per particle of about 4 − 6 MeV close to saturation density [41,58]. Since this repulsive contribution tends to increase with nuclear density, one obtains a smaller saturation density in the SCGF as compared to the BHF approach. A similar effect is also observed for the free energy per particle, since both the dynamical quasi-particle and the BHF entropies are quite close to each other [47]. It is interesting to compare the free energy and the chemical potential to check the fulfillment of TD consistency. Some properties of the system can be computed either microscopically (from Green's function theory) or macroscopically (from the TD properties of the system). A TD consistent many-body approximation will yield the same result for both of them. A very sensitive quantity to this test is the chemical potential. On the one hand, it can be computed microscopically from the normalization of the momentum distribution, Eq. (15), giving rise to the microscopic chemical potential,μ (diamonds in Figs. 1 and 2).
On the other hand, one can compute it from the derivative of the free energy with respect to the number of nucleons at constant temperature, µ = ∂F ∂N (dotted lines). The differences betweenμ and µ for the BHF approach can be larger than 15 MeV, showing its lack of consistency. Note that, in particular, the Hugenholtz-van Hove theorem is violated, i.e.μ does not coincide with F A at its minimum. The violation seems to be larger for CDBONN, ∼ 20 MeV, than for Argonne V18, ∼ 10 MeV. The SCGF results, however, fulfill TD consistency with less than one MeV error in a wide density range for both interactions. Note that µ has been computed by fitting a fourth-order polynomial to the free energy density and determining the derivatives from this polynomial fit. The lack of accuracy in the fit is responsible for the deviations at very low densities. In any case, the numerical implementation of the ladder approximation by means of the SCGF scheme leads to TD consistent results, independently of the NN potential under consideration [59].
The pressure is shown as a function of the density in Fig. 3 for several temperatures and for both the SCGF (left panels) and BHF (right panels) approaches. The upper panels, (a) and (b), correspond to Argonne V18 and the lower ones, (c) and (d), to CDBONN. The pressure is obtained from the TD relation, p = ρ(µ−F/A). Due to the conserving properties of the SCGF approach, we can compute the pressure at every density and temperature by using the microscopic chemical potential in the previous expression. For BHF, however, µ has to be computed as a numerical derivative of the free energy with respect to ρ. Note that numerical problems can appear in the low density limit, due to logarithmic density dependences in this region.
In general, one can say that the repulsive effect of the hole-hole propagation in the freeenergy, which we have already discussed above, is translated into larger pressures in the SCGF as compared to the BHF approach, especially at large densities. The SCGF method, therefore, yields a stiffer equation of state than the BHF approach. By construction, the TD chemical potential, µ, crosses the free-energy curve at its minimum, thus yielding a point of zero pressure. This defines the saturation density at each temperature. The repulsive effect induced by the propagation of holes in the SCGF approach is reflected in a smaller saturation density with respect to BHF, i.e. the minimum of the free energy is shifted to smaller densities in the SCGF approach. This effect has already been observed at zero temperature [58] and it appears to hold when thermal effects are taken into account. Above a certain temperature, the minimum of the free energy per particle disappears and the equation P (ρ) = 0 has no solution anymore. This defines the so-called flashing temperature, T f . The first column of Table I gives the flashing temperature for the different approaches and potentials. The SCGF results lead to T f 's which are about 3 MeV lower than the BHF ones.
The differences in pressures from the two many-body approaches are sizeable, but the differences due to the change of potentials are even larger. The pressure for Argonne V18, for instance, increases much more steeply with density than CDBONN does. There are also substantial differences in the saturation densities induced by the two potentials at all temperatures, with Argonne V18 leading to lower saturation densities than CDBONN. This is in agreement with the fact that at T = 0 the BHF saturation density is much larger for CDBONN than for Argonne V18 [60]. Note also that the temperature ranges explored in the upper and lower panels of Fig. 3 are not the same. The temperature dependences induced by the two potentials are therefore rather different. For instance, the flashing temperatures for Argonne V18 are about 4 − 5 MeV lower than for CDBONN with both SCGF and BHF.
This suggests that the different off-shell and tensor components of the NN forces do not only affect the properties of nuclear matter at zero temperature, but also its TD properties in an important way. It remains to be seen if the experimental knowledge gathered about the thermodynamical properties of nuclear systems can provide additional information to constrain the NN interaction and select the proper many-body approach to be used in their description.
For both approaches and potentials, the pressure decreases with density in a given range.
This signals the existence of a mechanical instability, which is associated to a first order liquid-gas phase transition. The properties of this transition are studied in the following.

IV. LIQUID-GAS PHASE TRANSITION
A physical interpretation of the TD unstable zone is customarily obtained by making use of the Maxwell construction. For each temperature, one should find the gas and liquid density, for which the equations µ(ρ g ) = µ(ρ l ) and p(ρ g ) = p(ρ l ) are simultaneously satisfied.
For a given temperature, the range ρ g − ρ l gives the coexistence region, where the gas and liquid phases coexist at constant pressure and chemical potentials. The spinodal region is defined by the violation of the TD stability criteria, ∂µ ∂ρ > 0 and ∂p ∂ρ > 0. For a one component system, both conditions are equivalent. This spinodal region lies within the liquid-gas coexistence region in the density-temperature plane, and the region between the two curves defines the so called metastable region. Finding the spinodal and coexistence densities at each temperature, one obtains the phase diagrams shown in Fig. 4 for Argonne V18 and in Fig. 5 for CDBONN. Note that the two figures have different vertical scales, because of the large differences in the TD properties induced by the two interactions. Once again, let us emphasize that one should be cautious when using fitting procedures for the free energy. These are necessary to obtain suitable analytical expressions for the derivatives of F (ρ, T ), needed to implement the liquid-gas coexistence conditions. It is also worth noting that, for the temperatures considered, both methods (BHF and SCGF) are numerically stable in a range of densities which covers from very low ones (gas phase, ρ = 0.01 fm −3 ) to relatively large ones (liquid phase, ρ = 0.30 fm −3 ). For the SCGF approach, numerical difficulties related to the pairing instability appear at lower temperatures [44,45].
The critical temperature for the liquid-gas phase transition of symmetric nuclear matter corresponds to the maximum of the spinodal and coexistence lines, which coincide with each other at the critical point. The flashing temperature, T f , always lies below T c and, as commented above, correspond to the maximum temperature at which the pressure presents a node. T f represents the highest temperature that finite nuclei can withstand without thermally dissociating. The critical properties as well as the limiting temperatures for the two approaches and interactions are listed in Table I. Both T f and T c depend strongly on the type of approximation used to study nuclear matter, and also on the NN interaction.
Note that, for a given potential, BHF results always lead to larger flashing and critical temperatures (about 3 − 4 MeV larger) than SCGF. Also, the results for Argonne V18 are, within each approximation, about 4 − 5 MeV lower than those of CDBONN.
At this point we are going to compare our predictions for the critical temperature with those obtained in other approximation schemes. We obtain the largest critical temperature for the CDBONN potential in the BHF approximation, with a value of T c = 23.3 MeV. The critical temperature for CDBONN using the SCGF approach is, however, T c = 18.5 MeV. This is in close agreement with the value of T c = 18 MeV obtained with the somewhat similar Bonn B potential in the Bloch-de Dominicis formalism [13]. The Argonne V18 interaction yields lower critical temperatures when compared to CDBONN, with T c = 18.1 MeV for BHF and T c = 11.6 MeV for SCGF. The first result can be compared with calculations performed with other approximations and the same NN potential. The LOCV results of Ref. [28] correspond to a somewhat larger critical temperature of T c = 22.2 MeV for V18, while a more recent variational calculation with frozen correlations and three-body forces leads to T c = 18 MeV [29], rather close to our BHF result. The authors of Ref.  [35] and [36], respectively, obtained using One-Boson-Exchange interactions fitted to NN data. The nonrelativistic semirealistic model of Ref. [61] found convergence problems, possibly associated to the liquid-gas transition, below a critical temperature as low as 9 MeV.
Can one explain in simple terms the origin of the large differences in critical temperatures between the two approaches and NN interactions? In fact, there are some simple models which try to relate the critical properties to the saturation properties of nuclear matter at zero temperature [14,62]. A particularly useful and simple estimate is obtained from the Kapusta model [15], which supposes that the temperature dependences are quadratic (as in the free Fermi gas close to the degenerate limit) and modulated by an effective mass, m * , which governs the density of states. Under the additional assumption that the zero temperature energy per particle can be characterized by the compressibility, K, one finds that: i.e. the critical temperature increases with K and the saturation density, ρ 0 . The presence of the pairing instability [44,45] [58]. Within the Kapusta model these two features would essentially compensate each other and predict a critical temperature for the SCGF approach which is about the same as the one obtained in the BHF approximation. Note, however, that the SCGF approach leads to an enhancement of the density of states at low excitation energies which could be described in terms of a larger effective mass m * as compared to BHF [63].
Independently of this analysis, it is generally true that the critical temperature is correlated and somewhat close to the binding energy at saturation, T c ∼ E/A [14].  Table I.
The Kapusta model also predicts the value of the critical density, ρc ρ 0 = 5 12 = 0.417. For the BHF calculations, we find ρc ρ 0 = 0.35 for Argonne V18 and ρc ρ 0 = 0.31 for CDBONN, which are somewhat closer to the empirical formula ρc ρ 0 = 1 3 [14]. Note also that the value for Argonne V18 in the BHF approximation lies within the range 0.07 − 0.09 fm −3 quoted in Ref. [22] for the V14 potential. The spread in critical densities is significant when we change from one potential to another, but it does not differ so much when we consider different approximations. The differences are however more drastic for the critical pressure, which changes by almost an order of magnitude when comparing different approximations and potentials. Yet, surprisingly, the discrepancies are somewhat less important for the adimensional parameter pc ρcTc . For a van der Waals equation of state, this parameter is  Table I are quite below this value. Intriguingly, the BHF results seem to lead to the same value, in spite of the large differences in each of their critical parameters.
Finally, let us recall again that we have not considered any three-body forces in the previous calculations. In terms of the phase diagram, one expects that the inclusion of three-body forces will shift the liquid coexistence branch to lower densities in the low temperature phase, but will presumably have a small effect at large temperatures. The effects on the gas phase, if any, would probably be very small in the homogeneous case. In fact, previous evaluations of the three-body effects on the critical temperature obtained within the BHF approach seem to indicate that these are rather small, about 1 − 3 MeV [22,53]. Such a decrease is much smaller than the discrepancies observed here when changing the two-body interaction or the many-body method.

V. SUMMARY AND CONCLUSIONS
We have calculated the TD properties of symmetric nuclear matter within the Self-Consistent Green's Function and the Brueckner-Hartree-Fock approaches for two realistic NN interactions, the CDBONN and the Argonne V18 potentials. The calculations cover a wide range of densities and temperatures. In the SCGF-LW approach, the entropy has been computed within the dynamical quasi-particle approximation, which takes into account the effects of correlations in the width of the quasi-particle peak. A very good agreement between the microscopic and macroscopic chemical potentials is found, highlighting the TD consistency of the SCGF-LW approach at the numerical level. This is in contrast to the BHF approximation, which yields a violation of the Hugenholtz-van Hove theorem by 10 MeV (20 MeV) for the Argonne (CDBONN) potential.
The essential difference between SCGF and BHF is the consistent inclusion of holehole propagation terms in the former approach, leading to non-trivial spectral distribution functions and partial occupation probabilities for states with momenta above and below the Fermi momentum. These hole-hole terms tend to provide some repulsion, which increases with density. This feature is also reflected in the calculated pressure, for which the SCGF approach yields larger values, in particular at higher densities. This implies that the equation of state derived within the SCGF approach tends to be stiffer than the corresponding one evaluated within the BHF approximation. The repulsive effect of the hole-hole terms also leads to a lower flashing temperatures for the SCGF approach as compared to BHF.
When comparing the results between two different NN interactions, we also find substantial differences, larger than those induced by the use of different many-body approaches.
In particular, the Argonne V18 interaction leads to a stiffer equation of state and a lower flashing temperature than CDBONN.
The liquid-gas phase diagram for nuclear matter has been studied for the first time in the framework of the microscopic SCGF approach for two realistic NN interactions and critically compared with the results obtained in the BHF approach. Substantial differences for the critical properties are found when changing the potential and the many-body approximation.
The SCGF leads to critical temperatures which are 5 − 7 MeV lower than those obtained with the BHF approximation. Within the same approximation, CDBONN leads to results which are 6 − 7 MeV larger. For BHF, where T = 0 calculations can be performed safely, we find that the critical density is about a third of saturation density and that the critical temperature can be well approximated by Eq. (21), in terms of the compressibility and the saturation density.
Of course, to have a proper estimation of the critical temperature for finite nuclei, one should also take into account Coulomb effects and the existence of a surface tension. These results would further reduce the critical temperature, by a factor of 1/2 − 1/4 [9,13]. In this paper we have found that realistic calculations allow for a large range of critical temperatures, in the same way that they predict different saturation properties. In particular, the important reduction in critical temperatures found for the Argonne interaction might be relevant when trying to connect the data of multifragmentation reactions with the liquid-gas phase transition for bulk nuclear matter.