1 Introduction

Ab initio calculations of nuclear structure [1-9] face the challenge of describing a complex, multiscale quantum many-body system. The goal is to directly solve the many-body problem for a system of protons and neutrons, with realistic internucleon interactions [10-13]. However, the nucleus is governed by a strong, short-range interaction. Short-range correlations, tightly-bound α clusters [14], and weakly-bound halo nucleons [15, 16] introduce dynamics over differing length scales and energy scales, within the same nucleus, which must be simultaneously described within the same many-body calculation.

Natural orbitals [17-22] provide a means of adapting the single-particle basis to better match the physical structure of the many-body wave function. Natural orbitals are obtained by diagonalizing the one-body density matrix, deduced from a preliminary many-body calculation using an initial reference single-particle basis. The Laguerre function basis [23, 24] has been used as the starting point for natural orbitals in atomic electron-structure calculations [18], while we start from the harmonic oscillator orbitals [25] more familiar to the nuclear structure context [26]. The natural orbital basis builds in important contributions from high-lying orbitals of the initial basis — for the present application, high-lying oscillator shells — thereby accelerating the convergence of wave functions, energies, and other observables.

In this letter, we present a framework for ab initio no-core configuration interaction (NCCI) [8] calculations with a natural orbital basis and demonstrate improved convergence for the lightest neutron halo nucleus ⁶He [15]. When used with recently-proposed infrared (IR) basis-extrapolation schemes [27, 28], we show that natural orbitals provide improved independence of basis parameters for predictions of energy and radius observables.

2 Natural orbitals

In NCCI calculations, the nuclear many-body Schrödinger equation is formulated as a matrix eigenproblem, where the Hamiltonian is represented within a basis of Slater determinants, i.e., antisymmetrized products of single-particle states. Conventionally, harmonic oscillator orbitals [25] are used, and the basis is truncated to a maximum allowed number N_max of oscillator excitations [8]. The calculated wave functions, energies, and observables depend upon both the truncation N_max and the oscillator length b of the basis (or, equivalently, the oscillator energy ħw∝ b^-2). The solution of the full, untruncated many-body problem could, in principle, be obtained to any desired accuracy, by retaining a sufficiently complete basis set. However, the dimension of the NCCI problem increases rapidly with the number of nucleons and included single-particle excitations, as shown in Fig. 1. Currently available computational resources therefore limit the convergence of calculated states and observables [29-31].

Figure 1:Full-screen View

Growth of the NCCI problem dimension as a function of the number of oscillator excitations N_max included in the basis, for selected nuclides, including ⁶He (red curve). The dimensions shown are for spaces with zero angular momentum projection (M=0) and positive parity.

We therefore seek a physically-adapted basis, in which the nuclear many-body wave function can be efficiently and accurately described, subject to the constraint of accessible problem dimensions. The natural orbital basis minimizes the mean occupation of states above the Fermi surface [21, 32], thus reducing the contribution of high-lying orbitals in describing the many-body wave function. Intuitively, the natural orbitals may be understood as representing an attempt to recover the basis in which the many-body wave function most resembles a single Slater determinant. Although we cannot expect to reduce the complex, highly-correlated nuclear wave function to a single Slater determinant, as assumed in the Hartree-Fock approximation, we may expect that transforming to a more natural single-particle basis could enhance the role of a comparatively small set of dominant Slater determinants, and thereby accelerate convergence within a Slater determinant expansion.

