Converging Many-Body Perturbation Theory for Ab Initio Nuclear-Structure: I. Brillouin-Wigner Perturbation Series for Closed-Shell Nuclei

Convergence aspects of nuclear many-body perturbation theory for ground states of closed-shell nuclei are explored using a Brillouin-Wigner formulation with a new vertex function enabling high-order calculations. A general formalism for Hamiltonian partitioning and a convergence criterion for the perturbation series are proposed. Analytical derivation shows that with optimal partitionings, the convergence criterion for ground states can always be satisﬁed. This feature attributes to the variational principle and does not depend on the choice of an internucleon interaction or a many-body basis. Numerical calculations of the ground state energies of 4 He and 16 O with Daejeon16 and a bare N 3 LO potential in both harmonic-oscillator and Hartree-Fock bases conﬁrm this ﬁnding.

Convergence aspects of nuclear many-body perturbation theory for ground states of closed-shell nuclei are explored using a Brillouin-Wigner formulation with a new vertex function enabling highorder calculations.A general formalism for Hamiltonian partitioning and a convergence criterion for the perturbation series are proposed.Analytical derivation shows that with optimal partitionings, the convergence criterion for ground states can always be satisfied.This feature attributes to the variational principle and does not depend on the choice of an internucleon interaction or a manybody basis.Numerical calculations of the ground state energies of 4 He and 16 O with Daejeon16 and a bare N 3 LO potential in both harmonic-oscillator and Hartree-Fock bases confirm this finding.
It is from the advent of quantum mechanics that theoretical description of many-fermion systems became a long-standing challenge.Originated in seminal works by Brueckner [1,2], Bethe [3,4] and Goldstone [5], on nuclear systems, many-body perturbation theory (MBPT) became the method of choice also in atomic and molecular physics [6].The majority of applications of MBPT in nuclear physics were of the Rayleigh-Schrödinger (RS) type via a diagrammatic representation [5,7] of the ground-state energy or of the effective interaction for the nuclear shell model [8][9][10][11][12][13].In those earlier studies, the so-called Brueckner G-matrix [4,[14][15][16] was used to deal with the strong short-range repulsion of the nuclear force.However, poor convergence of intermediatestate summations due to the strong tensor force [17] and weak evidence for order-by-order convergence of up to the 4th order expansions in G [18,19] inhibited extensive applications of MBPT in nuclear physics.Arrival of a new generation of internucleon potentials consistently obtained within chiral effective field theory [20,21] and development of novel renormalization techniques [22][23][24] to deal with the strong short-range forces revived MBPT for both ground-state energy calculation [25][26][27][28][29] and for effective interaction construction [30][31][32].
In spite of extensive work, a systematic study of the order-by-order convergence of MBPT for nuclear systems was impossible until the recursive method for the RS perturbation theory was adopted and results up to 30th order were obtained [27,28].The use of chiral interactions, softened within the similarity renormalization group approach, resulted in typical divergence of the perturbation series in a harmonic-oscillator (HO) basis against convergence in a Hartree-Fock (HF) basis for 4 He and 16,24 O.Note that a meaningful result can still be extracted from a divergent series by applying resummation techniques, such as Padé approximants [27,33] and eigenvector continuation method [34,35].
Although order-by-order convergence of MBPT has been observed for RS perturbation series in the HF basis, a comprehensive understanding is still missing.What de-termines the convergence of MBPT -a many-body basis, an interaction or a specific partitioning of the Hamiltonian?Can one assess a priori the nature of a perturbative expansion?In this letter we answer those questions for nuclear ground states with Brillouin-Wigner (BW) perturbation theory.A general partitioning scheme for a Hamiltonian and a convergence criterion are proposed.A new vertex function is devised to calculate high orders of BW perturbation series iteratively.We apply these developments to 4 He and 16 O using Daejeon16 [36] and bare N 3 LO [37] potentials in both HO and HF bases.The results are benchmarked with an exact matrix inversion method and with no-core shell model (NCSM) [38].
We start with a translationally-invariant intrinsic Hamiltonian for A point-like nucleons where m is the nucleon mass (approximated here as the average of the neutron and proton mass), p i is the ith nucleon's momentum, and V ij denotes the nucleon-nucleon interaction with the addition of Coulomb interaction for protons.Introducing an auxiliary spherically-symmetric one-body potential u, we can split the Hamiltonian into a solvable part, can be ideally solved by diagonalization of the Hamiltonian matrix computed in a complete set of orthonormal A-body basis, {|Φ α }, of eigenstates of H 0 : In practice, one has to truncate the basis to be finite.If the basis is large enough, the results are almost independent of its parameters.This is NCSM -a full configuration-interaction method, feasible for light nuclei.Because of rapidly growing basis dimensions with increasing number of nucleons, for heavier nuclei the method becomes prohibitive.In such a case, one can split the full model space into a smaller model space (Èspace) and its complement (É-space) using projection operators: P = α∈È |Φ α Φ α | and Q = 1−P .Then one projects Eq. ( 2) into È-space, getting for an E k -dependent effective Hamiltonian [39] H eff (E k ) = P HP + P HQ One may proceed within RS perturbation theory, expanding H eff (E k ) around an unperturbed energy [12].In this work, we stay within E k -dependent BW formulation.
A HO potential with the oscillator quantum ω=18 MeV and a spherical HF potential are used in this work.The latter is obtained by solving the spherical HF equation inside 12 major HO shells using Hamiltonian Eq. ( 1).
For both bases, we use N max -truncation, defined in the same way as in NCSM [38].Here, for the ground states of closed-shell nuclei, we choose the È-space to be onedimensional (N P max = 0), i.e., P = |Φ 0 Φ 0 |.In general, the matrix of the H eff (E) operator can be obtained as a function of energy E by taking the inverse of the (E − QHQ) matrix as seen from Eq. ( 4).Then, the solution of Eq. ( 3) can be found numerically by the Newton-Raphson method [40,41].As an example, Fig. 1 shows the curve of y=f 0 (E) which is the eigenvalue of H eff (E) for 4 He, obtained by matrix inversion using Daejeon16 in the HO basis with N max =2.The intersections of y=E and y=f 0 (E) give J π =0 + eigenvalues E k of H eff (E k ), which are in addition characterized by zero center-of-mass (CM) excitation.The latter is due to the fact that in the HO basis, H commutes with the CM Hamiltonian H CM (actually, both H 0 and H 1 commute with H CM ), and that in the È-space, the CM motion is in its ground state.
The fact that more than one solution appear from Eq. ( 3) is due to the energy dependence of H eff (E).In general, the true wave function contains two components: Solutions, appearing in Fig. 1, are those eigenstates of Eq. ( 2) which have a non-zero |Ψ È k component (that is equivalent to J π =0 + eigenstates with zero CM excitation in this example), and |Ψ k can be re- is seen to be small, which means that the ground state is dominated by the È-space configuration (95.245%), whereas other eigenstates are seen to be dominated by the É-space component.
The curve y=f 0 (E) exhibits singularities (in blue color in Fig. 1), slightly shifted from the É-space dominated eigenenergies E k≥1 .As seen from Eq. ( 4), these are the J π =0 + eigenvalues of the QHQ operator.Again, in the FIG. 1 (color online).The exact eigenvalue f0(E) of H eff (E) for 4 He from matrix inversion using Daejeon16 in the HO basis with ω=18 MeV, Nmax=2 (dashed line).The red dots are the solutions of Eq. (3) with J π =0 + and zero CM excitation, which are exactly reproduced by NCSM.The bottom blue vertical lines mark the J π =0 + eigenvalues of QHQ with zero CM excitation, which are the singularities of f0(E).
HO basis, these eigenvalues in addition correspond to zero CM excitation.
Before treating H eff (E) perturbatively, we introduce here a partition parameter ξ, which is diagonal in the full model space, so that the exact resolvent operator can formally be expanded into a BW perturbation series, provided the expansion converges, where is the E-dependent ratio parameter defined in the Éspace.The parameter ξ can be used to tune the convergence behavior, e.g., ξ= Φ 0 |H 1 |Φ 0 is similar to the Møller-Plesset (MP) partitioning with a normal-ordered Hamiltonian [42], and [43,44], in RS perturbation theory.Below, we treat ξ as a number.From Eq. ( 5), we observe that the BW perturbation series converges to the exact value if and only if ρ(R), the spectral radius of R in É-space, satisfies the inequality ρ(R) < 1, referred to as convergence criterion here (see Theorem 5.6.12 in Ref. [45], see also Ref. [40]).Let E qhq 1 be the lowest eigenvalue of QHQ.For E<E qhq 1 , choosing ξ to be large enough so that X is negative-definite and ρ(R−1)<2, we achieve that ρ(R)<1.This can be proved using Sylvester Theorem [45] and the fact that the eigenvalues of (−X) −1/2 (E − QHQ)(−X) −1/2 and of the (R−1) operator are the same (see Supplemental Material).When Eq. ( 5) is inserted into H eff (E), E qhq 1 is in addition characterized by the same values of quantum  4 He using Daejeon16 and ξ = 0 MeV, (e) 4 He using bare N 3 LO and ξ = 0, and (f) 16 O using Daejeon16 and ξ = −700 MeV, while panel (b) shows diagonal Padé approximants to various orders of (a).Inside each panel, ρ(R) is depicted on the top sub-panel, the details of the perturbative calculation near the exact ground-state energy are given on the right sub-panel, and the exact f0(E) obtained from matrix inversion is shown with black dash-dotted line.
numbers as the eigenvalue of H eff (E) and is therefore the lowest singularity of f 0 (E).The exact ground-state energy E 0 of H in the full model space satisfies E 0 <E qhq 1 as long as the variational principle holds, i.e., the groundstate energy of H in our finite full model space is smaller than the ground-state energy of QHQ in the É-space.
In other words, with a proper choice of ξ, the expansion (5) converges in the energy interval (−∞, E qhq 1 ), in which E 0 resides.Thus we can always make the BW perturbation series for ground state converge.Moreover, for open-shell nuclei where the È-space dimension is bigger than one, we can get the converged BW series for eigenenergies smaller than E qhq 1 , which is similar to the convergence condition proposed in Refs.[46,47] for an energy-independent effective interaction.
To check these findings numerically, we introduce a new vertex function, named K-box, in the ÈÉ-space, The effective Hamiltonian can be expressed via Then one can achieve a BW expansion of H eff (E) and calculate high-order terms by K-box iterations with In the following we apply this formalism for ground states of closed-shell nuclei.Figs.2(a,c-e) show the eigenvalues f sth 0 (E) of H eff (E) for 4 He, calculated within BW perturbation theory up to sth order for different choices of ξ in the HO basis with ω=18 MeV and N max =2, using K-box iterations.The values of ρ(R), corresponding to J π = 0 + and zero CM excitation, are obtained by diagonalizing QRQ.All calculations have been performed with Daejeon16, except for Fig. 2(e), where bare N 3 LO potential is used.In Fig. 2(a), we set ξ = Φ HO 0 |H 1 |Φ HO 0 = −132.927MeV.This is a MP partitioning with normal-ordered Hamiltonian, widely used in RS perturbation theory, which may induce divergent results for ground states in the HO basis [27,28].Here, we also get the BW expansion divergent.Indeed, we notice that the spectral radius ρ(R) for the present choice of ξ is below unity in the interval (−∞, −32.072MeV), while the exact ground-state energy E 0 = −26.822MeV belongs to the divergence interval [−32.072,54.218] where ρ(R) ≥ 1.We recall that the convergent result can be retrieved, as seen from Fig. 2(b), by using diagonal Padé approximants [27,48] to the expansion of H eff (E) with the 2nd order contribution being the first term.Starting from Padé [3/3] approximant (corresponding to the 8th order of H eff (E)), the difference between the Padé resummed value of the ground-state energy and the exact result is smaller than 1 keV.Singularities of the exact f 0 (E) can also be well restored by Padé approximants.Next, we demonstrate the role of ξ on the convergence behavior.As was explained, if ξ is large enough, the convergence criterion ρ(R)<1 can be satisfied for the ground-state calculation.For example, choosing a bit larger value of ξ= −110 MeV, we are able to make ρ(R)<1 in the interval (−∞, 1.081 MeV), where E qhq 1 = 1.081MeV is the lowest J π = 0 + eigenvalue of QHQ with zero CM excitation and hence the lowest singularity of the exact f 0 (E), resulting thus in a BW series converging to the exact ground-state energy E 0 = −26.822MeV (see Fig. 2(c)).Increasing ξ further, for example, choosing ξ = 0 MeV in Fig. 2(d), we continue to get a converging perturbation expansion, although the convergence rate slows down.The above results stem from the variational principle, as remarked before, and are independent of the internucleon interaction.For example, Fig. 2  .As is seen from the figure, BW perturbation series converge towards the exact ground-state energies in both cases.This outcome is similar to the case of RS perturbation theory studied in Ref. [28].
Calculations discussed above (Figs.2-3) are all performed at N max =2 for illustrative purpose.Table I summarizes ground-state energies from perturbative calculation for 4 He and 16 O using Daejeon16 in both HO and HF bases in larger model spaces.We observe again that with proper choices of ξ, the perturbation series converges in both bases, where the exact NCSM results serve as a benchmark.In addition, an advantageous choice of ξ can make the convergence rate in the HO basis competitive to that in the HF basis with MP partitioning.More results for different ω values can be found in the Supplemental Material.
To summarize, we report novel developments within BW type MBPT, applied to ground states of closed-shell nuclei.A general partitioning of the Hamiltonian and a convergence criterion for a BW perturbative expansion are proposed.In particular, we have shown that the convergence criterion for ground-state calculations can always be satisfied in both HO and HF bases with a proper choice of the Hamiltonian partitioning, which attributes to the variational principle and does not depend on the internucleon interactions.A new vertex K-box is designed to calculate high orders of BW perturbation series and applied to the ground-state calculation of 4 He and 16 O, using Daejeon16 and a bare N 3 LO potential.Our results confirm that BW perturbation series can be redefined to be converging in both, HO and HF bases, and for both potentials without any specific renormalization procedure.Study of excited states for open-shell nuclei and generalization of these findings to RS perturbation theory are under way.We believe this work paves a way towards ab-initio description of nuclear many-body TABLE I. NCSM and BW perturbative calculations (up to sth order) for ground-state energies of 4 He and 16 O in the HO basis (with ξ = −110 MeV for 4 He and ξ = −700 MeV for 16 O in the BW calculation) and in the HF basis (with ξ = Φ HF 0 |H1|Φ HF 0 in the BW calculation) using Daejeon16 at ω = 18 MeV.All results are in MeV.The authors acknowledge the financial support from CNRS/IN2P3, France, via ENFIA and ABICIA Master projects.Large-scale computations have been performed at MCIA, University of Bordeaux.

Supplemental Material for Converging Many-Body Perturbation Theory for Ab Initio Nuclear-Structure: I. Brillouin-Wigner Perturbation Series for Closed-Shell Nuclei
We denote the eigenvalue problem of H eff (E) at an arbitrary energy E as Since the effective Hamiltonian is one-dimensional in our case, we simply have |ψ È 0 (E) = |Φ 0 and H eff (E) = f 0 (E).We observe that the derivative of the eigenvalue f 0 (E) with respect to energy E is therefore at intersections of y = f 0 (E) and y = E, i.e., f 0 (E k ) = E k , we have The eigenfunction in the full model space can be restored by a wave operator Ω [6] where the corresponding É-space component is The occupation probability ratio of É-space to È-space for the kth solution is therefore

CONVERGENCE CRITERION FOR THE BRILLOUIN-WIGNER PERTURBATION SERIES
The BW perturbation series shown in the main text is simply an operator-valued geometric series with ratio R (see Eq. ( 6) in the main text), and is convergent if and only if lim R k tends to a finite number.We define R k , then we have We have used Theorem 5.6.12 of Ref. [45] i.e., lim k→∞ R k = 0 if and only if ρ(R) < 1, in the last step of the above derivation, where ρ(R) is the spectral radius (the maximum absolute magnitude of the eigenvalues) of R. It follows that the BW perturbation series converges to the exact resolvent if and only if ρ(R) < 1, referred to as convergence criterion of the BW perturbation series here.Now let us look at the role of ξ on the convergence behavior.Note that the parameter ξ is a scalar or an operator diagonal in the É-space.The E-dependent expansion ratio can be written as which means that the convergence criterion ρ(R) < 1 can be reduced to S being negative-definite and ρ(S) < 2, i.e., the eigenvalues λ i (S) of S satisfying where d q is the dimension of É-space.Let us denote the ordered eigenvalues of QHQ in the É-space as (S12) H 0 and ξ are diagonal in the É-space, hence, they can be expressed as Since we are interested in the ground state in this work, below we consider the interval E ∈ (−∞, E qhq 1 ) in which the ground-state energy resides because of the variational principle.Note that the operator S can be written as We can adjust {ξ α } to be large enough so that Then the characteristic equation of S can be written as which means S c ≡ (−X) −1/2 (E − QHQ)(−X) −1/2 and S have the same eigenvalues.The Hermitian matrices (E − QHQ) and S c are star-congruent and therefore they have the same number of positive eigenvalues and same number of negative eigenvalues, according to Sylvester Theorem [45].Since E < E qhq 1 , (E −QHQ) is negative-definite, and S c is therefore negative-definite.We immediately have negative-definite S. The condition ρ(S) = ρ(S c ) < 2 can be satisfied by choosing large ξ α .
In the following we roughly estimate how large the values {ξ α } we should choose to make ρ(S c ) < 2. Let λ k (S c ) and λ k (E − QHQ) be the eigenvalues of S c and (E − QHQ) respectively, in ascending order, The eigenvalues of S c and E − QHQ are related by Ostrowski theorem [45], i.e., The condition ρ(S c ) < 2 can therefore be guaranteed by simply setting i.e., Now we collect all the sufficient conditions to make ρ(R) < 1 for where the second condition can be guaranteed by the other two conditions and the above conditions are reduced to The above second condition can be achieved if we simply require To summarize, we can adjust {ξ α } to let it satisfy the condition Note that this condition is just a sufficient condition of ρ(R) < 1 and does not necessarily cover all the choice of {ξ α } to make ρ(R) < 1 for E < E qhq 1 .This conclusion is universal and is true for both closed-shell nuclei and open-shell nuclei.In the above analysis we have not mentioned the symmetries preserved by the intrinsic Hamiltonian H.In practice, when the resolvent operator is inserted into H eff (E), the Éspace and the corresponding projection operator Q in the above analysis should be replaced by its subspace É s and the corresponding projection operator Q s , characterized by the same values of good quantum numbers (e.g., angular momentum J, parity, and in certain cases of other quantum numbers, such as center-of-mass quantum numbers, etc.) as the È-space eigenstates of H eff (E) because of the presence of P HQ and QHP operators in H eff (E).

ADDITIONAL NUMERICAL RESULTS
Some additional calculations can be found in Tab.I-IX. ) and Padé resummation results for the ground-state energy of 4   for the BW calculation) using Daejeon16 at ω = 18 MeV.All results are in MeV.The BW calculation with Newton-Raphson method is accurate to five decimals.

