Hartree-Fock Many-Body Perturbation Theory for Nuclear Ground-States

We investigate the order-by-order convergence behavior of many-body perturbation theory (MBPT) as a simple and efficient tool to approximate the ground-state energy of closed-shell nuclei. To address the convergence properties directly, we explore perturbative corrections up to 30th order and highlight the role of the partitioning for convergence. The use of a simple Hartree-Fock solution to construct the unperturbed basis leads to a convergent MBPT series for soft interactions, in contrast to, e.g., a harmonic oscillator basis. For larger model spaces and heavier nuclei, where a direct high-order MBPT calculation in not feasible, we perform third-order calculation and compare to advanced ab initio coupled-cluster calculations for the same interactions and model spaces. We demonstrate that third-order MBPT provides ground-state energies for nuclei up into tin isotopic chain that are in excellent agreement with the best available coupled-cluster results at a fraction of the computational cost.

Introduction. The solution of the Schrödinger equation for atomic nuclei using realistic nuclear interactions is at the heart of ab initio nuclear structure theory. In practice this problem is solved by constructing approximate methods for a truncated, i.e., finite-dimensional Hilbert space. However, for the calculation of ground-state energies of heavy nuclei significant algorithmic and computational efforts are needed. There exist a plethora of different ab initio methods, e.g., coupled cluster (CC) theory [1][2][3][4][5][6], in-medium similarity renormalization group (IM-SRG) [7][8][9][10][11], or self-consistent Green's function methods [12][13][14]. However, it is desirable to have an alternative, light-weight framework available. A conceptually simple method to solve for the eigenenergies of a physical system is many-body perturbation theory (MBPT) [15][16][17]. A perturbative treatment is the standard approach for many problems from different fields of theoretical physics. The advantage of MBPT compared to other ab initio approaches is its simplicity, which also allows for straightforward generalizations to excited states and open-shell nuclei [18] without the need of sophisticated equation-of-motion techniques. The reason why MBPT usually is not considered as ab initio technique are convergence issues of the underlying perturbation series. Several studies of high-order MBPT based on Slater determinants constructed from harmonic oscillator (HO) single-particle states (HO-MBPT) have shown that the perturbation series is divergent in almost every case [18,19]. In such cases one heavily relies on the use of resummation techniques, e.g., Padé approximants, that enable a robust extraction of observables although the perturbative expansion diverges [19][20][21].
In this Letter, we formulate MBPT based on Hartree-Fock (HF) single-particle states (HF-MBPT), and, for the first time, investigate the convergence behavior of the perturbation series up to 30 th order. We compare the ground-state energies of 4 He and 16 O to results from exact diagonalizations in the configuration interaction (CI) approach using the same model space [22][23][24]. Based on the rapidly converging perturbation series resulting from the use of HF basis states, we study ground-state energies of selected closed-shell medium-mass and heavy nuclei at third-order MBPT, and compare to recent coupled cluster (CC) calculations [6].
The Nuclear Hamiltonian. For all following investigations we start from the chiral nucleon-nucleon (NN) interaction at next-to-next-to-next-to leading order (N 3 LO) by Entem and Machleidt [25] combined with the three-nucleon (3N) interaction at next-to-next-to leading order N 2 LO in its local form [26] with three-body cutoff Λ 3N = 400 MeV/c. Additionally, we use the similarity renormalization group (SRG) to soften the Hamiltonian through a continuous unitary transformation controlled by a flow parameter [27][28][29][30][31]. In principle this transformation induces beyond-3N operators up the mass number of the considered nucleus, which, however, we have to neglect. To avoid the complication of dealing with explicit 3N interactions, we make use of the normal-ordered two-body approximation (NO2B) of the 3N interaction that has been found to be very accurate for medium-mass nuclei, see Refs. [32,33]. For the matrix-element preparation we adopt the procedure we introduced in Ref. [6], i.e., in particular, we use the large SRG model-space and exploit the iterative scheme where necessary. Thus, the matrix elements and the treatment of the chiral NN+3N interaction are identical to Ref. [6] and we can compared directly to the CC results presented there.
Many-Body Perturbation Theory. The essence of Rayleigh-Schrödinger perturbation theory is the definition of an additive splitting, referred to as partitioning, of a given Hamiltonian H into an unperturbed part H 0 and a perturbation W. Introducing an auxiliary parameter λ yields a one-parameter family of operators, where the perturbation is defined by W = H − H 0 . As ansatz for the solution of the eigenvalue problem of H we take a power series expansion in terms of an auxiliary parameter λ, where the expansion coefficients are given by the energy corrections and state corrections, respectively. We choose H 0 to be the HF Hamiltonian arising from an initial NN+3N interaction. We have shown in Refs. [18,19]  allowing for detailed investigations of the convergence characteristics of the perturbation series. In general we cannot expect that a perturbation series is convergent [34][35][36], but one can exploit resummation-theory techniques to extract information on the observables of interest. There are different schemes and transformations that can be used to extract, e.g., the ground-state energy from a divergent expansion [37][38][39]. Padé approximants have proven to be particularly useful in the treatment of high-order HO-MBPT [18,19]. Additionally, they are well-known to mathematicians especially in the field of convergence acceleration [21,36,37]. However, the calculation of energy corrections up to sufficiently high orders is only feasible for light nuclei due to increasing computational requirements. When proceeding to the medium-mass region one must choose a different strategy. Depending on the rate of convergence, one might expect low-order partial sums of the perturbation series to be a reasonable approximation to the exact ground-state energy. Having only low-order information available, resummation methods are less effective, because one can only construct a small number of approximants that yield valid approximations only if the transformed sequence converges sufficiently fast [18]. However, one alternative is to exploit the freedom in the partitioning, i.e., the choice of the unperturbed basis, to improve the convergence of the perturbation series. For low-order calculations the corrections can be expressed in terms of the particle-hole formalism. Note that the Hartree-Fock energy corresponds to the the first-order partial sum, Therefore, the first contribution to the correlation energy appears in second-order HF-MBPT. The second-and third-order partial sum for the ground-state energy with respect to HF basis for a two-body operator are given by [40] In the third-order energy correction the first, second and third term are called particle-particle (pp), hole-hole (hh) and particle-hole (ph) correction, respectively. The i correspond to the HF single-particle energies and all matrix elements are taken to be antisymmetrized. Summation indices a, b, c, .. correspond to particle indices, i.e., above the Fermi level F , whereas i, j, k, ... correspond to hole indices ranging from lowest occupied single-particle states up to the Fermi level. The zero-and one-body parts of the normal-ordered Hamiltonian only enter in the first-order energy correction. Brillouin's theorem states that there is no mixing of the Hartree-Fock ground state with singly-excited determinants [17] and by orthogonality the zero-body part is only present in the expectation value of the perturbation. In principle, the derivation of energy corrections beyond third order is straight forward. However, considering a diagrammatic approach in terms of Hugenholtz diagrams, the number of contributing diagrams at a given perturbation order p increases rapidly [41] such that it becomes challenging to go beyond third-order in practice. Additionally, terms from higher-order corrections involve expressions that are notoriously hard to compute, because their effective implementation, e.g., by means of BLAS-enabled matrix operations, is not obvious. The computational power needed to perform third-order MBPT calculations up to the heavy-mass region can in principle be provided by a single computing node within 1 − 3% of the computing time needed for state-of-theart CC calculations.
Convergence Characteristics of Hartree-Fock Many-Body Perturbation Theory. We start with comparing perturbation series from HO-and HF-MBPT, and focus on their convergence characteristics and sensitivity to the SRG flow parameter. In Fig. 1 we present a direct comparison of the orderby-order behavior for the two partitionings up to 30th order for 16 O. For these high-order calculations we use an N maxtruncation of the many-body model space similar the no-core shell model (NCSM) [23]. The left column of Fig. 1 shows the high-order partial sums and the right column the individual energy corrections for each order. Panel (a) shows the partial sums from HO-MBPT for a sequence of model spaces with fixed SRG flow parameter α = 0.08 fm 4 . The partial sums are divergent for every model space. The divergence is also apparent from panel (c) which reveals the exponentially increasing energy corrections. In contrast, panel (b) shows the partial sums arising from HF-MBPT that are convergent for all model spaces. Furthermore, the converged values agree with direct CI results. As seen in panel (d), the energy corrections are exponentially suppressed for higher orders, giving rise to a robust convergence.
In Fig. 2 we show the high-order partial sums and energy corrections in HF-MBPT for different SRG flow parameters. Panels (a), (b) and (c) show the convergent perturbation series for 4 He, 16 O and 24 O, respectively. The calculations are performed for fixed N max = 6 for 4 He, 16 O and N max = 4 for 24 O and the flow-parameter dependence of the absolute energies results from the varying degree of convergence with respect to the many-body model space.
A more interesting flow-parameter dependence can be observed for the individual energy corrections in panels (d), (e), and (f). There is a clear dependence of the convergence rate on the flow parameter for the oxygen isotopes. For 16 O the series converges exponentially in all three cases and the larger the flow parameter, i.e. the softer the Hamiltonian, the more rapid the convergence-as might be naively expected. For 24 O the behavior is slightly more complicated. For the softest interaction with α = 0.08fm 4 there is still a clear exponential convergence. However, for the harder interactions, i.e., α = 0.02, 0.04fm 4 , we observe no systematic decrease of the high-order perturbative contributions anymore, they remain approximately constant and cause a small-amplitude os- cillatory behavior of the partial sums. However, even in these cases we can easily extract a robust estimate for the asymptotic value. In the case of 4 He the suppression is independent of α and we observe the same rapid convergence for all interactions.
The numerical values of the partial sums for selected orders of HF-MPBT for the three nuclei and the different flow parameters are summarized in Tab. I together with the results of direct CI calculations for the same Hamiltonians and model spaces. The higher-order partial sums are in good agreement with the CI results-in most cases the deviation of the groundstate energy is much smaller than 0.1%.
Based on our detailed analysis of high-order HF-MBPT and due to the exponential suppression of the energy corrections, we can take low-order partial sums as a reasonable approximation to the converged results. This motivates the investigation of third-order partial sums for selected medium-mass and heavy closed-shell nuclei in the following.
Explicit Summation for Heavy Nuclei. For heavier nuclei and larger model spaces we cannot compute the high-order perturbation series explicitly and, thus, we cannot investigate the convergence characteristics explicitly. We can, however, evaluate the perturbative contributions up to third order very efficiently. To demonstrate the validity of a low-order perturbative approximation, we need to compare our results to established ab initio techniques, in our case, coupled-cluster calculations with sophisticated triples corrections.
We consider a sequence of closed-shell nuclei ranging from 4 He to 132 Sn and perform calculations in second and thirdorder HF-MBPT in a large model space truncated with respect to the single-particle principal quantum number e max = 12. We restrict ourselves to SRG-evolved Hamiltonians with flow parameter α = 0.08 fm 4 , which was used extensively in previous calculations and showed favorable order-by-order convergence in our high-order studies. We cannot perform CI cal- culations for these large spaces, however, the coupled-cluster framework has proven to provide accurate results for groundstate energies of closed-shell nuclei [1][2][3][4]. We compare the HF-MBPT results to recent CC calculations at the CCSD and the CR-CC(2, 3) level [5,6,42]. Starting from a HF reference state this approach provides a complete inclusion of singly and doubly excited clusters on top of the reference state and, in the case of CR-CC(2, 3) an approximate non-iterative inclusion of triply excited clusters [43][44][45][46]. In Figs. 3 and 4 the ground-state energies per nucleon (a) as well as the correlation energy E corr = E − E HF per nucleon (b) from HF-MBPT and CR-CC(2, 3) are depicted for an initial chiral NN+3N and an initial chiral NN interaction. The SRGinduced three-nucleon contribution are taken into account in both cases, leading to the NN+3N-full and NN+3N-induced interactions, respectively.
These figures show a remarkable result: The binding energies in third-order HF-MBPT and CR-CC(2,3) are in excellent agreement with each other. The relative differences are in most cases much smaller than 1%. The same observation holds for the correlation energy, i.e., the corrections to the HF energy. The third-order energy corrections contribute approximately 0.2 MeV to the overall binding energy per nucleon and are, therefore, non negligible even though the third-order energy corrections in HF-MBPT are one order of magnitude smaller than the second-order correction. The third-order energy contribution (3) consists of three terms corresponding to three Hugenholtz diagrams. Figure 5 disentangles their individual contributions to the overall thirdorder energy correction. The contribution of the pp and hh corrections are almost constant over the entire mass range, whereas the energy correction arising from the ph term scales with increasing mass number in the case of a NN+3N-full interaction. For the tin isotopes the third-order energy correction contributes 3% of the overall binding energy in thirdorder HF-MBPT and is not negligible. In particular we see that most of the third-order energy correction arises from the ph diagram. In the case of a NN+3N-induced interaction we see that all three terms are suppressed with increasing mass number. These systematic dependencies of the individ-ual third-order contributions on the input Hamiltonian show that a partial inclusion of selected third-order terms may lead to wrong estimates.

