Towards accurate nuclear mass tables in covariant density functional theory

Nuclear binding energies (or masses) is one of the most fundamental properties of atomic nuclei. In experiment they are studied with extreme relative mass precision of $\delta m/m=10^{-9}$ and in some nuclei even with precision of $10^{-10}$ Mittig et al. (2003); Dilling et al. (2018). In absolute terms this corresponds to the accuracy of the measurements of the binding energies better than 1 keV. Accurate knowledge of binding energies is important for nuclear structure, neutrino physics and nuclear astrophysics. The measurements of nuclear masses far from stability provide a fundamental test for our understanding of nuclear structure Mittig et al. (2003) and the constraints on model predictions for the boundaries of nuclear landscape (see Ref. Erler et al. (2012); Afanasjev et al. (2013); Agbemava et al. (2014)). The measurements of mass differences ($Q$ values) of specific isotopes play a key role in neutrino physics: in order to determine the neutrino mass on a sub-eV level one should measure the $Q$ values of certain $\beta$ transitions better than 1 eV Dilling et al. (2018). Different nuclear astrophysical processes [such as r-process (see Refs. Arnould et al. (2007); Cowan et al. (2021)) or the processes in the crust of neutron stars (see Refs. Chamel and Haensel (2008); Lau et al. (2018))] also sensitively depend on nuclear masses since they affect nuclear reaction rates.

There was a substantial and continuous effort in non-relativistic models to improve the functionals and global reproduction of experimental masses (see, for example, Refs. Tondeur et al. (2000); Samyn et al. (2002); Goriely et al. (2013); Baldo et al. (2013); Möller et al. (2016); Kortelainen et al. (2010, 2012, 2014) and references quoted therein). Existing global fits of the model parameters to the ground state masses reach an accuracy of $\approx 0.5$ MeV in the global reproduction of experimental masses in the microscopic+macroscopic (mic+mac) approach Möller et al. (2016) and in Skyrme density functional theory (DFT) Goriely et al. (2013) and around 0.8 MeV in the Gogny DFT Goriely et al. (2009). These fits take into account the correlations beyond mean field by either adding rotational and vibrational corrections in a phenomenological way or by employing five-dimensional collective Hamiltonian (5DCH). Note that the global fits at the mean field level in the Skyrme DFT lead to $1.4-1.9$ MeV accuracy in the global description of experimental binding energies (see Refs. Kortelainen et al. (2010, 2012, 2014)).

In contrast, the attempts to create covariant energy density functionals (CEDF) based on the global set of data on binding energies are very limited in covariant density functional theory (CDFT). It is only in Refs. Peña-Arteaga et al. (2016); Taninah and Afanasjev (2023) that they were undertaken (see Sec. II for more detail). This is due to the fact that the calculations within the CDFT are significantly more numerically challenging and more time-consuming as compared with those carried out in nonrelativistic DFT due to relativistic nature of the CDFT. First, the wavefunctions of the single-particle states in the CDFT are represented by Dirac spinors the large and small components of which have opposite parities. As a consequence, the basis for a given maximum value of principal quantum number $N$ in the basis set expansion has approximately double size in the CDFT as compared with non-relativistic models. As a result, the numerical calculations of matrix elements and the diagonalization of the matrices requires substantially more time than in the case of nonrelativistic DFTs. A rough estimate provided in Ref. Reinhard et al. (1986) for spherical symmetry indicates that the calculations in the Skyrme DFT and CDFT for a given maximum value of principal quantum number $N$ in the basis set expansion differ by approximately an order of magnitude. Second, the most widely used classes of CEDFs contain the mesons in addition to fermions Vretenar et al. (2005); Agbemava et al. (2014). This requires the solution of Klein-Gordon equations which adds numerical cost to the calculations. Third, in the CDFT the nucleonic potential ($\approx-50$ MeV/nucleon) emerges as the sum of very large attractive scalar $S$ ($\approx-400$ MeV/nucleon) and repulsive vector $V$ ($V\approx 350$ MeV/nucleon) potentials (see Ref. Vretenar et al. (2005)). In the nucleus with mass $A$, these values are multiplied by $A$ and this leads to a cancellation of very large quantities. This may lead to a slower convergence of binding energies as a function of the size of the basis in the CDFT as compared with non-relativistic DFTs.

