Model nuclear energy density functionals derived from ab initio calculations

We separate the nuclear many-body Hamiltonian Ĥ = Ĥ₀ + Ĥ₁ into a noninteracting part $\hat{{H}_{0}}=\hat{T}+\hat{U}$, with the auxiliary one-body operator Û, and the interacting part defined as $\hat{{H}_{1}}=-\hat{U}+{\hat{V}}^{2\text{B}}+{\hat{V}}^{3\text{B}}$. We express the nuclear Hamiltonian in second quantisation as

$$\begin{align}\hfill \hat{H}=& \sum _{\alpha }{{\epsilon}}_{\alpha }^{0}{a}_{\alpha }^{{\dagger}}{a}_{\alpha }-\sum _{\alpha \beta }\langle \alpha \vert \hat{U}\vert \beta \rangle {a}_{\alpha }^{{\dagger}}{a}_{\beta }+\frac{1}{4}\sum _{\alpha \beta \gamma \delta }\langle \alpha \beta \vert {\hat{V}}^{2\text{B}}\vert \gamma \delta \rangle {a}_{\alpha }^{{\dagger}}{a}_{\beta }^{{\dagger}}{a}_{\delta }{a}_{\gamma }\hfill \\ \hfill & +\frac{1}{36}\sum _{\alpha \beta \mu \gamma \delta \nu }\langle \alpha \beta \mu \vert {\hat{V}}^{3\text{B}}\vert \gamma \delta \nu \rangle {a}_{\alpha }^{{\dagger}}{a}_{\beta }^{{\dagger}}{a}_{\mu }^{{\dagger}}{a}_{\nu }{a}_{\delta }{a}_{\gamma },\hfill \end{align} \tag{ 1 }$$

where ${{\epsilon}}_{\alpha }^{0}$ are the single-particle energies of Ĥ₀. By solving the corresponding Schrödinger equation,

$$\begin{equation}\hat{H}\vert {{\Psi}}_{0}^{A}\rangle ={E}_{0}^{A}\vert {{\Psi}}_{0}^{A}\rangle ,\end{equation} \tag{ 2 }$$

one obtains the ground-state energy ${E}_{0}^{A}$ and wave function $\vert {{\Psi}}_{0}^{A}\rangle$ of the nuclear system. The exact solution of equation (2) is very complicated and it is typically limited to systems with very few number of nucleons [31]. Rather than calculating the full many-body wave function, it is possible to expand the solution of the Schrödinger equation in terms of the propagation of single-particle excitations and the correlated density matrix of the system by using SCGF method [29]. These excitations represent basic collective degrees of freedom of the nucleus and they are described by the one-body Green's function, or propagator G_αβ. Using the Källén–Lehmann representation [32, 33], G_αβ takes the form

$$\begin{equation}{G}_{\alpha \beta }\left(E\right)=\sum _{n}\frac{\langle {{\Psi}}_{0}^{A}\vert {a}_{\alpha }\vert {{\Psi}}_{n}^{A+1}\rangle \langle {{\Psi}}_{n}^{A+1}\vert {a}_{\beta }^{{\dagger}}\vert {{\Psi}}_{0}^{A}\rangle }{E-\left({E}_{n}^{A+1}-{E}_{0}^{A}\right)+\mathrm{i}\eta }+\sum _{k}\frac{\langle {{\Psi}}_{0}^{A}\vert {a}_{\beta }^{{\dagger}}\vert {{\Psi}}_{k}^{A-1}\rangle \langle {{\Psi}}_{k}^{A-1}\vert {a}_{\alpha }\vert {{\Psi}}_{0}^{A}\rangle }{E-\left({E}_{0}^{A}-{E}_{k}^{A-1}\right)-\mathrm{i}\eta }.\end{equation} \tag{ 3 }$$

Here the complete set of eigenstates $\vert {{\Psi}}_{n}^{A+1}\rangle$, $\vert {{\Psi}}_{k}^{A-1}\rangle$ with eigenvalues ${E}_{n}^{A+1}$, ${E}_{k}^{A-1}$ introduce the intermediate (A ± 1)-body systems. Then, ${E}_{n}^{A+1}-{E}_{0}^{A}$ and ${E}_{0}^{A}-{E}_{k}^{A-1}$ are, respectively, the excitation energies of the propagating quasi-particle (index n) and quasi-hole (index k) states.

The propagator given in equation (3) satisfies the Dyson equation

$$\begin{equation}{G}_{\alpha \beta }\left(E\right)={G}_{\alpha \beta }^{\left(0\right)}\left(E\right)+\sum _{\gamma \delta }{G}_{\alpha \gamma }^{\left(0\right)}\left(E\right)\;{{\Sigma}}_{\gamma \delta }^{\star }\left(E\right)\;{G}_{\delta \beta }\left(E\right),\end{equation} \tag{ 4 }$$

where ${G}_{\alpha \beta }^{\left(0\right)}\left(E\right)$ is the propagator for the noninteracting system Ĥ₀, and ${{\Sigma}}_{\gamma \delta }^{\star }\left(E\right)$ is the irreducible self-energy. The latter represents the nonlocal and energy-dependent potential to which each nucleon is subjected when interacting within the nuclear medium. The nonlinearity of equation (4) in terms of G_αβ(E) requires an iterative procedure to reach convergence. The full self-consistency is required to satisfy fundamental symmetries and conservation laws [34].

A crucial ingredient of equation (4) is the self-energy. This is composed of three parts as

$$\begin{equation}{{\Sigma}}_{\gamma \delta }^{\star }\left(E\right)=-\langle \gamma \vert U\vert \delta \rangle +{{\Sigma}}_{\gamma \delta }^{\left(\infty \right)}+{\tilde {{\Sigma}}}_{\gamma \delta }\left(E\right),\end{equation} \tag{ 5 }$$

which, respectively, are the auxiliary potential, static mean-field, and energy-dependent component. To perform calculations of the energy-dependent component, in this work we employ the algebraic diagrammatic construction method [35, 36] up to third order, denoted by ADC(3). We refer to reference [37] for a more pedagogical introduction of the method. The accuracy of ADC(3) can be estimated to be of the fourth order in Ĥ₁, giving in practice an error of 1% of the total binding energy [37].

Another important aspect of the ab initio calculations is the presence of explicit three-body interaction ${\hat{V}}^{3\text{B}}$. By means of the modified Migdal–Galitski–Koltun sum rule [38], we can write the ground-state energy of the system as,

$$\begin{equation}{E}_{0}^{A}=\sum _{\alpha \beta }\frac{1}{2\pi }{\int }_{-\infty }^{{{\epsilon}}_{0}^{-}}\mathrm{d}E\left[\langle \alpha \vert \hat{T}\vert \beta \rangle +E{\delta }_{\alpha \beta }\right]\;\mathrm{Im}\left\{{G}_{\beta \alpha }\left(E\right)\right\}-\frac{1}{2}\langle {{\Psi}}_{0}^{A}\vert {\hat{V}}^{3\text{B}}\vert {{\Psi}}_{0}^{A}\rangle ,\end{equation} \tag{ 6 }$$

