Optimized nuclear energy density functionals including long-range pion contributions

Nuclear energy density functionals successfully reproduce properties of nuclei across almost the entire nuclear chart. However, nearly all available functionals are phenomenological in nature and lack a rigorous connection to systematically improvable nuclear forces. This issue might be solved with an energy density functional obtained from first principles. As an intermediate step towards this goal we construct the GUDE family of functionals that is obtained from a hybrid scheme consisting of long-range pion-exchange contributions derived from chiral effective field theory at the Hartree-Fock level and a phenomenological Skyrme part. When including pion contributions beyond next-to-leading order in the chiral expansion, we find significant improvements over a reference Skyrme functional constructed following the same protocol. We analyze the importance of different pion contributions and identify which terms drive the observed improvements. Since pions are incorporated without adding further optimization parameters to the functionals, the improvements can be attributed to the functional form of these terms. Our work therefore suggests that the considered chiral contributions constitute useful ingredients for true ab initio energy density functionals.


I. INTRODUCTION
Tremendous progress has been made in calculating nuclear structure from first principles [1,2], pushing descriptions toward heavy [3][4][5] and doubly open-shell [6][7][8][9][10][11] nuclei, and employing high-precision interaction models [12][13][14] and high-order many-body methods [15][16][17].However, due to their huge numerical cost, these microscopic approaches, usually generically referred to as ab initio methods [18], are not yet ready to be employed in large-scale, high-precision calculations of nuclear groundstate observables.Even if one could overcome this computational challenge, it is unclear whether ab initio calculations are going to be able to compete with less microscopic methods regarding the accuracy they can achieve.At present, they generally cannot [1,2,15,19].
Nuclear density functional theory (DFT) [20,21] is currently the most microscopic theoretical framework that can be used in global surveys thanks to its favorable computational scaling [22].It is rooted in the seminal work by Hohenberg and Kohn proving the existence of a universal functional of the density which, when minimized for fixed particle number, gives the ground state density and energy of a many-body system confined in an external potential [23].While this is most commonly employed for the description of electronic systems, later works extended the existence proof to self-bound systems as constituted by finite nuclei [24][25][26][27].In practice most calculations are carried out in the Kohn-Sham formulation of DFT [28], which allows for an efficient description of the kinetic energy of the system and of shell effects by expressing the density of interest in terms of auxiliary single-particle orbitals of an independent-particle system.
In nuclear physics different ansatze have been established for the form of the energy density functional (EDF).In the nonrelativistic sector, the Skyrme [29] and Gogny [30] EDFs are based on effective nucleon-nucleon interactions.Genuine energy functionals (not derived from an underlying potential) include the Fayans [31], the SeaLL1 [32], and the BCPM [33] functionals.Different forms are also available in covariant DFT; see, e.g., Refs.[34,35].Here we will limit ourselves to nonrelativistic functionals.
Significant progress in nuclear DFT has been achieved by using increasingly sophisticated parameter optimization protocols but it is widely believed that this avenue has been explored to such a degree that further improvements, necessary for instance for the description of r-process nucleosynthesis [36][37][38][39][40] or of singleparticle energies [41], need to come from elsewhere [41][42][43].The two most obvious routes are the explicit treatment of static correlations within a multireference framework [21,30,44,45] and the extension of the form of the employed EDFs.
In the latter direction, different empirical strategies have been pursued (see, e.g., Refs.[46][47][48][49][50][51][52]).They often consist in adding similar or higher-order terms to existing EDF structures and typically involve introducing additional adjustable parameters.Properly fitting such parameters is a nontrivial task since they cannot always be well constrained with available experimental data.This does not address the phenomenological nature of the EDFs, which is the root cause for potentially uncontrolled extrapolations outside the fitting regions [43,[53][54][55][56].
A unifying construction principle for nuclear EDFs might therefore be helpful.While different ideas to formulate an effective field theory (EFT) for EDFs have been discussed [57][58][59], none of them has been implemented yet.Alternatively, one can remain within the overall framework of nuclear DFT while seeking guidance from microscopic ab initio theories.By employing interactions derived from chiral EFT, which establishes a construction scheme based on a power counting estimating the importance of individual terms [60,61], ab initio calculations become systematically improvable by going to higher orders in the chiral expansion.At present the most accurate potentials are constructed at fifth order for nucleon-nucleon (NN) forces [12,13] and fourth order for three-nucleon (3N) forces [2].Different ideas exist for how to combine ab initio approaches and nuclear DFT [19,22,33,[62][63][64][65][66][67][68][69].They range from determining EDF parameters [64,67] and constraining the form of some functional terms [33,63,68] based on microscopic calculations to ideas for a full determination of the functional form from a chiral interaction model [22,62,65].
In this work, we follow a hybrid strategy first suggested in Refs.[19,70].It consists in adding terms arising from pion exchanges as described by chiral EFT interactions at the Hartree-Fock (HF) level on top of a Skyrme EDF structure.There are two motivations for this strategy.First, the form of the Skyrme EDF corresponds to calculating HF energies from contact interactions.Following chiral EFT, the first additional degree of freedom that appears when increasing the resolution of the description of the considered systems are the pions exchanged between the nucleons.Adding them explicitly should lead to a more accurate description of nuclear properties.Second, one notices that ab initio calculations with chiral EFT interactions often build correlations on top of an initial mean-field solution.In our approach, we employ the same interactions but instead of generating correlations via the many-body method, we adjust the short-range part of the interactions.This is because the dominant bulk correlations in nuclei, e.g., in expansions around HF, appear to be short range in nature [71] and could therefore be mimicked by contact interactions.
This semiphenomenological strategy was implemented in a series of papers, Refs.[70,[72][73][74][75][76].While improvements over EDFs without chiral terms were observed, the dependency of the results on the order of the chiral interactions showed large variability and puzzling systemat-ics [76].The goal of the present work is to carefully revisit the construction of EDFs incorporating chiral physics via a density-matrix expansion (DME).We study in detail the dependence of the results on the order of the employed chiral interaction and identify which terms are crucial to obtain improvement over EDFs without pionexchange terms.To perform these investigations we construct a new set of nuclear EDFs which we dub "Germany-USA DME EDFs" (GUDE1 for short).
We begin by laying out the theoretical framework of this study in Sec.II.In particular, we discuss the structure of the EDFs including the chiral contributions, the numerical setup used to determine nuclear properties from them, and the parameter optimization protocol.In Sec.III we present the obtained GUDE parametrizations and investigate their performance by comparing against experimental data.In particular, we construct a GUDE variant which reproduces the main improvements found in this work by adding only a minimal number of terms arising from pion exchanges.Section IV contains a detailed analysis of the order-by-order behavior of the functionals in the GUDE family.We end by summarizing our findings in Sec.V, where we also give an outlook on avenues for future work.

II. METHOD
The EDFs we construct in this work can be split into six parts according to (1) They are solved at the Hartree-Fock-Bogoliubov (HFB) level using the code HFBTHO [77], as detailed in Sec.II E. The conventional part of the EDFs consists of the latter four terms.The Skyrme part reads where and the isospin index t = 0 (t = 1) labels isoscalar (isovector) densities.Summations over spacial indices a, b are implied.In Eq. ( 2), we have suppressed the dependence on the position R of the (quasi)local densities, for which expressions can be found in Refs.[20,78].Since we only apply our EDFs for calculations of even-even nuclei, time-odd densities are not taken into account in the construction.Skyrme and pairing E pair (Sec.II D) contributions contain the parameters that are adjusted to data as described in Sec.II F. The Coulomb energy is obtained here as in Refs.[42,76,79,80]: the Hartree term is calculated exactly using the Gaussian substitution method [81,82] and the exchange term is calculated with the Slater approximation [83]; see Ref. [84] for an assessment of the accuracy of these methods.The kinetic energy is given by with ℏ 2 /(2m) = 20.73553MeV fm 2 .
In Sec.III we construct a conventional functional, below labeled as "no chiral", that contains only these four terms and serves as a reference functional for comparing the performance of the other EDFs that we construct following the same optimization protocol.These additionally contain the first two terms in Eq. ( 1), E χ H and E χ F , which represent the Hartree and Fock energy from pion exchanges, respectively.The expressions for the pion exchanges which enter the definitions of E χ H and E χ F are taken directly from interactions derived from chiral EFT at different orders; see Sec.II A. Because the low-energy constants of the pion exchanges are determined from fewbody data [85] and are not adjusted in the present work, the additional inclusion of these terms does not lead to an increase in the number of adjustable functional parameters.See Secs.II B and II C for details regarding the pion Hartree and Fock terms.
While the structure of the functionals constructed here agrees with the one from Ref. [76], we introduce several changes and improve various aspects in the construction and optimization of the functionals compared to that work.These changes, stated in detail in Secs.II A to II F, are mostly driven by the idea to enable a cleaner comparison of the functionals constructed at (different) chiral orders.

A. Chiral interactions
For the construction of the EDFs we consider pion exchanges at different orders in the chiral expansion up to next-to-next-to-leading order (N 2 LO) both with and without the explicit inclusion of intermediate ∆ isobars as well as with and without three-nucleon (3N) forces.Chiral EFT interactions contain pion exchanges and contact interactions.We take only the finite-range parts of the pion exchanges explicitly into account.Expressions for the corresponding interaction terms in coordinate space are given in Refs.[75,76].The low-energy constants that appear are taken from the determination of Ref. [85] (columns "Q 2 , no ∆" and "Q 2 , fit 1" of Table 1 therein).Note that we use g A = 1.27 and h A = 3g A / √ 2 as chosen in Ref. [85].The previous implementation [76] used the Fock coefficient functions derived in Ref. [75] for which the slightly inconsistent combination of g A = 1.29 with low-energy constants from Ref. [85] had been considered.The finite-range interactions are regularized by multiplying them with the local regulator function where we choose R c = 1.0 fm and n = 6 (cf.[86]).Investigating the choice of the regularization scheme is left for future work.Contact interactions as well as correlations involving pions beyond the HF level are assumed to be effectively captured by the EDFs by adjusting the parameters of E Skyrme and E pair to data from finite nuclei.

B. Long-range Hartree terms
The Hartree terms from the pion exchanges are included essentially exactly by evaluating the corresponding integrals.Since we consider only even-even nuclei, the spin density vanishes due to time-reversal symmetry so that only the central part of the NN interactions contribute: To make use of the capability of HFBTHO to solve the HFB equations for potentials given by sums of Gaussians [77], we approximate the central chiral potentials as A similar idea was implemented in Ref. [87].Together with B i = M i = 0 (which do not contribute here due to time-reversal invariance), Eqs. ( 7) and ( 8) correspond to a Gogny-like interaction, Note that in Eqs.(7) to (9) we correct several mistakes compared to Eqs. (30) to (33) of Ref. [76].The wrong equations in Ref. [76] led to an incorrect implementation of the Hartree terms in the functionals constructed therein.
To reproduce the behavior of the regulator [Eq.( 5)] at the origin, the conditions are imposed.The remaining free parameters W i , H i , µ i are obtained by a fitting routine.As in Ref. [76], N = 5 Gaussians are used here as a compromise between accuracy of the approximation and computational requirements for evaluating and storing the resulting integrals [88].The Gaussians used in Ref. [76] were obtained by simultaneously fitting all 13 parameters for the isoscalar V C and isovector W C potentials.Here, we fit first only the nine parameters for the isoscalar potential V C since it contributes significantly more to the energy of finite nuclei than its isovector counterpart.We keep the resulting Gaussian widths µ i fixed for the subsequent fitting of the remaining four parameters of the isovector potential W C .We obtain the parameters of the Gaussians by χ 2 minimizations where the loss functions are given by which are evaluated on an evenly spaced grid from r = 0 to 8 fm with step width 0.125 fm.We include the r 2 prefactor in the definition of the χ 2 to account for the increased importance of larger r due to the presence of the volume element in the Hartree energy, Eq. ( 6).This factor had not been included in the determination of the Gaussian parameters in Refs.[76,88].We provide the Gaussian parameters obtained in the new fit in the Supplemental Material [89].In Figs. 1 and 2 we plot r 2 [V t (r)− Ṽt (r)] including contributions up to including N 2 LO in the chiral expansion (without explicitly resolved ∆ excitations).The new fitting strategy improves the fit of V C without a significant degradation in fitting W C .When evaluating the Hartree energy expectation value in 208 Pb the difference between the value obtained with the exact and the approximated potential at N 2 LO is about 5 MeV (on a total Hartree energy of about 4000 MeV) with the Gaussian parameters obtained in this work.This is a significant improvement over the difference of 37 MeV obtained with the Gaussian parameters of Refs.[76,88].Similar improvements are obtained for the fits of the potentials at other chiral orders.For these comparisons the underlying single-particle orbitals were generated from a selfconsistent HF calculation with the SLy4 EDF [90] using the code HOSPHE [91].
Note that it is not clear if and how the observed improvements translate into improvements of the constructed EDFs.This is because the Skyrme parameters are fitted to data after adding the terms originating in the chiral potentials and this fitting can (partly) compensate the errors from the nonperfect Gaussian approximations.For the same reason it is also hard to gauge a priori the impact of other changes we introduced compared to Ref. [76].
For later reference we introduce a notation for contributions arising when performing a Taylor expansion of one density entering Eq. ( 6) in the relative coordinate r r 2 -weighted difference between isoscalar central potential at N 2 LO in the chiral expansion and its approximations by sums of five Gaussians according to Eq. (7).Both the approximation of Refs.[76,88] and the one obtained here are shown.FIG. 2. Same as Fig. 1 but for the isovector potential.
about the argument of the other density.We write with Finally, we recall that there are no Hartree contributions from the long-range parts of 3N forces at the orders we consider.
C. Long-range Fock terms The Fock energy arising from a local NN interaction V χ is given by A DME allows one to approximately rewrite the nonlocal one-body density matrix ρ as a sum of terms in which the nonlocality is factored out [92].After applying the DME and carrying out the traces and the integral in the nonlocality r, one obtains a quasilocal approximation for the Fock energy, which for the NN forces used here reads [ As before we consider only terms that contribute in timereversal invariant systems.Note J t,aa = 0 when axial symmetry is conserved [93], which is the case for all calculations performed in this work.Equation ( 15) looks very similar to the Skyrme part of the functional, Eq. ( 2).However, in Eq. ( 15) the prefactors of the density bilinears (the g coefficient functions g uv t ) are not constants but functions of the isoscalar density ρ 0 and are fixed once one picks a chiral interaction model and a DME variant.
In the actual HFB calculations with HFBTHO the g coefficients are approximated with interpolation functions of the form where N = 3 and the coefficients guv(w) t (0), a i , b i , c i are fitted separately for each g coefficient.For details on the interpolation see Ref. [76].Note that Eq. (47) therein contains an error which is corrected in Eqs.(17) and (18) above.
In this work we stick to the choice of Refs.[75,76] and use the (simplified) phase-space averaging (PSA) DME [70,74].The DME is applied to the isoscalar and isovector parts of the one-body density matrix using an isoscalar momentum scale, which works well for the former, but not for the latter [92].However, the isovector Fock contributions are small and again we expect the Skyrme parameter fitting to partly compensate the errors.We leave the investigation of the impact of choosing a different DME variant in the EDF construction for future work; see Ref. [92] for a study where similar tests are performed in a non-self-consistent scenario.In that work we found DMEs work well even for pion exchanges at leading order (LO) in the chiral expansion despite the long range of this interaction.Interaction terms at higher orders are of shorter range and therefore expected to be even more suited for a DME treatment.
Note that some of the 3N Fock terms used in Ref. [76] were incorrect; these have been corrected in the present work.We provide the resulting interpolation parameters entering Eqs. ( 17) and ( 18) in the Supplemental Material [89] and introduce the notation for the combination of Skyrme coefficient, Taylorexpanded Hartree contribution, as well as NN and 3N g coefficient functions of the same kind.

D. Pairing contribution
Within the HFB framework, the pairing contribution to our EDFs is given in the mixed-pairing prescription [94] as where ρq (R) are the pairing densities and ρ s = 0.16 fm −3 .The neutron and proton pairing strengths V n 0 and V p 0 are adjusted to data as described in Sec.II F. Because of the zero range of the underlying effective pairing force, a cutoff of E cut = 60 MeV to truncate the quasiparticle space is employed.This cutoff was missing in the implementation of Ref. [76].Thus, in that work the quasiparticle space was truncated implicitly only, via the finite size of the employed basis.
In Ref. [76], we approximated particle number projection with a variant of the Lipkin-Nogami (LN) prescription derived for a seniority-pairing interaction with an adjusted effective strength [95].In Ref. [96] it was shown that this scheme compared well against the numerically expensive variation-after-projection scheme in well-deformed nuclei, but not near closed shells; see also Ref. [97].In addition to the lack of consistency between the actual pairing interaction and the one used for the LN scheme, the LN scheme is not variational.For these reasons, we drop this prescription and work at the HFB level only.Future development of this work's EDFs should involve revisiting particle-number restoration.Note that the UNEDF1-HFB parametrization of the Skyrme EDF was also performed without the seniority-based LN scheme of its parent UNEDF1 and its performance was only slightly worse [98].

E. Hartree-Fock-Bogoliubov calculations
We obtain nuclear ground states based on the EDFs described in the previous subsections by performing HFB calculations.The HFB equations are solved with the DFT code HFBTHO, which expands the single-particle wave functions in a harmonic-oscillator (HO) basis in cylindrical coordinates [84].For calculations of ground states, bases without axial deformation are used.In all cases the basis consists of 20 HO shells and the spherical frequency ω 0 of the HO basis is set according to the empirical formula ω 0 = 1.2 × 41/A 1/3 MeV [82] unless noted otherwise.HFB solutions are obtained iteratively using the kickoff mode of HFBTHO in which an axial quadrupole deformation constraint is applied during at most the first ten HFB iteration steps to guide the solution towards the correct local minimum, then the constraint is lifted [77,82].
F. Optimization of Skyrme and pairing parameters E Skyrme and E pair contain in total 15 parameters C uv t , γ, and V q 0 which need to be determined from fitting to data.Note that E χ H and E χ F are free of adjustable parameters.Thus, the number of optimization parameters is the same for functionals constructed here with and without chiral terms.The volume parameters C ρρ t0 , C ρρ tD , C ρτ t , and γ can be related to properties of infinite nuclear matter (INM).Expressing the exponent γ in terms of INM parameters at saturation gives where quantities indexed "fr" represent the contributions from the finite-range Hartree terms to the INM parameters (see Ref. [99]).P denotes the pressure of symmetric matter at saturation density, C = (3/5)(3π 2 /2) 2/3 , and u c = (3π 2 ρ c /2) 1/3 /m π .The expressions for A γ (u c ) and B γ (u c ) are given in Appendix C of Ref. [73].The equations for the other volume parameters can easily be obtained from the ones given in Ref. [73] by adding the respective contributions from the finite-range Hartree terms [99].
Proceeding in this way we express the volume parameters via saturation density ρ c , saturation energy E sat , incompressibility of symmetric nuclear matter K, isoscalar effective mass M * s , symmetry energy at saturation density a sym , its slope L sym , and isovector effective mass M * v .As in previous works [42,76,79,80] we do not optimize the isovector effective mass but instead keep it fixed at its SLy4 value, M * −1 v = 1.249, which leaves 14 parameters to be optimized.
Using INM properties at saturation density as optimization parameters instead of EDF volume parameters allows us to impose physically motivated constraints on these parameters.The bounds that we impose are not allowed to be violated in our optimization procedure.We take the same bounds as in Refs.[42,76,79,80] except for K and L sym .For the incompressibility K we extend the allowed range to [180,260] MeV based on the analysis of Ref. [100] using different forces from chiral EFT, which obtained a range of [182,262] MeV, and a study assessing the nuclear matter properties of Skyrme EDFs, which used [200,260] MeV based on different experimental and empirical results [101].For the slope parameter L sym we use [30,80] MeV based on the overlapping region of different experimental and theoretical constraints; see Refs.[102,103].Collectively we denote our set of optimization parameters as x.The parameters and their ranges are summarized in Table I.
They are determined by minimizing a loss function, which is given by a weighted sum of squared errors: where s i,j (x) are the EDF predictions and d i,j the data.D T is the number of different data types.In this work we consider ground-state energies of spherical (E sph ) and deformed (E def ) nuclei, neutron (∆ n ) and proton (∆ p ) odd-even staggerings, proton point radii (R p ), and fission isomer excitation energies (E * ), therefore D T = 6.
For every data type i we employ a different inverse weight w i that represents the expected errors in describing the different observables [98].Rather than the somewhat arbitrary values set in Ref. [76], we choose for the weights the estimates determined from the Bayesian calibration of the UNEDF1 functional [104]; see as for UNEDF1.In addition, the form of the functionals (at least for our reference EDF without contributions from chiral EFT) as well as the employed optimization protocol are similar.Figure 3 shows in detail which data types are considered for which nuclei.The experimental data is similar to the data used in Refs.[42,76].However, we exclude single-particle level splittings from the data set.These were introduced in Ref. [42] together with removing the restriction of C JJ 0 = C JJ 1 = 0 for the tensor part of UN-EDF1 in an attempt to improve the description of nuclear shell structure.The reported standard deviations for the tensor coefficients were quite large and the observed improvement of the shell structure relatively small.Because the blocking calculations carried out to determine the single-particle structure are numerically expensive, we therefore decide to remove the single-particle level splittings from the data set.
With those exceptions, we consider the same data types for the same nuclei as in Refs.[42,76].The experimental binding energies -which determine E sph , E def , ∆ n , and ∆ p -are extracted from the 2020 Atomic Mass Evaluation (AME) [105] and the charge radii from Ref. [106]; see Refs.[42,76] for details.For 56 Ni, which had not been measured yet, we take the value determined in Ref. [107].The conversion from charge radius to proton point radius is based on the 2018 CO-DATA recommended value for the proton charge radius r p = 0.8414 fm [108] and the 2022 Particle Data Group average for the neutron charge radius square r2 n = −0.1155fm 2 [109].The fission isomer energies are taken from Ref. [110].The EDF predictions s i,j (x) are obtained for given values of the parameters x at every optimization step by solving HFB equations with the setup explained in Sec.II E. The value of the quadrupole moment used to initialize the kickoff mode is computed by assuming a ground-state deformation of β 2 = 0.3 for deformed nuclei and a fission isomer deformation of β 2 = 0.6 [104].In total, 81 HFB calculations are performed at every optimization step: 77 for the ground states of the nuclei in the data set, for which no axial basis deformation is used, and 4 for the fission isomers, which are calculated with an axial basis deformation parameter of β = 0.4.
We use the predicted average neutron (proton) HFB pairing gap as a proxy for neutron (proton) odd-even staggering.While this is an approximation [111], actually determining odd-even mass differences would require calculating ground states of odd nuclei for which additional EDF terms enter due to broken time-reversal invariance and the determination of odd ground states via blocking calculations is much more involved than calculating ground states of even-even nuclei [112].
To find the parameter set x for which χ 2 (x) is minimized within the bound constraints discussed above we employ the derivative-free optimization algorithm POUNDERS [113,114].It solves the nonlinear least squares problem by constructing a quadratic model for each term in the χ 2 .The resulting quadratic model for the χ 2 is assumed to be valid only within a certain trust region.Minimizing the model in this region yields a solution candidate point.Then new quadratic models are constructed around this point and the trust region is updated.In this way an iterative optimization procedure is obtained; see Ref. [114] for details on the algorithm.POUNDERS needs significantly fewer iteration steps to converge to a minimum than a conventional Nelder-Mead optimization routine [79,113].
At every iteration step, the trust region is essentially a hypersphere around the current candidate point (in a space where the different optimization parameters are scaled as described in Ref. [113]).The hypersphere's radius shrinks when getting closer to the minimum.Sometimes POUNDERS shrinks this radius too quickly despite the current candidate point not being sufficiently close to the optimum yet.In such scenarios, restarting POUNDERS from the current candidate point helps to accelerate the convergence and allows it to possibly jump to another valley in the parameter landscape.Therefore, we restart the optimization every 150 iteration steps and in doing so set the trust region radius back to its initial value of ∆ 0 = 0.1.
We use the parameter sets obtained at different orders in the chiral expansion in Ref. [76] as starting points for the optimization of the corresponding GUDE functionals constructed here.For the reference "no chiral" functional we start the optimization from the UNEDF2 parameters [42].For a few EDFs we carry out the optimizations more than once employing also other Skyrme parametrizations as starting points (e.g., SLy4 [90]).We find that if those optimization runs converge, they converge to the same solutions as the other optimizations.This gives us confidence that the parametrizations we obtain constitute global optima (within the employed bound constraints).

A. GUDE parametrizations
The parameter values obtained from the optimizations described in Sec.II F are given in Table III.Parameters that ended up at their bounds are underlined.We provide the EDF parameters with larger precision in the Supplemental Material [89], both in their explicit representation and equivalently in terms of INM properties.We refer to the Skyrme-type GUDE functional without any chiral terms as "no chiral".The other GUDE EDFs are labeled according to up to which order chiral terms are included and whether they include interaction terms with explicitly resolved intermediate ∆ excitations and 3N forces.We categorize the EDFs according to their properties discussed in the next paragraphs: we refer to the "no chiral" functional as class 0, to the LO and nextto-leading order (NLO) functionals collectively as class 1, and to the remaining functionals as class 2. This latter class contains also a functional labeled "min.chiral".It is constructed with the idea of adding as few terms as possible to the "no chiral" version while still obtaining an EDF that behaves like a member of class 2. Details of the construction of this functional are discussed in Sec.III B. In Table III, the different classes are indicated by vertical lines.
We start with a discussion of the INM parameters of the different GUDE variants.The saturation energy E sat ends up at its upper bound 2 for almost all optimized functionals.This also holds for the value of the incompressibility K for classes 0 and 1.For class 2 the incompressibility acquires lower values inside the allowed parameter range.All other considered nuclear matter parameters also indicate a qualitative difference between classes 0 and 1 on the one hand and class 2 on the other hand: the variation of the INM parameters within these groups is much smaller than the difference between them.The main parameter difference between class 0 and class 1 lies in an increased value of the slope parameter L sym for the chiral functionals.When going to the class 2 functionals, L sym gets significantly reduced and ends up at its lower bound for most of the EDFs, with a correspondingly lower a sym parameter.Note that for some of the EDFs the inverse isoscalar effective mass M * −1 s attains its lower bound, too.While M * −1 s = 0.9 is relatively low compared to typical values [101], this value was also obtained for UNEDF0 [79].
In Fig. 4 we show the energy per particle for pure neutron matter and symmetric nuclear matter for four functionals constructed in this work; one each from class 0 and 1 and two from class 2. The differences between the EDFs are very small up to about saturation density.This is not surprising since this region is probed by finite nuclei and hence strongly constrained by the fit to experimental data.The differences between the different classes become much more pronounced for ρ 0 ≳ ρ c , in particular for neutron matter.This region is not probed by finite nuclei, which is also why the deviation from the additionally given ab initio result observed for class 2 in this density regime for neutron matter is not surprising.The plotted uncertainty bands have been obtained by Drischler et al. [102] based on the MBPT calculations from Ref. [115] with a chiral NN+3N Hamiltonian at N 3 LO with momentum cutoff 500 MeV [116]  sults.Note that the curves for the two class 2 representatives, the N 2 LO∆+3N and the "min.chiral" variant, are very close to each other even for ρ 0 > ρ c .This holds analogously for other EDFs from the same class.
Overall, and in particular within the classes as defined above, the description of INM at saturation density and below shows a large consistency between the different functionals.This may be considered surprising given that the chiral contributions are quite different in size depending on the chiral order.However, it indicates that the optimization of the Skyrme and pairing coefficients to data can, to a large degree, wash out the effect of the additional terms.We return to this issue in Sec.IV.
In Table III we also provide the value of the γ exponent for the different EDFs.Compared to the "no chiral" variant it is larger for class 1, but smaller for class 2, indicating that the density-dependent terms absorb different physics for the two classes.Along the same lines we note that at every order γ is smaller by about 0.05 for functionals including chiral 3N contributions.
For all GUDE variants the generally observed hierarchy of pairing strengths 97,117] holds.The somewhat weaker strengths obtained for the class 2 EDFs when comparing to the other classes is in agreement with the lower M * −1 s values for class 2 [20].Note that a direct comparison of the surface parameters of the different GUDE variants makes little sense because the chiral contributions to the corresponding terms depend on the functional and are not included in the C uv t values given in Table III.
Based on starting optimization runs of the same GUDE variant from different initial points [118] we find that the parameters of the isovector part of the EDF are relatively ill-constrained with our optimization protocol.This is in agreement with observations made in other nuclear EDF optimizations [42,43,79,99,119].To better determine the isovector parameters the optimization data set has to be augmented; see also Sec.V. Also the C JJ 0 parameter seems poorly constrained.To quantify these statements a rigorous statistical analysis should be carried out in future work.
The last row of Table III contains the value of the loss function χ 2 at the optimum.For the "no chiral" EDF it is around 120. Adding the chiral terms at LO (and NLO) according to the construction described in Sec.II worsens the χ 2 at the minimum: it attains values around 145.This stems from a slightly worse description of groundstate and fission isomer energies.
However, the additional inclusion of chiral terms at N 2 LO or of the ∆ contributions at NLO reduces the χ 2 at the minimum to about 90.In particular experimental energies of spherical nuclei in the fitting set are better described by the class 2 functionals.The root-mean-square deviation (RMSD) for those is 2.5 MeV for the "no chiral" EDF, but only 1.6 MeV for the class 2 GUDE variants.The other data types in the χ 2 are typically either slightly improved or are equally well described when comparing to the "no chiral" functional.TABLE IV.Exact scalar Hartree energies and differences of scalar Hartree energies calculated with Taylor expansions of the densities up to a given order [cf.Eq. ( 12)] and the corresponding exact energies (all in MeV).The densities are generated from calculations with the SLy4 EDF.Results are given for the chiral pion exchanges considered here at N 2 LO and for the finite-range parts of the Gogny D1S functional [121].We note that the N 2 LO EDF constitutes a slight deviation to these general trends (which can also been seen from some of the parameter values listed in Table III): it describes the radii in the χ 2 worse than all other EDFs but proton odd-even staggerings are much improved.

B. Investigation of GUDE class 2 and construction of "min. chiral" functional
As discussed in the previous section and further in Sec.III C we observe an improvement over the "no chiral" functional when going to EDFs that include chiral terms entering at N 2 LO (or NLO when including interactions with explicit ∆ excitations).It turns out that only a small subset of the terms that contribute at these orders is actually necessary to achieve the improvement.
First, the inclusion of chiral isovector contributions is not required.This is hardly surprising given that the Skyrme part of the EDFs contains six parameters contributing solely to the isovector part, which is to be compared to seven parameters for the isoscalar terms, but the isoscalar energy contributions are at least an order of magnitude larger than the isovector ones [120].The similar number of parameters for the two EDF parts suggests one may expect a similar relative precision for the corresponding energy contributions.The resulting absolute deviations would then be much bigger for isoscalar energies.Thus, one can expect omitting chiral isovector contributions does not significantly impact the description of bulk properties of finite nuclei (after refitting the EDF parameters).Of course this is amplified by the inadequacy of the optimization data set to accurately fix the EDF isovector parameters.
Performing an optimization of an EDF as described by Eq. ( 1) but taking into account from the chiral side only Fock contributions up to N 2 LO yields a class-1-like functional which suggests that the switch to class 2 is due to the Hartree terms.Indeed N 2 LO (NLO with ∆s) is the first order which for even-even systems has isoscalar pion-exchange Hartree contributions.These are by far the largest chiral contributions to the energy.In Table IV we show the expectation values of the exact Hartree energy from pion-exchange contributions up to N 2 LO in the chiral expansion.They are obtained with densities generated from calculations performed with the code HOSPHE [91] employing the SLy4 EDF [90].Additionally, we provide the difference to these exact values for energies that we obtain when Taylor expanding one density entering the Hartree energy; see Eq. (12).For comparison we also provide the analogous numbers obtained with the finite-range parts of the Gogny EDF in the D1S parametrization; see Ref. [120] for a more extensive study.
One can see that the energies obtained with the Taylor series converge relatively slowly towards the exact values.In particular when going to second order in the Taylor expansion the approximated value is still off by about 40 MeV in 208 Pb.The second-order expression for the energy has a Skyrme-like structure (density bilinears consisting of up to second-order densities multiplied with constant prefactors).Therefore, one may expect that a Skyrme EDF cannot fully account for the chiral Hartree contributions at N 2 LO if they are left out (as is the case for classes 0 and 1).It is thus conceivable that class-2 GUDE variants behave differently from classes 0 and 1. 3  Carrying out the optimization of an EDF where in the chiral part only the isoscalar Hartree contributions entering at N 2 LO are included leads to a functional with χ 2 ≈ 112 at the minimum, which is clearly larger than the values observed for class 2. For this EDF the pairing strengths take a nonphysical value V q 0 ≈ 40 MeV fm 3 .These observations suggest another term is additionally needed to reproduce the class-2 behavior.
In Fig. 5 we show the contributions to the g ρρ 0 coefficient arising at different chiral orders, but the following discussion applies also similarly to other g coefficients.The total g ρρ 0 coefficient at a given order is the sum of all depicted contributions ∆g ρρ 0 up to that order.The LO contribution shows a strong density dependence with its value at ρ 0 = 0 being about five times as large as the value at ρ 0 = ρ c .The contributions at NLO and N 2 LO are much smaller and their density dependence is much weaker, which is why their effects can be easily captured by simply adjusting Skyrme coefficients.In principle even the strongly density-dependent LO coefficient could be quite well mimicked by a Skyrme EDF due the presence of the C ρρ tD ρ γ 0 term, but, since this term has to capture several different types of unresolved physics [122], one may expect that adding the LO g ρρ 0 contribution explicitly still has a relevant effect.Optimizing an EDF with both isoscalar chiral long-range Hartree contributions at N 2 LO and Fock contributions at LO yields a functional belonging to class 2 as desired. 3Note that the argument put forward above is not a direct proof because the fitting of the EDF parameters may shuffle around contributions among more terms than the ones technically entering the Taylor-expanded energy.
FIG. 5. Contributions to g ρρ 0 arising at different chiral orders.We show contributions at LO, NLO, and N 2 LO calculated from the interaction specified in Sec.II A using the PSA-DME.In addition, we the LO contribution when using the Slater approximation instead of the PSA-DME.
We showed in previous work [92] that Fock energies from a Yukawa interaction can be well approximated by using the Slater approximation instead of the more involved PSA-DME applied in this work so far.However, this comes at the price of a worse local reproduction of the Yukawa Fock energy density essentially everywhere in the nucleus.Using the Slater approximation instead of the PSA-DME reduces the amount of nonvanishing isoscalar NN g coefficient functions from five to one.We show the nonzero g ρρ 0 coefficient at LO in Fig. 5.We find that the resulting EDFs differ by similar amounts as other functionals in class 2 differ from each other.Therefore it seems safe to use the simpler Slater approximation in the present EDF construction, at least for bulk properties.
We refer to the EDF constructed according to Eq. ( 1) including for E χ H only the isoscalar NN pion-exchange Hartree contribution entering at N 2 LO and as E χ F the isoscalar NN pion-exchange Fock contribution at LO (described by the Slater approximation) as the "min.chiral" GUDE variant.The parameters obtained when optimizing this functional are given in Table III and with higher precision in the Supplemental Material [89], where we also provide the parameters used in the interpolations for the chiral Hartree and Fock contributions according to Eqs. ( 7) and (17).The INM parameters and the χ 2 value at the optimum are in the ranges of the other class-2 functionals (see Table III), indicating that the "min.chiral" variant indeed also belongs to this class.This explicitly demonstrates that the two identified terms are enough to achieve the improvement over classes 0 and 1.

C. Global comparison to experiment
We now investigate the performance of the different functional variants in the GUDE family obtained in Sec.III A by calculating the ground states of eveneven nuclei included in the 2020 AME [105].We include all 663 nuclei with actual measured masses, leaving out those for which only evaluated masses are available.Every nucleus is calculated five times with HFBTHO in kickoff mode setting the initial deformation constraint to β = −0.2,−0.1, 0, 0.1, 0.2.This is done so that oblate deformed, spherical, and prolate deformed solutions are considered as possible ground states for every nucleus.
The HFB calculations are carried out until they are converged (typically within at most about 100 HFB iteration steps) or until the amount of unconverged calculations for a given functional does not get further reduced for at least 800 HFB steps.For most GUDE variants only about a handful of the 3315 calculations end up unconverged at the end of this procedure.The N 2 LO EDF is the only exception from this rule: even after more than 3000 HFB steps, 111 calculations are still unconverged.Note, however, that only four of those constitute the calculation with lowest binding energy for the corresponding nucleus.
For every nucleus, we pick among the converged calculations the one with the lowest energy as a first groundstate candidate and apply on it two filters to exclude unphysical solutions.Whenever a filter is triggered, the calculation with the next-lowest energy for the same nucleus is considered instead.First we do not consider solutions with E/A < −11 MeV.This filter turns out to be triggered only a few times by calculations with EDFs that include interactions with explicit ∆ isobars in the chiral terms.Second we apply a filter to remove solutions with unphysically large deformations.This is done by applying the 1.5 interquartile range rule, which is a simple measure to detect outliers of a distribution, on the values of the deformation parameter β 2 of all remaining ground state candidates.The β 2 parameter is much less mass-number dependent than the axial quadrupole moment of the nucleus Q 20 and is related to it according to with the root-mean-square matter point radius R m .The deformation filter is in practice triggered at most for two nuclei per EDF.We compare the resulting ground-state energies against the values extracted from the 2020 AME.Table V contains the corresponding root-mean-square and mean deviations obtained for nuclei with Z ⩾ 8.We also give the deviations of the two-neutron (S 2n ) and two-proton (S 2p ) separation energies obtained from the same data set, and of the charge radii from Ref. [106].GUDE variants of the same class behave very similar for all these quantities with the only exception being somewhat larger mean deviations observed for separation energies for the N 2 LO functional compared to other class-2 EDFs.6. Distributions of ground-state energy differences between calculated and experimental results.They are shown for the "no chiral" and "min.chiral" GUDE functionals in bins with a width of 1 MeV each.Note that the last bin contains also values with an energy difference larger than 10.5 MeV.While classes 0 and 1 perform similarly, an improvement is observed for all observables when going to class 2. In particular, the ground-state energy RMSD is significantly reduced by roughly 30% from 2.1 MeV for classes 0 and 1 to about 1.5 MeV for the various class-2 EDFs.The mean deviation ⟨E theo −E exp ⟩ is almost halved down to 0.3 MeV, indicating that the energies are less biased towards underbinding for class 2. This can also be seen in Fig. 6,4 which shows the histogram of the quantity E theo − E exp .Calculations which produce extremely underbound nuclei (those at the very right of the distribution) occur much less often for the class-2 "min.chiral" functional than for the reference "no chiral" EDF.Such cases correspond mostly to very light nuclei.For the class-2 variants almost half of all nuclei are described with a mass error of less than 0.5 MeV.Note that while the binding energies included in the χ 2 are described better by class 0 than by class 1, the performance on all eveneven nuclei binding energies is very similar for these two classes.
In the upper row of Fig. 7 we show ground-state energy residuals for four GUDE variants.One can clearly see that the class-2 EDFs describe energies around the N = 82 and N = 126 shell closures much better than the class-0 and -1 variants.We note that, due to the parameter optimization involved in the construction of every functional, it is not clear if the additional chiral terms entering the class-2 functionals are actually directly improving the description of (near-)closed-shell nuclei or if they instead improve the open shells and indirectly allow the parameter optimization to yield a better reproduction of closed shells.In addition, the observed underbinding for light nuclei is reduced for the class-2 variants.
For both two-neutron and two-proton separation energies, class-2 EDFs give a small improvement over classes 0 and 1: the RMSD values are reduced by about 12%.In addition, the bias on S 2n values is almost completely gone while it is increased for S 2p .
The description of charge radii is least affected by the additional chiral terms added in class 2. This can also be seen in the lower row of panels of Fig. 7. Charge radii are only slightly better described for N ≈ 40 to 100 and their mean deviation is slightly closer to zero for class 2.

D. Shell structure and deformation properties
To investigate the quality of the GUDE family with respect to nuclear shell structure, we compute singleparticle levels using blocking calculations; see Refs.[80,112] for details on the procedure.Using blocking calculations at the HFB level is both logically consistent with the construction of the functionals at the HFB level and helps with reducing systematic errors when comparing with experiment [80].Calculations use the same setting for the HO basis as before, namely with 20 full, spherical shells.In this context one should be reminded that single-particle energies are not observables but extracted in a model-dependent way from experiment [123,124].
Here we compare to the values given in Ref. [125].Furthermore, it is well known that the single-particle shell structure depends strongly on beyond-mean-field effects such as particle-vibration couplings [126][127][128][129].As a consequence, blocking calculations should not be expected to perfectly match "experimental" single-particle data in closed shell nuclei.They are simply meant as a validation check to guarantee that basic features of the nuclear shell structure are properly reproduced.
As an illustrative example, we show in Fig. 8 the obtained neutron single-particle spectra of 208 Pb for selected GUDE EDFs representative of the different classes.One can make the following general observations.First, the single-particle levels turn out to be largely insensitive to the GUDE variant.Second, the obtained shell gaps in 208 Pb are in good agreement with the ones extracted from experiment and a little better reproduced than for the UNEDF1 functional.Third, the level ordering of the occupied neutron orbitals is also in slightly better agreement with experiment.These qualitative conclusions apply to other doubly closed shell nuclei and suggest a decent reproduction of the shell structure by the GUDE functionals.
Next, we test deformation properties of the EDFs on the standard fission benchmark case of 240 Pu.The HFB calculations are carried out in a deformed HO basis with 30 shells included and with the HO frequency and basis deformation optimized for that nucleus; see Ref. [130] for details.A constraint on the octupole moment is imposed during the first ten iterations to ensure the fission goes through the most likely pathway.Calculations assume axial symmetry.
In Fig. 9 we show the deformation energy, i.e., the energy difference between the configuration with given deformation and the ground state, as a function of the quadrupole moment for selected GUDE functionals as well as for UNEDF1 for comparison.Since including triaxiality typically reduces the height of the first fission barrier by about 2 MeV [81,130], the overall agreement with values extracted from experiment [131] is in fact very good for all considered GUDE variants.The energy of the fission isomer E * is predicted too low by about 1 MeV compared to the value used in the optimization set (2.8 MeV) [110].Seeing that the results for UN-EDF1, UNEDF2, and the DME EDFs of Ref. [76] agree very well with this experimental value, this is probably a consequence of the reduced weight of fission isomer energies in the present optimization protocol.Note that a newer experimental estimate for the fission isomer energy of 2.25 MeV [132] is closer to the GUDE values.
For values of Q 20 larger than the value at the fission isomer state a clear difference between results obtained for classes 0 and 1 and class 2 emerge as already observed for other quantities in this paper.We may speculate that such differences are the result of a competition between bulk and shell effects.Table III and Fig. 10 show that the symmetry energy a sym and the surface coupling function W surf (defined below), respectively, differ substantially for the class-0 and -1 and the class-2 parametrizations.For classes 0 and 1, the value of the symmetry energy is a sym ≈ 30 MeV while it is a sym ≈ 28.5 MeV for class-2 EDFs.The surface coupling function, which contains the full contribution to the isoscalar surface energy (Skyrme plus chiral terms), is given by where arises from integrating by parts: W surf is for intermediate densities much stronger for class-2 functionals than for classes 0 and 1. Together, a sym and W surf impact the surface and surface-symmetry contributions to the bulk energy, which are known to be key drivers of deformation properties [133,134].At the same time, Fig. 8 also shows a small but visible difference in the neutron shell structure between class 2 and the other GUDE variants functionals, with the N = 126 shell gap being a little smaller for class 2. Such differences will be amplified as deformation increases and this could play a role in the deformation energy.

IV. ANALYSIS OF CHIRAL CONTRIBUTIONS
In this section we analyze why the only significant effects we obtain from including chiral interactions explicitly into the GUDE functionals occur for the switch from class 1 to class 2, i.e., at N 2 LO (NLO when including ∆ isobars explicitly) in the chiral expansion.
As stated in Sec.III, only little change over the reference "no chiral" EDF is seen when going to LO in the present construction; see especially Table V.This is not surprising since one-pion exchange is known to largely average out for bulk properties [73,135] because at this order pions enter at the mean-field level only through Fock contributions, which are small.For nonbulk quantities such as behaviors along isotopic chains, small differences between the "no chiral" and LO EDFs are visible; see for instance the oxygen chain shown in Fig. 11.
At NLO pions enter at the HF level only through Fock and isovector Hartree contributions.Since these are very small and can be captured well by Skyrme terms due to the weak density dependence of the resulting g coefficients (see, e.g., Fig. 5), the almost identical performance 18 20  of the LO and NLO functionals is to be expected.When going to N 2 LO a significant improvement is achieved, in particular for the global description of ground-state energies.The detailed analysis of Sec.III B indicates that the interplay of two contributions is responsible for this.The attractive pion Hartree contribution at N 2 LO is large and apparently cannot be completely mimicked by Skyrme terms only.Its addition together with LO Fock terms leads to the improvement.
While the incompressibility is at its upper bound for classes 0 and 1, it is much smaller for the N 2 LO EDF (and the other class-2 ones); see Table III.This is probably a consequence of the strongly attractive central isoscalar two-pion exchange entering at N 2 LO in the chiral expansion [136].
This observation raises the question whether the additional chiral terms in class 2 lead to a better description of experiment by themselves or whether the improvement is realized indirectly by moving the unbounded optimum "closer" to the bound constraint region and thereby reducing the achievable χ 2 values within this region.To address this issue one could perform an unconstrained optimization for the different GUDE functionals.Preliminary unconstrained optimizations suggest that the latter mechanism is the dominant one because the difference of the obtained χ 2 values largely seems to vanish for the unbounded optima.Note, however, that these conclusions are preliminary, since for some of the EDFs competing minima seem to occur during the unbounded optimizations and sometimes the unconstrained optima seem to correspond to situations where some INM parameters attain values far away from physically expected regions (e.g., L sym ≈ 5 MeV).We leave the resolution of these issues for future work.
Similar improvement as for the N 2 LO EDF is observed for the NLO∆ EDF.This reflects the fact that in ∆-full chiral EFT the dominant two-pion-exchange contribution is promoted from N 2 LO to NLO [85].At N 2 LO∆ some additional attraction is brought in.For the interactions used here the additional contributions (which in ∆-less chiral EFT would occur in part at even higher orders) are similar in size as the difference between the chiral contributions at N 2 LO and NLO∆.The GUDE functionals are generally not sensitive to such differences on a qualitative level; see Table V.
All statements made above dealt with chiral NN interactions only.The inclusion of 3N forces does not seem to have a significant effect on the description of nuclei and INM at any considered order; see Sec.III.In ab initio calculations, 3N forces are important for a quantitative reproduction of nuclei, and are key for shell structure and for the limits of bound nuclei [137,138].For instance, for the oxygen isotopes, the additional repulsion from 3N forces moves the location of the predicted neutron drip line in agreement with experiment [7,[139][140][141][142].In Fig. 11 we show the ground-state energies of oxygen isotopes as predicted by a few GUDE functionals.Comparing the N 2 LO∆+3N results with the other EDFs shows that including 3N forces does not move the location of the neutron drip line for the EDFs.Similar conclusions hold for the other GUDE variants with 3N forces.In agreement with other EDF calculations [54], all EDFs constructed in this work predict 28 O to be the heaviest oxygen isotope stable against emitting two neutrons, while experimentally it is 24 O.
The crucial difference between the negligible role of 3N forces observed here and their relevant effects in ab initio calculations lies in the fact that the EDFs constructed here yield good saturation properties also without the presence of chiral 3N forces -see Table III and Fig. 4 -while they are absolutely necessary to achieve reasonable saturation in calculations of INM employing chiral interactions [143,144].In such ab initio calculations, the role of 3N forces is already visible at the HF level, so one could have expected an impact also here.The fact that this is not the case suggests the fitted EDF terms can compensate missing 3N pion exchanges in the density regime relevant for finite nuclei.
For the terms which depend only on ρ 0 this is illustrated in Fig. 12, which shows W ρρ 0 for different GUDE functionals.The curves for N 2 LO with and without 3N forces are basically on top of each other, signaling that for the EDF without 3N pion exchanges the Skyrme part of the EDF mostly takes over the role of the 3N terms (see also the different γ values in Table III).This observation correlates well with the original reason to introduce a density-dependent coefficient into nuclear EDFs, namely to replace a genuine 3N interaction [145].
The observation that fitting the EDF parameters can almost fully compensate missing 3N pion exchanges is in apparent contradiction with the wrong drip line position observed for the oxygen chain.In other words the question is, why does the GUDE family predict the wrong drip line location even though the functionals either explicitly contain or are essentially able to effectively encapsulate chiral 3N physics?One simple explanation is the lack of sufficiently neutron-rich nuclei in the experimental data set used in the optimization.Since chiral 3N contributions grow with increasing neutron number [139,146], the description of nuclei closer to stability might not be significantly altered but drip lines might be much improved when optimizing an EDF with chiral 3N contributions using an experimental data set containing more asymmetric nuclei.Another reason is the importance of beyond-mean-field effects that are known to significantly impact the nuclear structure in light nuclei [147,148].As alluded to above, the existence of strict bounds that we impose on some EDF parameters during their optimization somewhat complicates the analysis of the effect of different chiral contributions.Some conclusions drawn in the present section might thus not hold in other optimization settings.

V. CONCLUSIONS AND OUTLOOK
In this paper, we constructed semi-phenomenological EDFs, dubbed GUDE, consisting of pion exchanges taken from chiral EFT at different orders and a phenomenological Skyrme part.The long-range pion-exchange interactions are included at the Hartree-Fock level (using a DME for the Fock contributions) without adjustment and thereby do not change the number of free EDF parameters.The GUDE functionals with chiral terms perform significantly better than a reference Skryme functional without chiral terms constructed within the same protocol, especially in terms of accurately describing groundstate energies.These improvements can be traced back to the combination of two terms: Fock contributions from one-pion exchange at leading order in the chiral expansion and Hartree contributions from two-pion exchange at N 2 LO.This is demonstrated with the "min.chiral" variant of the GUDE EDFs which contains only those two terms in addition to the phenomenological part and achieves similar improvements as observed for the other class-2 GUDE functionals, which contain additional terms stemming from pion exchanges.
Conversely, adding only pion-exchange terms at LO or NLO does not give any improvement.While it might seem like a contradiction to the chiral EFT power counting -according to which the importance of additional terms is reduced with every higher order included -it may simply result from the fact that we include pion exchanges only at the HF level, i.e., beyond-mean-field effects from pions are not explicitly included and the structure of the contact interactions present in the EDFs does not change with increasing order, unlike in chiral EFT.Along similar lines, including long-range 3N forces does not yield significant improvement because the optimization procedure of the density-dependent contact terms in the traditional part of our EDFs allows for the approximate capture of their effects.
The order-by-order systematics of the GUDE functionals shows much less variability and surprising behavior compared to what was observed in Ref. [76], where functionals had been constructed following the same strategy as used here.In particular, we consider it promising that the inclusion of chiral long-range 3N forces does not lead to a worsening of the EDFs, unlike before.We attribute this to the different improvements, bugfixes, and other changes established in the present work.The analysis carried out in Sec.IV mostly explains the obtained order-by-order behavior.In some regards further insight is still needed.For instance, the detailed mechanism how the improvement is realized at N 2 LO (and why some LO terms are additionally needed which on their own do not provide improvement) is still unclear.We believe that insight might be gained from performing optimizations without imposing bound constraints on INM properties.It would also be of interest to investigate if adding pionexchange terms, in particular those included in the "min.chiral" variant, to other functionals, of Skyrme or other type, gives similar improvement as observed here.We have also left the study of the dependence of the EDFs on the chiral interactions including their regulators for future work.
Going beyond NLO in the present construction does not only improve the description of finite nuclei, it also considerably changes properties of INM as shown in Table III.The incompressibility K is significantly reduced and isovector parameters also change strongly.The decrease of the slope parameter L sym is particularly strong, with it typically ending up at our optimization protocol's lower bound of 30 MeV.
However, in current EDFs isovector terms are generally poorly constrained [43,149]; the present work is no exception.This is not of significant consequence when comparing to bulk properties of experimentally accessible nuclei as done here, but limits the predictive power for applications to extreme neutron-rich conditions in astrophysics.This is because the size of isovector contributions grows significantly when going to very neutron-rich systems.Including experimental data on neutron skins or dipole polarizabilities [149][150][151] in the optimization the EDF parameters, possibly combined with fitting to ab initio results for neutron drops [152][153][154][155], is expected to reduce the uncertainties on the isovector terms.
Extending the optimization data set could also be beneficial in other ways.Examples are the inclusion of ground-state information for nuclei close to the neutron drip line to better constrain isovector terms and to study the effect of chiral 3N forces, and the explicit inclusion of separation energies, which could help with their description and would therefore have significant impact on nucleosynthesis yields from r-process calculations [37,38,40].All GUDE variants underbind nuclei on average.This might be remedied by increasing the amount of data from open-shell nuclei in the fit or by adjusting the data weights in the optimization.
For practical applications, correlated uncertainties (or better, distributions) for the EDF parameters should be determined.They could be estimated using Bayesian inference; see Refs.[43,104,156] for example applications to EDFs.Such a scheme could also be extended to incorporate expectations for INM parameters via prior distributions in the optimization instead of imposing them as hard parameter bounds as done here.
The GUDE family may be plagued be self-interaction issues [157].For the chiral contributions this is because Fock contributions are included via a DME but the Hartree contributions are included quasiexactly by approximating the chiral potentials as sums of Gaussians.However, this could be remedied by also treating the Fock terms (at the same chiral order) quasiexactly, which does not lead to significant computational overhead.In this work, we used the DME because this simplifies the inclusion of 3N forces in EDF frameworks.However, their inclusion did not lead to significant improvement and they could thus be left out (like in the N 2 LO GUDE version).Treating self-pairing effects [157], that also occur for conventional functional parametrizations, would require larger adjustments of the EDF structure.
Our work shows that the explicit inclusion of longrange pion-exchange interactions from chiral EFT at the HF level into a Skyrme EDF improves the description of finite nuclei.This suggests that such terms will be relevant when generating an EDF completely from first principles.It might be necessary to account for effects of different types of correlations explicitly to create such an EDF.Collective correlations may be expected to be captured by going beyond the mean-field description.However, different instabilities and pathologies occur when EDFs not derived from actual Hamiltonians are used in those frameworks [45].Therefore, functionals of the GUDE form could not be used directly.Partially, these issues would be addressed by incorporating the pion exchanges consistently quasiexactly as discussed above.Including effects from short-distance correlations from resummed ladder diagrams as described by Brueckner-Hartree-Fock theory should be simpler: in Ref. [71] density-dependent Skyrme terms generated from a counterterm expansion capturing such correlations were computed.A next step towards ab initio EDFs could therefore be the inclusion of such terms.

FIG. 1 .
FIG. 1.r 2 -weighted difference between isoscalar central potential at N 2 LO in the chiral expansion and its approximations by sums of five Gaussians according to Eq. (7).Both the approximation of Refs.[76,88] and the one obtained here are shown.

FIG. 4 .
FIG. 4. Energy per particle in infinite nuclear matter for selected GUDE functionals constructed in this work.For each EDF, both pure neutron matter and symmetric nuclear matter energies are shown.The bound constraints on saturation density, saturation energy, and symmetry energy employed in the optimization of the EDFs are also depicted.For comparison, we show the 1σ uncertainty bands from a calculation employing a chiral Hamiltonian by Drischler et al. for ρ0 ⩾ 0.05 fm −3 [102].
FIG.6.Distributions of ground-state energy differences between calculated and experimental results.They are shown for the "no chiral" and "min.chiral" GUDE functionals in bins with a width of 1 MeV each.Note that the last bin contains also values with an energy difference larger than 10.5 MeV.

FIG. 7 .
FIG.7.Differences of ground-state energies (upper panels) and charge radii (lower panels) for even-even nuclei between values obtained with selected GUDE variants and experiment.See text for details on the experimental data.

FIG. 9 .
FIG. 9. Deformation energy of 240 Pu as a function of the axial quadrupole moment.Calculations assume axial symmetry.
)[J t,aa J t,bb + J t,ab J t,ba ] .

TABLE I .
Parameters optimized in this work and their bound constraints.

TABLE II .
Characteristics of the components of the loss function.ni is the number of data points for each data type i and wi is the inverse weight.For the latter, all units are MeV except Rp which is in fm.

TABLE III .
(22)meters of the different GUDE variants obtained in this work.Values that are underlined correspond to cases where the minimum was attained at a parameter bound.ρc is given in fm −3 , Esat, K, asym, and Lsym are in MeV, the surface coefficients C ρ∆ρ are in MeV fm5, and the pairing strengths V q 0 are in MeV fm3.The last row gives the value of the loss function(22)at the minimum.

TABLE V .
Deviations of ground-state energies, two-neutron and two-proton separation energies (all in MeV), and charge radii (in fm) calculated with the different GUDE variants and the corresponding experimental values.The upper half of the table contains root-mean-square deviations, the lower half lists mean deviations.The values are calculated from all even-even nuclei with Z ⩾ 8 included in the experimental data sets, see text for details on those.