# Many-Body Perturbation Theories for Finite Nuclei

^{1}Max-Planck-Institut für Kernphysik, Heidelberg, Germany^{2}Institut für Kernphysik, Technische Universität Darmstadt, Darmstadt, Germany^{3}ExtreMe Matter Institute EMMI, GSI Helmholtzzentrum für Schwerionenforschung GmbH, Darmstadt, Germany^{4}ESNT, Irfu, DPhN, CEA Saclay, Gif-sur-Yvette, France^{5}Irfu, DPhN, CEA Saclay, Gif-sure-Yvette, France^{6}Instituut voor Kern- en Stralingsfysica, KU Leuven, Leuven, Belgium

In recent years many-body perturbation theory encountered a renaissance in the field of *ab initio* nuclear structure theory. In various applications it was shown that perturbation theory, including novel flavors of it, constitutes a useful tool to describe atomic nuclei, either as a full-fledged many-body approach or as an auxiliary method to support more sophisticated non-perturbative many-body schemes. In this work the current status of many-body perturbation theory in the field of nuclear structure is discussed and novel results are provided that highlight its power as a efficient and yet accurate (pre-processing) approach to systematically investigate medium-mass nuclei. Eventually a new generation of chiral nuclear Hamiltonians is benchmarked using several state-of-the-art flavors of many-body perturbation theory.

## 1. Introduction

A major goal of quantum many-body theory is to provide accurate solutions of the stationary Schrödinger equation

for a given input Hamiltonian *H*, where $|{\Psi}_{k}^{\text{A}}\rangle $ denotes the *k*-th eigenstate of the *A*-body system with eigenvalue ${E}_{k}^{\text{A}}$. As for the study of the atomic nucleus at low energy, the starting point is a realistic Hamiltonian arising from the modeling of the strong interaction

in terms of nucleonic degrees of freedom. In Equation (2), *T* denotes the intrinsic kinetic energy, *V* the two-nucleon (2N) potential, *W* the three-nucleon (3N) potential, and so on. Nowadays, high-precision Hamiltonians are systematically constructed within the framework of chiral effective field theory (χEFT) [1–6]. Earlier on, more phenomenological models were employed that fitted a somewhat *ad hoc* parametrization to reproduce few-body observables, e.g., 2N scattering data. Examples of such potentials are the Argonne [7] or Bonn [8] potentials. It was recognized in various many-body studies that the inclusion of three-nucleon interactions is mandatory to reproduce nuclear phenomenology, e.g., nuclear saturation properties or the correct prediction of the oxygen neutron dripline [9, 10]. The relevance of three-nucleon interactions constitutes a key difference from other fields of many-body research like atomic physics, molecular structure, or solid-state physics, and poses significant challenges.

Over the past two decades, the range of applicability of *ab initio* nuclear many-body theory has been extended significantly. While 20 years ago first-principle solutions of the quantum many-body problem were restricted to nuclei lighter than *A* ≈ 12, formal and computational developments have made it possible to investigate a much wider range of masses. Initially, the *ab initio* treatment employed mostly large-scale diagonalization methods like the no-core shell model (NCSM) [11–13], or quantum Monte Carlo (QMC) [14–17] techniques. In addition, the few-body (A = 3, 4) solution could be constructed using the Fadeev approach or its Fadeev-Yakubowski extension [18, 19]. A major breakthrough occurred in the early 2000s when the re-import of so-called *expansion methods* from quantum chemistry provided systematically improvable many-body approximations for medium-mass closed-shell systems up to ^{132}Sn. In such approaches an initial guess for the exact wave function is taken as a reference state and corrections to this starting point are constructed through a chosen expansion scheme. Initially, this was done within the framework of self-consistent Green's function (SCGF) [20–25] and coupled cluster (CC) [26–31] theories that had proven to be extremely efficient at grasping *dynamical correlations* in electronic systems. Later on, the same strategy was transferred to other many-body expansion methods such as many-body perturbation theory (MBPT) [32–37], in-medium similarity renormalization group (IMSRG) [10, 38–44] or the unitary model operator approach (UMOA) [45, 46]. All of these methods provide a consistent description of ground-state energies of closed-shell mid-mass nuclei even though the rationales behind their expansions are not trivially related to one another. For sure, this consistency is a remarkable sign of success for *ab initio* nuclear many-body theory.

While closed-shell systems, dominated by so-called *dynamical correlations*, transparently allow for the use of single-reference techniques, the extension to open-shell systems requires a different strategy due to the degeneracy of single Slater-determinant reference states with respect to elementary excitations. Open-shell systems located in the vicinity of shell closures can be targeted via equation-of-motion (EOM) techniques where one or two nucleons are attached to the correlated ground state of a closed-shell nucleus [47, 48]. In nuclear systems, however, the strong coupling between spin and orbital angular momenta is such that long sequences of nuclei with open-shell character arise as the mass increases. Consequently, EOM approaches do not provide a viable option to tackle most of the open-shell systems that differ from closed-shell ones by more than one or two mass units. Two different routes have been followed in recent years to overcome this difficulty: (i) the construction of *ab initio*-rooted valence-space (VS) interactions used in a subsequent shell-model diagonalization and (ii) the use of correlated reference states capturing so-called *static correlations* and thus lifting from the outset the degeneracy with respect to elementary excitations. The design of VS interactions has been performed in various frameworks, going from simple low-order MBPT approaches [49] to more advanced non-perturbative schemes like IMSRG [39, 44] or CC [30]. While the design of the effective interaction can be performed at low polynomial cost, the final diagonalization, even though taking place in a limited valence space, still exhibits factorial scaling in the number of active nucleons pointing to the *hybrid scaling* of the approach.

In this article, the construction of VS interactions is not discussed and the focus is rather on the alternative strategy to overcome the limitations of single Slater-determinant-based expansions via the use of more general reference states. Reference states handling the bulk of *static correlations* from the outset re-introduce an *energy gap* in open-shell systems such that expanding the exact ground-state via elementary excitations of the reference state becomes well-defined again. In practice, this is done by resorting to either a *multi-determinantal* reference state or to a single *symmetry-broken* determinantal reference state. In an *ab initio* spirit, this was first done within the frame of SCGF theory, i.e., through the Gorkov extension of SCGF (GSCGF) [50] formulated on the basis of a Hartree-Fock-Bogoliubov (HFB) reference state breaking *U*(1) global-gauge symmetry associated with particle number conservation. Soon after, a multi-reference version of IMSRG (MR-IMSRG) was designed based on a particle-number-projected HFB state [51]. Around 2013, these two methods provided the first *ab initio* description of arbitrary mid-mass singly open-shell nuclei. Later on, expansion methods were merged with configuration interaction (CI) technology by using reference states from a prior NCSM calculation performed in a model space of limited size. In this way, one can systematically improve the many-body solution, either by increasing the size of the reference space or by relaxing the truncation of the many-body expansion. Within the framework of perturbative approaches this strategy yields multi-configurational perturbation theory (MCPT) [36] whereas for IMSRG it leads to the introduction of the in-medium no-core shell model (IM-NCSM) [52].

More recently, the use of particle-number-breaking reference states has been exploited in MBPT and CC theory, thus giving rise to Bogoliubov MBPT (BMBPT) and Bogoliubov CC (BCC). While BCC has only undergone proof-of-principle applications in limited model spaces so far [53], BMBPT has been applied successfully in large-scale applications up to medium-heavy isotopic chains [37]. Next, the additional or alternative breaking of *SU*(2) rotational symmetry associated with angular-momentum conservation will provide a systematic access to doubly open-shell nuclei [54]. While intensive efforts are already dedicated to this extension, no systematic result are available yet.

Whenever using a symmetry-broken reference state, an additional step must be envisioned to account for the quantum fluctuations eventually responsible for the lifting of the fictitious degeneracy associated with the broken symmetry. Indeed, the latter is only *emergent* in finite systems [55–57] such that it is mandatory to restore good symmetry quantum numbers. Doing so does not only change the energy (dramatically in certain situations) but also allows for the proper handling of transition operators characterized by symmetry selection rules. The formalism to achieve this step was proposed recently [54, 58, 59], but only applied to a schematic solvable model so far in the nuclear context [60] and will thus not be covered in the present article. Another topic of importance not presently covered relates to the benefit of applying so-called resummation methods to the Taylor expansion associated with perturbation theory. While traditional methods such as Padé resummation [32, 33, 61] can typically be employed with success, the newly formulated eigenvector continuation method [62, 63] was recently applied successfully to BMBPT [64] and show great promises for the future. Last but not least, and since the present document focuses on finite nuclei, the recent efforts made to describe nuclear matter, at zero and finite temperature, on the basis of MBPT are not discussed [65–68].

Eventually, the objective of the present article is to describe the on-going revival of MBPT, both in its basic form and in some of its novel extensions, to describe finite closed-shell and open-shell nuclei. Furthermore, the goal is to demonstrate how MBPT complement non-perturbative methods in two ways, i.e., (i) it can act as an inexpensive pre-processing method to accelerate non-perturbative techniques and (ii) as a post-processing tool to further improve upon non-perturbative many-body approaches.

While the present document discusses in detail many-body perturbation theory for (ground-state) energies, the discussion of other observables is intentionally left out. Targeting other observables via MBPT has never been investigated in depth in nuclear physics. Contrary to the energy that is evaluated projectively, other observables must be computed as the expectation value of the associated operator in the perturbatively-determined many-body state. Within CC theory, this is typically achieved via a linear-response treatment, giving rise to its so-called Λ-CC extension [69]. Employing perturbative CC amplitudes, the Λ-MBPT variant could be obtained to evaluate arbitrary observables. Most recently, various groups identified a resummation method called eigenvector continuation [70–73] (EC) as a powerful tool to robustly extract energies from the MBPT expansion, even in the case of a divergent perturbative expansion. Formally, the EC framework can be straightforwardly extended to deliver an expectation value estimate of arbitrary observables.

The present document is structured as follows. Before actually coming to perturbation theory, its possible sources of breakdown are discussed in section 2. The nuclear Hamiltonian and softening techniques are introduced in section 3. In section 4, formal perturbation theory is laid out. Section 5 is then dedicated to the standard Slater-determinant-based MBPT applicable to closed-shell systems. Multi-configurational perturbation theory and Bogoliubov many-body perturbation theory are discussed as open-shell extensions in sections 6 and 7, respectively. In the next two sections, MBPT is employed as a cheap and efficient *pre-processing* tool for non-perturbative many-body methods. In section 8 MBPT is used to pre-select important configurations, i.e., as a data compression tool, in a non-perturbative calculation whereas in section 9 it is used to pre-optimize single-particle states, thus accelerating the convergence of the non-perturbative calculation with respect to the one-body basis dimension. In section 10, MBPT is eventually employed as an inexpensive method to provide systematic tests of a newly designed family of χEFT nuclear Hamiltonians over a large set of nuclei. Conclusions and outlooks are provided in section 11.

## 2. Perturbative vs. Non-perturbative Problem

Before actually discussing MBPT, it is useful to consider the possible reasons for the failure of expansion methods. With these considerations at hand, various MBPT flavors can be better understood based on the interplay between the symmetries characterizing the reference state around which the exact eigenstate is expanded, its single- or multi-determinantal character as well as the resolution scale of the employed Hamiltonian.

### 2.1. Rationale

It has often been argued in the past that the nuclear many-body problem is “intrinsically non-perturbative” such that MBPT is bound to fail to describe correlations between nucleons. The statement is, in such generality, *not true* and the strong disfavor against MBPT techniques is, to a large extent, based on historically grown bias.

It is crucial to understand the basic fact that the “(non-)perturbative” character of a system characterized by a Hamiltonian *H* has no meaning in absolute and can only be stated with respect to a chosen starting point. This notion relates, at least implicitly, to an “unperturbed problem,” defined through an unperturbed Hamiltonian *H*_{0} and its associated eigenstates, with respect to which the targeted solution is meant to be expanded. If the expansion can be written as a converging powers series in *H*_{1} ≡ *H* − *H*_{0}, the problem is perturbative. Consequently, the more optimized *H*_{0}, the better the chances for the problem of actual interest to be *made* perturbative. For example, taking free particles as a reference point^{1}, the description of any bound state can only be *non* perturbative with respect to it. However, more optimized choices of *H*_{0}, making the description, perturbative might be accessible.

The critical point really resides in the cost required to find an appropriate *H*_{0} and its exact eigensolutions, i.e., if the cost to do so is similar to the one needed to employ a non-perturbative method, the perturbation theory built on top of it is not so appealing. In the end, the question is rather: can one find *H*_{0} and solve for its eigenstates at a moderate cost such that the eigenstates of *H* are obtained from them through a converging power series in *H*_{1} ≡ *H* − *H*_{0}? While success is certainly not guaranteed in general, the search for an optimal, yet simple enough, *H*_{0} must be performed in the most open-minded way. Typically, the statement that the nuclear many-body problem is “intrinsically non-perturbative” has been based on too restrictive assumptions of what *H*_{0} is allowed to be.

### 2.2. Ultra-Violet and Infra-Red Divergences

Whether a perturbative approach is viable or not certainly depends on the nature of the Hamiltonian, i.e., on the nature of the elementary degrees of freedom at play and of their interactions. As a matter of fact, two characteristics of the 2N interaction make the many-body problem hard to solve, e.g., *possibly* non-perturbative. The first one relates to the strong short-range central and tensor-forces between the nucleons inducing strong correlations in the ultra-violet (UV) regime. The second one relates to the large scattering lengths associated with the existence of a weakly bound proton-neutron state and of a virtual di-neutron state, which induce strong many-body correlations in the infra-red (IR) regime. These characteristics induce typical patterns, i.e., large matrix elements, of the residual interaction as schematically depicted in Figure 1: the upper-left corner corresponds to nuclear matrix elements between high-lying occupied and low-lying virtual states, i.e., to physics around the Fermi surface. Those matrix elements involve strong attractive IR couplings that fall off when moving further away from the Fermi surface. Complementary, repulsive UV couplings become more prominent if energetically higher virtual states are considered, independently of the particular pair of occupied states they interact with.

**Figure 1**. Schematic matrix representation of the residual nuclear two-body interaction at play in perturbation theory. The size of the matrix elements is proportional to the darkness of the pixel while the red (blue) color indicate schematically a negative (positive) sign. In the present case, typically associated to MBPT on top of a closed-shell Slater determinant, the IR (UV) couplings translate into large negative (positive) matrix elements between pairs of unoccupied and occupied states near the Fermi level (between any pairs of occupied and highly excited virtual states).

The potential occurrence of a UV-driven divergence is not nucleus specific and concerns all expansion methods, not just MBPT. Still, specific non-perturbative expansions may, at a given truncation order, resum a specific infinite subset of perturbation theory contributions that appropriately handles the divergence driven by the large UV couplings displayed in Figure 1. In the end, a pure power series in *H*_{1} is certainly the most sensitive expansion to strong low-to-high momentum couplings in the Hamiltonian that typically make the series to diverge after a few reasonable low-order contributions. It happens, however, that UV-driven divergences can be tamed to a large extent via a pre-processing of the Hamiltonian through renormalization group transformations. These transformations are briefly introduced in section 3.2 and their consequences on the MBPT expansion is illustrated in section 5.6.2.

