Ab Initio study of neutron drops with chiral Hamiltonians

We report ab initio calculations for neutron drops in a 10 MeV external harmonic-oscillator trap using chiral nucleon-nucleon plus three-nucleon interactions. We present total binding energies, internal energies, radii and odd-even energy differences for neutron numbers N = 2 - 18 using the no-core shell model with and without importance truncation. Furthermore, we present total binding energies for N = 8, 16, 20, 28, 40, 50 obtained in a coupled-cluster approach. Comparisons with Green's Function Monte Carlo results, where available, using Argonne v8' with three-nucleon interactions reveal important dependences on the chosen Hamiltonian.


Introduction
There is significant interest in ab initio solutions for systems of neutron drops trapped in external potentials aimed at providing insights into properties of unstable neutron-rich nuclei and neutron star matter [1,2,3]. At the same time, comparisons of neutron drop results using different microscopic interactions provide information on the isovector part of the nucleon-nucleon (NN) interaction and the T = 3/2 component of the threenucleon (3N) interaction. With these goals in mind, we present the first ab initio results for pure neutron systems using chiral NN+3N Hamiltonians in an external trap and compare with results previously obtained using high-precision phenomenological NN+3N Hamiltonians.
At the same time, significant theoretical advances regarding the underlying Hamiltonians, constructed within chiral effective field theory (EFT), provide a firm foundation for nuclear many-body calculations rooted in QCD [21,22] leading us to adopt a chiral EFT Hamiltonian here. We also make use of the similarity renormalization group (SRG) approach [23,24,25,26,27,28] that provides a straightforward and flexible framework for consistently evolving (softening) the Hamiltonian and other operators, including 3N interactions [11,29,30].
The goal of this paper is twofold. First, we aim to provide results for neutron drop systems in a 10 MeV harmonic-oscillator (HO) trap using realistic chiral NN+3N interactions with uncertainty estimates where feasible. Second, we present comparisons between our results and those of other high-quality NN+3N interactions. In particular, we compare with results obtained using Green's Function Monte Carlo (GFMC) [2,3] where the Argonne v 8 (AV8 ) NN interaction [31] was used in conjunction with the Urbana IX (UIX) 3N interaction [31] and with the Illinois-7 (IL7) 3N interaction [32].
We limit our investigations to a single form of the chiral NN+3N interaction. That is, we use the chiral NN interaction at N 3 LO with 500 MeV/c cutoff from Ref. [33] together with the 3N interaction at N 2 LO [34] in the local form of Ref. [35] for 500 MeV/c cutoff with low-energy constants determined entirely in the threenucleon sector [36]. This is also the Hamiltonian used in Refs. [11,29,30,37,38,39] for ab initio studies of nuclear properties. We evolve this Hamiltonian using the free-space SRG to two representative flow parameters or momentum scales to examine the scale-dependence of our results. As in the earlier applications, we retain the induced many-body interaction up to the threenucleon level and neglect induced four-nucleon (and beyond) interactions. Depending on whether the initial chiral 3N interaction is included or not we use the term NN+3N-full or NN+3N-induced, respectively, to characterize the SRG-evolved Hamiltonian.
In selected cases, we also compare our results with those obtained using the JISP16 NN interaction which lacks a 3N contribution [40,41]. The neutron drop results with JISP16 have also appeared previously in Ref. [3].
In Section 2, we briefly review the formalism and summarize related results from previous work. The results for our neutron drop observables are presented in Section 3. Section 4 summarizes our conclusions and provides perspectives on future efforts.