where ${{\epsilon}}_{0}^{-}$ is the highest quasi-hole energy, namely ${{\epsilon}}_{0}^{-}={\mathrm{max}}_{k}\left({E}_{0}^{A}-{E}_{k}^{A-1}\right)$. The determination of the expectation value of the three-body interaction would require calculation of many-body propagators. Instead, assuming that the magnitude of the three-body term is smaller than that of the two-body term, we take the lowest order approximation for the Hamiltonian, leading to

$$\begin{equation}\langle {{\Psi}}_{0}^{A}\vert {\hat{V}}^{3\text{B}}\vert {{\Psi}}_{0}^{A}\rangle \approx \frac{1}{6}\sum _{\alpha \beta \mu \gamma \delta \nu }\langle \alpha \beta \mu \vert {\hat{V}}^{3\text{B}}\vert \gamma \delta \nu \rangle \;{\rho }_{\gamma \alpha }{\rho }_{\delta \beta }{\rho }_{\nu \mu },\end{equation} \tag{ 7 }$$

where ρ_αβ is the one-body density matrix.

Having defined the main theoretical aspects of the ab initio method employed in this article, we now define the model energy density functional as

$$\begin{equation}\tilde {E}\left[\rho \right]={T}^{1\text{B}}\left[\rho \right]+{V}^{\text{Coul}}\left[\rho \right]+\sum _{j}{C}_{j}\;{V}_{j}^{\text{gen}}\left[\rho \right].\end{equation} \tag{ 8 }$$

The different terms correspond to average values, evaluated with respect to a Hartree–Fock (HF) state, of the kinetic energy ${T}^{1\text{B}}\left[\rho \right]\equiv \langle {\Phi}\vert {\hat{T}}^{1\text{B}}\vert {{\Phi}\rangle }_{\text{HF}}$, Coulomb potential ${V}^{\text{Coul}}\left[\rho \right]\equiv \langle {\Phi}\vert {\hat{V}}^{\text{Coul}}\vert {{\Phi}\rangle }_{\text{HF}}$, and interaction components ${V}_{j}^{\text{gen}}\left[\rho \right]\equiv \langle {\Phi}\vert {\hat{V}}_{j}^{\text{gen}}\vert {{\Phi}\rangle }_{\text{HF}}$. The operators ${\hat{V}}_{j}^{\text{gen}}$ represent our choice of the generators to build the functional, whereas C_j are the coupling constants we need to adjust on (meta)-data.

In the present article, we define generators ${\hat{V}}_{j}^{\text{gen}}$, based on ten individual terms ${\hat{T}}_{i}$ and ${\hat{T}}_{i}^{\sigma }$ for i = 0, 1, 2, and ${\hat{T}}_{e}$, ${\hat{T}}_{o}$, ${\hat{T}}_{W0}$, and ${\hat{T}}_{3}$ of the Skyrme functional generator [12], that is,

$$\begin{align}\hfill {\hat{V}}^{\text{Skyrme}}& ={t}_{0}\left(1+{x}_{0}{\hat{P}}^{\sigma }\right)\delta \left({\boldsymbol{r}}_{\boldsymbol{1}}-{\boldsymbol{r}}_{\boldsymbol{2}}\right)+\frac{1}{2}{t}_{1}\left(1+{x}_{1}{\hat{P}}^{\sigma }\right)\left[{\hat{\boldsymbol{k}}}^{\prime 2}\delta \left({\boldsymbol{r}}_{\boldsymbol{1}}-{\boldsymbol{r}}_{\boldsymbol{2}}\right)+\delta \left({\boldsymbol{r}}_{\boldsymbol{1}}-{\boldsymbol{r}}_{\boldsymbol{2}}\right){\hat{\boldsymbol{k}}}^{2}\right]\hfill \\ \hfill & \quad +\;{t}_{2}\left(1+{x}_{2}{\hat{P}}^{\sigma }\right){\hat{\boldsymbol{k}}}^{\prime }\cdot \delta \left({\boldsymbol{r}}_{\boldsymbol{1}}-{\boldsymbol{r}}_{\boldsymbol{2}}\right)\hat{\boldsymbol{k}}+i{W}_{0}\left(\hat{{\boldsymbol{\sigma }}_{\boldsymbol{1}}}+\hat{{\boldsymbol{\sigma }}_{\boldsymbol{2}}}\right)\cdot \left[{\hat{\boldsymbol{k}}}^{\prime }{\times}\delta \left({\boldsymbol{r}}_{\boldsymbol{1}}-{\boldsymbol{r}}_{\boldsymbol{2}}\right)\hat{\boldsymbol{k}}\right]\hfill \\ \hfill & \quad +\;\frac{{t}_{e}}{2}\left\{\left[3\left({\boldsymbol{\sigma }}_{\boldsymbol{1}}\cdot {\boldsymbol{k}}^{\prime }\right)\left({\boldsymbol{\sigma }}_{\boldsymbol{2}}\cdot {\boldsymbol{k}}^{\prime }\right)-\left({\boldsymbol{\sigma }}_{\boldsymbol{1}}\cdot {\boldsymbol{\sigma }}_{\boldsymbol{2}}\right){{\boldsymbol{k}}^{\prime }}^{2}\right]\delta \left({\boldsymbol{r}}_{\boldsymbol{1}}-{\boldsymbol{r}}_{\boldsymbol{2}}\right)\right.\hfill \\ \hfill & \left.\quad +\;\delta \left({\boldsymbol{r}}_{\boldsymbol{1}}-{\boldsymbol{r}}_{\boldsymbol{2}}\right)\left[3\left({\boldsymbol{\sigma }}_{\boldsymbol{1}}\cdot \boldsymbol{k}\right)\left({\boldsymbol{\sigma }}_{\boldsymbol{2}}\cdot \boldsymbol{k}\right)-\left({\boldsymbol{\sigma }}_{\boldsymbol{1}}\cdot {\boldsymbol{\sigma }}_{\boldsymbol{2}}\right){\boldsymbol{k}}^{2}\right]\right\}\hfill \\ \hfill & \quad +\;{t}_{o}\left\{3\left({\boldsymbol{\sigma }}_{\boldsymbol{1}}\cdot {\boldsymbol{k}}^{\prime }\right)\delta \left({\boldsymbol{r}}_{\boldsymbol{1}}-{\boldsymbol{r}}_{\boldsymbol{2}}\right)\left({\boldsymbol{\sigma }}_{\boldsymbol{2}}\cdot \boldsymbol{k}\right)-\left({\boldsymbol{\sigma }}_{\boldsymbol{1}}\cdot {\boldsymbol{\sigma }}_{\boldsymbol{2}}\right)\left[{\boldsymbol{k}}^{\prime }\cdot \delta \left({\boldsymbol{r}}_{\boldsymbol{1}}-{\boldsymbol{r}}_{\boldsymbol{2}}\right)\boldsymbol{k}\right]\right\}\hfill \\ \hfill & \quad +\;{t}_{3}\delta \left({\boldsymbol{r}}_{\boldsymbol{1}}-{\boldsymbol{r}}_{\boldsymbol{2}}\right)\delta \left({\boldsymbol{r}}_{\boldsymbol{2}}-{\boldsymbol{r}}_{\boldsymbol{3}}\right){\mathcal{A}}_{\mathrm{123}}\left({\hat{P}}^{\sigma },{\hat{P}}^{\tau }\right)\hfill \\ \hfill & ={t}_{0}{\hat{T}}_{0}+{t}_{0}{x}_{0}{\hat{T}}_{0}^{\sigma }+{t}_{1}{\hat{T}}_{1}+{t}_{1}{x}_{1}{\hat{T}}_{1}^{\sigma }+{t}_{2}{\hat{T}}_{2}+{t}_{2}{x}_{2}{\hat{T}}_{2}^{\sigma }\hfill \\ \hfill & \quad +\;{W}_{0}{\hat{T}}_{W0}+{t}_{e}{\hat{T}}_{e}+{t}_{o}{\hat{T}}_{o}+{t}_{3}{\hat{T}}_{3},\hfill \end{align} \tag{ 9 }$$