The potential occurrence of a IR-driven divergence is nucleus specific but concerns all expansion methods, not just MBPT. Infra-red-driven divergences occur whenever the *A*-body unperturbed reference state, i.e., the ground state of *H*_{0}, is (nearly) degenerate with respect to elementary, e.g., particle-hole, excitations. Considering a standard Slater-determinant reference state, and as illustrated in the top panel of Figure 2, this situation occurs whenever the number of constituents is such that the highest occupied shell is only partially filled (is too close to the first empty shell), thus defining so-called open-shell (closed-sub-shell) nuclei. Combined with large IR matrix elements (see Figure 1), this cancelation (reduction) of the particle-hole gap makes the perturbative expansion (nearly) singular from the outset such that even the few first terms (e.g., Equation 45) are not well behaved. A numerical illustration of the emergence of such a IR divergence is provided in the bottom panel of Figure 2. This major difficulty can be controlled to a large extent via the use of more general classes of unperturbed Hamiltonians *H*_{0} and reference states than typically considered in the past. This idea and the corresponding results are discussed in sections 6 and 7.

**Figure 2. (Top)** Schematic representation of neutron or proton energy shells and associated occupations corresponding to a two-particle/two-hole excitation on top of the reference Slater determinant, i.e., the ground state of *H*_{0}, appropriate to a ^{16}O-like nucleus (N = Z = 8). The last occupied shell in the reference state is the *Fermi level* and its energy separation to the first empty level is denoted as Δ*E*_{F}. Left: closed-shell nucleus for which the number of nucleons is such that (i) the Fermi level is fully occupied and (ii) Δ*E*_{F} ≫ 0. Center: sub-closed shell nucleus for which the number of nucleons is such that (i) the Fermi level is fully occupied and (ii) Δ*E*_{F} is small. Right: open-shell nucleus for which the number of nucleons is such that the Fermi level is only partly occupied such that Δ*E*_{F} = 0. **(Bottom)** Emergence of an infra-red divergence in the MBPT expansion of the ground-state energy of ^{16}O induced by a step-wise reduction (going from blue, to yellow, to green, to purple, and to red) of the size of the particle-hole gap in the spectrum of *H*_{0}.

Of course, even with the use of renormalization group transformations of the Hamiltonian and rather general classes of reference states, the two characteristics of the nucleon-nucleon interaction in the UV and IR regimes may eventually compromise the convergence of *any* practical perturbative expansion and call for resummation techniques or the use of explicitly non-perturbative methods.

## 3. The Nuclear Hamiltonian

### 3.1. The Bare Operator

As briefly mentioned in the introduction, chiral effective field theory (χEFT) provides a convenient framework to construct systematically improvable nuclear Hamiltonians valid in the low-energy regime relevant to nuclear structure [4, 6]. Starting from nucleons and pions as explicit degrees of freedom, the long- and mid-range parts of the interaction are mediated by multiple-pion exchanges whereas the unresolved short-range part is modeled via contact terms and derivatives of contact terms. In the early 1990's Weinberg paved the way for a systematic treatment of the strong interaction by introducing a power-counting scheme stipulating the *a priori* importance of the infinite number of allowed contributions in the operator expansion [1–3]. Operators with higher particle rank naturally arise at higher orders in this scheme. Eventually, the parameters, i.e., the low-energy constants (LEC's), entering the operator expansion are fitted to low-energy experimental data [74, 75].

Eventually, the second-quantized form of the many-body Hamiltonian takes, in an arbitrary basis $\left\{\right|p\rangle \equiv {c}_{p}^{\u2020}|0\rangle \}$ of the one-body Hilbert space ${{H}}_{1}$, the form

The Hamiltonian is, thus, represented via a set of one-, two-, and three-body matrix elements *t*_{pq}, ${\stackrel{\u0304}{v}}_{pqrs}$ and ${\stackrel{\u0304}{w}}_{pqrstu}$, respectively. In a modern language the above matrix elements define *tensors* of mode *n* = 2, 4, 6, respectively, where the *mode* specifies the number of indices.

### 3.2. Similarity Renormalization Group

While the tensors defining the Hamiltonian built within χEFT may display large low-to-high momentum couplings, pre-processing tools can be used to tame them. During the past decade the (free-space) similarity renormalization group (SRG) approach has become the standard technique to generate a “softened” basis representation of an operator more amenable to many-body calculations [76].

The SRG approach is based on a unitary transformation of the initial operator *O* parameterized by a continuous parameter α ∈ ℝ, i.e.,

Equation (4) can be re-cast into a first-order differential equation

involving an anti-Hermitian *generator* η(α) that can be chosen freely to achieve a desired decoupling pattern in the transformed operator. A convenient choice employed in many calculations is given by

such that the SRG evolution can be interpreted as a *pre-diagonalization* of the operator in momentum space, thus suppressing the coupling between high- and low-momentum modes. This procedure thus drives the Hamiltonian toward a band-diagonal form. Writing *H*(α) ≡ *T* + *V*(α) + *W*(α) + … in the same single-particle basis as the starting Hamiltonian, the SRG transformation corresponds to generating α-dependent tensors ${\stackrel{\u0304}{v}}_{pqrs}(\alpha )$, ${\stackrel{\u0304}{w}}_{pqrstu}(\alpha )\dots $ whose UV elements linking single-particle states corresponding to low and high momenta are strongly suppressed.

In many-body applications SRG-evolved operators display highly improved model-space convergence, thus facilitating studies of mid-mass nuclei. The impact on the convergence properties of the MBPT series will be illustrated in section 5.6.2. However, the numerical improvements come at the price of induced many-body operators, i.e., the unitary transformation shifts information to operators with higher particle ranks. For instance, employing an initial two-body operator *O*^{2B} leads to

In practice, Equation (7) must be truncated at a given operator rank, thus discarding higher-body operators. This approximation formally violates the unitarity of the transformation in Fock space and eventually induces a dependence of many-body observables on the SRG parameter α. A reasonable trade-off must be found for the value of α employed, i.e., it must improve the model-space convergence while keeping the effect of induced many-body operators at a minimum. The optimal parameter range may vary depending on the operator one starts from.

For the evaluation of nuclear properties, it is crucial to consistently transform all other operator representations to the same resolution scale as the Hamiltonian in order to provide a meaningful comparison in terms of a consistent operator basis smoothly connected in terms of the RG flow.

### 3.3. The “Standard” Hamiltonian