Theoretical Framework
We employ ab initio configuration interaction and coupled-cluster methods to solve for the properties of neutron drops. In the first approach, the no-core shell model (NCSM), we follow Refs. [4,5,6] where, for a given NN and 3N interaction we diagonalize the resulting many-body Hamiltonian in a sequence of truncated HO basis spaces. The basis spaces are characterized by two parameters: N max specifies the maximum number of total HO quanta above the lowest allowed HO Slater determinant and ω specifies the HO energy of the basis. This latter variable is distinct from the HO energy of the trap which is fixed to be 10 MeV in the present application. The goal is to achieve convergence as indicated by independence from these two basis parameters.
We also employ an extension of the NCSM, the importance truncated NCSM (IT-NCSM), for which we follow Refs. [10,11,12,30] where subspaces of the N max -truncated spaces are dynamically selected according to an importance measure derived from perturbation theory. The IT-NCSM uses this importance measure κ ν for the individual many-body basis states and retains only states with |κ ν | above a threshold κ min in the model space. Through a variation of this threshold and an a posteriori extrapolation κ min → 0 the contribution of discarded states is recovered. We use the sequential update scheme discussed in Refs. [10,30], which connects to the full NCSM model space and thus the exact NCSM results in the limit of vanishing threshold. We recently compared the IT-NCSM with the NCSM in basis spaces where calculations in both approaches are feasible [39].
For neutron drops that exhibit subshell closure we apply the coupled-cluster (CC) method which is capable of providing results for heavier systems. We use single-reference CC with singles and doubles excitations [42], in which the ground state |Ψ of a many-body Hamiltonian is parametrized by the exponential ansatz |Ψ = e T 1 +T 2 |Φ , where T n are n-particle-n-hole excitation operators acting on a single Slater-determinant reference state |Φ , which is the Hartree-Fock determinant in our calculations. Effects of the T 3 clusters are included through an a posteriori correction to the energy via the ΛCCSD(T) [20,43] method. The underlying single-particle basis is a HO basis truncated in the principal oscillator quantum number e = 2n + l ≤ e max . In this work we quote results from e max = 12 model spaces that are sufficiently well converged. Including explicit 3N interactions into CC calculations results in a significant increase of the computational expense [44,45,46]. To facilitate the calculations, we use the normal-ordered two-body approximation (NO2B) [44,45] to the 3N interaction, which was shown to be very accurate in calculations of atomic nuclei [44,45,46]. In the NO2B approximation, contributions of the 3N interaction are demoted to lower particle ranks through normal-ordering techniques, and the residual normal-ordered three-body operator is discarded [44,45]. Due to their enormous number, not all of the 3N matrix elements that would be required by the large model spaces employed by the CC method can be included in the normal-ordering procedure. For that reason we impose an energy truncation e 1 +e 2 +e 3 ≤ E 3 max = 14 on the 3N matrix elements. We have checked that our results are converged with respect to this truncation.
We adopt the chiral NN+3N interaction described above as this Hamiltonian has been applied in a range of ab initio calculations of light and medium-mass nuclei [11,29,30,37,38,39,47]. For a detailed discussion of the SRG evolution in the 3N sector adopted here, see Ref. [30]. Unless otherwise specified, we employ SRG-evolved interactions with α = 0.08 fm 4 corresponding to a momentum scale λ SRG = α −1/4 = 1.88 fm −1 .
In our (IT-)NCSM calculations, the size of the largest feasible model space is constrained by the total number of required 3N interaction matrix elements as well as by the number of many-body matrix elements that are com-puted and stored for the iterative Lanczos diagonalization algorithm. Through an efficient JT -coupled storage scheme and an on-the-fly decoupling during the calculation of the many-body Hamilton matrix [11,30,48], the limit arising from handling the 3N matrix elements has been pushed to significantly larger model spaces.
For the full NCSM calculations we employ the MFDn code [49,50,51,52,53] that is highly optimized for parallel computing. In order to exploit parallel architectures that include GPUs, such as the Titan facility at Oak Ridge National Laboratory, we have developed and implemented new algorithms that significantly speed up the required decoupling transformations [48]. The IT-NCSM calculations are performed with a dedicated code that has been developed to accommodate the specific demands of an importance-truncated calculation in a framework optimized for parallel performance. Due to the reduction of the model-space dimension resulting from the importance truncation, typically by two orders of magnitude, the many-neutron Hamiltonian matrix is significantly smaller and the memory needs are drastically reduced. The CC calculations are performed based on an highly-efficient angular-momentum coupled implementation [19]. The most time-consuming part, the calculation of the ΛCCSD(T) energy correction, is embarrassingly parallel and therefore exhibits a nearly perfect scaling. The evaluation of the ΛCCSD(T) energy correction of a N = 50 neutron drop can be performed within a few thousand CPU hours.

