# Converging Many-Body Perturbation Theory for Ab Initio Nuclear Structure:
II. Brillouin-Wigner Perturbation Series for Open-Shell Nuclei

Zhen Li<EMAIL_ADDRESS><EMAIL_ADDRESS>LP2IB (CNRS/IN2P3 –
Université de Bordeaux), 33170 Gradignan, France  Nadezda A. Smirnova
<EMAIL_ADDRESS>LP2IB (CNRS/IN2P3 – Université de Bordeaux),
33170 Gradignan, France

###### Abstract

Brillouin-Wigner (BW) perturbation theory is developed for both ground and
excited states of open-shell nuclei. We show that with optimal partitioning of
the many-body Hamiltonian proposed earlier by the authors [Z. Li and N.
Smirnova, arXiv:2306.13629], one can redefine the BW perturbation series for a
given state of the effective Hamiltonian in a small $\mathbbm{P}$-space to be
converging under the condition that the energy of this state is below the
lowest eigenvalue of the Hamiltonian matrix block belonging to the complement
of the $\mathbbm{P}$-space, characterized by the same good quantum numbers as
the state under consideration. Specifically, the BW perturbative calculations
for the lowest $J^{\pi}$ states are always converging due to the variational
principle. This property does hold for both soft and hard internucleon
interactions in the harmonic oscillator basis. To illustrate this method and
check the convergence behavior, we present numerical studies of low-energy
spectra of 5,6,7Li using the Daejeon16 and bare N3LO potentials.

##  I Introduction

Solving the many-body Schrödinger equation for atomic nuclei is the key to the
ab initio study of nuclear structure. The last two decades have witnessed huge
progress in the developments of nuclear ab initio many-body methods [1],
accompanied by the success of chiral effective field theory (EFT) potentials
[2, 3] and renormalization techniques [4]. For few-nucleon systems and very
light nuclei, Faddeev equations [5, 6], hyperspherical harmonics method [7, 8]
and quantum Monte Carlo methods [9] are able to provide very accurate
calculations. The full configuration interaction method, no-core shell model
(NCSM) [10], gives very detailed information on nuclear states and transitions
at low energies, by diagonalizing the Hamiltonian matrix calculated in a
sufficiently large model space spanned by many-body states. However, with the
current computing power NCSM can only reach nuclei to oxygen or so at the
moment, due to the rapid growth of the model-space dimension as the nucleon
number increases. For nuclei heavier than oxygen, specific approximations need
to be applied in the calculation, which has led to a variety of many-body
methods, such as many-body perturbation theory (MBPT) [11, 12], coupled
cluster method [13], in-medium similarity renormalization group [14], self-
consistent Green’s function [15], etc. These methods have achieved a lot of
success in the ab initio study of intermediate-mass nuclei [1] and even heavy
nuclei like lead [16].

Initiated from the 1950s [17, 18, 19, 20], MBPT has been extensively applied
to the calculation of the nuclear ground-state energy [12] and to the
construction of the effective Hamiltonian used in the valence-space nuclear
shell model calculation [21, 22, 11, 23]. Furthermore, MBPT has also become a
precious tool for the solution of atomic and molecular problems [24, 25, 26,
27, 28, 29, 30]. In spite of extensive applications, MBPT is brought into
doubt by the issue of order-by-order convergence in powers of interactions,
especially for atomic nuclear systems, where (strong) internucleon
interactions are used. Previously most attention has been attracted to the
Rayleigh-Schrödinger (RS) formulation of MBPT, because of its energy-
independence feature and a convenient implementation via a diagrammatic
approach [20, 30, 31, 32]. However, a huge number of diagrams required at
higher orders limited the application of this diagrammatic approach up to
usually 3rd order in nuclear structure calculations [33, 34, 35, 36, 37, 12,
22, 11], which makes it impossible for the study of order-by-order convergence
within this diagrammatic framework. The study of the convergence properties
became possible with the implementation of the algebraic recursive calculation
of the RS perturbation series terms [38, 39]. Extensive calculations for
molecular systems with various partitioning schemes reported cases of
converging or diverging series [40]. Similar calculations have also been
performed for atomic nuclei [41, 42, 43], to which we focus our attention in
the present paper. Although the recursive method was still restricted to
relatively light nuclei, some convergence properties were addressed in those
studies after numerical calculations. With the use of softened chiral EFT
potentials, the authors [41, 42, 43] concluded that the perturbative
calculation of nuclear ground and excited states typically diverges in the HO
basis, and only converges in the Hartree-Fock (HF) basis. However, harder
potentials can still spoil the convergence in the HF basis [43].

In our recent work [44], we chose to explore the energy-dependent BW
perturbation series for ground states of closed-shell nuclei. We showed that
the with an optimal Hamiltonian partitioning, the convergence criterion of a
perturbative expansion for the ground-state energy can always be satisfied.
This property holds due to the variational principle and does not depend on
the choice of the (HO or HF) basis or the choice of the (hard or soft)
internucleon interaction. In the present work, we generalize these ideas to
ground and excited states of open-shell nuclei and investigate the
corresponding convergence behavior. This paper is organized as follows. In
Sec. II we present the formalism of BW perturbation theory, the general
partitioning of Hamiltonian, the convergence criterion, and the $\hat{K}$-box
vertex to reach high order terms of the perturbation series. In Sec. III we
show the results of calculations for 5,6,7Li. The last section summarizes the
results and formulates perspectives.

##  II Formalism: Hamiltonian partitioning and convergence criterion

In this section, we generalize the BW formulation of MBPT given in Ref. [44]
to open-shell nuclei.

###  II.1 $\mathbbm{P}$-Space Eigenvalue Problem

Ideally, one searches to solve the eigenvalue problem for an intrinsic
Hamiltonian $H$ in a large, but finite-dimensional (denoted as
$d$-dimensional) model space

$\displaystyle H|\Psi_{k}\rangle=E_{k}|\Psi_{k}\rangle,\ k=0,1,2,\ldots,d-1.$
(1)

When the dimension $d$ becomes prohibitive for the exact diagonalization
method, we can resort to the projection technique, i.e., to project the full-
model-space eigenvalue problem into a much smaller model space (called
$\mathbbm{P}$-space) eigenvalue problem of an effective Hamiltonian. In this
context, given a complete set of basis states, $\\{|\Phi_{k}\rangle\\}$
($k=0,1,2,\ldots,d-1$), we introduce two projection operators $P$ and $Q$. $P$
projects the full $d$-dimensional model space on a smaller $d_{p}$-dimensional
$\mathbbm{P}$-space,

$\displaystyle P\equiv\sum_{k\in\mathbbm{P}}|\Phi_{k}\rangle\langle\Phi_{k}|$
(2)

while $Q$ projects the full model space on the $d_{q}$-dimensional
complementary space (called $\mathbbm{Q}$-space),

$\displaystyle Q=\sum_{k\in\mathbbm{Q}}|\Phi_{k}\rangle\langle\Phi_{k}|,$ (3)

where $d_{p}+d_{q}=d$. These projection operators satisfy

$\displaystyle P+Q=1,\,P^{2}=P,\,Q^{2}=Q,\,PQ=QP=0.$ (4)

To reduce the dimensionality of the eigenvalue problem, we project the
eigenvalue equation (1) for $H$ into the small $\mathbbm{P}$-space,

$\displaystyle H_{\rm
eff}(E_{k})|\Psi_{k}^{\mathbbm{P}}\rangle=E_{k}|\Psi_{k}^{\mathbbm{P}}\rangle,\
k=0,1,2,\ldots,d_{p}-1,$ (5)

for an energy-dependent effective Hamiltonian [45]

$\displaystyle H_{\rm eff}(E_{k})\equiv PHP+PHQ\frac{1}{E_{k}-QHQ}QHP,$ (6)

