Pade-resummed high-order perturbation theory for nuclear structure calculations

We apply high-order many-body perturbation theory for the calculation of ground-state energies of closed-shell nuclei using realistic nuclear interactions. Using a simple recursive formulation, we compute the perturbative energy contributions up to 30th order and compare to exact no-core shell model calculations for the same model space and Hamiltonian. Generally, finite partial sums of this perturbation series do not show convergence with increasing order, but tend to diverge exponentially. Nevertheless, through a simple resummation via Pade approximants it is possible to extract rapidly converging and highly accurate results for the ground state energy.


Introduction
The treatment of the nuclear many-body problem is a central and long-standing issue in nuclear structure theory. Ideally, we would like to solve the many-body problem ab initio, i.e., starting from a given nuclear Hamiltonian without any conceptual approximations. With the advent of high-precision nuclear potentials that are based systematically on Quantum Chromodynamics (QCD) through chiral effective field theory [1,2], the demand for exact ab initio solutions of the nuclear many-body problem has grown. Only these schemes establish a rigorous and quantitative connection between nuclear structure observables and the underlying QCD input.
The no-core shell model (NCSM) is one of the most universal exact ab initio methods, which gives access to all aspects of nuclear structure [3,4,5]. Other methods, are either restricted to certain classes of Hamiltonians, like the Green's Function Monte Carlo approach [6], or they are limited to certain nuclei and observables, like the coupled-cluster approach [7]. All of them are computationally demanding, which leads to a severe limitation regarding the number of nucleons that can be handled.
Therefore, approximate many-body schemes using the same Hamiltonians, i.e. approximate ab initio methods, also provide indispensable information. In particular approaches that use controlled and systematically improvable approximations are of great practical importance. In this category, many-body perturbation theory (MBPT) is one of the most powerful and widely used methods. On the one hand, the evaluation of loworders of perturbation theory is computationally simple and can be done for the whole nuclear mass range [8,9,10,11] as well as for infinite nuclear matter [12]. On the other hand, it is deemed systematically improvable, either by extending the MBPT calculations order-by-order or by using infinite partial summations, like ladder-or ring-type summations [13,14,15]. However, the accuracy of low-order perturbative estimates, e.g. for ground-state energies, or possible extensions of the MBPT series to higher orders and the resulting convergence pattern are rarely, if ever, addressed in the nuclear structure context.
In this paper, we apply MBPT for the calculation of the ground state energy of several closed-shell nuclei. We extend the order-by-order calculation of the perturbative energy contributions up to 30th order, study the convergence behavior, and compare to exact NCSM calculations for the same Hamiltonian and model space. We introduce Padé approximants as a highly efficient tool for the resummation of the divergent power-series of MBPT into a rapidly converging series and demonstrate their accuracy for the description of ground-state energies.