Many of the fits of the relativistic and non-relativistic functionals (including global ones) have been carried out using computer codes employing basis set expansion approach in which the wave functions are expanded into the basis of harmonic oscillator wave functions. This expansion is precise for the infinite size of the basis. However, due to numerical reasons the basis is truncated in the calculations. Unfortunately, a rigorous analysis of numerical errors in calculated binding energies emerging from the truncation of basis is missing in global calculations. However, there are some indications that they are not negligible. For example, the relativistic Hartree-Bogoliubov (RHB) calculations of the ground state in ²⁰⁸Pb with DD-ME2 functional indicate that binding energies calculated with $N_{F}=20$ fermionic shells deviate from asymptotic value of binding energy $B_{\infty}$ by approximately 250 keV (see discussion of Fig. 2 in Supplemental Material to Ref. Taninah and Afanasjev (2023)). Note that this difference represents only 0.015% error in the description of total binding energy. In a similar way, total binding energies of the ¹²⁰Sn and ^102,110Zr nuclei obtained in the Skyrme DFT calculations with the SLy4 force in the basis of $N=20$ and $N=25$ harmonic oscillator (HO) shells differ by 110-150 keV Dobaczewski et al. (2004); Pei et al. (2008). Ref. Pei et al. (2008) also indicates that the HO basis with $N=25$ is needed in order to describe the binding energies of the nuclei in this mass region with an accuracy of the couple of tens keV.

Thus, the present paper aims at the investigation of several issues related to the calculations of binding energies with high precision in the CDFT framework. One of the goals is to understand the numerical errors in the calculations of binding energies. In particular, their global dependence on the truncation of the basis both in the fermionic and bosonic sectors of the CDFT will be investigated for the first time. We will also study how these numerical errors depend on the type of the symmetry (spherical versus axially deformed) of the nuclear shape and how the convergence of binding energies as a function of the size of the basis depends on the type of functional. The second goal is to understand to which extent the procedures for the definition of asymptotic binding energies employed earlier in non-relativistic functionals are applicable in the CDFT. The third goal is to further understand the role and the features of the anchor based optimization method suggested in Ref. Taninah and Afanasjev (2023) as an alternative to other methods of global fitting of the functionals.

The paper is organized as follows. A short review of the approaches to global fitting of the energy density functionals (EDF) is presented in Sec. II. Sec.  III provides a brief outline of theoretical formalism and the discussion of major classes of CEDFs under study. A discussion of the anchor based approach to the optimization of energy density functionals and its comparison with alternative methods for global fitting of the functionals is given in Sec. IV. The impact of the truncation of basis in bosonic and fermionic sectors of the CDFT on binding energies as well as asymptotic behavior of the calculated binding energies are analysed in Sec. V. A new method of defining asymptotic values of calculated binding energies is presented in Sec. VI. Finally, Sec. VII summarizes the results of our paper.

The analysis of literature reveals that at present there are three different approaches to global fitting of energy density functionals. These are

• •

  Fully Global Approach (further FGA). In the FGA the optimization of the functional is performed using computer codes which include either axial or triaxial deformation. It allows to use in fitting protocol the quantities which depend on deformation such as fission barriers which are not accessible in alternative fitting protocols.

• •

  Reduced Global Approach (further RGA) has been outlined in Sec. II of Ref. Tondeur et al. (2000). The specific feature of this approach is the fact that the deformation effects are taken into account in deformed and transitional nuclei by calculating the deformation energy $E_{def}$

  $\displaystyle E_{def}(Z,N)=E_{sph}(Z,N)-E_{eq}(Z,N)$

  (1)

  where $E_{eq}$ is the energy at the equilibrium deformation and $E_{sph}$ is the energy of spherical configuration. The $E_{def}$ values are calculated using axially deformed code and then all experimental binding energies are renormalized to their ”equivalent spherical configuration” values. Then, the rms errors of the current iteration of the functional are calculated using a spherical code by comparing the binding energies it gives with these renormalized experimental binding energies. Thus, each iteration of optimization consists of defining the $E_{def}$ in deformed code and fitting the functional in spherical code. However, to reach an optimal fit one should repeat iterations until full convergence of $E_{def}(Z,N)$ is reached.

• •

  Anchor Based Optimization Approach (further ABOA). In the ABOA (see Ref. Taninah and Afanasjev (2023)), the optimization of the parameters of EDFs is carried out for a selected set of spherical anchor nuclei, the physical observables of which are modified by the correction function

  $\displaystyle E_{corr}(Z,N)=\alpha_{i}(N-Z)+\beta_{i}(N+Z)+\gamma_{i}$

  (2)

  which takes into account the global performance of EDFs. Here $i$ is the counter of the iteration in the anchor-based optimization. The parameters $\alpha_{i}$ and $\beta_{i}$ correct the deficiencies of EDF in isospin and in mass number while $\gamma_{i}$ is responsible for a global shift in binding energies¹¹1Note that similar in spirit approach is used in Ref. Gonzalez-Boquera et al. (2018) in which the quadrupole zero-point energy is replaced by a constant binding energy shift [similar to the $\gamma_{i}$ term in Eq. (2)] which is fixed by minimizing the global rms deviation..