where ${\mathcal{A}}_{\mathrm{123}}$ is the antisymmetric operator for combinations of ${\hat{P}}^{\sigma }=\frac{1}{2}\left(1+{\boldsymbol{\sigma }}_{\boldsymbol{1}}\cdot {\boldsymbol{\sigma }}_{\boldsymbol{2}}\right)$ and ${\hat{P}}^{\tau }=\frac{1}{2}\left(1+{\boldsymbol{\tau }}_{\boldsymbol{1}}\cdot {\boldsymbol{\tau }}_{\boldsymbol{2}}\right)$. We refer to reference [39] for more details.

To determine the coupling constants of the interaction using ab initio methods, it is convenient to transform the generators given in equation (9) to linear combinations that give specific isoscalar/isovector terms of the functional [28]⁹. This is done by means of the following matrix relation,

$$\begin{equation}\left[\begin{bmatrix}{c}\hfill {\hat{V}}_{0}^{\rho }\hfill \\ \hfill {\hat{V}}_{1}^{\rho }\hfill \\ \hfill {\hat{V}}_{0}^{{\Delta}\rho }\hfill \\ \hfill {\hat{V}}_{1}^{{\Delta}\rho }\hfill \\ \hfill {\hat{V}}_{0}^{\tau }\hfill \\ \hfill {\hat{V}}_{1}^{\tau }\hfill \\ \hfill {\hat{V}}_{0}^{J1}\hfill \\ \hfill {\hat{V}}_{1}^{J1}\hfill \\ \hfill {\hat{V}}_{W0}\hfill \\ \hfill {\hat{V}}_{t3}\hfill \end{bmatrix}\right]=\left[\begin{bmatrix}{cccccccccc}\hfill \frac{8}{3}\hfill & \hfill -\frac{4}{3}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -4\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -\frac{16}{3}\hfill & \hfill \frac{8}{3}\hfill & \hfill \frac{16}{3}\hfill & \hfill -\frac{8}{3}\hfill & \hfill -\frac{16}{15}\hfill & \hfill \frac{16}{15}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 8\hfill & \hfill -\frac{32}{3}\hfill & \hfill \frac{40}{3}\hfill & \hfill \frac{16}{5}\hfill & \hfill \frac{16}{15}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{4}{3}\hfill & \hfill -\frac{2}{3}\hfill & \hfill 4\hfill & \hfill -2\hfill & \hfill -\frac{8}{15}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -2\hfill & \hfill -8\hfill & \hfill 10\hfill & \hfill \frac{8}{5}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{8}{5}\hfill & \hfill \frac{8}{5}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -\frac{24}{5}\hfill & \hfill \frac{8}{5}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \end{bmatrix}\right]\left[\begin{bmatrix}{c}\hfill {\hat{T}}_{0}\hfill \\ \hfill {\hat{T}}_{0}^{\sigma }\hfill \\ \hfill {\hat{T}}_{1}\hfill \\ \hfill {\hat{T}}_{1}^{\sigma }\hfill \\ \hfill {\hat{T}}_{2}\hfill \\ \hfill {\hat{T}}_{2}^{\sigma }\hfill \\ \hfill {\hat{T}}_{e}\hfill \\ \hfill {\hat{T}}_{o}\hfill \\ \hfill {\hat{T}}_{W0}\hfill \\ \hfill {\hat{T}}_{3}\hfill \end{bmatrix}\right].\end{equation} \tag{ 10 }$$

We explicitly included only the vector part ${\hat{V}}_{T}^{J1}$ of the tensor term [40], which in the case of spherical symmetry considered here gives the only non-zero contribution. Note that the terms associated with the ${\hat{V}}_{W0}$ and ${\hat{V}}_{t3}$ generators do not allow for the separation of isoscalar/isovector terms, and therefore we keep them identical to the corresponding terms of the Skyrme generator in equation (9). The generators listed on the left-hand side of equation (10) are then used as generators, ${\hat{V}}_{j}^{\text{gen}}$, j = 1–10, that define the functional in equation (8).

In the electronic density functional theory [41], one uses the Levy–Lieb constrained variation [42–45] to obtain the ground state energy of the system. This procedure consists in a two-step minimisation,

$$\begin{equation}{E}_{\text{g.s.}}={\mathrm{min}}_{\rho }\left\{{\mathrm{min}}_{\vert {\Psi}\rangle \to \rho }\left[\langle {\Psi}\vert \hat{T}+\hat{V}\vert {\Psi}\rangle \right]\right\}\equiv {\mathrm{min}}_{\rho }E\left[\rho \right],\end{equation} \tag{ 11 }$$

where $\hat{V}$ stands for the Coulomb potential and symbol min_|Ψ⟩→ρ denotes an inner (first-step) minimisation over all correlated many-body states |Ψ⟩ that have a common fixed one-body density profile ρ(r). The outer (second-step) minimisation is then performed over all possible profiles ρ(r), and, in such a way, the global minimum of energy, and thus the exact ground-state energy E_g.s. is obtained.

The inner minimisation can be conveniently performed by an unconstrained minimisation of the Routhian $\hat{R}$ at fixed one-body external potential U(r) that plays a role of the Lagrange multiplier,

$$\begin{equation}R\left[U\right]={\mathrm{min}}_{\vert {\Psi}\rangle }\langle {\Psi}\vert \hat{R}\vert {\Psi}\rangle ={\mathrm{min}}_{\vert {\Psi}\rangle }\langle {\Psi}\vert \left[\hat{T}+\hat{V}+\int \mathrm{d}\boldsymbol{r}\;U\left(\boldsymbol{r}\right)\hat{\rho }\left(\boldsymbol{r}\right)\right]\vert {\Psi}\rangle .\end{equation} \tag{ 12 }$$

This gives the energy E[U] = R[U] − ∫drU(r)ρ(r) and density ρ[U] as functionals of the potential U(r). Assuming that the inverse functional U[ρ] can be found, we obtain the exact energy density functional,

$$\begin{equation}E\left[\rho \right]\equiv E\left[U\left[\rho \right]\right],\end{equation} \tag{ 13 }$$