An NCCI state of good total angular momentum can only, in general, be obtained as a superposition of several Slater determinants of nljm single particle states (here, n is the radial quantum number, l labels the orbital angular momentum, j labels the resultant angular momentum after coupling to the spin, and m labels its projection). Therefore, we cannot, in general, expect the nuclear eigenfunctions to resemble a single Slater determinant. Rather, we hope to recover wave functions which most resemble a single configuration of nucleons over nlj orbitals. For this purpose, we consider the scalar densities ${\rho }_{ab}^{(0)}\equiv \langle \Psi {[c_b^{†}{\tilde{c}}_a]}_{00}|\Psi \rangle$, where $c_a^{†}$ represents the creation operator for a nucleon in orbital a=(nalaja) and the brackets [...]₀₀ represent spherical tensor coupling to zero angular momentum. For a single configuration |Ψ〉, the scalar densities are diagonal, and the diagonal entries give the occupations of the contributing orbitals [${\mathcal{N}}_a=(2j_a+{1)}^{1/2}{\rho }_{aa}^{(0)}$] [33]. Otherwise, for the general case of a many-body state |Ψ〉, natural orbitals are defined by diagonalizing this scalar density matrix. The scalar density matrix only connects orbitals of the same l and j, i.e., differing at most in their radial quantum number n, so the transformation to natural orbitals induces a change of basis on the radial functions separately within each lj space [$|n^{\prime }lj\rangle ={\displaystyle {\sum }_n}{a}_{n^{\prime }\mathrm{,}n}^{(lj)}|nlj\rangle$].

An initial NCCI calculation is carried out in the oscillator basis. This provides a scalar density matrix for the ground-state wave function, which is then diagonalized to yield the natural orbitals. The eigenvalue associated with a natural orbital represents its mean occupation in the many-body wave function. We order the natural orbitals by decreasing eigenvalue of the density matrix [21], i.e., starting with n=0 for the natural orbital with highest eigenvalue [〈𝒪_0lj〉≥〈𝒪_1lj〉≥...], thereby providing an n quantum number for an N_max-type truncation scheme (see Ref. [34]).

The lowest p_3/2 natural orbital obtained from the ⁶He ground-state one-body densities is illustrated in Fig. 2, taking an example from the NCCI calculations presented below. In a traditional shell-model description, ⁶He consists of two protons and two neutrons in a filled s shell ($N\equiv 2n+l=0$), plus two neutrons in the valence p shell (N=1). In the oscillator-basis NCCI calculations, this general structure is reflected in a nearly-filled s shell. The next most heavily occupied orbital is then the neutron 0_p3/2. We observe that the corresponding natural orbital [Fig. 2] receives extended contributions from high-lying oscillator shells and thus acquires a substantial large-r tail compared to the oscillator orbital, as may be expected for a weakly-bound halo nucleon.

Figure 2:Full-screen View

Radial wave function for the neutron 0_p3/2 natural orbital (heavy curve) derived from the ⁶He ground state calculation in the harmonic oscillator basis, along with the contributions from individual oscillator basis functions (gray curves). The squared amplitudes P(N) of these contributions are shown in the inset. The initial oscillator basis for this calculation has N_max=16 and ħw=20 MeV.

The NCCI calculation using the natural orbital basis is no more computationally difficult than the original oscillator basis calculation. It is simply necessary, in preparation, to carry out a similarity transformation of the input Hamiltonian, as described in Refs. [35-36].

3 Results

Several experimental properties of the ground state of ⁶He support the interpretation that it consists of a weakly-bound two-neutron halo surrounding a tightly-bound α core [15, 16]. The two-neutron separation energy for ⁶He is only 0.97 MeV, out of a total binding energy of 29.27  MeV [37]. Experimentally, the onset of halo structure along the He isotopic chain is indicated by a jump in the measured charge and matter radii, from ⁴He to ⁶He. The root mean square (RMS) point-proton distribution radius rp, which may be deduced [38] from the measured charge radius, increases by ∼32% from ⁴He [rp=1.462(6)  fm ] to ⁶He [rp=1.934(9) fm] [39-41]. This increase may be understood as a consequence of halo structure, arising from the recoil of the charged α core against the halo neutrons (as well as possible contributions from swelling of the α core [41]).

The initial oscillator basis NCCI calculations from which we derive natural orbitals for ⁶He cover a range of oscillator basis parameters ħw=10 MeV to 40 MeV with truncations N_max≤16, as considered in Ref. [36]. The NCCI calculations are carried out using the code MFDn [42, 43], with the JISP16 two-body internucleon interaction [12] plus Coulomb interaction.