Each of these approaches has its own advantages and disadvantages which are discussed below.

The FGA potentially allows to avoid all uncertainties related to the selection of nuclei used in the fitting protocol and allows to use the physical quantities which depend on deformation (such as the height of fission barriers) in fitting protocol. However, the calculations in deformed codes are extremely numerically expensive and computational time does not depend drastically on the selection of the functional (see Table 1). For example, within the CDFT framework the numerical calculations in deformed RHB code are by more than two orders of magnitude more time consuming than those in spherical RHB one for the fermionic basis characterized by full $N_{F}=16$ fermionic shell (see columns 4 and 5 in Table 2). Further increase of fermionic basis leads to a drastic increase of the ratio of the calculational time in deformed and spherical RHB calculations: this ratio is approaching three orders of magnitude at $N_{F}=18$ (see Table 2).The simplification of the code by the transition to the relativistic mean field (RMF)+BCS framework does not change much this ratio. Similar situation holds also in non-relativistic Skyrme and Gogny DFT.

The FGA has been used in the number of studies but there are two possible ways of its application. In the first method, a large basis is used to reduce numerical errors in the calculations of binding energies but this restricts the number of the nuclei which can be included in the fitting protocol. However, the lack of strict mathematical procedure in the selection of nuclei leads to some subjectivity and possibly to some biased errors. This problem can be reduced by the inclusion of a larger set of nuclei into the fitting protocol. However, the computational cost raises drastically so that the increase of the number of nuclei leads to a substantial reduction of the basis size used in the calculations because of the limitations imposed on availability of computational resources. This, in turn requires the use of approximate relations to find the asymptotic value of the binding energy corresponding to infinite size of basis (see Sec. V for details). Thus, a very large number of the nuclei is used in the fitting protocol of the second method but this requires the use of moderate basis in numerical calculations and approximate relations to find the asymptotic value of the binding energy corresponding to infinite size of basis.

For example, the first method was used in Refs. Kortelainen et al. (2010, 2012, 2014) in Skyrme DFT for the development of the UNEDF class of the functionals. Note that fitting protocol of the UNEDF1 and UNEDF2 functionals has been supplemented by information on fission barriers and single-particle energies, respectively. In all optimizations, the large spherical basis with 20 harmonic oscillator shells ($N=20$) was used which led to quite substantial computational time in the calculations with deformed HFBTHO code. As a consequence, a limited number of nuclei (approximately 28 spherical and 46 deformed) were used in these optimizations.

The second method has been used in the fitting protocol of the BCPM (Barcelona-Catania-Paris-Madrid) functional (see Ref. Baldo et al. (2013)). It included 579 measured masses of even-even nuclei and the calculations have been carried out in a reduced basis which takes into account eleven harmonic oscillator shells for the $Z<50$ nuclei, 13 shells for the $50<Z<82$ nuclei, and 15 shells for $Z>82$ nuclei. Another example is the D1M and D1M* functionals which were fitted to 2149 and 620 nuclei, respectively, in axially deformed calculations corrected by quadrupole corrections to total binding energy and infinite-basis corrections $\Delta B_{\infty}$ [see Eq. (27) below] in Refs. Goriely et al. (2009); Gonzalez-Boquera et al. (2018). Note that only even-even nuclei were considered in the fit of the D1M* functional. The axially deformed calculations for these two functionals were performed in the basis with $N\leq 14$ major oscillator shells.

There is also a hybrid method which was used in the fitting of the DD-MEB1 and DD-MEB2 CEDFs in Refs. Peña-Arteaga et al. (2016): this fitting protocol contains 2353 experimental masses. The minimization of these functionals was performed in axially deformed RHB code with $N_{F}=16$ fermionic shells but with calculated binding energies corrected before and after minimization by the difference between the results obtained with $N_{F}=16$ and $N_{F}=20$ (see Ref. Goriely (2023a)).

The RGA has been used in the fitting protocol of the BSk1-32 series of the Skyrme functionals using very large bases (see, for example, Refs. Samyn et al. (2002); Goriely et al. (2013, 2016); Goriely (2023b)). For instance, it consisted of the $N=25$ HO shells in spherical calculations and $N=21$ HO shells in deformed calculations in Ref. Samyn et al. (2002). This allows to keep numerical error in calculated masses relatively low without involving the correction for infinite basis size. Note that fitting protocols of these functionals typically include all nuclei for which experimental data on binding energies is available.

