Many-body perturbation theories for finite nuclei

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 flavours of many-body perturbation theory.


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 | A k denotes the k-th eigenstate of the A-body system with eigenvalue E A k . 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][2][3][4][5][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][12][13], or quantum Monte Carlo (QMC) [14][15][16][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 reimport 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][21][22][23][24][25] and coupled cluster (CC) [26][27][28][29][30][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][33][34][35][36][37], inmedium similarity renormalization group (IMSRG) [10,[38][39][40][41][42][43][44] or the unitary model operator approach (UMOA) [45,46]. All of these methods provide a consistent description of groundstate 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 Slaterdeterminant 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 multiconfigurational 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 mediumheavy isotopic chains [37]. Next, the additional or alternative breaking of SU(2) rotational symmetry associated with angularmomentum 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][56][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][66][67][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 postprocessing tool to further improve upon non-perturbative manybody approaches.
While the present document discusses in detail manybody 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][71][72][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.

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 multideterminantal character as well as the resolution scale of the employed Hamiltonian.

Rationale
It has often been argued in the past that the nuclear manybody 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.

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 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).
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 upperleft 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.
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 preprocessing 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 (closedsub-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.
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 nucleonnucleon 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.

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 lowenergy 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][2][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 {|p ≡ c † p |0 } of the one-body Hilbert space H 1 , the form Frontiers in Physics | www.frontiersin.org  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 .
The Hamiltonian is, thus, represented via a set of one-, two-, and three-body matrix elements t pq ,v pqrs andw 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.

Similarity Renormalization Group
While the tensors defining the Hamiltonian built within χEFT may display large low-to-high momentum couplings, preprocessing 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 α ∈ R, 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 prediagonalization 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 tensorsv pqrs (α),w pqrstu (α) . . . whose UV elements linking singleparticle 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 modelspace convergence while keeping the effect of induced manybody 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.

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-leadingorder (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 SRGevolved 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 socalled normal-ordered two-body approximations (NO2B) are employed, depending on the particular nature of the A-body reference state [79,82,83].

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.

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 {| k , E (0) k ; k ∈ 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., 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 |H 1 |χ 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.

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.

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][85][86] 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.
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 , the first nontrivial 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 .

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 with the initial condition C (0) m0 = δ 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 orderp corrections, this approach is limited by the storage of the Hamiltonian due to the extensive size of the manybody 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.

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-2k 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 secondquantized 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 errorsafe way. Consequently, the last tool introduced to tackle this difficulty consists of an automatized generation and evaluation of diagrams [88][89][90][91][92][93][94][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.

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 {c † p } 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.

Normal Ordering
Applying Wick's theorem with respect to | , the Hamiltonian can be rewritten in terms of normal-ordered contributions pq : c † p c q : + 1 4 pqrs H [4] pqrs : c † p c † q c s c r : + . . . , (35) where : : denotes the normal order of the involved creation and annihilation operators. Thus, H [0] is the expectation value of H in | whereas H [2] pq and H [4] pqrs 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]).

Partitioning
To explicitly set up the partitioning of the Hamiltonian (Equation 8), one adds and subtracts a diagonal normal-ordered one-body operatorH [2] ≡ p e p : c † p c p : such that withH [2] ≡ H [2] −H [2] = p =q 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 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 ǫ ab··· ij··· > 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 takingH [2] = 0, i.e., H 1 = H [4] . This is obtained by using the reference Slater determinant | solution of the Hartree-Fock (HF) variational problem and by definingH [2] from the eigenvalues of the onebody HF Hamiltonian.

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.

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 noncanonical 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 secondorder 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.

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].

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.
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.

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   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 higherorder corrections (Figures 4E,F).

Realistic Calculations
The previous results reveal that using an optimized HF reference state combined with a sufficiently soft interaction defines a wellcontrolled 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.
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 7 Using the notion of tensors, the summation over a common index of several tensors defines a tensor contraction.
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 thirdorder 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.

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].

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 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 multideterminantal 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 ref are denoted by | ν , i.e., as a capital Greek letter carrying a lower-case Greek index, whereas (ii) Slater determinants outside M 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 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.

Partitioning
Formally, the unperturbed Hamiltonian is written in the spectral representation as where | k denote the NCSM eigenvectors within M 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 {| i / ∈ M 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 (0) i of the unperturbed Slater determinants making up the Q space, i.e., 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

Low Orders
With the partitioning of the Hamiltonian defined above, zerothand first-order MCPT contributions to the energy read as 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 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 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 | µ ′ 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 dim(M ref ) 2 · n 2 p · n h , where n p and n h denote the number of particle and hole states, respectively.

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.

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.
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.

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 .
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 ref max = 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 ref max = 2 diagonalization that are absent for N ref max = 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.

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 singleparticle 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  The NCSM-PT framework is thus highly valuable to perform light-weighted perturbative calculations in mediumlight 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.

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.

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 openshell 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 openshell 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.

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 {β † k , β k }, 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] β k ≡ p U * pk c p + V * pk c † p , (59a) 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 quasiparticle 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.

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., 11 In the present work, the Bogoliubov transformation is limited to treat the likeparticle pairing although it could be further generalized to address neutron-proton pairing as well.
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].

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 20 = 02 = 0 (65) and that 11 is in diagonal form, i.e., 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.

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 ref = (0) + (1) , according to 13 (2) = − 1 24 (3) = + 1 8 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 quasiparticle 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 (1) 0 = A 0 . Working with the HFB reference state, it can be shown that A (2) 0 = 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., 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. Figure 9 displays ground-state energies, two-neutron separation energies, particle-number variance and perturbative particlenumber corrections along oxygen, calcium, and nickel isotopic chains at the mean-field (i.e., HFB) level, as well as at secondand 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 thirdorder 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 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 σ ≡ A 2 − A 2 that grows with mass number. While the second-order contribution does not decrease yet the neutronnumber 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.