which after the outer (second-step) minimisation gives again the exact ground-state energy E_g.s..

In the case of a nuclear system, we consider the analogous first-step variation consisting in the minimisation of the Routhian ${\hat{R}}^{\text{ab}}$ as

$$\begin{equation}{\delta }_{{\Psi}}\langle {\Psi}\vert {\hat{R}}^{\text{ab}}\vert {\Psi}\rangle ={\delta }_{{\Psi}}\langle {\Psi}\vert \left[{\hat{H}}^{\text{ab}}+\int \mathrm{d}\boldsymbol{r}\;U\left(\boldsymbol{r}\right)\hat{\rho }\left(\boldsymbol{r}\right)\right]\vert {\Psi}\rangle =0,\end{equation} \tag{ 14 }$$

where Ĥ^ab is the Hamiltonian of the system, equation (1). We use the superscript ab to indicate that it is the Hamiltonian of the ab initio theory, distinguishing it from the Hamiltonian used to build the model functional. The minimisation gives us many-body states as functionals of the external potential, |Ψ(U)⟩.

The integrand on the right-hand side of equation (14) introduces a perturbation of the ground-state |Ψ_g.s.⟩. The response of the system to the perturbation causes a change in the density ρ. If we were able to probe the system with all possible potentials U(r), we would have obtained the functional E^ab[U] ≡ ⟨Ψ(U)|Ĥ^ab|Ψ(U)⟩, and then, as above, the exact energy density functional E^ab[ρ]. In the second step, the functional E^ab[ρ] is minimised with respect to ρ, which gives the exact ground-state energy ${E}_{\text{g.s.}}^{\text{ab}}$ and density ρ_g.s.(r).

Being unable to perturb the system with an infinite number of external potentials, let us introduce a discrete finite set of pre-defined external potentials u_i(r) and their corresponding strengths λ_i, whereby the first-step minimisation (14) becomes,

$$\begin{equation}{\delta }_{{\Psi}}\langle {\Psi}\vert {\hat{R}}^{\text{ab}}\vert {\Psi}\rangle ={\delta }_{{\Psi}}\langle {\Psi}\vert \left[{\hat{H}}^{\text{ab}}+\sum _{i}{\lambda }_{i}\int \mathrm{d}\boldsymbol{r}\;{u}_{i}\left(\boldsymbol{r}\right)\hat{\rho }\left(\boldsymbol{r}\right)\right]\vert {\Psi}\rangle =0.\end{equation} \tag{ 15 }$$

With this restriction, the ab initio energy and density, E^ab(λ_i) and ρ^ab(λ_i), become functions (not functionals) of the finite set of strengths λ_i. Obviously, we now do not know what is the full energy density functional E^ab[ρ] in the infinite-dimensional space of all possible one-body density profiles, however, we know it exactly on a finite-dimensional manifold of densities parametrised by strengths λ_i. Recall that this manifold still contains the point λ_i = 0 that corresponds to the exact ground state. The second-step variation would now correspond to the minimisation of function E^ab(λ_i) in finite dimensions.

As proposed in reference [28], we now conjecture that a meaningful manifold of ab initio densities ρ^ab(λ_i) can be obtained not by pre-defining external potentials u_i(r), but by perturbing the system with generators that are going to be used for modelling the functional, section 2.2, that is,

$$\begin{equation}{\delta }_{{\Psi}}\langle {\Psi}\vert {\hat{R}}^{\text{ab}}\vert {\Psi}\rangle ={\delta }_{{\Psi}}\langle {\Psi}\vert \left[{\hat{H}}^{\text{ab}}+\sum _{i}{\lambda }_{i}{\hat{V}}_{i}^{\text{gen}}\right]\vert {\Psi}\rangle =0.\end{equation} \tag{ 16 }$$

Since the unrestricted minimisation of the Routhian is equivalent to finding its exact ground state |Ψ^ab(λ_i)⟩ and eigenvalue R^ab(λ_i), we have

$$\begin{equation}\left[{\hat{H}}^{\text{ab}}+\sum _{i}{\lambda }_{i}{\hat{V}}_{i}^{\text{gen}}\right]\vert {{\Psi}}^{\text{ab}}\left({\lambda }_{i}\right)\rangle ={R}^{\text{ab}}\left({\lambda }_{i}\right)\vert {{\Psi}}^{\text{ab}}\left({\lambda }_{i}\right)\rangle .\end{equation} \tag{ 17 }$$

Formally, by replacing minimisation (15) with minimisation (16), we do not make any new assumptions. Indeed, the previously made assumption about the reversibility of the functional ρ[U] suffices. It stipulates that for any density ρ(r) generated by the latter minimisation there exists a potential U(r) that would have generated it by the former minimisation. However, when using minimisation (16) we do not have to pre-define any potential or, for that matter, we do not have to know it at all.

Since our goal is to use the ab initio energies with perturbation λ_i as metadata to determine the coupling constants of our model functional, we impose that the ab initio energy can be expressed in the form of functional (8), that is,

$$\begin{equation}{E}^{\text{ab}}\left({\lambda }_{i}\right)={T}^{1\text{B}}\left[{\rho }^{\text{ab}}\left({\lambda }_{i}\right)\right]+{V}^{\text{Coul}}\left[{\rho }^{\text{ab}}\left({\lambda }_{i}\right)\right]+\sum _{j}{C}_{j}\;{V}_{j}^{\text{gen}}\left[{\rho }^{\text{ab}}\left({\lambda }_{i}\right)\right].\end{equation} \tag{ 18 }$$

Equation (18) constitutes the central point of our approach. At first glance, it looks like an impossible miracle: it expresses the exact ab initio energy as a functional of the exact ab initio one-body density profile ρ(r) only. However, this is, in fact, exactly the main statement of the rigorous DFT [41, 46, 47]. Without repeating here the related arguments, at this point it is nevertheless worth summarizing the three assumptions we did make when formulating equation (18):

• (a)  
  We use a specific form of the density functional restricted to a few terms given by Skyrme generators (9). Obviously, the exact functional of equation (13) does not have to have such or a similar form. This is why in our approach we aim to derive a model functional and not the exact one. Nevertheless, a successful phenomenology of modelling nuclear phenomena by Skyrme-type functionals supports the idea of the model formulated in equation (18).
• (b)  
  Skyrme functional depends not only on the local density profiles ρ(r) for neutrons and protons, but also on a few second-order quasilocal densities [40, 48]. Although it can be rigorously proved that any arbitrary local density profile can be modelled by a density profile of a Kohn–Sham Slater determinant [49], a simultaneous modelling of local and quasilocal densities is only an approximation.
• (c)  
  When matching left- and right-hand sides of equation (18), we use a specific finite set of Lagrange multipliers λ_i. Had the form of functional on the right-hand side been exact, we could have used any set of λ_i. However, for the specific model functional we use, the results may depend on the specific values of the Lagrange multipliers. Our choice of using values around λ_i = 0 constitutes, therefore, another approximation.

In addition, considering that the Kohn–Sham kinetic energy represents a good approximation to the one-body kinetic energy, we assume that