The ABOA has been introduced in Ref. Taninah and Afanasjev (2023). It was shown that the use of this approach leads to a substantial improvement in the global description of binding energies for several classes of CEDFs (see also Sec. IV.1 below). This approach leads to the results which are comparable to those obtained in RGA (see Ref. Taninah and Afanasjev (2023) and the discussion in Sec. IV.3). The computational cost of defining a new functional within this approach is drastically lower as compared with FGA and RGA and this aspect is discussed in detail in Sec. IV.2 below. It turns out that the use of twelve anchor spherical nuclei distributed over the nuclear chart in the optimization is sufficient for this approach (see Ref. Taninah and Afanasjev (2023)). Note that existing optimizations of CEDFs within the ABOA have been obtained in the RHB calculations with $N_{F}=20$ fermionic and $N_{B}=20$ bosonic shells (see Ref. Taninah and Afanasjev (2023)). Such a truncation of the basis is used also in other recent CEDF optimizations²²2Unfortunately, the majority of the papers dedicated to the optimization of CEDFs do not explicitly provide information on the size of the basis used in the fitting protocols. (see, for example, Refs. Zhao et al. (2010); Taninah et al. (2020)). However, there are also the cases when smaller basis was used: for example, the DD-PC1 functional has been defined in deformed calculations with $N_{F}=16$ fermionic shells Nikšić et al. (2008).

It is necessary to mention that in these investigations the numerical errors in the binding energies associated with the truncation of basis or the definition of their asymptotic values have been either discussed only for a few selected nuclei or not mentioned at all. Thus, at present there is no published information on how these errors evolve across the nuclear landscape and how they depend on the type of the functional or theory employed. Thus, one of the goals of the present study is to close this gap in our knowledge in the CDFT framework.

The calculations were performed in the framework of relativistic Hartree-Bogoliubov (RHB) approach the technical details on the solution of which can be found in Refs. Vretenar et al. (2005); Agbemava et al. (2014); Nikšić et al. (2014). The calculations were carried out with two versions of the RHB code suitable for spherical and axially deformed nuclei, respectively. Both codes allow the calculations with a wide range of covariant energy density functionals. We employ a parallel version of the axial RHB computer code developed in Ref. Agbemava et al. (2014): it allows simultaneous calculations for a significant number of nuclei and deformation points in each nucleus. The unconstrained calculations are performed in axial RHB code using four initial deformations of the basis $\beta_{2}=-0.2,0.0,0.2$ and 0.4. As follows from the comparison with the calculations constrained on $\beta_{2}$ this procedure guarantees finding the lowest in energy minimum but it is substantially less numerically costly as compared with constrained calculations.

Both spherical and axial RHB computer codes are based on the expansion of the Dirac spinors and the meson fields in terms of harmonic oscillator wave functions with respective symmetry [so-called basis set expansion method] (see Ref. Gambhir et al. (1990); Nikšić et al. (2014)). We illustrate here this expansion for the case of axially deformed nuclei (see Refs. Gambhir et al. (1990); Nikšić et al. (2014) for more details and for the case of spherical symmetry).

The single-particle wave function $\psi_{i}(\vec{r},s,t)$ of the nucleon in the state $i$, which follows from the solution of the Dirac equation, is given by

$\displaystyle\psi_{i}(\vec{r},s,t)=\binom{f_{i}(\vec{r},s,t)}{ig_{i}(\vec{r},s,t)}$

(3)

where the large ($f_{i}(\vec{r},s,t)$) and small ($g_{i}(\vec{r},s,t)$) components of a Dirac spinor are expanded independently in terms of the HO eigenfunctions $\Phi_{\alpha}(\vec{r},s)$

	$\displaystyle f_{i}(\vec{r},s,t)=\sum_{\alpha}^{\alpha_{max}}f_{\alpha}^{(i)}\Phi_{\alpha}(\vec{r},s)\chi_{t_{i}}(t)$		(4)
	$\displaystyle g_{i}(\vec{r},s,t)=\sum_{\tilde{\alpha}}^{\tilde{\alpha}_{max}}f_{\tilde{\alpha}}^{(i)}\Phi_{\tilde{\alpha}}(\vec{r},s)\chi_{t_{i}}(t).$		(5)

The HO eigenfunctions have the form (see Ref. Nikšić et al. (2014))

$\displaystyle\Phi_{\alpha}(\vec{r},s)=\varphi_{n_{z}}(z_{,}b_{z})\varphi_{n_{r}}^{\Lambda}(r_{\perp},b_{\perp})\frac{e^{i\Lambda\phi}}{\sqrt{2\pi}}\chi(s)$

(6)