Low-Order Calculations in Mid-Mass Nuclei
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. 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.

Comparison to Non-perturbative Calculations
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 14 The strong deviation also concerns charge radii that are not presently discussed. 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.

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][106][107]. Historically, this idea was first applied in electronic structure theory to perform a pre-selection of multi-reference CI amplitudes [108,109].

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 (n) i ≡ φ i | n . 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 ref and κ i is used to assess a priori the importance of each basis state |φ i / ∈ M 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 importancetruncated CI approaches are discussed in Roth [98] and Stumpf et al. [106].

New Context: Coupled-Cluster Method
In the context of non-perturbative expansion methods, the general idea is to discard irrelevant entries of the largest moden 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.

Wave-Function Ansatz
To illustrate the concept, BMBPT is employed as a method to preprocess 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] | ≡ e T | , where the connected quasiparticle cluster operator T ≡ T 1 +T 2 + T 3 + . . . is defined through etc. The BCC amplitudes t 2n0 k 1 ···k 2n 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 (IT) 0 and a residual BMBPT correction δ (res) 0 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 loworder BMBPT estimates enter into play.

Importance-Truncated Tensor
Considering the n-tuple BCC amplitude tensor T n ≡ {t 2n0 k 1 ...k 2n }, the corresponding importance-truncated tensor based on the IT measure κ (p) (t 2n0 k 1 ...k 2n ) 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.

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 singleparticle 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 (eff) max ≤ 5, i.e., converged BCCSD results could be obtained using an effective model space including less than six major shells. 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 100 TB of memory. Performing the IT-compression of T 3 , the corresponding fourth-order contribution to the groundstate energy [4 T 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 [4 T ] 0 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.

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.

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 beyondmean-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 centerof-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-ofmass contamination and a slower convergence in large-scale NCSM diagonalizations compared to HO orbitals.

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 onebody 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.

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] ρ = ρ (00) + ρ (02) + ρ (20) + ρ (11) 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].

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. 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 tohω 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 loworder 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 configurationdriven methods like NCSM.

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].

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 ref max = 2 typically provide a good agreement with nonperturbative methods and that ground-state energies for larger reference spaces typically fall between N ref max = 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 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.
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 manybody and interaction uncertainties. The comparison reveals an excellent agreement between NCSM-PT and IM-NCSM, which themselves agree with experiment within estimated uncertainties.

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.
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 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 (q), second-order BMBPT ( ), and third-order BMBPT ( ). All calculations are performed using 15 oscillator shells and an oscillator frequency ofhω = 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. 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 ofhω = 20 MeV. The single-particle orbital angular-momentum quantum number was limited to l ≤ 10.
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., groundstate 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 nucleusindependent error of 3% is assumed. A recent study of highorder 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.
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.

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. Multiconfigurational perturbation theory and Bogoliubov MBPT were shown to provide complementary ways to tackle the particlehole degeneracy and the associated strong correlations in openshell 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 infrareddivergences, 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), socalled 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 mediummass 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.