Properties of spherical and deformed nuclei using regularized pseudopotentials in nuclear DFT

We developed new parameterizations of local regularized finite-range pseudopotentials up to next-to-next-to-next-to-leading order (N3LO), used as generators of nuclear density functionals. When supplemented with zero-range spin-orbit and density-dependent terms, they provide a correct single-reference description of binding energies and radii of spherical and deformed nuclei. We compared the obtained results to experimental data and discussed benchmarks against the standard well-established Gogny D1S functional.


Introduction
The nuclear density functional theory (DFT) offers one of the most flexible frameworks to microscopically describe structure of atomic nuclei [1,2]. A key element in the nuclear DFT is the energy density functional (EDF), which is usually obtained by employing effective forces as its generators. A long-standing goal of nuclear DFT is to construct an EDF with high precision of describing existing data and high predictive power.
The Skyrme and Gogny EDFs [1,3] are the most utilized non-relativistic EDFs in nuclear structure calculations. The Skyrme EDF is based on a zero-range generator, combined with a momentum expansion up to second order, whereas the Gogny EDF is based on the generator constructed with two Gaussian terms. While zero-range potentials are computationally simpler and less demanding, they lack in flexibility of their exchange terms. In addition, in the pairing channel they manifest the wellknown problem of nonconvergent pairing energy, which needs to be regularized, see Refs. [4,5] and references cited therein. While Skyrme-type EDFs can reproduce various nuclear bulk properties relatively well, their limits have been reached [6], and proposed extensions of zero-range generators [7,8] did not prove efficient enough [9]. New approaches are, therefore, required.
To improve present EDFs, a possible route is to use EDFs based on regularized finite-range pseudopotentials [10]. Such EDFs stem from a momentum expansion around a finite-range regulator and thus have a form compatible with powerful effective-theory methods [11,12]. Here, as well as in our earlier studies [13,14], we chose a Gaussian regulator, which offers numerically simple treatment, particularly when combined with the harmonic oscillator basis. The momentum expansion can be built order-by-order, resulting in an EDF with increasing precision. Due to its finite-range nature, treatment of the pairing channel does not require any particular regularization or renormalization.
The ultimate goal of building EDFs based on regularized finite-range pseudopotentials is to apply them to beyond mean-field multi-reference calculations. However, before that, to evaluate expected performance and detect possible pitfalls, their predictive power should be benchmarked at the single-reference level. The goal of this work is to adjust the single-reference parameters of pseudopotentials up to next-to-next-tonext-to-leading order (N 3 LO) and to compare the obtained results to experimental data on the one hand and to those obtained for the Gogny D1S EDF [15] on the other. The D1S EDF offers an excellent reference to compare to, because it contains finite-range terms of a similar nature, although its possible extensions to more than two Gaussians [16], cannot be cast in the form of an effective-theory expansion. Because EDFs adjusted in present work are intended to be used solely at the single-reference level, they include a density-dependent term. This term significantly improves infinite nuclear matter properties, with the drawback that such EDFs become unsuitable for multi-reference calculations, see, e.g., Refs. [17,18].
This article is organized as follows. In Sec. 2, we briefly recall the formalism of the regularized finite-range pseudopotential and in Sec. 3 we present details of adjusting its parameters. Then, in Secs. 4 and 5, we present results and conclusions of our study, respectively. In Appendices A-D, we give specific details of our approach and in the supplemental material (https://arxiv.org/e-print/2003.10990) we collected files with numerical results given in a machine readable format.

Pseudopotential
In this study, we use the local regularized pseudopotential with terms at nth order introduced in [13], where the Gaussian regulator is defined as and1 σ and1 τ are respectively the identity operators in spin and isospin space andP σ andP τ the spin and isospin exchange operators. Standard relative-momentum operators are defined as k ij = 1 2i (∇ i − ∇ j ) and relative positions as r ij = r i − r j . Up to the nth order (N n LO), this pseudopotential depends on the following parameters, • 8 parameters up to the next-to-leading-order (NLO): 1 , H 1 , H 1 and M (6) 1 .
In the present study, we determined coupling constants of pseudopotentials that are meant to be used at the single-reference level. Therefore, we complemented pseudopotentials (1) with standard zero-range spin-orbit and density-dependent terms, which carry two additional parameters W SO and t 3 . The density-dependent term, which has the same form as in the Gogny D1S interaction [15], represents a convenient way to adjust the nucleon effective mass in infinite nuclear matter to any reasonable value in the interval 0.70 m * m 0.90 [19]. This term contributes neither to the binding of the neutron matter nor to nuclear pairing in time-even systems. To avoid using a zero-range term in the pairing channel, we neglect contribution of the spin-orbit term to pairing.
Expressions giving the contributions to the EDFs of the local regularized pseudopotential (1) can be found in Ref. [14], whereas those of the zero-range spinorbit (3) and density-dependent term (4) can be found, for example, in Refs. [20,21].

Adjustments of parameters
As explained in Sec. 2, pseudopotentials considered here contain 10 parameters at NLO, 14 at N 2 LO, and 18 at N 3 LO. In this study, we adjusted 15 series of parameters with effective masses m * /m equal to 0.70, 0.75, 0.80, 0.85, and 0.90 at NLO, N 2 LO, and N 3 LO. For each series, the range a of the regulator was varied between 0.8 and 1.6 fm.
Our previous experience shows that the use of a penalty function only containing data on finite nuclei is not sufficient to efficiently constrain parameters of pseudopotentials, or even to constrain them at all. Typical reasons for these difficulties are (i) appearance of finite-size instabilities, (ii) phase transitions to unphysical states (for example those with very large vector pairing), or (iii) numerical problems due to compensations of large coupling constants with opposite signs. To avoid these unwanted situations, the penalty function must contain specially designed empirical constrains. Before performing actual fits, such constrains cannot be easily defined; therefore, to design the final penalty function, we went through the steps summarized below.
• Step 1: We made some preliminary fits so as to detect possible pitfalls and devise ways to avoid them. The main resulting observation was that it seems to be very difficult, if possible at all, to adjust parameters leading to a value of the slope of the symmetry energy coefficient L in the range of the commonly accepted values, which is roughly between 40 and 70 MeV [22,23,24]. Therefore, for all adjustments performed in this study, we set its value to L = 15 MeV. This value is rather low, although it is at a similar lower side as those corresponding to various Gogny parameterizations: L = 18.4 MeV for D1 [25], L = 22.4 MeV for D1S [15], L = 24.8 MeV for D1M [26], or L = 43.2 MeV for D1M* [27].
• Step 2: With the fixed value of L = 15 MeV, we performed a series of exploratory adjustments with fixed values of other infinite-nuclear-matter properties, that is, for the saturation density of ρ sat = 0.16 fm −3 , binding energy per nucleon in infinite symmetric matter of E/A = −16 MeV, compression modulus of K ∞ = 230 MeV, and symmetry energy coefficient of J = 32 MeV. These initial values were the same as for the Skyrme interactions of the SLy family [28,29]. The conclusion drawn from this step was that the favoured values for ρ sat and J were slightly lower than the initial ones. Therefore, we decided to fix ρ sat and J at the results corresponding to pseudopotentials giving the lowest values of the penalty function χ 2 , see Fig. 1 and Table 1.
In a consistent effective theory, with increasing order of expansion, the dependence of observables on the range a of the regulator should become weaker and weaker. In our previous work [14], where all terms of the pseudopotential were regulated with the same range, such a behaviour was clearly visible. In the present work, the regulated part of the pseudopotential is combined with two zero-range terms. As a result, even at N 3 LO, there remains a significant dependence of the penalty functions on a, see Fig. 1. Therefore, in step 3 we picked for further analyses the parameterizations of pseudopotentials that correspond to the minimum values of penalty functions. Then, for each of the five values of the effective mass and for each of the three orders of expansion, we optimized the corresponding parameters of the pseudopotential, but this time with the infinite-matter properties not rigidly fixed but allowed to change within small tolerance intervals, see Table 1.
In the supplemental material (https://arxiv.org/e-print/2003.10990), the corresponding 15 sets of parameters are listed in a machine readable format. Following the naming convention adopted in Ref. [30], these final sets are named as   For brevity, in the remaining of this paper, we omit the date of the final adjustment, denoted by 190617, which otherwise is an inherent part of the name.
We are now in a position to list all contributions to the penalty function χ 2 , which come from the empirical constrains used in step 3 of the adjustment and from those corresponding to the nuclear data and pseudo-data that we used.
(i) Empirical properties of the symmetric infinite nuclear matter. These correspond to: saturation density ρ sat , binding energy per nucleon E/A, compression modulus K ∞ , isoscalar effective mass m * /m, symmetry energy coefficient J, and its slope L.
The target values and the corresponding tolerance intervals are listed in Table 1.
(ii) Potential energies per nucleon in symmetric infinite nuclear matter. We used values in four spin-isospin channels (S, T ) determined in theoretical calculations of Refs. [31,32]. Although it is not clear if these constraints have any significant impact on the observables calculated in finite nuclei, we observed that they seem to prevent the aforementioned numerical instabilities due to compensations of large coupling constants with opposite signs. Explicit formulas for the decomposition of the potential energy in the (S, T ) channels are given in Appendix A.
(iii) Energy per nucleon in infinite neutron matter. We used values calculated for potentials UV14 plus UVII (see Table III in [33]) at densities below 0.4 fm −3 with a tolerance interval of 25 %.
(iv) Energy per nucleon in polarized infinite nuclear matter. Adjustment of parameters often leads to the appearance of a bound state in symmetric polarized matter. To avoid this type of result, we used the constraint of E/A = 12.52 MeV at density 0.1 fm −3 (taken from Ref. [34]) with a large tolerance interval of 25 %.
(v) Average pairing gap in infinite nuclear matter. Our goal was to obtain a reasonable profile for the average gap in symmetric infinite nuclear matter and to avoid too frequent collapse of pairing for deformed minima (especially for protons). Therefore, we used as targets the values calculated for the D1S functional at k F = 0.4, 0.8, and 1.2 fm −1 with the tolerance intervals of 0.1 MeV.
(vii) Proton rms radii. We used values taken from Ref. [36] for 40 Ca, 48 Ca, 208 Pb, and 214 Pb with the tolerance intervals of 0.02 fm and that for 56 Ni (which is extrapolated from systematics) with the tolerance interval of 0.03 fm.
(viii) Isovector and isoscalar central densities. To avoid finite-size scalar-isovector (i.e. S = 0, T = 1) instabilities, we used isovector density at the center of 208 Pb and isoscalar density at the center of 40 Ca. A use of the linear response methodology (such as in Ref. [37] for zero-range interactions) would lead to too much timeconsuming calculations. As a proxy, we used the two empirical constraints on central densities, which are known to grow uncontrollably when the scalar-isovector instabilities develop. We used the empirical values of ρ 1 (0) < 0 fm −3 in 208 Pb and ρ 0 (0) < 0.187 fm −3 in 40 Ca with asymmetric tolerance intervals as described in Ref. [14]. For ρ 0 (0) in 40 Ca, we have used the central density obtained with SLy5 [29] as an upper limit. In the parameter adjustments performed in this study, possible instabilities in the vector channels (S = 1) are still not under control.
(ix) Surface energy coefficient. As it was recently shown [38,39], a constraint on the surface energy coefficient is an efficient way to improve properties of EDFs. For the regularized pseudopotentials considered here, we calculate a simple estimate of the (x) Coupling constants corresponding to vector pairing. Terms of the EDF that correspond to this channel are given in Eq. (36) of Ref. [14]. To avoid transitions to unphysical regions of unrealistically large vector pairing, we constrain them to be equal to 0 ± 5 MeV fm 3 .

Parameters and statistical uncertainties
For the purpose of presenting observables calculated in finite nuclei, we decided to use a criterion of binding energies of spherical nuclei, see Sec. 4.3. It then appears that optimal results are obtained for m * /m = 0.85 at N 3 LO [40] and a = 1.50 fm, that is, for the pseudopotential named REG6d. Following this guidance, below we also present some results corresponding to the same effective mass of m * /m = 0.85 and lower orders: REG2d (NLO and a = 0.85 fm) and REG4d (N 2 LO and a = 1.15 fm). For an extended comparison with the Gogny D1S parameterization [15], which corresponds to m * /m = 0.697, we also show results for m * /m = 0.70, that is, for REG6a (N 3 LO and a = 1.60 fm). Parameters of the four selected pseudopotentials are tabulated in Appendix C. In the supplemental material (https://arxiv.org/e-print/2003.10990) they are collected in a machine readable format. We performed the standard analysis of statistical uncertainties as presented in Ref. [41]. For REG2d, REG4d and REG6d, eigenvalues of the Hessian matrices corresponding to penalty functions scaled to χ 2 = 1 are shown in Fig, 2(a). The numbers of eigenvalues correspond to the numbers of parameters optimized during the adjustments, and, therefore, vary from 10 (NLO) to 18 (N 3 LO).
The magnitude of the eigenvalues of the Hessian matrices reveals how well the penalty functions are constrained in the directions of the corresponding eigenvectors in the parameter space. We observe that for the three pseudopotentials considered here, there is a rapid decrease of magnitude from the first to the third eigenvalue and then a slower and almost regular decrease, where no clear gap can be identified. This suggests that all parameters of the pseudopotentials are important.
For three tin isotopes of different nature: 100 Sn (closed-shell, isospin symmetric, unpaired), 120 Sn (open-shell, isospin asymmetric, paired) and 132 Sn (closed-shell, isospin asymmetric, unpaired), we calculated the propagated statistical uncertainties of the total binding energies as functions of the number of kept eigenvalues of the Hessian matrices, Figs 2(b)-(d) for REG2d-REG6d, respectively. For each of the considered parameterizations, after a given number of kept eigenvalues (denoted in Figs 2(b)-(d) by vertical lines), we observe a saturation of the propagated statistical uncertainties. Therefore, we performed the final determination of the statistical uncertainties by keeping these minimal numbers of eigenvalues, i.e. 6 eigenvalues for REG2d (NLO) and 7 for REG4d (N 2 LO) and REG6d (N 3 LO).

Infinite nuclear matter
In Table 2, we list quantities characterizing the properties of infinite nuclear matter. We present results for pseudopotentials REG2d, REG4d, REG6d, and REG6a compared to those characterizing the D1S interaction [15]. For the two strongly constrained quantities, m * /m and L, the target values are almost perfectly met, whereas, for the other ones, we observe some deviations, which, nevertheless, are well within the tolerance intervals allowed in the penalty function.
For pseudopotentials REG6a and REG6d, the isoscalar effective mass in symmetric matter and energies per particle (equations of state) for symmetric, neutron, polarized, and polarized neutron matter are plotted in Fig. 3 along with the same quantities for D1S [15]. The plotted equations of state can be obtained from those calculated in four Table 2. Infinite nuclear matter properties corresponding to pseudopotentials REG2d, REG4d, REG6d, and REG6a, compared to those of the Gogny D1S interaction [15].  spin-isospin (S, T ) channels, see Appendix A. For these two N 3 LO pseudopotentials, equations of state of symmetric matter are somewhat stiffer than that obtained for D1S. This is because of its slightly larger compression modulus K ∞ . We also can see that for polarized symmetric matter, a shallow bound state appears at low density. This feature also affects D1S. The constraint on the equation of state of polarized symmetric matter introduced in the penalty function has probably limited the development of this state, but did not totally avoid its appearance. Further studies are needed to analyze to what extent it could impact observables calculated in time-odd nuclei and how this possible flaw might be corrected. The two main differences that appear when we compare the properties in infinite nuclear matter of REG6a and REG6d on one hand and those of D1S on the other hand relate to the equation of state of the neutron matter and isoscalar effective masses. First, near saturation, the regularized pseudopotentials give equations of state of neutron matter slightly lower than D1S, which can be attributed to its lower symmetry energy. Second, for the N 3 LO pseudopotentials, dependence of the effective on density is less regular than for D1S. We note, however, that the N 3 LO effective masses are monotonically decreasing functions of the density, and thus the pseudopotentials obtained in this study do not lead to a surface-peaked effective mass, a feature which was expected to improve the description of the density of states around the Fermi energy [42].

Binding energies, radii, and pairing gaps of spherical nuclei
In this section, we present results of systematic calculations performed for spherical nuclei and compared with experimental data. For the purpose of such a comparison, we have selected a set of 214 nuclei that were identified as spherical in the systematic calculations performed for the D1S functional in Refs. [43,44]. In Fig. 4, we present an overview of the binding-energy residuals obtained for the D1S, REG6a, and REG6d functionals. Experimental values were taken from the 2016 atomic mass evaluation [35]. The obtained root-mean-square (RMS) binding-energy residuals are equal to 2.582 MeV for D1S, 1.717 MeV for REG6a, and 1.458 MeV for REG6d. We also see that for REG6d, the trends of binding-energy residuals along isotopic chains in heavy nuclei become much better reproduced. As a reference, we have also determined the analogous RMS value corresponding to the UNEDF0 functional [36,45], which turns out to be equal to 1.900 MeV. In Figs. 5 and 6, we show detailed values of binding-energy residuals along the isotopic or isotonic chains of semi-magic nuclei. In most chains one can see a clear improvement of the isospin dependence of masses. In particular, in almost all semimagic chains, kinks of energy residuals at doubly magic nuclei either decreased or even vanished completely, like at N = 82 and 126, see Figs. 6(b) and (c), respectively.
In Fig. 7, for the same set of EDFs and nuclei as those used in Fig. 4, we show the analogous residuals of the charge radii of spherical nuclei. The experimental values were taken from Ref. [46]. Again, the N 3 LO EDFs provide the smallest deviations from data. We note that the residuals of the order of 0.02 fm are typical for many Skyrme-like EDFs, for example, for the UNEDF family of EDFs [6]. Figures 8(a) and (b) present summary of the RMS residuals of binding energies and charge radii, respectively, which were obtained in this study. We see that a decrease of the penalty functions when going from NLO to N 2 LO, see Fig. 1, is often accompanied by an increase of the RMS residuals. This indicates that the data for 17 spherical nuclei, which are included in

Single particle energies
In Figs. 11 and 12, we show comparison of single-particle energies calculated in semimagic nuclei for the D1S [15], REG6a (m * /m = 0.70), and REG6d (m * /m = 0.85) functionals with the empirical values taken from the compilation published in the supplemental material of Ref. [47], which contains the single-particle energies collected within three data sets. In all panels of Figs represent spreads between the minimum and maximum values. Quantum numbers in parentheses indicate single-particle states with corresponding attributed spectroscopic factors smaller than 0.8 or unknown. The spin-orbit interaction corresponding to functional REG6a (m * /m = 0.70) is smaller than that of D1S, which may explain differences between the single-particle energies of states with large orbital angular momenta. Differences between the results   obtained for functionals with m * /m = 0.70 and m * /m = 0.85 mostly amount to a global compression. Generally speaking, the calculated single-particle energy spacings are larger than those of the empirical ones, which is typical for the effective masses being smaller than one. We note that the comparison between the calculated and empirical single-particle energies is given here only for the purpose of illustration. Indeed, both are subjected to uncertainties of definition and meaning. The calculated ones, which are here determined as the eigenenergies of the mean-field Hamiltonian, could also be evaluated from calculated odd-even mass differences. Similarly, determination of the empirical ones is always uncertain from the point of view of the fragmentation of the single-particle strengths. For these reasons, we did not include single-particle energies in the definition of our penalty function, see Sec. 3. Nevertheless, positions and ordering of single-particle energies are crucial for a correct description of other observables, such as, for example, ground-state deformations or deformation energies. Therefore, we consider comparisons presented in Figs. 11 and 12 to be very useful illustrations of properties of the underlying EDFs.

Deformed nuclei
Using the methodology of extrapolating results calculated for N 0 = 16 harmonicoscillator shells to infinite N 0 , presented in Appendix D, for a set of nine nuclei from 54 Cr to 252 Cf we determined deformation energies, Fig. 13, and binding-energy deformed Figure 14. Same as in Fig. 13 but for the binding-energy residuals calculated at spherical shapes (a) and deformed minima (b).
residuals, Fig. 14. The figures compare results obtained for the D1S [15] and N 3 LO REG6d functionals. For N 3 LO REG6d, in Fig. 14(a) we also show propagated uncertainties [41] of spherical energies determined using the covariance matrix available in the supplemental material (https://arxiv.org/e-print/2003.10990). For illustration, the same propagated uncertainties are also plotted in Fig. 14(b), whereas determination of full propagated uncertainties of deformed minima is left to a future analysis of deformed solutions. In general, the pattern of deformations obtained for both functionals is very similar. This is gratifying, because deformed nuclei were not included in the adjustment of either one of the two EDFs. For this admittedly fairly limited sample of nuclei, the pattern of RMS binding-energy residuals is fairly analogous to what we have observed in spherical nuclei, see Sec. 4.3, with REG6d giving values that are about 30% smaller than those for D1S. It is interesting to see that in several instances, the two functionals generate absolute energies of the deformed minima that are more similar to one another than those of the spherical shapes. In the present study we limit ourselves to presenting results only for these very few nuclei, whereas attempts of using deformed nuclei in penalty functions [48] and systematic mass-table calculations [49] are left to forthcoming publications.

Conclusions
In this article, we reported on the next step in adjusting parameters of regularized finiterange functional generators to data. We have shown that an order-by-order improvement of agreement with data is possible, and that the sixth order (N 3 LO) functional describes data similarly or better than the standard Gogny or Skyrme functionals.
We implemented adjustments of parameters based on minimizing fairly complicated penalty function. Our experience shows that a blind optimization to selected experimental data seldom works. Instead, one has to implement sophisticated constraints, which prevent wandering of parameters towards regions where different kinds of instabilities loom.
We consider the process of developing new functionals and adjusting their parameters a continuous effort to better their precision and predictive power. At the expense of introducing single density-dependent generator, here we were able to raise the values of the effective mass, obtained in our previous study [14], well above those that are achievable with purely two-body density-independent generators [19]. Such a solution is perfectly satisfactory at the single-reference level. However, for multireference implementations, the density-dependent term must be replaced by secondorder three-body zero-range generators [50], or otherwise entirely new yet unknown approach would be required.
Although a definitive conclusion about usefulness of EDFs obtained in this study can only be drawn after a comparison of observables of more diverse nature, this class of pseudopotentials looks promising, even if it can clearly be further improved. In the future, we plan on continuing novel developments by implementing non-local regularized pseudopotentials along with their spin-orbit and tensor terms [13]. This may allow us to fine-tune spectroscopic properties of functionals and facilitate precise description of deformed and odd nuclei.
Appendix A. Decomposition of the potential energy in (S, T ) channels The techniques to derive decomposition of the potential energy per particle, E (S,T ) pot /A, into four spin-isospin (S, T ) channels are the same for finite-range and zero-range pseudopotentials. Therefore, we do not repeat here the details of the derivation, which can be found, for example, in Ref. [51].
First, we recall the expression for the auxiliary function F 0 (ξ), already introduced in Ref. [30], Then, in the symmetric infinite nuclear matter with Fermi momentum k F and density ρ 0 = 2k 3 F /3π 2 , contributions of the finite-range local pseudopotential (1) at order zero (n = 0) to (S, T ) channels can be expressed as: 5) and those at higher orders n as: Appendix B. Estimate of the surface energy coefficient In section 3, we introduced a constraint on the estimate of the surface energy coefficient a LDM surf , calculated with a liquid-drop-type formula. In the case of local functionals (such as Skyrme functionals), to calculate the surface energy coefficient [38], several approaches can be considered, such as the Hartree-Fock (HF) calculation [52], approximation of the Extended Thomas Fermi (ETF) type [53] or Modified Thomas Fermi (MTF) type [54], or within a leptodermous protocol, which is based on an analysis of calculations performed for very large fictitious nuclei [55]. Some of these approaches are not usable with the regularized pseudopotentials considered in this article. Indeed, the ETF and MTF methods can only be used for functionals that depend on local densities. In principle, the leptodermous protocol could be used, but it would require a significant expense in CPU time. Moreover, the HF calculations cannot be considered because the Friedel oscillations of the density make the extraction of a stable and precise value of the surface energy coefficient very difficult (see discussion in Ref. [38] and references therein).
Therefore, for the purpose of performing parameter adjustments, we decided to use a very simple estimate of the surface energy coefficient, which is usable with any kind of functional. After determining the self-consistent total binding energy E of a fictitious symmetric, spin-saturated, and unpaired N = Z = 40 nucleus without Coulomb interaction, we used a simple liquid-drop formula to calculate the surface energy coefficient, where A = 80 and a v is the volume energy coefficient in symmetric infinite nuclear matter at the saturation point.
Values of a LDM surf obtained in this way do depends on A, but, at least in the case Skyrme functionals, they vary linearly with the surface energy coefficients obtained using full HF calculations. Detailed study of the usability of estimates (B.1) will be the subject of a forthcoming publication [56].
In section 3, we used the value of a LDM surf = 18.5 MeV as the target value of the parameter adjustments. This value is only slightly below the value obtained for the Skyrme functional SLy5s1 (18.6 MeV), which is an improved version of the SLy5 functional with optimized surface properties [38,39]. This target value we used was only an educated guess, and it may require fine-tuning after a systematic study of the properties of deformed nuclei will have been performed.

Appendix C. Parameters of the pseudopotentials
Parameters of the pseudopotentials used in Sec. 4, REG2d at NLO, REG4d at N 2 LO, and REG6d at N 3 LO with m * /m ≃ 0.85, and REG6a at N 3 LO with m * /m ≃ 0.70 are reported in Table C1 along with their statistical uncertainties. As it turns out, values of parameters rounded to the significant digits, which would be consistent with the statistical uncertainties, give results visibly different than those corresponding to unrounded values. Therefore, in the Table we give all parameters up to the sixth decimal figure. Moreover, the statistical uncertainties of parameters are only given for illustration, whereas the propagated uncertainties of observables have to evaluated using the full covariance matrices [41]. Parameters of other pseudopotentials derived in this study along with the covariance matrix corresponding to REG6d are listed in the supplemental material (https://arxiv.org/e-print/2003.10990).

Appendix D. Extrapolation of binding energies of deformed nuclei to infinite harmonic-oscillator basis
In this study, results for spherical nuclei were obtained using code finres 4 [57], which solves HFB equations for finite-range generators on a mesh of points in spherical space coordinates. Because of the spherical symmetry, it is perfectly possible to perform calculations with a mesh dense enough and a number of partial waves high enough for the results to be stable with respect to any change of the numerical conditions. Results for deformed nuclei were obtained using the 3D code hfodd (v2.92a) [58,59] or axialsymmetry code hfbtemp [60]. These two codes solve HFB equations by expanding single-particle wave functions on harmonic-oscillator bases. Since for deformed nuclei the amount of CPU time and memory is much larger than for spherical ones, it is not practically possible to use enough major harmonic-oscillator shells to reach the asymptotic regime, especially for heavy nuclei.
In order to estimate what would be the converged asymptotic value of the total binding energy of a given deformed nucleus, we proceeded in the following way. First, using code finres 4 , we determined the total binding energy E sph of a given nucleus at the spherical point. Second, using codes hfodd or hfbtemp, for the same nucleus and for a given number of shells N 0 , we determined total binding energies E sph (N 0 ) (constrained to sphericity) and E def (N 0 ) (constrained to a non-zero deformation). Third, we assumed that with N 0 → ∞, the deformation energy ∆E def (N 0 ) = E def (N 0 ) −E sph (N 0 ) converges much faster to its asymptotic value than either of the two energies. Fourth, within this assumption, we estimate the asymptotic energies of deformed nuclei as E def (N 0 = ∞) = E sph + ∆E def (N 0 ) = E def (N 0 ) + E sph − E sph (N 0 ). (D.1) In Fig. D1, we present typical convergence pattern that supports the main assumption leading to estimate (D.1). Fig. D1(a) shows deformation energies ∆E of 166 Er calculated for the numbers of shells between N 0 = 10 and 16 using the D1S functional [15]. It is clear that in the scale of the figure, differences between the four curves are hardly discernible. In a magnified scale, in Fig. D1(b) we show total energies E def (N 0 ) − E def (16) relative to that obtained for N 0 = 16. We see that at N 0 = 10 and 12, the relative energies are fairly flat; however, they also exhibit significant changes, including sudden jumps related to individual orbitals entering and leaving the space of harmonic-oscillator wave functions that are included in the basis with changing deformation. Nevertheless, already at N 0 = 14, the relative energy becomes very smooth and almost perfectly flat. This behaviour gives strong support to applying estimate (D.1) already at N 0 = 16. Such a method was indeed used to present all total binding energies of deformed nuclei discussed in this article.
Finally, in Fig. D1(c), in the absolute energy scale we show total binding energies in 166 Er obtained using codes hfodd (up to N 0 = 20, large full circles) and hfbtemp (up to N 0 = 30, small empty circles). Calculations were constrained to the spherical point and thus the results are directly comparable with the value of E sph = −1339.069768 obtained using the spherical code finres 4 (horizontal line). These results constitute a very strong benchmark of our implementations of the N 3 LO pseudopotentials in three very different codes. For N 0 ≤ 20, differences between the hfodd and hfbtemp total energies (not visible in the scale of the figure) do not exceed 3 keV. By fitting an exponential curve to the hfbtemp results (thin line) we obtained the extrapolated asymptotic value of energy E sph (N 0 = ∞) = −1339.097 (34), which within the extrapolation error of 34 keV (shown in the figure by the shaded band) perfectly agrees with the finres 4 value.