Spurious finite-size instabilities with Gogny-type interactions

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. 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 properties of nuclei, like neutron and proton densities, if this D1M* force is used. The conclusions from this study are then, for the first time, tested against mean-field calculations in a coordinate representation for several nuclei. Unphysical results for several 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 led to a new parameter 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. This latter value makes D1M * the a e-mail: martini.marco@gmail.com b e-mail: depace@to.infn.it c e-mail: bennaceur@ipnl.in2p3.fr first Gogny-type interaction able to predict neutron star masses in agreement with the already mentioned recent observations.
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 finite-size instabilities can easily be revealed with mean-field calculations on a mesh [14][15][16], a series of such HFB calculations is then performed to possibly confirm their 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, R (S,T ) (q, ω = 0), exhibits a pole. In the left panel of 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 finite-size instabilities established in ref. [11] is satisfied for the new D1M * force. The behavior of the critical 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 the right panel of fig. 1.
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 finite-size instabilities in nuclei calculated with the D1M * force and to remove the ambiguity concerning D1N we have performed a series of mean-field calculations using the recently developed code FINRES 4 [16], which solves the selfconsistent Hartree-Fock-Bogolyubov equations on a mesh for finite-range interactions using the method described in [17]. This series of calculations is done for three very different nuclei: the doubly-magic stable nucleus 208 Pb, the semi-magic nucleus 120 Sn for which pairing is active, and the very light symmetric nucleus 4 He. This latter is not necessarily supposed to be correctly described by a meanfield 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 * . In fig. 2 the proton and neutron densities from an Hartree-Fock calculation for 208 Pb with the D1M, D1N and D1M * Gogny interactions are represented. With the D1M * interaction, the calculation did not converge and led to oscillations of the isovector density with very large amplitude. The proton and neutron densities for different numbers of HFB iterations are represented in 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 large difference is due to significant oscillations of the densities in the bulk of the nucleus and 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 in fig. 3. 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 in fig. 2, 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).
In the right panel of fig. 1, 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 transferred 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 in fig. 4. As for 208 Pb and 120 Sn, the interaction D1M leads to a convergent result. As expected 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 not all nuclei lead to finite-size instabilities with D1M * . For example, we obtained convergent results with D1M * for 16 O, 100 Sn and 176 Sn (and with all 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 * (or D1N) interaction should only be used with the basis employed to fit its parameters. This, of course, prevents its use for calculations on a mesh, but also leads to results that depend on the size of the basis when the equations are solved by expanding the solutions on harmonic oscillator wave functions. This is exemplified in fig. 5, where the binding energy E(N sh ) for 48 Ca is plotted as a function of the number of shells N sh for the interactions D1S and D1M * using the code HFBTHO [18]. On can clearly see that E(N sh ) rapidly reaches a plateau with D1S while, with D1M * , its change is significant for each increase of N sh up to N sh = 24, the last value with which we were able to obtain a converged results using HFBTHO.
Second, 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 a robust correspondence exists between the finite-size instabilities which can be inferred 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 instabilities. 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 [19][20][21].

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: All data generated during this study are contained in this published article.] Publisher's Note The EPJ Publishers remain neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.