Formalism
We aim at a perturbative expansion of the many-nucleon Schrödinger equation for the translational invariant nuclear Hamiltonian H = T − T cm + V, where we assume V to be a two-body interaction for simplicity. In a first step we have to chose the unperturbed basis, which in turn defines the unperturbed Hamiltonian. From the practical point of view, a basis of Slater-determinants constructed from a set of single-particle states is most convenient. The underlying single-particle basis will typically be a Hartree-Fock or a harmonic oscillator basis-for simplicity we assume the latter. The unperturbed Hamiltonian H 0 is a one-body operator containing the kinetic energy T and a harmonic oscillator potential. The unperturbed Slater determinants |Φ n fulfill the eigenvalue relation with eigenvalues ǫ n being the sum of the single-particle energies of the occupied states. After the unperturbed Hamiltonian is fixed, the perturbation is defined through W = H − H 0 . This partitioning leads to the Møller-Plesset formulation of MBPT and obviously other partitionings of the Hamiltonian are possible [16,17]. For ease of presentation, we assume that the unperturbed state corresponding to the eigenstate we are interested in is non-degenerate, as it is the case for the ground state of closed shell nuclei. In the case of degeneracy, as e.g. for the excited states of closed shell nuclei, one would have to diagonalize the full Hamiltonian in the degenerate subspace and pick the eigenstates with the desired quantum numbers as unperturbed states. The standard Rayleigh-Schrödinger perturbation series can now be constructed based on a Hamiltonian (using the notation from Ref. [17]) containing an auxiliary expansion parameter λ that continuously connects the unperturbed Hamiltonian H 0 = H(λ = 0) with the full Hamiltonian H = H(λ = 1). The energy eigenvalues E n (λ) and the corresponding eigenvectors |Ψ n (λ) of H(λ) are formulated as a power series in λ In the absence of degeneracy the lowest-order contributions are simply given by the unperturbed quantities, i.e., Inserting the Hamiltonian (3) and the power series (4) into the Schrödinger equation (1) leads to the fundamental equation Assuming that the unperturbed states form an orthonormal basis and using the intermediate normalization n |Ψ (p) n = 0 for p > 0, which allows us to project-out all required information on the individual contributions in the power series. By multiplying Eq. (6) with Ψ (0) n | and matching same orders of λ on both sides, we immediately obtain a simple expression for the pth-order energy contribution By multiplying Eq. (6) with Ψ (0) m | with m n and matching λ-orders, we obtain an expression for the amplitudes which characterize the perturbative corrections to the eigenstates |Ψ (p) n expanded in the unperturbed basis with C (p) n,n = 0 for p > 0 and C (0) n,m = δ n,m . We can cast Eqs. (7) and (8) into a more transparent form by systematically introducing the amplitudes C (p) n,m and formulating all matrix elements in terms of the unperturbed states. For the pth-order energy contribution we obtain Similarly we obtain for the pth-order amplitudes Together with C (0) n,m = δ n,m and E (0) n = ǫ n these relations form a recursive set of equations which uniquely determines the perturbative corrections for all energies and states to all orders.
Usually one would use these general expressions to derive explicit formulae for the lowest-order corrections. The matrix elements of the perturbation in the unperturbed Slaterdeterminant states can be evaluated explicitly and the summations over the many-body basis set can be replaced by summations over single-particle states. In this way we would recover the standard expressions for, e.g., the second-and third-order energy corrections [8,11,16].

Evaluation to high orders
When attempting to evaluate the perturbative corrections beyond third-or forth-order the explicit formulae for the energy corrections become impractical because of the large number of nested summations. A much more elegant way to evaluate highorder contributions makes use of the recursive structure of Eqs. (10) and (11). The only ingredients needed are the many-body matrix elements of the full Hamiltonian H with respect to the unperturbed basis |Φ n . Starting from the zeroth-order coefficients C (0) n,m = δ n,m we can readily evaluate the first-order energy contribution E (1) n from (10). This in turn allows us to compute the first-order coefficients C (1) n,m via (11). Generally, for the evaluation of the energy contribution E (p) n only the coefficients n,m of the previous order are required. For the evaluation of the coefficients C (p) n,m all energy contributions up to order p and all coefficients up to order (p − 1) need to be known.
Technically, the recursive evaluation of the perturbation series bears some resemblance to the Lanczos algorithm for the iterative solution of the eigenvalue problem for a few extremal eigenvalues as it is used in the NCSM. The most significant operation is a matrix-vector multiplication of the Hamiltonian matrix with the coefficient vector from the previous order, which constitutes the first term in the evaluation of the coefficients (11). Because the second term in (11) involves the coefficient vectors from all previous orders, we store them for simplicity. These computational elements are the same as for a simple Lanczos algorithm in the NCSM or in corresponding configuration interaction (CI) approaches, therefore, an implementation of high-order MBPT using NCSM or CI technologies is straight forward. However, since the computational elements are the same, so are the computational limitations: This direct implementation of high-order MBPT is limited to the same model spaces and particle numbers as the full NCSM. This is not a concern for the present study, but for an application of MBPT beyond the domain of the NCSM one has to resort to other evaluation schemes.