$$\begin{equation}\langle {{\Psi}}^{\text{ab}}\left({\lambda }_{i}\right)\vert {\hat{T}}^{1\text{B}}\vert {{\Psi}}^{\text{ab}}\left({\lambda }_{i}\right)\rangle \approx {T}^{1\text{B}}\left[{\rho }^{\text{ab}}\left({\lambda }_{i}\right)\right].\end{equation} \tag{ 19 }$$

If in the ab initio Hamiltonian the two-body center-of-mass correction ${\hat{T}}^{2\text{B}}$ is included, it needs to be removed from the left-hand side of equation (18). Moreover, we suppose that the Coulomb contribution in the ab initio Hamiltonian is close to the Hartree–Fock average in the functional, namely,

$$\begin{equation}\langle {{\Psi}}^{\text{ab}}\left({\lambda }_{i}\right)\vert {\hat{V}}^{\text{Coul}}\vert {{\Psi}}^{\text{ab}}\left({\lambda }_{i}\right)\rangle \approx {V}^{\text{Coul}}\left[{\rho }^{\text{ab}}\left({\lambda }_{i}\right)\right].\end{equation} \tag{ 20 }$$

Based on these two additional assumptions, we subtract the kinetic and the Coulomb energies from both sides of equation (18), rewriting it as

$$\begin{equation}\langle {{\Psi}}^{\text{ab}}\left({\lambda }_{i}\right)\vert {\hat{V}}^{\text{ab}}\vert {{\Psi}}^{\text{ab}}\left({\lambda }_{i}\right)\rangle =\langle {{\Psi}}^{\text{ab}}\left({\lambda }_{i}\right)\vert {\hat{H}}^{\text{ab}}-\hat{T}-{\hat{V}}^{\text{Coul}}\vert {{\Psi}}^{\text{ab}}\left({\lambda }_{i}\right)\rangle =\sum _{j}{C}_{j}\;{V}_{j}^{\text{gen}}\left[{\rho }^{\text{ab}}\left({\lambda }_{i}\right)\right],\end{equation} \tag{ 21 }$$

where ${\hat{V}}^{\text{ab}}$ is the ab initio potential. The coupling constants C_j can now be obtained by linear regression analysis to best match the ab initio results appearing on both sides of equation (21) and determined for a meaningful set the strength parameters λ_i.

To evaluate the right-hand side of equation (21), we used standard expressions for Hartree–Fock expectation values of two- and three-body generators, that is,

$$\begin{equation}\langle {\Psi}\left({\lambda }_{i}\right)\vert {\hat{V}}_{2\text{B}}^{\text{gen}}\vert {\Psi}{\left({\lambda }_{i}\right)\rangle }_{\text{HF}}=\frac{1}{2}\sum _{\alpha \beta \gamma \delta }\langle \alpha \beta \vert {\hat{V}}_{2\text{B}}^{\text{gen}}\vert \gamma \delta \rangle {\rho }_{\gamma \alpha }^{\text{ab}}\left({\lambda }_{i}\right){\rho }_{\delta \beta }^{\text{ab}}\left({\lambda }_{i}\right),\end{equation} \tag{ 22 }$$

$$\begin{equation}\langle {\Psi}\left({\lambda }_{i}\right)\vert {\hat{V}}_{3\text{B}}^{\text{gen}}\vert {\Psi}{\left({\lambda }_{i}\right)\rangle }_{\text{HF}}=\frac{1}{6}\sum _{\alpha \beta \mu \gamma \delta \nu }\langle \alpha \beta \mu \vert {\hat{V}}_{3\text{B}}^{\text{gen}}\vert \gamma \delta \nu \rangle \;{\rho }_{\gamma \alpha }^{\text{ab}}\left({\lambda }_{i}\right){\rho }_{\delta \beta }^{\text{ab}}\left({\lambda }_{i}\right){\rho }_{\nu \mu }^{\text{ab}}\left({\lambda }_{i}\right),\end{equation} \tag{ 23 }$$

where ${\rho }_{\alpha \beta }^{\text{ab}}$ is the ab initio one-body density matrix, and $\langle \alpha \beta \vert {\hat{V}}_{2\text{B}}^{\text{gen}}\vert \gamma \delta \rangle$ and $\langle \alpha \beta \mu \vert {\hat{V}}_{3\text{B}}^{\text{gen}}\vert \gamma \delta \nu \rangle$ are the antisymmetrized two- and three-body matrix elements.

To evaluate the left-hand side of equation (17), we have to take into account the fact that the SCGF solver is not able to separate the different potential contributions to the Routhian in equation (16), but it provides us with the total interaction energy, defined as

$$\begin{equation}{V}^{\text{tot}}\left({\lambda }_{i}\right)=\langle {{\Psi}}^{\text{ab}}\left({\lambda }_{i}\right)\vert {\hat{V}}^{\text{ab}}+{\hat{V}}^{\text{Coul}}+\sum _{i}{\lambda }_{i}{\hat{V}}_{i}^{\text{gen}}\vert {{\Psi}}^{\text{ab}}\left({\lambda }_{i}\right)\rangle .\end{equation} \tag{ 24 }$$

This gives equation (21) in the form

$$\begin{align}\hfill {V}^{\text{ab}}\left({\lambda }_{i}\right)& \equiv {V}^{\text{tot}}\left({\lambda }_{i}\right)-\langle {{\Psi}}^{\text{ab}}\left({\lambda }_{i}\right)\vert {\hat{V}}^{\text{Coul}}+\sum _{i}{\lambda }_{i}{\hat{V}}_{i}^{\text{gen}}\vert {{\Psi}}^{\text{ab}}\left({\lambda }_{i}\right)\rangle =\sum _{j}{C}_{j}\;{V}_{j}^{\text{gen}}\left[{\rho }^{\text{ab}}\left({\lambda }_{i}\right)\right],\hfill \end{align} \tag{ 25 }$$

that is, we have to evaluate the ab initio expectation values of the Coulomb potential and functional generators.

For a generic two-body operator, the exact expectation value ⟨Ô^2B⟩ ≡ ⟨Ψ(λ_i)|Ô^2B|Ψ(λ_i)⟩ is given as infinite expansion in terms of the effective interaction and dressed propagators, as sketched in figure 1 (first line). From a practical point of view, such a summation becomes computationally difficult, because the dressed propagators contain many poles that multiply matrix elements of every interaction line. Therefore, we approximate the single-particle propagator G_αβ(E) with an optimised reference state (OpRS) propagator [50, 51], defined as

$$\begin{equation}{G}_{\alpha \beta }^{\mathrm{O}\mathrm{p}\mathrm{R}\mathrm{S}}\left(E\right)=\sum _{n\notin F}\frac{{\left({\phi }_{\alpha }^{n}\right)}^{{\ast}}{\phi }_{\beta }^{n}}{E-{{\epsilon}}_{n}^{\mathrm{O}\mathrm{p}\mathrm{R}\mathrm{S}}+\mathrm{i}\eta }+\sum _{k\in F}\frac{{\phi }_{\alpha }^{k}{\left({\phi }_{\beta }^{k}\right)}^{{\ast}}}{E-{{\epsilon}}_{k}^{\mathrm{O}\mathrm{p}\mathrm{R}\mathrm{S}}-\mathrm{i}\eta }.\end{equation} \tag{ 26 }$$