where $|\Psi_{k}^{\mathbbm{P}}\rangle\equiv P|\Psi_{k}\rangle$ is the
$\mathbbm{P}$-space component of the $k$th full model space eigenstate
$|\Psi_{k}\rangle$. The full model space eigenstates of Eq. (1) can be
restored via

$\displaystyle|\Psi_{k}\rangle=|\Psi_{k}^{\mathbbm{P}}\rangle+\frac{1}{E_{k}-QHQ}QHP|\Psi_{k}^{\mathbbm{P}}\rangle,$
(7)

and the corresponding $\mathbbm{Q}$-space component is thus

$\displaystyle|\Psi_{k}^{\mathbbm{Q}}\rangle\equiv
Q|\Psi_{k}\rangle=\frac{1}{E_{k}-QHQ}QHP|\Psi_{k}^{\mathbbm{P}}\rangle.$ (8)

We note that the effective Hamiltonian depends on the exact eigenenergies
$E_{k}$, and hence the $\mathbbm{P}$-space eigenvalue problem Eq. (5) needs to
be solved self-consistently.

Let us define the eigenvalue problem for $H_{\rm eff}(E)$ at an arbitrary
energy $E$ as

$\displaystyle H_{\rm eff}(E)|\psi_{n}(E)\rangle=f_{n}(E)|\psi_{n}(E)\rangle,\
n=0,1,\ldots,d_{p}{-}1.$ (9)

Here energy $E$ does not necessarily coincide with the eigenenergy and
therefore the eigenvalues $f_{n}(E)$ and eigenvectors $|\psi_{n}(E)\rangle$
are functions of energy $E$. Remark that each $f_{n}(E)$ exhibits
singularities because of the presence of the resolvent operator $(E-QHQ)^{-1}$
in $H_{\rm eff}(E)$, as shown in Eq. (6). Obviously, the singularities are the
eigenvalues of $QHQ$ characterized by the same good quantum numbers (angular
momentum, parity, etc) as carried by $f_{n}(E)$, due to the presence of $PHQ$
and $QHP$ operators in Eq. (6). The eigenvalues $E_{k}$ of Eq. (5) are the
energies satisfying $f_{n}(E)=E$. In practice these energies can be found by
finding the intersections of $y=E$ and $y=f_{n}(E)$, which can be quickly
located with the Newton-Raphson method. At each intersection $E=E_{k}$, only
one of the $d_{p}$ eigenpairs of Eq. (9) satisfies $f_{n}(E_{k})=E_{k}$ and
$|\psi_{n}(E_{k})\rangle=|\Psi_{k}^{\mathbbm{P}}\rangle$, assuming that the
eigenenergies of Eq. (5) are nondegenerate. Examples with one-dimensional
$\mathbbm{P}$-space can be found in Ref. [44] and are also given in Fig. 1 of
this paper. The cases with multi-dimensional $\mathbbm{P}$-space are given in
Figs. 4, 9, and 12.

The first derivative of $f_{n}(E)$ can be easily obtained

$\displaystyle
f_{n}^{\prime}(E)=-\frac{\langle\psi_{n}(E)|PHQ\frac{1}{(E-QHQ)^{2}}QHP|\psi_{n}(E)\rangle}{\langle\psi_{n}(E)|\psi_{n}(E)\rangle}.$
(10)

We note that $f_{n}^{\prime}(E)\leq 0$ for all $E$. At the intersections,
i.e., $f_{n}(E_{k})=E_{k}$,
$|\psi_{n}(E_{k})\rangle=|\Psi_{k}^{\mathbbm{P}}\rangle$, we have

$\displaystyle
f_{n}^{\prime}(E_{k})=-\frac{\langle\Psi_{k}^{\mathbbm{P}}|PHQ\frac{1}{(E_{k}-QHQ)^{2}}QHP|\Psi_{k}^{\mathbbm{P}}\rangle}{\langle\Psi_{k}^{\mathbbm{P}}|\Psi_{k}^{\mathbbm{P}}\rangle}\leq
0,$ (11)

which is related to the occupation probability ratio of $\mathbbm{Q}$-space to
$\mathbbm{P}$-space components of the $k$-th eigenstate of Eq. (5) via

$\displaystyle\frac{\langle\Psi_{k}|Q|\Psi_{k}\rangle}{\langle\Psi_{k}|P|\Psi_{k}\rangle}$
$\displaystyle=$
$\displaystyle\frac{\langle\Psi_{k}^{\mathbbm{P}}|PHQ\frac{1}{(E_{k}-QHQ)^{2}}QHP|\Psi_{k}^{\mathbbm{P}}\rangle}{\langle\Psi_{k}^{\mathbbm{P}}|\Psi_{k}^{\mathbbm{P}}\rangle}$
(12) $\displaystyle=$ $\displaystyle-f_{n}^{\prime}(E_{k})$
$\displaystyle\geq$ $\displaystyle 0.$

It follows that for $\mathbbm{P}$-space dominated eigenstates, we have
$|f_{n}^{\prime}(E_{k})|<1$, while for $\mathbbm{Q}$-space dominated
eigenstates we have $|f_{n}^{\prime}(E_{k})|>1$.

Because of the energy-dependence of the effective Hamiltonian, all the
eigenenergies $E_{k}$ (carried by the same good quantum numbers as carried by
the $\mathbbm{P}$-space) of the full-model-space eigenvalue problem Eq. (1)
can be found with the method described above. The corresponding
$\mathbbm{P}$-space eigenvectors $|\Psi_{k}^{\mathbbm{P}}\rangle$ for these
eigenenergies are non-orthogonal, since they are the eigenstates of the
Hermitian effective Hamiltonian $H_{\rm eff}(E)$ at different energies.

Until now, we have not mentioned the way to construct the effective
Hamiltonian $H_{\rm eff}(E)$. The straightforward and exact way to construct
$H_{\rm eff}(E)$ is by taking the inverse of the $(E-QHQ)$ matrix in the
$\mathbbm{Q}$-space directly. This matrix inversion method, however, cannot be
applied to large model spaces, and hence it serves mainly as a benchmark in
this work. In the following subsections, we show how to compute $H_{\rm
eff}(E)$ perturbatively within the BW framework.

###  II.2 Convergence Criterion of the Brillouin-Wigner Perturbation Series
for the Effective Hamiltonian

We use in Eq. (1) the $A$-body intrinsic Hamiltonian

$\displaystyle H$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{A}\left(\frac{\bm{p}_{i}^{2}}{2m}+u_{i}\right)+\sum_{\begin{subarray}{c}i<j\end{subarray}}^{A}\left(V_{ij}-\frac{\bm{p}_{i}\cdot\bm{p}_{j}}{mA}\right)$
(13)
$\displaystyle-\sum_{i=1}^{A}\left(u_{i}+\frac{\bm{p}_{i}^{2}}{2mA}\right),$

where $m$ is the nucleon mass (approximated here as the average of the neutron
and proton mass), $\bm{p}_{i}$ is the $i$th nucleon’s momentum, and $V_{ij}$
is the nucleon-nucleon interaction with additional Coulomb interaction for
protons. We have introduced an auxiliary potential $u$ so that the eigenvalue
problem for the first term
$H_{0}\equiv\sum_{i=1}^{A}\left(\bm{p}_{i}^{2}/2m+u_{i}\right)$ (i.e., the
unperturbed Hamiltonian) can be easily solved. The rest part of the intrinsic
Hamiltonian $H_{1}\equiv H-H_{0}$ is the residual interaction. In this work,
we use the spherically-symmetric HO potential as the auxiliary potential,
i.e., $u=\frac{1}{2}m\omega^{2}\bm{r}^{2}$, where $\omega$ is the HO angular
frequency. We denote the eigenvalue equation of $H_{0}$ as

$\displaystyle H_{0}|\Phi_{k}\rangle=\mathcal{E}_{k}|\Phi_{k}\rangle,$ (14)