Applications: 4 He, 16 O and 40 Ca
As examples for a direct application of high-order MBPT we consider the ground-state energies of 4 He, 16 O, and 40 Ca. Throughout this work we use an intrinsic Hamiltonian with a soft two-nucleon interaction that is derived from the chiral N3LO potential [1] via a Similarity Renormalization Group (SRG) transformation [18,19,20]. The final flow parameter for the SRG evolution of the interaction is α = 0.02 fm 4 which corresponds to a momentum scale of Λ = 2.66 fm −1 . This choice for the flow parameters leads to a unitarily transformed interaction which is sufficiently soft to warrant excellent convergence properties with respect to model space size in the NCSM but at the same time yields ground-state energies which are in reasonable correspondence with experiment in the mass-range considered here. To allow for a direct comparison with exact NCSM calculations for the same Hamiltonian and the same model space, we use an N max Ω model space also for the MBPT calculations. We have confirmed, however, that all conclusions regarding the performance and limitations of the MBPT approach do not depend on this particular choice.
In Fig. 1 we summarize the results of an order-by-order MBPT calculation up to 30th order for the ground state energy of the three nuclei. For 4 He the calculations were performed in a 12 Ω model space, for 16 O in 6 Ω, and for 40 Ca in 4 Ω, each with two different oscillator frequencies Ω. The partial sum of the perturbative energy contributions up to order p, is depicted in upper row and the modulus of the individual pthorder contributions, |E (p) |, on a logarithmic scale in the lower row. Here and in the following we omit the index n = 0 for convenience. Already the first glance at Fig. 1 reveals a fundamental problem with the convergence behavior of the perturbation series. For 4 He, as depicted in Fig. 1(a), we observe two different patterns depending on the oscillator frequency. For Ω = 20 MeV the partial sum E sum (p) shows an alternating behavior with a systematically decreasing amplitude. Beyond 10th order one might consider the perturbation series converged and the resulting energy is in excellent agreement with the result of an NCSM calculation with the same Hamiltonian in the same model space. However, a change of the oscillator frequency destroys this picture. For Ω = 32 MeV, where the NCSM provides a lower ground-state energy, we again observe an alternating sequence of energy contributions E (p) , but this time without any sign of convergence. The absolute value of the individual energy corrections does not decrease with increasing order, it even shows a slightly increasing trend. Hence even in the simplest case, the 4 He ground state, the convergence of the perturbation series is not guaranteed.
The situation is even more dramatic for 16 O or 40 Ca as depicted in Fig. 1(b) and (c). In all cases the size of the perturbative energy contributions |E (p) | grows exponentially with p. The partial sum E sum (p) exhibits a strong oscillatory behavior with increasing amplitude. Only the lower orders, typically up to 10th order for 16 O and up to 5th order for 40 Ca, lead to binding energies in a physically meaningful energy range. At the 30th order the perturbative contributions are in the order of 10 4 MeV for 16 O and 10 9 MeV for 40 Ca-this is beyond any physical energy scale present in the nuclear many-body problem. We were not able to find a convergent scenario by varying the oscillator frequency or the model space size or truncation for these nuclei.
The explosion of the perturbative corrections beyond any meaningful energy scale suggests a principal mathematical problem in the representation of the energy eigenvalue E(λ) as a partial sum of a simple power series (4) .