The ground state energy eigenvalues obtained for ⁶He using the harmonic oscillator (dashed lines) and natural orbital (solid lines) bases are compared in Fig. 3(a). The energies from the natural orbital calculations are lower (thus, by the variational principle, closer to the true value) than those from the harmonic oscillator calculations and are also less dependent upon ħw. The improvement in convergence afforded by the natural orbitals ranges from approximately one step in N_max in the vicinity of the variational minimum (ħw≈15–20 MeV) to several steps in N_max towards the ends of the calculated ħw range [44].

Figure 3:Full-screen View

Comparison of ⁶He ground state properties calculated using harmonic oscillator (HO, dashed lines) and natural orbital (NO, solid lines) bases: (a) energy and (b) ground state proton radii. These are shown as functions of the oscillator basis ħw, for N_max=10 to 16 (as labeled).

For the ⁶He proton radius, shown in Fig. 3(b), the natural orbital basis NCCI calculations lead the harmonic oscillator basis calculations in convergence by more than one step in N_max at ħw≈20 MeV and by several steps in N_max at the high end of the ħw range. These radii obtained from the natural orbital calculations are also less dependent upon ħw than those obtained from the harmonic oscillator calculations.

The goal we set out to achieve is to find the true results for observables as they would be obtained in the full, infinite-dimensional space. Although full convergence is not achieved, even with the natural orbitals, the improved convergence motivates us to attempt to obtain estimates of the converged results via basis extrapolation methods [27, 29, 45, 46]. Infrared extrapolation schemes [27, 28, 47-49] are based on the premise that the solution of the many-body problem in a truncated space effectively imposes infrared (long-range) and ultraviolet (short-range) cutoffs. For bases with high enough ħw (and N_max) ultraviolet convergence is assumed, and any remaining incomplete convergence is attributed to the failure of the basis to reproduce the long-range tail of the many-body wave function.

A basis consisting of harmonic oscillator orbitals with no more than N quanta cannot fully resolve long-range physics beyond the classical turning point L(N,ħw)=[2(N+3/2)]^1/2 b(ħw) [27],³ where b(ħw)=(ħc)/[(mNc²)(ħw)]^1/2 is again the oscillator length, with mN the nucleon mass. The calculated energy and observables are expected to depend only on this cutoff L, approaching the true converged values as L→∞. For energy eigenvalues, it is expected that [47, 48]

where E_∞, a₀, and k∞ are to be deduced as fitting parameters from the results of calculations in truncated spaces. Taking →∞, we extract E_∞ as an estimate for the true energy. For mean square radii, it is expected that, letting $\beta \equiv 2k_{\infty }L$,

for β≫1, where r∞, c₀, and c₁ are similarly deduced from calculations in truncated spaces, and r∞ provides an estimate of the true RMS radius.

The extrapolated values for the ⁶He ground state energy and proton radius are shown in Fig. 4. We restrict ourselves to a straightforward application of (1) and (2), based on three-point extrapolation in N_max at fixed ħw. Calculations at low ħw may not provide the assumed ultraviolet convergence, while poor infrared convergence at high ħw leads to an excessively large correction and thus poor extrapolation.

Figure 4:Full-screen View

Infrared basis extrapolations for the ⁶He ground state energy (top) and point proton radius (bottom), based on calculations in the harmonic oscillator basis (left) and natural orbital basis (right). The extrapolations (diamonds) are shown along with the underlying calculated results (plain lines) as functions of ħw at fixed N_max (as indicated). Experimental values (circles) are shown with uncertainties. The shaded bands reflect the mean values and standard deviations of the extrapolated results, at the highest N_max, over the ħw range considered.