where $\mathcal{E}_{k}$ are the sum of $A$ single-particle HO energies, and
$|\Phi_{k}\rangle$ are the $A$-body HO Slater determinants. In this work, the
full model space is spanned by the $M$-scheme HO Slater determinants
(characterized by good quantum numbers: $z$-component of total angular
momentum $M$, parity $\pi$ and $z$-component of total isospin $M_{T}$) and
truncated by $N_{\rm max}$ (the total HO excitation quantum above the lowest
configuration), as commonly used in NCSM calculations [10]. We choose the
$M$-scheme basis state(s) with the lowest energy (the so-called $0\hbar\omega$
model space) as the $\mathbbm{P}$-space. Therefore the $\mathbbm{P}$-space is
one-dimensional for closed-shell nuclei and multi-dimensional for open-shell
nuclei.

Now let us consider the exact resolvent operator $(E-QHQ)^{-1}$ in $H_{\rm
eff}(E)$ and introduce an additional partitioning parameter $\xi$, which is an
operator, diagonal in the $\mathbbm{Q}$-space. The resolvent operator can be
expanded to a BW perturbation series

$\displaystyle\frac{1}{E-QHQ}=\frac{1}{E-QH_{0}Q-QH_{1}Q}$ (15)
$\displaystyle=$ $\displaystyle\frac{1}{\underbrace{(E-QH_{0}Q-Q\xi
Q)}_{X}-\underbrace{(QH_{1}Q-Q\xi Q)}_{Y}}$ $\displaystyle=$
$\displaystyle\frac{1}{X}+\frac{1}{X}Y\frac{1}{X}+\frac{1}{X}Y\frac{1}{X}Y\frac{1}{X}+\frac{1}{X}Y\frac{1}{X}Y\frac{1}{X}Y\frac{1}{X}+\cdots$
$\displaystyle=$
$\displaystyle\lim_{n\rightarrow\infty}\sum_{k=0}^{n}R^{k}\frac{1}{X},$

where

$\displaystyle R\equiv\frac{1}{X}Y=\frac{1}{E-QH_{0}Q-Q\xi Q}(QH_{1}Q-Q\xi Q)$
(16)

is the expansion ratio depending on the energy $E$, and $X$ is diagonal in the
$\mathbbm{Q}$-space so that its inverse can be easily calculated. In the
derivation of Eq. (15), we repeatedly used the following identity:

$\displaystyle\frac{1}{X-Y}=\frac{1}{X}+\frac{1}{X-Y}Y\frac{1}{X}.$ (17)

Inserting Eq. (15) back to the effective Hamiltonian $H_{\rm eff}(E)$, we are
able to construct $H_{\rm eff}(E)$ perturbatively, providing the BW
perturbation series is converging. The BW series in Eq. (15) is nothing but an
operator-valued geometric series, which converges to the exact resolvent
operator if and only if the spectral radius (the maximum absolute eigenvalues)
of $R$ is smaller than unity, i.e., $\rho(R)<1$ [46, 47]. We refer to this
condition as the convergence criterion of BW perturbation series [44].
Noticeably, the smaller the value of $\rho(R)$ is, the faster the speed of
convergence will be. The partitioning parameter $\xi$ we introduced here can
be used to tune the convergence behavior of the BW series.

Note that the expansion ratio $R$ depends on the energy $E$. Let $E_{\rm
min}^{qhq}$ and $E_{\rm max}^{qhq}$ be the lowest and the highest eigenvalues
of the operator $QHQ$ in the $\mathbbm{Q}$-space. As we have shown in Ref.
[44], for energy $E<E_{\rm min}^{qhq}$ the convergence criterion $\rho(R)<1$
for $f_{n}(E)$ can always be satisfied with a specific choice of the
partitioning parameter $\xi$. We refer to this energy interval
$-\infty<E<E_{\rm min}^{qhq}$ as the convergence interval. This convergence
interval was also discussed in Refs. [48, 49] in the context of the
convergence of a Rayleigh-Schrödinger perturbative calculation. Recall that
$H_{0}$ and $\xi$ are diagonal operators in the $\mathbbm{Q}$-space. We denote
their diagonal entries as
$\\{\mathcal{E}_{1}^{q}\leq\mathcal{E}_{2}^{q}\leq\mathcal{E}_{3}^{q}\leq\cdots\leq\mathcal{E}_{d_{q}}^{q}\\}$
and
$\\{\xi_{1}^{q}\leq\xi_{2}^{q}\leq\xi_{3}^{q}\leq\cdots\leq\xi_{d_{q}}^{q}\\}$,
respectively. We have shown in the supplemental material of Ref. [44] that the
convergence criterion $\rho(R)<1$ for $E<E_{\rm min}^{qhq}$ can be always
satisfied by adjusting $\xi_{k}^{q}$ ($k=1,2,\ldots,d_{q}$) such that

$\mathcal{E}_{k}^{q}+\xi_{k}^{q}>\dfrac{1}{2}\left(E_{\rm min}^{qhq}+E_{\rm
max}^{qhq}\right),\ k=1,2,\ldots,d_{q}.$ (18)

This means that the diagonal entries $\xi_{k}^{q}$ should be large enough to
satisfy the above condition. Note that the above inequality (18) is just a
sufficient condition and does not necessarily cover all values of $\xi$ which
make the perturbation series converging.

Inserting the BW perturbation series Eq. (15) back into the effective
Hamiltonian $H_{\rm eff}(E)$, we have to consider the symmetries preserved by
the intrinsic Hamiltonian $H$ because of the presence of the $PHQ$ and $QHP$
operators in $H_{\rm eff}(E)$, as shown in Eq. (6). Therefore, the
$\mathbbm{Q}$-space and its projection operator $Q$ used in the above analysis
should be replaced with its subspace $\mathbbm{Q}_{s}$ and the corresponding
projection operator $Q_{s}$, characterized by the same good quantum numbers
(e.g. angular momentum $J$, parity $\pi$, center-of-mass quantum numbers
$N_{\rm cm}$, $L_{\rm cm}$, $M_{\rm cm}$, etc) as the eigenvalues $f_{n}(E)$.
Therefore $E_{\rm min}^{qhq}$ becomes the lowest eigenvalue of the
$Q_{s}HQ_{s}$ operator, which is exactly the lowest singularity of $f_{n}(E)$.
Note that we choose the $0\hbar\omega$ model space as our $\mathbbm{P}$-space,
whose center-of-mass motion is in the ground state, i.e., $N_{\rm cm}=L_{\rm
cm}=M_{\rm cm}=0$. Therefore, the subspace $\mathbbm{Q}_{s}$ does not contain
spurious center-of-mass motions.

In Ref. [44] we concluded that for ground states of closed-shell nuclei, the
convergence criterion $\rho(R)<1$ is always satisfied by choosing large enough
$\xi$ due to the variational principle. For open-shell nuclei, the
$0\hbar\omega$ model space is not one-dimensional anymore, and hence excited
states can also be calculated. The conclusion for the ground states of closed-
shell nuclei can now be generalized to the lowest $J^{\pi}$ states for open-
shell nuclei, i.e., the convergence criterion can always be satisfied for each
lowest $J^{\pi}$ state by choosing large enough $\xi$ (namely large enough
diagonal entries of $\xi$) due to the variational principle and the conserved
quantum numbers $J$ and $\pi$. For excited states of each $J^{\pi}$, only
those states with eigenenergies in the convergence interval $E_{k}<E_{\rm
min}^{qhq}$ can get converged. The presence of intruder
states111$\mathbbm{Q}$-space dominated states which occur in the energy region
of $\mathbbm{P}$-space dominated states. may reduce the size of the
convergence interval and spoil the convergence of higher excited states of
each $J^{\pi}$. We will discuss this with numerical results in Sec. III.

###  II.3 $\hat{K}$-Box Iterative Calculation of High-Order Terms of
Brillouin-Wigner Perturbation Series

