Microscopically-based energy density functionals for nuclei using the density matrix expansion: Full optimization and validation
R. Navarro Perez N. Schunck A. Dyhdalo R. J. Furnstahl S. K. Bogner
nucl-th hep-ph nucl-ex







APS/123-QED















[]navarrop@ohio.edu



Institute of Nuclear and Particle Physics and Department
  of Physics and Astronomy, Ohio University, Athens, OH 45701, USA


[]schunk1@llnl.gov



Nuclear and Chemical Sciences Division, Lawrence Livermore
  National Laboratory, Livermore, CA 94551, USA


[]dyhdalo.2@osu.edu



Department of Physics, The Ohio State University,
  Columbus, OH 43210, USA


[]furnstahl.1@osu.edu



Department of Physics, The Ohio State University,
  Columbus, OH 43210, USA


[]bogner@nscl.msu.edu



National Superconducting Cyclotron Laboratory and
  Department of Physics and Astronomy, Michigan State University, East
  Lansing, MI 48824, USA







  Background
Energy density functional methods provide a generic framework to compute 
properties of atomic nuclei starting from models of nuclear potentials and the 
rules of quantum mechanics. Until now, the overwhelming majority of functionals 
have been constructed either from empirical nuclear potentials such as the 
Skyrme or Gogny forces, or from systematic gradient-like expansions in the 
spirit of the density functional theory for atoms.

  Purpose
We seek to obtain a usable form of the nuclear energy density 
functional that is rooted in the modern theory of nuclear forces. We thus 
consider a functional obtained from the density matrix expansion of local nuclear 
potentials from chiral effective field theory. We propose a parametrization of 
this functional carefully calibrated and validated on selected ground-state properties that is suitable 
for large-scale calculations of nuclear properties.

  Methods
Our energy functional comprises two main components. The first component is a non-local
functional of the density and corresponds to the direct part (Hartree 
term) of the expectation value of local chiral potentials on a Slater 
determinant. Contributions to the mean field and the energy of this term are 
computed by expanding the spatial, finite-range components of the chiral 
potential onto Gaussian functions. The second component is a local functional of the 
density and is obtained by applying the density matrix expansion to the 
exchange part (Fock term) of the expectation value of the local chiral 
potential. We apply the  optimization protocol to determine the 
coupling constants of this energy functional.

  Results
We obtain a set of microscopically-constrained functionals for local chiral 
potentials from leading-order up to next-to-next-to-leading order with and 
without three-body forces and contributions from Δ excitations. These 
functionals are validated on the calculation of nuclear and neutron matter, 
nuclear mass tables, single-particle shell structure in closed-shell nuclei and 
the fission barrier of ^240Pu. Quantitatively, they perform noticeable 
better than the more phenomenological Skyrme functionals.

  Conclusions
The inclusion of higher-order terms in the chiral perturbation expansion seems 
to produce a systematic improvement in predicting nuclear binding energies 
while the impact on other observables is not really significant. This result is 
especially promising since all the fits have been performed at the 
single-reference level of the energy density functional approach, where 
important collective correlations such as center-of-mass correction, rotational 
correction or zero-point vibrational energies have not been taken into account 
yet.









Microscopically-based energy density functionals for nuclei using the density matrix expansion: Full optimization and validation
    S.K. Bogner
    December 30, 2023
================================================================================================================================








§ INTRODUCTION


chiEFT provides the framework for the modern theory of nuclear forces 
<cit.>. It allows the systematic construction of 
nuclear interaction potentials from first principles by introducing an expansion 
of the momentum and pion mass over the chiral symmetry breaking scale (of the order of 1 
GeV). Using chiral interactions to compute properties of heavy nuclei relevant 
to applications such as fission, nucleosynthesis or superheavy science poses a 
number of challenges. These interactions are presumed to 
represent realistic in-medium nuclear forces. Therefore, they should only be 
used in the framework of many-body methods that fully incorporate all many-body 
correlations induced by these potentials. In light nuclei, the no-core shell 
model <cit.> or Quantum Monte-Carlo methods <cit.> are 
popular examples of such direct approaches; in heavier nuclei, alternative 
methods such as the coupled-cluster <cit.> or in-medium similarity 
renormalization group <cit.> can provide good approximations of the 
exact many-body solution for nuclei near closed-shell. In spite of 
very impressive recent success, the majority of nuclei remain out of reach to 
ab initio methods, and the most microscopic approach available relies on 
the nuclear EDF formalism <cit.>.

The EDF approach stands in contrast to ab initio approaches in that 
it is based on enforcing that the wave-function of the nucleus take a simple 
form such as a Slater determinant in the HF theory or a quasiparticle 
vacuum in the HFB theory[This statement applies in the 
SR-EDF version of the EDF approach. In its MR-EDF version, 
one needs to consider two different reference states.]. By definition, such an ansatz for the many-body 
wavefunction cannot be compatible with the use of realistic potentials, and 
most energy functionals are instead derived from phenomenological NN potentials 
such as the zero-range, Skyrme potential or the  
finite-range Gogny force <cit.>. The parameters of these 
effective NN forces are typically adjusted on properties of nuclear 
matter or finite nuclei. Because of their phenomenological nature, these 
EDFs do not offer a way to systematically improve their predictive power.

