Natural orbitals for many-body expansion methods
NNatural orbitals for many-body expansion methods
J. Hoppe,
1, 2, ∗∗ A. Tichai,
3, 1, 2, †† M. Heinz,
1, 2, ‡‡ K. Hebeler,
1, 2, §§ and A. Schwenk
1, 2, 3, ¶¶ Technische Universität Darmstadt, Department of Physics, 64289 Darmstadt, Germany ExtreMe Matter Institute EMMI, GSI Helmholtzzentrum für Schwerionenforschung GmbH, 64291 Darmstadt, Germany Max-Planck-Institut für Kernphysik, Saupfercheckweg 1, 69117 Heidelberg, Germany
The nuclear many-body problem for medium-mass systems is commonly addressed using wave-function ex-pansion methods that build upon a second-quantized representation of many-body operators with respect to achosen computational basis. While various options for the computational basis are available, perturbativelyconstructed natural orbitals recently have been shown to lead to significant improvement in many-body ap-plications yielding faster model-space convergence and lower sensitivity to basis set parameters in large-scaleno-core shell model diagonalizations. This work provides a detailed comparison of single-particle basis sets anda systematic benchmark of natural orbitals in non-perturbative many-body calculations using the in-mediumsimilarity renormalization group approach. As a key outcome we find that the construction of natural orbitalsin a large single-particle basis enables for performing the many-body calculation in a reduced space of muchlower dimension, thus o ff ering significant computational savings in practice that help extend the reach of abinitio methods towards heavier masses and higher accuracy. I. INTRODUCTION
Nuclear many-body theory has witnessed major develop-ments over the last two decades, extending the reach of the ab initio solution of the stationary Schrödinger equation overa wide range of mass numbers in the nuclear chart, cover-ing closed- and open-shell nuclei and including exotic nu-clei [11, 22]. This progress is mainly based on i) the constructionof improved nucleon-nucleon (NN) and three-nucleon (3N)interactions based on chiral e ff ective field theory (EFT) [33–1414] and ii) the extension of many-body theories applicable tomedium-mass nuclei [1515–1818]. The advances of many-bodycalculations are intimately linked to the use of wave-functionexpansion methods, which exhibit mild computational scal-ing in mass number, instead of the exponential scaling re-quired by exact methods, resulting in the recent milestone abinitio calculation of Sn [1919]. In practice, various many-body approaches exist for medium-mass nuclei, e.g., many-body perturbation theory (MBPT) [1818, 2020–2222], coupled clus-ter (CC) theory [1616, 2323], the in-medium similarity renormal-ization group (IMSRG) [1717, 2424, 2525], self-consistent Green’sfunction (SCGF) theory [1515, 2626], and nuclear lattice simula-tions [2727]. In particular, the recent use of non-perturbativemany-body approaches has generated an unprecedented levelof accuracy in medium-mass applications for a diverse set ofnuclear observables, see, e.g., [2828–3030].All these frameworks require the introduction of a compu-tational basis for the representation of the (second-quantized)many-body operators. In the limit of a one-body Hilbert spaceof infinite dimension, di ff erent choices of the computationalbasis yield identical results. However, due to computationallimitations, in practice one is always restricted to using a finite ∗ [email protected]@theorie.ikp.physik.tu-darmstadt.de † [email protected]@physik.tu-darmstadt.de ‡ [email protected]@theorie.ikp.physik.tu-darmstadt.de § [email protected]@physik.tu-darmstadt.de ¶ [email protected]@physik.tu-darmstadt.de basis size and, consequently, the resulting observables will de-pend on the underlying computational basis.It has been realized only recently in nuclear physics thatthe optimization of the single-particle basis provides a pow-erful tool to stabilize many-body calculations and enables amore reliable extraction of physical observables from large-scale calculations [3131, 3232]. In other fields of many-body re-search, like quantum chemistry, this is much more exploredand the construction of suitable single-particle basis functionshas been an integral part of the ab initio endeavor, yielding arich variety of basis sets. In this work, we investigate bene-fits and limitations of various single-particle bases used in thesolution of the nuclear many-body problem.Choosing the single-particle basis in nuclear many-bodytheory primarily requires addressing the following questions:i) What is the best choice for obtaining rapid convergencewith respect to the model-space size?ii) What is the best strategy to minimize the dependence ofphysical observables on basis set parameters?iii) To what extent is the factorization of center-of-mass andintrinsic motion contaminated?In practice, optimizing with respect to all of the above simulta-neously is not possible. Historically, most many-body calcula-tions either employ harmonic oscillator (HO) or Hartree-Fock(HF) single-particle states. Harmonic oscillator basis statesrigorously ensure factorization of center-of-mass and intrin-sic degrees of freedom of the many-body wave function whencombined with an N max -truncation, as in no-core shell model(NCSM) approaches [3333, 3434]. However, in practice a strongdependence on the basis set parameters such as the oscillatorfrequency of the confining potential is observed, especiallyfor heavier nuclei or for observables that are more sensitive tothe long-range part of the nuclear wave function. This makesthe extraction of such observables challenging. Using HF or-bitals based on a prior mean-field solution typically lowersthe frequency dependence, while numerically still leading toa factorization of the center-of-mass and intrinsic wave func-tion in large enough model spaces [3535]. However, selected a r X i v : . [ nu c l - t h ] S e p nuclear observables may still show sensitivity to the oscillatorfrequency in HF basis as observed, e.g., for charge radii ofmedium-mass nuclei in IMSRG calculations [1010].Recently, applications of natural orbitals (NAT), defined aseigenvectors of the one-body density matrix, revealed fastermodel-space convergence and significantly reduced sensitiv-ity to basis parameters in large-scale NCSM calculations [3232].Furthermore, they have been shown to drastically reduce therequired amount of three-particle–three-hole amplitudes inCC applications [3636], allowing for novel calculations withleading triples corrections, e.g., for deformed nuclei [3636] andnuclear matrix elements of the neutrinoless double-beta de-cay [3737]. In this work, the natural orbital basis is inspected indetail and systematically applied in non-perturbative medium-mass studies using modern chiral interactions.This paper is organized as follows. In Sec. IIII various op-tions of single-particle bases are shown and discussed. Sec-tion IIIIII introduces the concept of normal ordering as well asthe impact of the reference state on the many-body formal-ism. In Sec. IVIV the IMSRG approach is briefly introduced.The di ff erent single-particle bases considered are comparedin detail in Sec. VV and applied to medium-mass systems usingthe IMSRG formalism in Sec. VIVI. Finally, we summarize andconclude in Sec. VIIVII. II. BASIS OPTIMIZATIONA. Rationale
While HO basis sets have been used extensively for a longtime in various many-body frameworks, they constitute anagnostic choice with respect to any specific properties ofthe target system, e.g., in terms of mass number or mean-field e ff ects. This can be addressed by using HF orbitals in-stead. Hartree-Fock orbitals account for bulk properties ofthe nucleus stemming from a variational minimization of theground-state energy. Observables like the energy or the radiusare therefore well captured at the HF level as long as the nu-clear interaction is soft enough, and the single-particle wavefunctions possess an improved radial dependence as opposedto the Gaussian fall o ff of HO eigenfunctions. Proton and neu-tron single-particle potentials in general di ff er in the HF ap-proach, thus accounting for mean-field contributions inducedby Coulomb and isospin-breaking e ff ects.Still, the HF procedure by construction only provides anoptimization of occupied single-particle states (holes) whileleaving virtual single-particle states (particles) untouched be-yond fixing the normalization. However, wave-function ex-pansion methods aim at capturing dynamic correlations linkedto particle-hole excitations, which also involve single-particleorbitals that are not optimized by the HF approach. Therefore,incorporating such e ff ects in the construction of the computa-tional basis is key when trying to robustly determine observ-ables to high precision. B. Notation
In the following, we denote second-quantized n -body ( n B)operators via O ( n B) ≡ n !) (cid:88) k ... k n o ( n B) k ... k n a † k · · · a † k n a k n · · · a k n + , (1)where the lower-case letters o ( n B) k ... k n represent their matrix el-ements and a † k ( a k ) denote the single-particle creation (anni-hilation) operators. Normal-ordering techniques are exploitedto re-express the operator with respect to an A -body referencestate,˜ O ( n B) ≡ n !) (cid:88) k ... k n ˜ o ( n B) k ... k n : a † k · · · a † k n a k n · · · a k n + : , (2)where strings of normal-ordered creation and annihilation op-erators are denoted by colons and we use the tilde to distin-guish the reference-state normal-ordered operator and its ma-trix elements from the initial one. While the specific vacuumis absent in this notation, it will be clear from the context whatreference state we are referring to.For the particular case of the nuclear Hamiltonian the fol-lowing notation is employed to denote its normal-ordered con-tributions H = E + (cid:88) pq f pq : a † p a q : + (cid:88) pqrs Γ pqrs : a † p a † q a s a r : + (cid:88) pqrstu W pqrstu : a † p a † q a † r a u a t a s : , (3)where E , f , Γ , and W denote the zero-, one-, two-, and three-body matrix elements. Because the operator in Eq. (33) is inreference-state normal order the expectation value is given byits normal-ordered zero-body part (cid:104) Φ | H | Φ (cid:105) = E , (4)where | Φ (cid:105) denotes the reference state. In the case of a Hartree-Fock reference state | Φ (cid:105) = | HF (cid:105) this corresponds to theHartree-Fock mean-field energy, E = E HF . C. Natural orbitals
Natural orbitals are defined as the eigenbasis of the one-body density matrix with its matrix elements given by γ pq ≡ (cid:104) Ψ | : a † p a q : | Ψ (cid:105)(cid:104) Ψ | Ψ (cid:105) , (5)where | Ψ (cid:105) denotes the exact ground state. The calculation ofthe exact one-body density matrix requires the full solution ofthe Schrödinger equation, which is out of reach beyond thelightest systems. However, early attempts in quantum chem-istry revealed that approximate natural orbitals can be veryuseful [3838, 3939]. Such basis sets are obtained by using an ap-proximate many-body state | Ψ approx. (cid:105) to obtain an approximateone-body density matrix, γ approx. pq ≡ (cid:104) Ψ approx. | : a † p a q : | Ψ approx. (cid:105)(cid:104) Ψ approx. | Ψ approx. (cid:105) . (6)In regions of the nuclear chart where the exact wave functionis computationally inaccessible, this provides an alternativeoption for defining a basis for the many-body calculation. Inpractice, a reasonable trade-o ff between the accuracy of themany-body truncation for the construction of the approximatewave function and the associated computational cost needs tobe found.In the case where the approximate wave function is an HFSlater determinant | HF (cid:105) the density matrix γ HF pq ≡ (cid:104) HF | : a † p a q : | HF (cid:105)(cid:104) HF | HF (cid:105) (7)has the particularly simple form γ HF = (cid:32) hh
00 0 (cid:33) , (8)where hh denotes the identity in the sub-block of hole states.As the HF state is normalized to unity, the presence of themean-field overlap (cid:104) HF | HF (cid:105) does not a ff ect the HF densitymatrix. Furthermore, the density matrix corresponds to a nor-malized many-body state, meaning that its trace yields the par-ticle number of the state, i.e., tr( γ HF ) = A .For an HF reference state, the canonical orbitals , definedas the eigenbasis of the one-body HF Hamiltonian, and thenatural orbitals based on the HF density matrix in Eq. (88) co-incide. Therefore, one must include correlations beyond meanfield in the construction of the density matrix to gain a benefitfrom the natural orbitals. D. Perturbatively improved density matrix
As discussed in the previous section, accounting forparticle-hole couplings in the density matrix is essential forproviding a more refined computational basis. The simplestapproach to including such e ff ects is by employing a pertur-batively corrected one-body density matrix. Following thedescription in Ref. [4040], the one-body density matrix up tosecond order in the interaction ( λ ), based on expanding theeigenstate of the approximate wave function up to second or-der in MBPT (MP2), can be written as γ MP2 ≡ γ HF + γ (02) + γ (20) + γ (11) + O ( λ ) , (9)where γ ( mn ) pq ≡ (cid:104) Φ ( m ) | : a † p a q : | Φ ( n ) (cid:105) (10)is the MBPT contribution for the density matrix arising fromthe bra and ket wave function corrections at order m and n , respectively. Terms of order λ or higher in the interaction arediscarded. Moreover, terms of the form γ (01) / (10) are absentwhen using a canonical HF reference state due to Brillouin’stheorem [4141]. Explicit expressions for the various contribu-tions in terms of single-particle orbitals are given by D (hp) i (cid:48) a (cid:48) = (cid:88) abi Γ i (cid:48) iab Γ aba (cid:48) i (cid:15) a (cid:48) i (cid:48) (cid:15) abi (cid:48) i , (11a) D (hp) i (cid:48) a (cid:48) = − (cid:88) ai j Γ i (cid:48) ai j Γ i ja (cid:48) a (cid:15) a (cid:48) i (cid:48) (cid:15) a (cid:48) ai j , (11b) D (hh) i (cid:48) j (cid:48) = − (cid:88) abi Γ i (cid:48) iab Γ ab j (cid:48) i (cid:15) abi (cid:48) i (cid:15) abj (cid:48) i , (11c) D (pp) a (cid:48) b (cid:48) = (cid:88) ai j Γ a (cid:48) ai j Γ i jb (cid:48) a (cid:15) a (cid:48) ai j (cid:15) b (cid:48) ai j , (11d)where the labels i , j , k , ... ( a , b , c , ... ) correspond to single-particle states occupied (unoccupied) in the reference determi-nant, i.e., the HF state in our case. The matrix elements Γ pqrs are given in the HF basis, thus corresponding to an HF par-titioning in the MBPT expansion of the density matrix [2121].Furthermore, the short-hand notation (cid:15) abi j ≡ (cid:15) i + (cid:15) j − (cid:15) a − (cid:15) b (12)is used, with (cid:15) p ≡ f pp denoting the HF single-particle energyof orbital p . Consequently, the MP2 density matrix is givenby γ MP2 = (cid:32) γ hh γ hp γ ph γ pp (cid:33) , (13)where the hole-particle and particle-hole blocks are non-zeroand given by γ hp = D (hp) + D (hp) = (cid:16) γ ph (cid:17) (cid:124) , (14)and the hole-hole and particle-particle blocks by γ hh = γ HF + D (hh) , (15a) γ pp = D (pp) , (15b)respectively. Note that in contrast to the HF density matrix,the MP2 density matrix contains particle-particle and particle-hole couplings as shown for O in Fig. 11. Since (cid:88) i D (hh) ii + (cid:88) a D (pp) aa = , (16)the second-order density matrix still fulfills the trace normal-ization condition tr( γ MP2 ) = A as in the HF case.In practice, the construction of the MP2 density matrix isrealized using a spherically constrained scheme, i.e., enforc-ing angular-momentum conservation throughout the initial HFsolution and the following MBPT calculation. Specifically,the single-particle orbitals are then characterized by the quan-tum numbers n , l , j , t , which are (2 j + n is the radial quantum number, l the orbital angular mo-mentum, j the total angular momentum, and t the isospin pro-jection. In actual calculations, we truncate the single-particlestates at e ≤ e max , with quantum numbers e = n + l . s p p d d s f f p p g g d d s s p p d d s f f p p g g d d s FIG. 1. Correlated one-body proton density matrix γ MP2 in the HFbasis for O in an e max = LO 450 in-teraction based on the N LO NN potential from Ref. [88] with N LO3N forces constructed in Ref. [99]. The first three proton orbitals areoccupied in the O reference state, while the remaining ones are un-occupied. The perturbative corrections beyond HF can be seen in thediagonal particle-particle contributions and the o ff -diagonal particle-hole and hole-particle contributions. Note that for this N = Z nucleusthe neutron density matrix is very similar. Consequently, the resulting MP2 density matrix is block di-agonal in the quantum numbers l jt as only states with di ff erentradial quantum number n couple. The diagonalization of theMP2 density matrix is performed in sub-blocks to ensure sym-metry conservation. The resulting eigenvectors and eigenval-ues correspond to the transformation coe ffi cients from the HFto the natural orbital (NAT) basis and the occupation numbersof the natural orbitals, respectively. E. Basis transformation
The natural orbital states are obtained as linear combina-tions of the HF states, mixing radial excitations only | n α p (cid:105) NAT = (cid:88) n (cid:48) NAT C α p nn (cid:48) | n (cid:48) α p (cid:105) HF , (17)where α p is a collective index for the quantum numbers l p , j p , t p and NAT C α p nn (cid:48) denotes the expansion coe ffi cients in theHF basis obtained by the diagonalization, i.e., HF (cid:104) n (cid:48) α p | n α p (cid:105) NAT = NAT C α p nn (cid:48) , (18)where the m -projection of the total angular momentum issuppressed since the transformation coe ffi cients and single-particle states do not depend on it as long as rotational sym-metry is enforced. By expanding the HF states in the HO ba-sis, we can also express the natural orbital states in the HO s p p d d s f f p p g g d d s Orbitals n p O HFNAT e max = 2 NAT e max = 10 FIG. 2. Occupation numbers n p of the single-particle proton orbitalsfor the HF and NAT basis in O, using the 1.8 / (cid:126) ω =
16 MeV. We show results for twomodel-space truncations e max = e max =
10 in the NAT basisconstruction. As for Fig. 11, the occupations of the neutron orbitalsare nearly identical. basis | n α p (cid:105) NAT = (cid:88) n (cid:48) n (cid:48)(cid:48) NAT C α p nn (cid:48) HF C α p n (cid:48) n (cid:48)(cid:48) | n (cid:48)(cid:48) α p (cid:105) HO = (cid:88) n (cid:48)(cid:48) NAT / HF C α p nn (cid:48)(cid:48) | n (cid:48)(cid:48) α p (cid:105) HO , (19)where the coe ffi cients NAT / HF C α p nn (cid:48)(cid:48) now combine the transfor-mation from the HO to the HF and from the HF to the NATbasis.Note that the set of occupation numbers for the natural or-bitals n p ∈ [0 ,
1] obtained from the eigenvalues now leads toa fractional filling of all orbitals, in contrast to the occupa-tion numbers n p ∈ { , } obtained from the HF solution. Thisfeature is illustrated in Fig. 22 comparing the NAT and HF oc-cupation numbers for an O reference state.Since the reference state for the MP2 density matrix is nota single Slater determinant due to mixing of particle-hole ex-citations the occupation numbers must di ff er from the mean-field picture. As will be discussed in the following, this alsoa ff ects the normal-ordering procedure with respect to natural-orbital basis states.While the employed MP2 density matrix provides a sim-ple approximation to the exact one-body density matrix, non-perturbative many-body schemes can be used to refine the ap-proximation, e.g., a Λ -approach in CC theory [4242], dressedpropagators from Green’s function theory [4343], or a fully cor-related CI calculation [4444]. A balance between accuracy andcomputational complexity needs to be found, and a low-orderMBPT approach provides a reasonable approximation to theone-body density matrix at low computational cost. F. Intrinsic kinetic energy
The intrinsic Hamiltonian, here considered up to three-bodycontributions, can be split into a kinetic part and an interactionpart H = T − T cm + V (2) + V (3) = T int + V (2) + V (3) , (20)with the intrinsic kinetic energy T int and the two- and three-body potentials V (2) and V (3) , respectively. The intrinsic ki-netic energy is obtained by subtracting the center-of-mass ki-netic energy T cm from the full kinetic energy T . The intrinsickinetic energy can be represented either as a sum of one- andtwo-body operators T (1 + = (cid:32) − A (cid:33) (cid:88) i p i m − A (cid:88) i < j p i · p j m (21)or as a pure two-body operator T (2)int = A (cid:88) i < j ( p i − p j ) m . (22)Of course, both cases are equal and can be transformed intoeach other by (cid:88) i < j ( p i − p j ) m = (cid:88) i < j ( p i + p j − p i · p j )2 m = ( A − (cid:88) i p i m − (cid:88) i < j p i · p j m . (23)The one- and two-body matrix elements of the Hamiltonianobviously di ff er depending on the choice of T int . Neverthelessboth cases result in the same HF determinant with identicaltotal HF energy, as studied in detail in Refs. [4545, 4646]. The HFsingle-particle energies are di ff erent for both choices and canbe related by a unitary transformation of the occupied single-particle states [4545]. These findings are based on the assump-tion of a reference state with well-defined particle number A .For a discussion of particle-number breaking theories, e.g., theHartree-Fock-Bogoliubov approach, see Ref. [4747].The partitioning of the kinetic energy operator also a ff ectsthe construction of the natural orbital basis. By employ-ing T (2)int the initial Hamiltonian (before normal ordering) nolonger has a one-body part, and the two-body matrix elementsin the construction of γ MP2 , see Eqs. (11a11a)-(11d11d), di ff er fromthe ones obtained by using the one- plus two-body form ofthe kinetic energy, resulting in altered transformation coe ffi -cients and NAT occupation numbers. The partitioning alsochanges the single-particle energies, further changing the re-sulting γ MP2 .In general, we apply the intrinsic kinetic energy operator ofEq. (2121) with a one- and two-body part for the IMSRG calcu-lations performed in this work. However, in the following wewill additionally study the impact of using a pure two-bodykinetic energy operator T (2)int . III. NORMAL ORDERING
The concept of normal ordering facilitates the formulationof Wick’s theorem [4848] and defines an in-medium optimizedrepresentation of the operator [4949]. For this reason, it is com-monly employed in many-body frameworks applied to nucleiand nuclear matter. In the following, we first address the for-mal details of working with a correlated reference state.
A. Multi-reference formulation
When employing a perturbatively improved one-body den-sity matrix the associated reference state is no longer a sin-gle Slater determinant. Consequently, the notion of nor-mal ordering needs to be extended to cope with the multi-configurational character of the vacuum. Such an exten-sion can be naturally addressed in terms of the generalizedMukherjee-Kutzelnigg normal ordering [5050]. Even thoughthis scheme is not numerically benchmarked in this work, it isstill worth anticipating the additional complications that arisefrom a multi-reference treatment of the MP2 density matrix.For simplicity, we neglect three-body contributions in thefollowing analysis and start from an arbitrary many-body op-erator O containing up to two-body contributions O ≡ O (0B) + O (1B) + O (2B) . (24)Performing the normal ordering of the operator O in Eq. (2424)with respect to a non-product-type vacuum yields [5151]˜ o (0B) = o (0B) + (cid:88) pq o (1B) pq γ pq + (cid:88) pqrs o (2B) pqrs γ pqrs , (25a)˜ o (1B) pq = o (1B) pq + (cid:88) rs o (2B) prqs γ rs , (25b)˜ o (2B) pqrs = o (2B) pqrs , (25c)involving one- and two-body density matrices γ pq and γ pqrs ,respectively. The two-body density matrix, which contributesto the zero-body part of the normal-ordered operator, is givenby γ pqrs ≡ (cid:104) Ψ | : a † p a † q a s a r : | Ψ (cid:105)(cid:104) Ψ | Ψ (cid:105) (26)and can be decomposed into a factorized part of products ofone-body operators and an irreducible two-body part λ pqrs , γ pqrs = (cid:16) γ pr γ qs − γ ps γ qr (cid:17) + λ pqrs . (27)The appearance of λ pqrs is a consequence of the reference statebeing no longer of mean-field character . In the following, In practice, such states are obtained, e.g., from particle-number-broken and-restored Hartree-Fock-Bogoliubov vacua [5252] or small-scale CI diagonal-izations [5353]. the irreducible two-body part is discarded for simplicity and amean-field-like approximation is employed γ pqrs ≈ γ pr γ qs − γ ps γ qr . (28)Equation (2828) is exact as long as a many-body state ofproduct-type is used . While in principle it is straightfor-ward to derive two-body density matrices in MBPT, the factor-ized approximation is expected to provide a reasonably goodchoice for basis optimization.Considerable simplifications are obtained by working in thenatural orbital basis, i.e., γ pq = n p δ pq , n p ∈ [0 , , (29)where the lack of a well-defined particle-hole picture meansthe occupation numbers are no longer zero or one. Expres-sions for the normal-ordered matrix elements in Eq. (3333) are˜ o (0B) = o (0B) + (cid:88) p o (1B) pp n p + (cid:88) pq o (2B) pqpq n p n q , (30a)˜ o (1B) pq = o (1B) pq + (cid:88) r o (2B) prqr n r , (30b)˜ o (2B) pqrs = o (2B) pqrs , (30c)now involving single-particle summations running over thefull one-body Hilbert space for the summation indices p , q ,and r instead of hole orbitals only, a consequence of thesmeared-out Fermi distributions in the occupation numbers n p , as shown in Fig 22. B. Single-reference case
In the simplest case, a single Slater-determinant referencestate is employed in the many-body expansion, | Φ (cid:105) = A (cid:89) i = a † i | (cid:105) , (31)where { a † i } denotes the single-particle creation operators in thecomputational basis. For the occupation numbers of the indi-vidual orbitals one has n p = , if p is a hole state0 , if p is a particle state . (32)Performing the single-reference normal ordering with respectto the reference state in Eq. (3131), the corresponding normal-ordered matrix elements of the operator are obtained as [4949] Approximating a two-body density matrix from one-body density matri-ces is closely related to the approach followed in SCGF theory, wherehigher-order Green’s function are typically factorized products of one-bodyGreen’s functions, thus neglecting their irreducible higher-body contribu-tions [4343]. ˜ o (0B) = o (0B) + (cid:88) i o (1B) ii + (cid:88) i j o (2B) i ji j + (cid:88) i jk o (3B) i jki jk , (33a)˜ o (1B) pq = o (1B) pq + (cid:88) i o (2B) piqi + (cid:88) i j o (3B) pi jqi j , (33b)˜ o (2B) pqrs = o (2B) pqrs + (cid:88) i o (3B) pqirsi , (33c)˜ o (3B) pqrstu = o (3B) pqrstu , (33d)where the labels i , j indicate hole states occupied in the ref-erence state | Φ (cid:105) . In Eqs. (33a33a)-(33d33d) three-body contributionsare explicitly included. In practice, the normal-ordered two-body approximation (NO2B) is employed [5454, 5555], where theresidual three-body part, Eq. (33d33d), is discarded to lower thecomputational complexity.Because the MP2 density matrix does not correspond to asingle Slater determinant, an auxiliary many-body state | NAT (cid:105) is constructed by filling the first A states with the highest oc-cupation numbers. Similar to Eq. (3131) these orbitals are filledwith updated occupations n i ∈ { , } to conserve the particle-number expectation value, thus establishing a well-definedparticle-hole picture. Consequently, in the following appli-cations standard Slater-determinant-based codes can be usedfor the many-body expansion. Note that, even though this ref-erence state has product-type character, the information aboutthe correlated density matrix is encoded in the transformationmatrix from the HO to the NAT basis [see Eq. (1919)] for theone- and two-body parts of the intrinsic Hamiltonian. By us-ing such an auxiliary vacuum the reference-state expectationvalue is larger than the HF expectation value since there is nounderlying variational principle, i.e., (cid:104) NAT | H | NAT (cid:105) > (cid:104) HF | H | HF (cid:105) . (34) IV. IN-MEDIUM SIMILARITY RENORMALIZATIONGROUP
For the medium-mass applications in this work we use thenon-perturbative in-medium similarity renormalization groupapproach. For a detailed discussion of the many-body formal-ism the reader is referred to Refs. [1717, 2424, 5656].
A. Formalism
In the IMSRG framework the many-body Schrödingerequation is solved by performing a decoupling of particle-holeexcitations from the reference state by a continuous unitarytransformation U ( s ) parametrized in terms of a real-valuedflow parameter s , H ( s ) = U ( s ) H U † ( s ) , (35)where H ( s =
0) is the initial, i.e., unevolved, Hamiltonian.Equation (3535) can be rewritten as a first-order ordinary di ff er-ential equation (ODE) in s ,d H ( s )d s = (cid:2) η ( s ) , H ( s ) (cid:3) , (36)with the anti-Hermitian generator η ( s ) defined by the unitarytransformation.The in-medium character of the decoupling condition isachieved by performing the renormalization group evolutionon the normal-ordered representation of the operator. In thesimplest case this is done with respect to a single Slater de-terminant, as indicated in Sec. III BIII B. At s =
0, the normal-ordered zero-body part of the Hamiltonian is the reference-state energy expectation value, e.g., the HF energy in the caseof the HF reference state. Over the course of the evolution,the o ff -diagonal matrix elements of the Hamiltonian are sup-pressed, and the exact ground-state energy of the fully inter-acting system is given by the normal-ordered zero-body partof the evolved Hamiltonian E ( s → ∞ ) = lim s →∞ (cid:104) Φ | H ( s ) | Φ (cid:105) . (37)Solving the IMSRG flow equation absorbs the dynamic cor-relations non-perturbatively and smoothly decouples particle-hole excitations from the reference state.In practice, the expansion in Eq. (3636) is computation-ally intractable since the repeated commutator evaluation in-duces higher-body operators in each integration step, yield-ing up to A -body operators when applied to an A -body sys-tem. In the following, the IMSRG(2) approximation is used,where all operators are truncated at the normal-ordered two-body level. The IMSRG(2) provides an accurate truncationscheme that incorporates all MBPT corrections up to order λ while resumming many higher-order contributions in a non-perturbative way through the IMSRG flow [1717]. In medium-mass applications, the many-body uncertainty related to theIMSRG(2) approximation is estimated to be ≈
2% for ground-state energies (see, e.g., [1111]), thus providing a versatile andprecise many-body scheme at moderate computational cost.In all subsequent calculations the modified normal order-ing as discussed in Sec. III BIII B is employed, thus enabling theuse of the standard Slater-determinant-based IMSRG. Eventhough various open-shell extensions have been designed andapplied within the IMSRG approach [5252, 5353, 5757], this workis restricted to closed-shell applications to provide a simpletestbed for the natural orbital basis.
B. Magnus reformulation
In this work, the Magnus formulation of the IMSRG [5858,5959] is employed, which provides direct access to the unitarytransformation U ( s ) in a computationally e ffi cient way usingthe parametrization U ( s ) = e Ω ( s ) , (38)where Ω ( s ) is the anti-Hermitian Magnus operator. Analo-gously to Eq. (3636), an ODE for the Magnus operator can be derived. Having the transformation matrix at hand allowsfor the IMSRG-evolution of any other operator by the Baker-Campbell-Hausdor ff formula for the many-body operator ac-cording to O ( s ) = e Ω ( s ) O ( s = e − Ω ( s ) (39)instead of solving an additional set of ODEs for each opera-tor. In practice, the Baker-Campbell-Hausdor ff expansion isperformed using nested commutator evaluations to some trun-cation level. Moreover, solving the Magnus expansion has theadvantage of allowing the use of a simpler ODE solver for theMagnus operator without loss of accuracy. C. Correlation e ff ects from MP2 natural orbitals Before discussing IMSRG applications using natural or-bitals, it is worth addressing the interplay of the correlationsbuilt into the MP2 density matrix and the correlations that areresummed within the IMSRG flow.Using a natural orbital reference determinant yields ahigher ground-state energy at s = s = ∞ usingthe MP2 density matrix does not improve upon the IMSRG(2)results obtained in any other single-particle basis. The MP2density matrix only incorporates correlations to one-particle-one-hole and two-particle-two-hole excitations. Within theIMSRG(2) approximation such e ff ects are resummed to all or-ders [1717] such that no improvement on the final observable isexpected. Once higher-body excitations are included, addi-tional correlations will enter the description which are absentin the IMSRG(2) scheme. Practically, this is achieved by in-cluding third-order terms in the MBPT expansion, i.e., λ , orallowing three-body operators in the normal-ordered Hamilto-nian, thus generating additional contributions in the first-orderstate correction. Both options will generate the leading con-tributions to three-particle-three-hole excitations. V. DIAGNOSTICS FOR THE DENSITY MATRIX
We begin by investigating the MP2 density matrix and theassociated NAT basis. These results allow us to gain a betterunderstanding of the relationship between the various compu-tational bases and their sensitivity to the nuclear Hamiltonianused in their construction. We focus on two sets of chiral in-teractions, the N LO NN potential from Ref. [88] with N LO3N forces constructed in Ref. [99], which in the following is re-ferred to as “N LO” with the corresponding cuto ff value, andthe “1.8 / A. The “softness” of the interaction “Soft” interactions are low-resolution interactions thatshow weak coupling between low- and high-energy states.The softness of an interaction, i.e., the degree of decouplingbetween low and high momenta in the Hamiltonian, can bevaried by changing the regulator scale for Hamiltonians con-structed from an EFT as well as by applying (S)RG methodsto decouple or integrate out high-momentum degrees of free-dom [4949, 6060]. Soft interactions applied in many-body methodshave been shown to improve convergence with respect to basistruncation and order in the many-body expansion. In particu-lar, the use of an SRG-evolved Hamiltonian is required to en-able a perturbative solution even for closed-shell systems [2121].A Weinberg-eigenvalue analysis, which provides a metric ofthe perturbativeness of an interaction, shows that the softnessis intimately linked to the SRG resolution scale [6161, 6262].Because the one-body density matrix is constructed froman MBPT expansion, we expect the density matrix and the re-sulting NAT basis to be more sensitive to the basis frequencyand truncation for hard interactions. For unevolved chiral po-tentials the mean-field wave function may exhibit unphysicalproperties, giving rise to an unbound HF solution. With such apoor reference state, the many-body expansion is significantlymore complicated, in particular if perturbative techniques areemployed. The key idea of a many-body expansion is to startfrom a qualitatively correct reference state while residual dy-namic correlation e ff ects are brought in as (small) corrections.This rationale is obviously broken once the mean-field refer-ence is unbound or not under control, manifesting in final re-sults via, e.g., strong frequency dependence. B. Single-particle wave functions
While the HF approach targets the optimization of the oc-cupied single-particle states from a variational approach, theunoccupied orbitals are left unmodified up to normalization.Therefore, the HF basis is expected to properly describe oc-cupied orbitals while failing for unoccupied ones. The naturalorbital basis, however, accounts for particle-hole admixturesand therefore may qualitatively improve the description of un-occupied states as will be tested in the following calculations.In the following, a single-particle basis is employed includingstates up to a principal quantum number e max . Additionally,we introduce a truncation in three-body space keeping onlyconfigurations with e + e + e ≤ E < e max due to theextensive size of three-body matrix elements.In Fig. 33, we show the squared absolute value of the radialwave functions for di ff erent oscillator frequencies using theHO, HF, and NAT bases. Di ff erent rows correspond to dif-ferent single-particle orbitals; only the first row (0 p / ) corre-sponds to an occupied orbital. Clearly, using a HO basis leadsto strong frequency dependence in all cases, even for the occu-pied 0 p / state. Hence, HO wave functions are ruled out as areliable computational basis and will not be considered furtherin this work. While the 0 p / orbitals are more robust in theHF case as expected, unoccupied HF orbitals show frequencydependence comparable to that of HO orbitals, a consequenceof the fact that unoccupied orbitals are not optimized in theHF approach.Switching to natural orbitals nicely resolves many of the re- p HO N LO 450 e max = HF O e max = NAT d = 16 MeV = 20 MeV = 24 MeV = 36 MeV | u ( r ) | [ f m ] s f p d r [fm] FIG. 3. Squared absolute value of the radial wave function u ( r ) of O as a function of r for di ff erent proton orbitals in the HO, HF, andNAT bases in the first, second, and third columns, respectively. Weshow results for the occupied 0 p / and some of the first unoccupiedorbitals for the N LO 450 interaction and oscillator frequencies (cid:126) ω = −
36 MeV. The HF and NAT orbitals include single-particle HOstates up to e max =
10 and E = maining artifacts, revealing only minor frequency dependencefor both occupied and unoccupied states. As the softer 1.8 / e max =
10 for this interaction, while additionallybenchmarking the e ff ect of natural orbitals when going to alarger basis size of e max =
14 (right). While the frequencydependence for this softer interaction is much milder in theHF case, high-lying single-particle states still significantly de-pend on (cid:126) ω . A residual frequency dependence is still seen inthe 1 p / orbital at e max =
10 in the natural orbital basis, butthis fully vanishes when going to larger spaces of e max = p e max = 10 HF e max = 10 NAT O e max = 14 NAT d = 16 MeV = 20 MeV = 24 MeV = 36 MeV | u ( r ) | [ f m ] s f p d r [fm] FIG. 4. Same as Fig. 33 but for the 1.8 / e max =
10 in the left andmiddle columns, respectively, as well as e max =
14 for the NAT basisin the right column.
C. Positive definiteness as diagnostic tool
The density matrix is a positive definite operator and thusits eigenvalues, the occupation numbers, are non-negative.Therefore, unphysical negative occupations or occupationslarger than one should not show up during the diagonalization.Previous investigations in quantum chemistry showed that theappearance of negative occupation numbers can be linked toa breakdown of a single-reference description and hint at theonset of strong static correlations [6464]. Therefore, we aim toutilize occupation numbers as diagnostic and investigate theirsensitivity to the softness of the nuclear interaction. As theHF ground-state energy is directly related to the softness ofthe interaction, a correlation between the HF energy and thesize of negative occupations is expected.Figure 55 depicts the magnitude of the negative occupationsusing N LO interactions for various cuto ff values in , O. Inboth nuclei, we observe a decrease in size for softer interac-tions, as indicated by going from the harder potentials withcuto ff Λ =
500 MeV to
Λ =
400 MeV, in both the NN-onlyand the NN +
3N cases. Consequently, an unbound HF solutionstrongly a ff ects the appearance of unphysical negative occu-pations. In general, using the two-body form of the kinetic
40 30 20 10 0 10 20 300.150.100.050.00 ( j i + ) n N A T O p -orbitals N LO 400N LO 420N LO 450N LO 500
40 30 20 10 0 10 20 30 E HF [MeV]0.150.100.050.00 ( j i + ) n N A T O p -orbitals NN-only ( T (1+2)int )NN+3N ( T (1+2)int ) NN-only ( T (2)int )NN+3N ( T (2)int ) FIG. 5. Negative occupations of the p -orbitals scaled by (2 j i +
1) inthe NAT basis for O (top) and O (bottom) as a function of the HFenergy. We show results for various cuto ff s with the NN-only N LOEMN and NN +
3N N LO interactions indicated by triangles and cir-cles, respectively. We apply both choices for the kinetic energy oper-ator, T (1 + (filled symbols) and T (2)int (open symbols), and use a modelspace of e max / E = /
14 with (cid:126) ω =
20 MeV. All negative occu-pations arise only for high radial quantum number. Note that thereare no negative occupations for the softer NN-only EMN 400 and420 interaction. energy operator T (2)kin results in smaller negative occupationsfor both nuclei. We also verified that softening the interac-tion by a consistent SRG evolution of NN and 3N contribu-tions [1010, 6565] significantly reduces the magnitude of the neg-ative occupations, eventually letting them vanish completely.Increasing the model-space size seems to increase the magni-tude of these occupations. Moreover, the e ff ect is generallyless pronounced for heavier nuclei, e.g., in Ni.In addition, we investigate the size of negative occupationsin the case of C. Due to the cluster structures and weak shellclosure in C the quality of single-reference many-body ap-proaches is expected to deteriorate in comparison to the dou-bly magic nucleus O. An analysis of the single-particle spec-trum revealed only a small shell gap in the single-particlespectrum, thus significantly enhancing the size of perturba-tive corrections to the MP2 density matrix in the particle-particle and hole-hole channel, see Eqs. (11c11c) and (11d11d). Con-sequently, highly erratic occupation numbers were observed(not shown). Empirically, we found that the use of T (2)int witha slightly larger shell gap was superior to T (1 + , significantlyreducing, though not fully resolving, the large negative occu-pations. The results of this analysis for the occupation num-bers is also evidence for the challenges of the single reference-state starting point for a description of C.0 E [ M e V ] N LO 450 O HF NAT16 20 24 28 32 36 [MeV]2.22.42.62.83.0 R c h [ f m ] NN+3N NN-only
16 20 24 28 32 36 [MeV] e max =6e max =8e max =10 e max =12Expt. E [ M e V ] N LO 450 Ca HF NAT16 20 24 28 32 36 [MeV]2.62.83.03.23.43.63.8 R c h [ f m ] NN+3N NN-only
16 20 24 28 32 36 [MeV] e max =6e max =8e max =10 e max =12Expt. FIG. 6. Ground-state energies (upper rows) and charge radii (lower rows) of O and Ca in the left and right plots, respectively, as a functionof the oscillator frequency, for the NN-only N LO EMN 450 (triangles) and NN +
3N N LO 450 (circles) interactions. We show results for theHF and NAT bases in the left and right panels of each plot, respectively, using various single-particle truncations e max with E = VI. IMSRG RESULTS
After addressing in detail properties of the single-particlebasis itself, the various choices are benchmarked for medium-mass closed-shell systems using the IMSRG framework, fo-cusing on O, Ca, and Ni. All many-body calcu-lations employed the publicly available IMSRG-solver byR. Stroberg [6666].
A. Comparing the HF and NAT basis
We compare results for ground-state energies and chargeradii of O and Ca in the HF and NAT bases in Fig. 66 fora large range of oscillator frequencies for the NN-only andNN +
3N N LO 450 interactions. For the NN-only potential,we observe nearly no change when going from the HF to theNAT basis on this scale. Since the HF solution is bound, bulkproperties are well captured at the mean-field level and ap-plying the NAT basis does not yield an improvement in finalresults. Both energies and radii are almost flat as a functionof (cid:126) ω for the largest model space and rapidly converge withmodel-space size in both the HF and natural orbital bases. When 3N forces are included, the NN +
3N results similarly tothe NN-only case show almost no change from the HF to theNAT basis, but the (cid:126) ω -dependence becomes more pronouncedfor the radii.In order to systematically understand the di ff erence be-tween the basis sets, we examine the converged IMSRG(2)ground-state energies in greater detail. In Fig. 77, we show thedi ff erence of the results in the HF and NAT bases as a functionof the SRG evolution scale for three closed-shell nuclei O, Ca, and Ni for the NN-only interaction. The analysis isperformed in absence of three-body interactions to eliminatethe sensitivity of the di ff erent reference states to the NO2B ap-proximation. For harder interactions (larger λ ) the di ff erenceis of the order of 1 MeV with the natural orbitals yieldingstronger binding for O and Ca and slightly weaker bind-ing for Ni. Softening the potential (small λ ) significantlyreduces the e ff ect such that, eventually, only di ff erences of theorder of tens of keV remain at λ = . − . These di ff erencesare marginally enhanced when including 3N forces, i.e., nat-ural orbitals provide slightly more binding compared to theHF basis and lead to a minor decrease of the (cid:126) ω -dependenceof charge radii. We emphasize again that this NN +
3N inter-action leads to an unbound HF solution, such that the totalbinding has to be produced by correlation e ff ects during the1 ]0.60.40.20.00.20.40.60.81.0 E N A T E H F [ M e V ] N LO 450 O Ca Ni FIG. 7. Di ff erence of the ground-state energies in the NAT and HFbases for O, Ca, and Ni as a function of the SRG resolutionscale λ using the NN-only N LO EMN 450 interaction and a modelspace of e max =
14 with (cid:126) ω =
20 MeV.
IMSRG flow and the mean field provides a poor referencestate as discussed in Sec. V AV A. The minor di ff erences in con-verged energies are assumed to be driven by induced many-body contributions that di ff er in HF and natural orbital bases.A further investigation requires systematic evaluation of lead-ing three-body contributions beyond the IMSRG(2) approxi-mation which is beyond the scope of the present work.In summary, we do not observe the desired independenceof the oscillator frequency in the smaller model spaces, whichone could have guessed from Fig. 33, and do not improve onthe frequency dependence in the largest model spaces shownhere compared to the HF basis. B. Di ff erences between NCSM and IMSRG Given the great performance of MP2 natural orbitals inNCSM results as shown in Ref. [3232], the above results seemsurprising at first since no substantial improvement over HForbitals is obtained. The key di ff erence between the IMSRGcalculations performed in this work so far and the NCSM cal-culations is the model space in which the many-body solutionis obtained.In the NCSM one conveniently employs an N max -truncationwhere only many-body configurations up to a given relativeexcitation level are included [3434]. In this case, construct-ing the correlated one-body density matrix in a large single-particle basis includes excitations that are absent from theNCSM configuration space and, therefore, improves the fre-quency dependence. On the other hand, this is inherentlydi ff erent to an IMSRG application where both the reference-state construction and IMSRG flow typically take place in thesame model space , parametrized by e max . Since the high-lyingstates are included already in the initial single-particle basisfor the HF calculation, we cannot expect significant improve-ment for the simplest natural-orbital-based IMSRG calcula- tions over HF-based IMSRG calculations.Consequently, the key idea in the following will be the con-struction of the MP2 density matrix in a large space whilesolving the many-body problem in a reduced space in pres-ence of the full-space correlations embedded into the basistransformation. C. Reduced-basis calculations from NAT / HF constructed infull space
In the following, the initial MP2 density matrix is built ina large single-particle basis M full , while a smaller subspace M reduced ⊆ M full is used for performing the IMSRG evolu-tion. We restrict the natural orbital states to orbitals includedin M reduced and discard higher-lying orbitals in the transfor-mation. The transformation of the single-particle states forboth model-space choices is contrasted in Eqs. (4040) and (4141): | n α p (cid:105) (cid:103) NAT = n reduced (cid:88) n (cid:48) NAT / HF (cid:101) C α p nn (cid:48) | n (cid:48) α p (cid:105) HO , (40) | n α p (cid:105) NAT = n full (cid:88) n (cid:48) NAT / HF C α p nn (cid:48) | n (cid:48) α p (cid:105) HO . (41)The tilde for the state indicates that it is constructed in thereduced space and the upper bounds n reduced and n full specifythe maximal radial quantum number contained in M reduced and M full for the corresponding | α p (cid:105) state, respectively. The tildeon top of the transformation coe ffi cients NAT / HF (cid:101) C indicates thatthe HF calculation and the construction of the natural orbitalsis performed in M reduced , while untilded coe ffi cients originatefrom the full model space M full . Even though parts of theinformation contained in the natural orbital basis are lost dur-ing this reduction process, the resulting matrix representationof operators in M reduced still contains information about thelarge space due to the mixing of radial excitations up to n full included in M full that are otherwise not contained in M reduced .Excluding higher-lying states from M reduced is also motivated,because for low-resolution Hamiltonians we expect the many-body expansion to be dominated by excitations to low-lyingstates.In practice, M reduced and M full will be parametrized bytwo values, e max for the IMSRG evolution (in M reduced ) and e HF / NATmax for the basis construction (in M full ). For the follow-ing results, we employ e HF / NATmax =
14 for the 1 . / . e max = , , and 10 for O, Ca, and Ni in Figs. 88 and 99. Construct-ing the NAT basis in the full space leads to a significant re-duction of the (cid:126) ω -dependence for both ground-state energiesand charge radii as well as improved convergence behaviorwith respect to e max . The resulting improvement is similar towhat was seen in NCSM calculations for O [3232] with nearlyfrequency-independent energies and radii, shown in the rightcolumn of the first plot in Fig. 88.2 E [ M e V ] HF 1.8/2.0 EM O e HF/NATmax = 14
NAT12 16 20 24 28 32 [MeV]2.42.52.6 R c h [ f m ]
12 16 20 24 28 32 [MeV] e max =6e max =8e max =10 e max =14Expt. E [ M e V ] HF 1.8/2.0 EM Ca e HF/NATmax = 14
NAT12 16 20 24 28 32 [MeV]3.03.13.23.33.4 R c h [ f m ]
12 16 20 24 28 32 [MeV] e max =6e max =8e max =10 e max =14Expt. FIG. 8. Ground-state energies (upper rows) and charge radii (lower rows) of O and Ca in the left and right plots, respectively, as a functionof the oscillator frequency in the HF and NAT bases for the 1.8 / M full with e HF / NATmax =
14 to constructthe NAT basis, whereas the IMSRG calculations are performed for e max =
6, 8, 10, and 14, with E =
16 in both cases.
Analogous conclusions hold for heavier nuclei, where theconvergence pattern is improved and we obtain converged re-sults already in smaller model spaces e max . Although we can-not improve the results beyond the model space of e HF / NATmax employed for the initial transformation, we only have to solvefor the natural orbital basis once in the largest possible e HF / NATmax space without having to solve the computationally more ex-pensive IMSRG equations in the full space. Assuming, we canobtain comparable results in e max =
10 (1140 single-particlestates) using an MP2 density matrix constructed in e max = R ≈ − ffi cient alter-native to the full-space IMSRG-calculations. Very advancedtruncation schemes such as IMSRG(3) scale as N , where N is a measure of the size of the single-particle basis. Therefore,naive speedups of the order R ≈ − can be anticipated.Consequently, further improving the construction of single-particle basis sets will significantly help to advance to heaviernuclei and higher accuracies in ab initio applications. VII. SUMMARY AND CONCLUSIONS
In this work, we performed an extensive study of single-particle bases in nuclear ab initio applications. We focusedon a set of natural orbitals, Hartree-Fock, and harmonic-oscillator basis states, with natural orbitals based on a per-turbatively improved one-body density matrix. The single-particle wave function and its dependence on the oscillatorfrequency as well as the potential occurrence of negative oc-cupations in the NAT basis have been investigated in detail.Going to su ffi ciently large model spaces, the natural orbitalsprovide frequency-independent wave functions for both occu-pied and unoccupied states. A reasonable mean-field solution,preferably bound, and a lower resolution Hamiltonian are keyfactors to generate a reasonable correlated one-body densitymatrix and resulting natural orbital basis. When these con-ditions are not met, the construction of the natural orbitalsdoes not completely lead to the desired frequency indepen-dence and produces states with unphysical negative occupa-tion numbers. Using the two-body form of the kinetic energydecreased the size of the negative occupations compared tothe one- plus two-body form of the kinetic energy.Hartree-Fock and natural orbital basis states have been3 E [ M e V ] HF 1.8/2.0 EM Ni e HF/NATmax = 14
NAT12 16 20 24 28 32 [MeV]3.13.23.33.43.53.63.73.8 R c h [ f m ]
12 16 20 24 28 32 [MeV] e max =6e max =8 e max =10e max =14 FIG. 9. Same as Fig. 66 but for Ni. benchmarked in medium-mass applications using the IMSRG as a non-perturbative many-body framework. Comparingground-state energies and charge radii of medium-mass nu-clei, we observed only small di ff erences between the HF andNAT bases for Hamiltonians with only two-body interactionsand slightly enhanced variations when three-body forces wereincluded. In both cases, the results came closer in agreementby a consistent SRG evolution of NN and 3N interactions tosmaller resolution scales. Even though we did not obtain thedesired frequency independent results by applying the NATbasis directly, significantly improved frequency independenceand faster model-space convergence has been found by con-structing the natural orbitals in a large model space and eval-uating a subsequent IMSRG evolution in a reduced space.This strategy presents a promising improvement to advancethe reach of ab initio methods to heavier nuclei and demon-strates the benefits of investigating the computational basis inmore detail.One possibility for further exploration is the investiga-tion of higher-body contributions and specifically how thedi ff erence between the HF and natural orbital basis resultsarises. Another direction is to incorporate the natural or-bitals in a multi-reference approach to avoid the single Slater-determinant approximation and fully capitalize on the dy-namic correlations included in the perturbatively improveddensity matrix. ACKNOWLEDGEMENT
We thank P. Arthuis, G. Hagen, T. Papenbrock, andR. Stroberg for useful discussions on this manuscript. Thiswork is supported in part by the Deutsche Forschungsgemein-schaft (DFG, German Research Foundation) – Projektnummer279384907 – SFB 1245, and the Max Planck Society. [1] K. Hebeler, J. D. Holt, J. Menéndez, and A. Schwenk, “Nuclearforces and their impact on neutron-rich nuclei and neutron-richmatter,” Annu. Rev. Nucl. Part. Sci. , 457 (2015)Annu. Rev. Nucl. Part. Sci. , 457 (2015).[2] H. Hergert, “A Guided Tour of Ab Initio Nuclear Many-BodyTheory,” (2020), arXiv:2008.05061arXiv:2008.05061.[3] E. Epelbaum, H.-W. Hammer, and U.-G. Meißner, “ModernTheory of Nuclear Forces,” Rev. Mod. Phys. , 1773 (2009)Rev. Mod. Phys. , 1773 (2009).[4] R. Machleidt and D. R. Entem, “Chiral e ff ective field theoryand nuclear forces,” Phys. Rep. , 1 (2011)Phys. Rep. , 1 (2011).[5] K. Hebeler, S. K. Bogner, R. J. Furnstahl, A. Nogga,and A. Schwenk, “Improved nuclear matter cal-culations from chiral low-momentum interactions,”Phys. Rev. C , 031301(R) (2011)Phys. Rev. C , 031301(R) (2011).[6] A. Ekström, G. R. Jansen, K. A. Wendt, G. Hagen,T. Papenbrock, B. D. Carlsson, C. Forssén, M. Hjorth-Jensen, P. Navrátil, and W. Nazarewicz, “Accurate nu-clear radii and binding energies from a chiral interaction,”Phys. Rev. C , 051301(R) (2015)Phys. Rev. C , 051301(R) (2015).[7] E. Epelbaum, H. Krebs, and U.-G. Meißner, “PrecisionNucleon-Nucleon Potential at Fifth Order in the Chiral Expan-sion,” Phys. Rev. Lett. , 122301 (2015)Phys. Rev. Lett. , 122301 (2015).[8] D. R. Entem, R. Machleidt, and Y. Nosyk, “High-quality two-nucleon potentials up to fifth order of the chiral expansion,” Phys. Rev. C , 024004 (2017)Phys. Rev. C , 024004 (2017).[9] C. Drischler, K. Hebeler, and A. Schwenk, “Chiral interactionsup to next-to-next-to-next-to-leading order and nuclear satura-tion,” Phys. Rev. Lett. , 042501 (2019)Phys. Rev. Lett. , 042501 (2019).[10] J. Hoppe, C. Drischler, K. Hebeler, A. Schwenk, andJ. Simonis, “Probing chiral interactions up to next-to-next-to-next-to-leading order in medium-mass nuclei,”Phys. Rev. C , 024318 (2019)Phys. Rev. C , 024318 (2019).[11] T. Hüther, K. Vobig, K. Hebeler, R. Machleidt, andR. Roth, “Family of chiral two- plus three-nucleoninteractions for accurate nuclear structure studies,”Phys. Lett. B , 135651 (2020)Phys. Lett. B , 135651 (2020).[12] E. Epelbaum, H. Krebs, and P Reinert, “High-precision nu-clear forces from chiral EFT: State-of-the-art, challenges andoutlook,” Front. Phys. , 98 (2020)Front. Phys. , 98 (2020).[13] K. Hebeler, “Three-Nucleon Forces: Implementation and Ap-plications to Atomic Nuclei and Dense Matter,” (2020),arXiv:2002.09548arXiv:2002.09548.[14] W. G. Jiang, A. Ekström, C. Forssén, G. Hagen, G. R.Jansen, and T. Papenbrock, “Accurate bulk properties of nu-clei from A = ∞ from potentials with ∆ isobars,” (2020),arXiv:2006.16774arXiv:2006.16774. [15] W. H. Dickho ff and C. Barbieri, “Self-consistentGreen’s function method for nuclei and nuclear matter,”Prog. Part. Nucl. Phys. , 377 (2004)Prog. Part. Nucl. Phys. , 377 (2004).[16] G. Hagen, T. Papenbrock, M. Hjorth-Jensen, and D. J.Dean, “Coupled-cluster computations of atomic nuclei,”Rep. Prog. Phys. , 096302 (2014)Rep. Prog. Phys. , 096302 (2014).[17] H. Hergert, S. K. Bogner, T. D. Morris, A. Schwenk,and K. Tsukiyama, “The In-Medium Similarity Renormal-ization Group: A Novel Ab Initio Method for Nuclei,”Phys. Rept. , 165 (2016)Phys. Rept. , 165 (2016).[18] A. Tichai, R Roth, and T. Duguet, “Many-body perturbationtheories for finite nuclei,” Front. Phys. , 164 (2020)Front. Phys. , 164 (2020).[19] T. D. Morris, J. Simonis, S. R. Stroberg, C. Stumpf, G. Ha-gen, J. D. Holt, G. R. Jansen, T. Papenbrock, R. Roth,and A. Schwenk, “Structure of the lightest tin isotopes,”Phys. Rev. Lett. , 152503 (2018)Phys. Rev. Lett. , 152503 (2018).[20] J. D. Holt, J. Menéndez, J. Simonis, and A. Schwenk, “Three-nucleon forces and spectroscopy of neutron-rich calcium iso-topes,” Phys. Rev. C , 024312 (2014)Phys. Rev. C , 024312 (2014).[21] A. Tichai, J. Langhammer, S. Binder, and R. Roth, “Hartree-Fock many-body perturbation theory for nuclear ground-states,” Phys. Lett. B , 283 (2016)Phys. Lett. B , 283 (2016).[22] A. Tichai, P. Arthuis, T. Duguet, H. Hergert, V. Somá, andR. Roth, “Bogoliubov Many-Body Perturbation Theory forOpen-Shell Nuclei,” Phys. Lett. B , 195 (2018)Phys. Lett. B , 195 (2018).[23] S. Binder, J. Langhammer, A. Calci, and R. Roth, “Ab initiopath to heavy nuclei,” Phys. Lett. B , 119 (2014)Phys. Lett. B , 119 (2014).[24] K. Tsukiyama, S. K. Bogner, and A. Schwenk, “In-medium Similarity Renormalization Group for Nuclei,”Phys. Rev. Lett. , 222502 (2011)Phys. Rev. Lett. , 222502 (2011).[25] S. R. Stroberg, H. Hergert, S. K. Bogner, and J. D. Holt,“Nonempirical interactions for the nuclear shell model: An up-date,” Annu. Rev. Nucl. Part. Sci. , 307 (2019)Annu. Rev. Nucl. Part. Sci. , 307 (2019).[26] V. Somà, P. Navrátil, F. Raimondi, C. Barbieri, and T. Duguet,“Novel chiral Hamiltonian and observables in light andmedium-mass nuclei,” Phys. Rev. C , 014318 (2020)Phys. Rev. C , 014318 (2020).[27] T. A. Lähde, E. Epelbaum, H. Krebs, D. Lee, U.-G. Meißner,and G. Rupak, “Lattice E ff ective Field Theory for Medium-Mass Nuclei,” Phys. Lett. B , 110 (2014)Phys. Lett. B , 110 (2014).[28] G. Hagen, A. Ekström, C. Forssén, G. R. Jansen,W. Nazarewicz, T. Papenbrock, K. A. Wendt, S. Bacca,N. Barnea, B. Carlsson, C. Drischler, K. Hebeler, M. Hjorth-Jensen, M. Miorelli, G. Orlandini, A. Schwenk, and J. Simonis,“Neutron and weak-charge distributions of the Ca nucleus,”Nat. Phys. , 186 (2016)Nat. Phys. , 186 (2016).[29] P. Gysbers, G. Hagen, J. D. Holt, G. R. Jansen, T. D. Morris,P. Navrátil, T. Papenbrock, S. Quaglioni, A. Schwenk, S. R.Stroberg, and K. A. Wendt, “Discrepancy between experimen-tal and theoretical β -decay rates resolved from first principles,”Nat. Phys. , 428 (2019)Nat. Phys. , 428 (2019).[30] J. M. Yao, B. Bally, J. Engel, R. Wirth, T. R. Rodríguez,and H. Hergert, “Ab initio treatment of collective corre-lations and the neutrinoless double beta decay of Ca,”Phys. Rev. Lett. , 232501 (2020)Phys. Rev. Lett. , 232501 (2020).[31] M. A. Caprio, P. Maris, and J. P. Vary, “Coulomb-sturmian basis for the nuclear many-body problem,”Phys. Rev. C , 034312 (2012)Phys. Rev. C , 034312 (2012).[32] A. Tichai, J. Müller, K. Vobig, and R. Roth, “Natu-ral orbitals for ab initio no-core shell model calculations,”Phys. Rev. C , 034321 (2019)Phys. Rev. C , 034321 (2019).[33] P. Navratil, S. Quaglioni, I. Stetcu, and B. R. Barrett,“Recent developments in no-core shell-model calculations,”J. Phys. G , 083101 (2009)J. Phys. G , 083101 (2009). [34] B. R. Barrett, P. Navrátil, and J. P. Vary, “Ab initio no core shellmodel,” Prog. Part. Nucl. Phys. , 131 (2013)Prog. Part. Nucl. Phys. , 131 (2013).[35] G. Hagen, T. Papenbrock, and D. J. Dean, “Solution ofthe center-of-mass problem in nuclear structure calculations,”Phys. Rev. Lett. , 062503 (2009)Phys. Rev. Lett. , 062503 (2009).[36] S. J. Novario, G. Hagen, G. R. Jansen, and T. Papenbrock,“Charge radii of exotic neon and magnesium isotopes,” (2020),arXiv:2007.06684arXiv:2007.06684.[37] S. J. Novario, P. Gysbers, J. Engel, G. Hagen, G. R. Jansen,T. D. Morris, P. Navrátil, T. Papenbrock, and S. Quaglioni,“Coupled-cluster calculations of neutrinoless double-beta de-cay in Ca,” (2020), arXiv:2008.09696arXiv:2008.09696.[38] P. J. Hay, “On the calculation of natural orbitals by perturbationtheory,” J. Chem. Phys. , 2468 (1973)J. Chem. Phys. , 2468 (1973).[39] A. K. Q. Siu and E. F. Hayes, “Configuration interac-tion procedure based on the calculation of perturbationtheory natural orbitals: Applications to H and LiH,”J. Chem. Phys. , 37 (1974)J. Chem. Phys. , 37 (1974).[40] M. R. Strayer, W. H. Bassichis, and A. K. Kerman, “CorrelationE ff ects in Nuclear Densities,” Phys. Rev. C , 1269 (1973)Phys. Rev. C , 1269 (1973).[41] A. Szabo and N.S. Ostlund, Modern Quantum Chemistry: In-troduction to Advanced Electronic Structure Theory , DoverBooks on Chemistry (Dover Publications, 1989).[42] Isaiah Shavitt and Rodney J. Bartlett,
Many-Body Methods inChemistry and Physics: MBPT and Coupled-Cluster Theory ,Cambridge Molecular Science (Cambridge University Press,2009).[43] A. Carbone, A. Cipollone, C. Barbieri, A. Rios, and A. Polls,“Self-consistent Green’s functions formalism with three-bodyinteractions,” Phys. Rev. C , 054326 (2013)Phys. Rev. C , 054326 (2013).[44] C. Constantinou, M. A. Caprio, J. P. Vary, and P. Maris,“Natural orbital description of the halo nucleus He,”Nucl. Sci. Tech. , 179 (2017)Nucl. Sci. Tech. , 179 (2017).[45] S. B. Khadkikar and V. B. Kamble, “Centre of mass motionand single particle separation energies in Hartree-Fock theory,”Nucl. Phys. A , 352 (1974)Nucl. Phys. A , 352 (1974).[46] L. Jaqua, M. A. Hasan, J. P. Vary, and B. R. Barrett,“Kinetic-energy operator in the e ff ective shell-model interac-tion,” Phys. Rev. C , 2333 (1992)Phys. Rev. C , 2333 (1992).[47] H. Hergert and R. Roth, “Treatment of the IntrinsicHamiltonian in Particle-Number Nonconserving Theories,”Phys. Lett. B , 27 (2009)Phys. Lett. B , 27 (2009).[48] G. C. Wick, “The Evaluation of the Collision Matrix,”Phys. Rev. , 268 (1950)Phys. Rev. , 268 (1950).[49] S. K. Bogner, R. J. Furnstahl, and A. Schwenk,“From low-momentum interactions to nuclear structure,”Prog. Part. Nucl. Phys. , 94 (2010)Prog. Part. Nucl. Phys. , 94 (2010).[50] W. Kutzelnigg and D. Mukherjee, “Normal order and extendedWick theorem for a multiconfiguration reference wave func-tion,” J. Chem. Phys. , 432 (1997)J. Chem. Phys. , 432 (1997).[51] E. Gebrerufael, A. Calci, and R. Roth, “Open-shell nuclei andexcited states from multireference normal-ordered Hamiltoni-ans,” Phys. Rev. C , 031301(R) (2016)Phys. Rev. C , 031301(R) (2016).[52] H. Hergert, S. Binder, A. Calci, J. Langhammer, andR. Roth, “Ab Initio Calculations of Even Oxygen Iso-topes with Chiral Two-Plus-Three-Nucleon Interactions,”Phys. Rev. Lett. , 242501 (2013)Phys. Rev. Lett. , 242501 (2013).[53] E. Gebrerufael, K. Vobig, H. Hergert, and R. Roth, “Ab Ini-tio Description of Open-Shell Nuclei: Merging No-Core ShellModel and In-Medium Similarity Renormalization Group,”Phys. Rev. Lett. , 152503 (2017)Phys. Rev. Lett. , 152503 (2017).[54] G. Hagen, T. Papenbrock, D. J. Dean, A. Schwenk, A. Nogga,M. Włoch, and P. Piecuch, “Coupled-cluster theory for three-body Hamiltonians,” Phys. Rev. C , 034302 (2007)Phys. Rev. C , 034302 (2007). [55] R. Roth, S. Binder, K. Vobig, A. Calci, J. Langhammer, andP. Navrátil, “Medium-Mass Nuclei with Normal-Ordered Chi-ral NN +
3N Interactions,” Phys. Rev. Lett. , 052501 (2012)Phys. Rev. Lett. , 052501 (2012).[56] H. Hergert, “In-Medium Similarity Renormalization Group forClosed and Open-Shell Nuclei,” Phys. Scr. , 023002 (2017)Phys. Scr. , 023002 (2017).[57] S. R. Stroberg, A. Calci, H. Hergert, J. D. Holt,S. K. Bogner, R. Roth, and A. Schwenk, “Nucleus-dependent valence-space approach to nuclear structure,”Phys. Rev. Lett. , 032502 (2017)Phys. Rev. Lett. , 032502 (2017).[58] W. Magnus, “On the exponential solution ofdi ff erential equations for a linear operator,”Commun. Pure Appl. Math. , 649 (1954)Commun. Pure Appl. Math. , 649 (1954).[59] T. D. Morris, N. M. Parzuchowski, and S. K. Bogner, “Magnusexpansion and in-medium similarity renormalization group,”Phys. Rev. C , 034331 (2015)Phys. Rev. C , 034331 (2015).[60] R. J. Furnstahl, “The Renormalization Group in NuclearPhysics,” Nucl. Phys. Proc. Suppl. , 139 (2012)Nucl. Phys. Proc. Suppl. , 139 (2012).[61] S. K. Bogner, R. J. Furnstahl, S. Ramanan, and A. Schwenk,“Convergence of the Born series with low-momentum interac- tions,” Nucl. Phys. A , 203 (2006)Nucl. Phys. A , 203 (2006).[62] J. Hoppe, C. Drischler, R. J. Furnstahl, K. Hebeler, andA. Schwenk, “Weinberg eigenvalues for chiral nucleon-nucleoninteractions,” Phys. Rev. C , 054002 (2017)Phys. Rev. C , 054002 (2017).[63] J. Simonis, S. R. Stroberg, K. Hebeler, J. D. Holt, andA. Schwenk, “Saturation with chiral interactions and conse-quences for finite nuclei,” Phys. Rev. C , 014303 (2017)Phys. Rev. C , 014303 (2017).[64] M. S. Gordon, M. W. Schmidt, G. M. Chaban, K. R. Glaese-mann, W. J. Stevens, and C. Gonzalez, “A natural orbital di-agnostic for multiconfigurational character in correlated wavefunctions,” J. Chem. Phys. , 4199 (1999)J. Chem. Phys. , 4199 (1999).[65] K. Hebeler, “Momentum-space evolution of chiral three-nucleon forces,” Phys. Rev. C , 021002(R) (2012)Phys. Rev. C , 021002(R) (2012).[66] S. R. Stroberg, https: // github.com / ragnarstroberg //