Before going into the numerical calculations, we here outline the
$\hat{K}$-box iterative method introduced in Ref. [44] to efficiently
calculate the high order terms of the BW perturbation series for the effective
Hamiltonian $H_{\rm eff}(E)$. Applying Eq. (17) into the exact resolvent
operator in $H_{\rm eff}(E)$ we obtain

$\displaystyle H_{\rm eff}(E)=PHP+P\hat{K}(E)Q\frac{1}{E-QH_{0}Q-Q\xi Q}QHP,$

where we have defined a special vertex function $\hat{K}$-box in
$\mathbbm{P}\mathbbm{Q}$-space, introduced first in Ref. [44] :

$\displaystyle\hat{K}(E)\equiv PHQ+PHQ\frac{1}{E-QHQ}(QH_{1}Q-Q\xi Q).$ (20)

Applying again Eq. (17) to the resolvent operator in the above equation, we
obtain the following recursive equation for $\hat{K}(E)$:

$\displaystyle\hat{K}(E)=PHQ+P\hat{K}(E)Q\frac{1}{E-(H_{0}+\xi)}Q(H_{1}-\xi)Q.$

The value of $\hat{K}(E)$ can be calculated iteratively via

$\displaystyle\hat{K}^{(s)}(E)$ $\displaystyle=$ $\displaystyle PHQ$
$\displaystyle+P\hat{K}^{(s-1)}(E)Q\frac{1}{E-(H_{0}+\xi)}Q(H_{1}-\xi)Q,$

starting with $\hat{K}^{(0)}(E)=0$, where $s=1,2,3,\cdots$. The $(s-1)$th
iterative result $\hat{K}^{(s-1)}(E)$ corresponds to the accumulated BW
perturbation series of $H_{\rm eff}(E)$ up to the $s$th order, $H_{\rm
eff}^{s{\rm th}}(E)$. We denote the corresponding eigenvalues as $f_{n}^{s{\rm
th}}(E)$. With the above iterative calculations of $\hat{K}$-box, we can
easily reach high orders of BW perturbation series of $H_{\rm eff}(E)$. The
advantage of using $\hat{K}$-box is that we only need to store a $(d_{p}\times
d_{q})$ matrix $\hat{K}$ in the memory during iterations.

The full model space eigenvectors can be calculated from $\hat{K}$-box via

$\displaystyle|\Psi_{k}\rangle=|\Psi_{k}^{\mathbbm{P}}\rangle+\frac{1}{E_{k}-QH_{0}Q-Q\xi
Q}\hat{K}^{\dagger}(E_{k})|\Psi_{k}^{\mathbbm{P}}\rangle.$ (23)

where we made use of another equality, namely,

$\displaystyle\frac{1}{X-Y}=\frac{1}{X}+\frac{1}{X}Y\frac{1}{X-Y}.$ (24)

Therefore other physical observables apart from energies can also be
calculated with the full-model-space eigenvectors $|\Psi_{k}\rangle$.

(a)  (b)

Figure 1: The eigenvalue functions $f_{0}(E)$ ($J^{\pi}=1/2^{-},3/2^{-}$) of
$H_{\rm eff}(E)$ for 5Li calculated from the matrix inversion method using (a)
the Daejeon16 potential and (b) the bare N3LO potential with $\hbar\omega=18$
MeV at $N_{\rm max}=2$. The blue dots and red stars are the solutions of Eq.
(5), which can be exactly reproduced by NCSM. The bottom red vertical lines
mark the positions of the lowest eigenvalues of $QHQ$ with
$J^{\pi}=1/2^{-},3/2^{-}$ (zero CM excitation), which are the lowest
singularities of $f_{0}(E)$ at $J^{\pi}=1/2^{-},3/2^{-}$, respectively.

Figure 2: The eigenvalue functions $f_{0}^{s{\rm th}}(E)$
($J^{\pi}=1/2^{-},3/2^{-}$) of $H_{\rm eff}^{s{\rm th}}(E)$ from BW
perturbative calculations up to various orders $s$ for 5Li using the Daejeon16
potential with $\xi=\langle{\rm C}|H_{1}|{\rm C}\rangle=-130.227$ MeV at
$\hbar\omega=18$ MeV, $N_{\rm max}=2$. Inside each panel, the value of
$\rho(R)$ with the same $J^{\pi}$ is depicted on the top, and the exact
$f_{0}(E)$ obtained from matrix inversion is plotted with black dash-dotted
lines. The blue dots mark the exact result from NCSM calculations below the
lowest singularities.

(a)

(b)

Figure 3: The two lowest energies ($J^{\pi}=1/2^{-},3/2^{-}$) of 5Li from the
BW perturbative calculation up to various orders $s$ using (a) the Daejeon16
potential with $\xi=\langle{\rm C}|H_{1}|{\rm C}\rangle=-130.227$ MeV and (b)
the bare N3LO potential with $\xi=\langle{\rm C}|H_{1}|{\rm
C}\rangle=-103.136$ MeV, at $\hbar\omega=18$ MeV, $N_{\rm max}=2,4,6$. Inside
each panel, the energy difference in absolute value between the BW
perturbative calculation and the NCSM calculation is also plotted.

Figure 4: The eigenvalue functions $f_{n}(E)$
($J^{\pi}=0^{+},1^{+},2^{+},3^{+}$) of $H_{\rm eff}(E)$ for 6Li calculated
from the exact matrix inversion method using the Daejeon16 potential at
$\hbar\omega=18$ MeV, $N_{\rm max}=2$. The blue dots ($\mathbbm{P}$-space
dominated) and red stars ($\mathbbm{Q}$-space dominated) are the results from
NCSM calculations, which are also the intersections of $y=f_{n}(E)$ and $y=E$
(i.e., solutions of the $\mathbbm{P}$-space eigenvalue problem Eq. (5)). The
bottom vertical lines mark the position of the lowest eigenvalues of $QHQ$
with $J^{\pi}=0^{+},1^{+},2^{+},3^{+}$ and zero CM excitation, which are the
lowest singularities of $f_{n}(E)$ at $J^{\pi}=0^{+},1^{+},2^{+},3^{+}$
respectively. The details of the eigenvalue “crossing” marked by the gray
rectangle shadow in the panel of $J^{\pi}=1^{+}$ is plotted in Fig. 5. Figure
5: The details of the eigenvalue “crossing” marked by gray rectangle shadow in
the $J^{\pi}=1^{+}$ panel of Fig. 4. The eigenvalue functions $f_{n}(E)$ avoid
crossing.

Figure 6: The eigenvalue functions $f_{n}^{s{\rm th}}(E)$ ($J^{\pi}=0^{+}$) of
$H_{\rm eff}^{s{\rm th}}(E)$ from the BW perturbative calculation up to
various orders $s$ for 6Li using the Daejeon16 potential with $\xi=\langle{\rm
C}|H_{1}|{\rm C}\rangle=-128.427$ MeV, at $\hbar\omega=18$ MeV, $N_{\rm
max}=2$. For clarity, different states of perturbative calculations are
depicted separately. Inside each panel, the value of $\rho(R)$ with the same
$J^{\pi}$ is depicted on the top, and the exact $f_{n}(E)$ obtained from the
matrix inversion method is shown with black dash-dotted lines. The blue dots
mark the exact result from NCSM calculations below the lowest singularities.

Figure 7: The energies of low-lying states of 6Li
($J^{\pi}=0^{+},1^{+},2^{+},3^{+}$) from the BW perturbative calculation up to
various orders $s$ using the Daejeon16 potential with $\xi=\langle{\rm
C}|H_{1}|{\rm C}\rangle=-128.427$ MeV, at $\hbar\omega=18$ MeV, $N_{\rm
max}=2,4,6$. Inside each panel, the energy difference in absolute value
between the BW perturbative calculation and the NCSM calculation is also
plotted.