This is a model propagator for independent-particle states with energy ^OpRS and wave function ϕ. Such a propagator contains a reduced number of poles compared to the dressed one, while the energies and wave functions of the OpRS propagator are constrained to give the same momenta of the dressed propagators. We estimate ⟨Ô^2B⟩ using the dressed propagators in the leading order (LO), or Hartree–Fock average, and the OpRS propagators in the next-to-leading order (NLO) and in the all-orders re-summation diagrams, as depicted in the second line of figure 1. We account for higher orders with the ph-, pp- and hh-RPA (random phase approximation) insertions, respectively, in the ring and ladder re-summation diagrams. At the cost of introducing a small error of the density, the use of the OpRS propagator allows us to remarkably speed up calculation of $\langle {\hat{V}}^{\text{Coul}}\rangle$ and $\langle {\hat{V}}_{i}^{\text{gen}}\rangle$, which, otherwise, would not have been possible.

Zoom In Zoom Out Reset image size

Figure 3 shows the ab initio energies in function of λ_i calculated for different perturbations ${\hat{V}}_{i}^{\text{gen}}$ in ³⁴Si. We had expected that the value at λ = 0 (black triangle) would be the one with the lowest energy E^ab(0), because for any state different than the exact ground state |Ψ(0)⟩, the variational principle stipulates that

$$\begin{equation}{E}^{\text{ab}}\left({\lambda }_{i}\right){\geqslant}{E}^{\text{ab}}\left(0\right),\end{equation} \tag{ 27 }$$

assuming, of course, that the ab initio energies are calculated exactly. From the plot, it is evident that there are cases of energies E^ab(λ_i) smaller than E^ab(0) in violation of the variational principle. For the perturbations induced by the three-body generator ${\hat{V}}_{t3}$ (hexagons), the energy does not present a minimum, but increases monotonically with λ_i. This effect can be partially related to the way the contribution of the three-body interaction is extracted. In fact, we can estimate only the leading order as in equation (7).

Zoom In Zoom Out Reset image size

The violation of the variational principle must be traced back to a series of approximations used to calculate the SCGF results. To make any sensible use of these metadata when fitting the functional coupling constants, we thus have to estimate the uncertainties associated with the SCGF calculations.

The major source of error in our calculation is probably related to the subtraction procedure of the perturbation energies from solutions given by the minimisation of the Routhian. For simplicity, we further assume that only one perturbing potential contributes to this uncertainty. The resulting subtraction error, ${\Delta}{E}_{\text{S}}^{\text{ab}}\left({\lambda }_{i}\right)$, associated with the value of E^ab(λ_i) can be estimated as

$$\begin{equation}{\Delta}{E}_{\text{S}}^{\text{ab}}\left({\lambda }_{i}\right)=\left\vert \right.{\langle {\lambda }_{i}{\hat{V}}_{i}^{\text{gen}}\left({\lambda }_{i}\right)\rangle }_{\text{RL}}-{\langle {\lambda }_{i}{\hat{V}}_{i}^{\text{gen}}\left({\lambda }_{i}\right)\rangle }_{\text{NLO}}\left\vert \right.,\end{equation} \tag{ 28 }$$

where ⟨⟩_RL, ⟨⟩_NLO, indicates the truncation, respectively, up to ring and ladder RPA, to NLO approximation. Such error can be viewed as a relative error between the perturbed solutions and the unperturbed one. In figure 4, we show the subtraction errors ${\Delta}{E}_{\text{S}}^{\text{ab}}\left({\lambda }_{i}\right)$ calculated for perturbations induced by the potential ${\hat{V}}_{0}^{J1}$ in ³⁴Si. As expected, this error is zero for λ_i = 0 and grows rapidly with increasing values of the perturbation strength parameter λ_i.

Zoom In Zoom Out Reset image size

We also found another way to estimate uncertainties of the SCGF calculations, which is independent of the approximately calculated average values of the perturbation potentials ${\hat{V}}_{i}^{\text{gen}}$. Indeed, we can determine such estimates by using the well-known Hellmann–Feynman theorem [62]. For any Hamiltonian ${\hat{H}}_{\lambda }={\hat{H}}_{0}+\lambda {\hat{V}}^{\text{pert}}$ that depends explicitly on the parameter λ, the theorem states that

$$\begin{equation}\frac{\mathrm{d}{E}_{\lambda }}{\mathrm{d}\lambda }\equiv \frac{\mathrm{d}}{\mathrm{d}\lambda }\langle {\Psi}\left(\lambda \right)\vert {\hat{H}}_{\lambda }\vert {\Psi}\left(\lambda \right)\rangle =\langle {\Psi}\left(\lambda \right)\vert \frac{\partial {\hat{H}}_{\lambda }}{\partial \lambda }\vert {\Psi}\left(\lambda \right)\rangle =\langle {\Psi}\left(\lambda \right)\vert {\hat{V}}^{\text{pert}}\vert {\Psi}\left(\lambda \right)\rangle .\end{equation} \tag{ 29 }$$

Equation (29) is valid under the condition that |Ψ(λ)⟩ is an eigenstate of Ĥ_λ, Ĥ_λ|Ψ(λ)⟩ = E_λ|Ψ(λ)⟩, or an Hartree–Fock wave function [63], or a variational wave function [64]. However, when the wave function is expanded in a truncated basis [65], or it is a solution of a perturbative expansion [64], the Hellman–Feynman theorem is violated. The ground-state wave function in the SCGF method is not variational because the ADC(3) approximation is a truncated expansion. The degree to which the Hellman–Feynman theorem is violated then illustrates the degree of violation of the variational principle.

The derivative at λ = 0 in the left-hand side of equation (29) can be, for small values of λ, determined by the finite difference method,

$$\begin{equation}\frac{\mathrm{d}{E}_{\lambda }}{\mathrm{d}\lambda }\left(0\right)\approx \frac{E\left(\lambda \right)-E\left(-\lambda \right)}{2\lambda }.\end{equation} \tag{ 30 }$$

In our numerical test, we studied the case of the perturbation given by ${\hat{V}}^{\text{pert}}={\hat{V}}^{\text{ab}}+{\hat{V}}^{\text{Coul}}$, that is, by the full potential V^tot(0), equation (24), that defines Ĥ^ab at λ_i = 0. This compares the finite difference of the energy calculated for the perturbed cases (λ and −λ) with the average value of the interaction energy in the unperturbed case (λ = 0). Such a comparison offers an estimate of the difference between the approximated energy in the ADC(3) method and the exact energy. Consequently, we define the error of the ab initio energy as

$$\begin{equation}{\Delta}{E}_{\text{H--F}}^{\text{ab}}=\left\vert \right.\frac{\mathrm{d}{E}_{\lambda }}{\mathrm{d}\lambda }\left(0\right)-\langle {\Psi}\left(0\right)\vert {\hat{V}}^{\text{tot}}\left(0\right)\vert {\Psi}\left(0\right)\rangle \left\vert \right.,\end{equation} \tag{ 31 }$$

