Model nuclear energy density functionals derived from ab initio calculations

We present the first application of a new approach, proposed in [Journal of Physics G: Nuclear and Particle Physics, 43, 04LT01 (2016)] to derive coupling constants of the Skyrme energy density functional (EDF) from ab initio Hamiltonian. By perturbing the ab initio Hamiltonian with several functional generators defining the Skyrme EDF, we create a set of metadata that is then used to constrain the coupling constants of the functional. We use statistical analysis to obtain such an ab initio-equivalent Skyrme EDF. We find that the resulting functional describes properties of atomic nuclei and infinite nuclear matter quite poorly. This may point out to the necessity of building up the ab initio-equivalent functionals from more sophisticated generators. However, we also indicate that the current precision of the ab initio calculations may be insufficient for deriving meaningful nuclear EDFs.


Introduction
Different energy scales that appear in nuclear systems suggest a theoretical approach based on effective field theories (EFT), which use relevant degrees of freedom adapted to a given energy scale [1]. A remarkable example is the chiral effective field theory (χ-EFT) [2,3]: by using neutrons, protons, and pions as degrees of freedom, χ-EFT is able to provide consistent description of numerous observables in atomic nuclei. Chiral interactions are used in the framework of ab initio methods to solve the many-body Schrödinger equations, employing controlled approximations. The numerical solutions require huge computational resources; for this reason, actual state-of-the-art ab initio calculations are limited to nuclei that are not too heavy and/or close to (semi)magic systems [4,5,6,7]. In particular, the description of open-shells nuclei has a recent history, started with the Gorkov-Green's Function approach [8].
Another example of successful effective theory is the approach based on nuclear energy density functionals (EDFs) [9, 10]. Nuclear EDF is a versatile tool that, at low computational cost, allows us to describe properties of atomic nuclei across the entire nuclear chart from drip-line to drip-line and from light to super-heavy nuclei. The underlying non-relativistic functionals are usually obtained from phenomenological two-body potentials called functional generators, which have radial form factors of zerorange, for Skyrme [11,12], or finite-range, for Gogny [13] implementations, with coupling constants adjusted to reproduce selected nuclear observables. In addition to finite nuclei, also infinite nuclear matter properties can be addressed within the EDFs formalism [14, 15,16], and included among the observables.
In Ref. [17], it was shown that Skyrme EDFs have reached their best in terms of reproducing the experimental observables. As a consequence, several groups are exploring new forms of functional generators to improve precision of describing data [18,19,20,21,22]. The resulting functionals are then adjusted following standard procedures based on constraining the coupling constants to a given set of nuclear observables [23,24].
In this article, we aim to explore a complementary approach, which is based on explicitly bridging the ab initio methods with the nuclear EDFs. Formal connections between the two approaches are the subject of recent studies [25,26]. Furthermore, in Ref. [27] the density matrix expansion formalism has been used to obtain a set of functionals microscopically constrained from chiral effective field theory interactions. Here, we proceed by fitting parameters that enter the energy functional directly to metadata generated by the ab initio calculations, using the method suggested in Ref. [28]. In Ref. [28] this method was applied adjusting the Skyrme coupling constants to metadata generated by the Gogny functional. In this work, we perform a realistic application of the method by employing full-fledged ab initio calculations that use Self-Consistent Green's Function (SCGF) [29] methodology with the chiral interaction NNLO sat [30].
The paper is organised as follows: Section 2 briefly recalls the SCGF (2.1) and EDF (2.2) methods, and explains the formalism to derive the model functionals (2.3). Results are discussed in Section 3 and conclusions are given in Section 4.