Figure 8: Same as Fig. 7 but with the bare N3LO potential and $\xi=\langle{\rm
C}|H_{1}|{\rm C}\rangle=-101.336$ MeV.

Figure 9: The eigenvalue functions $f_{n}(E)$
($J^{\pi}=1/2^{-},3/2^{-},5/2^{-},7/2^{-}$) of $H_{\rm eff}(E)$ for 7Li
calculated from the exact matrix inversion method using the Daejeon16
potential at $\hbar\omega=18$ MeV, $N_{\rm max}=2$. The blue dots
($\mathbbm{P}$-space dominated) and red stars ($\mathbbm{Q}$-space dominated)
are the eigenenergies of the $\mathbbm{P}$-space eigenvalue problem Eq. (5).
The red bottom vertical lines mark the position of the lowest eigenvalues of
$QHQ$ with corresponding $J^{\pi}$ and zero CM excitation, which are the
lowest singularities of $f_{n}(E)$.

Figure 10: The energies of low-lying states of 7Li ($J^{\pi}=1/2^{-}$) from
the BW perturbative calculation up to various orders using the Daejeon16
potential with $\xi=\langle{\rm C}|H_{1}|{\rm C}\rangle=-127.141$ MeV, at
$\hbar\omega=18$ MeV, $N_{\rm max}=2,4$. Inside each panel, the energy
difference in absolute value between the BW perturbative calculation and the
NCSM calculation is also plotted.

Figure 11: Same as Fig. 10 but with the bare N3LO potential with
$\xi=\langle{\rm C}|H_{1}|{\rm C}\rangle=-100.050$ MeV Figure 12: Same as Fig.
9 for $J^{\pi}=1/2^{-}$ of 7Li but with the bare N3LO potential.

Figure 13: NCSM calculations for seven lowest $J^{\pi}=1/2^{-}$ states of 7Li
using the Daejeon16 potential and the bare N3LO potential, at $N_{\rm max}=2$.
The lowest $\mathbbm{Q}$-space dominated state is marked in red.

##  III Results

In this section we apply the above BW formulation of MBPT to open-shell nuclei
5,6,7Li to illustrate the convergence criterion at work. The high order terms
of the BW perturbation series of the effective Hamiltonian are calculated from
the $\hat{K}$-box iterative method. We use the exact matrix inversion method
and the NCSM to benchmark the results from perturbative calculations. Two
nucleon-nucleon potentials, Daejeon16 [50] and bare N3LO [51], with additional
Coulomb potential for protons, are used in the numerical calculation. The
Daejeon16 potential is based on a similarity renormalization group (SRG) [4]
evolved N3LO potential [52, 51] and hence is a soft potential. The $M$-scheme
HO many-body basis and $N_{\rm max}$-truncation are used, the same as the
basis used in NCSM. We take the $0\hbar\omega$ model space as the
$\mathbbm{P}$-space. In this way, the center-of-mass motion is in the ground-
state for $\mathbbm{P}$-space basis, and therefore the eigenvalues of $H_{\rm
eff}(E)$ are also characterized by zero CM excitation quanta. This is due to
the fact that the intrinsic Hamiltonian $H$ commutes with the center-of-mass
Hamiltonian, and the presence of $PHQ$ and $QHP$ operators in Eq. (5) makes
sure that nonzero contributions from the $\mathbbm{Q}$-space have the same CM
motion quantum numbers as the $\mathbbm{P}$-space.

###  III.1 5Li

Figure 1(a) depicts the eigenvalue functions $f_{0}(E)$
($J^{\pi}=1/2^{-},3/2^{-}$) of $H_{\rm eff}(E)$ for 5Li calculated from the
exact matrix inversion method using the Daejeon16 potential at
$\hbar\omega=18$ MeV, $N_{\rm max}=2$. For 5Li, the $\mathbbm{P}$-space is
two-dimensional characterized with $J^{\pi}=1/2^{-},3/2^{-}$. For each
$J^{\pi}$, the $\mathbbm{P}$-space is one-dimensional. We see that the first
derivative of the eigenvalue $f_{0}(E)$ is negative, as discussed in Sec.
II.1. The intersections between $y=f_{n}(E)$ and $y=E$, marked with blue dots
and red stars in the figure, are the solutions of the $\mathbbm{P}$-space
eigenvalue problem Eq. (5), which should reproduce the results from NCSM
calculations. Indeed, we find that these intersections can be exactly
reproduced with NCSM calculations. The NCSM calculation shows that the states
of blue dots are $\mathbbm{P}$-space dominated and red stars
$\mathbbm{Q}$-space dominated, which coincides with the property of the
eigenvalues $f_{n}(E)$ addressed in Sec. II.1, i.e.,
$|f_{n}^{\prime}(E_{k})|<1$ and $|f_{n}^{\prime}(E_{k})|>1$ correspond to
$\mathbbm{P}$\- and $\mathbbm{Q}$-space dominated states respectively. Indeed,
we observe at blue dots the slopes of $f_{0}(E)$ are small and at red stars
the slopes are large. Another feather of $f_{0}(E)$ is that it exhibits
singularities, which are nothing but the eigenvalues of the $QHQ$ operator in
the $\mathbbm{Q}$-space. The lowest eigenvalue $E_{\rm min}^{qhq}$ of $QHQ$
for each $J^{\pi}$ is marked with the bottom red vertical line. Recall that
for energy $E<E_{\rm min}^{qhq}$, we can always make the BW perturbative
calculation of $H_{\rm eff}(E)$ converging by choosing large enough diagonal
entries of $\xi$. As shown in the figure, the two lowest $J^{\pi}$ states are
both blow the corresponding first singularities, which is guaranteed by the
variational principle. This means that the BW perturbative calculation of
these two states can get converged. We will see this later. A similar
calculation with the bare N3LO potential is plotted in Fig. 1(b). The
properties of $f_{0}(E)$ are similar to the case of Daejeon16. Since the bare
N3LO potential is not softened, the eigenenergies of 5Li calculated with the
N3LO potential are higher that those calculated with the Daejeon16 potential.
Nevertheless, for each $J^{\pi}$ the lowest state is still lower than the
first singularity. As discussed before, this comes from the variational
principle, and is independent of the softness of the potential. Therefore, we
expect converging BW perturbative calculations for these two lowest $J^{\pi}$
states.

Now let us go to the BW perturbative calculation. As we expected, the BW
perturbation series for the effective Hamiltonian $H_{\rm eff}(E)$ is
converging in the energy interval $E<E_{\rm min}^{qhq}$ with large enough
diagonal entries of $\xi$. Equation (18) suggests that all the diagonal
entries of $\xi$ (i.e., $\xi_{k}^{q}$) in the $\mathbbm{Q}$-space should be
larger than $-144.782$ MeV and $-149.936$ MeV for $J^{\pi}=1/2^{-}$ and
$J^{\pi}=3/2^{-}$ respectively. We here simply choose $\xi_{k}=\langle{\rm
C}|H_{1}|{\rm C}\rangle=-130.227$ MeV, i.e., $\xi$ is a scalar, where $|{\rm
C}\rangle$ is the Slater determinant with the lowest energy of the core
nucleus 4He. This choice of $\xi$ is similar to the Møller-Plesset
partitioning [53] with a normal-ordered Hamiltonian, commonly used in the RS
formalism of MBPT [12]. In Fig. 2, we show the eigenvalues $f_{0}^{s{\rm
th}}(E)$ ($J^{\pi}=1/2^{-},3/2^{-}$) of $H_{\rm eff}^{s{\rm th}}(E)$
constructed from the BW perturbation series up to various orders $s$ for 5Li
using the Daejeon16 potential at $\hbar\omega=18$ MeV, $N_{\rm max}=2$. The
spectral radius $\rho(R)$ is also plotted for each $J^{\pi}$. We observe that
in the energy intervals $E<10.152$ MeV for $J^{\pi}=1/2^{-}$ and $E<4.299$ MeV
for $J^{\pi}=3/2^{-}$, the convergence criterion $\rho(R)<1$ are satisfied.
The energies $10.152$ MeV and $4.299$ MeV are nothing but the lowest
eigenvalues of $QHQ$ operator with $J^{\pi}=1/2^{-}$ and $J^{\pi}=3/2^{-}$
respectively, i.e., the lowest singularities of $f_{0}(E)$. Therefore the
perturbative calculations inside the convergence intervals should get
converged. Indeed, we observe for each $J^{\pi}$ the eigenvalue $f_{0}^{s{\rm
th}}(E)$ is converging to the exact $f_{0}(E)$ calculated from the matrix
inversion method in the interval $E<E_{\rm min}^{qhq}$. The lowest states
(intersections) of $J^{\pi}=1/2^{-}$ and $J^{\pi}=3/2^{-}$ are both in the
convergence interval, and thus can be obtained via converging BW perturbative
calculations. The intersections residing outside the convergence interval
cannot be obtained via BW perturbative calculations. The calculation with the
bare N3LO potential is similar to the case of Daejeon16.