where subscript H–F stands for the error extracted from the Hellmann–Feynman theorem. This error represents an absolute error of the total energy that is due to the approximated solution of the SCGF method. It depends on the nucleus, but we assume it is independent of the perturbation, namely the value calculated at λ = 0 is attributed to all perturbed and unperturbed total energies of a given nucleus.

In the case of N_max = 9 (1bodycm), the violation of the Hellmann–Feynman theorem, as defined in equation (31), is for lighter (heavier) nuclei studied in this paper equal to about 1% (3%–4%) of the total energy. This error is larger than the estimate provided in reference [37], because it gives a cumulative effect resulting from the ADC(3) truncation and reduced model space.

In figure 4, we also show the Hellmann–Feynman error ${\Delta}{E}_{\text{H--F}}^{\text{ab}}$ determined in ³⁴Si. We see that for this particular generator, ${\hat{V}}_{0}^{J1}$, ${\Delta}{E}_{\text{H--F}}^{\text{ab}}$ is more than twice larger than ${\Delta}{E}_{\text{S}}^{\text{ab}}$. We checked that in most cases considered in this paper ${\Delta}{E}_{\text{H--F}}^{\text{ab}}{ >}{\Delta}{E}_{\text{S}}^{\text{ab}}$, even in very light nuclei. Given the very different magnitude and sources of the subtraction and Hellmann–Feynman errors, we decided to keep them separate and perform two independent analyses of the coupling constants.

In view of the identified uncertainties, we can now conclude that the explicit violation of the variational principle obtained for E^ab(λ_i ≠ 0), see figure 3, can be considered acceptable consequences of the inherent imprecision encountered in the SCGF approach.

From the set of ab initio calculations described in section 2.3, we obtained d = 284 equations (25) whose left-hand sides represent regression dependent variable y_i, i = 1–284, and whose p expectation values on the right-hand sides form the matrix of features, that is,

$$\begin{equation}{y}_{i}=\sum _{j=1}^{p}{\mathcal{J}}_{ij}{C}_{j},\end{equation} \tag{ 32 }$$

where each coupling constant C_j corresponds to a generator of the model functional given in equation (8).

Introducing for compactness the vector notation: C = (C₁, ..., C_p) and y = (y₁, ..., y_d), we build the penalty function by of the least-square method as

$$\begin{equation}{\chi }^{2}\left(\boldsymbol{C}\right)=\frac{1}{d-p}{\left(\mathcal{J}\boldsymbol{C}-\boldsymbol{y}\right)}^{T}W\left(\mathcal{J}\boldsymbol{C}-\boldsymbol{y}\right),\end{equation} \tag{ 33 }$$

where the weight matrix W is a diagonal matrix with elements W_ii = w_i. We define the weight of each data point as the inverse of the estimated error squared,

$$\begin{equation}{w}_{i}=\frac{1}{{\left({\Delta}{y}_{i}\right)}^{2}},\end{equation} \tag{ 34 }$$

where Δy_i is composed of three contributions [66],

$$\begin{equation}{\left({\Delta}{y}_{i}\right)}^{2}={\left({\Delta}{y}_{i}^{\text{ab}}\right)}^{2}+{\left({\Delta}{y}^{\text{num}}\right)}^{2}+{\left({\Delta}{y}^{\text{mod}}\right)}^{2}.\end{equation} \tag{ 35 }$$

${\Delta}{y}_{i}^{\text{ab}}$ is the error attributed to the SCGF approach, namely, we take it as ${\Delta}{E}_{\text{S}}^{\text{ab}}$ (28) or ${\Delta}{E}_{\text{H--F}}^{\text{ab}}$ (31). Δy^num is the numerical precision of the SCGF calculations, which is smaller that 5 × 10⁻⁵  MeV and can be neglected. Δy^mod represents the error associated with the model itself, which is entirely unknown, and thus it has to be tuned to normalise the penalty function χ². Starting from an arbitrary value, Δy^mod can be increased iteratively up to the value at which the χ² approaches the value of 1. Then, the penalty function in equation (33) satisfies the typical statistical normalisation condition χ²(C) → 1 at the minimum, and Δy^mod acquires interpretation of the Birge factor [67].

Minimisation of χ² gives the solution C_min, covariance matrix $\mathcal{K}$, statistical error associated with the parameters, ΔC_min, and propagated errors of observables [66, 68, 69].

We minimise the χ² defined in equation (33) using the two types of errors discussed in previous section. We thus obtain two sets of results: the one labelled D_S(10), that is, obtained using the errors defined in equation (28) and that labelled D_H–F(10)—using the errors defined in equation (31).

Given the large uncertainties of the metadata and the fact that they could be obtained only in fairly light nuclei with small isospin asymmetry, the resulting coupling constants are poorly determined with quite large error bars. This is particularly evident for the D_H–F(10) set of parameters. The relative errors of the isovector coupling constants are often of the order of 100%, meaning that the obtained values are compatible with 0. The very large errors mean that the χ² surface is fairly flat in these particular directions of the parameter space. As a consequence, when used to calculate atomic nuclei, the obtained coupling constants often immediately lead to finite-size instabilities [70].

The quality of the fit can be easily judged by inspecting relative residuals of equation (25), defined as

$$\begin{equation}\mathrm{R}\mathrm{e}\mathrm{s}\mathrm{i}\mathrm{d}\mathrm{u}\mathrm{a}\mathrm{l}\mathrm{s}\equiv \frac{{\sum }_{j}{C}_{j}\;{V}_{j}^{\text{gen}}\left[{\rho }^{\text{ab}}\left({\lambda }_{i}\right)\right]-{V}^{\text{ab}}\left({\lambda }_{i}\right)}{{V}^{\text{ab}}\left({\lambda }_{i}\right)},\end{equation} \tag{ 36 }$$

and for ⁵⁶Ni shown in figure 5(a) for D_S(10) and (b) for D_H-F(10). For a good-quality fit, we would expect the residuals lie close to the dashed horizontal line, which indicates zero values. Rapid departures of residuals from zero, especially for the second-order isovector coupling constants, illustrate the fact that a reasonable description of the ab initio results in terms of Skyrme functional generators could not be obtained.

Zoom In Zoom Out Reset image size

For all nuclei and generators studied here, the residuals corresponding to the unperturbed systems (λ = 0) are around zero, whereas when λ_i moves away from zero, the residuals increase. In addition, they are not normally distributed, but they exhibit clear trends, meaning that the hypothesis formulated in equation (25) is not correct. Shadows shown in figure 5(a) represent the propagated errors of the corresponding observables obtained in the fit D_S(10). For D_H–F(10), the propagated errors are much larger, so in figure 5(b), for clarity we only show the one corresponding to generator ${\hat{V}}_{W0}$. We note that the propagated errors related to the violation of the Hellmann Feynman theorem are much larger than those corresponding to the subtraction errors, and thus they make it very hard to draw reasonable conclusions about the structure of the residuals.

