Spurious finite-size instabilities of a new Gogny interaction suitable for astrophysical applications

Recently, a new parameterization of the Gogny interaction suitable for astrophysical applications, named D1M*, has been presented. We investigate the possible existence of spurious finite-size instabilities of this new Gogny force by repeating a study that we have already performed for the most commonly used parameterizations (D1, D1S, D1N, D1M) of the Gogny force. This study is based on a fully-antisymmetrized random phase approximation (RPA) calculation of the nuclear matter response functions employing the continued fraction technique. The conclusions from this study are then tested against mean-field calculations in a coordinate representation for several nuclei. It turns out that this new Gogny interaction is affected by spurious finite-size instabilities in the scalar isovector channel; hence, unphysical results are expected in the calculation of bulk properties, like neutron and proton densities, of some nuclei if this D1M* force is used. Unphysical results for some nuclei are also obtained with the D1N parameterization of the Gogny force. These observations strongly advocate for the use of the linear response formalism to detect and avoid finite-size instabilities during the fit of the parameters of Gogny interactions as it is already done for some Skyrme forces.

Recently, a new parameterization of the Gogny interaction suitable for astrophysical applications has been presented [1] since it has been found [2,3] that the most successful parameterizations of this force for describing finite nuclei, namely D1S [4], D1N [5] and D1M [6], commonly suffers from a rather soft neutron matter equation of state and, as a consequence, are unable to reach neutron star masses of about 2M ⊙ , as required by recent astrophysical observations [7,8].
The authors of Ref. [1] proposed a reparameterization scheme that preserves the main properties of the Gogny force but allows for tuning the density dependence of the symmetry energy which, in turn, modifies the predictions for the maximum stellar mass. This scheme works well with D1M as a starting point, and leads to a new parameter Email addresses: martini.marco@gmail.com (M. Martini), depace@to.infn.it (A. De Pace), bennaceur@ipnl.in2p3.fr (K. Bennaceur) set, dubbed D1M * . The parameters of this new interaction are constrained by requiring the same saturation density, energy per particle, compressibility, effective mass in symmetric nuclear matter, symmetry energy at density 0.1 fm −3 as in the original D1M force. The change in the slope is from L = 24.83 MeV, the value of D1M, to L = 43.18 MeV.
In order to preserve the good performance of D1M in describing nuclear structure features of finite nuclei the authors of Ref. [1] checked that the basic bulk properties of D1M * , such as binding energies and charge radii of even-even nuclei, remain globally unaltered as compared to D1M. For this purpose, they carried out Hartree-Fock-Bogolyubov (HFB) calculations for 620 even-even nuclei of the 2012 AME [9] database using the HFBaxial code [10] which solves the HFB equations in a harmonic oscillator basis.
Beyond these evaluations performed by the authors, it is interesting to check the behavior of this newly adjusted interaction with respect to the development of finite-size instabilities, since it was pointed out that they can be hidden by the use of a representation of the quasiparticle wave functions on a limited number of harmonic oscillator shells [11] whereas these instabilities can develop when the calculation is performed on a mesh.
To study the possible onset of finite-size instabilities, we start by repeating here the same study as the one published in Ref. [12] where we developed a fully-antisymmetrized random phase approximation (RPA) calculation of the nuclear matter response functions based on the continued fraction (CF) technique to investigate the possible existence of spurious finite-size instabilities of the Gogny forces. In Ref. [12] we considered the most commonly used parameterizations of the Gogny force (D1 [13], D1S, D1N, D1M), as well as recent generalizations that include tensor terms. Here we complete the study by considering D1M * . Since finitesize instabilities can easily be revealed with meanfield calculations on a mesh [14,15,16], a series of such HFB calculations is then performed to possibly confirm there appearance for some of the considered interactions.
The key quantities to investigate unphysical finite-size instabilities are the critical densities ρ c , i.e. the lowest densities at which the nuclear response calculated at zero transferred energy (ω = 0), R (S,T ) (q, ω = 0), exhibits a pole. In Fig. 1 we show the critical densities in the three spin-isospin (S, T ) channels (0, 1), (1, 1) and (1, 0) as a function of the transferred momentum q for the D1M * force. We disregard the (S, T ) = (0, 0) channel and its corresponding physical spinodal instability. All the details of the calculations are here omitted since they can be found in Ref. [12]. The only difference with respect to this article being the parameters of the considered Gogny forces. The result that we obtain in the scalar-isovector channel (S, T ) = (0, 1) is particularly important: the critical density rapidly decreases with q reaching values around ρ c ≃ 1.2ρ 0 , where ρ 0 = 0.16 fm −3 is the empirical saturation density. Furthermore the values of ρ c only slowly increase with q after the minimum. In Ref. [11], a systematic quantitative analysis of the connection between finite nuclei and nuclear matter instabilities in the (S, T ) = (0, 1) channel was performed finding that a functional is stable if the lowest critical density at which a pole occurs in nuclear matter calculations is larger than the typical central density often obtained for 40 Ca which is, in practice, around 1.2 times the saturation density. In addition, one also has to verify that this pole represents a distinct global minimum in the (ρ c , q) plane. It therefore appears that neither of the two stability criteria to avoid spurious finitesize instabilities established in Ref. [11] is satisfied for the new D1M * force. The behavior of the criti- cal density in the (S, T ) = (0, 1) channel for D1M * is worse with respect to the corresponding behavior for D1M and D1N, as one can observe on Fig. 2.
One of the conclusions of the analysis of Ref. [12] was that the D1N parameterization of the Gogny force should be treated with some caution since the stability criteria are not strictly respected by this force.
To possibly confirm the appearance of finitesize instabilities in nuclei calculated with the D1M* force and to remove the ambiguity concerning D1N we have performed a series of meanfield calculations using the recently developed code FINRES 4 [16] which solves the self-consistent Hartree-Fock-Bogolyubov equations on a mesh for finite-range interactions. This series of calculations is done for three very different nuclei: the doublymagic stable nucleus 208 Pb, the semi-magic nucleus 120 Sn and the very light symmetric nucleus 4 He. This latter is not necessarily supposed to be correctly described by a mean-field calculation, but it represents an interesting test case for the study of instabilities in the (S, T ) = (0, 1) channel because of its small number of constituents and its approximate proton-neutron symmetry only weakly broken by the small Coulomb field. The calculations were done in a spherical box with radius R = 20 fm on a mesh with a spacing δr = 0.1 fm and with an expansion up to ℓ max = 19 for the densities. The Coulomb exchange contribution and the center-of-mass correction were treated exactly. With the D1M interaction, the self-consistent calculations were initialized from a schematic Woods-Saxon potential. The final converged densities obtained with D1M were then used to initialize the calculations for D1N and D1M*.
On Figure 3 are represented the proton and neutron densities from an Hartree-Fock calculation with the D1M, D1N and D1M* Gogny interactions for 208 Pb. With the D1M* interaction, the calculation did not converge and lead to oscillations of the isovector density with very large amplitude. The proton and neutron densities for different numbers of HFB iterations are represented on the figure. The calculations converged with D1M and D1N, for which only the final densities are represented. We can nonetheless observe that the difference between the proton and neutron densities at r = 0 is much larger with D1N than with D1M. This can be interpreted as a warning signal for a finite-size instability close to show up.
A series of HFB calculations were also done for 120 Sn with the same interactions. The corresponding densities are represented on Figure 4. In this case, we observe that the calculation converges with D1M, but neither with D1M* nor D1N. For the latter two interactions, in the same manner as on Fig. 3, we have plotted the proton and neutron densities for different number of iterations during the calculations.
From these calculations we can conclude that D1N does not always lead to convergent results when it is used to calculate finite nuclei on a mesh. The non-convergence observed here with 120 Sn is not the only one we encountered. Similar instabilities were observed for nuclei with moderate or large neutron number and for calculations done with the HF or HFB approximations (for example 60 Ca calculated with the HF approximation).
On Fig. 2, we observe that for D1N, the critical densities in the (S, T ) = (0, 1) channel approach the saturation density for a relatively narrow interval of transfered momenta q. Since the typical wave-length of the oscillations leading to an instability is λ ∼ 2π/q, we can expect that D1N cannot lead to the appearance of instabilities in very small nuclei. The situation is rather different with D1M* which shows a plateau of critical densities at approximately 1.2 × ρ sat for 2.8 fm −1 q 4 fm −1 . This means that with D1M*, finite-size instabilities may well appear in relatively small nuclei.
To test this conjecture, we have performed a series of calculations for 4 He. The results are presented on Fig. 5. As for 208 Pb and 120 Sn, the interaction D1M leads to a convergent result. As ex- pected from the aforementioned argument, the calculation did converge with D1N but not with D1M* for which, despite the very limited number of nucleons and the tiny asymmetry of the system, an instability develops. Let us mention that the only nuclei for which we obtained convergent results with D1M* were N = Z nuclei calculated without Coulomb interaction. Let us also mention that random tests with D1S for nuclei with mass number above 40 did not show any sign of instabilities.
Before we conclude, we can comment on the results obtained in Ref. [1] using the HFBaxial code. The convergence of the calculations, with results in good agreement with data, obtained by the authors for finite nuclei are owed to the use of a basis with a given number of oscillator shells. This representation strongly renormalizes the interaction and inhibits the development of instabilities and therefore makes the obtained results sound. But this renormalization has several important drawbacks. First, the D1M * interaction should only be used with the basis employed to fit its parameters. This, of course, prevents its use for calculations on a mesh. This also questions the fact that one usually uses different harmonic oscillator parameter values for different nuclei and shapes and, therefore, implicitely uses a different interaction. Second, calculations beyond mean-field, such as RPA, for finite nuclei may lead to spurious results [17] which means that this interaction should be used at the mean-field level only. Finally, the strong renormalization on the harmonic oscillator basis breaks the link with the properties of the interaction in infinite nuclear matter where the properties of the saturation point and the various equations of state are calculated without accounting for this renormalization but with the assumption that local densities are constant in space.
In this article, we have shown that it exists a robust correspondance between the finite-size instabilities which can be infered from RPA calculations in infinite nuclear matter and the one observed in nuclei for Gogny forces as it is the case for Skyrme interactions [11]. Therefore we encourage to use infinite nuclear matter RPA calculations as a tool to avoid these finite-size instabilites. This would allow to use the obtained parameters sets on any basis and would maintain the link between infinite nuclear matter and finite nuclei. Such a tool already exists [12] and could be used in the fit protocol for the construction of new Gogny forces, as already done for the Skyrme functionals [18,19,20].