FIG. 2 (
FIG. 2 (color online).The eigenvalue f sth0 (E) of H eff (E) up to various orders s (labeled by "Os") of BW perturbation series in HO basis at ω = 18 MeV, Nmax = 2 for (a)4 He using Daejeon16 and ξ = Φ HO 0 |H1|Φ HO 0 = −132.927MeV, (c)4 He using Daejeon16 and ξ = −110 MeV, (d)4 He using Daejeon16 and ξ = 0 MeV, (e)4 He using bare N 3 LO and ξ = 0, and (f)16 O using Daejeon16 and ξ = −700 MeV, while panel (b) shows diagonal Padé approximants to various orders of (a).Inside each panel, ρ(R) is depicted on the top sub-panel, the details of the perturbative calculation near the exact ground-state energy are given on the right sub-panel, and the exact f0(E) obtained from matrix inversion is shown with black dash-dotted line.
(e) depicts the perturbation series obtained with the bare N 3 LO potential assuming ξ=0 MeV.The BW series clearly converges to the exact value E 0 = −0.6527MeV, which is smaller than the lowest singularity of f 0 (E) at 25.77 MeV.The convergence aspects of16 O in the HO basis are similar.Increasing ξ from Φ HO 0 |H 1 |Φ HO 0 = −753.239MeV to −700 MeV, the diverging series becomes converging, as seen in Fig.2(f), for energies E smaller than −89.001MeV (the lowest singularity of f 0 (E)), to which the exact ground-state energy belongs.
He in the HO basis using Daejeon16 at ω = 18 MeV.The Padé resummation f Padé[N/N] 0 (E) is performed on the BW perturbation series of the energy-dependent effective Hamiltonian H eff (E) starting from the 2nd order contribution, and the corresponding ground-state energy is obtained by finding the intersection of y = E and y = f Padé[N/N] 0 (E) + f s=1 0 (E), where f s=1 0 (E) is the eigenvalue of the first order H eff (E).All results are in MeV.

TABLE II .
NCSM and BW perturbative calculations (up to s-th order) for the ground-state energy of4He in the HO basis (with ξ = −110 MeV for the BW calculation) and in the HF basis (with ξ = Φ HF 0 |H1|Φ HF 0

TABLE III .
NCSM and BW perturbative calculations (up to s-th order) for the ground-state energy of4He in the HO basis (with ξ = −150 MeV for the BW calculation) and in the HF basis (with ξ = Φ HF 0 |H1|Φ HF 0 for the BW calculation) using Daejeon16 at ω = 25 MeV.All results are in MeV.The BW calculation with Newton-Raphson method is accurate to five decimals.

TABLE IV .
NCSM and BW perturbative calculations (up to s-th order) for the ground-state energy of4He in the HO basis (with ξ = −100 MeV for the BW calculation) and in the HF basis (with ξ = Φ HF 0 |H1|Φ HF 0 for the BW calculation) using bare N 3 LO at ω = 18 MeV.All results are in MeV.The BW calculation with Newton-Raphson method is accurate to five decimals.

TABLE V .
NCSM and BW perturbative calculations (up to s-th order) for the ground-state energy of4He in the HO basis (with ξ − 130 MeV for the BW calculation) and in the HF basis (with ξ = Φ HF 0 |H1|Φ HF 0 for the BW calculation) using bare N 3 LO at ω = 25 MeV.All results are in MeV.The BW calculation with Newton-Raphson method is accurate to five decimals.

TABLE VI .
NCSM and BW perturbative calculations (up to s-th order) for the ground-state energy of16O in the HO basis (with ξ = −700 MeV for the BW calculation) and in the HF basis (with ξ = Φ HF 0 |H1|Φ HF 0 for the BW calculation) using Daejeon16 at ω = 18 MeV.All results are in MeV.The BW calculation with Newton-Raphson method is accurate to five decimals.

TABLE VII .
NCSM and BW perturbative calculations (up to s-th order) for the ground-state energy of16O in the HO basis (with ξ = −950 MeV for the BW calculation) and in the HF basis (with ξ = Φ HF 0 |H1|Φ HF 0 for the BW calculation) using Daejeon16 at ω = 25 MeV.All results are in MeV.The BW calculation with Newton-Raphson method is accurate to five decimals.

TABLE VIII .
NCSM and BW perturbative calculations (up to s-th order) for the ground-state energy of16O in the HO basis (with ξ = −560 MeV for the BW calculation) and in the HF basis (with ξ = Φ HF 0 |H1|Φ HF 0 for the BW calculation) using bare N 3 LO at ω = 18 MeV.All results are in MeV.The BW calculation with Newton-Raphson method is accurate to five decimals.

TABLE IX .
NCSM and BW perturbative calculations (up to s-th order) for the ground-state energy of16O in the HO basis (with ξ = −750 MeV for the BW calculation) and in the HF basis (with ξ = Φ HF 0 |H1|Φ HF 0 for the BW calculation) using bare N 3 LO at ω = 25 MeV.All results are in MeV.The BW calculation with Newton-Raphson method is accurate to five decimals.