where $n_{z}$ and $n_{r}$ are the number of nodes in the $z$- and $r_{\perp}$- directions, respectively, and $\Lambda$ and $m_{s}$ are projections of the orbital angular momentum and spin on the intrinsic $z$-axis, respectively. $b_{z}$ and $b_{\perp}$ are oscillator lengths in respective directions.

The quantum numbers $\left|\alpha\right>=\left|n_{z}n_{r}\Lambda s\right>$ and $\left|\tilde{\alpha}\right>=\left|\tilde{n_{z}}\tilde{n_{r}}\tilde{\Lambda}\tilde{s}\right>$ in the sums of Eqs.  (4) and (5) run over all permitted combinations of $n_{z}$ and $n_{r}$, $\Lambda$ and $m_{s}$ (and their tilted counterparts) defined by the truncation of the basis. In absolute majority of the CDFT calculations, full fermionic HO shells are used in the basis set expansion and we follow this prescription in the present paper. Thus, $\alpha_{max}$ is selected in such a way that all basis states corresponding to full $N_{F}=2n_{r}+|\Lambda|+n_{z}$ fermionic shells are included into the basis of large components. The basis of small components includes $N_{F}+1$ full fermionic shells to avoid spurios contributions to the RHB equation (see Ref. Gambhir et al. (1990) for details). Note that large and small components of the Dirac spinor have opposite parities: this leads to an approximate doubling of the basis in the CDFT as compared with non-relativistic theories.

The solution of the Klein-Gordon equation is obtained by expanding mesonic fields in a complete set of basis states [HO eigenfunctions] along [$\varphi_{n_{z}}(z_{,}b_{z})$] and perpendicular [$\varphi_{n_{r}}^{0}(r_{\perp},b_{\perp})]$ the axis of symmetry of the nucleus

$\displaystyle\phi(z,r_{\perp})=\sum_{n_{z}n_{r}}^{N_{B}}\phi_{n_{z}n_{r}}\varphi_{n_{z}}(z_{,}b_{z})\varphi_{n_{r}}^{0}(r_{\perp},b_{\perp})$

(7)

This expansion is truncated at full $N_{B}$ bosonic shells.

The spherical RHB code is not parallelized but it allows the calculations to a very large value of $N_{F}=38$ (see Sec. V.2 below). Thus, the convergence of the results of the calculations in spherical nuclei as a function of $N_{F}$ has been tested using this code. In contrast, the axial RHB code works only up to $N_{F}=30$ due to memory allocation problems at larger values of $N_{F}$. Note that similar to Fig. 2 of Supplemental Material to Ref. Taninah and Afanasjev (2023) we verified for selected set of spherical and deformed nuclei that the codes developed in our group and those existing in the DIRHB package of the RHB codes (see Ref. Nikšić et al. (2014)) provide almost the same (within a few keVs) results for binding energies as a function of $N_{F}$ and $N_{B}$. Note that only DD-ME* and DD-PC* types of the functionals are included into the DIRHB codes of Ref. Nikšić et al. (2014).

In order to understand the dependence of the convergence properties on the functional we consider three major classes of CEDFs used at present. These are (i) those based on meson exchange with non-linear meson couplings (NLME) (see Ref. Boguta and Bodmer (1977)), (ii) those based on meson exchange with density dependent meson-nucleon couplings (DDME) (see Ref. Typel and Wolter (1999)), and finally (iii) those based on point coupling (PC) models containing various zero-range interactions in the Lagrangian (see Ref. Bürvenich et al. (2002)).

The Lagrangians of these three major classes of CEDFs can be written as: $\mathcal{L}=\mathcal{L}_{common}+\mathcal{L}_{model-specific}$ where the $\mathcal{L}_{common}$ consist of the Lagrangian of the free nucleons and the electromagnetic interaction. It is identical for all three classes of functionals and is written as

$\displaystyle\mathcal{L}_{common}=\mathcal{L}^{free}+\mathcal{L}^{em}$

(8)

with

$\displaystyle\mathcal{L}^{free}=\bar{\psi}(i\gamma_{{\mu}}\partial^{\mu}-m)\psi$

(9)

and

$\displaystyle\mathcal{L}^{em}=-\frac{1}{4}F^{\mu\nu}F_{\mu\nu}-e\frac{1-\tau_{3}}{2}\bar{\psi}\gamma^{\mu}\psi A_{\mu}.$

(10)