Formalism
Prompted by the drastic failure of a partial sum of a simple power series to describe the energy E(λ) at the physical point λ = 1 one might consider more general expansions of this function. A next step would be an expansion of the energy E(λ) in terms of a rational function composed of separate power series for numerator and denominator Obviously we will not attempt to re-derive perturbation theory for this type of expansion. The above is useful only, if we could use the information contained in the standard MBPT energy contributions E (p) to construct this rational expansion. Exactly this is achieved through the Padé approximants [21,22]. Given a power series (4) of the function E(λ), then the Padé approximant with numerator being a polynomial of order M and the denominator a polynomial of order N is constructed such that its Taylor expansion reproduces the first M + N orders of the initial power series, i.e.
From this definition one can immediately construct a coupled system of equations that determines the coefficients a n and b m of the Padé approximant from a given set of E (p) with p = 0, ..., N + M. An alternative and more elegant form [21,22] relates the Padé approximants to determinants of two (N + 1) × (N + 1) matrices containing directly the power-series coefficients E (p) where we set E (p) ≡ 0 for p < 0. We will use this form to evaluate various Padé approximants in the following. Before considering numerical results, we should like to mention a few formal properties of the Padé approximants that are of importance for the present application. The mathematical foundation for using Padé approximants for our purpose in the first place is provided by the Padé conjecture (simplified) [21,22]: Let E(λ) be a continuous function for |λ| ≤ 1, then there is an infinite subsequence of diagonal Padé approximants [N/N](λ) that for N → ∞ converges locally uniformly against E(λ) for |λ| ≤ 1. For our application the continuity requirements for the function E(λ) are always fulfilled, thus we expect the diagonal Padé approximants to show a convergence behavior-unlike the simple power series.
Additionally the Padé approximants have a number of specific properties that would be extremely valuable in the present context. If the power series expansion of E(λ) is a Stieltjes series, then the Padé approximants fulfill the condition for λ ≥ 0 as well as a whole set of related inequalities [21,22]. Thus the diagonal and the super-diagonal Padé approximants provide upper and lower bounds for the full energy E(λ), respectively. Furthermore, these bounds improve monotonically with increasing order M of the Padé approximant. Unfortunately, it turns out that the power series we start from is not a Stieltjes series in general.

Applications: 4 He, 16 O and 40 Ca
Using the results of the order-by-order calculation of the energy corrections E (p) up to 30th order of MBPT we can construct all Padé approximants with N + M ≤ 30 from Eq. (16). Evaluating the approximant at λ = 1 yields an estimate for the ground-state energy of the perturbed system A collection of all diagonal as well as sub-and superdiagonal approximants with N + M ≤ 30 for 4 He, 16 O, and 40 Ca using the oscillator frequencies that yield the lowest groundstate energy is provided in Fig. 2. The first remarkable observation is that the Padé approximants converge very quickly for sufficiently large order-we have observed this behavior in all cases we considered. For M + N 10 essentially all Padé approximants provide the same ground-state energy. We emphasize that the input for the construction of those Padé approximants are the exponentially diverging coefficients from the power-series formulation of MBPT discussed in Fig. 1. The Padé resummation of these coefficients efficiently regularizes these divergencies and leads to exceptionally stable results for all orders M + N 10. The second important observation results from the comparison of the converged Padé approximants 1 with the exact energy eigenvalue obtained from a solution of the matrix eigenvalue problem for the Hamiltonian in the same model space-i.e., from the corresponding NCSM calculation-as indicated by the dashed horizontal line in Fig. 2. The converged Padé approximants exactly reproduce the corresponding energy eigenvalues, i.e., Padé resummed perturbation theory and the exact solution of the eigenvalue problem become equivalent.
A quantitative comparison is presented in Tab. 1, where the difference of the partial sums E sum (p) and the Padé approximants E Padé (M/N) to the exact NCSM eigenvalues E NCSM are shown. The latter were obtained using the Antoine shell-model code [23]. Starting from M + N ≈ 10 the deviations of the Padé approximants from the exact result are getting very small and starting from M + N ≈ 20 the Padé approximants are numerically identical to the exact result for all nuclei. In this regime the individual MBPT contributions E (p) are already increasing exponentially for 16 O and 40 Ca (cf. Fig. 1) and the partial sum E sum (p) does not provide any sensible estimate of the groundstate energy. The Padé approximants prove to be a highly efficient tool to extract a virtually exact and stable result for the energy from the first 10 or more coefficients E (p) of the strongly fluctuating and non-converging power series. Considering the scale of the order-to-order fluctuations and the absolute size of the perturbative contributions E (p) the stability and the precision of the converged Padé approximants is truly remarkable.
For application purposes, the behavior at low orders is also of interest. As shown in Tab. 1, the deviations ∆E sum (p) and ∆E Padé (M/N) are of comparable magnitude up to about fifth order, both showing sizable fluctuations. Hence, in this low-order domain the Padé approximants do not improve on the results obtained from a simple partial sum. Only beyond M + N ≈ 5 do the Padé approximants start to converge, i.e., the variations within a set of approximants of neighboring order reduce systematically. At the same time the deviations of the partial sums, ∆E sum (p), start to increase exponentially.
The stability of the Padé approximants E Padé (M/N) across various neighboring orders M and N is an important intrinsic criterion for convergence and for the accuracy of the Padé approximants as compared to the exact result. Therefore, it seems advisable to always consider sets of several approximants. Moreover, there is always the possibility that individual approximants completely escape the general trend, such as the E Padé (2/2) approximant for 40 Ca at Ω = 20 MeV that has a large positive and thus unphysical value. These cases are a reminder that the convergence theorems for Padé approximants, e.g., the Padé conjecture, only cover subsequences of approximants. Finally we note that the MBPT power series in the present examples turns out not to be a Stieltjes series. As the Padé approximants of Tab. 1 show, the inequality (17) as well as related inequalities are not fulfilled. The Padé approximants for the ground state energy on nuclei in the present MBPT framework do not provide rigorous bounds for the exact eigenvalues.