The above calculations are done at $N_{\rm max}=2$ for illustrative purpose.
In practice, we use the Newton-Raphson method to quickly locate the position
of the intersections between $y=f_{n}(E)$ and $y=E$. Figure 3(a) shows the BW
perturbative calculation of the two lowest states ($J^{\pi}=1/2^{-},3/2^{-}$)
of 5Li using the Daejeon16 potential with $\xi=\langle{\rm C}|H_{1}|{\rm
C}\rangle=-130.227$ MeV, at $\hbar\omega=18$ MeV, $N_{\rm max}=2,4,6$. The
NCSM calculation with the same potential and model space is also performed to
provide benchmarks. As discussed previously, the BW perturbative calculation
of these two lowest $J^{\pi}$ states is converging with large enough $\xi$ due
to the variational principle. We see that the choice of $\xi=\langle{\rm
C}|H_{1}|{\rm C}\rangle=-130.227$ MeV is enough to make the calculations
converging to the NCSM results. Varying the value of $\xi$ is able to tune the
convergence speed of the BW perturbative calculation. Here we notice that as
$N_{\rm max}$ increases, the speed of convergence becomes lower and lower.
This may be explained as follows. The size of the $\mathbbm{Q}$-space
increases as we increase the size of the full model space (characterized by
$N_{\rm max}$), since we choose the $0\hbar\omega$ model space as the
$\mathbbm{P}$-space. Larger $\mathbbm{Q}$-space gives lower eigenvalues of
$QHQ$, and hence pushes the lowest eigenvalue characterized by $J^{\pi}$ of
$QHQ$ closer to the lowest $J^{\pi}$ eigenvalue of the Hamiltonian in the full
model space (i.e., pushes the lowest singularity of $f_{n}(E)$ closer to the
lowest intersection between $y=E$ and $y=f_{n}(E)$ for a given $J^{\pi}$).
Therefore the value of $\rho(R)$ at the lowest intersection gets larger but
cannot exceed one, which slows down the speed of convergence. If we continue
to increase $N_{\rm max}$, $\rho(R)<1$ still holds at the lowest intersection.
This is assured by the variational principle and is reflected by the feature
of negative slope of $f_{n}(E)$. The calculations with the bare N3LO potential
and $\xi=\langle{\rm C}|H_{1}|{\rm C}\rangle=-103.136$ MeV are plotted in Fig.
3(b). As we expected, even with the unsoftened N3LO potential, the BW
perturbative calculation of each lowest $J^{\pi}$ state can still get
converged with a proper $\xi$. Furthermore, the speed of convergence is even
faster than the case of Daejeon16, with the same choice of $\xi=\langle{\rm
C}|H_{1}|{\rm C}\rangle$. This is more obvious in the calculations of 6Li and
7Li, which will be given later.

###  III.2 6Li

The case of 6Li in the $p$-shell provides us with the first example of
$\mathbbm{P}$-space with $d_{p}>1$ for a given $J^{\pi}$. The eigenvalue
functions $f_{n}(E)$ of $H_{\rm eff}(E)$ for $J^{\pi}=0^{+},1^{+},2^{+},3^{+}$
calculated from the exact matrix inversion method using the Daejeon16
potential at $\hbar\omega=18$ MeV, $N_{\rm max}=2$ are plotted in Fig. 4. The
properties of these eigenvalues $f_{n}(E)$ are similar to the case of 5Li
discussed above, except that now we have multi-dimensional $\mathbbm{P}$-space
for $J^{\pi}=0^{+},1^{+},2^{+}$. We observe eigenvalue “crossings” for multi-
dimensional $\mathbbm{P}$-space in this figure, as marked with the gray
rectangle shadow, for example. These “crossings” are actually avoided level
crossings, as shown in Fig. 5 for the marked “crossing”, for example. The
intersections between the eigenvalue functions $y=f_{n}(E)$ and $y=E$ are the
solutions of the $\mathbbm{P}$-space eigenvalue problem Eq. (5). These
intersections can exactly reproduce the corresponding NCSM calculations. The
energies of $\mathbbm{P}$\- and $\mathbbm{Q}$-space dominated states are
marked with blue dots and red stars respectively. Correspondingly, the slopes
of $f_{n}(E)$ are small at blue dots and are large at red stars. From the
discussions in Sec. II.2, we expect that the perturbative calculation for
these marked states below the first singularities can get converged with large
enough $\xi$. As an example, Fig. 6 shows the eigenvalue functions
$f_{n}^{s{\rm th}}(E)$ of $H_{\rm eff}^{s{\rm th}}(E)$ for $J^{\pi}=0^{+}$
from the BW perturbative calculation up to various orders $s$. The
corresponding values of $\rho(R)$ at various energies are plotted in the top
panel. We see that the perturbative calculations in the convergence interval
$E<E_{\rm min}^{qhq}=4.621$ MeV (where $\rho(R)<1$) can indeed converge to the
results from the exact matrix inversion calculation. There are two
intersections residing in the convergence interval, which correspond to the
two lowest $0^{+}$ states and the BW perturbative calculation for these two
states can get converged. For the intersections (states) residing outside the
convergence interval, the BW perturbative calculation cannot get converged.

The BW perturbative calculations using the Daejeon16 potential for the states
($J^{\pi}=0^{+},1^{+},2^{+},3^{+}$) residing inside the convergence intervals
at $N_{\rm max}=2,4,6$ are plotted in Fig. 7. $\xi=\langle{\rm C}|H_{1}|{\rm
C}\rangle=-128.427$ MeV is used in the calculation. The corresponding NCSM
calculations are also performed to provide benchmarks. We see that the
convergence behavior of the perturbative calculations is similar to the case
of 5Li. We have now more than one state for $J^{\pi}=0^{+},1^{+},2^{+}$. The
BW perturbative calculation for the lowest state of each $J^{\pi}$ is
converging due to the variational principle. Although the perturbative
calculations of the excited states for each $J^{\pi}$ are all converging to
the NCSM results in this figure, the convergence of the excited state
calculation cannot be guaranteed for higher $N_{\rm max}$ or other potentials,
since an intruder state might occur and spoil the converging series for these
excited states. Figure 8 plots the similar calculations for 6Li with the bare
(unsoftened) N3LO potential and $\xi=\langle{\rm C}|H_{1}|{\rm
C}\rangle=-101.336$ MeV. We observe that all the calculations with the bare
N3LO converge to the NCSM results, and the speed of convergence is even faster
than in the case of calculations with the Daejeon16 potential with the same
choice of $\xi=\langle{\rm C}|H_{1}|{\rm C}\rangle$. We will also see this
feature in the case of 7Li in the next subsection.

###  III.3 7Li