We tested the performance of the obtained coupling constants in the description of infinite nuclear matter. For D_S(10), we obtain a fairly reasonable energy per particle of E/A = −14.6 ± 0.2 MeV and saturation density of ρ₀ = 0.132 ± 0.002 fm⁻³. The D_H–F(10) provides slightly different values for the binding energy E/A = −16.0 ± 1.2 MeV and saturation density ρ₀ = 0.128 ± 0.01 fm⁻³, but with much larger error bars than in D_S(10). For both sets, the nuclear incompressibility K is also in the range of acceptable values [71]; in particular, we have K = 380 ± 17 MeV for D_S(10) and K = 393 ± 143 MeV for D_H–F(10). Other properties of infinite matter, such as the symmetry energy and its first derivative, are poorly determined, probably because the set of nuclei used for the fit is not rich enough to describe isovector properties of the nuclear medium sufficiently well. In particular, we find that the symmetry energy has a value compatible with zero within the error bars.

Another major drawback of the derived coupling constants is an unrealistic value of the effective mass. Although the effective mass is not strictly speaking an observable, we can extract information about its value from other many-body methods [72]. For the isoscalar effective mass, an acceptable range of values is m^*/m ∈ [0.7–0.9], although it is worth mentioning that also other values are found in the literature. Both sets of coupling constants, D_S(10) and D_H–F(10), give effective masses that are off by roughly an order of magnitude (cf figure 7), and this is probably the main reason why they do not lead to realistic results when used in calculations of finite nuclei.

Since a simple least-square minimisation does not provide us with satisfactory values of the effective mass, we also performed a constrained linear regression to drag the coupling constants towards reasonable values of m*/m. In this way, we want to test if the poor determination of the coupling constants, reflected in their large errors, can be exploited to improve values of the effective mass. Linear regression with constraints is a procedure in the spirit of the Bayesian inference, where a prior information about the parameters is known. Here we consider it in the form of the Tikhonov regularisation [73] or ridge regression [74], which consists of minimising the penalty function of the form

$$\begin{equation}{\chi }_{\text{T}}^{2}\left(\boldsymbol{C}\right)=\frac{1}{d-p+f}\left[{\left(\boldsymbol{y}-\mathcal{J}\boldsymbol{C}\right)}^{T}W\left(\boldsymbol{y}-\mathcal{J}\boldsymbol{C}\right)+{\lambda }_{\text{T}}{\left(\boldsymbol{b}-\boldsymbol{Q}\left[\boldsymbol{C}\right]\right)}^{T}\left(\boldsymbol{b}-\boldsymbol{Q}\left[\boldsymbol{C}\right]\right)\right].\end{equation} \tag{ 37 }$$

The Tikhonov parameter λ_T is a real positive number and b = Q[C] represents a system of f linear equations in parameters C with constant terms b. The final values of C will depend on λ_T. As one can see, equation (37) is defined as a sum of two terms: ${\chi }_{\text{data}}^{2}$, which depends on weights W, and ${\chi }_{\text{constr}}^{2}$, which depends on λ_T. Increasing the value of the Tikhonov parameter gives more importance to ${\chi }_{\text{constr}}^{2}$ and eventually may deteriorate the description of data represented by ${\chi }_{\text{data}}^{2}$.

We used one constraint, namely the definition of the in-medium isoscalar effective mass,

$$\begin{equation}\frac{{m}^{{\ast}}}{m}\left({\rho }_{\text{sat}}\right)={\left[1+\frac{2\;m}{{\hslash }^{2}}{C}_{0}^{\tau }{\rho }_{\text{sat}}\right]}^{-1}.\end{equation} \tag{ 38 }$$

It is important to note that $\frac{{m}^{{\ast}}}{m}$ depends explicitly only on one coupling constant, ${C}_{0}^{\tau }$, however, implicitly, through the saturation density ρ_sat it non-linearly depends on all other coupling constants. We set the target value of the effective mass to b₁ ≡ m*/m = 0.70 and in the Tikhonov term we varied the value of log₁₀λ_T from −4 to 2. As an example, in ${\chi }_{\text{data}}^{2}$ we took the inputs used for D_H–F(10). By implementing the Tikhonov regularisation, we aim to distinguish between two possible scenarios, namely, whether at an expense of a small deterioration of ${\chi }_{\text{data}}^{2}$, description of effective mass can or cannot be improved.

In figure 6(a), we show the obtained evolution of the constrained parameters C in function of the Tikhonov parameter. As one can see, changes in the parameters happen around log₁₀λ_T ≈ −2. Beyond this region, their values are quite stable. ${C}_{1}^{{\Delta}\rho }$, ${C}_{0}^{J1}$, ${C}_{1}^{J1}$ and C_t3 are the coupling constants to which the constraint makes the largest impact. We already noted that these coupling constants were poorly determined by the unconstrained regression. In figure 6(b), we show relative contributions to the penalty function coming from the data points, ${\chi }_{\text{data}}^{2}$, and from the constraint, ${\chi }_{\text{constr}}^{2}$, (see equation (37)). We see that the contribution of the constraint increases with λ_T up to log₁₀λ_T ≈ −1.8 and after that it quickly drops to zero. At the peak of ${\chi }_{\text{constr}}^{2}$, the coupling constants begin to adjust to the requested value of m*/m. A further increase of λ_T decreases ${\chi }_{\text{constr}}^{2}$, because b₁ − Q₁[C] → 0, whereas the reproduction of the data points deteriorates and the differences $\boldsymbol{y}-\mathcal{J}\boldsymbol{C}$ increase.

Zoom In Zoom Out Reset image size

In figure 7, we present nuclear matter properties ρ_sat, E/A, m*/m, and symmetry energy J obtained by the Tikhonov regularisation of the D_H–F(10) coupling constants with m*/m constrained to 0.70. For small λ_T, the obtained values are equal to those corresponding to the original D_H–F(10) results. In the region of log₁₀λ_T between −2 and 0, nuclear-matter properties change abruptly, and effective mass is dragged towards the target value already at log₁₀λ_T ≃ −1. In the region of log₁₀λ_T ⩽ 0, nuclear-matter properties ρ_sat, E/A, and J cross their empirical values of ρ_sat = 0.17 ± 0.03  fm⁻³, E/A = −16 ± 1  MeV, and J = 32.5 ± 2.5  MeV (shown in grey) [71]. Such a crossing occurs at slightly different values of λ_T for the three quantities and before the effective mass reaches its empirical range. We note here, that for the NNLO_sat chiral interaction, the corresponding nuclear-matter values read ρ_sat = 0.166  fm⁻³ and E/A = −14.5  MeV [30]. Beyond log₁₀λ_T ≃ 0, nuclear-matter properties settle at values far away from the standard nuclear-matter values. When the effective mass approaches the target value of 0.70, the energy per particle E/A becomes too low and the symmetry energy becomes very large. We conclude that there is no region of the Tikhonov parameter where all four nuclear-matter quantities would be in their expected domains, even when the propagated errors (represented by shaded areas) are taken into account. In other words, the obtained incorrect values of the effective mass cannot be corrected by a small deterioration of the penalty function ${\chi }_{\text{data}}^{2}$.

Zoom In Zoom Out Reset image size