For each model there is a specific term in the Lagrangian: for the DDME models we have

	$\displaystyle\mathcal{L}_{DDME}$	$\displaystyle=$	$\displaystyle\frac{1}{2}(\partial{\sigma})^{2}-\frac{1}{2}m_{\sigma}^{2}\sigma^{2}-\frac{1}{4}\Omega_{\mu\nu}\Omega^{\mu\nu}+\frac{1}{2}m_{\omega}^{2}\omega^{2}$		(11)
		$\displaystyle-$	$\displaystyle\frac{1}{4}\vec{R}_{\mu\nu}\vec{R}^{\mu\nu}+\frac{1}{2}m_{\rho}^{2}\vec{\rho}^{\,2}-g_{\sigma}(\bar{\psi}\psi)\sigma$
		$\displaystyle-$	$\displaystyle g_{\omega}(\bar{\psi}\gamma_{\mu}\psi)\omega^{\mu}-g_{\rho}(\bar{\psi}\vec{\tau}\gamma_{\mu}\psi)\vec{\rho}^{\,\mu}$

with the density dependence of the coupling constants given by

	$\displaystyle g_{i}(\rho)$	$\displaystyle=$	$\displaystyle g_{i}(\rho_{0})f_{i}(x)~{}~{}~{}{\rm for}~{}i=\sigma,\omega$		(12)
	$\displaystyle g_{\rho}(\rho)$	$\displaystyle=$	$\displaystyle g_{\rho}(\rho_{0})\exp[-a_{\rho}(x-1)]$		(13)

where $\rho_{0}$ denotes the saturation density of symmetric nuclear matter and $x=\rho/\rho_{0}$. The functions $f_{i}(x)$ are given by the Typel-Wolter ansatz Typel and Wolter (1999)

$\displaystyle f_{i}(x)=a_{i}\frac{1+b_{i}(x+d_{i})}{1+c_{i}(x+d_{i})}.$

(14)

Because of the five conditions $f_{i}(1)=1$, $f^{\prime\prime}_{i}(1)=0$, and $f^{\prime\prime}_{\sigma}(1)=f^{\prime\prime}_{\omega}(1)$, only three of the eight parameters $a_{i}$, $b_{i}$, $c_{i}$, and $d_{i}$ are independent and we finally have the four parameters $b_{\sigma}$, $c_{\sigma}$, $c_{\omega}$, and $a_{\rho}$ characterizing the density dependence. In addition we have the four parameters of the Lagrangian $\mathcal{L}_{DDME}$ $m_{\sigma}$, $g_{\sigma}$, $g_{\omega}$, and $g_{\rho}$. As usual the masses of the $\omega$- and the $\rho$-meson are kept fixed at the values $m_{\omega}=783$ MeV and $m_{\rho}=763$ MeV Lalazissis et al. (2005); Roca-Maza et al. (2011). We therefore have $N_{par}=8$ parameters in the DDME class of the models. The DD-MEX Taninah et al. (2020), DD-ME2 Lalazissis et al. (2005) and DD-MEY Taninah and Afanasjev (2023) functionals have been used in the present study.

The NLME class of the functionals has the same model specific Lagrangian as the DDME class except that the coupling constants $g_{\sigma}$, $g_{\omega}$, and $g_{\rho}$ are constants and there are extra terms for a non-linear $\sigma$ meson coupling. These couplings are important for the description of surface properties of finite nuclei, especially the incompressibility Boguta and Bodmer (1977) and for nuclear deformations Gambhir et al. (1990).

$\displaystyle\mathcal{L}_{NLME}=\mathcal{L}_{DDME}-\frac{1}{3}g_{2}\sigma^{3}-\frac{1}{4}g_{3}\sigma^{4}$

(15)

For the NLME class we have $N_{par}=6$ parameters $m_{\sigma}$, $g_{\sigma}$, $g_{\omega}$, $g_{\rho}$, $g_{2}$, and $g_{3}$. In the present study, we use the NL5(E) CEDF from Ref. Agbemava et al. (2019a).

In general, the Lagrangian of the PC models contains three parts:
(i) the four-fermion point coupling terms:

$\displaystyle\begin{aligned} \mathcal{L}^{4f}=&-\frac{1}{2}\alpha_{S}(\bar{\psi}\psi)(\bar{\psi}\psi)-\frac{1}{2}\alpha_{V}(\bar{\psi}\gamma_{\mu}\psi)(\bar{\psi}\gamma^{\mu}\psi)&\\ &-\frac{1}{2}\alpha_{TS}(\bar{\psi}\vec{\tau}\psi)(\bar{\psi}\vec{\tau}\psi)-\frac{1}{2}\alpha_{TV}(\bar{\psi}\vec{\tau}\gamma_{\mu}\psi)(\bar{\psi}\vec{\tau}\gamma^{\mu}\psi),&\end{aligned}$

(16)

(ii) the gradient terms which are important to simulate the effects of finite range:

