Accurate bulk properties of nuclei from $A = 2$ to $\infty$ from potentials with $\Delta$ isobars

We optimize $\Delta$-full nuclear interactions from chiral effective field theory. The low-energy constants of the contact potentials are constrained by two-body scattering phase shifts, and by properties of bound-state of $A=2$ to $4$ nucleon systems and nuclear matter. The pion-nucleon couplings are taken from a Roy-Steiner analysis. The resulting interactions yield accurate binding energies and radii for a range of nuclei from $A=16$ to $A=132$, and provide accurate equations of state for nuclear matter and realistic symmetry energies. Selected excited states are also in agreement with data.

Introduction.-Ideas from chiral effective field theory (EFT) and the renormalization group [1][2][3][4][5], advances in computing power, and developments of many-body methods that scale polynomially with mass number, have propelled ab initio calculations of atomic nuclei from light [6,7] to medium-mass isotopes [8][9][10][11][12][13][14][15][16]. Such approaches are now starting to explore ever-increasing fractions of the nuclear chart [17]-a task once thought to be reserved for computationally less expensive mean-field methods [18][19][20]. Ab initio computations make controlled approximations that give a quantifiable precision in the solution of the quantum many-body problem. All these calculations are only as accurate as the employed nuclear Hamiltonians.
In the past two decades, our understanding of the nuclear interactions has evolved from phenomenological models [21,22] to potentials whose improvement is guided by ideas from EFT [3,5,23]. The quest to link such potentials to quantum chromodynamics, the microscopic theory of the strong nuclear interaction, is ongoing [24][25][26][27]. Recent progress in nuclear potentials from chiral EFT include (i) the identification [28,29] of a redundant term at next-to-next-to-next-to leading order, (ii) the systematic and simultaneous optimization of nucleon-nucleon and three-nucleon potentials [30], (iii) the construction of local potentials for use with quantum Monte Carlo methods [31,32], and (iv) the development of high-order interactions [33][34][35].
In spite of these advances, nuclear potentials from chiral EFT have long struggled to accurately reproduce bulk nuclear properties such as charge radii and binding energies of finite nuclei, and the saturation point and the symmetry energy of infinite nuclear matter. Notable exceptions are the 1.8/2.0(EM) [36] and the N N +3N (lnl) interactions [37] (which accurately reproduce binding energies and spectra of selected nuclei up to tin isotopes [15] but yield too small radii), and the NNLO sat potential [38] (which accurately describes binding energies and radii up to calcium isotopes [39] but is less accurate for spectra). The novel family of interactions [40] also seems promising but is not much explored yet. While it is not clear what distinguishes these particular interactions from their many peers that become inaccurate beyond oxygen isotopes [41], some key ingredients to nuclear binding have been uncovered: Nuclear lattice EFT computations revealed that non-locality is essential for lowmomentum interactions (with momentum cutoffs below about 600 MeV) [42], and that elastic alpha-alpha scattering, for example, is very sensitive to the degree of nonlocality [43]. In our interpretation, these results indicate that the finite size of nucleons plays an important role in nuclear binding [44]. Consistent with this view is the finding that the inclusion of ∆ isobar degrees of freedom, i.e. excited states of the nucleon that reflect its finite size, considerably improve the saturation properties of chiral potentials [45]. This suggests that it might be profitable to include ∆-isobar degrees of freedom and nuclear matter properties into the construction and optimization, respectively, of nuclear potentials [46,47].
In this Letter we report on chiral potentials that accurately describe bulk properties of finite nuclei and nuclear matter. We optimized ∆-full interactions and included nuclear matter properties as calibration data. We used the coupled-cluster (CC) method to compute groundstate energies, charge radii, and spectra of nuclei up to tin, plus the equation of state for nuclear matter.
Optimization.-The ∆-full interactions are based on Refs. [23,[48][49][50][51], and its specific form was used in Ref. [45]. We employed standard nonlocal regulator functions f (p) = exp (p/Λ) 2n that act on relative momenta p, see, e.g., Refs. [52,53]. We constructed three potentials: two of them with Λ = 450 MeV (and power n = 3 in the regulator) at next-to-leading order (NLO) and one at next-to-next-to-leading (NNLO) order with a cutoff Λ = 394 MeV (and power n = 4). The softer interaction has a momentum cutoff and regulator power exactly as the 1.8/2.0(EM) potential. At NNLO there are 17 low-energy coefficients (LECs) that parametrize the interaction. The pion-nucleon LECs c 1,2,3,4 were held fixed during the optimization and taken as the central values from the recent Roy-Steiner analysis [54]-see Table I. The LECs of the nucleon-nucleon and three-nucleon po-tentials were simultaneously constrained by the following data: low-energy nucleon-nucleon scattering data from the Granada phase shift analysis [55] up to 200 MeV scattering energy in the laboratory system; observables of few-nucleon systems (with mass numbers A ≤ 4) as listed in Table II; the saturation energy and density, and the symmetry energy of nuclear matter. The inclusion of the latter is in contrast to the construction of the potential NNLO sat [38] which exhibits deficiencies for neutronrich nuclei and neutron matter. The minimization of the objective function was performed with the algorithm POUNDerS [56]. During this process we periodically calculated selected medium-mass nuclei to guide the optimization. The LECs of Ref. [45] were taken as a starting point for the optimization.
The inclusion of nuclear matter properties into the optimization procedure is not without challenges. Here, we used coupled-cluster calculations [57], which are based on a discrete lattice in momentum space. Nondegenerate reference states are "closed shell" configurations, and we used periodic boundary conditions in position space. Systems of 132 nucleons and 66 neutrons exhibit one of the smallest finite-size corrections for symmetric nuclear matter and neutron matter, respectively [57,58]. Unfortunately, such large particle numbers are numerically too expensive to be used in the optimization. However, systems of 28 nucleons and 14 neutrons exhibit predictable differences (about 10%) from systems consisting of 132 and 66 particles, respectively [57]. This allowed us to use the lower-precision computations with smaller system sizes in the optimization. We checked periodically that our estimates for the finite-size corrections were accurate.
To explore the resolution-scale dependence, we optimized the interaction using two cutoffs, namely Λ = 450 and Λ = 394 MeV at the NNLO level. For Λ = 450 MeV we also optimized an NLO interaction. Table I shows the optimized LECs values for the three interactions of this work. In what follows, we label these "Gothenburg-Oak Ridge" potentials as ∆NLO GO (450), ∆NNLO GO (450) and ∆NNLO GO (394). Most of the LECs of the newly constructed potentials are close to the starting point of Ref. [45], with a few exceptions. In particular, the c D and c E for the short-ranged three-nucleon forces have different signs for ∆NNLO GO (450). Figure 1 shows the phase shifts of the new potentials for a few representative neutron-proton channels and compares them to the Granada partial wave analysis [55]. Overall, the phase shifts of ∆NNLO GO are improved compared to NNLO sat , and they are close to the data for laboratory energies below about 125 MeV. We note that the phase shifts of ∆NLO GO (450) are within the uncertainty estimates expected at this order [34,45,59]. Table II summarizes the results of bound-state observables for light nuclei with A ≤ 4. The theoretical results were obtained with the no-core shell model   Table II are obtained from the computed point-proton radii with standard nucleonsize and relativistic corrections, see, e.g., Ref [38]. The two ∆NNLO GO potentials reproduce the experimental energies within less than 0.5% and other observables within 2%. We note that it was important to include the deuteron quadrupole moment as a calibration datum in the optimization. Note that the Q = 0.27 value we aimed is obtained from a theoretical calculation based on the high-precision meson-exchange N N model CD-Bonn. Figure 2 shows the energy per nucleon in symmetric nuclear matter (top) and pure neutron matter (bottom) as a function of density, using 132 nucleons and 66 neutrons, respectively. The black rectangle indicates the empirical saturation region with E/A = −16 ± 0.5 MeV and ρ = 0.16 ± 0.01 fm −3 [36,63]. The CC calculations were performed in the CCD(T) approximation, i.e. with 2p-2h excitations and perturbative 3p-3h corrections as done in Refs. [45,57]. We find the saturation density ρ 0 = 0.169 fm −3 , the symmetry energy S 0 = 32.0 MeV and its slope L = 65.2 for the ∆NNLO GO (450) potential, and ρ 0 = 0.163 fm −3 , S 0 = 31.5 MeV and L = 58.4 for ∆NNLO GO (394). These nuclear-matter properties are more accurate than those reported in Ref. [45].
Results.-Our CC computations of heavier nuclei started from a spherical Hartree-Fock state built from a model-space consisting of 15 major oscillator shells with frequency ω = 16 MeV. The three-nucleon force had the additional energy cut of E 3max = 16 ω. Our calculations employed the coupled-cluster singles, doubles, and leading triples approximation CCSDT-1 [64]. No further truncations were imposed on the three-particle-three-hole amplitudes. This computational achievement was enabled by the Nuclear Tensor Contraction Library  Table III shows the binding energies of selected closedshell nuclei up to 132 Sn. We found that ∆NNLO GO (394) converges faster, especially for heavier nuclei, than ∆NNLO GO (450). This is expected due to its lower momentum cutoff. A dash indicates that the employed model space was too small to achieve reasonably converged energies. Uncertainties from the coupled-cluster method (about 30% of the difference between doubles and triples energies) and from the model space are given in subsequent parenthesis, respectively. For lighter nuclei, the uncertainties from the model-space are omitted because they are much smaller than those of the method. The model-space uncertainties combine the truncation of single-particle model space and the employed E 3max = 16 cut of the three-nucleon interaction. Reference [15] found that for the 1.8/2.0(EM) interaction the binding energy of 100 Sn changes by less than 1% by increasing E 3max from E 3max = 16 to E 3max = 18. This finding guided our estimated model-space uncertainty as the 1.8/2.0(EM) has identical three-nucleon regulator and momentum cut- off as the ∆NNLO GO (394) potential. It is non-trivial to estimate the EFT truncation errors for bound nuclear states since the relevant momentum scale is unknown and the lack of a spin-orbit (LS) force at leading order (LO) give energy degeneracies that hamper CC calculations of non-LS-closed nuclei. Nevertheless, based on the observed order-by-order convergence in Ref. [45], we estimate the EFT truncation errors for the ∆NNLO GO interactions to 1 MeV and 7 MeV in the ground state energies of oxygen and calcium isotopes respectively. We also expect the truncation error to be the dominating source of uncertainty in heavier nuclei. For symmetric nuclear matter and pure neutron matter we expect a truncation error of ±1 MeV and ±2 MeV per nucleon, respectively, with ∆NNLO GO , see Ref. [45] for details. Figure 3 shows the ground-state energies of selected calcium isotopes and compares them to the 1.8/2.0(EM) interaction [66] and to data. The lower borders of the bands are results from CCSDT-1 while the upper borders are from Λ-CCSD(T) [67] which treats triples excitations perturbatively. For 53,55 Ca we employed the particleattached equation-of-motion coupled-cluster (EOM-CC) method with perturbative three-particle-two-hole excitations from Ref. [15], while for 56 Ca we employed the twoparticle attached EOM-CC method from Refs. [68,69]. All interactions accurately describe isotopes from 40 Ca to 56 Ca. We note that the ground-state of 60 Ca is bound by about 10 MeV with respect to 54 Ca, consistent with recent data [70][71][72]. Figure 4 shows that the interactions of this work yield significantly larger charge radii than the 1.8/2.0(EM) potential. Nevertheless, the ∆NNLO GO potentials still fail to explain the unusually large charge radii of neutronrich calcium isotopes. For speculations about the origin of these large charge radii in calcium isotopes we refer the reader to Ref. [73].
We employed EOM-CC methods with singles, doubles, and leading-order triples excitations, EOM-CCSDT-1 [74], to compute excited states of 16,22,24 O, and 48 Ca. The EOM-CCSDT-1 method is computationally demanding, and we therefore limited the number of three-  particle-three-hole excitations by employing the energy cutẼ pqr =ẽ p +ẽ q +ẽ r <Ẽ 3max , whereẽ p = |N p − N F | is the energy difference of the single-particle energies with respect to the Fermi surface N F . This energy cut improved the convergence of the EOM-CCSDT-1 method with respect to the number of three-particle-three-hole excitations [75,76]. In this work we also employed the method developed in Ref. [76] to correct perturbatively for three-particle-three-hole excitations abovẽ E 3max . Here we employed the energy cutẼ 3max = 6 which was sufficient to converge all excited states to within approximately 100 keV. The results are summarized in Table IV   flects that the charge radius is well reproduced in this nucleus (see Ref. [38] for a more detailed discussion on this point). The uncertainties are estimated based on Refs. [76,77]. Here excited states were computed in different truncations within the EOM-CC approach, and it was found that the triples correction to the excited states were about 20% of the EOM-CCSD correlation energy. Using this, we give a conservative error estimate in the first parenthesis which amounts to 6% and 15% of the total excitation energies for the interactions with cutoffs 394 MeV and 450 MeV, respectively. The uncertainties from the truncated model space are given in the second parenthesis.
Summary.-We developed chiral interactions with ∆ degrees of freedom by calibrating LECs to reproduce nucleon-nucleon scattering phase shifts, bound-state observables of few-nucleon systems, and properties of infinite nuclear matter. The resulting ∆NNLO GO potentials yield accurate (within about 2%) binding energies of nuclei up to mass numbers A = 132, and improved radii for medium-mass nuclei. The description of neutron-rich calcium isotopes is improved by including the symmetry energy in the optimization. Selected excited states are also accurately reproduced. This shows that key nuclear properties can be obtained by chiral interactions at nextto-next-to-leading order.
We thank Ragnar Stroberg for fruitful discussions and for providing us with the 1.