Figure 9 presents the eigenvalue functions $f_{n}(E)$
($J^{\pi}=1/2^{-},3/2^{-},5/2^{-},7/2^{-}$) of $H_{\rm eff}(E)$ for 7Li
calculated by the exact matrix inversion method using the Daejeon16 potential
at $\hbar\omega=18$ MeV, $N_{\rm max}=2$. The result is similar to the case of
6Li. However, for $J^{\pi}=1/2^{-}$ only 5 intersections between $y=f_{n}(E)$
and $y=E$ are below the lowest singularity, whereas the $\mathbbm{P}$-space is
6-dimensional. We notice that there is an intruder state ($\mathbbm{Q}$-space
dominated) below the highest $\mathbbm{P}$-space dominated state. The presence
of this intruder state (or the slightly lower singularity) hinders the
convergence in the BW perturbative calculation, which will be seen later. The
corresponding BW perturbative calculation for $J^{\pi}=1/2^{-}$ is shown in
Fig. 10 with $\xi=\langle{\rm C}|H_{1}|{\rm C}\rangle=-127.141$ MeV. Indeed,
we observe that the BW perturbative calculations are converging to the NCSM
results expect for the highest $J^{\pi}=1/2^{-}$ state, at $N_{\rm max}=2$.
The difference between the result of MBPT and NCSM changes the sign at 60th
order for this highest $J^{\pi}=1/2^{-}$ state at $N_{\rm max}=2$. At $N_{\rm
max}=4$, this change of sign will be present at higher orders because of the
slower speed of convergence. We also observe that for a given $J^{\pi}$, as
the energy of the state increases, the convergence speed slows down until the
divergence is present, which corresponds to the fact that the value of
$\rho(R)$ becomes larger and larger and eventually exceeds 1. However there is
an anomaly for the $J_{k}^{\pi}=1/2_{3}^{-}$ state, whose speed of convergence
is faster than that of $J_{k}^{\pi}=1/2_{2}^{-}$. This is due to the
approximate isospin symmetry, namely $T\simeq 3/2$ for the
$J_{k}^{\pi}=1/2_{3}^{-}$ state, whereas $T\simeq 1/2$ for the other
$J^{\pi}=1/2^{-}$ states. The same perturbative calculation with the bare N3LO
potential is shown in Fig. 11. Surprisingly, the perturbative calculations of
all six $J^{\pi}=1/2^{-}$ states are quickly converging to the NCSM result,
compared to the case of Daejeon16. The corresponding exact eigenvalue
functions $f_{n}(E)$ ($J^{\pi}=1/2^{-}$) of $H_{\rm eff}(E)$ as a function of
$E$ is depicted in Fig. 12. From this figure, we see that the gap between the
$\mathbbm{P}$-space dominated states and the $\mathbbm{Q}$-space dominated
states (and also the lowest singularity) is pretty large and no intruder
states are present, in contrast to the case of Daejeon16. This may be
explained as follows. The SRG evolution makes the NCSM calculation converging
faster with respect to the size of the model space, and therefore the intruder
state can occur at lower energy even in a small model space. To illustrate
this idea, we plot in Fig. 13 seven lowest $J^{\pi}=1/2^{-}$ states of 7Li,
calculated with NCSM at $N_{\rm max}=2$, using the Daejeon16 and bare N3LO
potentials (upper and lower panels, respectively). As seen from the figure,
the lowest $\mathbbm{Q}$-space dominated state becomes an intruder state at
$\hbar\omega<20$ MeV using the Daejeon16 potential (it appears at lower
energies than the highest $\mathbbm{P}$-space dominated state). At the same
time a large gap between $\mathbbm{P}$-space dominated states and the lowest
$\mathbbm{Q}$-space dominated state is observed using the bare N3LO potential.
The presence of intruder states at low energy slows down the speed of, or even
prevents, the order-by-order convergence of the BW perturbative calculation.

##  IV Conclusions and Perspectives

In summary, we generalize the novel developments of BW MBPT for the ground
states of closed-shell nuclei given in Ref. [44] to ground and excited states
of open-shell nuclei. The general Hamiltonian partitioning scheme, the
convergence criterion, and the $\hat{K}$-box iterative method for the BW
perturbative calculation are introduced in the calculation of ground and
excited states of open-shell nuclei. Analytical derivations show that with a
large enough partitioning parameter $\xi$ (namely large enough diagonal
entries of $\xi$) the BW perturbation series can always be made convergent for
the states lower than the lowest eigenvalue of the Hamiltonian matrix in the
excluded $\mathbbm{Q}$-space for a given $J^{\pi}$. It follows that it is
always converging for each lowest $J^{\pi}$ state with large enough $\xi$, due
to the variational principle. This conclusion is independent of the choice of
the basis or the choice of the internucleon interaction.

To numerically check the above conclusions, we perform proof-of-principle
calculations for the open-shell nuclei 5,6,7Li using the (soft) Daejeon16
potential and the (hard) bare N3LO potential in the HO basis, benchmarked with
the results from the exact matrix inversion method and NCSM. We simply choose
$\xi=\langle{\rm C}|H_{1}|{\rm C}\rangle$ in these calculations. From these
calculations we can conclude:

  1. (a)

The BW perturbative calculation for each lowest $J^{\pi}$ state is always
converging with both soft and hard potentials in the HO basis with
$\xi=\langle{\rm C}|H_{1}|{\rm C}\rangle$. This conclusion is different from
the RS perturbative calculations given in Refs. [41, 42], where the
perturbation series is diverging in the HO basis even with soft potentials.

  2. (b)

The speed of convergence of BW perturbative calculation becomes lower as the
size of full model space increases, and for excited $J^{\pi}$ states the
convergent BW perturbation series can become divergent due to the presence of
intruder states, especially in large model space.

  3. (c)

The BW perturbative calculation with hard (bare N3LO) potential converges
faster than with soft (Daejeon16) potential using the same choice of
$\xi=\langle{\rm C}|H_{1}|{\rm C}\rangle$ at same $N_{\rm max}$.

With the developments in this paper, many works can be done in the future. For
example, with the iterative $\hat{K}$-box method, the energy levels and other
physical observables for other open-shell nuclei can also be calculated.
However, for open-shell nuclei not very close to the closed-shell, the
$\mathbbm{P}$-space dimension is large, which makes this kind of calculations
less efficient. To overcome this problem, we can use the calculations for one-
and two-valence nucleon systems to construct effective Hamiltonians used in
valence-space shell model calculations, so that the open-shell nuclei can be
calculated efficiently. Another problem is that the use of the $M$-scheme
many-body basis in the $\hat{K}$-box iterative calculation may prohibit
reaching heavier nuclei, similar to the case of NCSM. Instead of this
algebraic iterative $\hat{K}$-box method, we can formulate the $\hat{K}$-box
iterative equation (II.3) in a (Feynman) diagrammatic form so that particle-
hole excitations can be symmetrically truncated. This diagrammatic method will
be a brand new nonperturbative ab initio many-body method with the ability to
reach intermediate-mass and even heavy nuclei. The three-body forces can be
introduced in the algebraic iterative $\hat{K}$-box method in the way as NCSM
does. In the diagrammatic $\hat{K}$-box calculation, the normal-ordered
approximation up to two-body level for the three-body forces can be used.
Besides, the partitioning parameter $\xi$ can also be introduced to the RS
perturbation theory, especially in the diagrammatic formalism, to tune the
convergence behavior. These developments are under way and the results will be
published separately.

###### Acknowledgements.

The authors acknowledge the financial support from CNRS/IN2P3, France, via
ENFIA and ABI-CONFI Master projects. Large-scale computations have been
performed at MCIA, University of Bordeaux.