The extrapolated ⁶He ground state energies from the natural orbital NCCI calculations [Fig. 4(b)] are considerably less ħw-dependent than the extrapolated energies from the harmonic oscillator NCCI calculations [Fig. 4(a)]. The extrapolations obtained for different N_max are also considerably more consistent (in the figure, N_max refers to the highest N_max in the three-point extrapolation). The extrapolated ground state energies obtained with the harmonic oscillator and natural orbital bases at ħw=20 MeV (chosen close to the variational energy minimum) and N_max=16 are consistent with each other to within their respective variations, giving E≈-28.79 MeV and E≈-28.80 MeV, respectively.

Once the many-body calculation is under control, any remaining deviation of calculated values from nature may be attributed to deficiencies in the internucleon interaction. Comparing to the experimental binding energy of 29.27 MeV thus indicates that the JISP16 interaction underbinds ⁶He by ∼0.5 MeV⁴. (For comparison, the binding of ⁴He obtained with JISP16 matches experiment to within ∼0.003 MeV [29].)

The extrapolated proton radii extracted from the NCCI calculations with the natural orbital basis [Fig. 4(d)] similarly demonstrate a reduced ħw dependence and N_max dependence, as compared to the extrapolations from the oscillator-basis calculations [Fig. 4(c)]. At the highest calculated N_max (N_max=16), the extrapolated rp varies only by ∼0.02 fm across the range of ħw values shown (ħw≈14 MeV to 40 MeV), and the N_max dependence is comparable. We must emphasize that the variations in extrapolated values at best provide a rough guide to how well we can trust these extrapolated values as reflecting the true radius which would be obtained in an untruncated many-body calculation. Nonetheless, the ħw-independence and N_max-independence of the calculations at the ∼0.02 fm level is reassuring.

Taking the extrapolated proton radius at ħw=20 MeV and N_max=16 as representative gives rp≈1.82 fm.⁵ Thus, it would appear that the ab initio NCCI calculations with the JISP16 interaction, while qualitatively reproducing the increase in proton radius with the onset of halo structure in ⁶He, do yield a quantitative shortfall of ∼0.12 fm (or ∼6%) for the proton radius of ⁶He.

4 Conclusion

Describing the nuclear many-body wave function within truncated spaces is challenging due to the need to describe, simultaneously, long-range asymptotics and short-range correlations. Natural orbitals, obtained here by diagonalizing one-body density matrices from initial NCCI calculations using the harmonic oscillator basis, build in contributions from high-lying oscillator shells, thus accelerating convergence.

In the present application to the halo nucleus ⁶He, improvement is by about one step in N_max near the variational minimum in ħw, and significantly more for other ħw values (Fig. 3). To put these gains in perspective, we note that an increment in N_max results in an increase in matrix dimension of about a factor of 3–5, as seen in Fig. 1, with much larger increase in the computational costs [48]. Although full convergence is still not achieved, the calculations using natural orbitals provide improved basis parameter independence for extrapolations with respect to the infrared cutoff of the basis (Fig. 4).

The successful application of natural orbitals to ab initio nuclear NCCI calculations presented here provide a starting point for exploring ideas (some taken from electron structure theory) which may more fully realize the potential of the natural orbital approach:

(1) NCCI calculations based on natural orbitals yield improved one-body densities which can, in turn, be diagonalized to yield new natural orbitals. Natural orbitals constructed through such an iterative method can rapidly build in additional contributions from high-lying shells, thereby potentially further accelerating convergence [20, 49].

(2) An improved reference basis for the initial NCCI calculation may also boost the convergence of the subsequent natural orbital calculations. For instance, the Laguerre functions, commonly used as the starting point for natural orbitals in electron-structure calculations, also have the correct exponential asymptotics for nucleons bound by a finite-range potential.

(3) The structure of nuclear excited states can vary markedly from that of the ground state. Natural orbitals constructed by diagonalizing the density matrices from excited states, rather than from the ground state, may more effectively accelerate convergence of those excited states [20].

(4) Finally, natural orbitals are conducive to a more efficient many-body truncation scheme than the conventional oscillator N_max scheme. The eigenvalues associated with the natural orbitals, by providing an estimate of the mean occupation of each orbital in the many-body wave function, also suggest a means of estimating the relative importance of Slater determinants involving these orbitals.