All many-body applications discussed below, except for the novel ones presented in section 10, employ a chiral Hamiltonian containing a 2N interaction at next-to-next-to-next-to-leading-order (N3LO) with a cutoff value of Λ_{2N} = 500 MeV/*c* [77, 78]. Three-body forces are included up to next-to-next-to-leading order (N2LO) with a local regulator [78] based on a cutoff value of Λ_{3N} = 400 MeV [79]. This constitutes a “standard” Hamiltonian used in many recent *ab initio* studies of light and medium-mass nuclei.

Additionally, the intrinsic Hamiltonian is consistently SRG-evolved in the two- and three-body sectors [80, 81]. The particular value of the SRG parameter is specified in each individual application. To avoid the complication of dealing with genuine three-body operators various forms of so-called normal-ordered two-body approximations (NO2B) are employed, depending on the particular nature of the *A*-body reference state [79, 82, 83].

## 4. Formal Perturbation Theory

The presentation of perturbation theory can be separated into *formal perturbation theory* and *many-body perturbation theory* [69]. Formal perturbation theory allows one to understand the general rationale and most relevant properties of the formalism. This is done by employing abstract Dirac notations and by specifying the initial assumptions via the action of Hilbert or Fock space operators on basis vectors. In particular, many key results can be obtained without specifying the content of the Hamiltonian (e.g., the rank of the operators it contains), the nature of the partitioning (e.g., the symmetries characterizing each contribution) and the associated reference state.

### 4.1. Partitioning

The starting point of perturbation theory relates to a *partitioning of the Hamiltonian*

into an unperturbed part *H*_{0} and a perturbation *H*_{1} ≡ *H* − *H*_{0}. The main assumption relies on the fact that the eigenvalue equation for *H*_{0} is numerically accessible, i.e.,

delivering the set of unperturbed eigenstates and eigenergies $\left\{\right|{\Phi}_{k}\rangle \text{},{E}_{k}^{(0)};k\in \mathbb{N}\}$ making up an orthonormal, i.e.,

basis of the many-body Hilbert space.

**Remark:**A large part of this document is dedicated to the description of nuclear ground states, i.e., *k* = 0. Consequently, the corresponding index is dropped in the following whenever targeting the ground state, e.g., $|{\Psi}_{0}^{\text{A}}\rangle =|{\Psi}^{\text{A}}\rangle ,|{\Phi}_{0}\rangle =|\Phi \rangle \text{or}{E}_{0}^{(0)}={E}^{(0)}$.

One typically employs *intermediate normalization*, i.e., the ground state |Ψ^{A}〉 of *H* is connected^{2} to the unperturbed ground-state |Φ〉 of *H*_{0} such that

Associated with the above partitioning are the *projection operators*

where *P*|Φ〉 = |Φ〉 and *Q*|Φ〉 = 0 by orthonormality. It can be shown that *P* and *Q* do meet the requirements of projection operators, i.e., Hermiticity and idempotency [69]. The operator *Q* can be explicitly written as

where the primed sum indicates the exclusion of the reference state from the summation. With these operators at hand, the exact ground-state can be written as

where the *correlated part* |χ〉 ≡ *Q*|Ψ^{A}〉, which is the unknown to be solved for, denotes the *orthogonal complement* of |Φ〉.

Eventually, the exact ground-state energy is typically accessed in a *projective* way^{3} by left-multiplying Equation (1) with the reference state 〈Φ| such that

where ${E}_{\text{ref}}\equiv \langle \Phi \left|H\right|\Phi \rangle ={E}^{(0)}+\langle \Phi \left|{H}_{1}\right|\Phi \rangle $ and $\Delta E\equiv {E}^{\text{A}}-{E}_{\text{ref}}=\langle \Phi \left|{H}_{1}\right|\chi \rangle $ denote reference and correlation energies, respectively. When using a reference state of product type, e.g., a Slater determinant, Δ*E* accounts for correlations between the nucleons beyond the mean-field approximation.

### 4.2. Resolvent Operator

The complete derivation of formal perturbation theory is best performed in terms of the (Rayleigh-Schrödinger) many-body *resolvent operator*

which, due to orthonormality of the employed many-body basis, annihilates the reference state

It is possible to employ alternative choices, such as the *Brillouin-Wigner* resolvent

which differs from *R*^{RS} by the presence of the *exact energy* in the denominator instead of the unperturbed energy *E*^{(0)}. In practice, Brillouin-Wigner perturbation theory requires an (computationally intensive) iterative solution and, additionally, suffers from a lack of size-extensivity^{4}. Therefore, this choice is only scarcely used in many-body applications. All of the subsequent results are obtained using a Rayleigh-Schrödinger resolvent. Consequently, the upper-case label “RS” is dropped to avoid notational clutter.

### 4.3. Power-Series Expansion

After a long but straightforward derivation [69], one obtains the correlated part of many-body ground-state and associated energy under the form [84–86]

The lower index “c” stipulates the *connected* character of the expansion ensuring its size-extensivity, i.e., proper scaling of observables with system size^{5}. Combining Equations (15) and (19b), one obtains the total ground-state wave-function and energy as power series in *H*_{1}

such that |Ψ^{(0)}〉 = |Φ〉 and ${E}_{\text{ref}}={E}^{(0)}+{E}^{(1)}$, i.e., the first non-trivial correction contributing to the ground-state correlation energy corresponds to the second-order term of the power series. It reads as

and can be re-written more explicitly by expanding the resolvent as

Equation (23) provides a prototypical example of a PT expression associated with ground-state energy corrections involving a resolvent operator connecting the unperturbed reference state (i.e., the *P* space) to excited states of *H*_{0} (i.e., the *Q*-space), and then going back to the reference state through the perturbation *H*_{1}. As will be seen with explicit MBPT, the nature of the elementary excitations of the reference state effectively involved at a given order depend on the rank, i.e., the *k*-body character, of the perturbation *H*_{1}.

### 4.4. Recursive Formulation

Equations (19a) and (19b) conveniently provide explicit expressions for the energy and state corrections whenever working at rather low orders. To go to high orders and study the convergence properties of perturbation theory as a power series, a different scheme becomes more useful. It relates to (i) making more explicit that the perturbative expansion relates to the power series expansion of a mathematical function taken at a particular value of its variable and to (ii) computing the coefficient of the series in a recursive way.

On the basis of the partitioning introduced in Equation (8), one defines a one-parameter family of Hamiltonians

such that *H*(0) = *H*_{0} and *H*(1) = *H*. Perturbation theory assumes that the exact eigenstates and eigenenergies of *H*(λ) can be parameterized through the *power series ansatz*

Setting λ = 0, the problem becomes equivalent to Equation (9) while setting λ = 1, one recovers the expansions of Equations (20) and (21) associated with the fully interacting problem.

Inserting the power series ansatz into the stationary Schrödinger equation for *H*(λ) and grouping together the terms proportional to λ^{p} leads to

Left multiplying Equation (26) with 〈Φ| and using intermediate normalization yields

which allows one to write the *p*-order ground-state energy correction as

Left multiplying Equation (26) with 〈Φ_{m}|, *m* ≠ 0, and matching the terms proportional to λ^{p} provides the relation

Introducing the coefficients

allows one to expand the *p*-order ground-state correction |Ψ^{(p)}〉 on the unperturbed basis {|Φ_{m}〉} according to

such that Equation (28) becomes

Inserting Equation (31) into Equation (30) further provides a recursive scheme to compute ${C}_{m0}^{(p)}$

with the initial condition ${C}_{m0}^{(0)}={\delta}_{m0}$.

Eventually, Equations (32) and (33) form a set of *recursive relations* from which the ground-state energy and state corrections can be obtained to all orders. In practice, the eigenbasis of *H*_{0} and the associated matrix elements of *H*_{1} are built such that the latter is stored in memory. The recursive steps are then identified as large-scale matrix-vector multiplications, thus, using the same technology as configuration-interaction approaches like the NCSM. While being formally convenient and obviating the explicit algebraic computation of order-*p* corrections, this approach is limited by the storage of the Hamiltonian due to the extensive size of the many-body basis. Therefore, only proof-of-principle studies in model spaces of limited dimensionality can be performed and realistic calculations of mid-mass nuclei are out of reach in this way.

## 5. Closed-Shell Many-Body Perturbation Theory

With the results of formal perturbation theory at hand, one can envision to apply them to specific many-body systems. To do so, one must further specify the nature of the partitioning and of the associated reference state. In particular, the maximum rank and symmetries of *H*_{0} must be characterized. Additionally, the goal of MBPT is to express all quantities, e.g., many-body matrix elements and unperturbed eigenenergies, entering the formulae at play in terms of the actual inputs to the many-body problem, i.e., the mode-2*k* tensors defining the *k*-body contributions to the Hamiltonian (Equation 3).

A series of tools exists to compute the expectation value of products of (many) operators in a vacuum state in an incrementally faster, more flexible and less error-prone way. The first step in this series corresponds to using the second-quantized representation of many-body operators in a chosen single-particle basis and to performing canonical commutations of fermionic operators. Next comes Wick's theorem [87], which is nothing but a procedure to capture the result in a condensed and systematic fashion. Still, the combinatorics associated with the application of Wick's theorem quickly becomes cumbersome whenever a long string of creation and annihilation operators is involved. Furthermore, many terms thus generated give identical contributions to the end result. Many-body diagrams address this issue [69] by providing a pictorial representation of the contributions and, even more importantly, by capturing at once all identical contributions, thus reducing the combinatorics tremendously. While being incredibly useful, the number of diagrams itself grows tremendously when applying MBPT beyond the lowest orders, thus leading to yet another combinatorial challenge. This translates into the difficulty to both generate all allowed diagrams at a given order without missing any and to evaluate their expression in a quick and error-safe way. Consequently, the last tool introduced to tackle this difficulty consists of an automatized generation and evaluation of diagrams [88–95]. All these technical, yet crucial, aspects of MBPT are not addressed in the present article and the interested reader is referred to the references.

### 5.1. Reference State

The present chapter is dedicated to the simplest form of MBPT appropriate to closed-shell systems. This first version relies on the use of a symmetry-conserving Slater determinant reference state

where the set of single-particle creation operators $\left\{{c}_{p}^{\u2020}\right\}$ acts on the physical vacuum |0〉. This constitutes an appropriate starting point of the perturbative expansion as long as |Φ〉 denotes a closed-shell Slater determinant in agreement with the left-hand case in Figure 2. While, in principle, the single-particle basis is completely arbitrary, applications will reveal its significant impact on the qualitative behavior of the perturbative expansion.

### 5.2. Normal Ordering

Applying Wick's theorem with respect to |Φ〉, the Hamiltonian can be rewritten in terms of normal-ordered contributions

where :: denotes the normal order of the involved creation and annihilation operators. Thus, *H*^{[0]} is the expectation value of *H* in |Φ〉 whereas ${H}_{pq}^{\left[2\right]}$ and ${H}_{pqrs}^{\left[4\right]}$ define matrix elements of *effective*, i.e., normal-ordered, one-body and two-body operators, respectively. The dots denote normal-ordered operators of higher ranks, up to the maximum rank *k*_{max} characterizing the initial Hamiltonian (Equation 3). Through the application of Wick's theorem, an effective operator of rank *k*_{eff} receives contributions from all initial operators with rank *k*, where *k*_{eff} ≤ *k* ≤ *k*_{max}. Using an initial Hamiltonian with up to three-nucleon interactions and working in the normal-ordered two-body approximation (NO2B) [79, 82], the residual three-body part *H*^{[6]} is presently discarded. For explicit expressions of the matrix elements defining the normal-ordered operator (see [79, 82, 83]).

### 5.3. Partitioning

To explicitly set up the partitioning of the Hamiltonian (Equation 8), one adds and subtracts a *diagonal* normal-ordered one-body operator

such that

with

Introducing the set of Slater determinants obtained from |Φ〉 via *n*-particle/*n*-hole excitations

one obtains an orthonormal basis of the *A*-body Hilbert space

which is nothing but the eigenbasis of *H*_{0}

where

sums (subtracts) the *n* one-body energies of the particle (hole) states the nucleons are excited into (from). Equation (41) corresponds to the explicit form of Equation (9) in the case of a Slater determinant reference state.

**Convention:** One-body states occupied (unoccupied) in the reference determinant are labeled by *i, j, k*, ... (*a, b, c*, ...) and are referred to as *hole* (*particle*) states. Generic one-body states are denoted by *p, q, r*, ....

The single-particle energies {ϵ_{p}} are parameters of the theory that are fixed by the partitioning, which itself defines the reference state. They can be chosen arbitrarily as long as the *A* occupied hole states have lower energies than the remaining particle states, such that ${\u03f5}_{ij\cdots}^{ab\cdots}>0$. A simple choice employed in nuclear physics consists of building |Φ〉 by filling up the *A* lowest single-particle eigenstates of the spherical *harmonic oscillator* Hamiltonian [32, 33], i.e., setting

where the oscillator frequency ω specifies the width of the potential. A more standard choice throughout various fields of many-body physics and chemistry relates to the so-called Møller-Plesset partitioning that corresponds to taking ${\stackrel{\u02d8}{H}}^{[2]}=0$, i.e., ${H}_{1}={H}^{\left[4\right]}$. This is obtained by using the reference Slater determinant |Φ〉 solution of the Hartree-Fock (HF) variational problem and by defining ${\stackrel{\u0304}{H}}^{\left[2\right]}$ from the eigenvalues of the one-body HF Hamiltonian.

### 5.4. Perturbative Expansion

Given the eigenbasis of *H*_{0} characterized by Equations (34), (39), and (41), the many-body resolvent (Equation 16) takes the form

and is to be fed into Equation (19b) that, once truncated at a given power in *H*_{1}, provides the correlation energy at the corresponding perturbative order.

### 5.5. Low-Order Formulas

As alluded to above, the evaluation of low-order corrections is facilitated by representing the MBPT expansion diagrammatically. This is typically done using either Hugenholtz or (anti-)symmetrized Goldstone diagrams, i.e., the time-ordered counterpart of Feynman diagrams that are used to compute matrix elements in quantum field theory. The interested reader is referred to the literature, e.g., Shavitt and Bartlett [69], for an elaborate discussion of the diagrammatic rules and their relation to Wick's theorem.

Focusing on the first non-trivial correction to the reference energy (Equation 23), the second-order correction takes the algebraic form

and is, thus, expressed in terms of the tensors defining the residual interaction *H*_{1} (Equation 37b) in normal-ordered form. The first contribution in Equation (45) relates to a so-called *non-canonical diagram* that vanishes if the reference state is taken to be the HF Slater determinant. The second term constitutes the genuine and dominant second-order correction that contributes for any Slater determinant reference state. Using the HF Slater determinant reference state has the practical benefit of lowering the number of many-body diagrams to be considered. While this feature is not relevant at second-order, the proliferation of non-canonical diagrams at higher order [95] makes the writing of numerical codes more cumbersome. Still, at a given order non-canonical diagrams are always of sub-leading complexity from a computational point of view, i.e., they involve fewer single-particle summations, such that they do not drive the computational cost.

Since *E*^{(2)} provides the leading contribution to the perturbative expansion, one observes that dynamical correlations are dominated by low-lying 2p2h-contributions. Most importantly, it is clear from Equation (45) that the second-order correction is manifestly negative, i.e., it increases the binding energy. This stems from the fact that the numerators are squared norms of matrix elements contributing to *H*_{1} and that the denominators are positive as long as the Slater determinant reference state displays a non-zero shell gap between occupied and unoccupied states, i.e., as long as one deals with a closed-shell nucleus.

### 5.6. Results

The first goal of the present analysis is to study the convergence characteristics of the perturbative expansion. In absence of analytical knowledge, this study must be based on empirical observations of high-order corrections, which is achieved through the recursive formulation of section 4.4 in small model spaces. Following this analysis, results of low-order MBPT calculations in realistic model spaces are presented to illustrate state-of-the-art *ab initio* applications to doubly closed-shell nuclei [34].

#### 5.6.1. Impact of Partitioning

While perturbation theory defines a general framework to access nuclear observables, the performance strongly depends on the choice of the partitioning *H* = *H*_{0} + *H*_{1} or, equivalently, on the underlying vacuum fixing the starting point for the expansion. Subsequently, two choices for *H*_{0} are presently compared in the calculation of the ground-state energy of ^{16}O, i.e., the one-body (i) spherical harmonic oscillator (HO) and (ii) self-consistent HF^{6} Hamiltonians (cf. section 5.3). The model space is truncated employing the *N*_{max}-truncation similar to the NCSM. Figure 3A shows the sequence of partials sums using a HO partitioning for a set of model spaces. The partial sums are divergent in all cases, which can equally be seen from the exponential divergence of high-order energy corrections in Figure 3C. On the other hand, using a HF reference state yields a rapidly converging perturbation series (Figure 3B) and the energy corrections are exponentially suppressed as a function of the perturbative order (Figure 3D), indicating robust convergence. In all cases the converged results agree up to numerical accuracy with the exact CI diagonalization.

**Figure 3**. Partial sums **(A,B)** and order-by-order corrections **(C,D)** for the ground-state energy of ^{16}O using the Hamiltonian described in section 3. Top (bottom) panels correspond to the HO (HF) partitioning. All calculations are performed using *N*_{max} = 2, 4, 6 (, , ) and an oscillator frequency of ℏω = 24 MeV. The SRG parameter is set to α = 0.08 fm^{4}. Reproduced from Tichai et al. [34] under the Creative Commons CC-BY license.

Obviously, the reference state heavily affects the performance of MBPT. In the above case, this can be understood by the poor quality of the HO reference, e.g., the wrong asymptotic radial dependence of single-particle HO eigenstates (Gaussian instead of exponential suppression). Consequently, in the following a HF determinant is used as a reference state when results are reported for closed-shell nuclei.

#### 5.6.2. SRG Dependence

Using a HF partitioning, the impact of the SRG transformation of the Hamiltonian on the perturbative series is now illustrated. In Figure 4, the ground-state energy of ^{4}He, ^{16}O, and ^{24}O is displayed while varying the value of the flow parameter α defining the SRG transformation. The left-hand panels show that the perturbative series converge in all cases thus demonstrating the reliability of HF-MBPT. For the light ^{4}He, the results are independent of the flow parameter and the MBPT expansion converges rapidly in all cases. For ^{16}O and ^{24}O, the rate of convergence is slower for harder interactions, i.e., for lower values of α. Furthermore, the partial sums admit a damped oscillatory behavior in the oxygen isotopes for α = 0.02 fm^{4}. These features can be better seen from the right-hand panels, where lower values of α induce a slower suppression of higher-order corrections (Figures 4E,F).

**Figure 4**. Partial sums **(A–C)** and order-by-order corrections **(D–F)** corresponding to the ground-state energy obtained from the Hamiltonian described in section 3. Top, middle, and bottom panels correspond to ^{4}He (*N*_{max} = 6), ^{16}O (*N*_{max} = 6), and ^{24}O (*N*_{max} = 4), respectively. All calculations are performed setting the oscillator frequency to ℏω = 20 MeV. The different plot markers correspond to SRG flow parameter α = 0.02(), 0.0625(), and 0.08 fm^{4}(). Reproduced from Tichai et al. [34] under the Creative Commons CC-BY license.

#### 5.6.3. Realistic Calculations

The previous results reveal that using an optimized HF reference state combined with a sufficiently soft interaction defines a well-controlled regime where perturbation theory can be robustly applied to closed-shell nuclei. For heavier nuclei and larger model spaces the high-order perturbative series cannot be computed recursively and, rather, low-order expressions are evaluated via explicit single-particle summations^{7}.

In Figure 5, the ground-state energy per particle (top panel) and the correlation energy per particle (bottom panel) of a selection of doubly closed-shell nuclei ranging from ^{4}He to ^{132}Sn is displayed at second- and third-order in MBPT [34]. A model space built out of 13 major harmonic oscillator shells is employed. Since the target nuclei are out of reach of exact diagonalization, MBPT results are compared to state-of-the CC calculations employing the same input Hamiltonian and the same HF determinant reference state.

**Figure 5**. Ground-state binding energies per particle **(A)** and correlation energy per particle **(B)** of closed-shell nuclei from second-order MBPT (), third-order MBPT (), CCSD (), and CR-CC(2,3) () using the Hamiltonian described in section 3. All calculations are performed using 13 oscillator shells and an oscillator frequency of ℏω = 20 MeV. The SRG parameter is set to α = 0.08 fm^{4}. Reproduced from Tichai et al. [34] under the Creative Commons CC-BY license.

The top panel demonstrates that third-order calculations fully capture the bulk part of the ground-state energy and are in remarkable agreement with more sophisticated non-perturbative CC results, i.e., the deviation with CR-CC(2,3) is less than 1% in all cases. A more refined analysis can be deduced from the bottom panel where the mean-field binding energy has been subtracted. The correlation energy accounts in most mid-mass systems for about 1 − 2 MeV/*A* such that the reference HF results capture roughly 60 − 70% of the overall binding energy. The CCSD correlation energy lies between second- and third-order results even though the CCSD wave function resums correlation effects beyond third order, thus indicating a repulsive effect on the binding from 2-particle/2-hole-excitations at 4th order and beyond. When (approximately) including triple excitations through CR-CC(2,3), slightly stronger binding than in third-order MBPT is generated.

Of course, the enormous benefit of low-order MBPT is that it excellently reproduces highly sophisticated CC results at a computational cost that is two orders of magnitude lower due to its non-iterative character.

## 6. Multi-Configurational Perturbation Theory

Driven by the capacity of HF-MBPT to grasp dynamical correlations induced in closed-shell systems by soft chiral Hamiltonians, extensions to genuine open-shell systems were envisioned. As alluded to in section 2, the presence of degeneracies with respect to particle-hole excitations does not allow the use of a single symmetry-conserving Slater-determinant reference state. A first possibility is to use a multi-determinantal reference state, which was shown to be very useful in electronic structure applications [96, 97] and as a pre-processing tool in NCSM calculations [98].

### 6.1. Rationale

The starting point of multi-configurational perturbation theory (MCPT) is the definition of an initial reference space

built from a set of orthonormal many-body Slater determinants |Φ_{ν}〉. The extended reference space redefines the nature of the *P*-space (and thus of the *Q*-space) introduced in Equation (12). The multi-configurational reference state |Ψ_{ref}〉 is chosen to be a normalized vector obtained from a diagonalization in ${{M}}_{\text{ref}}$, i.e.,

where *c*_{ν} denotes the expansion coefficients, typically obtained from a NCSM calculation. The initial diagonalization provides a set of *non-degenerate* but multi-determinantal reference states carrying good symmetry quantum numbers. At the price of giving up the product-type character of the reference state, the degeneracy is lifted and a well-defined perturbative expansion can be designed [36].

**Convention:** The symbol |Ψ〉 is used to emphasize the multi-determinantal character of the reference state and distinguish it from a single product state. Furthermore, two different notations are employed to designate the Slater determinants spanning the complete Hilbert space: (i) Slater determinants belonging to ${{M}}_{\text{ref}}$ are denoted by |Φ_{ν}〉, i.e., as a capital Greek letter carrying a lower-case Greek index, whereas (ii) Slater determinants outside ${{M}}_{\text{ref}}$ are denoted by |ϕ_{i}〉, i.e., as a lower-case Greek letter carrying a lower-case Roman index.

It is worth noting that MCPT naturally accesses excited states by building the perturbation theory on top of the various vectors produced through the prior diagonalization in ${{M}}_{\text{ref}}$. Of course, it is not guaranteed that the energetically lowest state in the initial NCSM calculation eventually corresponds to the ground state in the fully correlated limit, i.e., perturbative corrections may induce *level crossings* among the various states.

### 6.2. Partitioning

Formally, the unperturbed Hamiltonian is written in the spectral representation as

where |Ψ_{k}〉 denote the NCSM eigenvectors within ${{M}}_{\text{ref}}$, including the particular one, e.g., the lowest state, playing the role of the reference state. Consequently, one obtains the following eigenvalue relation for the unperturbed Hamiltonian

since the set of determinants $\left\{\right|{\Phi}_{i}\rangle \notin {{M}}_{\text{ref}}\}$ is orthogonal to |Ψ_{k}〉. In principle, the construction of *H*_{0} requires a *full diagonalization* in the reference space, i.e., the solution of all eigenvectors and eigenvalues making the construction of the unperturbed solution rapidly unfeasible if the dimension of the reference state grows significantly. As will become clear below, the computation of the lowest-order correction only require to access the reference state, which is thus the only the eigenvector that needs to be solved for explicitly in this case. In that case one may resort to Lanczos algorithms, thus targeting a limited number of extremal eigenstates.

Zeroth-order energies ${E}_{i}^{(0)}$ of the unperturbed Slater determinants making up the *Q* space, i.e., $\left\{\right|{\Phi}_{i}\rangle \notin {{M}}_{\text{ref}}\}$, are given by the sum of occupied single-particle energies defined as diagonal matrix elements of the one-body Hamiltonian

where the one-body density matrix of the reference state^{8} is introduced

In principle, an explicit three-body term can be included as well at the price of invoking the two-body density matrix. However, for the sake of computational simplicity a normal-ordered two-body (NO2B) approximation is employed to approximately account for the inclusion of 3N interactions [82].

The zeroth-order energy of the reference state is also defined via the single-particle energies defined in Equation (50) while taking into account the multi-determinantal character of the reference state through the mean occupation of single-particle states, i.e., the diagonal elements of the one-body density matrix ${\rho}_{pp}^{(0)}$, so that

### 6.3. Low Orders

With the partitioning of the Hamiltonian defined above, zeroth- and first-order MCPT contributions to the energy read as

such that their sum reproduces the full reference energy ${E}_{\text{ref}}^{\text{NCSM}}$ obtained via the diagonalization of the full Hamiltonian *H* in ${{M}}_{\text{ref}}$.

The second-order energy correction reads similarly to the one at play in standard MBPT, i.e.,

where the sum runs over states outside of the reference space and where the contribution from *H*_{0} vanishes by orthogonality 〈Ψ_{ref}|Φ_{i}〉 = 0. To explicitly evaluate *E*^{(2)} the reference state is expanded according to Equation (47)

All many-body matrix elements appearing in the algebraic expressions of the perturbative corrections involve Slater determinants only and can be readily evaluated using standard NCSM technology. As an efficient alternative, normal-ordering techniques and standard Wick's theorem are employed such that an associated diagrammatic can be designed. It is worth noting that intermediate states from within ${{M}}_{\text{ref}}$ only start contributing at fourth order [36] such that they do not appear in the evaluation of *E*^{(2)}. In Equation (56), the Hamiltonian is normal ordered with respect to the rightmost determinant $|{\Phi}_{{\mu}^{\prime}}\rangle $ for each term in the sum over μ′ and the two matrix elements are evaluated using the associated Wick's theorem. Similar techniques have been applied in quantum chemistry [101, 102]. The computational scaling of the second-order correction for large reference spaces is given by $\text{dim}{({{M}}_{\text{ref}})}^{2}\xb7{n}_{p}^{2}\xb7{n}_{h}$, where *n*_{p} and *n*_{h} denote the number of particle and hole states, respectively.

### 6.4. Results

In the following the performance of MCPT, specifically denoted as NCSM-PT in the present case, is gauged in a similar spirit as for HF-MBPT in section 5.

#### 6.4.1. High-Order Corrections in Light Systems

The recursive treatment of HF-MBPT laid out in section 4.4 has proven invaluable to understand the convergence characteristic of the perturbative expansion. While NCSM-PT does not employ a Slater-determinant reference, the recursive formulation can be extended in a straightforward way [36]. Figure 6 displays the convergence behavior of the perturbative series for^{6}^{,}^{7} Li built on top of the four lowest states of the *N*_{max} = 0 NCSM diagonalization. The left-hand panels show that the perturbative series is convergent for both systems and all target states with a slight overbinding of the second-order partial sum for many states. In all cases the converged results agree with exact diagonalization. Furthermore, right panels reveal an (almost) exponential suppression of higher-order energy corrections indicating rapid convergence of the expansion. The rate of convergence is mostly independent of the target state or the nucleus, except for the ground and first 2^{+} states in ^{6}Li that both converge slightly slower.

**Figure 6**. High-order binding energies from NCSM-PT employing ${N}_{\text{max}}^{\text{ref}}=0$ reference states using the Hamiltonian described in section 3. All calculations are performed including configurations up to *N*_{max} = 8 and employ an oscillator frequency of ℏω = 20 MeV. The SRG parameter is set to α = 0.08 fm^{4}. Reproduced from Tichai et al. [36] under the Creative Commons CC-BY license.

The high-order benchmarks strongly motivate the use of low-order partial sums as good approximations to the binding energies of heavier systems. Subsequently, systems with mass number *A* ≈ 30 are investigated through realistic MCPT calculations in large model spaces.

#### 6.4.2. Low Orders in *A*≈30 Nuclei

In Figure 7, ground-state energies of carbon, oxygen, and fluorine isotopes are shown and compared to large-scale importance truncated NCSM (IT-NCSM) diagonalization whenever available. The reference states are obtained from a diagonalization in a *N*_{max} = 0 or 2 space. In all cases a HF single-particle basis is employed in order to minimize the dependence on the oscillator frequency^{9}.

**Figure 7**. Reference (/) and second-order NCSM-PT (/) energies with ${N}_{\text{max}}^{\text{ref}}=0$ and 2, respectively, for the ground states of ^{11−20}C, ^{16−26}O, and ^{17−31}F using the Hamiltonian described in section 3. All calculations are performed using 13 oscillator shells and an oscillator frequency of ℏω = 20 MeV. The SRG parameter is set to α = 0.08 fm^{4}. Importance-truncated NCSM calculations () are shown for comparison. Experimental values are indicated by black bars. Reproduced from Tichai et al. [36] under the Creative Commons CC-BY license.

For all nuclei the reference energies and second-order NCSM-PT results show a sizeable dependence on the size of the reference space. Throughout all investigated isotopic chains though, results from reference states built within a ${N}_{\text{max}}^{\text{ref}}=2$ space almost perfectly reproduce the exact (IT-)NCSM ones. This significant improvement in ground-state energies hints at important correlations incorporated through the reference states obtained from a ${N}_{\text{max}}^{\text{ref}}=2$ diagonalization that are absent for ${N}_{\text{max}}^{\text{ref}}=0$, thus providing an ideal compromise between computational efficiency and accuracy. In particular, neutron-rich fluorine isotopes are out of reach of conventional NCSM calculations and NCSM-PT provides an efficient *ab initio* approach to investigate the neutron drip line. A single NCSM-PT calculation requires typically two to three orders of magnitude less computational resources than the corresponding IT-NCSM calculation.

In practice, the reference states employed in the above calculations contain between several hundreds of thousands up to a few million determinants, thus, providing an excellent account of static correlations. Note that in most cases the reference state accounts for up to 80% of the overall binding energy such that residual dynamical correlations can indeed be grasped efficiently from low-order perturbation theory.

#### 6.4.3. Low-Lying Spectroscopy

As already exemplified in the high-order investigation, excited states can be straightforwardly accessed through NCSM-PT by targeting different reference states from the NCSM spectrum. From absolute NCSM-PT binding energies, excitation energies are obtained by subtracting the correlated ground-state energy. In Figure 8, the associated NCSM-PT spectra are compared to bare NCSM calculations for a selection of open-shell carbon and oxygen isotopes. All calculations employ a HO single-particle basis to separate center-of-mass degrees of freedom in the many-body wave function. It is well-known that NCSM excitation energies of states with identical parity display a much faster convergence than absolute binding energies, thus yielding stable results in the right-hand columns of each panel. When a level re-ordering appears with increasing *N*_{max} in the NCSM calculation, e.g. for the two lowest states in ^{12}C or the third and fourth states in ^{19}O, NCSM-PT reproduces the correct level ordering at small values of ${N}_{\text{max}}^{\text{ref}}$. As for ground-state energies, NCSM-PT results based on a reference state with ${N}_{\text{max}}^{\text{ref}}=2$ reproduce well the NCSM spectra, while going to ${N}_{\text{max}}^{\text{ref}}=4$ only refines the quality of a subset of excitation energies.

**Figure 8**. Spectra obtained via second-order NCSM-PT for selected carbon and oxygen isotopes using the Hamiltonian described in section 3. All calculations are performed using 13 oscillator shells and an oscillator frequency of ℏω = 20 MeV. The SRG parameter is set to α = 0.08 fm^{4}. A HO single-particle basis was used to separate center-of-mass contaminations. Large-scale IT-NCSM calculations are shown for comparison. Reproduced from Tichai et al. [36] under the Creative Commons CC-BY license.

The NCSM-PT framework is thus highly valuable to perform light-weighted perturbative calculations in medium-light systems up to mass numbers *A* ≲ 40. Due to its versatility and conceptual simplicity the low-lying spectrum of genuine open-shell nuclei can be described at low computational cost.

## 7. Bogoliubov Many-Body Perturbation Theory

While the use of multi-configurational reference states can efficiently resolve situations of strong static correlations, it displays several limitations, i.e., (i) the physical origin of the underlying correlations in the reference state is unclear, (ii) it does not easily ensure size-extensivity, and (iii) it is numerically prohibitive in heavy nuclei. The objective is thus to present an alternative based on single-reference product states that bypasses these limitations.

### 7.1. Rationale

An alternative route to lift the particle-hole degeneracy of the reference state in open-shell systems is to authorize the reference state to break a symmetry of the underlying Hamiltonian. In semi-magic nuclei, the relevant symmetry is *U*(1) global gauge symmetry associated with particle-number conservation^{10}. Breaking *U*(1) symmetry permits to efficiently deal with Cooper pair instability associated with the superfluid character of open-shell nuclei. The degeneracy of a Slater determinant with respect to particle-hole excitations is lifted via the use of a Bogoliubov reference state and transferred into a degeneracy with respect to transformations of the symmetry group. As a consequence, the ill-defined (i.e., singular) expansion of exact quantities with respect to a symmetry-conserving Slater determinant is replaced by a well-behaved one. Extending the treatment to doubly open-shell nuclei requires a similar treatment of *SU*(2) symmetry associated with the conservation of angular momentum.

Eventually, the degeneracy with respect to *U*(1) transformations must also be lifted by restoring the symmetry. However, BMBPT only restores the symmetry in the limit of an all-order resummation, and thus displays a symmetry contamination at any finite order. While BMBPT can still be used as a stand-alone approach as is done in the present work, it eventually provides the first step toward the implementation of the so-called particle-number projected BMBPT (PNP-BMBPT) [58] that restores good particle number at any truncation order.

### 7.2. Bogoliubov Algebra

The BMBPT formalism is based on the introduction of the Bogoliubov reference state

where ${C}$ is a complex normalization constant and |0〉 denotes the physical vacuum. The Bogoliubov product state presently defining the *P*-space of perturbation theory is a vacuum for the quasi-particle operators $\left\{{\beta}_{k}^{\u2020},{\beta}_{k}\right\}$, i.e.,

obtained from the creation and annihilation operators associated with a basis of the one-body Hilbert space via the unitary Bogoliubov transformation [103]

Generically speaking, the Bogoliubov transformation is only constrained by unitarity such that a large manifold of Bogoliubov states is at hand. To actually set up the perturbation theory, a particular Bogoliubov reference state must be specified. Typically, the transformation matrices (*U, V*) are obtained by solving the Hartree-Fock-Bogoliubov (HFB) variational problem that naturally extends the simpler HF approximation to treat pairing correlations while sticking to a single reference product state. The columns of the transformation matrices (*U, V*) correspond to the eigenvectors of the HFB eigenvalue equation [103] whereas the associated eigenvalues {*E*_{k} > 0} deliver the so-called quasi-particle energies^{11}. Since Bogoliubov states are not eigenstates of the particle-number operator *A*, the expectation value of *A* is constrained to match a specific number of particles, e.g., the particle number *A*_{0} of the target system. In the HFB method for example, the constraint is enforced via the use of a Lagrange multiplier λ in the minimization of the expectation value of the grand potential

In actual applications, separate Lagrange multipliers λ_{N} and λ_{Z} are used to constrain proton and neutron numbers *N* and *Z*, respectively. In the subsequent formalism *A* stands for either one of them.

### 7.3. Quasi-Particle Normal Ordering

In the next step, Wick's theorem is employed to normal order the grand potential Ω with respect to the Bogoliubov reference state

where Ω^{ij} denotes the normal-ordered component involving *i* (*j*) quasi-particle creation (annihilation) operators, e.g.,

The tensors defining each normal-ordered term display antisymmetry properties, i.e.,

where σ(*P*) refers to the signature of the permutation *P*. The notation *P*(…|…) denotes a separation into the *i* quasiparticle creation operators and the *j* quasiparticle annihilation operators such that permutations are only considered among members of the same group. Thus, Ω^{00} is the expectation value of Ω in |Φ〉 whereas Ω^{[2]} and Ω^{[4]} define effective, i.e., normal-ordered, one-body and two-body operators, respectively. Working in the particle-number-conserving normal-ordered two-body approximation (PNO2B) [83], the effective three-body part Ω^{[6]} is presently discarded^{12}. Details on the normal-ordering procedure as well as expressions of the matrix elements of each operator Ω^{ij} in terms of the original matrix elements of the Hamiltonian and of the (*U, V*) matrices can be found in Signoracci et al. [53] and Ripoche et al. [83].

### 7.4. Partitioning

To set up the perturbation theory, the Hamiltonian (i.e., grand potential) must be partitioned into an one-body unperturbed part Ω_{0} and a residual part Ω_{1}, i.e.,

Focusing on the case where the Bogoliubov reference state is the solution of the HFB variational problem, i.e., using a Møller-Plesset scheme, Ω appearing in Equation (61) is naturally partitioned given that

and that Ω^{11} is in diagonal form, i.e.,

with *E*_{k} > 0 for all *k*. Introducing all many-body states obtained via an even number of quasi-particle excitations of the vacuum

the unperturbed system is fully characterized by its complete set of orthonormal eigenstates in Fock space

where the strict positivity of unperturbed excitations

characterizes the lifting of the particle-hole degeneracy authorized by the spontaneous breaking of *U*(1) symmetry in open-shell nuclei at the mean-field level.

With these ingredients at hand, the perturbation theory can be entirely worked out algebraically or diagrammatically. This can be done on the basis of a (imaginary) time-dependent formalism or of a time-independent formalism. While the former framework leads to working with Feynman (time-dependent) diagrams, the latter makes use of Goldstone (time-ordered) diagrams. Recently, the complete Rayleigh-Schrödinger BMBPT formalism, including the automatic generation and algebraic evaluation of all possible diagrams appearing at an arbitrary order *n* on the basis of 2N and full 3N interactions has been published in Arthuis et al. [95].

Eventually, the BMBPT expansion of the correlation energy can be written in compact form as a *superfluid extension* of the Goldstone formula (Equation 19b)

where the Hamiltonian is replaced by the grand potential.

### 7.5. Low Orders

As a result of Wick's theorem with respect to |Φ〉, the first few orders contribute to Equation (70), with Ω^{(p)} ≡ *E*^{(p)} − λ*A*^{(p)} and ${\Omega}_{\text{ref}}={\Omega}^{(0)}+{\Omega}^{(1)}$, according to^{13}

The lifting of the degeneracy with respect to particle-hole excitations is embodied in the fact that the energy denominators in Equation (71) are non-singular and well behaved. Indeed, quasi-particle energies are bound from below by the superfluid *pairing gap* at the Fermi energy, i.e.,

This would not be true in standard MBPT based on a Slater determinant reference state, where energy denominators associated with particle-hole excitations within the open shell would be zero in Equation (45). Of course, BMBPT does strictly reduce to standard MBPT in a closed-shell system [95]. In particular, the single third-order diagram whose algebraic expression is given in Equation (71c) generates the three, i.e., particle-particle, hole-hole and particle-hole, third-order HF-MBPT diagrams [95]. This reduction of the number of diagrams at any order *p* is a consequence of working in a quasi-particle representation that does not distinguish particle and hole states. Conversely, all summations over quasi-particle labels run over the entire dimension of the one-body Hilbert space, which significantly increases the computational cost compared to standard MBPT. In any case, low-order BMBPT corrections only induce low polynomial scaling with respect to quasi-particle summation and do not suffer from the storage of large tensors as in more sophisticated all-order many-body approaches such as (B)CC or IMSRG.

Extracting the *p*-order contribution to the binding energy from Equation (70) requires the subtraction of the Lagrange term λ*A*^{(p)}. Computing *A*^{(p)} can be done straightforwardly by replacing the leftmost operator Ω by *A* in Equation (70) [95]. As the reference state is constrained to have the correct particle number on average, it implies that ${A}_{0}^{(1)}={\text{A}}_{0}$. Working with the HFB reference state, it can be shown that ${A}_{0}^{(2)}=0$ due to the fact that Ω^{20} = Ω^{02} = 0. Consequently, the first correction to the average particle number appears at third order such that

i.e., the computed average particle number does not match the targeted number *A*_{0} of the physical system. This feature requires an iterative BMBPT scheme in order for the particle number to be correct at the working order, e.g., *p* ≥ 3, of interest [104]. To do so, one needs to rerun the HFB calculation with a *p*-dependent chemical potential such that, through a series of iterations, one eventually obtains, e.g., ${A}_{0}^{(1)}+\dots +{A}_{0}^{(p)}={\text{A}}_{0}$. Such a costly algorithm can fortunately be very well approximated by an *a posteriori* correction scheme that entirely bypasses the iterative scheme [104]. The third-order results presented in section 7.6.1 have been computed without any adjustment of the average particle number whereas the novel ones discussed in section 10.2 have been obtained on the basis of the *a posteriori* correction scheme.

While the BMBPT expansion efficiently grasps static correlation effects associated to nuclear superfluidity, the breaking of a continuous symmetry in a finite quantum system is always fictitious. Consequently, a full-fledged many-body formalism requires the additional restoration of the broken symmetry, i.e., a mixing of gauge-rotated Bogoliubov vacua which are connected to each other via (highly non-perturbative) symmetry transformations. While the formalism has already been laid out [58], realistic calculations remain yet to be performed. Still, proof-of-principle applications to the model pairing Hamiltonian employing particle-number-projected BCC theory [60] revealed the significant impact of the symmetry projection in the weakly broken regime corresponding to open-shell nuclei in the vicinity of shell closures.

### 7.6. Results

#### 7.6.1. Low-Order Calculations in Mid-Mass Nuclei

Figure 9 displays ground-state energies, two-neutron separation energies, particle-number variance and perturbative particle-number corrections along oxygen, calcium, and nickel isotopic chains at the mean-field (i.e., HFB) level, as well as at second- and third-order in the BMBPT expansion. Top panels reveal that, while static correlations have partially been accounted for by employing a symmetry-broken reference state, the bulk of dynamical correlations are efficiently grasped via low-order BMBPT corrections. For (sub-)closed-shell nuclei the third-order correction is consistently suppressed indicating rapid convergence. In open-shell systems third-order partial sums are strongly contaminated due to a significant excess of neutrons brought by the third-order contribution to the average neutron number. As discussed earlier, calculations at third-order and beyond must eventually be done while constraining the average particle number to match the physical value [104]. This is reported on in section 10.2.

**Figure 9**. BMBPT systematics along O, Ca, and Ni isotopic chains: **(A)** absolute binding energies, **(B)** two-neutron separation energies, **(C)** neutron-number dispersion, and **(D)** perturbative correction to the average neutron number. The employed Hamiltonian is defined in section 3. Plot markers correspond to HFB (), second-order BMBPT (), and third-order BMBPT without particle-number adjustment (). All calculations are performed using 13 oscillator shells and an oscillator frequency of ℏω = 20 MeV. The SRG parameter is set to α = 0.08 fm^{4}. Reproduced from Tichai et al. [37] under the Creative Commons CC-BY license.

Figure 9B exhibits a qualitative reproduction of two-neutron separation energies already at the mean-field level. Results are quantitatively improved once second-order effects are incorporated. Figure 9C shows the neutron-number dispersion $\sigma \equiv \sqrt{\langle {A}^{2}\rangle -{\langle A\rangle}^{2}}$ that grows with mass number. While the second-order contribution does not decrease yet the neutron-number dispersion, one expects higher orders to do so [104]. In closed-shell systems, the particle-number dispersion is zero as a hallmark of the particle-number-conserving character of the wave function throughout the expansion.

A detailed study of the convergence characteristics of the BMBPT expansion via the calculation of high-order corrections similar to the ones presented in section 5 (section 6.4.1) for HF-MBPT (NCSM-PT) calculations of closed-shell (open-shell) nuclei has just be completed [104] and is thus not reported on here.

#### 7.6.2. Comparison to Non-perturbative Calculations

Figure 10 compares second-order BMBPT ground-state and two-neutron separation energies with results obtained from state-of-the-art non-perturbative many-body frameworks along oxygen, calcium, and nickel isotopic chains. In particular, IT-NCSM provides essentially exact reference results along the oxygen chain. In heavier closed-shell nuclei, advanced CR-CC(2,3) calculations also provide reference results.

**Figure 10**. Absolute ground-state binding energies **(top)** and two-neutron separation energies **(bottom)** along O, Ca, and Ni isotopic chains using the Hamiltonian defined in section 3. Plot markers correspond to second-order BMBPT (), second-order NCSM-PT (), large-scale IT-NCSM (), GSCGF-ADC(2) (), MR-IMSRG(2) (), and CR-CC(2,3) (). All calculations are performed using 13 oscillator shells and an oscillator frequency of ℏω = 20 MeV. The SRG parameter is set to α = 0.08 fm^{4}. Reproduced from Tichai et al. [37] under the Creative Commons CC-BY license.

In all cases, second-order BMBPT is in excellent agreement with other methods, displaying deviations of <2%. In particular, IT-NCSM results are very well reproduced in oxygen isotopes. While NCSM-PT (see section 6) and MR-IMSRG systematically generate stronger binding, GSCGF at the ADC(2) level is very close for all systems investigated. In the case of closed-shell nuclei the stronger binding obtained in the CR-CC(2,3) calculation highlights the importance of an (approximate) incorporation of 3-particle/3-hole excitations.

As seen from two-neutron separation energies, all *ab initio* methods consistently predict the same shell structure and location of the neutron trip line. This feature highlights both the great success of recently developed many-body methods and the excellent performance of perturbative techniques such as NCSM-PT and BMBPT in particular. On the other hand, the strong deviation of absolute binding energies^{14} from experimental data, most pronounced in neutron-rich calcium and nickel isotopes, is common to all employed frameworks and clearly points to defects of the employed Hamiltonian. A detailed discussion of this crucial issue is postponed to section 10, where a new family of chiral Hamiltonians is tested with the goal to cure the poor agreement presently obtained for nuclear ground-state observables of mid-mass nuclei.

## 8. Importance Truncation

In previous sections, MBPT was considered as a full-fledged standalone framework to access the solution of the quantum many-body problem. However, MBPT techniques can also be used to support non-perturbative approaches by pre-processing the many-body configuration space. Subsequently, the concept of *importance truncation* (IT) is presented as a procedure to a priori remove a set of tensor entries to be solved for on the basis of an importance measure [98, 105–107]. Historically, this idea was first applied in electronic structure theory to perform a pre-selection of multi-reference CI amplitudes [108, 109].

### 8.1. Original Context: Configuration Interaction Methods

The IT concept is ideally suited for basis expansion methods, such as general configuration interaction (CI) approaches or, more specifically, the NCSM. Given the expansion of an eigenstate |Ψ_{n}〉 of *H* in some many-body basis |ϕ_{i}〉, the importance of each basis state can be quantified through the associated expansion coefficient ${C}_{i}^{(n)}\equiv \langle {\varphi}_{i}|{\text{\Psi}}_{n}\rangle $. This importance measure can be estimated within first-order MCPT, discussed in section 6, starting from a multi-configurational reference state |Ψ_{ref}〉

The reference state provides an initial approximation of |Ψ_{n}〉 within the reference space ${{M}}_{\text{ref}}$ and κ_{i} is used to assess *a priori* the importance of each basis state $\text{|}{\varphi}_{i}\rangle \notin {{M}}_{\text{ref}}$. Only states with importance measures above an importance threshold κ_{min} are included into the model space for the subsequent CI calculation.

Since the eigenvalue problem is solved exactly in the importance-truncated model space, IT provides a variational approximation to the full solution. Moreover, the effect of discarded configurations can be estimated via the second-order MCPT energy correction (56). This together with extrapolations to vanishing importance thresholds κ_{min} → 0 can be used to reconstruct the energies in the full model space. As testified by the results already presented in Figures 7, 8, 10, the IT concept has allowed one to push NCSM calculations beyond p-shell nuclei, up to neutron-rich oxygen isotopes. More details on importance-truncated CI approaches are discussed in Roth [98] and Stumpf et al. [106].

### 8.2. New Context: Coupled-Cluster Method

In the context of non-perturbative expansion methods, the general idea is to discard irrelevant entries of the largest mode-*n* tensor that needs to be solved for and that drives the numerical scaling as well as the storage cost of the method. This is done by *a priori* estimating the importance of each of its entries on the basis of a less costly many-body method, e.g., MBPT. While the underlying formalism is only sketched in this section, the reader is referred to Tichai et al. [107] for a detailed discussion.

#### 8.2.1. Wave-Function Ansatz

To illustrate the concept, BMBPT is employed as a method to pre-process the cluster amplitudes constituting the unknown tensors to be solved for within the BCC method. This non-perturbative approach is based on the wave-function ansatz [53]

where the connected quasiparticle cluster operator ${T}\equiv {{T}}_{1}+{{T}}_{2}+{{T}}_{3}+\dots $ is defined through

etc. The BCC amplitudes ${t}_{{k}_{1}\cdots {k}_{2n}}^{2n0}$ constitute the unknown mode-*n* tensors of present interest.

Instead of solving for all unknown cluster amplitudes, IT selects a subset of tensor entries that are determined non-perturbatively by solving BCC amplitude equations while residual corrections from the discarded tensor entries are treated in low-order BMBPT. Consequently, the ground-state energy ${{E}}_{0}$ consists of a contribution from BCC in the IT subspace ${{E}}_{0}^{(\text{IT})}$ and a residual BMBPT correction ${\delta}_{0}^{(\text{res})}$ from the complementary subspace

Of course, the selection of the IT subspace must be performed *without the knowledge* of the full solution, which is where low-order BMBPT estimates enter into play.

#### 8.2.2. Importance-Truncated Tensor

Considering the *n*-tuple BCC amplitude tensor ${{T}}_{n}\equiv \left\{{t}_{{k}_{1}\dots {k}_{2n}}^{2n0}\right\}$, the corresponding *importance-truncated* tensor based on the IT measure ${\kappa}^{(p)}({t}_{{k}_{1}\dots {k}_{2n}}^{2n0})$ from BMBPT at order *p* is obtained as

where κ_{min} defines the IT threshold. The original tensor is recovered in the limit of κ_{min} → 0, i.e.,

With an importance-truncated tensor comes the associated *data compression ratio*

which relates the initial amount of data to the compressed one resulting from the IT process. Whenever *R*_{c} > 1, the compressed tensor requires less storage than the original one.

#### 8.2.3. Results

Applying IT to ${{T}}_{3}$, the goal is eventually (not done here) to solve BCCSDT equations for the retained entries and correct for the omitted ones in perturbation. The feasibility of the approach directly depends on the reduction offered by the IT for a desired accuracy given that a full BCCSDT calculation is currently undoable in realistic model spaces, even when employing an angular-momentum-coupled scheme.

Figure 11 displays the impact of the data compression obtained by applying IT for ^{18}O as a function of the HO single-particle basis size. To provide a reference, the compression of ${{T}}_{2}$ is illustrated first in the left panel before coming to the real challenge constituted by the treatment of ${{T}}_{3}$. Three different storage schemes are displayed corresponding to a full storage of all tensor entries (diamonds), the tensor entries allowed by fundamental symmetries of the interaction (squares) and the IT-compressed entries from perturbative estimate (circle) for different values of κ_{min}. The percentage associated with the IT-data indicates the corresponding error on the second-order contribution to the ground-state energy

Results demonstrate that even for ${{T}}_{2}$, i.e., a mode-4 tensor, the storage of all index combinations requires tremendous amounts of memory (> 200 GB per working copy) that is at the edge of supercomputing facilities^{15}. Employing a symmetry-adapted storage scheme, several orders of magnitude can be saved. More importantly, discarding irrelevant tensor entries through IT additionally compresses the data by three orders of magnitude at the price of inducing a systematic error of <1%, i.e., a few hundreds of keV, on the second-order correlation energy. Correspondingly, the storage of the ${{T}}_{2}$ amplitudes is lowered to less than 1 MB per working copy. One can assign an “effective one-body dimension” to this number of tensor entries yielding ${e}_{\text{max}}^{\text{(eff)}}\le 5$, i.e., converged BCCSD results could be obtained using an effective model space including less than six major shells.

**Figure 11**. Number of tensor entries of ${{T}}_{2}$ amplitudes **(left)** and ${{T}}_{3}$ amplitudes **(right)** as a function of model space using the Hamiltonian described in section 3. The oscillator frequency is given by ℏω = 20 MeV. Different plot markers indicate different storage scheme including all entries (), symmetry-allowed entries (), and IT-compressed entries (). For *e*_{max} ≥ 10 an additional *E*_{3max}-truncation () for the ${{T}}_{3}$ amplitudes is introduced. All calculations employ a canonical HFB reference state. Reproduced from Tichai et al. [107] with permission from European Physical Journal (EPJ).

While the storage of ${{T}}_{2}$ amplitudes is actually within reach of state-of-the-art capacities, the extension to ${{T}}_{3}$ amplitudes poses a severe computational problem. Performing the same analysis as above, even the symmetry-adapted tensor requires more than 100TB of memory. Performing the IT-compression of ${{T}}_{3}$, the corresponding fourth-order contribution to the ground-state energy

can be evaluated with 1% accuracy while discarding 99.99% of all ${{T}}_{3}$ tensor entries. Because triples correction to the ground-state energy first appear at fourth order, a 1% relative error on $\Delta {\Omega}_{0}^{\left[{4}_{T}\right]}$ relates to an absolute error of a few keV, which is negligible in mid-mass *ab initio* studies.

Eventually, the use of IT shifts the boundary of what is computationally unfeasible such that, in the following years, the implementation of IT-BCCSDT defines an ambitious and yet reachable goal.

## 9. Basis Optimization

The present section discusses how MBPT can be used as an auxiliary tool to construct an alternative single-particle basis accelerating the convergence of non-perturbative many-body methods as compared to commonly employed bases.

### 9.1. Rationale

In nuclear structure calculations, tensors defining *H* are typically expressed in the eigenbasis of the spherical HO Hamiltonian. Next, a self-consistent mean-field, e.g., HF, calculation is typically performed to generate the reference state such that the beyond-mean-field step and the associated tensors are transformed into the corresponding mean-field, e.g., HF, one-body basis.

While the use of HO orbitals authorizes the rigorous center-of-mass factorization in NCSM calculations, they are plagued at long distances with a pathological Gaussian falloff and a strong dependence of the many-body results on the width of the confining HO potential. Contrarily, HF orbitals display a proper exponential falloff but may induce a sizeable center-of-mass contamination and a slower convergence in large-scale NCSM diagonalizations compared to HO orbitals.

### 9.2. Natural Orbitals

Naturally, one wonders about the existence of a single-particle basis that admit only little sensitivity to the oscillator frequency while maintaining rapid model space convergence. Because the exact one-body density matrix

contains information about the fully correlated wave function, its eigenbasis, defining so-called *natural orbitals*, is expected to deliver an optimal choice.

This expectation was confirmed by employing the one-body density matrix obtained from a large-scale NCSM calculation [110]. For mid-mass nuclei, this strategy has the downside to require a full CI solution to obtain the optimized single-particle basis. However, as already observed several decades ago in quantum chemistry calculations [111, 112], *approximate* natural orbitals perform surprisingly well in applications. Instead of using the exact CI wave function, a correlated density matrix is constructed within a polynomially scaling expansion method, MBPT providing a particularly simple example to do so [113]. Similar approaches can be followed by using the Λ-extension of CC theory or by diagonalizing the dressed one-body propagator in SCGF theory. Note, however, that the construction of an auxiliary basis already involves a computationally non-trivial solution within a non-perturbative many-body approach.

### 9.3. Density Matrix in Closed-Shell MBPT

Starting from Equation (25b) and following Strayer et al. [114], the one-body density matrix up to second order in the residual interaction can be written as [113]

where

denotes the zeroth-order HF density matrix whereas

denote its MBPT corrections whose matrix elements are expressed in the HF basis. Explicit expressions for the individual corrections to the mean-field density in terms of two-body matrix elements can be found in Strayer et al. [114] and Tichai et al. [113].

### 9.4. Results

The impact of using natural orbitals is now gauged by performing a systematic variation of the underlying single-particle basis in NCSM calculations. Figure 12 displays ^{16}O ground-state energy as a function of the frequency characterizing the underling HO basis for various values of *N*_{max}. Results obtained from the NCSM diagonalization performed in the HO and HF bases as well as in the natural-orbital basis extracted from second-order HF-MBPT are displayed.

**Figure 12**. NCSM ground-state energy of ^{16}O as a function of the frequency ℏω (denoted as ℏΩ on the figure) of the underling HO basis for *N*_{max} = 2, 4, 6, 8, 10 spaces. The NCSM diagonalization itself is performed either in the initial HO basis **(left)**, in the HF basis **(center)**, or in the natural-orbital basis from second-order HF-MBPT **(right)**. The Hamiltonian described in section 3 is employed. Reproduced from Tichai et al. [113] with permission from the American Physical Society.

The left panel shows the typical strong dependence of NCSM results on the underlying frequency when using HO one-body basis states, thus yielding a pronounced parabolic shape around the variational minimum for all values of *N*_{max}. While the use of HF orbitals (middle panel) lowers the sensitivity to ℏω the model space convergence is not improved and the ground-state energy has a linear falloff at higher frequencies. Finally, the NCSM ground-state energy obtained using the natural-orbital basis (right panel) is almost independent of the oscillator frequency and displays a faster model-space convergence as a function of *N*_{max}. Similar conclusions hold for other observables, e.g., charge radii, low-lying excitation energies and electro-magnetic transition strengths, as discussed in Tichai et al. [113].

Eventually, the use of one-body density matrices from low-order MBPT to construct of an optimal computational basis is a simple yet powerful tool to account for high-lying particle-hole excitations that are otherwise hard to grasp in configuration-driven methods like NCSM.

## 10. Testing Novel _{χ}EFT Nuclear Hamiltonians

The computational efficiency of perturbative approaches makes them ideally suited for survey calculations over a large range of nuclei. A specific scenario for such survey calculations is the development of novel chiral EFT interactions, where large numbers of many-body calculations are necessary to constrain or validate the choice of low-energy constants as well as to characterize the effect of different regulator formulations. Recently, a new family of consistent chiral 2N plus 3N interactions up to N3LO based on non-local regulators was developed [115]. This new set of interactions is able to simultaneously reproduce ground-state energies and charge radii up into the medium-mass regime. Moreover, the full sequence of chiral orders from LO to N3LO is available for different cutoff values, which facilitates a rigorous uncertainty quantification. In the applications discussed below, the new interactions are consistently SRG-evolved in the 2N and 3N sectors down to α = 0.04 fm^{4} and tested through non-expensive perturbative methods over a larger set of nuclei than those addressed in Hüther et al. [115].

### 10.1. NCSM-PT

First, NCSM-PT is employed to calculate ground-state energies along the oxygen chain using interactions from Hüther et al. [115]. In Figure 13, results obtained from N2LO and N3LO interactions are compared to experiment including a quantification of systematic theoretical uncertainties resulting from both the many-body and interaction truncations. The benchmarks discussed in section 6 showed that calculations with the ${N}_{max}^{\text{ref}}=2$ typically provide a good agreement with non-perturbative methods and that ground-state energies for larger reference spaces typically fall between ${N}_{max}^{\text{ref}}=0$ and 2 results. Therefore, the difference between these two calculations offers a convenient way to estimate many-body uncertainties in NCSM-PT calculations. Uncertainties resulting from truncating the chiral expansion of nuclear interactions are estimated from the order-by-order variation of the energies using the prescription discussed in Epelbaum et al. [116], Binder et al.[117, 118]. The interaction-induced uncertainty of an observable *X*_{N3LO} at order N3LO, e.g., is given by $max(Q|{X}_{\text{N3LO}}-{X}_{\text{N2LO}}|,{Q}^{2}|{X}_{\text{N2LO}}-{X}_{\text{NLO}}|,{Q}^{3}|{X}_{\text{NLO}}-{X}_{\text{LO}}|,{Q}^{5}|{X}_{\text{LO}}|)$ with the expansion parameter *Q* estimated as the ratio of a typical momentum scale characterizing mid-mass systems over the breakdown scale, which results in *Q* ≈ 1/3 [118]. The uncertainty band shown in Figure 13 combines both uncertainties, where many-body uncertainties are typically larger than interaction uncertainties.

**Figure 13**. Even oxygen isotopes ground-state energies from reference (empty symbols) and second-order (full symbols) NCSM-PT with ${N}_{max}^{\text{ref}}=2$ **(left)** and from IM-NCSM **(right)**. The consistent non-local 2N plus 3N chiral interactions at N^{2}LO () and N^{3}LO () presented in Hüther et al. [115] are employed with a cutoff Λ = 500 MeV. The uncertainty bands represent both many-body and interaction uncertainties (see text). Both many-body calculations make use a natural-orbital basis obtained in 13 oscillator shells with ℏω = 20 MeV.

The NCSM-PT ground-state energies can be compared with results from the in-medium no-core shell model (IM-NCSM) displayed on the right-hand side of Figure 13. The non-perturbative IM-NCSM method combines a multi-reference in-medium SRG transformation of the Hamiltonian that decouples a small NCSM model space, with a subsequent NCSM calculation using the transformed Hamiltonian [52]. The associated uncertainty bands include once again both many-body and interaction uncertainties. The comparison reveals an excellent agreement between NCSM-PT and IM-NCSM, which themselves agree with experiment within estimated uncertainties.

### 10.2. BMBPT

In addition to low-order NCSM-PT benchmarks in oxygen isotopes, the novel set of interactions is tested on a large set of mid-mass semi-magic nuclei up to third-order in BMBPT^{16}.

In Figure 14, ground-state energies of oxygen and calcium isotopes are displayed for two different values of *E*_{3max} = 14, 16. The first striking result is that the overbinding obtained in Ca isotopes with the ‘standard’ Hamiltonian (see Figure 10) is resolved, i.e., the systematic trend throughout O and Ca isotopes is now consistent with experimental data. In particular, low-order BMBPT results provide a much better reproduction of experimental ground-state energies of very neutron-rich Ca isotopes than with the “standard” Hamiltonian for which the degrading increases systematically with neutron number.

**Figure 14**. BMBPT systematics along O **(left)** and Ca **(right)** isotopic chains using the non-local chiral two- plus three-nucleon interactions at N^{3}LO presented in Hüther et al. [115]. Top panel: absolute binding energies. Bottom panel: two-neutron separation energies. Plot markers correspond to HFB (), second-order BMBPT (), and third-order BMBPT (). All calculations are performed using 15 oscillator shells and an oscillator frequency of ℏω = 20 MeV. The single-particle orbital angular-momentum quantum number was limited to *l* ≤ 10. Open and closed symbols correspond to different truncations *E*_{3max} = 14, 16, respectively.

Going in more details, one observes that the reference HFB energy, i.e., the lowest order in the BMBPT expansion, accounts for much less binding (i.e., less than half of the total binding energy) than with the “standard” Hamiltonian, thus indicating a “harder” interaction^{17}. Moreover, the third-order correction is typically larger in the present case. While the first three BMBPT orders are systematically and strongly suppressed, this indicates that the new interactions are indeed less perturbative. One can thus expect higher orders to contribute non-negligibly.

Focusing now on two-neutron separation energies, the end results are very satisfactory and of similar or even greater quality than with the “standard” Hamiltonian. Interestingly, HFB results are further away in the present case, indicating that not only absolute values but also the trend with neutron number is different. Still, second and third order contributions consistently compensate for this apparent, but in fact fictitious, degrading at the HFB level. One eventually observes that results are still less accurate near major shell closures, which is consistent with the expectation that restoring particle-number symmetry through PBMBPT [58] will have the largest impact near shell closures.

While the results in oxygen isotopes are virtually identical for both *E*_{3max} values, the truncation in *E*_{3max} plays an increasingly important role in neutron-rich calcium isotopes, i.e., ground-state energies beyond *A* = 50 obtained for *E*_{3max} = 14 and 16 start deviating more strongly, pointing to the importance of truncated three-body matrix elements as already shown in Hüther et al. [115]. This constitutes a technical challenge for the future when *ab initio* calculations move to even heavier and more neutron-rich nuclei.

Let us now further discuss these novel results by adding both many-body and interaction uncertainties to BMBPT(3) calculations. Interaction uncertainties are estimated as discussed previously. To estimate the many-body uncertainty arising from the finite truncation of the BMBPT expansion, a nucleus-independent error of 3% is assumed. A recent study of high-order BMBPT calculations revealed that the third-order partial sum typically induces a many-body error of about 2% [104]. Because the employed Hamiltonian was softer and further SRG-evolved to lower resolution scale, a more conservative 3%-error in the current setting seems plausible. The overall uncertainty is then given by summing this many-body error and the interaction uncertainty. Additional uncertainties due to the finite single-particle basis are presently not incorporated. Clearly, the construction of reasonable models of many-body uncertainties is highly nontrivial and employing simplistic parametric models, e.g., geometric sums for perturbation series, can lead to wrong estimates if the perturbation series does not obey the underlying assumptions.

In Figure 15, ground-state energies of oxygen and calcium isotopes are displayed for both N2LO and N3LO Hamiltonians. In each case, results are accompanied by a band combining the two sources of uncertainty according to the method outlined above. While the mean value of the results only slightly changes when going from N2LO to N3LO, the uncertainty is consistently reduced. Furthermore, N3LO results are essentially consistent with experimental data within the error band, although the narrowing of that uncertainty allows one to stipulate that this agreement is deteriorating when going to neutron-rich isotopes.

**Figure 15**. BMBPT ground-state energies for oxygen and calcium isotopes at third order of the perturbative expansion using consistent non-local chiral two- plus three-nucleon interactions at N^{2}LO (blue) and N^{3}LO (red). Open symbols show the corresponding HFB ground-state energies. The uncertainty bands represent both, many-body and interaction uncertainties (see text). All calculations are performed using 15 oscillator shells and an oscillator frequency of ℏω = 20 MeV. The single-particle orbital angular-momentum quantum number was limited to *l* ≤ 10.

Clearly, the above features points toward a systematically improved quality of the novel set of interactions compared to prior generations of chiral Hamiltonians that tend to fail in mid-mass systems. In the future, this statement will be further benchmarked by validating the consistency of low-lying spectroscopy and charge radii with experimental data in the mid-mass region.

## 11. Conclusions and Outlook

The present paper reviewed the status of MBPT techniques in the field of *ab initio* many-body calculations of finite nuclei. After discussing formal properties of the power-series ansatz, the main goal was to provide an in-depth understanding of the MBPT expansion and its range of applicability. Most importantly, two ways of extending Slater-determinant-based MBPT toward open-shell systems were discussed. Multi-configurational perturbation theory and Bogoliubov MBPT were shown to provide complementary ways to tackle the particle-hole degeneracy and the associated strong correlations in open-shell nuclei.

Each of the subparts dedicated to the various flavors of many-body perturbation theory displayed results of large-scale calculations touching upon the very frontier of current *ab initio* studies dedicated to mid-mass (open-shell) nuclei. While highly accurate non-perturbative alternatives are available, the formal and computational simplicity is at the very heart of perturbative approaches. Indeed, these features enable large surveys along complete isotopic chains at a small fraction of the computational cost necessary to utilize more sophisticated approaches. In view of the increasing activity related to the construction of improved generations of nuclear Hamiltonians, simple and yet accurate many-body schemes constitute a very valuable tool to quickly characterize their overall performance.

Furthermore, MBPT provides not only a standalone many-body framework, but can also be used as an inexpensive pre-processing tool to accelerate non-perturbative methods. Two facets discussed in this work are data compression from MBPT-based importance truncation and basis optimization from low-order MBPT density matrices. Both tools have been used with great success to either tame the curse of dimensionality or cure pathological defects in the computational basis at low computational price.

This (by far not exhaustive) list of MBPT-related applications in nuclear theory highlights the start of a *renaissance of perturbation theory* fostered by renormalization group techniques used to soften nuclear Hamiltonians. In parallel, the deepened understanding of infrared-divergences, their origin in open-shell systems and systematic ways to overcome them will hopefully help cure the (unfortunately still present) disfavor against MBPT for nuclear structure research.

Future developments will expand many of the ideas laid out in this work. Better control on uncertainties in the nuclear Hamiltonian pave the way for targeting more exotic systems in the upcoming years. In particular, a proper account of static correlations associated with nuclear deformation has been identified as a critical goal by various groups worldwide. Consequently, the implementation of frameworks employing a systematic breaking (and restoration) of *SU*(2) symmetry will be of high relevance. However, due to the increasing number of basis states these efforts will require significant computational resources as well as extensive formal developments. To control the increase of computational requirements (memory and runtime), so-called *tensor factorization techniques* have been proposed recently where high-mode tensors are decomposed into sums of product of lower-rank ones [107, 119]. While initial proof-of-principle applications have demonstrated the high potential, extensive additional research is required in this direction.

On a short timescale the benchmarking of importance truncation and natural orbitals in non-perturbative medium-mass applications will deepen the understanding of the role of the computational basis and distribution of correlations in Hilbert space. Furthermore, by employing data compression tools will greatly help relax the many-body truncation at play such that first-principle calculations will witness a level of accuracy that has been out of reach so far.

## Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

## Author Contributions

All authors contributed to the writing of the text. RR performed *ab initio* calculations of NCSM-PT and IM-NCSM. AT performed calculations of BMBPT.

## Funding

This publication was based on work supported in part by the framework of the Espace de Structure et de réactions Nucléaires Théorique (ESNT) at CEA, the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—Projektnummer 279384907—SFB 1245, and the BMBF through contract 05P18RDFN1 and Verbundprojekt 05P2018 (ErUM-FSP T07).

## Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Acknowledgments

We would like to thank M. Frosini for generating the data used in Figure 2. Parts of this work were performed at the Lichtenberg cluster at the computing center of the TU Darmstadt and the JURECA supercomputing facility at Forschungszentrum Jülich.

## Footnotes

1. ^The unperturbed Hamiltonian *H*_{0} sums individual kinetic energies in this case.

2. ^Both states are supposed to be adiabatically connected when the perturbation *H*_{1} is switched on.

3. ^The projective character of standard MBPT or CC is to be distinguished from an *expectation value* approach, where the correlated state appears both as the bra and the ket in the evaluation of the energy.

4. ^A quantum-mechanical method is coined as *size-extensive* if the energy of a systems computed with this method scales linearly in the number of particles.

5. ^The first disconnected contribution originally appears at fourth order. It can be shown that the renormalization terms cancel such disconnected contributions at every order [69, 85, 86], thus providing the final connected form.

6. ^The HF problem is solved in a symmetry-restricted way, enforcing rotational invariance of the resulting single-particle basis.

7. ^Using the notion of tensors, the summation over a common index of several tensors defines a *tensor contraction*.

8. ^If the density matrix involved in Equation (50) were the one of the fully correlated eigenstate of *H*, the one-body operator *h* would be nothing but the *Baranger Hamiltonian* [99, 100]. Whenever the reference state reduces to the HF Slater determinant, *h* identifies with the HF one-body Hamiltonian.

9. ^The HF problem is solved employing a so-called *equal-filling approximation* to ensure rotational invariance of the mean-field density. The underlying NCSM-PT formulation is based on *m*-scheme quantities and does not employ angular-momentum coupling techniques as most correlation expansions do. Consequently, even- and odd-mass systems can be described on equal footing.

10. ^In fact the relevant symmetry group is the direct product *U*_{N}(1) × *U*_{Z}(1) since both proton and neutron number are conserved separately.

11. ^In the present work, the Bogoliubov transformation is limited to treat the like-particle pairing although it could be further generalized to address neutron-proton pairing as well.

12. ^The particle-number-conserving nature of the PNO2B approximation further requires to drop specific contributions to Ω^{[4]} as well; see Ripoche et al. [83] for a detailed discussion.

13. ^Non-canonical contributions are set to zero here given that we use the HFB reference state. When using an arbitrary Bogoliubov reference state, additional non-canonical contributions arise [95].

14. ^The strong deviation also concerns charge radii that are not presently discussed.

15. ^Solving non-linear problems requires the storage of various copies of the coupled-cluster amplitudes.

16. ^As compared to Figure 9, third-order are presently corrected for the particle-number contamination via the so-called *a posteriori* correction studied and validated in Demol et al. [104].

17. ^This feature reflects both the different characters of the chiral Hamiltonians themselves and the fact that the “standard” Hamiltonian was SRG-evolved down to α = 0.08 fm^{4} instead of α = 0.04 fm^{4} for the new ones.

## References

1. Weinberg S. Nuclear forces from chiral lagrangians. *Phys Lett B*. (1990) **251**:288. doi: 10.1016/0370-2693(90)90938-3

2. Weinberg S. Effective chiral lagrangians for nucleon-pion interactions and nuclear forces. *Nucl Phys Sect B*. (1991) **363**:3–18. doi: 10.1016/0550-3213(91)90231-L

3. Weinberg S. Three-body interactions among nucleons and pions. *Phys Lett B*. (1992) **295**:114–21. doi: 10.1016/0370-2693(92)90099-P

4. Entem DR, Machleidt R. Accurate charge-dependent nucleon-nucleon potential at fourth order of chiral perturbation theory. *Phys Rev C*. (2003) **68**:041001. doi: 10.1103/PhysRevC.68.041001

5. Epelbaum E. Few-nucleon forces and systems in chiral effective field theory. *Prog Part Nucl Phys*. (2006) **57**:654–741. doi: 10.1016/j.ppnp.2005.09.002

6. Machleidt R, Entem DR. Chiral effective field theory and nuclear forces. *Phys Rep*. (2011) **503**:1–75. doi: 10.1016/j.physrep.2011.02.001

7. Wiringa RB, Stoks VGJ, Schiavilla R. Accurate nucleon-nucleon potential with charge-independence breaking. *Phys Rev C*. (1995) **51**:38–51. doi: 10.1103/PhysRevC.51.38

8. Machleidt R. The high-precision, charge-dependent Bonn nucleon-nucleon potential (CD-Bonn). *Phys Rev C*. (2000) **63**:024001. doi: 10.1103/PhysRevC.63.024001

9. Otsuka T, Suzuki T, Holt JD, Schwenk A, Akaishi Y. Three-body forces and the limit of oxygen isotopes. *Phys Rev Lett*. (2010) **105**:32501. doi: 10.1103/PhysRevLett.105.032501

10. Hergert H, Bogner SK, Binder S, Calci A, Langhammer J, Roth R, et al. In-medium similarity renormalization group with chiral two- plus three-nucleon interactions. *Phys Rev C*. (2013) **87**:34307. doi: 10.1103/PhysRevC.87.034307

11. Navrátil P, Quaglioni S, Stetcu I, Barrett B. Recent developments in no-core shell-model calculations. *J Phys G Nucl Part Phys*. (2009) **36**:83101. doi: 10.1088/0954-3899/36/8/083101

12. Roth R, Langhammer J, Calci A, Binder S, Navrátil P. Similarity-transformed chiral NN+3N interactions for the *ab initio* description of 12-C and 16-O. *Phys Rev Lett*. (2011) **107**:72501. doi: 10.1103/PhysRevLett.107.072501

13. Barrett BR, Navrátil P, Vary JP. *Ab initio* no core shell model. *Prog Part Nucl Phys*. (2013) **69**:131–81. doi: 10.1016/j.ppnp.2012.10.003

14. Gezerlis A, Tews I, Epelbaum E, Gandolfi S, Hebeler K, Nogga A, et al. Quantum Monte Carlo calculations with chiral effective field theory interactions. *Phys Rev Lett*. (2013) **111**:032501. doi: 10.1103/PhysRevLett.111.032501

15. Carlson J, Gandolfi S, Pederiva F, Pieper SC, Schiavilla R, Schmidt KE, et al. Quantum Monte Carlo methods for nuclear physics. *Rev Mod Phys*. (2015) **87**:1067. doi: 10.1103/RevModPhys.87.1067

16. Lynn JE, Tews I, Carlson J, Gandolfi S, Gezerlis A, Schmidt KE, et al. Quantum Monte Carlo calculations of light nuclei with local chiral two- and three-nucleon interactions. *Phys Rev C*. (2017) **96**:054007. doi: 10.1103/PhysRevC.96.054007

17. Lynn JE, Tews I, Gandolfi S, Lovato A. Quantum Monte Carlo methods in nuclear physics: recent advances. *Annu Rev Nucl Part Sci*. (2019) **69**:279–305. doi: 10.1146/annurev-nucl-101918-023600

18. Nogga A, Kamada H, Glöckle W. Solution of the Faddeev-Yakubovsky equations using realistic NN and 3N interactions. *Nuclear Phys A*. (2001) **689**:357–60. doi: 10.1016/S0375-9474(01)00854-5

19. Lazauskas R, Carbonell J. The Faddeev-Yakubovsky symphony. *Few Body Syst*. (2019) **60**:62. doi: 10.1007/s00601-019-1529-5

20. Dickhoff WH, Barbieri C. Selfconsistent Green's function method for nuclei and nuclear matter. *Prog Part Nucl Phys*. (2004) **52**:377–496. doi: 10.1016/j.ppnp.2004.02.038

21. Cipollone A, Barbieri C, Navrátil P. Isotopic chains around oxygen from evolved chiral two- and three-nucleon interactions. *Phys Rev Lett*. (2013) **111**:062501. doi: 10.1103/PhysRevLett.111.062501

22. Carbone A, Cipollone A, Barbieri C, Rios A, Polls A. Self-consistent Green's functions formalism with three-body interactions. *Phys Rev C*. (2013) **88**:54326. doi: 10.1103/PhysRevC.88.054326

23. Somá V, Cipollone A, Barbieri C, Navrátil P, Duguet T. Leading chiral three-nucleon forces along isotope chains in the calcium region. *Phys Rev C*. (2014) **89**:61301. doi: 10.1103/PhysRevC.89.061301

24. Raimondi F, Barbieri C. Nuclear electromagnetic dipole response with the Self-Consistent Green's Function formalism. *Phys Rev C*. (2019) **99**:054327. doi: 10.1103/PhysRevC.99.054327

25. Somá V, Navrátil P, Raimondi F, Barbieri C, Duguet T. Novel chiral Hamiltonian and observables in light and medium-mass nuclei. *Phys. Rev. C* (2019) **101**:014318. doi: 10.1103/PhysRevC.101.014318

26. Kowalski K, Dean DJ, Hjorth-Jensen M, Papenbrock T, Piecuch P. Coupled cluster calculations of ground and excited states of nuclei. *Phys Rev Lett*. (2004) **92**:132501. doi: 10.1103/PhysRevLett.92.132501

27. Hagen G, Papenbrock T, Dean DJ, Hjorth-Jensen M. *Ab initio* coupled-cluster approach to nuclear structure with modern nucleon-nucleon interactions. *Phys Rev C*. (2010) **82**:34330. doi: 10.1103/PhysRevC.82.034330

28. Binder S, Langhammer J, Calci A, Roth R. *Ab initio* path to heavy nuclei. *Phys Lett B*. (2014) **736**:119–23. doi: 10.1016/j.physletb.2014.07.010

29. Henderson TM, Dukelsky J, Scuseria GE, Signoracci A, Duguet T. Quasiparticle coupled cluster theory for pairing interactions. *Phys Rev*. (2014) **C89**:54305. doi: 10.1103/PhysRevC.89.054305

30. Jansen GR, Engel J, Hagen G, Navratil P, Signoracci A. *Ab initio* coupled-cluster effective interactions for the shell model: application to neutron-rich oxygen and carbon isotopes. *Phys Rev Lett*. (2014) **113**:142502. doi: 10.1103/PhysRevLett.113.142502

31. Hagen G, Papenbrock T, Hjorth-Jensen M, Dean DJ. Coupled-cluster computations of atomic nuclei. *Rep Prog Phys*. (2014) **77**:096302. doi: 10.1088/0034-4885/77/9/096302

32. Roth R, Langhammer J. Padé-resummed high-order perturbation theory for nuclear structure calculations. *Phys Lett B*. (2010) **683**:272–7. doi: 10.1016/j.physletb.2009.12.046

33. Langhammer J, Roth R, Stumpf C. Spectra of open-shell nuclei with Padé-resummed degenerate perturbation theory. *Phys Rev C*. (2012) **86**:054315. doi: 10.1103/PhysRevC.86.054315

34. Tichai A, Langhammer J, Binder S, Roth R. Hartree-Fock many-body perturbation theory for nuclear ground-states. *Phys Lett B*. (2016) **756**:283–8. doi: 10.1016/j.physletb.2016.03.029

35. Hu BS, Xu FR, Sun ZH, Vary JP, Li T. *Ab initio* nuclear many-body perturbation calculations in the Hartree-Fock basis. *Phys Rev C*. (2016) **94**:14303. doi: 10.1103/PhysRevC.94.014303

36. Tichai A, Gebrerufael E, Vobig K, Roth R. Open-shell nuclei from no-core shell model with perturbative improvement. *Phys Lett B*. (2018) **786**:448–52. doi: 10.1016/j.physletb.2018.10.029

37. Tichai A, Arthuis P, Duguet T, Hergert H, Somá V, Roth R. Bogoliubov many-body perturbation theory for open-shell nuclei. *Phys Lett B*. (2018) **786**:195. doi: 10.1016/j.physletb.2018.09.044

38. Tsukiyama K, Bogner SK, Schwenk A. In-medium similarity renormalization group for nuclei. *Phys Rev Lett*. (2011) **106**:222502. doi: 10.1103/PhysRevLett.106.222502

39. Tsukiyama K, Bogner SK, Schwenk A. In-medium similarity renormalization group for open-shell nuclei. *Phys Rev C*. (2012) **85**:061304. doi: 10.1103/PhysRevC.85.061304

40. Bogner SK, Hergert H, Holt JD, Schwenk A, Binder S, Calci A, et al. Nonperturbative shell-model interactions from the in-medium similarity renormalization group. *Phys Rev Lett*. (2014) **113**:142501. doi: 10.1103/PhysRevLett.113.142501

41. Hergert H, Bogner SK, Morris TD, Schwenk A, Tsukiyama K. The in-medium similarity renormalization group: a novel *Ab initio* method for nuclei. *Phys Rep*. (2016) **621**:165. doi: 10.1016/j.physrep.2015.12.007

42. Parzuchowski NM, Morris TD, Bogner SK. *Ab initio* excited states from the in-medium similarity renormalization group. *Phys Rev C*. (2017) **95**:044304. doi: 10.1103/PhysRevC.95.044304

43. Morris TD, Simonis J, Stroberg SR, Stumpf C, Hagen G, Holt JD, et al. Structure of the lightest tin isotopes. *Phys Rev Lett*. (2018) **120**:152503. doi: 10.1103/PhysRevLett.120.152503

44. Stroberg SR, Calci A, Hergert H, Holt JD, Bogner SK, Roth R, et al. Nucleus-dependent valence-space approach to nuclear structure. *Phys Rev Lett*. (2017) **118**:032502. doi: 10.1103/PhysRevLett.118.032502

45. Miyagi T, Abe T, Okamoto R, Otsuka T. Ground-state energies and charge radii of 4He, 16O, 40Ca, and 56Ni in the unitary-model-operator approach. *Prog Theor Exp Phys*. (2015) **2015**:041D01. doi: 10.1093/ptep/ptv042

46. Miyagi T, Abe T, Okamoto R, Otsuka T. Introduction of the one-body correlation operator in the unitary-model-operator approach. *Phys Rev C*. (2017) **96**:054312. doi: 10.1103/PhysRevC.96.054312

47. Piecuch P, Gour JR, Włoch M. Left-eigenstate completely renormalized equation-of-motion coupled-cluster methods: review of key concepts, extension to excited states of open-shell systems, and comparison with electron-attached and ionized approaches. *Int J Quant Chem*. (2009) **109**:3268. doi: 10.1002/qua.22367

48. Jansen GR. Spherical coupled-cluster theory for open-shell nuclei. *Phys Rev C*. (2013) **88**:024305. doi: 10.1103/PhysRevC.88.024305

49. Holt JD, Menéndez J, Simonis J, Schwenk A. Three-nucleon forces and spectroscopy of neutron-rich calcium isotopes. *Phys Rev C*. (2014) **90**:024312. doi: 10.1103/PhysRevC.90.024312

50. Somá V, Duguet T, Barbieri C. *Ab-initio* self-consistent Gorkov-Green's function calculations of semi-magic nuclei. I. Formalism at second order with a two-nucleon interaction. *Phys Rev C*. (2011) **84**:64317. doi: 10.1103/PhysRevC.84.064317

51. Hergert H, Bogner SK, Morris TD, Binder S, Calci A, Langhammer J, et al. *Ab initio* multireference in-medium similarity renormalization group calculations of even calcium and nickel isotopes. *Phys Rev C*. (2014) **90**:041302. doi: 10.1103/PhysRevC.90.041302

52. Gebrerufael E, Vobig K, Hergert H, Roth R. *Ab initio* description of open-shell nuclei: merging no-core shell model and in-medium similarity renormalization group. *Phys Rev Lett*. (2017) **118**:152503. doi: 10.1103/PhysRevLett.118.152503

53. Signoracci A, Duguet T, Hagen G, Jansen GR. *Ab initio* Bogoliubov coupled cluster theory for open-shell nuclei. *Phys Rev C*. (2015) **91**:64320. doi: 10.1103/PhysRevC.91.064320

54. Duguet T. Symmetry broken and restored coupled-cluster theory: I. Rotational symmetry and angular momentum. *J Phys G Nucl Part Phys*. (2015) **42**:25107. doi: 10.1088/0954-3899/42/2/025107

55. Ui H, Takeda G. Equivalence of stochastic quantization method to conventional field theories through super transformation invariance. *Prog Theor Phys*. (1983) **70**:176. doi: 10.1143/PTP.70.176

56. Yannouleas C, Landman U. Symmetry breaking and quantum correlations in finite systems: studies of quantum dots and ultracold Bose gases and related nuclear and chemical methods. *Rep Prog Phys*. (2007) **70**:2067. doi: 10.1088/0034-4885/70/12/R02

57. Papenbrock T, Weidenmueller HA. Effective field theory for finite systems with spontaneously broken symmetry. *Phys Rev C*. (2014) **89**:014334. doi: 10.1103/PhysRevC.89.014334

58. Duguet T, Signoracci A. Symmetry broken and restored coupled-cluster theory. II. Global gauge symmetry and particle number. *J Phys*. (2017) **G44**:15103. doi: 10.1088/0954-3899/44/1/015103

59. Qiu Y, Henderson TM, Zhao J, Scuseria GE. Projected coupled cluster theory. *J Chem Phys*. (2017) **147**:064111. doi: 10.1063/1.4991020

60. Qiu Y, Henderson TM, Duguet T, Scuseria GE. Particle-number projected Bogoliubov coupled cluster theory. Application to the pairing Hamiltonian. *Phys Rev C*. (2019) **99**:044301. doi: 10.1103/PhysRevC.99.044301

62. Frame D, He R, Ipsen I, Lee D, Lee D, Rrapaj E. Eigenvector continuation with subspace learning. *Phys Rev Lett*. (2018) **121**:032501. doi: 10.1103/PhysRevLett.121.032501

63. Frame DK. *Ab initio* simulations of light nuclear systems using eigenvector continuation and auxiliary field Monte Carlo. *arXiv [Preprint]. 2019. arXiv:1905.02782.*

64. Demol P, Duguet T, Ekström A, Frosini M, Hebeler K, König S, et al. Improved many-body expansions from eigenvector continuation. *arXiv [Preprint]. 2019. arXiv:1911.12578*. doi: 10.1103/PhysRevC.101.041302

65. Hebeler K, Bogner SK, Furnstahl RJ, Nogga A, Schwenk A. Improved nuclear matter calculations from chiral low-momentum interactions. *Phys Rev C*. (2011) **83**:031301. doi: 10.1103/PhysRevC.83.031301

66. Drischler C, Carbone A, Hebeler K, Schwenk A. Neutron matter from chiral two- and three-nucleon calculations up to N3LO. *Phys Rev C*. (2016) **94**:054307. doi: 10.1103/PhysRevC.94.054307

67. Wellenhofer C, Drischler C, Schwenk A. Dilute Fermi gas at fourth order in effective field theory. *Phys Lett B.* (2020) **802**:135247.

68. Drischler C, Hebeler K, Schwenk A. Chiral interactions up to next-to-next-to-next-to-leading order and nuclear saturation. *Phys Rev Lett*. (2019) **122**:042501. doi: 10.1103/PhysRevLett.122.042501

69. Shavitt I, Bartlett RJ. *Many-Body Methods in Chemistry and Physics: MBPT and Coupled-Cluster Theory (Cambridge Molecular Science)*. Cambridge University Press (2009). Available online at: http://www.amazon.com/Many-Body-Methods-Chemistry-Physics-Coupled-Cluster/dp/052181832X. doi: 10.1017/CBO9780511596834

70. Frame D, He R, Ipsen I, Lee D, Lee D, Rrapaj E. Eigenvector continuation with subspace learning. *Phys Rev Lett*. (2018) **121**:032501. doi: 10.1103/PhysRevLett.121.032501

71. König S, Ekström A, Hebeler K, Lee D, Schwenk A. Eigenvector continuation as an efficient and accurate emulator for uncertainty quantification. *arXiv [Preprint]. 2019. arXiv: 1909.08446*.

72. Ekström A, Hagen G. Global sensitivity analysis of bulk properties of an atomic nucleus. *Phys Rev Lett*. (2019) **123**:252501. doi: 10.1103/PhysRevLett.123.252501

73. Demol P, Duguet T, Ekström A, Frosini M, Hebeler K, König S, et al. Improved many-body expansions from eigenvector continuation. *Phys Rev C*. (2019) **101**:041302(R). doi: 10.1103/PhysRevC.101.041302

74. Ekström A, Baardsen G, Forssén C, Hagen G, Hjorth-Jensen M, Jansen GR, et al. Optimized chiral nucleon-nucleon interaction at next-to-next-to-leading order. *Phys Rev Lett*. (2013) **110**:192502. doi: 10.1103/PhysRevLett.110.192502

75. Ekström A, Jansen GR, Wendt KA, Hagen G, Papenbrock T, Carlsson BD, et al. Accurate nuclear radii and binding energies from a chiral interaction. *Phys Rev C*. (2015) **91**:051301. doi: 10.1103/PhysRevC.91.051301

76. Jurgenson E, Navratil P, Furnstahl R. Evolution of nuclear many-body forces with the similarity renormalization group. *Phys Rev Lett*. (2009) **103**:2. doi: 10.1103/PhysRevLett.103.082501

77. Entem DR, Machleidt R. Accurate charge-dependent nucleon-nucleon potential at fourth order of chiral perturbation theory. *Phys Rev C*. (2003) **68**:041001(R). doi: 10.1103/PhysRevC.68.041001

78. Navrátil P. Local three-nucleon interaction from chiral effective field theory. *Few-Body Syst*. (2007) **41**:117–40. doi: 10.1007/s00601-007-0193-3

79. Roth R, Binder S, Vobig K, Calci A, Langhammer J, Navrátil P. Medium-mass nuclei with normal-ordered chiral NN+3N interactions. *Phys Rev Lett*. (2012) **109**:052501. doi: 10.1103/PhysRevLett.109.052501

80. Roth R, Calci A, Langhammer J, Binder S. Evolved chiral *NN*+3*N* Hamiltonians for ab initio nuclear structure calculations. *Phys Rev C*. (2014) **90**:024325. doi: 10.1103/PhysRevC.90.024325.

81. Maris P, Vary JP, Calci A, Langhammer J, Binder S, Roth R. ^{1}2C properties with evolved chiral three-nucleon interactions. *Phys Rev C*. (2014) **90**:014314. doi: 10.1103/PhysRevC.90.014314

82. Gebrerufael E, Calci A, Roth R. Open-shell nuclei and excited states from multireference normal-ordered Hamiltonians. *Phys Rev C*. (2016) **93**:031301. doi: 10.1103/PhysRevC.93.031301

83. Ripoche J, Tichai A, Duguet T. Normal-ordered k-body approximation in particle-number-breaking theories. *arXiv [Preprint]. 2019. arXiv:1908.00765*. doi: 10.1140/epja/s10050-020-00045-8

84. Hugenholtz NM. Perturbation theory of large quantum systems. *Physica*. (1957) **23**:481–532. doi: 10.1016/S0031-8914(57)92950-6

85. Goldstone J. Derivation of the Brueckner many-body theory. *Proc R Soc Lond*. (1957) **239**:267. doi: 10.1098/rspa.1957.0037

86. Bloch C. Sur la détermination de l'état fondamental d'un systéme de particules. *Nucl Phys*. (1958) **7**:451. doi: 10.1016/0029-5582(58)90284-0

87. Wick GC. The evaluation of the collision matrix. *Phys Rev*. (1950) **80**:268. doi: 10.1103/PhysRev.80.268

88. Paldus J, Wong HC. Computer generation of Feynman diagrams for perturbation theory I. General algorithm. *Comput Phys Commun*. (1973) **6**:1–7. doi: 10.1016/0010-4655(73)90016-7

89. Wong HC, Paldus J. Computer generation of Feynman diagrams for perturbation theory II. Program description. *Comput Phys Commun*. (1973) **6**:9–16. doi: 10.1016/0010-4655(73)90017-9

90. Csépes Z, Pipek J. An effective recursive algorithm for generating many-body Hugenholtz and Goldstone diagrams. *J Comput Phys*. (1988) **77**:1–17. doi: 10.1016/0021-9991(88)90153-2

91. Lyons J, Moncrieff D, Wilson S. Diagrammatic many body perturbation expansion for atoms and molecules: automatic generation & analysis of 5th order Hugenholtz energy diagrams. *Comput Phys Commun*. (1994) **84**:91–101. doi: 10.1016/0010-4655(94)90205-4

92. Kállay M, Surján PR. Higher excitations in coupled-cluster theory. *J Chem Phys*. (2001) **115**:2945. doi: 10.1063/1.1383290

93. Kállay M, Szalay PG, Surján PR. A general state-selective multireference coupled-cluster algorithm. *J Chem Phys*. (2002) **117**:980. doi: 10.1063/1.1483856

94. Stevenson PD. Automatic generation of vacuum amplitude many-body perturbation series. *Int J Mod Phys C*. (2003) **14**:1135. doi: 10.1142/S0129183103005236

95. Arthuis P, Duguet T, Tichai A, Lasseri RD, Ebran JP. ADG: Automated generation and evaluation of many-body diagrams I. Bogoliubov many-body perturbation theory. *Comput Phys Commun*. (2019) **240**:202–27. doi: 10.1016/j.cpc.2018.11.023

96. Rolik Z, Szabados Á, Surján PR. On the perturbation of multiconfiguration wave functions. *J Chem Phys*. (2003) **119**:1922. doi: 10.1063/1.1584424

97. Surján PR, Rolik Z, Szabados Á, Köhalmi D. Partitioning in multiconfiguration perturbation theory. *Annalen der Physik (Leipzig)*. (2004) **13**:223–31. doi: 10.1002/andp.200310074

98. Roth R. Importance truncation for large-scale configuration interaction approaches. *Phys Rev C*. (2009) **79**:064324. doi: 10.1103/PhysRevC.79.064324

99. Baranger M. A definition of the single-nucleon potential. *Nucl Phys A*. (1970) **149**:225. doi: 10.1016/0375-9474(70)90692-5

100. Duguet T, Hergert H, Holt JD, Somá V. Nonobservable nature of the nuclear shell structure: meaning, illustrations, and consequences. *Phys Rev C*. (2015) **92**:034313. doi: 10.1103/PhysRevC.92.034313

101. Hose G, Kaldor U. Diagrammatic many-body perturbation theory for general model spaces. *J Phys B Mol Phys*. (1979) **12**:3827. doi: 10.1088/0022-3700/12/23/012

102. Hose G, Kaldor U. A general-model-space diagrammatic perturbation theory. *Phys Scripta*. (1980) **21**:357. doi: 10.1088/0031-8949/21/3-4/019

103. Ring P, Schuck P. *The Nuclear Many-Body Problem*. New York, NY: Springer Verlag (1980). doi: 10.1007/978-3-642-61852-9

104. Demol P, Frosini M, Tichai A, Somá V, Duguet T. Bogoliubov many-body perturbation theory under constraint as a series. Unpublished (2020).

105. Roth R, Navrátil P. *Ab initio* study of ^{4}0Ca with an importance-truncated no-core shell model. *Phys Rev Lett*. (2007) **99**:092501. doi: 10.1103/PhysRevLett.99.092501

106. Stumpf C, Braun J, Roth R. Importance-truncated large-scale shell model. *Phys Rev C*. (2016) **93**:021301. doi: 10.1103/PhysRevC.93.021301

107. Tichai A, Ripoche J, Duguet T. Pre-processing the nuclear many-body problem. *Eur Phys J A*. (2019) **55**:90. doi: 10.1140/epja/i2019-12758-6

108. Buenker RJ, Peyerimhoff SD. Individualized configuration selection in CI calculations with subsequent energy extrapolation. *Theor Chim Acta*. (1974) **35**:33–58. doi: 10.1007/BF02394557

109. Buenker RJ, Peyerimhoff SD. Energy extrapolation in CI calculations. *Theor Chim Acta*. (1975) **39**:217–28. doi: 10.1007/BF00555301

110. Constantinou C, Caprio MA, Vary JP, Maris P. Natural orbital description of the halo nucleus 6He. *Nucl Sci Techn*. (2017) **28**:179. doi: 10.1007/s41365-017-0332-6

111. Hay PJ. On the calculation of natural orbitals by perturbation theory. *J Chem Phys*. (1973) **59**:2468. doi: 10.1063/1.1680359

112. Siu AKQ, Hayes EF. Configuration interaction procedure based on the calculation of perturbation theory natural orbitals: applications to H2 and LiH. *J Chem Phys*. (1974) **61**:37–40. doi: 10.1063/1.1681646

113. Tichai A, Müller J, Vobig K, Roth R. Natural orbitals for *Ab initio* no-core shell model calculations. *Phys Rev C*. (2019) **99**:034321. doi: 10.1103/PhysRevC.99.034321

114. Strayer MR, Bassichis WH, Kerman AK. Correlation effects in nuclear densities. *Phys Rev C*. (1973) **8**:1269. doi: 10.1103/PhysRevC.8.1269

115. Hüther T, Vobig K, Hebeler K, Machleidt R, Roth R. Family of chiral two- plus three-nucleon interactions for accurate nuclear structure studies. *arXiv [Preprint]. 2019. arXiv:1911.04955*.

116. Epelbaum E, Krebs H, Meißner UG. Improved chiral nucleon-nucleon potential up to next-to-next-to-next-to-leading order. *Eur Phys J A*. (2015) **51**:53. doi: 10.1140/epja/i2015-15053-8

117. Binder S, Calci A, Epelbaum E, Furnstahl RJ, Golak J, Hebeler K, et al. Few-nucleon systems with state-of-the-art chiral nucleon-nucleon forces. *Phys Rev C*. (2016) **93**:044002. doi: 10.1103/PhysRevC.93.044002

118. Binder S, Calci A, Epelbaum E, Furnstahl RJ, Golak J, Hebeler K, et al. Few-nucleon and many-nucleon systems with semilocal coordinate-space regularized chiral nucleon-nucleon forces. *Phys Rev C*. (2018) **98**:014002. doi: 10.1103/PhysRevC.98.014002

Keywords: many-body theory, *ab initio*, perturbation theory, correlation expansion, open-shell nuclei

Citation: Tichai A, Roth R and Duguet T (2020) Many-Body Perturbation Theories for Finite Nuclei. *Front. Phys.* 8:164. doi: 10.3389/fphy.2020.00164

Received: 28 January 2020; Accepted: 20 April 2020;

Published: 05 June 2020.

Edited by:

Carlo Barbieri, University of Surrey, United KingdomCopyright © 2020 Tichai, Roth and Duguet. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Alexander Tichai, alexander.tichai@physik.tu-darmstadt.de