Self-Consistent Green's Function method
We separate the nuclear many-body HamiltonianĤ =Ĥ 0 +Ĥ 1 into a noninteracting partĤ 0 =T +Û, with the auxiliary one-body operatorÛ , and the interacting part defined asĤ 1 = −Û +V 2B +V 3B . We express the nuclear Hamiltonian in second quantisation aŝ where ǫ 0 α are the single-particle energies ofĤ 0 . By solving the corresponding Schrödinger equation,Ĥ one obtains the ground-state energy E A 0 and wave function |Ψ A 0 of the nuclear system. The exact solution of Eq. (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 propagator of one-body Green's function G αβ . Using the Källén-Lehmann representation [32,33], G αβ takes the form Here the complete set of eigenstates |Ψ A+1 are respectively the excitation energies of the propagating quasi-particle (index n) and quasi-hole (index k) states. The propagator given in Eq. (3) satisfies the Dyson equation where G αβ (E) is the propagator for the noninteracting systemĤ 0 , and Σ ⋆ γδ (E) 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 Eq. (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 Eq. (4) is the self-energy. This is composed of three parts as 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 Ref. [37] for a more pedagogical introduction of the method. The accuracy of ADC (3) can be estimated to be of the fourth order inĤ 1 , 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 interactionV 3B . By means of the modified Migdal-Galitski-Koltun sum rule [38], we can write the ground-state energy of the system as, where ǫ − 0 is the highest quasi-hole energy, namely ). 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 where ρ αβ is the one-body density matrix.

Model Energy Density Functionals and generators
Having defined the main theoretical aspects of the ab initio method employed in this article, we now define the model energy density functional as [ρ] ≡ Φ|V gen j |Φ HF . The operatorsV gen j 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 generatorsV gen j , based on ten individual termsT i andT σ i for i = 0, 1, 2, andT e ,T o ,T W 0 , andT 3 of the Skyrme functional generator [12], that is,V where A 123 is the antisymmetric operator for combinations ofP σ = 1 2 (1 + σ 1 · σ 2 ) and P τ = 1 2 (1 + τ 1 · τ 2 ). We refer to Ref. [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 Eqs. (9) to linear combinations that give specific isoscalar/isovector terms of the functional [28]. This is done by means of the following matrix relation, We explicitly included only the vector partV J1 T 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 theV W 0 andV 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 Eq. (9). The generators listed on the left-hand side of Eq. (10) are then used as generators,V gen j , j = 1 − 10, that define the functional in Eq. (8).

Derivation of the functionals
In the electronic density functional theory [41], one uses the Levy-Lieb constrained variation [42,43,44,45] to obtain the ground state energy of the system. This procedure consists in a two-step minimisation, whereV stands for the Coulomb potential and symbol min |Ψ →ρ denotes an inner (firststep) 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 RouthianR at fixed one-body external potential U(r) that plays a role of the Lagrange multiplier, This gives the energy as functionals of the potential U(r). Assuming that the inverse functional U[ρ] can be found, we obtain the exact energy density functional, 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 RouthianR ab as whereĤ ab is the Hamiltonian of the system, Eq. (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 Eq. (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 ab g.s. 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, 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 finitedimensional 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 secondstep variation would now correspond to the minimisation of function E ab (λ i ) in finite dimensions.
As proposed in Ref. [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, Sec. 2.2, that is, Since the unrestricted minimisation of the Routhian is equivalent to finding its exact ground state |Ψ ab (λ i ) and eigenvalue R ab (λ i ), we have Formally, by replacing minimisation (15) with minimisation (16), we do not make any new assumption. 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, Considering that the Kohn-Sham kinetic energy represents a good approximation to the one-body kinetic energy, we assume that If in the ab initio Hamiltonian the two-body center-of-mass correctionT 2B is included, it needs to be removed from the left-hand-side of Eq. (18). Furthermore, we suppose that the Coulomb contribution in the ab initio Hamiltonian is close to the Hartree-Fock average in the functional, namely, With these two additional assumptions, we subtract the kinetic and the Coulomb energies from both sides of Eq. (18), rewriting it as whereV 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 Eq. (21) and determined for a meaningful set the strength parameters λ i . To evaluate the right-hand side of Eq. (21), we used standard expressions for Hartree-Fock expectation values of two-and three-body generators, that is, where ρ ab αβ is the ab initio one-body density matrix, and αβ|V gen 2B |γδ and αβµ|V gen 3B |γδν are the antisymmetrized two-and three-body matrix elements. To evaluate the left-hand side of Eq. (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 Eq. (16), but it provides us with the total interaction energy, defined as This gives Eq. (21) in the form 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 Fig. 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 [46,47], defined as 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 resummation diagrams, as depicted in the second line of Fig. 1. We account for higher orders with the ph-, pp-and hh-RPA (Random Phase Approximation) insertions respectively in the ring and ladder resummation 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 V Coul and V gen i , which, otherwise, would not have been possible.

Results
In this section, we present our results obtained using SCGF with ADC(3) approximation. We chose the interaction NNLO sat [30] for the two-body and three-body sector. We made this choice since this interaction has been optimised to reproduce ground state energies and radii for isotopes up to mass A = 24 and it has been shown to predict accurate saturation properties also for larger isotopes and for infinite nuclear matter [48,49,50,51]. 16  The SCGF calculations were performed using a basis of spherical harmonic oscillator wave functions with oscillator energy ω=20 MeV. The model space was limited to states with principal quantum numbers smaller or equal than N max . In terms of computational resources, even if the present technology allows calculations up to N max = 13, this would require a too large amount of CPU hours to complete the full set of perturbations required in the current analysis. We thus decided to limit the model space to N max = 9. This model space may converge poorly to the experimental values, but, in the first attempt at this approach, our attention focus on the validity of the method rather than on the accurate prediction of the experimental masses.
In Fig. 2 Values of E/A are reported for three choices of the size of model space: N max = 7, 9, and 11, to be compared with the experimental values taken from Ref. [52]. The convergence in terms of N max is not fully achieved; nevertheless, with increasing N max , most binding energies decrease towards the experimental values. The model space of N max = 7 appears to be too small to reproduce experimental results. For N max = 11, the maximum difference with experiment appears for 56 Ni, with the difference of around 3% of the total energy. This discrepancy is a combination of the 1% error on the correlation energy that is associated with the ADC(3) truncation [37] and of the accuracy of stateof-the-art chiral nuclear forces [53,51]. In absolute terms, the calculated binding energy of 56 Ni is more than 10 MeV above the measured value. Such a deviation is much larger than the standard deviation of typical EDFs, which is of the order of 1 MeV [24,54]. Clearly, with the model space truncation discussed below, the overall accuracy with respect to predicting the experiment is too poor to obtain novel functionals capable of reducing the discrepancies between the EDF approach and the experiment. However, our goal is to reproduce the ab initio energies, irrespective of their detailed agreement with experiment. Typical ab initio computations subtract the kinetic energy of the center of mass to directly access the intrinsic ground state energy [55]. This implies adding a one-and a two-body correction terms to the nuclear Hamiltonian [56,37]. The result of retaining only the one-body term is shown in Fig. 2 for N max = 9 as "1bodycm" and it amounts to an overcorrection. However, this additional discrepancy will not affect our goal of investigating the consistence between microscopic results and the functionals that they generate. For our purposes, it is technically simpler to drop the two-body center-ofmass correction because it avoids further approximations in the way the kinetic energy is treated in the SCGF and EDF approaches. Since the EDF results also depend on the number of oscillator shells included in the model space [57], we fixed the model spaces of both approaches to be N max = 9 (1bodycm) for the following analysis.
In the model functional, Eq. (8), we consider generators of zero-range contact interaction (Skyrme-like) defined in Eq. (10), namely,V ρ T ,V ∆ρ T ,V τ T ,V J1 T ,V W 0 , and V t3 . Subscripts T = 0 and 1 denote generators inducing terms of the functional that depend on isoscalar and isovector densities, respectively. Such a link between generators and densities is valid only on the EDF level; for the chiral interactions used in our SCGF computations, already at the next-to-leading order in powers of the interaction, expectation values of generators contain contributions from both isoscalar and isovector densities.
We perturbed ground states of each of the seven selected nuclei using four different intensities of the perturbation strength λ i , separately for each of the 10 generators. In this way, we obtained 284 converged results ‡, which represent our full database of the perturbed and unperturbed ground-state energies. The choice of using non-zero values of λ i separately for eachV gen i represents a compromise between the volume of calculations and a coverage of the full manifold in the space of perturbations λ i . Obviously, the ‡ For 36 S, 3 data points did not converge.
larger the λ i 's the wider is the density space probed, however, if perturbations are too strong, the numerical SCGF solutions may diverge. In practice, for each generator, we have found a most suitable range of values of λ i 's that were used in the final calculations.    34 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

Estimated errors on the SCGF calculations
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 generatorV 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 Eq. (7). 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, ∆E ab S (λ i ), associated with the value of E ab (λ i ) can be estimated as 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 ∆E ab S (λ i ) calculated for perturbations induced by the potentialV J1 0 in 34 Si. As expected, this error is zero for λ i = 0 and grows rapidly with increasing values of the perturbation strength parameter λ i .
We also found another way to estimate uncertainties of the SCGF calculations, which is independent of the approximately calculated average values of the perturbation potentialsV gen i . Indeed, we can determine such estimates by using the well-known Hellmann-Feynman theorem [58]. For any HamiltonianĤ λ =Ĥ 0 + λV pert that depends explicitly on the parameter λ, the theorem states that Eq. (29) is valid under the condition that |Ψ(λ) is an eigenstate ofĤ λ ,Ĥ λ |Ψ(λ) = E λ |Ψ(λ) , or an Hartree-Fock wave function [59], or a variational wave function [60]. However, when the wave function is expanded in a truncated basis [61], or it is a solution of a perturbative expansion [60], 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 Eq. (29) can be, for small values of λ, determined by the finite difference method, In our numerical test, we studied the case of the perturbation given byV pert = V ab +V Coul , that is, by the full potential V tot (0), Eq. (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 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 Eq. (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 Ref. [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 ∆E ab H-F determined in 34 Si. We see that for this particular generator,V J1 0 , ∆E ab H-F is more than twice larger than ∆E ab S . We checked that in most cases considered in this paper ∆E ab H-F > ∆E ab S , 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 Fig. 3, can be considered acceptable consequences of the inherent imprecision encountered in the SCGF approach.

Linear regression analysis to determine the coupling constants
From the set of ab initio calculations described in Sec. 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, where each coupling constant C j corresponds to a generator of the model functional given in Eq. (8).
Introducing for compactness the vector notation: C = (C 1 , . . . , C p ) and y = (y 1 , . . . , y d ), we build the penalty function by of the least-square method as 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, where ∆y i is composed of three contributions [62], ∆y ab i is the error attributed to the SCGF approach, namely, we take it as ∆E ab S (28) or ∆E ab H-F (31). ∆y num is the numerical precision of the SCGF calculations, which is smaller that 5×10 −5 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 χ 2 . Starting from an arbitrary value, ∆y mod can be increased iteratively up to the value at which the χ 2 approaches the value of 1. Then, the penalty function in Eq. (33) satisfies the typical statistical normalisation condition χ 2 (C) → 1 at the minimum, and ∆y mod acquires interpretation of the Birge factor [63].
Minimisation of χ 2 gives the solution C min , covariance matrix K, statistical error associated with the parameters, ∆C min , and propagated errors of observables [64,62,65].

Fitted coupling constants
We minimise the χ 2 defined in Eq. (33) using the two types of errors discussed in previous section. We thus obtain two sets of results: the one labeled D S (10), that is, obtained using the errors defined in Eq. (28) and that labeled D H-F (10) -using the errors defined in Eq. (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 χ 2 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 [66].
The quality of the fit can be easily judged by inspecting relative residuals of Eq. (25), defined as and for 56 Ni shown in Fig. 5(a) for D S (10) and 5(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. 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 Eq. (25) is not correct. Shadows shown in Fig. (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 Fig. (5)(b), for clarity we only show the one corresponding to generatorV W 0 . 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 saturation properties of infinite nuclear matter, with an energy per particle E/A = −14.6 ± 0.  a value for the saturation density of ρ 0 = 0.132 ± 0.002 fm −3 . The D H-F (10) provides slightly different values for the binding energy E/A = −16.0 ± 1.2 MeV and saturation densityρ 0 = 0.128 ± 0.01 fm −3 , 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 [67]; 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 [68]. 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. Fig. 7), and this is probably the main reason why they do not lead to realistic results when used in calculations of finite nuclei.

Constraints on the nuclear matter properties
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 toward 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 [69] or ridge regression [70], which consists of minimising a penalty function of the form 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 , which is unique for the f independent constraining equations. Eq. (37) is defined as a sum of two terms: χ 2 data , which depends on weights W , and χ 2 constr , which depends on λ T . Increasing the value of the Tikhonov parameter gives more importance to χ 2 constr and eventually may deteriorate the description of data represented by χ 2 data . We use one constraint, namely the definition of the in-medium isoscalar effective mass, It is important to note that m * m depends explicitly only on one coupling constant, C τ 0 , 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 1 ≡ m * /m=0.70 and in the Tikhonov term we varied the value of log 10 λ T from -4 to 2. As an example, in χ 2 data we took the inputs used for D H-F (10).  In Fig. 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 10 λ T ≈ −2. Beyond this region, their values are quite stable. C ∆ρ 1 , C J1 0 , C J1 1 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 Fig. 6(b), we show separate contributions from the data points, χ 2 data , and from the constraint, χ 2 constr , (see Eq. (37)). We see that the contribution of the constraint increases with λ T up to log 10 λ T ≈ -1.8 and after that it quickly drops to zero. At the peak of χ 2 constr , the coupling constants begin to adjust to the requested value of m * /m. A further increase of λ T decreases χ 2 constr , because b 1 − Q 1 [C] → 0, whereas the reproduction of the data points deteriorates and the differences y − J C increase.
In Fig. 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 10 λ T between -2 and 0, nuclear-matter properties change abruptly, and effective mass is dragged towards the target value already at log 10 λ T ≃ −1. Beyond log 10 λ T ≃ 0 nuclear-matter properties do not vary much. The curves for ρ sat , E/A, and symmetry energy cross their regions of empirical values (shown in grey) before setting to values far away from the standard nuclear-matter values. Such crossing occurs at slightly different values of λ T for the three quantities and before the effective mass reaches its empirical range. 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. A possible reason for this failure can probably be traced back to strong correlations between the Skyrme nuclear-matter properties, cf. Ref. [71].

Conclusions
Applying the methodology suggested in Ref. [28], we studied the link between the nuclear Skyrme functional and the NNLO sat chiral interaction used within the ab initio Self-Consistent Green's Function calculations in ADC(3) approximation. We performed the ab initio calculations in seven light closed shell nuclei by perturbing their ground states with ten functional generators that define Skyrme functional. By employing the linear regression method, the obtained metadata were used to derive the functional coupling constants. We analysed two possible sources of uncertainties of the ab initio calculations: the first one related to approximate determination of average values of two-body potentials and the second one to an imprecise determination of nuclear ground states, for which we employed the Hellmann-Feynman theorem.
The obtained values of Skyrme coupling constants were very different than those typically obtained in phenomenological adjustments to nuclear observables. We have identified several possible reasons of such a result: First, it appears that a relatively high level of uncertainties arising in the ab initio calculations induces large uncertainties of the derived coupling constants, which then propagate to large uncertainties of the nuclear matter properties and to instabilities when solving the self-consistent equations.
Second, it appears that the ab initio energies are poorly reproduced by the terms in the functional generated by the Skyrme zero-range potentials. Third, it appears that the information contents within the perturbed ground-state energies of light semi-magic nuclei is insufficient to properly determine Skyrme coupling constants, especially those corresponding to second-order terms depending on isovector densities.
Certainly, future research may be focused on applications of ab initio technologies with improved overall precision, which would better correspond to the ambition of reducing discrepancies between the phenomenological EDF results and experiment. Another promising avenue would be to repeat present analyses by using finite-range functional generators. However, the most challenging aspect of the methodology proposed in Ref. [28] is the fact that it is based, similarly as most other generic ab initio DFT approaches, on the Levy-Lieb construction that essentially pertains to variational studies of ground-state energies. This is in opposition to the methodology of adjusting functionals directly to experimental data, where one uses not only ground-state energies, but also other essential observables like radii, deformations, or transition probabilities.