For these reasons, one would like to combine the predictive power, systematic 
power-counting scheme, and connection to QCD of chiral potentials with the 
computational scaling, versatility and physics intuition of phenomenological 
potentials. One route to achieving this is based on the DME of 
expectation values <cit.>. Thanks to the Wick theorem, the expectation value of an 
arbitrary potential on a product state can be expressed as a functional of the 
one-body density matrix (or generalized density if pairing correlations are 
present). In the general case, the density matrix is fully non-local, that is, 
of the form ρ(rστ,r'σ'τ') with σ, 
τ the spin and isospin projections, respectively. The basic idea of the 
DME  is to expand ρ around the local density ρ(r) in 
order to turn the expectation value of the potential into a functional of the 
local density and gradient-like corrections. 

The method was first outlined by Negele and Vautherin in <cit.>. Several refinements to the original method to increase its 
accuracy were proposed in <cit.>. In <cit.>, the DME was 
applied to both the direct and exchange terms of the HF expectation value 
for unregulated momentum-space chiral potentials. The parameters of the resulting EDF were 
adjusted approximately using the singular value decomposition algorithm, and 
were tested in calculations of radii, single-particle spectra in doubly-closed 
shell nuclei and deformation energy. While promising, the authors reported 
numerical instabilities in the practical implementation of the DME, and 
emphasized that direct terms were treated in the local density approximation 
and that tensor contributions to the EDF had been neglected. 

The goal of this paper is to remedy some of these limitations, 
in particular by taking into account recent developments in 
chiEFT and employing high performance computing tools. Precisely, we want
to fully calibrate and validate an energy functional constrained by local 
chiral nuclear potentials. To this end, we build the EDF (in the ph 
channel) by computing the expectation value of chiral potentials on a Slater 
determinant. We use the DME of local chiral potentials presented in 
<cit.> to recast the exchange contribution in the form of a 
local functional of the density. We adjust the coupling constants of this 
EDF to ground-state properties of finite nuclei by solving the HFB 
equation. We provide results for EDFs corresponding to different orders 
in the chiral expansion up to N2LO order. Δ excitations and 
3N forces are included.

The paper is organized as follows. In section II, we recall the expressions for 
NN and 3N chiral potentials in coordinate space. Section III 
describes how these potentials can be transformed into an EDF with the 
DME (more details are given in <cit.>), 
and how these EDFs are implemented in current DFT solvers. We also 
provide in that section the result of our calibration process. In section IV, 
we test the predictive power of these EDFs on the equation of state of
nuclear matter, mass tables, sp energies of doubly-closed shell nuclei 
and the fission barrier of ^240Pu. Finally, we present some conclusions and 
perspectives in section V.






§ LOCAL CHIRAL POTENTIAL IN COORDINATE SPACE


We consider local chiral potentials up to N2LO with and without Δ 
excitations including also 3N forces <cit.>. 
Chiral interactions contain finite- and zero-range contributions, with the 
zero-range couplings usually fine-tuned to reproduce low energy π N and 
NN scattering data along with selected few-body properties and sometimes
properties of
nuclei up to Oxygen <cit.>. In the present work, we implement the finite-range 
contributions `as is', since they correspond to the long-range pion physics, 
which is well described by chiEFT. In contrast, the zero-range 
contribution will be replaced by a Skyrme-like potential, and we will take the 
contact coupling constants as adjustable parameters to be determined on 
selected properties of finite nuclei.

The finite-range contributions depend on a few set of parameters
including the pion mass m_π, the Δ-N mass splitting
M_Δ-N, the pion decay constant f_π, the nucleon axial
vector coupling g_A, the N-to-Δ axial vector coupling
h_A, as well as the LECs c_1, c_2, c_3,
c_4 and the LECs linear combination b_3 + b_8[For simplicity we retain 
the b_3 + b_8 combination even though it has been shown to be redundant <cit.>.]. 
Although the values for most of these parameters are well determined, the LECs have been
determined through different analyses of low energy π N and NN
scattering observables yielding different results
<cit.>. For this particular work, we use the
determination of Ref. <cit.> and leave the study of the
impact of the value of the LECs for a future work.

Since the finite-range potentials obtained from their corresponding diagrams 
diverge as r goes to zero, in practice a short range regulator is used in 
order to make scattering and structure calculations feasible. While several 
arguments have been made about the effects of the regulator on the statistical 
and systematic uncertainties <cit.>, for this work we restrict ourselves to the particular 
regulator

    f(r) = [1 - exp(-r^2/R_c^2) ]^n,

with R_c = 1.0 fm and n=6. We also leave the study of the dependence on the 
strength of the regulator controlled by the R_c and n parameters for a 
future work. The inclusion of the regulator is one of the improvements over the 
early work of <cit.>. In the two following subsections we recall the expressions for the 
long-range part only of the NN and 3N chiral potentials. 






 §.§ Two-Body Potential


The finite-range contribution to the local chiral potential in coordinate space 
is given by

    V_χ(r) 
        = 
    ( V_C(r) + W_C(r) τ_1 ·τ_2 ) 
        + 
    ( V_S(r) + W_S(r) τ_1 ·τ_2 ) σ_1 ·σ_2 
        + 
    ( V_T(r) + W_T(r) τ_1 ·τ_2 ) Ŝ_12(r̂),

where r ≡ |r|, Ŝ_12(r̂) is the usual tensor 
operator

    Ŝ_12(r̂) = 
    3 (σ_1·r̂)
      (σ_2·r̂) -σ_1 ·σ_2,

and σ_i (τ_i) is the spin (isospin) operator for
the ith particle. The potential components at LO
correspond to the well-known one-pion exchange and are given by

    W_S^ LO(r)    = m_π^3/12 π(g_A/2 f_π)^2 Y(r) 
    
    W_T^ LO(r)    = m_π^3/12 π(g_A/2 f_π)^2 Y(r) T(r)

where Y(r), U(r) and T(r) are the usual Yukawa, scalar and tensor 
functions, respectively,

    Y(r) = e^-x/x,  
    U(r) = 1 + 1/x,  
    T(r) = 1 + 3/xU(r),

with x = m_π r. The potential components at NLO including only nucleons and 
pions are given by


    W_C^ NLO(r) 
        = + m_π^5/8π^3(2f_π)^41/x^4{ 
    x [1+10g_A^2 -g_A^4(23+4x^2)]K_0(2x) 
    + [1+2g_A^2(5+2x^2) -g_A^4(23+12x^2) ]K_1(2x) 
    }, 
    
    
    V_S^ NLO(r) 
        = +m_π^5/2π^3( g_A/2fπ)^4 
    1/x^4[ 3 x K_0(2x) 
     + (3 +2x^2) K_1(2x) ], 
    
    
    V_T^ NLO(r) 
        = -m_π^5/8π^3( g_A/2fπ)^4 
    1/x^4[12 x K_0(2x)
    + (15 +4x^2) K_1(2x) ],


where K_0(x) and K_1(x) are the modified Bessel functions of the second 
kind. The potential components at N2LO including only nucleons and pions 
are given by


    V_C^ N2LO(r)     = 
    +3/2g_A^2 m_π^6/(2f_π)^4π^2e^-2x/x^6[2 c_1 x^2(1+x)^2 + c_3 (6+12x+10x^2+4x^3+x^4 ) ], 
    
    W_S^ N2LO(r)     = +1/3g_A^2 m_π^6/(2f_π)^4π^2e^-2x/x^6 c_4 (1+x)(3+3x+2x^2), 
    
    W_T^ N2LO(r)     = -1/3g_A^2 m_π^6/(2f_π)^4π^2e^-2x/x^6 c_4 (1+x)(3+3x+x^2).


The expressions for the potential components coming from one- and two-Δ 
excitations can be found in the supplemental material of 
Ref. <cit.>.






 §.§ Three-Body Potential


A general, local, three-body potential consists of all permutations with 
respect to the three two-body subsystems,

    V^ 3N = V_12 + V_23 + V_13,

where the V_ij potential depends on two of the relative coordinates and the 
spin and isospin of the three particles,

    V_ij = V(r_ik,r_jk,σ_1,τ_1,σ_2,τ_2,σ_3,τ_3),
      r_ij = r_i - r_j.

Because of the symmetry of the potential under subscript interchange
in the DME implementation (see Ref. <cit.>) only one
of the three terms in (<ref>) is necessary to express
the full three-body potential. For definiteness, we choose
V_23(r_21,r_31,{στ}).

The different terms appearing in the 3N chiral potential can be 
classified as (i) long-range, which are vertices of c_i or h_A and have no 
Dirac delta functions, (ii) intermediate-range, which are vertices of c_i, 
h_A or c_D and have one Dirac delta function, and (iii) short-range, which 
are vertices of c_i, h_A, c_D or c_E with two Dirac delta functions. In 
our implementation of the 3N chiral interaction in an EDF, we 
include only the long-range terms along with the c_i and h_A 
intermediate-range terms. We assume that all short-range terms and c_D 
vertices can be effectively absorbed by the optimization of the contact 
couplings on nuclear properties.

There is no contribution from 3N potentials up to N2LO, unless 
Δ excitations are included. The 3N potentials at NLO with 
Δ and N2LO have a very similar structure that can be summarized as

    V_3N^ NLOΔ    = ∑_i=1^3 α_i^ NLOΔ V_C,i
    
    V_3N^ N2LO    = ∑_i=1^3 α_i^ N2LO V_C,i + V_D + V_E.

We employ the NLOΔ label to emphasize that this contribution is 
only present when the Δ contributions are included. As mentioned before, 
the V_D (short-range term controlled by c_D) and V_E (short-range term 
controlled by c_E) terms are in fact not included in the present 
implementation. The α_i prefactors are given by

    α_1^ NLOΔ    = 0,   α_2^ NLOΔ   = -h_A^2m_π^6g_A^2/2592f_π^4π^2M_Δ-N, 
    
    α_3^ NLOΔ    =  h_A^2m_π^6g_A^2/10368f_π^4π^2M_Δ-N,

and

    α_1^ N2LO    = c_1m_π^6g_A^2/16f_π^4π^2,   α_2^ N2LO   = c_3m_π^6g_A^2/288f_π^4π^2, 
    α_3^ N2LO    = c_4m_π^6g_A^2/576f_π^4π^2,

and the V_C,i potentials are given by


    V_C,1   = 
    (τ_2 ·τ_3) 
    (σ_2 ·r̂_21)
    (σ_3 ·r̂_31) U(r_21) Y(r_21) U(r_31) Y(r_31) 
    
    
    V_C,2   = 
    (τ_2 ·τ_3)
    {16π^2/m_π^6 (σ_2 ·σ_3) δ^3(r_21) δ^3(r_31)
    - 
    4π/m_π^3[S_23 (r̂_21) T(r_21)
    + 
    (σ_2 ·σ_3) ] Y (r_21) δ^3(r_31) 
    .
       .
    - 
    4π/m_π^3[S_23 (r̂_31) T(r_31) 
    + 
    (σ_2 ·σ_3) ] Y (r_31) δ^3(r_21) 
    . 
       . 
    + 
    [ 
    9 (σ_2 ·r̂_21)(σ_3 ·r̂_31)(r̂_21·r̂_31) 
    - 
    3 (σ_2 ·r̂_21)(σ_3 ·r̂_21) 
    - 
    3 (σ_2 ·r̂_31)(σ_3 ·r̂_31)
    + 
    (σ_2 ·σ_3)  
    ]
    T(r_21) Y(r_21)  T(r_31) Y(r_31) 
    .
       .
    + (σ_2 ·σ_3) Y(r_21) Y(r_31) 
    + S_23(r̂_21) T(r_21)Y(r_21) Y(r_31),
    + S_23(r̂_31) T(r_31)Y(r_31) Y(r_21) }
    
    V_C,3    = 
    τ_2 · (τ_3 ×τ_1) 
    {16π^2/m_π^6δ^3(r_21) δ^3(r_31) σ_2 · (σ_3 ×σ_1) 
     .
       . 
    -12π/m_π^3 (σ_2 ·r̂_21)
    r̂_21· (σ_3 ×σ_1) T(r_21) Y(r_21) δ^3(r_31)
    +4π/m_π^3σ_2 · (σ_3 ×σ_1)
    3/m_π r_21 U(r_21)) Y(r_21) δ^3(r_31) 
    .
       . 
    -12π/m_π^3 (σ_3 ·r̂_31)
    r̂_31· (σ_1 ×σ_2) T(r_31) Y(r_31) δ^3(r_21)
    +4π/m_π^3σ_2 · (σ_3 ×σ_1)
    3/m_π r_31 U(r_31)) Y(r_31) δ^3(r_21) 
    .
       . 
    +9 (σ_2 ·r̂_21) (σ_3 ·r̂_31)
    σ_1 · (r̂_21×r̂_31) T(r_21)Y(r_21) T(r_31) Y(r_31)
    .
       . 
    -3 (σ_2 ·r̂_21)r̂_21· (σ_3 ×σ_1)
    T(r_21) 3/m_π r_31 U(r_31)Y(r_21)Y(r_31) 
    .
       . 
    -3 (σ_3 ·r̂_31) r̂_31· (σ_1 ×σ_2)
    T(r_31) 3/m_π r_21U(r_21)) Y(r_21) Y(r_31) 
    .
       . 
    -σ_1 · (σ_2 ×σ_3)
    3/m_π r_21U(r_21)3/m_π r_31U(r_31) Y(r_21) Y(r_31) 
    },


where the Yukawa Y(r), scalar U(r) and tensor T(r) functions are defined 
in (<ref>). As mentioned above, short-range terms with two Dirac 
delta functions are not included in the DME implementation of the 
3N interaction as we expect their effect to be absorbed by the 
calibration of the EDF contact couplings.






§ IMPLEMENTATION OF ΧEFT IN DFT



In our implementation of chiral interactions in the DFT framework, we 
will seek to write the total energy of a nucleus in the following form,

    E = 
    E_χ^ Har.
    +
    E_χ^ Skyrme
    +
    E^ Coul. + E^ pair.

The first term is the Hartree (direct) contribution to the expectation value of 
the long-range part of the local chiral potentials on Slater determinants 
reference states. In practice, we will see below that only the two-body chiral 
potential contributes to it because we only work with time-even systems. The 
second term is formally identical to an extended Skyrme-like functional and 
contains both the effects of the short-range part of the chiral 
potentials (in an effective way) and of the exchange contribution from the 
long-range part (through the DME). Finally, the terms 
E^ Cou and E^ pair are the usual Coulomb and pairing energy, 
obtained here by following exactly the same recipes as in <cit.>.

Two remarks are in order at this point:


  * It should be clear from Eq.(<ref>) that we restrict 
ourselves to the SR-EDF level. In other words, we seek to calibrate a 
functional built out of a single, HFB reference state in complete analogy 
with, e.g., the  family of functionals <cit.>, the  functional <cit.> or the  functional <cit.>. As a result, 
we expect our functional to be limited in its description of the fine structure 
of N=Z nuclei, for instance mirror displacement energies where isospin mixing 
and restoration are essential <cit.>) or the arc-like structure of 
binding energies near closed shell nuclei caused by quadrupole correlation 
energies <cit.>. We should also expect limitations in 
describing the shell structure of closed shell nuclei, where effects such as 
particle-vibration couplings should be taken into account <cit.>.

  * In this work, we take pairing functionals derived from a zero-range, 
surface-volume, two-body force as in the  functionals. The primary 
motivation for this choice is to focus on the effect of the DME on the 
ph channel only before considering its application to the pp 
channel. Early studies of pairing observables in finite nuclei with chiral 
potentials also suggest that it is mostly the short-range part of the latter 
that affect the pp channel <cit.>.
Finally, since we work at the SR-EDF level, the consistency of the 
generating kernels between the two channels is not really an issue.


To achieve the decomposition (<ref>), we express the 
expectation value of chiral potentials on a Slater determinant reference state. 
In Section <ref>, we briefly recall how this works for the two-body 
channel (<ref>): the Hartree term is expanded as a sum of 
Gaussians, while the Fock term is transformed into a generalized Skyrme 
functional with the DME. In Section <ref>, we 
give the expressions for the three-body channel, where only the Fock term
contributes; detailed derivations can be found in <cit.>.






 §.§ Two-Body Potentials



In configuration space, the contribution to the energy from a two-body 
potential reads

    E^ NN = 
    1/2∑_ij⟨ ij|  V_χ^ NN |kl ⟩ρ_kiρ_lj,

with ρ_ij the matrix elements of the one-body density matrix on an 
arbitrary basis of the single-particle Hilbert space. The two-body potential is 
antisymmetrized,

    V_χ^ NN = V_χ^ NN(1-P_σ P_τ P_r),

with the usual spin-, isospin- and space-exchange operators P_σ, 
P_τ and P_r

    P_σ≡1/2(1+σ_1·σ_2),     
     P_τ≡1/2(1+τ_1  ·τ_2  ).

The antisymmetrization operator results in direct and exchange contributions, 
also referred to as the Hartree and Fock energies respectively. 

After transforming (<ref>) to coordinate space by inserting 
resolutions of the identity, changing to relative ( r) and center of 
mass ( R) coordinates and assuming translational invariance along with 
a local potential, the two-body interaction energy term becomes

    E^ NN   = 1/2 Tr_1  Tr_2 ∫ d  R∫ d r ⟨ rσ_1 τ_1 σ_2 τ_2 | V_χ( r) |  rσ_3 τ_3 σ_4 τ_4 ⟩
       ×[ ρ_1 (  R +  r/2) ρ_2 ( R -  r/2) .
    
       .
    -ρ_1 ( R -  r/2,  R +  r/2)
     ρ_2 ( R +  r/2,  R -  r/2) P_12^στ],

where the traces refer to summation over spin and isospin quantum numbers 
and the local density matrix is

    ρ(x) ≡ρ(x,x).

The first term in Eq. (<ref>) corresponds to the Hartree 
energy, while the second one to the Fock energy. The following subsections 
describe our implementation of each of these two terms.  





  §.§.§ Hartree Term


The one-body density matrix in Eq. (<ref>) can be
decomposed into scalar-isoscalar, scalar-isovector, vector-isoscalar
and vector-isovector components <cit.>

    ρ(xσ_1τ_1,yσ_2τ_2) 
        = 
    1/4[ 
    ρ_0(x,y) + ρ_1(x,y)τ_z 
    
        +
    S_0(x,y) ·σ +
    S_1(x,y) ·στ_z 
    ].

Inserting this decomposition into the first term of 
Eq. (<ref>) and performing the traces one obtains

    E_ H = 
    1/2∑_t=0,1∫ d r∫ d R [
    ρ_t(R^+) ρ_t(R^-) Γ_ρρ^t 
    
      + 
    S_t(R^+) ·S_t(R^-) Γ_SS^t 
    
     + 
    (S_t(R^+) ·r̂ ) (S_t(R^-) ·r̂ ) Γ_Sr^t ],

where R^± = R± r/2 and

  
    Γ_ρρ^t    =
        
          V_C      for   t = 0 
    
          W_C      for   t = 1
        
    Γ_SS^t    =
        
          V_S - V_T      for   t = 0 
    
          W_S - W_T      for   t = 1
        
    Γ_Sr^t    =
        
          3   V_T      for   t = 0 
    
          3   W_T      for   t = 1

  

with t=0 (t=1) indicating the isoscalar (isovector) case. Note that for 
systems with time-reversal symmetry, all terms diagonal in the spin density 
vanish, i.e., S(x) = 0. Hence for even-even nuclei, only the terms 
proportional to the central part of the potential in Γ_ρρ^t 
contribute to the Hartree energy.



At this stage the DME could be applied to the calculation of the Hartree 
energy. However, it has been established that computing this term exactly 
provides a more precise description of the density fluctuations and energy 
contribution <cit.>. Furthermore, the inclusion of 
the DME approximation in the calculation of the Hartree field introduces 
numerical instabilities <cit.>. For these reasons, we choose to 
compute the direct term exactly. To compute the matrix element of the finite 
range of the chiral potential, we expand the Yukawa form factors on a series of 
Gaussian functions <cit.>. This allows us to take full advantage of the 
many analytic properties of Gaussian matrix elements in the Harmonic Oscillator basis 
<cit.> and of the existing implementation of the Gogny force in 
the latest version of  <cit.>. As we show in 
figure <ref>, five Gaussian functions already give an excellent 
approximation to the spatial part of the potential.  

The chiral potential in Eq. (<ref>) is expressed in a 
spin-isospin operator basis. In the code  <cit.>, 
the finite range part of the Gogny functional is implemented in a exchange 
operator basis, that is,

    V_ G = ∑_i=1^N (W_i + B_i P_σ + H_i P_σ P_τ + M_i P_σ P_τ)
    e^-r^2/μ_i^2.

After inserting the definitions of Eq. (<ref>) and 
rewriting the terms in the Gogny functional, we find

    V_ G = ∑_i=1^N 
       [ W_i + B_i/2 + H_i/2 + M_i/4. 
       . + ( B_i/2 + M_i/4) σ_1 ·σ_2  
           + ( H_i/2 + M_i/4) τ_1 ·τ_2 
    .
    
       .
    + M_i/4 (σ_1 ·σ_2)(τ_1 ·τ_2) ]e^-r^2/μ_i^2.

Since only systems with time-reversal symmetry are being considered for 
this work, only the central components of the Chiral potential are considered 
in Eq.(<ref>). Therefore, we can set B_i = M_i = 0 and use 
the approximations

    V_C(r) →Ṽ_C(r)    = ∑_i=1^N ( W_i + H_i/2) e^-r^2/μ_i^2, 
    
    W_C(r) →W̃_C(r)    = ∑_i=1^N H_i/2 e^-r^2/μ_i^2.

In order to reproduce the behavior of the regulator in Eq. (<ref>) 
the conditions

    H_N = - ∑_i=1^N-1 H_i,     
      W_N = - ∑_i=1^N-1 W_i,

are imposed. This conditions ensure that the potentials vanish as 
r→ 0. The remaining free parameters W_i, H_i and μ_i are 
adjusted numerically to reproduce the central components of the chiral 
potential at different orders. In figure <ref> we show the 
difference between the approximations as a sum of five Gaussian functions and 
the corresponding chiral potential up to a certain order in the chiral 
expansion. Note that while the scale in figure <ref> is
10^-2MeV the potentials have an order of magnitude, at their highest values,
of 10 MeV.






  §.§.§ Fock Term


Inserting Eq. (<ref>) into the second term of 
Eq. (<ref>), the two-body Fock term becomes

    E_ F   = -1/2∑_t=0,1∫ d r∫ d R[ 
    ρ_t^2(R^+,R^-) Ξ_ρρ^t . 
       . 
    - S_t^2(R^+,R^-) Ξ_SS^t 
    -(S_t(R^+,R^-) ·r̂ )^2  Ξ_Sr^t ],

where the symmetries ρ_t(x,y)=ρ_t(y,x) and
S_t(x,y)=-S_t(y,x) for time-reversal invariant 
systems have been used and the Ξ^t functions are given by

  
    Ξ_ρρ^t    =
        1/4 V_C + 3/4 W_C + 3/4 V_S + 9/4 W_S     for   t = 0 
    1/4 V_C - 1/4 W_C + 3/4 V_S - 3/4 W_S     for   t = 1
        
    Ξ_SS^t    =
        1/4 V_C + 3/4 W_C - 1/4 V_S - 3/4 W_S 
    - 1/2 V_T - 3/2 W_T     for   t = 0 
    1/4 V_C - 1/4 W_C - 1/4 V_S + 1/4 W_S 
    - 1/2 V_T + 1/2 W_T     for   t = 1
        
    Ξ_Sr^t    =
        3/2 V_T + 9/2 W_T     for   t = 0 
    3/2 V_T - 3/2 W_T     for   t = 1.


The DME consists in expanding the non-diagonal density matrices in such manner 
that the non-locality is factorized using the following formula,

    ρ_t   (R^+,R^-)    ≈∑_n=0^n_ maxΠ_n^ρ (kr) 𝒫_n (R), 
    S_t (R^+,R^-)    ≈∑_m=0^m_ maxΠ_m^s (kr) 𝒬_m (R)

where the Π functions are specified by the DME variant and 
𝒫_n (R), 𝒬_m (R) denote various local 
densities. The arbitrary momentum scale k in the Π functions sets the 
scale for the fall-off in the off-diagonal direction. In this work, we follow 
common practice and truncate the expansion at n_ max = 2 and 
m_ max = 1 such that

    ρ_t(R^+,R^-)    ≈Π_0^ρ(k_ Fr) ρ_t(R) 
        + r^2/6Π_2^ρ[1/4Δρ_t(R) - τ_t(R) + 3/5 k_ F^2 ρ_t(R) ], 
    
      S_t,b(R^+,R^-)    ≈ i Π_1^s (k_ F r) ∑_a=x^z r_a J_t,ab(R) ,

where the kinetic density τ_t and spin current density J_t are defined as

    τ_t(r)     = ∇·∇' ρ_t(r,r')|_r=r', 
    
    J_t,ab(r)     = -i/2(∇_a - ∇_a')S_t,b(r,r')|_r=r'.

In Gebremariam's improved phase-space-averaging DME variant
<cit.>, the momentum scale k is 
chosen to be the Fermi momentum k_ F with the Π functions given by

    Π_0^ρ (k_ F r) 
    = Π_2^ρ (k_ F r) 
    = Π_1^s (k_ F r) 
    = 3 j_1 (k_ F r)/k_ F r

where j_1 is a spherical Bessel function of the first kind and k_ F is 
related to the isoscalar density in the usual way,

    k_ F = ( 3 π^2/2ρ_0 (R) )^1/3.

By inserting the expansions of 
Eqs. (<ref>)-(<ref>) into the exact 
Fock energy of Eq. (<ref>), the Fock energy can be approximated 
by expressions involving only products of local densities. Terms beyond 
second-order in the density expansions are dropped e.g., 
Π_2^ρ(k_ Fr) Π_2^ρ(k_  Fr). After performing the DME and 
organizing the different terms in Eq. (<ref>) by the different 
densities, we find the more compact expression

    E_ F≈   ∑_t=0,1∫ d R [ 
    g_t^ρρρ_t ρ_t + g_t^ρτρ_t τ_t + g_t^ρΔρρ_t Δρ_t .
    
       . + g_t^JJ,1 J_t,aa J_t,bb + g_t^JJ,2 J_t,ab
                 J_t,ab + g_t^JJ,3 J_t,ab J_t,ba]

where the R-dependence of the local densities and couplings g has been 
omitted for simplicity and the coupling functions are given by


    g^ρρ_t(ρ_0)     = 
        [                       -4π/2∫ dr r^2 [ Π_0^ρ(k_ Fr)^2; + r^2 k_ F^2/5Π_0^ρ(k_ Fr) Π_2^ρ(k_ Fr) ] Ξ_ρρ^t(r), ]
    
    g^ρτ_t(ρ_0)     = 4π/2∫ dr r^2 [ 
      r^2/3Π_0^ρ(k_ Fr) Π_2^ρ(k_ F) ] Ξ_ρρ^t(r),
    
    
    g^ρΔρ_t(ρ_0)    = -4π/2∫ drr^2 [ 
      r^2/12Π_0^ρ(k_ Fr) Π_2^ρ(k_ F) ] Ξ_ρρ^t(r),
    
    
    g^JJ,1_t(ρ_0)    = -4π/2∫ dr r^2 [ 
      r^2/15Π_1^S(k_ Fr)^2 ] Ξ_Sr^t(r),
    
    
    g^JJ,2_t(ρ_0)    = -4π/2∫ dr r^2 [ 
      r^2/15Π_1^S(k_ Fr)^2 ] [ 5 Ξ_SS^t(r) + Ξ_Sr^t(r) ],
    
    
    g^JJ,3_t(ρ_0)    = -4π/2∫ dr r^2 [ 
      r^2/15Π_1^S(k_ Fr)^2 ] Ξ_Sr^t(r).



As already highlighted in <cit.>, one of the practical differences 
between Skyrme and DME-based functionals is that each Skyrme coupling constant 
becomes a coupling function, which is dependent on the isoscalar density. The 
calculation of these density-dependent couplings requires performing several
multidimensional numerical integrals, some of them converging slowly at small 
values of ρ_0. However, these couplings are completely independent of the 
system being calculated or any other characteristic of the HFB simulation like 
the basis size or oscillator length. Therefore, we can tabulate all the 
relevant coupling functions for different values of the density ρ_0. In 
the actual HFB calculation the couplings are approximated via the interpolating function




    g̃^uv(ρ) = g^uv(0) + ∑_i=1^N a_i tan^-1(b_i ρ_0^c_i).


The parameters of the interpolating functions were adjusted to reproduce
the tabulated values using N=3. While other forms of interpolating
functions were considered, this one gave a better description of the
coupling functions while avoiding numerically unstable behavior at
small and large values of ρ_0.  An additional advantage of using
interpolating functions is that the inclusion of Δ excitations
and 3N forces does not imply any increase on computational cost
since the same type of interpolating function is used for all
cases. Figure <ref> shows the numerical
precision of these interpolating functions for two-body
couplings. While the numerical precision of the interpolation is of the order 
of 10^-3 MeV fm^-3, this is 4 orders of magnitude smaller than the 
scale of the couplings themselves, which are ∝ 10 MeV fm^-3.






 §.§ Three-Body Term



The contribution of the V_χ^ 3N three-body chiral potential to the 
total energy is given by

    E^ 3N = 
    1/6∑_ijk⟨ ijk | V_χ^ 3N𝒜_123 | lmn ⟩ρ_liρ_mjρ_nk,

with 𝒜_123 = (1+P_13 P_12 + P_23 P_12)(1-P_12) the full 
three-body antisymmetrization operator. Since all the terms in the 3N 
Hartree energy contain at least one spin density matrix, which vanishes in 
time-reversal invariant systems, there is in fact no contribution from the 
3N channel to the Hartree energy. The application of the DME to the 
Fock term results in 23 trilinears of local densities, each one with its 
corresponding density-dependent coupling,

  
    E_ F^ 3N   ≈∫ dR(
        g^ρ_0^3ρ_0^3 
        + g^ρ_0^2τ_0ρ_0^2τ_0 +
        g^ρ_0^2Δρ_0ρ_0^2Δρ_0 
        + g^ρ_0(∇ρ_0)^2ρ_0∇ρ_0·∇ρ_0 
        + g^ρ_0ρ_1^2ρ_0ρ_1^2 
        + g^ρ_1^2τ_0ρ_1^2τ_0 
        + g^ρ_1^2Δρ_0ρ_1^2Δρ_0 
        . 
        + g^ρ_0ρ_1τ_1ρ_0ρ_1τ_1
        + g^ρ_0ρ_1Δρ_1ρ_0ρ_1Δρ_1 
        + g^ρ_0(∇ρ_1)^2ρ_0∇ρ_1·∇ρ_1 
        + ρ_0 ϵ_ijk[
          g^ρ_0∇ρ_0J_0∇_iρ_0J_0,jk + g^ρ_0∇ρ_1J_1∇_iρ_1J_1,jk] 
        
        + ρ_1 ϵ_ijk[
          g^ρ_1∇ρ_1J_0∇_iρ_1J_0,jk + g^ρ_1∇ρ_0J_1∇_iρ_0J_1,jk]
        + ρ_0 [
          g^ρ_0J_0^2,1J_0,aaJ_0,bb + g^ρ_0J_0^2,2J_0,abJ_0,ab + g^ρ_0J_0^2,3J_0,abJ_0,ba] 
        
        + ρ_0 [
          g^ρ_0J_1^2,1J_1,aaJ_1,bb + g^ρ_0J_1^2,2J_1,abJ_1,ab + g^ρ_0J_1^2,3J_1,abJ_1,ba] .
        
        + ρ_1 . [
          g^ρ_0J_0J_1,1J_1,aaJ_0,bb + g^ρ_0J_0J_1,2J_1,abJ_0,ab + g^ρ_0J_0J_1,3J_1,abJ_0,ba]
        ).





We refer the reader to <cit.> and its supplemental material for 
an in-depth derivation of the 3N energy Fock term and complete expressions for 
the density-dependent couplings. Similarly to the two-body case, the 
calculation of the density functions in terms of the density ρ_0 requires 
several multidimensional numerical integrals with slow convergence. To avoid 
calculating these coupling functions at every iteration of the HFB 
calculation, we employ the same type of interpolating function as in 
Eq. (<ref>) and adjust its parameters to reproduce 
tabulated couplings. Figure <ref> shows the accuracy of 
these interpolations at different chiral orders for 2 representative functions, 
g^ρ_0^3 and g^ρ_0ρ_1^2. The accuracy of the interpolation is 
comparable for all coupling functions. The irregularities in the curves show 
that the interpolation has an accuracy similar to the numerical multidimensional 
integral. As with two-body couplings, the numerical error of the interpolation 
is about 4 orders of magnitude smaller than the couplings themselves.










 §.§ Optimization of Contact Couplings


Unlike ab initio methods in which many-body correlations are encoded into 
the nucleon wave-functions, DFT assumes independent (quasi-)particles through 
uncorrelated, product wave-functions. Therefore, many-body correlations have to 
be included both through the form and parameters of the EDF (single- or 
multi-reference) and through the symmetry-breaking mechanism. Recall that in 
our decomposition (<ref>) of the total energy, the Skyrme-like 
part is the sum of three terms: an effective two-body Skyrme functional that 
mocks up the effects of the short range part of the chiral potentials; a 
generalized Skyrme functional of the form (<ref>) that contains the 
exchange contribution of the two-body chiral potential; and a generalized 
Skyrme functional of the form (<ref>) for the exchange 
contribution of the 3N channel. The two-body part thus reads

    E_ NN^ Skyrme = ∑_t=0,1∫ d R[ 
    ( 
    U_t0^ρρ + 
    U_tD^ρρρ_0^γ) ρ_t^2 + 
    U_t^ρτρ_t τ_t .
    . + 
    U_t^ρΔρρ_t Δρ_t +
    U_t^ρ∇ Jρ_t ∇·J_t + 
    U_t^JJ J_t,ab J_t,ab]

where each coupling function is given by 

    U_t^uu'(ρ) = C_t^uu' + g_t^uu'[ρ_0(R)]

with the functions g_t^uu'[ρ_0(R)] listed in 
Eqs. (<ref>)-(<ref>). The pairing energy is given by

    E^ pair = 1/4∑_q=n,p∫ d R  
    V_0^q [ 1 - 1/2ρ_0(R)/] ρ̃^2(R)

where ρ̃(R) is the pairing density and =0.16 
fm^-3. The coefficients C_t^uu' and V_0^q are the unknown 
parameters that we will determine in the calibration process. 

In practice, note that we are only fitting parameters in the NN exchange and 
pairing channel of the functional. The two-body Hartree term, 
Eq. (<ref>), is computed “exactly” (that is, without any 
adjustment of parameters) and so is the three-body Fock term, 
Eq. (<ref>). 

For the optimization of the contact couplings, we follow the same prescription 
as for the  parametrization of the Skyrme functional 
<cit.>. Among the fourteen parameters in 
Eqs. (<ref>) and (<ref>), C_t0^ρρ, 
C_tD^ρρ, C_t^ρτ and γ are volume couplings and, 
therefore, can be directly related to INM properties, which allows using 
tighter, physically-motivated bounds; see Section IV.A-C and Appendix C in 
<cit.> for actual expressions relating INM properties with 
the coupling functions (<ref>). In practice, we thus optimize the 
following INM parameters: , , , ,  and 
. As with all  parametrizations of the Skyrme functional, we do 
not optimize the vector effective mass  but instead keep it fixed at its 
SLy4 value, = 1.249. The remaining eight parameters are fitted directly 
with the same bounds as for .

Note that, in contrast to Skyrme EDFs, our DME EDFs contain a 
finite-range term that contributes to INM properties. The contribution of 
this finite-range term to the energy per particle for asymmetric nuclear matter 
reads 
<cit.>


    e(ρ,β) = 1/2∑_i=1^N { A_0^i ρ + A_1^i ρβ^2  + B_nn^i [1+β/2 g(μ_i k_F^n) + 1-β/2 g(μ_i k_F^p)  ] + B_np^i h(μ_i k_F^n,μ_i k_F^p) },


where ρ = ρ_n + ρ_p is the total density, β=(ρ_n-ρ_p)/ρ the isospin 
asymmetry parameter, k_F^n = k_F(1+β)^1/3 and 
k_F^p = k_F(1-β)^1/3 are the Fermi momenta of the corresponding
isospin symmetric system, the parameters A and B are 

    A_0^i    = π^3/2μ_i^3/4( 4 W_i + 2 B_i - 2 H_i - M_i),
    
      A_1^i    = π^3/2μ_i^3/4( - 2 H_i - M_i), 
    
      B_nn^i    = - 1/√(π)(W_i + 2B_i- H_i -2M_i ), 
    
      B_np^i    =   1/√(π)(H_i  + 2M_i ),

and the functions functions g(q) and h(q_1,q_2) are the result of a double 
integration of the exchange matrix elements over the same Fermi surface for 
g, and two different surfaces for h, which result in

    g(q)    = 2/q^3 - 3/q - (2/q^3 - 1/q) e^-q^2 + √(π) erf(q) 
    
      h(q_1,q_2)    = 2 q_1^2 -q_1q_2+q_2^2-2/q_1^3+q_2^3 e^-(q_1+q_2)^2/4
       - 2  q_1^2 +q_1q_2+q_2^2-2/q_1^3+q_2^3 e^-(q_1-q_2)^2/4
        -√(π)q_1^3-q_2^3/q_1^3+q_2^3 erf(q_1-q_2/2) + √(π) erf(q_1+q_2/2)




The calculation of the INM properties simply requires performing a Taylor 
expansion around the saturation density  which yields

    e(ρ,β)    = e(ρ) + S_2(ρ) β^2 + S_4(ρ) β^4 + …,
    
      e(ρ)    = /A + /δρ  + /18^2(δρ)^2 + …, 
    
      S_2(ρ)    =  + /3 + δρ + ….

For detailed derivations of the finite-range contributions to the infinite 
nuclear matter properties see <cit.>.

The new EDFs were optimized from LO to N2LO. The inclusion of 
a Δ excitation produces two additional versions of the functional at 
NLO and N2LO, denoted as NLOΔ and N2LOΔ. 
Finally the incorporation of 3N forces at the appropriate orders adds three 
more versions denoted as N2LO+3N, NLOΔ+3N and 
N2LOΔ+3N. Each of these 8 EDFs has its own set of 
Gaussian functions to represent the finite-range contribution to the Hartree 
field, density-dependent couplings for the two- and three-body Fock fields and 
calibrated two-body contact couplings to recover many-body correlations.

As already highlighted in <cit.>, the parameter space of DME 
EDFs could be significantly different from that of traditional Skyrme 
EDFs. To avoid possible difficulties during the optimization process 
(which was initialized with the  parameter set), we took advantage 
of the built-in regulator (<ref>): for large values of R_c, 
f(r) → 0 and the finite-range contributions vanish. Therefore, the 
EDF reduces to a traditional Skyrme EDF. Starting from the 
 parameter set, we thus produced intermediate parametrizations of 
all the DME EDFs at R_c = 2.0 fm. As mentioned earlier, the final 
parametrizations were obtained with R_c = 1.0 fm. 

It is important to note that the Δ-less and Δ-full version of the 
chiral potential employ different sets of low-energy constants. In this work we 
derive the EDFs from the local potentials in coordinate space as 
presented in <cit.>, for which the DME approximation has 
been applied in <cit.>. The values for the LECs and other 
physical parameters used for the chiral potentials in this work are listed in 
Table I of <cit.>.

The actual parameter sets of all 8 DME EDFs are listed in 
Table <ref>. The optimization was carried out with the  
optimization package from Argonne National Laboratory, with all HFB 
calculations performed with the  solver with the exact same basis 
characteristics as in <cit.>. We notice that for most DME
EDFs, the scalar effective mass ends up at its bound. As expected, there 
are substantial variations among the different parametrizations, in particular 
when it comes to γ, C_0^ρΔρ and C_1^JJ. Large 
fluctuations in the isovector channel are not surprising, since it is widely 
believed that the lack of constraints on these parameters comes from a lack 
experimental data in very neutron-rich nuclei. The observed fluctuations in the 
power of the density-dependence and the isoscalar surface terms are indicative 
of the strong non-linearity of the optimization process. With the exception of 
the EDF at LO, the value of the objective function is similar for all 
EDF at around 164±5, but we will see in the next section that there 
are significant differences in predictive powers.








§ VALIDATION AGAINST EXPERIMENTAL DATA







 §.§ Infinite nuclear matter properties




We calculated the EOS with each of the new microscopically constrained 
EDFs. As mentioned in the previous section, the inclusion of 
density-dependent couplings and finite-range contributions in the density 
functional brings additional terms to the corresponding EOS. These terms 
were included following the derivations in <cit.>. In Fig. <ref> we show the energy per 
nucleon E/A as a function of the density ρ for SNM and PNM.

Since the value of the saturation density and other INM properties at 
saturation were used to constrain the contact couplings for the EDFs, it 
is not surprising that all curves exhibit very similar behavior around the 
saturation point. The curves for the different EDFs start to deviate from 
one another at large values of ρ, specially for cases in which the 
3N terms are not included even though the corresponding diagrams are 
present at such order. For the Δ-less implementation (left panels), a 
convergence pattern can be seen when including the 3N terms, i.e., the 
difference between LO and NLO is larger than the difference between 
NLO and N2LO+3N. Unfortunately such a convergence pattern can 
not be found in the Δ-full implementation.

As a reference point we include the recent calculation by Logoteta, Bombaci and 
Kievsky (LBK) for a local chiral potential with Δ 
isobar <cit.>. While a direct comparison can not be done with 
our current results since the LBK calculations correspond to a Δ-full 
implementation of the 2-body force at N3LO and the 3N terms at 
N2LO, it still provides a useful reference. Also for comparison, the recent
auxiliary field diffusion Monte Carlo simulations <cit.> using a
chiral interaction at N2LO give an energy per nucleon for SNM between
12 and 16 MeV at saturation density depending on the value of the regulator. These
simulations are in agreement with our calculations. Overall, our parametrization 
seems to yield a stiffer EOS than fully microscopic calculations.






 §.§ Nuclear Mass Tables


For each of the EDF listed in Table <ref>, we computed the 
binding energies of all even-even nuclei from Z=8 to Z=120. The driplines 
were identified by the requirement that the two-neutron separation energy 
change sign and become negative. For each even-even nucleus, we considered 11 
different configurations characterized by their axial quadrupole deformation 
β_2 = -0.25, …, +0.25 (by steps of 0.05). For each configuration, we 
used the small-deformation approximation of the quadrupole moment, 
Q_2 ≈β_2 √(5/π)(Z+N)^5/3/100 (in barns), to impose a 
constraint on Q_2 for the HFB solution. The constraint was only active 
during the first 20 iterations of the self-consistent loop and was then 
automatically released. The binding energy retained for the even-even nucleus 
is then the lowest energy of these 11 configurations.



Binding energies of odd-even and odd-odd nuclei were not computed explicitly, 
as it would require performing numerous blocking calculations. There are at 
least two reasons why such an explicit calculation is not mandated here: (i) 
the main focus of this work is a global assessment of microscopically-constrained 
EDFs, not the production of a mass model, and (ii) our optimization protocol, 
also used for the Skyrme  functional, does not put special emphasis 
on nuclear masses, which are 1 of 5 different types of observables. For these 
reasons, we rely instead on a popular approximation, where the energy of an odd-even nucleus 
(Z-1,N) (with both Z and N even-even) is given by

    E(Z-1,N) = 
    1/2[ E(Z,N) + E(Z-2,N) ] 
    
    + 
    1/2[ Δ_p(Z,N) + Δ_p(Z-2,N) ],

with Δ_p(Z,N) the average proton pairing gap (obtained as 
Δ(Z,N) = 1/2TrΔρ_p, where Δ is the 
pairing field of the HFB matrix). Similar formula hold for even-odd 
nuclei. For odd-odd systems, we first compute

    Δ_p(Z,N-1) = 1/2[ Δ_p(Z,N) + Δ_p(Z,N-2) ]

and, similarly, Δ_p(Z-2,N-1), and combine them to get

    E(Z-1,N-1) = 
    1/2[ E(Z,N-1) + E(Z-2,N-1) ] 
    
    + 
    1/2[ Δ_p(Z,N-1) + Δ_p(Z-2,N-1) ]

All calculations were performed with the code  in a deformed (stretched) 
basis of 20 shells with an axial deformation β = β_2.



Table <ref> summarizes the characteristics of the mass tables for 
each EDF. It lists the rms deviation between theoretical and 
experimental nuclear binding energies as well as the number of experimental 
measurements. Experimental atomic masses are taken from the 2016 Atomic Mass 
Evaluation <cit.>. Nuclear binding energies are obtained 
after taking into account the binding energy of the electrons. Following 
<cit.>, we adopt the following empirical formula

    B_e(Z) = 1.44381× 10^-5Z^2.39 + 1.55468× 10^-12 Z^5.35

with the energy given in MeV. We only included true experimental measurements 
and did not take into account evaluated masses. Further details on how nuclear 
binding energies are extracted from the mass evaluation can be found in 
<cit.>.

Perhaps the most surprising (and promising) result is the relatively large 
variation of the results, with a rms ranging from 1.26 MeV for 
N2LOΔ to 2.02 MeV for NLO – a 60% difference in predictive 
power. It is also very encouraging to note that the EDFs seem to perform 
better and better overall as we go from LO to NLO to N2LO and 
add Δ excitations. In fact, the quality of the N2LO+Δ 
EDF is rather spectacular. Without any “beyond mean-field” corrections 
such as the rotational or vibrational corrections, Wigner energy, etc., this 
EDF does markedly better than . Remember that (i) the 
determination of the EDF parameters was made with the exact same protocol 
and optimizer (ii) all mass tables were computed with the exact same code, 
basis characteristics and overall algorithms to identify driplines. Therefore, 
the origin of all differences listed in Table <ref> can be 
attributed to the form of the EDF only.

A visual representation of the difference between theory and experiment 
highlights a few additional interesting features of these EDFs. For light 
nuclei, Fig. <ref> shows that the LO and NLO EDFs 
behave like the  (and older ) EDF: discrepancies 
with experimental masses are larger, which was explained in 
<cit.> as resulting from neglecting the 
center-of-mass correction in the EDF – a choice that we also made for 
all the DME EDFs. Surprisingly, this feature is much attenuated for 
EDFs based on higher-order chiral potentials. 

We also notice that both EDFs including the effect of three-body force 
and, to a lesser extent, that of the Δ excitations have more pronounced 
spikes near closed-shell nuclei, as shown by comparing, e.g., the mass tables 
for N2LO and N2LO+3N, or N2LOΔ and 
N2LOΔ+3N. Overall, we also notice that the effect of the 
three-body force seems to be the largest near closed shells, in particular near 
^208Pb.



Table <ref> completes the picture by showing the mean value 
and standard deviations computed from the residuals of nuclear binding 
energies. Compared with  and LO (which is not much different 
from  by construction), DME functionals have a larger 
systematic bias – which also tends to decrease as we go to higher order in the 
chiEFT expansion. Conversely, the standard deviation for the DME 
functionals is much smaller than for Skyrme, and the trend is also towards 
smaller standard deviations. Recall that for a random variable with mean μ 
and standard deviation σ, we have rms^2 = σ^2 + μ^2. 



Another indicator of the global quality of a EDF is the proton radius. We 
show in Fig.<ref> the residuals for proton radii for the 339 nuclei 
listed in <cit.> in the two particular case of the  and 
the NLOΔ+3N functionals. We extracted experimental proton 
radii from the table of <cit.> by using the formula

    r_ ch^2 = ⟨ r_p^2⟩ + r_p^2 + N/Zr_n^2

with r_ ch^2 the charge radius of the nucleus, 
⟨ r_p^2⟩ the expectation value of the proton radius on the 
HFB solution, and r_n^2 = -0.1161 fm^2 and r_p^2 = 0.769 
fm^2 the charge mean square radii of the free neutron and proton, 
respectively. Table <ref> lists the rms, mean value and 
standard deviation of the residuals for all 8 functionals listed in Table 
<ref>. Overall, the prediction of proton radii is on par with competing 
functionals, see, e.g., <cit.>, although it is slightly 
worse than for the Skyrme . We also observe a similar effect as for 
masses: DME functionals have a larger systematic bias than the pure 
Skyrme . However, this bias is very small (<0.01 fm) and may not 
be very significant.








 §.§ Shell Structure


We turn to the sp shell structure of closed shell nuclei. As a reminder, 
we extract sp energies of the nucleus (Z,N) by performing blocking 
calculations <cit.> at the equal filling approximation <cit.> 
in the neighboring odd nuclei, e.g., (Z,N± 1) for neutrons sp states; 
see, e.g., <cit.> for studies of the blocking prescription on the ground-state 
properties of odd nuclei. Specifically, we define 


    E_ s.p.^(part.)    =     E_ bl(A+1) - E(A), 
    
    E_ s.p.^(hole)    =     E(A) - E_ bl(A-1),


where A is the particle number of the reference, doubly-magic, nucleus of 
interest and E_ bl is the energy of the blocked configuration in the
neighboring odd nucleus. The labels “hole” and “particle” refer to whether
the corresponding sp levels would be, respectively, fully occupied or 
empty in the corresponding HF calculation of the doubly-magic nucleus. 

This method presents two advantages. First, it ensures the consistency of the 
calculations for all observables. Whether we consider masses, sp energies 
or fission barriers, we always perform computations in the same HFB framework 
with the Lipkin-Nogami correction. Second, we automatically include the small 
shape polarization induced by the blocking calculation – even though this 
polarization is restricted here to axial shapes owing to the built-in 
symmetries of  <cit.>.

We recall that in , blocking configurations can only be specified by the 
Nilsson quantum numbers [Nn_zΛ]Ω of the requested sp state; 
see <cit.> for details. Since these quantum numbers are only valid 
approximately (Ω corresponds to a conserved symmetry of the mean field, 
but not the others; see <cit.> for a discussion), the convergence of 
the blocking calculations can sometimes fail. In particular, we found that for 
low-j orbitals, we had to introduce tiny constraints either on the expectation 
value of Q̂_2 or Q̂_4 in order to converge the blocking 
calculations. Since the effective sp energy is defined as an energy 
difference, the numerical error introduced is very small – less than 50 keV 
overall. Note that similar difficulties were experienced with the  
family of EDFs presented in <cit.>.

The figure <ref> shows the example of neutron sp states in 
^208Pb for the various EDFs listed in Table <ref>. Contrary 
to binding energies, we do not observe a very clear improvement or degradation 
of the shell structure as a function of the EDF used. The 
sp spectrum in other closed-shell nuclei yields similar conclusions. This 
could be attributed to the fact that the optimization protocol is the same for 
all EDFs and explicitly include a constraint on a few spin-orbit 
splittings. In addition, work with either Skyrme EDF or covariant 
DFT suggest that correlations such as particle-vibration couplings play a 
major role in improving the shell structure in closed shell nuclei. It is  
unlikely that the DME functionals we consider have this type of 
correlations built-in.








 §.§ Deformation Properties


As mentioned in the introduction, one of the primary applications of DFT 
is the study of nuclear structure and excited states for deformed nuclei, 
including rotational and vibrational spectra and fission. Although the shell 
structure of closed-shell nuclei discussed in the previous section and the 
large-scale mass tables shown in Figure <ref> are indicative of an 
EDF with a good overall predictive power, it is important to also test 
the behavior of the EDF at large deformations. For this reason, the 
 optimization protocol includes the excitation energy of a few 
fission isomers in actinides. As discussed in <cit.>, this provides constraints on both the shell structure – 
inasmuch as deformation properties are partly driven by the particular ordering 
and level density of sp j-shells in spherical nuclei – and on surface 
properties of the EDF, which are related in particular to .





We report in Table <ref> the excitation energy of the fission 
isomer and height of the first and second barriers extracted from calculations 
of the potential energy curve in ^240Pu. Across all 8 DME EDFs, 
the rms deviation for the excitation energy of the fission isomer is 
0.29 MeV, which is comparable to the predictive power of the  Skyrme 
functionals across all actinides; see <cit.>. However, fission 
barriers tend to be too high – even when taking into account the extra 
≈ 2.5 MeV caused by the lack of triaxiality in calculations of the 
first barrier, column marked E_A (est.) in the table. For example, the 
rms deviation for the second barrier is 1.55 MeV, compared with 1.39 MeV 
for  and 0.69 MeV for  (across all actinides). We show in 
the two panels of Figure <ref> the potential energy curve of 
^240Pu as a function of the axial quadrupole moment for all the 8 
EDFs considered here and listed in Table <ref>, together with 
that of the ,  and  functionals for 
comparison. 








§ CONCLUSIONS


In this work, we have calibrated and validated a set of energy density
functionals derived from local chiral potentials through the DME.
We have provided the parametrization of these EDFs at LO, 
NLO, N2LO, N2LO+3N, NLOΔ, 
NLOΔ+3N, N2LOΔ, N2LOΔ+3N for 
a value of R_c = 1.0 fm and n = 6 in the regulator cutoff. The optimization 
was performed with the  protocol, and results were validated on the 
EOS of infinite nuclear matter and pure neutron matter, nuclear mass 
tables, the shell structure of doubly-closed shell nuclei and the deformation 
energy of ^240Pu.

The overall predictive power of these EDFs is better than the Skyrme 
EDFs obtained with the same optimization protocol. The relatively large 
variations among the 8 considered parametrizations is very encouraging as it 
suggests that even without “beyond mean-field effects” such as zero-point 
correlations energies, particle-vibration couplings, etc., observables such as 
binding energies are sensitive to the details of the EDF. It is 
remarkable that, on average, the quality of the prediction increases noticeably 
as one goes further up in the chiral expansion. The exception is the 
N2LOΔ+3N EDF, for which predictions of binding 
energies and even the sp structure degrade relatively to other 
EDFs. However, one should keep in mind that, from a statistical 
perspective, calibrating functionals at the single-reference level implies that 
the model has “defects”, i.e., it is not designed to accurately reproduce 
specific observables; see discussion on section <ref>. In practice, 
trying to fit both binding energies in closed shell nuclei like ^208Pb and 
well-deformed nuclei in the rare-earth region could lead to overfitting issues. 

In this work, we have left out the estimate of uncertainties – only 
quantifying numerical errors induced by approximating Yukawa form factors by a 
sum of Gaussian functions and by interpolating coupling functions. In 
particular, it could be worth studying in more details the exchange 
contribution to the energy by (i) calibrating a Hartree-only functional where 
the exchange contribution of the long-range chiral potential would be dropped 
entirely, and (ii) conversely, calibrating a functional where the Fock 
contribution from the chiral potential would be computed exactly by expanding 
it onto a sum of Gaussians, like the Hartree term. If we restrict ourselves to 
NN potentials only, the computational effort is not significantly larger.

Since the DME EDFs originate from NN and 3N potentials 
from chiral perturbation theory and their coupling constants have been 
determined with the exact same optimization protocol, these EDF lend 
themselves particularly well to studies of systematic uncertainties. Together 
with the well-established machinery to quantify statistical uncertainties with 
either covariance or Bayesian methods, such studies could shed more light on 
the true predictive power of these EDFs. Since the DME functionals 
seem to encode some effects traditionally associated with beyond mean-field 
physics, it would also be natural to explore fits at the MR-EDF level.







Support for this work was partly provided through the Scientific
Discovery through Advanced Computing (SciDAC) program funded by
U.S. Department of Energy, Office of Science, Advanced Scientific
Computing Research and Nuclear Physics. It was partly performed under
the auspices of the US Department of Energy by the Lawrence Livermore
National Laboratory under Contract DE-AC52-07NA27344. This work was
supported by the US Department of Energy under grant number
DE-FG02-93ER-40756 and the National Science Foundation under
grant number PHY–1614460. Computing support for this work came from the
Lawrence Livermore National Laboratory (LLNL) Institutional Computing
Grand Challenge program.

 