## References

  * Hergert [2020] H. Hergert, Front. Phys. 8, 379 (2020).
  * Epelbaum  _et al._ [2009] E. Epelbaum, H.-W. Hammer, and U.-G. Meißner, Rev. Mod. Phys. 81, 1773 (2009).
  * Machleidt and Entem [2011] R. Machleidt and D. Entem, Phys. Rep. 503, 1 (2011).
  * Bogner  _et al._ [2010] S. Bogner, R. Furnstahl, and A. Schwenk, Prog. Part. Nucl. Phys. 65, 94 (2010).
  * Faddeev [1961] L. D. Faddeev, Sov. Phys. JEPT 12, 1014 (1961).
  * Yakubowsky [1967] O. Yakubowsky, Sov. J. Nucl. Phys. 5, 937 (1967).
  * Ballot and Fabre de la Ripelle [1980] J. L. Ballot and M. Fabre de la Ripelle, Ann. Phys. (N.Y.) 127, 62 (1980).
  * Barnea  _et al._ [1999] N. Barnea, W. Leidemann, and G. Orlandini, Nucl. Phys. A 650, 427 (1999).
  * Pieper and Wiringa [2001] S. C. Pieper and R. B. Wiringa, Annu. Rev. Nucl. Part. Sci. 51, 53 (2001).
  * Barrett  _et al._ [2013] B. R. Barrett, P. Navrátil, and J. P. Vary, Prog. Part. Nucl. Phys. 69, 131 (2013).
  * Coraggio  _et al._ [2009] L. Coraggio, A. Covello, A. Gargano, N. Itaco, and T. Kuo, Prog. Part. Nucl. Phys. 62, 135 (2009).
  * Tichai  _et al._ [2020] A. Tichai, R. Roth, and T. Duguet, Front. Phys. 8, 164 (2020).
  * Hagen  _et al._ [2014] G. Hagen, T. Papenbrock, M. Hjorth-Jensen, and D. J. Dean, Rep. Prog. Phys. 77, 096302 (2014).
  * Hergert  _et al._ [2016] H. Hergert, S. Bogner, T. Morris, A. Schwenk, and K. Tsukiyama, Phys. Rep. 621, 165 (2016), memorial Volume in Honor of Gerald E. Brown.
  * Dickhoff and Barbieri [2004] W. Dickhoff and C. Barbieri, Prog. Part. Nucl. Phys. 52, 377 (2004).
  * Hu  _et al._ [2022] B. Hu, W. Jiang, T. Miyagi, Z. Sun, A. Ekström, C. Forssén, G. Hagen, J. D. Holt, T. Papenbrock, S. R. Stroberg, _et al._ , Nat. Phys. 18, 1196 (2022).
  * Brueckner [1955] K. A. Brueckner, Phys. Rev. 100, 36 (1955).
  * Brueckner and Levinson [1955] K. A. Brueckner and C. A. Levinson, Phys. Rev. 97, 1344 (1955).
  * Bethe [1956] H. A. Bethe, Phys. Rev. 103, 1353 (1956).
  * Goldstone [1957] J. Goldstone, Proc. Roy. Soc. A 239, 267 (1957).
  * Kuo and Osnes [1990] T. T. S. Kuo and E. Osnes,  _Folded-Diagram Theory of the Effective Interaction in Nuclei, Atoms and Molecules_ (Springer-Verlag, Berlin/Heidelberg, 1990).
  * Hjorth-Jensen  _et al._ [1995] M. Hjorth-Jensen, T. T. Kuo, and E. Osnes, Phys. Rep. 261, 125 (1995).
  * Coraggio  _et al._ [2023] L. Coraggio, G. De Gregorio, T. Fukui, A. Gargano, Y. Ma, Z. Cheng, and F. Xu, Prog. Part. Nucl. Phys. , 104079 (2023).
  * Kelly [1963] H. P. Kelly, Phys. Rev. 131, 684 (1963).
  * Kelly [1964a] H. P. Kelly, Phys. Rev. 134, A1450 (1964a).
  * Kelly [1964b] H. P. Kelly, Phys. Rev. 136, B896 (1964b).
  * Bartlett and Silver [1974a] R. J. Bartlett and D. M. Silver, Chem. Phys. Lett. 29, 199 (1974a).
  * Bartlett and Silver [1974b] R. J. Bartlett and D. M. Silver, Int. J. Quantum Chem. 8, 271 (1974b).
  * Bartlett [1981] R. J. Bartlett, Ann. Rev. Phys. Chem. 32, 359 (1981).
  * Shavitt and Bartlett [2009] I. Shavitt and R. Bartlett,  _Many-Body Methods in Chemistry and Physics: MBPT and Coupled-Cluster Theory_ (Cambridge University Press, Cambridge, 2009).
  * Kuo  _et al._ [1971] T. Kuo, S. Lee, and K. Ratcliff, Nucl. Phys. A 176, 65 (1971).
  * Ellis and Osnes [1977] P. J. Ellis and E. Osnes, Rev. Mod. Phys. 49, 777 (1977).
  * Coraggio  _et al._ [2003] L. Coraggio, N. Itaco, A. Covello, A. Gargano, and T. T. S. Kuo, Phys. Rev. C 68, 034320 (2003).
  * Coraggio  _et al._ [2005] L. Coraggio, A. Covello, A. Gargano, N. Itaco, T. T. S. Kuo, and R. Machleidt, Phys. Rev. C 71, 014307 (2005).
  * Roth  _et al._ [2006] R. Roth, P. Papakonstantinou, N. Paar, H. Hergert, T. Neff, and H. Feldmeier, Phys. Rev. C 73, 044312 (2006).
  * Günther  _et al._ [2010] A. Günther, R. Roth, H. Hergert, and S. Reinhardt, Phys. Rev. C 82, 024319 (2010).
  * Hu  _et al._ [2016] B. S. Hu, F. R. Xu, Z. H. Sun, J. P. Vary, and T. Li, Phys. Rev. C 94, 014303 (2016).
  * Löwdin [1965] P.-O. Löwdin, Phys. Rev. 139, A357 (1965).
  * Laidig  _et al._ [1985] W. D. Laidig, G. Fitzgerald, and R. J. Bartlett, Chem, Phys. Lett. 113, 151 (1985).
  * Perrine  _et al._ [2005] T. M. Perrine, R. K. Chaudhuri, and K. F. Freed, Int. J. Quantum Chem. 105, 18 (2005).
  * Roth and Langhammer [2010] R. Roth and J. Langhammer, Phys. Lett. B 683, 272 (2010).
  * Langhammer  _et al._ [2012] J. Langhammer, R. Roth, and C. Stumpf, Phys. Rev. C 86, 054315 (2012).
  * Tichai  _et al._ [2016] A. Tichai, J. Langhammer, S. Binder, and R. Roth, Phys. Lett. B 756, 283 (2016).
  * Li and Smirnova [2023] Z. Li and N. A. Smirnova, arXiv preprint arXiv:2306.13629 (2023).
  * Bloch and Horowitz [1958] C. Bloch and J. Horowitz, Nucl. Phys. 8, 91 (1958).
  * Löwdin [1962] P.-O. Löwdin, J. Math. Phys. 3, 969 (1962).
  * Horn and Johnson [2013] R. A. Horn and C. R. Johnson,  _Matrix Analysis_ , 2nd ed. (Cambridge University Press, 2013).
  * Schucan and Weidenmüller [1972] T. Schucan and H. Weidenmüller, Ann. Phys. (N.Y.) 73, 108 (1972).
  * Schucan and Weidenmüller [1973] T. H. Schucan and H. A. Weidenmüller, Ann. Phys. (N.Y.) 76, 483 (1973).
  * Shirokov  _et al._ [2016] A. Shirokov, I. Shin, Y. Kim, M. Sosonkina, P. Maris, and J. Vary, Phys. Lett. B 761, 87 (2016).
  * Entem and Machleidt [2003] D. R. Entem and R. Machleidt, Phys. Rev. C 68, 041001 (2003).
  * Entem and Machleidt [2002] D. Entem and R. Machleidt, Phys. Lett. B 524, 93 (2002).
  * Møller and Plesset [1934] C. Møller and M. S. Plesset, Phys. Rev. 46, 618 (1934).