Spinodal Instabilities in Asymmetric Nuclear Matter Based on Realistic $NN$ Interactions

A density dependent relativistic mean-field model is determined to reproduce the components of the nucleon self-energy at low densities. This model is used to investigate spinodal instabilities in isospin asymmetric nuclear matter at finite temperatures. The inhomogeneous density distributions in the spinodal region are investigated through calculations in a cubic Wigner-Seitz cell. Compared to results obtained in phenomenological calculations the spinodal region is large, i.e. the spinodal region at zero temperature can reach densities above 0.12 fm$^{-3}$. The predicted spinodal region is concentrated around isospin symmetric nuclear matter and the critical temperature is considerably lower than in the previous microscopic based investigation within a non-relativistic Brueckner-Hartree-Fock approach.


I. INTRODUCTION
The properties of isospin asymmetric nuclear matter systems like their instabilities and phase transitions are presently a topic of great interest [1][2][3][4][5][6][7]. At either high energy or high baryon density one explores the phase transition to quark matter, i.e. the hadrons dissolve into a phase of de-confined quarks. Below a critical temperature or density, the homogeneous fluid consisting of nucleons may condense into nuclei. This transition can be regarded to be of the liquid-gas type. Heavy ion collisions are well suited to investigate these transitions, since the multi-fragmentation phenomenon occurs in these experiments. In the multi-fragmentation phenomenon, highly excited composed nuclei are formed in a gas of evaporated particles. Therefore, this multi-fragmentation phenomenon can be interpreted as the coexistence of a liquid and a gas phase. This liquid-gas phase transition is important in the description of the nucleation processes. Furthermore, it plays an important role in the collapse of supernovae into neutron stars.
The hadrons involved in this transition can be either neutrons or protons. Therefore, one should consider nuclear matter as composed of two different fluids. These phase transitions in nuclear matter are related to the thermodynamic instability regions, which are limited by so-called spinodals. This spinodal instability is dominated by mechanical or isoscalar density fluctuations which lead to a liquid-gas phase separation with restoration of the isospin symmetry in the dense phase. This phenomenon is known as the isospin distillation or fragmentation effect [8], i.e. the formation of large droplets of high density symmetric nuclear matter within a neutron gas with only a very small fraction of protons.
In the past systematic analysis of the stability conditions of asymmetric nuclear matter against the liquid-gas phase transition have been performed using different methods. These approaches can roughly be divided into phenomenological methods and microscopic based methods. Phenomenological studies of spinodal instabilities involve mean-field calculations with effective forces of the Skyrme or Gogny type [1,[9][10][11] and relativistic mean-field calculations using constant and density dependent coupling parameters [3,4,7,12,13]. Some of these models predict that the total density at which spinodal instability appears will decrease if the asymmetry increases, whereas for other models it will increase up to large asymmetry and will finally decrease. However, a common feature of these models is that the spinodal region is in principle situated at densities below ρ B < 0.1 fm −3 . The microscopic based studies are a calculation with low-density functionals [5] based on the Dirac-Brueckner-Hartree-Fock (DBHF) approach [14] with a realistic Bonn A potential [15] and a calculation within a microscopic Brueckner-Hartree-Fock (BHF) approach [6] using the realistic Argonne V18 potential [16] plus a three-body force of Urbana type. Both calculations predict larger spinodal regions compared to those obtained in the phenomenological calculations, i.e. the spinodal region can reach densities above 0.12 fm −3 close to symmetric nuclear matter. The critical temperature in the BHF calculation is T c = 17.5 MeV, whereas finite temperature effects are not considered in the work of Ref. [5].
In the present work, we construct a density dependent relativistic mean-field (DDRMF) model, i.e. a density dependent relativistic Hartree (DDRH) model, which is based on microscopic DBHF calculations [17] using the realistic Bonn A potential [15]. Apart from the application of the results of this DBHF calculation to spinodal instabilities [5], it has also already been used for other applications such as investigations of the non-homogeneous structures in the neutron star crust [18], investigations of hybrid stars [19], and studies of nucleon-nucleus scattering [20]. In the previous DBHF based investigation of the spinodal instability [5], a functional is constructed on a density expansion of the relativistic mean-field theory with fixed coupling constants. However, in the present work the approximation of an expansion is not made and the effective coupling functions, which are obtained from the DBHF self-energy components, are density dependent. Special attention is devoted to reproduce the DBHF self-energy also at very low densities, i.e. ρ B < 0.10f m −3 . In the framework of this DDRH theory based on a DBHF approach, the equation of state at zero and finite temperature can then be determined from an iteration procedure until convergence is reached. Within this model, the temperature effects and the nature of spinodal instabilities are investigated. In addition to these investigations in infinite nuclear matter calculations, the density distribution in the spinodal region is studied in cubic Wigner-Seitz (WS) cell calculations.
The paper is organized as follows. The density dependent relativistic Hartree (DDRH) theory is treated in Sec. II. The parameterization for the DDRH theory based on the DBHF results are discussed in Sec. III. Section IV is devoted to the theory of the spinodal instabilities. The results for the spinodal instabilities are presented and discussed in Sec. V. Finally, we end with a summary and the conclusion in Sec. VI.

II. DENSITY DEPENDENT RELATIVISTIC HARTREE THEORY.
A Lagrangian density of an interacting many-particle system consisting of nucleons and mesons is the starting point of a DDRH theory. The Lagrangian density of the DDRH theory used for our investigation includes as well the isoscalar mesons σ and ω as the isovector mesons δ and ρ. Therefore, the Lagrangian density has the following parts: the free baryon Lagrangian density L B , the free meson Lagrangian density L M , and the interaction Lagrangian density L int : which can be written as In the above Lagrangian density the field strength tensor for the vector mesons is defined as The nucleon field is Ψ and the nucleon rest mass is M . The meson fields are denoted by Φ σ ,Φ δ , A (ω) , and A (ρ) . Moreover, the bold symbols denote vectors in the isospin space. The mesons couple to the nucleons with the strength of the coupling constants g σ , g δ , g ω , and g ρ . The rest masses of these mesons are m σ , m ω , m δ , and m ρ . By minimizing the action for variations of the fieldsΨ included in the Lagrangian density of Eq. (1), one obtains the field equations for the nucleons, Finally the following field equation can be derived Without density dependent meson-baryon vertices the third term vanishes in Eq. (6) and the normal Euler-Lagrange field equation is obtained. From the evaluation of Eq. (6), we obtain the Dirac equation, with Σ the nucleon self-energy and Σ (r) the so-called rearrangement contribution to the nucleon self-energy. The self-energy Σ for the Lagrangian density in Eqs. (1)-(4) can be written as behavior under Lorentz transformations, The Σ s , Σ o , and Σ v components are Lorentz scalar functions. Hence, one can define the following effective quantities and In the on-shell case, the effective energy can also be given by The contributions to the self-energy of the Lagrangian density presented in Eqs.
(1)-(4) can be written as and The thermal occupation factors in Eqs. (14) and (15) are given by in the "no-sea" approximation. The Σ (r) term in Eq. (7) is generated by the third term in Eq. (6). It will only give a nonzero contribution, if meson-baryon vertices are density dependent. This rearrangement contribution Σ (r) can be written as Rearrangement contributions are essential to provide a symmetry conserving approach [21,22], i.e. an approach satisfying energy-momentum conservation and thermodynamic consistency like the Hugenholtz -van Hove theorem.

III. PARAMETERIZATION
We have constructed a density dependent relativistic mean-field (DDRMF) model, i.e. a density dependent relativistic Hartree (DDRH) model. The effective density dependent coupling functions of this DDRH model are determined by the requirement that the results of a DBHF calculation at zero temperature using Bonn A [17] for the components   of the nucleon self-energy at the Fermi momentum are reproduced. The rearrangement terms of DDRH are not included in the fit, since the DBHF approach has no rearrangement contributions.
The density dependence of the coupling functions are parameterized by with x = ρ B /ρ 0 and ρ 0 = 0.16 fm −3 . The DBHF data at densities of ρ B = 0.017, 0.100, 0.197, 0.313, and 0.467 fm −3 at a proton fraction of Y p = 0.4 are used to obtain these parameters of our DDRH model. The meson masses are chosen to be identical to those of the Bonn A potential. All parameters are given in table I. The density dependence of these coupling functions is displayed in Fig. 1. The density dependence is rather weak except at very low densities. The strength for the isovector mesons (ρ and δ) is significantly smaller than for the isoscalar mesons (σ and ω). All coupling functions tend to decrease with increasing density. This may reflect the quenching of the DBHF G-matrix underlying this fit due to Pauli and dispersive effects in the Bethe-Goldstone equation. In Fig. 2 the binding energy of the DDRH model is shown for neutron matter and symmetric nuclear matter. The results are in reasonable agreement with the DBHF results. Note, that the DDRH parameterization has solely been based on the DBHF self-energies at the corresponding Fermi momenta for a proton fraction Y p = 0.4. No attempt has been made to account for the momentum and further asymmetry dependence of the DBHF self-energies, which enter the calculated energies for symmetric nuclear matter and neutron matter. Both approaches show a nonlinear convergence to zero when the density decreases. which is in qualitative agreement with a study of low-density nuclear matter [23] based on a virial expansion that includes protons, neutrons, and α-particle degrees of freedom. In particular at densities below ρ < 0.1fm −3 , the DDRH model presented here is in better agreement with the DBHF result than those models constructed in Ref. [24].
From table II, one can see that the saturation density ρ sat is shifted to lower densities and the binding energy E/A of the saturation point is weaker in the DDRH model compared to the original DBHF results. The compression modulus K is somewhat smaller and symmetry energy a s is in a fairly good agreement with the DBHF results.   [14,17]. The quantities listed include the saturation density ρsat, the binding energy E/A, the compressibility modulus K, and the symmetry energy as. The symmetry energy as has been obtained assuming a quadratic dependence of the EOS on the asymmetry parameter β = (ρn − ρp)/(ρn + ρp).

IV. SPINODAL INSTABILITIES
The stability conditions for isospin asymmetric nuclear matter at finite temperatures are obtained from the free energy density F . This free energy density is defined as where E is the energy density and S denotes the entropy density. The energy density E is obtained from the energymomentum tensor where f i (k) is the thermal occupation probability. The rearrangement term Σ (r) does not appear in Eq. (20), since rearrangement does not contribute to the energy density. However, note that the rearrangement does modify other quantities such as the pressure. Furthermore, the entropy density S is given by The necessary ingredients to analyze the stability of isospin asymmetric nuclear matter against separation into two phases are now introduced. The nuclear system will be stable if the free energy of a single phase is lower than the free energy in any two-phase configurations. This stability condition implies that the free energy density is a convex function of the neutron and proton densities ρ n and ρ p . This means that the curvature matrix, is positive-definite. The requirement that the local curvature is positive demands that both the trace and the determinant are positive, The [F ij ] is a 2 * 2 symmetric matrix, since one considers a two-fluids system. Therefore, the stability matrix [F ij ] has two real eigenvalues λ ± : associated to eigenvectors δρ ± given by (i = j) Therefore, the requirements in Eq. (23) can equivalently be written as This means that both eigenvalues should be positive to guarantee the stability of the system. Due to the essentially quadratic dependence of the symmetry energy on the asymmetry parameter and a positive increasing symmetry energy with the total density, it turns out that only one eigenvalue can become negative [1]. In principle λ + will always be positive in isospin asymmetric nuclear. Only λ − can become negative and the eigenvector associated with this negative eigenvalue indicates the direction of the instability. Therefore, if Eq. (26) is violated, or equivalently λ − is negative for isospin asymmetric nuclear matter, the system will be in the unstable region of a phase transition.

V. RESULTS AND DISCUSSION
In the present work, we have investigated the spinodal region in the case of isospin asymmetric nuclear matter at finite temperature. The projection of the spinodal contour in the neutron-proton density plane for several temperatures ranging from T=0 to T=10 MeV is depicted in Fig. 3. Keep in mind that the isospin symmetry of the nucleon-nucleon interaction guarantees the symmetry of the spinodal contours with respect to the interchange of the neutron density ρ n and the proton density ρ p .
Let us first consider the results for zero temperature (T= 0 MeV). One finds that the spinodal region extends up to a maximum density of 0.127 fm −3 . This is considerably larger than the maximum density obtained in phenomenological models [1-4, 7, 9-13], which indicate that the spinodal region is situated at densities below ρ B < 0.1 fm −3 . The large spinodal region, however, is in agreement with microscopic calculations [5,6]. This difference between the phenomenological and the microscopic based models is of course related to the larger saturation densities, which are obtained in the microscopic calculations. Note, however, that the present study yields a larger maximal density than the non-relativistic BHF calculations of Vidana and Polls [6], although the saturation density for the present approach (ρ 0 = 0.176 fm −3 ) is smaller than the one obtained in the BHF calculation (ρ 0 = 0.182 fm −3 ).
The shape of the spinodal region for T=0 projected in the neutron-proton density plane, is similar to the shapes obtained in calculations before. In detail, however, there is a difference: The region of instability is more confined to isospin symmetric matter.
Also the dependence of the region of instability on temperature in Fig. 3 is very similar to the one observed in previous calculations [1-4, 6, 7, 9-13]. The spinodal region shrinks with increasing temperature, until the critical temperature T c is reached. In the present approach no spinodal region can be obtained above T = 11 MeV. Therefore, the critical temperature is considerably lower than the value T c = 17.5 MeV obtained in the other microscopic based calculation [6] or the values for T c reported in the phenomenological studies [4]. It is worth mentioning that the spinodal region obtained here at temperatures close to T c occurs at slightly higher densities than in the other studies.
It has been shown that the spinodal instabilities in isospin asymmetric nuclear appear as a mixture of isoscalar and isovector instabilities [1,2]. Isoscalar or mechanical instabilities are related to density fluctuations. Pure isoscalar move in phase δρ − p /δρ − n > 0, or of the isovector nature, if they move out of phase δρ − p /δρ − n < 0. In Fig. 4 the results obtained from our model for this ratio δρ − p /δρ − n is plotted as a function of the asymmetry parameter β = (ρ n − ρ p )/(ρ n + ρ p ) for several temperatures ranging from T = 0 up to T = 8 MeV and for two different densities. It can be observed from Fig. 4 that the ratio δρ − p /δρ − n is positive in all considered cases. This result implies that the spinodal instability is dominated by mechanical or isoscalar density fluctuations, which is in agreement with previous investigations [1,2,5,6]. Since δρ − p /δρ − n is larger than the ratio ρ p /ρ n presented by the black dotted line in Fig. 4, the dense liquid phase of the system is driven towards a more isospin symmetric composition and the gas phase will be more isospin asymmetric. These effects are less strong for higher temperatures. The liquid phase is in general also less strongly driven towards a isospin symmetric composition, when the density ρ B or the asymmetry parameter β increase.
Up to this point we have determined the spinodal zone from the evaluation of the curvature matrix [F ij ] in (22) for the free energy of the homogeneous system of nuclear matter. Having identified the regions of instability we will now try to verify these instabilities and determine the density profile in this area.
For that purpose we consider the very same DDRH and perform variational calculations in a cubic Wigner-Seitz (WS) cell. The size and shape of the WS cell should be determined in a variational manner, however, in a kind of explorative study we have chosen the length of the WS cell to be 20 fm in each direction with the origin of the coordinate system in the center of the box. Also we have restricted the study to isospin symmetric systems. In order to be consistent with our infinite nuclear matter calculations, the Coulomb energy is not included in the WS cell calculations. Therefore the proton density distributions are identical to the corresponding neutron ones displayed in Fig. 5.
The transition from isolated nuclei to a phase of homogeneous baryon matter can be observed in the various panels Fig. 5 for zero temperature. The left top panel shows the density distribution at a density of ρ = 0.005 fm −3 , i.e. only 20 protons and 20 neutrons ( 40 Ca) are present in the cubic WS cell. At these low densities, isolated nuclei will be present. The density profiles are identical in all three Cartesian directions, which corresponds to crystal of spherical nuclei located at the centers of the WS cell.
However, at higher densities this is not the case anymore. At a density of ρ = 0.05 fm −3 , the density profile in z-direction deviates a little from that in the x and y direction. Furthermore, one has a high-density region at the center of the cubic WS cell and a low-density region along the connection lines to the nearest neighbor at the boundary of the box.
However, the density distribution is not as simple as suggested by Fig. 5. This can be seen from the complete density distribution in the z = 0 plane as displayed in Fig. 6. Apart from the low density regions at the boundary of the box along the connection lines between the centers of the WS cells, one has also regions in which the density is almost zero. Therefore, the geometrical structure at this density can be described by large quasi-nuclei located in the centers of the WS cell, which are connected to each other by low-density bridges. At even higher densities, i.e. around ρ = 0.10 fm −3 , the opposite density distribution can be observed, i.e isolated low-density regions completely surrounded by high-density matter. At a density of ρ = 0.125 fm −3 , so just a little bit below the critical density (ρ c = 0.127 fm −3 ) predicted from the study of the curvature matrix in homogeneous one obtains a grid structure with small fluctuations, which resembles homogeneous nuclear matter.
Therefore, within the numerical accuracy we regard this result to be in agreement with the result of the infinite nuclear matter calculation in Fig. 3.
Performing DDRH calculations in the same kind of WS cell at the finite temperature T = 8 MeV, we obtain the neutron density distributions displayed in Fig. 7. It is clear that finite temperature effects tend to dissolve the geometrical structures observed for T = 0 MeV. The structures at T = 8 MeV are less pronounced than those for T = 0 MeV in Fig. 6 and the transition density is lower at T = 8 MeV. These results are in agreement with the results presented in Fig. 3 of the infinite nuclear matter calculations.

VI. SUMMARY AND CONCLUSION
In the present work, we construct a Density Dependent Relativistic Mean Field model (DDRH) to parameterize the results of microscopic Dirac Brueckner Hartree Fock (DBHF) calculations [17] using the realistic Bonn A potential [15]. Special attention is devoted to reproduce the DBHF self-energy in particular at densities below the saturation density. The equation of state at zero and finite temperature can then be determined in the framework of this DDRH theory for homogeneous matter as well as inhomogeneous structures. The spinodal instabilities of homogeneous matter including their temperature effects are investigated within this model.
Our investigations confirm the main results of earlier studies based on phenomenological [1,2] as well as realistic models of the NN interaction [6], that the spinodal instability is dominated by mechanical or isoscalar density fluctuations which lead to a liquid-gas phase separation with restoration of the isospin symmetry in the liquid phase. However, this restoration is less strong for higher temperatures.
Furthermore, the DDRH model predicts some details, which are different from earlier studies: A large critical density at zero temperature, a smaller critical temperature T c and details in the shape of the spinodal region. At zero temperature the spinodal region reaches densities above 0.12 fm −3 . This is different from the predictions of phenomenological calculations [1-4, 7, 9-13], but in line with the microscopic studies of Vidana and Polls [6]. It can to some extent be attributed to the differences in the predicted saturation point of nuclear matter. The present DDRH approach yields a smaller critical temperature (T c less than 11 MeV) and higher densities for the spinodal region at temperatures just below T c . The predicted spinodal region obtained from our calculation is more strongly concentrated around isospin symmetric nuclear matter.
The studies of the spinodal regions which are based on the curvature of the free energy of homogeneous matter are supplemented by variational calculations allowing for density fluctuations explicitly. Calculations employing the same DDRH model are performed in a cubic Wigner Seitz cell by minimizing the total energy to obtain density distributions in isospin symmetric nuclear matter. At a density of 0.125 fm −3 , a geometrical structure resembling homogeneous nuclear matter can be found. At lower densities, non-homogeneous nuclear matter structures occur. Finite temperature effects tend to dissolve the geometrical structures observed, i.e. the structures at finite temperature are less pronounced and the transition density to homogeneous nuclear matter is lower.