Conclusions & Outlook
We have formulated and applied many-body perturbation theory up to high orders for the description of ground-state energies of closed-shell nuclei using realistic Hamiltonians. In contrast to typical applications of MBPT in nuclear physics that are limited to second or third order, we extend the order-byorder evaluation of the perturbation series up to 30th order using a simple recursive scheme. In order to facilitate the comparison with exact eigenvalues obtained in NCSM calculations, we have limited ourselves to an harmonic-oscillator single-particle basis and a N max Ω space. However, results for other single-particle bases and model-space truncations are qualitatively similar.
Our major results and conclusions are: First, a simple partial sum of the perturbative series in general does not converge. On the contrary, we typically observe an exponential increase of the individual perturbative contributions |E (p) | and an oscilla- tory behavior of the partial sum E sum (p) as function of p. Thus, finite partial summations, even if they are extended to high orders, do not provide a stable and systematically improvable approximation for the exact energy eigenvalue. Second, Padé approximants offer a computationally simple yet powerful tool to extract a convergent series from a finite set of perturbative energy contributions E (p) . Solely through a resummation of the finite pth-order power series to a rational function, whose Taylor expansion up to order p is identical to the initial power series, we are able to extract a highly stable and convergent approximation for the energy. The information entering these Padé approximants of order p = M + N is identical to the power series of order p and so is the computational effort. However, whereas a finite partial sum E sum (p) explodes with increasing order p, the Padé approximants E Padé (M/N) converge.
Third, beyond a sufficiently large order, typically M + N 10, the different Padé approximants become very stable and converge rapidly to a unique value for the energy. This energy is identical to the exact eigenvalue obtained in the corresponding NCSM calculation. In this sense, Padé resummed MBPT and the Lanczos diagonalization become numerically equivalent. Unfortunately, the computational effort is also comparable at least for the implementation we adopted for the MBPT here.
Fourth, at low orders, i.e. M + N 4, the Padé approx-imants generally do not yield an improved approximation for the energy. Their deviations from the exact eigenvalue fluctuate and are of the same order of magnitude as the errors of the corresponding partial sums. If one is limited, for computational reasons, to very low orders, then the version of MBPT used here can only provide a rough estimate that might differ significantly from the exact eigenvalue.
These initial studies open a number of new avenues for the study and application of MBPT in nuclear structure. Beyond the MBPT calculations presented here, one can use optimized single-particle bases and improved partitionings of the Hamiltonian to influence the convergence behavior of a finite order-byorder MBPT calculation. Furthermore, one can exploit infinite partial summations, e.g., ladder-or ring-type summations including MBPT contributions from all orders, and compare their quality to Padé resummed finite-order calculations. These improvements and extensions will be the subject of future studies.