$\displaystyle\begin{aligned} \mathcal{L}^{der}=&-\frac{1}{2}\delta_{S}\partial_{\nu}{(\bar{\psi}\psi)}\partial^{\nu}{(\bar{\psi}\psi)}&\\ &-\frac{1}{2}\delta_{V}\partial_{\nu}{(\bar{\psi}\gamma_{\mu}\psi)}\partial^{\nu}{(\bar{\psi}\gamma^{\mu}\psi)}&\\ &-\frac{1}{2}\delta_{TS}\partial_{\nu}{(\bar{\psi}\vec{\tau}\psi)}\partial^{\nu}{(\bar{\psi}\vec{\tau}\psi)}&\\ &-\frac{1}{2}\delta_{TV}\partial_{\nu}{(\bar{\psi}\vec{\tau}\gamma_{\mu}\psi)}\partial^{\nu}{(\bar{\psi}\vec{\tau}\gamma^{\mu}\psi)},&\end{aligned}$

(17)

(iii) The higher order terms which are responsible for the effects of medium dependence

$\displaystyle\begin{aligned} \mathcal{L}^{hot}=&-\frac{1}{3}\beta_{S}(\bar{\psi}\psi)^{3}-\frac{1}{4}\gamma_{S}(\bar{\psi}\psi)^{4}&\\ &-\frac{1}{4}\gamma_{V}[(\bar{\psi}\gamma_{\mu}\psi)(\bar{\psi}\gamma^{\mu}\psi)]^{2}.&\end{aligned}$

(18)

However, the structure of two PC functionals, namely, PC-PK1 and DD-PC1, employed in the present paper depends on particular realization of this scheme. The PC-PK1 CEDF closely follows this scheme but it neglects the scalar-isovector channel, i.e. $\alpha_{TS}=\delta_{TS}=0$. This is because the information on masses and charge radii of finite nuclei does not allow one to distinguish the effects of the two isovector mesons $\delta$ and $\rho$ (see Ref. Roca-Maza et al. (2011)). Thus, for this PC CEDF we have $N_{par}=9$ parameters $\alpha_{S}$, $\alpha_{V}$, $\alpha_{TV}$, $\delta_{S}$, $\delta_{V}$, $\delta_{TV}$, $\beta_{S}$, $\gamma_{S}$, $\gamma_{V}$. The DD-PC1 functional Nikšić et al. (2008) introduces the density dependence of the $\alpha_{S}$, $\alpha_{V}$ and $\alpha_{TV}$ constants via

$\displaystyle\alpha_{i}(\rho)=a_{i}+(b_{i}+c_{i}x)e^{-d_{i}x}\qquad({\rm for}\,\,\,i=S,V,TV),$

and eliminates the terms proportional to $\alpha_{TS}$, $\delta_{TS}$ [this is similar to the PC-PK1 functional] and $\delta_{V}$, $\delta_{TV}$, $\beta_{S}$, $\gamma_{S}$ and $\gamma_{V}$ in Eqs. (16-18). The introduced density dependence of the $\alpha_{S}$, $\alpha_{V}$ and $\alpha_{TV}$ constants partially takes care of the impact of omitted terms. Note that in the isovector channel a pure exponential dependence is used, i.e. $a_{TV}=c_{TV}=0$. Thus, the DD-PC1 functionals contains the $N_{par}=10$ parameters, i.e. $a_{S}$, $b_{S}$, $c_{S}$, $d_{S}$, $a_{V}$, $b_{V}$, $c_{V}$, $d_{V}$, $b_{TV}$ and $d_{TV}$.

Note that in the deformed RHB calculations we use $ngh=30$ and $ngl=30$ integration points in the Gauss-Hermite and Gauss-Laguerre quadratures in the $z>0$ and $r_{\perp}>0$ directions, respectively, in order to properly take into account the details of the single-particle wave functions. The same $ngh=30$ is used in the Gauss-Hermite quadratures for $r>0$ in spherical RHB calculations.

The experimental binding energies of twelve anchor spherical nuclei (¹⁶O, ^40,48Ca, ⁷²Ni, ⁹⁰Zr, ^116,124,132Sn, ^204,208,214Pb and ²¹⁰Po) and charge radii of nine of them (namely, ¹⁶O, ^40,48Ca, ⁹⁰Zr, ^116,124Sn and ^204,208,214Pb) are used in the anchor based optimization approach (see Ref. Taninah and Afanasjev (2023) for more details). The correction function of Eq. (2) is defined in ABOA using 855 even-even nuclei for which experimental binding energies are available in Ref. Wang et al. (2017).