Conclusions.
We have discussed Rayleigh-Schrödinger MBPT as an efficient approach to compute ground-state energies for closed-shell nuclei up to the heavy-mass region. The use of a HF basis has enabled us to overcome convergence problems that generally arise in HO-MBPT. Investigating 16 O in different model spaces showed convergent partial sums when using HF-MBPT and the limit of the perturbation series coincides with the results from explicit CI calculations. Additionally, we found systematic dependencies of the convergence rate on the SRG parameter in the case of 16 O and 24 O. Thus, in HF-MBPT we can improve the conver-  gence behavior of the perturbation series by further evolving the Hamiltonian, whereas the divergent HO-MBPT series are unaffected.
We can identify a hierarchy of elements influencing the convergence properties of perturbation series. Defining a par-titioning, or equivalently, defining a starting point of the recursive calculation is the most important part. We have seen from the radically different behavior of the perturbation series in HF-MBPT and HO-MBPT that the partial sums are very sensitive to the partitioning. When using HF-MBPT we can improve the order-by-order convergence by using softer interactions corresponding, e.g., to larger SRG flow parameters. Even for HF basis sets harder interactions can spoil convergence. The 'softness' of the interaction has been characterized in terms of Weinberg eigenvalues, which are connected to the spectrum of two-body Green's functions [47][48][49]. Similar expressions also appear in the formulas for the first-order state correction. Though the general connection seems obvious, one should be careful with conclusions about the convergence of MBPT for a finite nucleus based on the softness of the interaction. Our work has shown that the partitioning is key for convergence. Our observation that the convergence of HF-MBPT deteriorates for harder interactions could simply be explained by the fact that the unperturbed HF solution becomes a much worse approximation for the ground state in these cases.
The superior convergence properties of HF-MBPT has motivated the use of low-order approximations to investigate nuclei in the medium-mass region. We have validated these loworder approximations to the most sophisticated CC calculations and found excellent agreement of third-order HF-MBPT and CR-CC(2,3) at the level of better than 1%. The consistency of high-order partial summations with exact CI diagonalizations as well as the agreement of low-order summations with coupled-cluster results may qualify HF-MBPT as an ab initio approach. However, the strong dependence of convergence on the partitioning should be a reason for caution. The HF partitioning seems to be robust for sufficiently soft interactions, but there is no formal guarantee for convergence.
The great advantage of low-order HF-MBPT is its simplicity: Computationally, the third-order calculations are much cheaper than CC or IM-SRG calculations. It is, therefore, ideal for survey calculations over a large range of mediummass nuclei, e.g., to explore the ground-state systematics for new interactions. Formally, the underlying equations and algorithms are trivial compared to CC or IM-SRG. As a result of this formal simplicity, extensions to the description exited states and open-shell nuclei are straight-forward. We have demonstrated this already for light nuclei using high-order degenerate HO-MBPT [18]. Alternative multi-configurational formulations for open-shell nuclei are under investigation.