Results
We begin with a demonstration of the ground-state energy convergence with increasing basis space truncation N max and with a range of HO basis ω values as shown in Fig. 1 for N = 10 neutrons using the NCSM and the IT-NCSM approaches. We observe an excellent agreement of the NCSM and the IT-NCSM energies wherever both are available. As guaranteed by the variational principle, the results converge uniformly from above for all values of ω. The convergence is fastest for basis ω values slightly above the HO trap strength of 10 MeV as may be expected. Since we find the convergence pattern for all other neutron numbers very similar to that shown in Fig. 1, we will quote the lowest energy for fixed neutron number at the largest N max obtained as our final result for the total ground state energy. From the convergence pattern in Fig. 1 we deduce that it is reasonable to take our uncertainty as the difference between the quoted energy and the energy at the next smaller value of N max at the same ω value. NCSM IT-NCSM 3N-full 3N-ind. 3N-full 3N-ind. N J π α = 0.04 α = 0.08 α = 0.08 α = 0.08 α = 0.08 2 0 + 23.90(1) 23.900 (5) (7) 283.47 (7) Table 2: Total ground-state energies E tot and correlation energies ∆E in units of MeV at closed subshells obtained with coupled-cluster theory at the ΛCCSD(T) level for a 10 MeV HO trap. We use the SRGevolved NN+3N-full and NN+3N-induced Hamiltonian at SRG evolution scales α = 0.04 and 0.08 fm 4 as in Table 1. We use Hartree-Fock reference states obtained in an HO single-particle basis with e max = 12 and ω = 16 MeV. Uncertainties are quoted in parenthesis for the last significant figure and discussed in the text. We then portray our results for the total energy in Fig. 2 scaled by the Thomas-Fermi N-dependence (N 4/3 ) and by the HO well strength following the practice of Refs. [1,2,3]. Remarkably, our results with chiral Hamiltonians are rather insensitive to the presence or absence of the initial 3N interaction indicating that the contribution of the chiral 3N interaction is very small in the T = 3/2 channel. This contrasts the results with AV8 plus 3N interactions (UIX and IL7) from Refs. [2,3] in Fig. 2. To highlight the difference in the sensitivity to the 3N interaction, we provide an inset in Fig. 2 depicting the ratio of the results with the initial NN+3N interaction to those with the initial NN interaction alone. For the chiral interactions, this ratio fluctuates at a level of less than 1% over this range of N. However, for the other interactions, it is either a steadily increasing (AV8 +UIX) or decreasing (AV8 +IL7) function of N.
We find a pronounced dip in the total energy for the chiral interactions due to expected shell closure at N = 8 that is similar to dips seen in the results with the other interactions. For N ≥ 10 we observe an apparent coincidence where the chiral results follow closely those of AV8 without 3N interactions.
We tabulate ground-state energies (unscaled) in Table 1 for two values of the SRG evolution parameter α and compare the NCSM and IT-NCSM results. These results show only a weak dependence on the SRG evo- The results for the chiral interactions are obtained in the IT-NCSM or NCSM with the largest accessible N max (cf. Table 1) and for closed subshells in ΛCCSD(T) (cf. Table 2). The results labelled AV8 , AV8 +UIX and AV8 +IL7 are adopted from Ref. [3]. The inset shows the ratio of the total energies obtained with the initial 3N interaction and without for the largest N max . Characteristic uncertainties are shown in the inset for the largest N only. We discuss the uncertainties for our results in the text.
lution scale and are consistent between the NCSM and IT-NCSM to within quoted uncertainties. Note that the IT-NCSM energies are always obtained for larger N max and are thus systematically below the NCSM results, as expected from the variational principle. The uncertainties quoted in the table are estimated from the difference of the total energies in the two largest model spaces.
In the case of the IT-NCSM the additional uncertainties from the threshold extrapolation are typically an order of magnitude smaller than the quoted values.
In Table 2 we summarize the total ground-state energies as well as the correlation energies beyond Hartree-Fock obtained from ΛCCSD(T) calculations for neutrons with closed subshells up to N = 50. The quoted uncertainties are again estimated from the energy difference in the largest two model spaces. In addition they include the effects of the E 3 max -truncation of the 3N matrix elements and of contributions beyond the triply excited clusters in a very conservative estimate. We observe an excellent quantitative agreement of ΛCCSD(T), IT-NCSM, and NCSM results well within the estimated uncertainties.
As mentioned above, all of these results show that the inclusion of the chiral 3N interaction affects the total energies by less than 1%, even in the heaviest neutron drops considered here. Furthermore, the dependence   on the SRG flow parameter, indicative for the contribution of SRG-induced multi-nucleon interactions beyond the 3N level, is extremely small. Note that the NN+3N-full Hamiltonian adopted here leads to a significant flow-parameter dependence in medium-mass nuclei [45]. This indicates that the main origin of the flowparameter dependence in symmetric nuclear systems are the T = 1/2 components of the chiral 3N interaction. We present the single-particle root-mean-square  (rms) radii in Fig. 3 for a selection of Hamiltonians, including the nonlocal NN interaction JISP16. Our results for the chiral Hamiltonian show a smooth trend with increasing N in contrast to the other Hamiltonians that yield significant odd-even effects below N = 8. Above N = 9 our rms radii follow closely those obtained with JISP16, whereas the radii obtained with AV8 +UIX are significantly larger. It is noteworthy that the chiral results for the rms radii are again rather insensitive to the presence or absence of the chiral 3N interaction.
The energy differences between neighboring systems in N are shown in Fig. 4 where, without interactions, we expect the differences to be simple multiples of the HO trap energy as indicated by the solid horizontal lines in Fig. 4. The interactions produce the strong odd-even effect conventionally characterized as "pairing energy". Within the indicated uncertainties, the pairing effects are in rough agreement between the different interactions up to N = 8. For 11 ≤ N ≤ 14, our chiral results are in agreement with the JISP16 results (within error estimates), but the GFMC results for the AV8 +UIX interaction indicate less pairing.
We portray the internal energies (scaled) in Fig. 5 which are defined as the total energies of Fig. 1 less the expectation value of the HO trap potential computed with the ground-state wavefunction. Here one observes the odd-even staggering due to pairing and, again, the lack of significant sensitivity of the chiral results to the inclusion of the 3N interaction. We note that the role of pairing appears to be more pronounced in the results with the chiral interactions than with the other interactions depicted in Fig. 5.

Summary and Conclusions
We have presented NCSM, IT-NCSM, and CC results for neutron drops in a 10 MeV external HO trap using chiral NN+3N interactions. We examined total binding energies, odd-even energy differences, internal energies, and radii. By comparing with GFMC results using AV8 plus 3N interactions we found significant dependences on the selected Hamiltonian which should have an impact on phenomenological energy-density functionals that may be derived from them. Furthermore, we found surprisingly weak contributions in these neutrondrop observables from the inclusion the chiral 3N interaction. Based on systematic trends shown in previous neutron-drop investigations [1,2,3], we anticipate these conclusions will persist over a range of HO well strengths from 5 MeV to 20 MeV and will follow the trends seen here for larger values of the neutron number N as well.
Furthermore, as more than one many-body method is applied to the same problem, we establish that our NCSM, IT-NCSM and CC results are consistent with each other within the quoted uncertainties. We find little sensitivity of these results to the SRG evolution of the chiral interactions over the range of the SRG flow parameter we investigated.
Next steps will include the study of excited states and of the sensitivity of our observables to the underlying chiral Hamiltonian. Furthermore, the comparison of our ab initio results to predictions from state-of-the-art energy-density functionals will help constrain the latter in regimes of extreme isospin.