The anchor-based optimization of energy density functionals is defined in Ref. Taninah and Afanasjev (2023) and the basic results for this approach have been discussed in this reference and in Supplemental Material to it. It is based on the introduction of the correction function $E_{corr}(Z,N)$ given in Eq. (2) by which the binding energies of spherical anchor nuclei are redefined for global performance of the EDFs (see Ref. Taninah and Afanasjev (2023)) for full details). The behavior of these parameters during the iterative procedure has been illustrated for the DD-MEX1, DD-MEY and NL5(Y) functionals in Tables I, II and III of the Supplemental Material to Ref. Taninah and Afanasjev (2023). In addition, the evolution of the rms deviations $\Delta E_{rms}$ between calculated and experimental binding energies $E$ for these functionals during the iterative procedure was also illustrated in Fig. 1 of Supplemental Material to Ref. Taninah and Afanasjev (2023).

In this section we extend the discussion of the ABOA to other functionals, introduce new improved strategy for the fitting of the functionals within ABOA, provide an estimate of the computation time for the ABOA and RGA approaches and compare the results obtained for the DD-MEY type of the functionals in the RGA and ABOA.

The values of $N_{F}$ and $N_{B}$ employed in the RMF+BCS or RHB calculations depend on available computational power [which increases with time] and on the approach of theory group to handling the numerical errors in physical observables. In absolute majority of such calculations, $N_{B}$ is fixed at 20 since the calculations in bosonic sector of the models are relatively cheap. In contrast, the size of the fermionic basis shows substantial variations and dependence on time of publication, the sophistication of theoretical approach and theory group since the cost of numerical calculations in the RMF+BCS and RHB approaches is mostly defined by this sector.

Let us briefly review truncation schemes of fermionic basis of existing global mass calculations in the CDFT framework. Note that we specify $N_{B}$ below only if it differs from 20. The $N_{F}=12$ value has been used in global RMF+BCS mass calculations with NL3 CEDF in Ref. Lalazissis et al. (1999) with a larger number of fermionic shells used in very heavy nuclei. The RMF+BCS calculations of Ref. Geng et al. (2005) have been carried with $N_{F}=N_{B}=14$ and TMA CEDF. Global RMF+BCS calculations with CEDFs NL3, TM1, FSUGold and BSR4 were performed in Ref. Reinhard and Agrawal (2011) with $N_{F}=18$. $N_{F}=20$ was used in global RHB calculations with the DD-ME2, DD-PC1, NL3* and DD-ME$\delta$ functionals in Ref.  Agbemava et al. (2014). Ref. Zhang et al. (2014) presents the results of the global RMF+BCS calculations with phenomenological treatment of the dynamical correlation energy and CEDF PC-PK1: it uses $N_{F}=14$ for $Z<60$ nuclei and $N_{F}=18$ for heavier ones. The transition to microscopic beyond mean field approaches leads to a significant increase of computational time. As a consequence, the fermionic basis is reduced in such studies. For example, the first global investigation within five-dimensional collective Hamiltonian (5DCH) based on triaxial RMF+BCS calculations with PC-PK1 CEDF employs $N_{F}=12$, 14 and 16 for the nuclei with $Z<20$, $20\leq Z\leq 82$, and $Z\geq 82$, respectively (see Ref. Lu et al. (2015)). The same fermionic basis is used in the 5DCH studies based on triaxial RHB calculations with PC-PK1 CEDF presented in Ref. Yang et al. (2021).

While the majority of these schemes of the truncation of the fermionic basis provide a sufficient accuracy for absolute majority of physical observables (such as radii, deformations, etc) it is still an open question on how accurate they are for the calculations of the binding energies. This is because the assessment of the accuracy of the truncation of the basis in these publications has been either not carried or performed by comparing the solutions obtained with $N_{F}$ and $N_{F}+2$ (or $N_{F}+6$ as in Ref. Agbemava et al. (2014)) for very restricted set of nuclei. As a consequence, a number of the issues has not even been considered in earlier publications. For example, how the truncation errors in binding energies depend on the type of the functional? Or what is the evolution of truncation errors with proton and neutron numbers? Are those errors gradually increase with particle numbers or there are some local fluctuations caused by the shell effects?

To address these questions one should consider an asymptotic limit for the binding energies $B(Z,N)$⁹⁹9$B(Z,N)$ is the (negative) binding energy of the nucleus with $Z$ protons and $N$ neutrons. which corresponds to $N_{F}\rightarrow\infty$ and $N_{B}\rightarrow\infty$. The asymptotic (negative) binding energy corresponding to these limits is defined as $B_{\infty}(Z,N)=B(N_{F}\rightarrow\infty,N_{B}\rightarrow\infty)$ and it serves as a benchmark with respect of which the truncation errors in binding energies in the basis with ($N_{F},N_{B}$) have to be defined. The consideration of this limit and related truncation errors are presented below for bosonic and fermionic sectors separately.