A Guided Tour of ab initio Nuclear Many-Body Theory

Over the last decade, new developments in Similarity Renormalization Group techniques and nuclear many-body methods have dramatically increased the capabilities of ab initio nuclear structure and reaction theory. Ground and excited-state properties can be computed up to the tin region, and from the proton to the presumptive neutron drip lines, providing unprecedented opportunities to confront two- plus three-nucleon interactions from chiral Effective Field Theory with experimental data. In this contribution, I will give a broad survey of the current status of nuclear many-body approaches, and I will use selected results to discuss both achievements and open issues that need to be addressed in the coming decade.


INTRODUCTION
Over the past decade, the reach and capabilities of ab initio nuclear many-body theory have grown exponentially. The widespread adoption of Renormalization Group (RG) techniques, in particular the Similarity Renormalization Group (SRG) [1], and Effective Field Theory (EFT) [2,3,4] in the 2000s laid the foundation for these developments. Consistent two-nucleon (NN) and three-nucleon (3N) interactions from chiral EFT were quickly established as a new "standard" inputs for a variety of approaches, which made true multi-method benchmarks possible. The SRG equipped us with the ability to dial the resolution scale of nuclear interactions, accelerating model-space and many-body convergence alike. Suddenly, even (high-order) Many-Body Perturbation Theory (MBPT) became a viable tool for rapid benchmarking [5,6], and exact diagonalization approaches were able to extend their reach into the lower sd-shell [7,8,9]. A variety of of computationally efficient techniques with controlled truncations were readied, like the Self-Consistent Green's Function method (SCGF) [10], the In-Medium SRG (IMSRG) [11] and Coupled Cluster (CC) [12], the prodigal son [13,14] who returned home after finding success in foreign lands, i.e., quantum chemistry and solid state physics.
At the start of the last decade the race was on, and Fig. 1 documents the progress that ensued. Calculations started at closed-shell nuclei [15,16,17,18,19] and their vicinity before extending to semi-magic isotopic chains with the development of the Multi-Reference IMSRG [20,21] and Gor'kov SCGF [22,23] techniques, and just a couple of years later, the use of CC [24,25] and IMSRG [26,27] techniques to construct valence-space interactions opened all nuclei that were amenable to Shell Model calculations for exploration. Owing to very recent developments that extend these combined approaches to multi-shell valence spaces, the open region between the nickel and tin isotopic chain is poised to be filled in rapidly [28]. Development of the no-core versions of these methods has continued as well, and made direct calculations for intrinsically deformed nuclei possible [29].
The growing reach of ab initio many-body methods made it possible to confront chiral NN+3N forces with a wealth of experimental data, revealing shortcomings of those interactions and sparking new efforts toward their improvement. There were other surprises along the way, some good, some bad. Due to the benchmarking capabilities and further developments in many-body theory, we are now often able to understand the reasons for the failure of certain calculations (see, e.g., Ref. [27])hindsight is 2020, as they say 1 .
The present collection of Frontiers in Physics contributions provides us with a timely and welcome opportunity to attempt a look back at some of the impressive results from the past decade and the developments that brought us here, as well as a look ahead at the challenges to come as we enter a new decade.
Let us conclude this section with a brief outline of the main body of this work. In Section 2, I will discuss the main ingredients of modern nuclear many-body calculations: The input interactions from chiral EFT, the application of the SRG to process Hamiltonians and operators, and eventually a variety of many-body methods that are used to solve the Schrödinger equation. I will review key ideas but keep technical details to a minimum, touching only upon aspects that will become relevant again later on. Section 3 presents selected applications from the past decade, and discusses both Hergert A Guided Tour of Ab Initio Nuclear Many-Body Theory the advances they represent as well as open issues. This will provide a starting point for Section 4, which presents ideas for addressing the aforementioned issues and highlights important directions for the next decade.
Naturally, the discussion in Sections 3 and 4 is highly subjective. While this work grew from a more restricted scope into a rambling, albeit not random, walk through the landscape of modern nuclear many-body theory, it still cannot encompass the field in its entirety. The upside is that this reflects the breadth of ideas that are being pursued by the ab initio nuclear theory community, including those with cross-disciplinary impact, as well as our community's ability to attract junior researchers. The downside is that the present work can only scratch the tip of the iceberg of impressive results from the past decade. I hope that the readers will use it as a jumping-off point for delving into the cited literature, including the contributions to this volume.

Interactions from Chiral Effective Field Theory
Quantum Chromodynamics (QCD) is the fundamental theory of the strong interaction between quarks and gluons. One of its characteristic features is that the strong coupling, which governs the strength of interaction processes, is sufficiently small to allow perturbative expansions at high energies, but large in the low-energy domain relevant for nuclear structure and dynamics [30,31]. This makes the description of all but the lightest nuclei at the QCD level inefficient at best, and impossible at worst. However, strongly interacting matter undergoes a phase transition that leads to the confinement of quarks in composite hadronic particles, like nucleons and pions. These particles can be used as the degrees of freedom for a hierarchy of EFTs that describe the strong interaction across multiple scales. Following Weinberg [32,33], one can construct effective Lagrangians that consist of interactions that are consistent with the symmetries of QCD and organized by an expansion in (Q/Λ). Here, Q is a typical momentum of the interacting system, and Λ is the breakdown scale of the theory, which is associated with physics that is not explicitly resolved. In chiral EFT with explicit nucleons and pions, Λ = Λ χ is traditionally considered to be in the range 700 − 1000 MeV, although newer analyses of observable truncation errors using Bayesian methods favor slightly lower values [34,35,36]. From a chiral EFT Lagrangian, one can then construct a systematic low-momentum expansion of nuclear interactions, as shown in Fig. 2 (see Refs. [32,37,2,3,38]). These interactions consist of (multi-)pion exchanges between nucleons, indicated by dashed lines, as well as nucleon contact interactions. The different types of vertices are proportional to the low-energy constants (LECs) of chiral EFT, which encode physics that is not explicitly resolved because it involves either a high momentum scale or excluded degrees of freedom. Eventually, one hopes to calculate these LECs directly from the underlying QCD either through matching or renormalization group evolution of the couplings (see Section 2.2), but at present, the LECs are fit to experimental data [3,38,4,39,40].
The power counting scheme shown in Fig. 2 yields consistent two-, three-and higher many-nucleon interactions, and explains their empirical hierarchy, i.e., V NN > V 3N > V 4N > . . .. Moreover, one can readily extend the chiral Lagrangian with couplings to the electroweak sector by gauging the derivatives. In this way, nuclear interactions and electroweak currents depend on the same LECs, and one can use electroweak observables to constrain their values [42,43,44,45]. Last but not least, the existence of a power counting scheme offers inherent diagnostics for assessing the theoretical  . Chiral two-, three-and four nucleon forces through next-to-next-to-next-to-leading order (N 3 LO) (see, e.g., [37,41,2] ). Dashed lines represent pion exchanges between nucleons. The large solid circles, boxes and diamonds represent vertices that are proportional to low-energy constants (LECs) of the theory (see text).

The Similarity Renormalization Group
Renormalization group methods are a natural companion to the hierarchy of EFTs for the strong interaction. They provide the means to systematically dial the resolution scales and cutoffs of these theories, and this makes it possible, at least in principle, to connect the different levels in our hierarchy of EFTs. The RGs also expand the diagnostic toolkit for assessing the inherent consistency of EFT power counting schemes, e.g., by tracing the enhancement or suppression of specific operators, or by identifying important missing operators.
In nuclear many-body theory, the SRG has become the method of choice. In contrast to Wilsonian RG [52], which is based on decimation, i.e., integrating out high-momentum degrees of freedom, SRGs decouple low-and high-momentum physics using continuous unitary transformations. Note that this concept is not limited to RG applications: we can construct transformations that adapt a Hergert A Guided Tour of Ab Initio Nuclear Many-Body Theory many-body Hamiltonian or other observables of interest to our needs, e.g., to extract eigenvalues [11,53], or impose specific structures on the operator [1,26,54,27,55].
We define the flowing Hamiltonian where H(s = 0) is the starting Hamiltonian, and the flow parameter s parameterizes the unitary transformation. Instead of making an ansatz for U (s), we take the derivative of Eq. (1) and obtain the operator flow equation where the anti-Hermitian generator η(s) is related to U (s) by We can choose η(s) to achieve the desired transformation of the Hamiltonian as we integrate the flow equation (2) for s → ∞. Wegner [56] originally proposed a class of generators of the form that is widely used in applications, although it gives rise to stiff flow equations, and more efficient alternatives exist for specific applications [1,11,53]. Wegner generators are constructed by splitting the Hamiltonian into suitably chosen diagonal (H d (s)) and off-diagonal (H od (s)) parts. These labels are a legacy of applying this generator to drive finite-dimensional matrices towards diagonality. For our purposes, they reflect the desired structure of the operator in the limit s → ∞: We want to keep the diagonal part and drive H od (s) to zero by evolving it via Eq. (2) (see Refs. [56,57,1,11,53]).
To implement the operator flow equation (23), we need to express η(s) and H(s) in a basis of suitable operators {O i } i∈N , where η i (s) and H i (s) are the running couplings of the operators. If the algebra of the operators O i is closed naturally or with some truncation, we have and Eq. (2) becomes a system of flow equations for the coupling coefficients:

A Guided Tour of Ab Initio Nuclear Many-Body Theory
where the bold quantities collect the algebra's structure constants and the running couplings, respectively. From this discussion, it is clear that the choice of the O i can have a significant effect on the size of the system of flow equations, as well as the quality of any introduced truncations.
An important application of the SRG in nuclear many-body theory is the dialing of the operators' resolution scales. This is achieved by using the Wegner-type generator to band-diagonalize the Hamiltonian in momentum space, and thereby decouple low-and highmomentum physics in the operators and eigenstates. As indicated in Eq. (9) the flow is typically re-parameterized by λ = s −1/4 , which characterizes the width of the band in momentum space and controls the magnitude of the momentum transferred in an interaction process. For example, |k i − k f | λ in a two-nucleon system [1,58].
Nowadays, the momentum space evolution is regularly performed for two-and three-nucleon forces [59,1,60,61,62]. In light of the previous discussion, it can be understood as choosing the operator basis with creation and annihilation operators referring to (discretized) single-particle momentum modes, and truncating four-and higher-body terms that appear when the commutators of the basis operators are evaluated. Since the commutator of an M -body and an N -body operator in the basis (10) acts at least on K = max(M, N ) particles, the SRG evolution is exact for A ≤ 3 systems under this truncation [59,61]. It is implemented by working with the matrix representations of H(s) in two-and three-nucleon systems, whose entries correspond to the coupling constants in our chosen operator basis (cf. Eq. (6)). For efficiency, an additional basis change is made to center-of-mass and relative coordinates.
In principle, the strategy for evolving nuclear interactions towards some form of "diagonality" could be used to determine eigenvalues of many-body Hamiltonians, but the computational cost for dealing either with exponentially growing matrix representations or induced terms of high particle rank is prohibitive. This motivates the implementation of the flow equation with a different choice of basis operators in the In-Medium SRG (see Section 2.3.3).

Many-Body Methods
Let us now discuss commonly used many-body methods for solving the nuclear Schrödinger equation. Roughly speaking, they fall into two categories: configuration space methods that expand the nuclear eigenstates on a basis of known many-body states, or coordinate-space methods that work directly with the wave function and optimize them in some fashion. Our goal is to use approaches that systematically converge to an exact result, e.g., by adding more and more particle-hole excitations of a selected reference state to the many-body basis of a configuration space, or by exhausting the distribution of meaningful wave function parameters.
The discussion in the following sections will be light on mathematical details, which can be found in more specialized articles and reviews, including other contributions to the present volume. The goal is to review only certain ideas that will become relevant later on.

The Many-Body Problem in Configuration Space
Let us briefly discuss the general setup of the configuration-space approaches. We choose a single-particle basis, e.g., the eigenstates of a harmonic oscillator, and use it to construct a basis of Slater determinants for the many-body Hilbert space. Usually, the many-body basis is organized by selecting a reference state |Φ and constructing its particle-hole excitations in order to account for the natural energy scales of the system under consideration. For further use, we define where particle (a, b, . . .) and hole (i, j, . . .) indices run over unoccupied and occupied single-particle states, respectively 2 . The parentheses indicate that the strings of creation and annihilation operators are normal ordered with respect to the reference state. They are related to the original operators by where the indices p, q, . . . run over all single-particle states, and the contractions are defined as (see, e.g., Refs. [11,53] for more details).
Let us now consider a Hamiltonian containing up to two-body interactions, for simplicity. In normal-ordered form, it is given by

A Guided Tour of Ab Initio Nuclear Many-Body Theory
where E is the energy expectation value of the reference state, while f and Γ are the mean-field Hamiltonian and residual two-body interaction, respectively [11,53]. Our task is to solve the many-body Schrödinger equation for this Hamiltonian to determine its eigenvalues and eigenstates, either in an approximate fashion or by exactly diagonalizing its matrix representation, which is shown in Fig. 3(a).

Many-Body Perturbation Theory
Many-Body Perturbation Theory (MBPT) is the simplest configuration-space approach for capturing correlations in interacting quantum many-body systems. It has enjoyed widespread popularity in treatments of the many-electron system since the early days of quantum mechanics, and it comes in a myriad of flavors (see, e.g., Ref. [64] and references therein). A major factor in its success is that the Coulomb interaction is sufficiently weak to make perturbative treatments feasible. Applications in nuclear physics had long been hindered by the strong short-range repulsion and tensor interactions in realistic nuclear forces, despite the introduction of techniques like Brueckner's G matrix formalism that were meant to resum the strong correlations from these contributions [65,66,67,68]. These issues were overcome with the introduction of the SRG evolution to low resolution scales, which makes nuclear interactions genuinely perturbative, albeit at the cost of inducing three-and higher many-body interactions [1]. As a consequence, MBPT has undergone a renaissance in nuclear physics in the past decade [69], leading to efficient applications for the computation of ground-state properties [5,70,6] and the construction of effective Shell Model interactions and operators (see, e.g., Refs. [71,72,73,74], or the reviews [75,76] and references therein). These successes have also motivated the development of novel types of MBPTs [77,78,69].
In a nutshell, MBPT assumes that the Hamiltonian can be partitioned into a solvable part H 0 and a perturbation H I , which then allows an order-by-order expansion of its eigenvalues and eigenstates in powers of H I , usually starting from a mean-field solution. In the Rayleigh-Schrödinger formulation of MBPT, which is widely used for its convenience, where E (0) is the unperturbed energy. If we assume that the reference Slater determinant |Φ has been variationally optimized by solving the Hartree-Fock equations, E 0 in Eq. (15) is the Hartree-Fock energy and f is diagonal. Then we can introduce the so-called Møller-Plesset partitioning, and note that the Slater determinants of the basis introduced in Section 2.3.1 are eigenstates of H 0 :

A Guided Tour of Ab Initio Nuclear Many-Body Theory
The eigenvalues of H 0 then become the unperturbed energies appearing in Eqs. (17), (18), and the energy including a finite number of correction terms can be evaluated straightforwardly. For example, the ground-state energy through second order is given by For a more detailed discussion, we refer to Ref. [69] and references therein.
The expression (21) can serve to illustrate both advantages and drawbacks of an MBPT treatment of nuclei. We see that the second-order energy can be evaluated very efficiently, since it requires a non-iterative calculation whose computational effort scales polynomially in the single-particle basis size N , namely as O(N 4 ). The reason is that the construction of the Hamiltonian matrix ( Fig. 3(a)) can be avoided. In fact, the computational scaling is even more favorable, because we can distinguish particle and hole states and achieve O(N 2 p N 2 h ), and we typically have N h ∼ A N p . Although there is a proliferation of terms with increasing order [79,63,69], MBPT is still fundamentally polynomial and therefore more efficient than an exact diagonalization, whose cost scales exponentially with N . It is also clear from Eq. (21) that the expansion of the exact eigenvalue will break down if one (or more) of the energy denominators become small due to (near-)degeneracies of the unperturbed energies. Thus, MBPT works best for ground states in systems with a strong energy gap, i.e., closed-shell nuclei, although extensions for more complex scenarios exist (see Refs. [68,63,69] and references therein). A noteworthy new development is Bogoliubov MBPT, in which particle number symmetry is broken and eventually restored [77,80,81].
As mentioned at the beginning of this section, MBPT can be used to derive effective interactions and operators. The primary tool for such efforts is theQ-box or folded-diagram resummation of the perturbative series (see Refs. [82,75,76] and references therein).

In-Medium Similarity Renormalization Group
As already mentioned in our discussion of the SRG in Section 2.2, we could envision applying SRG techniques not only to preprocess the nuclear interactions, but also to compute eigenvalues and eigenstates. For all but the lightest nuclei, applying the SRG to the Hamiltonian matrix is hopeless, so we work with the operators instead.
Let us again consider the matrix representation shown in Fig. 3 (a). We want to design a transformation that will decouple the one-dimensional 0p0h block in the Hamiltonian matrix, spanned by a reference state Slater determinant |Φ , from all excitations as the flow equation (2) is integrated. The matrix element in this block will then be driven towards an eigenvalue (up to truncation errors), and the unitary transformation becomes a mapping between the reference Slater determinant and the exact eigenstate (see below). In principle, we could use a suitably chosen reference to target different eigenstates, e.g., by taking references which are expected to have a large overlap with the target state (see Section 10.3 in Ref. [58]). In practice, we usually target the ground state by using a Hartree-Fock Slater determinant as our reference.
To implement the operator flow, we need to choose an operator basis to express H(s) and the generator η(s). Instead of using the basis (10), we switch to operators that are normal ordered with Commutators of these operators can feed into terms of lower particle rank: For instance, a commutator of M -body and N −body operators generates |M − N |-body through (M + N − 1)-body operators, while the lower bound for the basis (10) is max(M, N ) (cf. Section 2.2). As a result, the complexity of the flow equations for the operators' coupling coefficients increases due to the appearance of additional terms that depend on the contractions introduced in Eqs. (12) and (13). These contractions translate into density matrices (or occupation numbers) -hence the name In-Medium SRG. At the same time, we achieve a reduction of the truncation error because only the residual, contractionindependent parts of the operators (12) and (13) are omitted. In the majority of applications to date, we truncate all operators and their commutators at the two-body level, defining the IMSRG(2) truncation scheme. More details can be found in Refs. [11,53,58,76].
In the chosen basis we now identify the parts of the Hamiltonian that are responsible for coupling the reference state to 1p1h and 2p2h excitations, and define the off-diagonal Hamiltonian (cf. 2.2) as We use this H od to construct a generator, either using Wegner's ansatz (4) or an alternative choice [11,53]. Plugging the generator into the operator flow equation (2), we obtain a system of flow equations for the energy E(s) and the coefficients f pq (s), Γ pqrs (s), . . . (cf. Eq. (8) and Refs. [11,53,76]). By integrating these flow equations, we evolve the Hamiltonian operator so that its matrix representation assumes the shape shown in Fig. 3 (b). We note that the suppression of H od not only leads to the desired ground-state decoupling, but also eliminates the outermost band in the Hamiltonian matrix. This simplification makes the evolved Hamiltonian an attractive input for other approaches, e.g., configuration interaction (CI) or equation-of-motion methods (see Refs. [83,27,76,84,85,86,29] and discussion below).
Valence-Space IMSRG. Soon after introducing the IMSRG in nuclear physics [87], Tsukiyama, Bogner and Schwenk proposed the use of the IMSRG flow to derive Hamiltonians (and other effective operators) for use in nuclear Shell Model calculations [88]. This is achieved by partitioning the single-particle basis into core, valence, and beyond-valence states, normal ordering all operators with respect to a Slater determinant describing the closed-shell core, and extending the definition of the off-diagonal Hamiltonian (23) to include all terms that couple valence and non-valence states. The eigenvalue problem for the evolved Hamiltonian can then be solved in the valence space with widely available Shell model codes [89,90,91,92,93]. After a study of the oxygen isotopic chain revealed an increasing overbinding away from the chosen core [26], we adopted a normal-ordering scheme that uses an ensemble of Slater determinants to account for partially filled shells in open-shell nuclei [54,27]. This improved operator basis, along with the valence decoupling procedure and subsequent Shell Model diagonalization defines what is nowadays called the valence-space IMSRG (VS-IMSRG) -see Ref. [76] for a recent review.
The transformation shifts correlations from the wave function into the evolved, RG-improved Hamiltonian H(s) = U (s)HU † (s), and any many-body method that uses this Hamiltonian as input now needs to describe U (s) |Ψ k , which should be less correlated than the exact eigenstate |Ψ k . In the extreme cases, U (s) = 1 and the wave function carries all correlations, or U (s) has shifted all correlations into the Hamiltonian and |Φ = U (s) |Ψ is a simple Slater determinant.
Correlated reference states can be particularly useful for the description of systems with strong static or collective correlations, like open-shell nuclei with strong intrinsic deformation or shape coexistence. Reference states that describe these types of correlations efficiently, e.g., through symmetry breaking and restoration (also see Section 2.3.4), are an ideal complement to the IMSRG transformation, which excels at capturing dynamic correlations, involving the excitation of a few particles up to high energies. This complementarity is schematically illustrated in Fig. 4: Collective correlations that would require as much as an IMSRG(A) calculation in the conventional approach are built into the reference state, and an MR-IMSRG(2) calculation is sufficient to treat the bulk of the dynamical correlations in the system.
Reference state correlations are built into the MR-IMSRG framework by using a generalized normal ordering [94,95,53] that is extended with contractions of higher rank, namely the irreducible k-body density matrices λ (k) : λ pqrs ≡ ρ pqrs − ρ pr ρ qs + ρ qr ρ ps , Hergert A Guided Tour of Ab Initio Nuclear Many-Body Theory etc. The irreducible densities matrices encode the correlation content of an arbitrary reference state |Φ , hence they vanish for Slater determinants. While the basis of normal-ordered operators superficially is the same as in the conventional IMSRG, shown in Eq. (22), the inclusion of the irreducible densities (cf. Eqs. (12) and (13)) equips the basis with the capability to describe the correlations that are present in the reference state, which in turn should help to reduce MR-IMSRG truncation errors. To understand this, let us assume that we know the ground state of our system, and we normal order the Hamiltonian with respect to this correlated state. Then the zero-body part of the normal ordered Hamiltonian already is the exact ground-state energy, and the normal-ordered one-, two-and higher-body parts do not matter at all for our result, and neither does their evolution under an exact or truncated MR-IMSRG flow. Thus, the better the reference state matches the ground state, the less work the MR-IMSRG evolution and any subsequent many-body method have to do to obtain the correct ground-state energy.
Computational Scaling and Magnus Expansion. The computational scaling of all three IMSRG flavors discussed here -traditional, VS-IMSRG, and MR-IMSRG -is governed by the truncation scheme. If we truncate operators and commutators at the two-body level, as briefly mentioned above, the number of flow equations scales as O(N 4 ) with the single-particle basis size N , and the computational effort for evaluating the right-hand sides as O(N 6 ). This holds despite the greater complexity of the MR-IMSRG flow equations, which contain terms containing irreducible two-and higher-body density matrices.
Any observables of interest must, in principle, be evolved alongside the Hamiltonian for consistency, which would create a significant overhead. In practice, we can address this issue by using the so-called Magnus formulation of the IMSRG [96,83,58,76]: Assuming that the IMSRG transformation can be written as an explicit exponential, U (s) = exp Ω(s), we can solve a single set of flow equations for the anti-Hermitian operator Ω(s) instead of evolving observables separately. All operators of interest can then be computed by applying the Baker-Campbell-Hausdorff expansion to IMSRG Hybrid Methods. As noted earlier in this section, the conventional IMSRG evolution makes the matrix representation of the Hamiltonian more diagonal by suppressing couplings between the npnh excitations of the reference state. This implies a decoupling of energy scales of the manybody system, analogous to the decoupling of momentum scales by the free-space SRG, although there are differences in detail that are associated with the operator bases in which the flow is expressed (cf. Eqs. (10), and (22)).
From this realization, it is not a big step to consider using the IMSRG to construct RG-improved Hamiltonians for applications in other methods, defining novel hybrid approaches. In fact, even the original IMSRG formulation can be understood from this perspective: The evolution generates a Hamiltonian that yields the exact ground-state energy (up to truncations) in a Hartree-Fock calculation, except the HF equations are automatically satisfied for the evolved H, and we can read off the ground-state energy directly. The same Hamiltonian can then be used as input for EOM methods to compute excitation spectra [83]. Likewise, the VS-IMSRG produces an RG-improved Hamiltonian that serves as input for a Shell Model diagonalization.
Applying the same logic as in the VS-IMSRG case, the IMSRG has been merged with the No-Core Shell Model (NCSM, see Section 2.3.6) into the In-Medium NCSM [84,97]. In this approach, the IMSRG improves the Hamiltonian with dynamical correlations from high-energy few-nucleon Hergert A Guided Tour of Ab Initio Nuclear Many-Body Theory excitations that would require enormously large model spaces in the conventional NCSM, and the exact diagonalization in a small model space describes the dynamics of many-nucleon excitations. The NCSM as the "host" method is rooted in the same particle-hole expansion picture as the IMSRG itself, but this is not a requirement. Another new hybrid method is the In-Medium Generator Coordinate Method (IM-GCM), which relies on the GCM as a host method to capture collective correlations [85,86,29]. In this approach, a many-body basis is generated by restoring the symmetries of mean field solutions with various types of shape and gauge configuration constraints, which is very different from the particle-hole excitation basis discussed so far.

Coupled Cluster Methods
The Coupled Cluster (CC) method [63,12] is an older cousin of the IMSRG approach. It can also be understood as a decoupling transformation of the Hamiltonian, but in contrast to the IMSRG, it relies on a non-unitary similarity transformation (see Fig. 3). Traditionally, CC is motivated by an exponential ansatz for the exact wave function of a system, where |Φ is a reference Slater determinant, and T is the so-called cluster operator. This operator is expanded on particle-hole excitations, with the cluster amplitudes t ai , t abij , . . .. In practical applications, the T is truncated to include up to 2p2h (CC with Singles and Doubles, or CCSD) or 3p3h terms (CCSDT, including Triples). Various schemes exist for iteratively or non-iteratively including subsets of Triples [98,99,63,100,12]. When it acts on the reference state |Φ , e T admixes arbitrary powers of few-particle, few-hole excitations. Note, however, that the cluster operator T is not anti-Hermitian because it lacks de-excitation operators, and therefore e T is not unitary.
The cluster amplitudes are determined by demanding that the transformed Hamiltonian, does not couple the reference to 1p1h and 2p2h states (see Fig. 3). Using notation introduced in Section 2.3.1, the decoupling conditions lead to the following system of non-linear equations: Here, E CC is the CC ground-state energy, which corresponds to the one-dimensional block in the upper left of Fig. 3 (c) and is analogous to the zero-body part of the IMSRG-evolved Hamiltonian, as discussed in the previous section. The other blocks in the first column of the matrix vanish because of the CC equations (30)-(32).

A Guided Tour of Ab Initio Nuclear Many-Body Theory
Since the CC transformation is non-unitary, one needs to be careful when one evaluates observables using the CC wave function, or uses H CC as input for equation-of-motion calculations or other applications [63,12]. For instance, the non-Hermiticity of H CC forces us to consider left and right eigenstates separately. This is a drawback compared to unitary transformation methods like the IMSRG. Coupled Cluster also has advantages, though: For instance, the Baker-Campbell-Hausdorff expansion appearing in Eqs. (30)-(32) automatically terminates at finite order because the cluster operator only contains excitation operators. For the same reason, Eq. (31) will automatically solve the Hartree-Fock equations, so any Slater determinant is equally well suited as a reference state, while MBPT, IMSRG and even exact diagonalization approaches exhibit (some) reference-state dependence.
Symmetry Breaking and Collective Correlations. While most applications of CC theory in nuclear physics have enforced and exploited spherical symmetry, the capabilities for performing M -scheme calculations that allow nuclei to develop intrinsic deformation have existed for more than a decade. This is a more natural approach for capturing collective correlations than the construction of Triples, Quadruples (4p4h) and ever higher particle-hole excitations of a spherical reference (cf. Section 2.3.3). Converging such calculations is challenging because the single-particle basis typically grows by an order of magnitude or more, and the broken symmetries must eventually be restored. The formalism for symmetry restoration in CC has been developed in Refs. [101,102,103,104]. In fact, the work of Duguet et al. forms the basis of recent works on symmetry breaking and restoration in MBPT [77,80,81]. Applications are currently underway.
Shell-Model CC. Like the IMSRG, the CC framework can be used to construct effective interactions and operators for Shell model calculations. Initial work in that direction applied Hilbert space projection techniques (cf. Section 2.3.6) to construct a so-called CC effective interaction (CCEI) [24,105], but the construction of the model spaces via Equation-of-Motion CC methods proved to be computationally expensive. The CCEI approach is now superseded by the Shell Model CC method [25], which applies a second similarity transformation to H CC in Fock space, similar to VS-IMSRG decoupling (cf. Section 2.3.3).
Unitary CC. While almost all applications of CC in nuclear physics use the traditional ansatz (27), unitary CC (UCC) approaches that parameterize the wave function as |Ψ U CC = e T −T † |Φ have been used in numerous studies in quantum chemistry (see, e.g., [106,107]). Unitary CC wave functions have also become a popular ansatz for the Variational Quantum Eigensolver (VQE) algorithm on current and near-term quantum devices [108,109]. It is also worth noting that the recently revived Unitary Model Operator Approach (UMOA) is closely related to UCC [110,111].

Self-Consistent Green's Functions
Self-Consistent Green's Function (SCGF) theory is another prominent approach for solving the nuclear many-body problem with systematic approximations [112,113,114,115]. The Green's Functions in question are correlation functions of the form which describe the propagation of nucleons in the exact ground state |Ψ A 0 of the system. Using Wick's theorem, the exact A-body propagator (33) can be factorized into products of irreducible one-, two-, etc. propagators, similar to the decomposition of density matrices briefly touched upon Hergert A Guided Tour of Ab Initio Nuclear Many-Body Theory in Section 2.3.3. One can then formulate coupled equations of motion for propagators, and introduce truncations to obtain polynomially scaling methods, again somewhat analogous to IMSRG and CC. We must remain aware that the propagators of SCGF, the induced operators of IMSRG, and the CC amplitudes are all different objects, and while their definitions may make the seem complementary to each other, there are subtle distinctions. One of these is that the g (k) are formally defined with respect to the exact wave function, while IMSRG and CC use definitions with respect to a reference state.
Practical implementations of the SCGF technique usually work with the Fourier transforms of the propagators to the energy domain. One needs to solve integral equations of motion of the form where g 0 is the propagator of the non-interacting system and Σ a kernel that encodes the particles' interactions, which is constructed using diagrammatic techniques. For example, the one-body propagator is obtained by solving the so-called Dyson equation. From this propagator, one can compute the one-body density matrix where C + indicates an integration contour in the complex upper half plane. Higher-body density matrices are connected to the corresponding higher-body propagators in analogous fashion. Using the density matrices, one can then evaluate any operator expectation values of interest. For more details, we refer to the contributions [10,115] to the present volume, and the works cited therein.
Current applications of SCGF techniques in nuclear physics make use of the so-called Algebraic Diagrammatic Construction (ADC) scheme, with increasing orders, denoted by ADC(n), converging to an exact solution. For closed-shell nuclei, calculations up to ADC(3) are be performed regularly, which contain correlations that are roughly comparable to IMSRG(2) with a perturbative 3p3h correction (see Section 2.3.3 and Refs. [116,83,86]) and CCSD(T) (cf. Section 2.3.4). Somà and collaborators have extended the ADC scheme to open-shell nuclei by using Gor'kov Green's Functions with explicitly broken particle number symmetry [117,118]. Applications of this framework have used a self-consistent second-order scheme, denoted Gor'kov-ADC (2), and the extension to Gor'kov-ADC(3) as well the integration of particle-number projection to restore the broken number symmetry are in progress [80,114].
While the computation of the Green's Functions tends to be a more involved task than solving the IMSRG flow equations or CC amplitude equations, the propagator contains more information from a single computation than these other methods. For instance, one can immediately extract spectral information about the neighboring nuclei and the response of the system [119,120], which requires the application of additional techniques in the IMSRG [83] and CC approaches [121,122,12], or, indeed, the computation of the Green's Function using similarity-transformed operators. Furthermore, the kernels of the equations of motion (34) are energy-dependent effective interactions that govern the dynamics of (few-)nucleon-nucleus interactions. For example, the one-nucleon self-energy in Eq. (35)

A Guided Tour of Ab Initio Nuclear Many-Body Theory
is an ab initio version of an optical potential, as used in reaction theory [123,124,125]. We will return to this discussion in Section 4.5.

Configuration Interaction Approaches
No-Core Configuration Interaction Methods. The most straightforward but also most computationally expensive approach to solving the many-body Schrödinger equation is to exactly diagonalize the Hamiltonian in a basis of many-body states. In general, we refer to such approaches as No-Core Configuration Interaction (NCCI). "No core" makes it explicitly clear that all nucleons are treated as active degrees of freedom, in contrast to the nuclear Shell model discussed below.
In light nuclei, the exact diagonalization can be directly formulated in Jacobi coordinates, using translationally invariant harmonic oscillator [126] or hyperspherical harmonic wave functions [127,128]. Since the construction of the basis states themselves and the matrix representation of the Hamiltonian becomes increasingly complicated and computationally expensive as the particle number grows, one eventually has to switch to Slater determinants in the laboratory system, using a construction along the lines discussed in Section 2.3.1.
A common choice for the single-particle basis in the laboratory system are spherical harmonic oscillator (SHO) states, because they allow an exact separation of center-of-mass and intrinsic degrees of freedom provided one uses an energy-based truncation for the model space [129,130].
These choices define what we specifically call the No-Core Shell Model (NCSM). A disadvantage of using SHO orbitals is that they are not optimized to the energy scales of specific nuclei, and they are poorly suited for describing physical features like extended exponential wave function tails. Other popular choices are Hartree-Fock single-particle states, and perturbatively [131] or nonperturbatively enhanced natural orbitals [132,133,134]. Model spaces built on these bases no longer guarantee the separation of center-of-mass and intrinsic coordinates, but fortunately, center-of-mass contaminations either remain small automatically [135], or they can be suppressed using techniques like the Lawson method [136].
Importance Truncation and Symmetry Adaptation. As indicated above, the main issue with exact diagonalization approaches is the exponential (or greater) growth of the Hilbert space dimension, which is proportional to N A with single particle basis size N and particle number A. A variety of strategies can be used to address this often-quoted "explosion" of the basis size. One direction is to avoid the construction of the full model space basis by applying importance-based truncation or sampling methods, leading to the Importance-Truncated NCSM [9] or Monte-Carlo (No-Core) CI approaches [137,138]).
Another important research program is the exploration of many-body states that are constructed from the irreducible representations (irreps) of the symplectic group Sp(3,R), which describes an approximate emergent symmetry of finite nuclei [139,140]. An exact diagonalization in such a symmetry-adapted basis will offer a much more efficient description of nuclear states with intrinsic deformation than the conventional NCSM, which would need to use massive model spaces with many-particle-many-hole excitations. This reduction of the model space dimensions also allows such symmetry-adapted NCSM [139,140] and NCCI approaches [141] to reach heavier nuclei than the conventional versions.
Interacting Nuclear Shell Model with a Core (Valence CI). Instead of treating all of the nucleons as active, one can also factorize the nuclear wave function by introducing an inert core and This, of course, is the traditional nuclear Shell model approach. Even with the substantial reduction of the single-particle basis to a relatively small number of valence orbitals, the numerical cost for an exact diagonalization quickly becomes unfeasible for many medium-mass and heavy nuclei, especially if one needs multi-shell valence-spaces to capture complex nuclear structure features like coexisting intrinsic shapes.
In previous sections, we have discussed how a variety of many-body methods can be used to derive valence-space interactions, hence it is not a surprise that this is possible in NCCI approaches as well. One strategy is to project solutions of no-core calculations for the core and its neighboring nuclei onto a valence-configuration space to extract the effective Hamiltonian. The viability of this approach has been demonstrated in several publications [142,143,144,145], although there are ambiguities in the extraction of the valence-space Hamiltonian, and the initial NCCI calculations that serve as input for the projection rapidly become expensive.

Description of Continuum Effects and Nuclear
Dynamics. An important breakthrough in ab initio calculations for light nuclei has been the merging of the NCSM with resonating group method (RGM) techniques [130,146]. This makes it possible to describe clustered states as well as reactions between light projectile(s) and targets. In the original NCSM/RGM approach, compact clusters of nucleons are described by NCSM states, which are then used to construct a basis of configurations |χ i that place such clusters at different relative distances. In this basis, one can then solve the generalized eigenvalue problem, known as the Griffin-Hill-Wheeler equation [147] in the RGM context: where H and N are the so-called Hamiltonian and norm kernels. The latter appears because the chosen basis configurations are not orthogonal in general. The dimension of Eq. (38) is typically small, certainly compared to the NCSM model space, but the computation of the kernels is computationally expensive since it relies on the construction of up to three-body transition density matrices. In recent years, the NCSM/RGM has been extended to the NCSM with Continuum (NCSMC), which accounts for the coupling between the NCSM and RGM sectors of the many-body basis [130]. It requires solving the generalized eigenvalue problem where h and 1 are the Hamiltonian and norm kernel in the NCSM sector (the latter being diagonal), H and N the corresponding kernels in the RGM sector (cf. Eq. (38)), andh andn encode the coupling between the sectors of the basis.
Alternative approaches to the description of continuum effects in the NCSM are the Single-State HORSE (Harmonic Oscillator Representation of Scattering Equations) method [148,149,150], for which the nomen is omen, as well as the No-Core Gamow Shell Model (GSM), a no-core CI approach that constructs Slater determinants from a single-particle Berggren basis [151] consisting of bound, resonant and scattering states [152,153,154,155].

Quantum Monte Carlo
The most commonly used Quantum Monte Carlo (QMC) techniques in nuclear physics make use of many-body wave functions in coordinate space representation [156,157,158,159]. As such, they are well suited for the description of nuclear states with complex intrinsic structures, and they can readily use interactions with a high momentum cutoff, as opposed to the configuration space methods which would exhibit poor convergence in such cases. This allows QMC calculations to explore physics across the interfaces of the hierarchy of EFTs for the strong interaction (cf. Sections 2.1 and 4.4), e.g., for processes that explore energies approaching the breakdown scale of chiral EFT [160,161,162,163].
A typical ansatz for a QMC trial state is where F(a) is an operator that explicitly imprints correlations on the mean-field like state |Φ(b) , and a, b are vectors of tunable parameters. The first step of most QMC calculations is a variational minimization of the energy in the trial state , followed by an imaginary-time evolution to project out the true ground state in a quasi-exact fashion: This projection can be implemented using Monte Carlo techniques in a variety of ways, which gives rise to different approaches like Green's Function Monte Carlo (GFMC) or Auxiliary-Field Diffusion Monte Carlo (AFDMC) [156,158].
A major challenge in QMC calculations is that most commonly used algorithms suffer from some form of sign problem [156,158]. Many quantities of interest like the wave functions or local operator expectation values in these wave functions are not positive definite across their entire domain, which means that they cannot be immediately interpreted as probability distributions that the algorithms sample. This is one of the main reasons why QMC methods can only be used with Hamiltonians that are either completely local, or have a nonlocality that is at most quadratic in the momenta, e.g., p 2 or l 2 .
While QMC applications in ab initio nuclear structure have been focused on coordinate space, there are a wide variety of approaches that merge QMC techniques with the configuration space approaches discussed in previous sections. Examples include sampling the intermediate-state summations in MBPT [164], diagrammatic expansions [165,166,167], or the coefficients of correlated CC [168] or (No-Core) CI wave functions [137,138,169,170,171].

Lattice Effective Field Theory
Lattice methods are nowadays widely used to simulate the dynamics of nonperturbative field theories on finite space-time lattices. The most prominent example is Lattice QCD, but implementations of various Effective Field Theories on the Lattice have been developed and applied with impressive Hergert A Guided Tour of Ab Initio Nuclear Many-Body Theory outcomes in the past two decades -see, for example, Refs. [172,173,174,175,173] and references therein, which also provide pedagogical introductions to Lattice EFT for nuclear systems.
Lattice EFT simulations are built around the partition function, which is defined for a pure state |Ψ as Here, H is an EFT Hamiltonian, typically truncated at a given order of the EFT's power counting scheme. In practice, the partition function is evaluated as a path integral in which field configurations are sampled using Monte Carlo techniques. At large τ , one can extract information about the ground state and low-lying excited states of the system directly from Z (cf. Section 2.3.7), and general expectation values can be evaluated using The use of discretized spatial lattices makes Lattice EFT particularly suited for the description of nuclear states with complex geometries like cluster structures [176,177,178]. Depending on the size of the lattices, it will also typically require less computational effort than the imaginary-time evolution of states that are formulated in continuum coordinates, as in AFDMC or GFMC (see Section 2.3.7). Moreover, the development of the so-called adiabatic projection method (APM) [179,180] in recent years has made it possible to compute scattering cross sections for reactions of (light) clusters on the lattice. Conceptually, the APM is reminiscent of the resonating-group method used to describe reactions in the NCSMC framework discussed in Section 2.3.6.
Of course, Lattice EFT is not free of disadvantages, which are usually caused by the discretization of space(time). The finite size and lattice spacing are related to infrared (long-range, low-momentum) and ultraviolet (short-range, high-momentum) cutoffs of a calculation, which need to be carefully considered. Since the recognition of cutoff scales is an inherent aspect of EFTs, one can systematically correct for these effects [181,182]. The discrete lattice also breaks continuous spatial symmetries that may need to be restored approximately or exactly before comparisons with experimental data are made [172,182].

THE PAST IS PROLOGUE: ACHIEVEMENTS IN THE LAST DECADE
In this section, I will discuss selected achievements of the ab initio nuclear many-body community in the past decade, and the issues that were encountered in the process. As stated in the introduction, this selection is subjective, and giving full justice to the breadth of research accomplishments is beyond the scope of this work. I hope that the present discussion will serve as an invitation for further exploration, for which the cited literature may serve as a useful starting point.

Benchmarking Nuclear Forces
One of the biggest issues in nuclear theory was the lack of comparability between different approaches for describing the structure of medium-mass or heavy nuclei. These nuclei were well in reach of the Shell Model and nuclear Density Functional Theory (DFT), but whenever issues emerged, it was unclear whether they resulted from approximations in the many-body method, or deficiencies in the effective interactions, i.e., the valence-space Hamiltonians or energy density  Figure 5. Ground-state energies of the oxygen isotopes for various many-body approaches, using the chiral NN+3N(400) interaction at λ = 1.88 fm −1 [183]. Details on the Lattice EFT calculation can be found in Ref. [177]. Gray bars indicate experimental data [184].
functionals (EDF). Moreover, one cannot simply perform a valence CI calculation with an EDF, or a DFT calculation with a Shell Model interaction, because the interactions are tailored to their specific many-body method.
The development of the RG/EFT and many-body methods discussed in Section 2 has opened up a new era for benchmarking the same nuclear interactions across multiple approaches, and on top of that, these methods provide a systematic framework for analyzing, and eventually quantifying, the reasons for differences between the obtained results.
One of the earliest testing grounds for ab initio calculations of medium-mass nuclei was the oxygen isotopic chain, which was accessible to all of the approaches that emerged at the beginning of the past decade. Figure 5 shows the ground-state energies of even oxygen isotopes for the same chiral NN+3N interaction, obtained with several of the configuration space approaches introduced in Section 2.3. In addition, results for applying various types of MBPT to the same interaction and nuclei are presented in Ref. [69] -I only refrained from including them here to avoid overloading the figure. As we can see, the ground-state energies obtained from the different approaches are in good agreement with each other and with experiment. Since our results include quasi-exact IT-NCSM values, the deviation of the other methods' energies from these values provide us with an estimate of the theoretical uncertainties due to any employed truncations, which is on the order of 1-2%. As we can see from Fig. 5, essentially all of the used many-body methods place the drip line in the oxygen isotopic chain at 24 O, although the signal is exaggerated. Continuum effects that have been omitted in these calculations would lower the energy of the 26 O resonance, which is experimentally constrained to be a mere 18(7) keV above the two-neutron threshold [185], and produce a very flat trend in the energies towards 28 O. Similar features were found in calculations for other isotopic chains and other chiral interactions [21,118,186,114]. The 16 O ground state energies obtained for the employed chiral NN+3N Hamiltonian are also compatible with a Lattice EFT result that was obtained at a similar resolution scale [177].
This last comparison shows that some obstacles to the ideal cross-validation scenario still remain. Since coordinate-space approaches like Lattice EFT or QMC are truly complementary Hergert A Guided Tour of Ab Initio Nuclear Many-Body Theory to configuration-space methods, it would be highly desirable to test the same chiral NN+3N Hamiltonians in both types of calculations. However, the Hamiltonians used in configuration space are typically given in terms of harmonic oscillator matrix elements (especially if SRG evolved) instead of the coordinate-space operators required by Lattice EFT or QMC calculations. Furthermore, Lattice EFT and QMC cannot handle all possible types of nonlocality in the Hamiltonian (cf. Section 2.3.7), including the forms generated by the nonlocal regulators that are favored for configuration-space Hamiltonians. Conversely, local chiral interactions that have been constructed explicitly for QMC applications [187,188,189,190,4,158] exhibit slow model-space convergence in configuration-space calculations because they still tend to require a significant repulsive core at short distance to describe nucleon-nucleon scattering data, albeit a far weaker one than interactions like Argonne V18 [191].

Extending the Reach of Ab Initio Theory
The reach of ab initio many-body theory has increased dramatically over the past decade. Figure  1 illustrates this growing coverage of the nuclear chart, but it tells only part of the story. The expansion has happened in many "dimensions" besides the mass number A, namely by pushing towards exotic nuclei via improved treatments of the continuum degrees of freedom, filling in gaps in the coverage that are occupied by doubly open-shell nuclei with strong intrinsic deformation, and expanding the types of observables that can be computed from first principles. Recalling Section 3.1, the ongoing push against the limitations of our many-body approaches will continue to grow the opportunities for benchmarking current-and next-generation chiral Hamiltonians.

Pushing the Mass Boundaries
First calculations for selected nuclei and semi-magic isotopic chains up to tin were already published in the first half of the last decade [19,21,23]. For the most part, they were using a family of chiral NN+3N interactions that gave a good description of the oxygen ground-state energies (cf. Fig. 5) as well as the spectroscopy of the lower sd-shell region [26,24]. However, the same interactions underpredict nuclear charge radii [192], and start to overbind as we approached the calcium chain (cf. Fig. 7), eventually leading to an overbinding of 1 MeV per nucleon in tin. While model-space convergence in CC, IMSRG and SCGF calculations suggested that calculations for heavier nuclei would have been technically possible, it made little sense to pursue them.
The growing number of results for medium-mass nuclei and the problems they revealed motivated a new wave of efforts to refine chiral interactions. One direction of research aimed to achieve a simultaneous description of nuclear energies and radii up to 48 Ca by including selected many-body data in the optimization protocol of the chiral LECs. This work resulted in the so-called NNLO sat interaction [193]. While NNLO sat definitely improved radii [194], its model-space convergence was found to become problematically slow already in lower pf -shell nuclei [195,196,114].
Simultaneously with the efforts to develop new interactions, attention also turned towards an older, less consistently constructed family of chiral NN+3N interactions that exhibited reasonable saturation properties in nuclear matter calculations [199,200]. These forces are referred to as EMλ/Λ, where λ indicates the resolution scale of the NN interaction, the SRG-evolved N 3 LO potential of Entem and Machleidt [201], and Λ is the cutoff of an NNLO three-nucleon interaction whose low-energy constants have been adjusted to fit the triton binding energy and 4 He charge radius [199,200]. Ni Figure 6. Energies of the first excited 2 + states from VS-IMSRG [197] and Equation-of-Motion CC [195] calculations for several chiral two-plus three-nucleon interactions. Experimental values [198,197] are indicated as black bars. Data courtesy of J. D. Holt, J. Menéndez, and G. Hagen. [195]. As shown in Fig. 6, these findings have been reinforced by subsequent VS-IMSRG calculations, as well as the experimental observation of the first excited 2 + state in this nucleus [197].
Since this initial application in medium-mass nuclei, the EMλ/Λ family has seen widespread use in ab initio calculations due to its empirical quality, although the Hamiltonian's theoretical uncertainties are less well defined than for interactions that obey the chiral power counting more rigorously. Indeed the EM1.8/2.0 interaction was used in VS-IMSRG calculations to produce what is to my knowledge the first attempt at producing an ab initio mass table for nuclei up to the iron isotopes [186]. For selected nuclei up to the tin region, it also yields converged energies for ground and low-lying states that are in good agreement with experimental data [202,203]. It also yields slightly larger radii than previous interactions, although the underprediction is not eliminated entirely (see Refs. [194,202] and Section 3.2.3).
Multiple applications of the EMλ/Λ Hamiltonians in support of spectroscopy experiments have been published in recent years (see, e.g., [204,196,205,206,207]), and additional studies are underway, including an effort to better understand what makes the EM1.8/2.0 Hamiltonian so successful. Furthermore, a new generation of chiral NN+3N interactions is now available for applications in medium-mass and heavy nuclei [46,208,209,210,114].

Towards the Drip Lines
Neutron-rich nuclei are excellent laboratories for disentangling the interplay of nuclear interactions, many-body correlations and the continuum. Thus, data from the experimental push towards the drip line can offer important constraints for the refinement of chiral interactions if the many-body truncations and continuum effects are under control.
In practice, ab initio results for observables like the absolute energies of states still exhibit significant scale and scheme dependence due to truncations that are made in the EFT, the potential implementation of SRG evolutions, and the many-body methods. Since such variations tend to be systematic within families of interactions (and sometimes even across multiple families), differential quantities like separation and excitation energies or transition matrix elements often exhibit a weaker scale and scheme dependence -note, for example, the small systematic variation of the first excited 2 + states of the neutron-rich nickel isotopes for EMλ/Λ interactions. This makes energy differences an ideal observable for confronting ab initio results with experimental data.
Let us consider two-neutron separation energies as a concrete example. Sudden drops in these observables are a signal of (sub)shell closures (albeit not universally [194]) and in the neutron-rich domain, they are important indicators for the proximity of the drip line. Figure 7 shows MR-IMSRG ground-state and two-neutron separation energies of the calcium isotopes, obtained with the NN+3N(400) interaction used in Fig. 5, as well as the NNLO sat and EM1.8/2.0 interactions briefly discussed in the previous section. We note the overbinding produced by NN+3N(400) and the baffling accuracy of the EM1.8/2.0 results, given the approximations that went into the construction of this force, as well as the MR-IMSRG truncation. Common to all three interactions is the emergence of a very flat trend in the ground-state and separation energies in neutron-rich calcium isotopes, which will likely be further enhanced by the inclusion of continuum effects, and extended beyond the shown mass range. Similar flat trends emerge in many isotopic chains, as shown both in ab initio surveys based on chiral interactions [114,10,186] as well as a sophisticated Bayesian analysis of empirical EDF models [212]. Naturally, this will make the precise determination of the neutron drip line in the medium-mass region a challenging task, but also suggests that interesting features like alternating patterns of unbound odd nuclei and weakly-bound even nuclei with multi-neutron halos could emerge. This is an exciting prospect for the experimental programs at rare-isotope facilities.
With the exception of the NCSMC and HORSE methods discussed in Section 2.3.6, the inclusion of continuum degrees in configuration-space techniques has been focused on the use of the Berggren basis [151]. While such calculations are challenging due to the significantly increased single-particle basis size and the difficulties of handling the resulting complex symmetric Hamiltonians, applications in CC (see Refs. [12,213] and references therein), both valence and No-Core Gamow Shell Model [153,154,155,214,215] and IMSRG [216] calculations have been published. Common to all these approaches is that a configuration space interaction that is given in terms of SHO matrix elements is expanded on a basis containing SHO and Berggren states, hence it is still an open question how a direct implementation of the interactions in a basis with continuum degrees of freedom might modify  existing results. It is worth noting that such a construction has been achieved for phenomenological GSM interactions that have been tuned for light nuclei [217,218,219,220,221,222].
In light nuclei, the NCSMC has been applied with impressive success to describe a variety of exotic nuclei with up to three-cluster structures. For example, Calci et al. [223] carried out NCSMC calculations for 11 Be with several chiral NN+3N interactions to investigate the parity inversion of the ground and first-excited states in this nucleus from first principles. The authors found that the coupling between the NCSM and RGM sectors of the generalized eigenvalue has strong effects, but that among the tested interactions, only NNLO sat can produce the experimentally observed ordering of the states (see Fig. 8). However, it still underpredicts the splitting of these levels and as a result, overestimates the cross section for the photodisintegration 11 Be(γ, n) 10 Be. Additional applications of the NCSMC for exotic nuclei can be found in the review [130] and references therein, as well as the more recent works [224,225,226].

Accessing More Observables
The capabilities of ab initio approaches have also significantly expanded when it comes to the evaluation of observables other than the energies.
Nuclear Radii. Figure 9 shows MR-IMSRG results for the charge radii of calcium isotopes. The left panel illustrates the reasonable reproduction of the 40 Ca and 48 Ca charge radii that can be obtained for NNLO sat . The MR-IMSRG(2) results are slightly smaller than the experimental data due to differences in the truncations from the CCSD charge radius calculations that were used in the NNLO sat optimization protocol [193]. Note the steep increase in the experimental charge radii beyond 48 Ca: At the time of the measurement, NNLO sat was the only chiral NN+3N interaction exhibiting this feature, although other more recent interactions can replicate this trend as well [114,10]. Also note that none of the calculations are able to reproduce the inverted arc of the charge radii between 40 Ca and 48 Ca. In a CI picture, it is caused by strong mixing with 4p4h excitations into the pf -shell [227]. Since the MR-IMSRG(2) calculations shown here included only up to (generalized) 2p2h excitations and used particle-number projected Hartree-Fock Bogoliubov vacua as reference states that do not contain collective correlations (cf. Section 2.3.3), it is not surprising that the inverted arc cannot be reproduced. We will return to this issue of missing collectivity later.
While the EMλ/Λ interactions underpredict the absolute charge radii, they fare quite well in the description of radius differences, as suggested in the previous section. Figure 9(b) is adapted from a recent study that suggests a correlation between the charge radius difference of mirror nuclei, ∆R ch , and the slope of the symmetry energy in the nuclear matter equation of state [229]. We see that the MR-IMSRG results for ∆R ch are actually compatible with results from a multitude of Skyrme EDFs, and the value for the magic EM1.8/2.0 interaction falls into the uncertainty band of the experimental result.
Electromagnetic Transitions. Since the second half of the past decade, ab initio calculations for transitions in medium-mass nuclei have become more frequent, owing to the appropriate extensions of the IMSRG, CC and SCGF methods [230,204,231]. While results for transitions that are dominated by a few nucleons, e.g., M 1 transitions [230] or β decays (see Ref. [232] and the discussion below) can be quite good, the description of collective transitions is hampered by inherent truncations of these many-body methods, which are better suited for dynamical, few-particle correlations (see Sections 2. 3.3 and 2.3.4). Results from the SA-NCSM [139,140] and the IM-GCM discussed in Section 2.3.3 show that the modern chiral interactions themselves adequately support the emergence of nuclear collectivity.
Consider for example Fig. 10, which shows VS-IMSRG(2) results for the quadrupole transition from the first excited 2 + state to the ground state in 14 C, 22 O and 32 S [230]. The picture is fairly consistent for all four chiral NN+3N interactions that were used in the study: The 2 + energies are described quite well, but energies are not very sensitive to the details of the nuclear wave functions. In 14 C, the E2 transition is weakly collective, so the E2 matrix element is reasonably reproduced, while the matrix element for the collective transition in 32 S is underpredicted by 25-50%. The NN+3N(400) interaction gives a particularly poor result, but this is also related to the significant underestimation of the point-proton radius we obtain for this Hamiltonian, as discussed earlier.
The result for 22 O deserves special attention. The E2 transition matrix element is only a third of the experimental value, although the transition is only weakly collective. However, 22 O only has neutrons in an sd valence space, so the E2 matrix element would vanish in a conventional Shell Model calculation unless the neutrons have an effective charge. Such effective charges must be introduced by hand and fit to data in phenomenological Shell Model calculations. Here, we see that the VS-IMSRG decoupling naturally induces a non-vanishing quadrupole moment through an effective neutron charge in the one-body transition operator as well as an induced two-body contribution (see Ref. [230], and Ref. [231] for an analogous effort in SCGF theory). It is likely that the E2 strength could be improved by performing the VS-IMSRG calculation in a psd valence space, so that the proton dynamics is treated explicitly instead of implicitly by valence-space decoupling. Until recently, we were unable to perform such a multi-shell decoupling because of the IMSRG version of the intruder-state problem, but a promising workaround was introduced in Ref. [28].
Gamow-Teller Transitions. In recent years, there have been concerted efforts to understand the mechanisms behind the empirically observed quenching of Gamow-Teller (GT) transitions in medium-mass nuclei, in part due to its relevance to neutrinoless double-beta decay searches (see below). In Ref. [232], the authors show that this issue is largely resolved by properly accounting for the scale and scheme dependence of configuration-space calculations. By dialing the resolution scale to typical values favored by approaches like NCSM, CC and VS-IMSRG, correlations are shifted from the wave functions into induced two-and higher-body contributions to the renormalized transition operator, just as in the quadrupole case discussed above.
The transition operator, including two-body currents, is consistently evolved to lower resolution scale alongside the nuclear interactions, keeping the induced contributions. The transition matrix  elements of the evolved operator are then computed with the NCSM in light nuclei, and VS-IMSRG in sd-and pf -shell nuclei, leading to agreement with experimental GT strengths within a few %. In contrast, the bare GT operator must be quenched by 20-25% via the introduction of an effective axial coupling, g eff A < g A , to yield agreement with experimental beta decay rates. The GT transitions in light nuclei have also been evaluated in the GFMC, most recently with consistently constructed local chiral interactions and currents [234,235]. Interestingly, the inclusion of two-body currents seems to consistently enhance the GT matrix elements, while it tends to quench the matrix element in NCSM calculations. Since this is almost certainly related to the differences in the resolution scale and calculation scheme, the disentanglement of these observables might yield further insights into the interplay of wave function correlations and the renormalization of the transition operators.
Neutrinoless Double Beta Decay. Due to the high impact the observation of neutrinoless double beta decay (or lack thereof) would have on particle physics and cosmology, the computation of nuclear matrix elements (NMEs) for neutrinoless double beta decay is a high priority for nuclear structure theory. Precise knowledge of the NMEs for various candidate nuclei is required to extract key observables like the absolute neutrino mass scale from the measured lifetimes (or at least, any new bounds that would be provided by experiment). Most calculations of the NME to date were subject to the lack of comparability between phenomenological nuclear structure results that was discussed in Section 3.1, hence a new generation of ab initio calculations with quantified uncertainties is required.
A major step in that direction was the first calculation of the NME for the decay 48 Ca → 48 Ti based on chiral interactions [29]. The IM-GCM approach discussed in Section 2.3.3 was used to describe the structure of the intrinsically deformed daughter nucleus 48 Ti, achieving a satisfactory reproduction of the low-lying states and their quadrupole transitions (see Fig. 11). Since the initial publication (blue spectra in Fig. 11(a)), the description of the excited states has been improved further through the admixing of cranked configurations (red spectra), without affecting the NME Hergert A Guided Tour of Ab Initio Nuclear Many-Body Theory ( Fig. 11(b)). Work on quantifying the uncertainties due to the many-body method, the Hamiltonian, and the transition operator is underway, in preparation for the computation of the NMEs of more realistic candidate nuclei like 76 Ge and 136 Xe.

Response and Scattering
From the computation of transitions between low-lying levels, it is only a small step to the computation of nuclear response functions and cross sections, although the implementation can be challenging and the applications are often computationally expensive.
Nuclear Response Functions. In light nuclei, GFMC is a powerful yet numerically heavy tool for computing exact nuclear response functions (see, e.g., Refs. [236,237]). In medium-mass nuclei, applications of SCGF and CC techniques to the computation of the nuclear response have been published in recent years. As mentioned in Section 2.3.5, the Green's functions computed in the standard or Gor'kov ADC Green's function schemes inherently contain information about the nuclear response that has been used to study both electromagnetic and weak processes of medium-mass nuclei [120,238,239,119,240].
In the Coupled Cluster framework, response functions have been computed by merging CC with up to Triples excitations with the Lorentz Integral Transformation (LIT) technique [241,242,243,244,245]. Immediately after its inception, this approach was used to for the first ab initio calculations of dipole response and the related photodisassociation cross section of medium-mass closed-shell nuclei [241,242]. More recently, it was used to compute the electric dipole polarizability α D of nuclei like 48 Ca [243,246,244] and 68 Ni [247]. Together with measurements of the charge radius, this quantity can be used to constrain ab initio calculations that will in turn allow the theoretical extraction of the neutron point radius as well as the thickness of the neutron skin.
An important application for nuclear response calculations is to map out the neutrino response of 40 Ar, the primary target material in detectors for the short-baseline [248] and long-baseline neutrino experiments, like the Deep Underground Neutrino Experiment (DUNE) Far Detector [249,250]. At low energies, the cross section for coherent neutrino elastic scattering is essentially determined by the weak form factor of 40 Ar, which has recently been computed using CC techniques [251]. This work is complementary to SCGF calculations of the neutrino response in the region of the quasi-elastic peak by Barbieri et al. [238].
Nuclear Reactions. As discussed in Section 2.3, there has been enormous progress in the development of unified treatments of ab initio nuclear structure and reactions. Here, I want to highlight two among a bevy of impressive results. Figure 12(a) shows S− and D−wave phase shifts for α − α scattering, computed order by order in Lattice EFT [179,180]. These calculations are made possible by the lattice's capability to describe clustered states (also see Refs. [176,177,178]), as well the development of the APM and associated algorithms. The results for the phase shifts show the desired order-by-order improvement, and the inclusion of higher-order terms of the chiral expansion is expected to improve agreement with experimental data. The near identical NLO and NNLO phase shifts in the S−wave appear to be the result of an accidental cancellation that is not occurring in the D−wave phase shifts.
In Ref. [252], the authors studied deuterium-tritium (D-T) fusion using the NCSMC. One of the main results of this work is shown in Fig. 12(b), which compares the NCSMC D-T reaction rates for polarized and unpolarized fuels to each other, as well as rates obtained with several widely used part of the nuclear data required by these applications. A xample is the fusion of deuterium (D) with tritium ( 3 H or nerate a 4 He nucleus (α-particle), a neutron, and 17.6 energy released in the form of kinetic energy of the s. This reaction, used at facilities such as ITER 1 and NIF 2 ursuit of sustained fusion energy production, is charby a pronounced resonance at the center-of-mass (c.m.) f 65 keV above the free D and T nuclei due to the forof the J π = 3/2 + resonance of the unbound 5 He nucleus. ars ago, it was estimated 3 that, in the ideal scenario in e spins of the reactants are perfectly aligned in a totalconfiguration and assuming that the reaction is isotropic, ld achieve an enhancement of the cross section by a factor .5, thus improving the economics of fusion energy gen-. However, while the unpolarized cross section and some g-power data exist, no correlation coefficients have been d yet to confirm this prediction 5 . More generally, what known about the properties of the polarized DT fusion rred from measurements of the D 3 He reaction 6 . DT fusion is a primary example of a thermonuclear in which the conversion of two lighter elements to a one occurs through the transfer of a nucleon from the e (D) to the target (T). Despite the fairly small number of s involved in this process, arriving at a comprehensive nding-in terms of the laws of quantum mechanics and erlying theory of the strong force (quantum chromody--of the interweaving of nuclear shell structure and dynamics giving rise to the DT fusion already represents able challenge for nuclear theory. rds this goal, a pioneering ab initio study of the DT fusion ied out in ref. 7 , using a nucleon-nucleon (NN) interact accurately describes two-nucleon data and representing function on a basis of continuous "microscopic-cluster" ade of D+T and n+ 4 He pairs in relative motion with o each other. However, this approach was unable to yield f adequate fidelity, due to the omission of the three-(3N) force-disregarded for technical reasons. Numerous ave shown that this component of the nuclear interacssential for the reproduction of single-particle propermasses [13][14][15] , and spin properties 10,16 , all impactful in the case. Besides the 3N force, the approach of ref. 7 also complete treatment of short-range five-nucleon correwhich are crucial to arrive at the accurate description of resonance. The formation of this rather long-lived ce as a correlated, localized system of five nucleons built g the fusion process is integral to the reaction mechanally, no polarization observables were calculated in the ref. 7 . e following, we report on ab initio predictions for the d DT fusion using validated NN and 3N forces derived in ework of chiral effective field theory (EFT) 17,18 , a powl that enables the organization of the interactions among and neutrons in a systematically improvable expansion o the fundamental theory of quantum chromodynamics. ntum-mechanical five-nucleon problem is solved using ore shell model with continuum (NCSMC) 10,19 , where the model space includes D+T and n+ 4 He microscopic-cluster states, plus conventional static solutions for the aggregate 5 He system 20 . This enables a fully integrated description of the reaction in the incoming (outgoing) channel, where the reactants (products) are far apart, as well as when all five nucleons are close together. We show that this approach yields an accurate reproduction of the DT cross section for unpolarized reactants, discriminating among reaction rates from phenomenological evaluations and demonstrating in detail the small contribution of anisotropies in the vicinity of the 3/2 + resonance. The maximum enhancement of the polarized cross section varies as a function of the deuterium incident energy, dropping significantly above 0.8 MeV. However, such variation is slow in the narrow range of optimal energies for the reaction, resulting in a rather constant enhancement of the rate of fusion, compatible with the historic approximate estimate.

Results
Validation of model for unpolarized reaction observables. We begin our study with a validation of our ab initio reaction method on existing experimental data for the unpolarized DT  36 , which is intended for applications in astrophysics simulations. Overall, we find that they agree well even at energies above the resonance. In more detail, our calculation agrees best with the phenomenological R-matrix evaluation, particularly at higher energies where data are typically scarcer. In our case, the uncertainties due to the finiteness of the model space are indisti initio Supple analys with t stemm tainty prohib enolog thresh practic to ma peak. albeit with This enhan the en energi the M probab nifican polariz the un that i rates p illustra of a st its ope compr fuel in in ref.
Angul the de small a observ they p Computed maximum enhancement factor (over all possible values of p z q z and p zz ) of energy. The maximum enhancement is always found for p z q z ; p zz ¼ 1. The "NCSMCincluding '≠0 partial waves after a phenomenological correction of −5 keV to the posit enhancement factor obtained without such phenomenological correction (that is, the Also shown is the maximum enhancement factor obtained by retaining only the ' ¼ wave Comparison between the computed DT reaction rate (N A σν h i, with N A the Avogadro number) for unpolarized and polarized fuel with aligned spins as a function of the temperature, T. The "Polarized" and "Unpolarized" labels stand for the present results obtained with the phenomenological correction of −5 keV to the position of the 3/2 + resonance (dubbed NCSMC-pheno). We use reactants' polarization parameters achievable in the laboratory, that is p z ; p zz ¼ 0:8 and q z ¼ 0:8. Also shown for comparison are the unpolarized reaction rates obtained from the widely adopted parametrization of the DT fusion cross section of Bosch and Hale 34 (labeled as "Bosh and Hale"), from the R-matrix fit of Descouvemont 35 (labeled as "Descouvemont") and from the NACRE compilation 36  parameterizations of the D-T fusion cross section. The NCSMC calculations indicate that for an experimentally realizable polarized fuel with aligned spins, a reaction rate of the same magnitude as for unpolarized fuel can be achieved at about half the temperature. Naturally, this suggests that polarized D-T fuels will allow a more efficient power generation in thermonuclear reactors.

Emergence of Empirical Nuclear Structure Models from Ab Initio Calculations
The progress in ab initio calculations over the past decade has not only led to impressive results for nuclear observables, but also revealed the long-surmised underpinnings of empirical models of nuclear structure. In many cases, the ideas that led to the formulation of such models were shown to be correct, but they could not be verified at the time because RG and EFT techniques or sufficient computing power for a more thorough exploration were not available.  The Nuclear Shell Model. The first prominent example I want to discuss is the nuclear Shell Model and some of the "folklore" surrounding it. We can immediately make the observation that the Shell model picture is inherently a low-momentum description of nuclear structure. It is based on the assumption that nucleons are able to move (almost) independently in a mean field potential, and that nuclear spectra can be explained by the mixing of a few valence configurations above an inert core via the residual interaction. As we know now, the existence of a bound mean-field solution and a weak, possibly perturbative residual interaction relies on the decoupling of low and high momenta in the nuclear Hamiltonian [1,6,253], e.g., by an SRG transformation. Historical approaches to exploit this connection to construct the Shell model from realistic nuclear forces [254,255,256] failed in part because the decoupling of the momentum scales via Brueckner's G−matrix formalism [65,66,67] was not as good as believed [1].
In addition to the momentum-space decoupling, one must also decouple the valence space configurations from the excluded space. This can be achieved using a variety of techniques (cf. Sections 2.3.3-2.3.6), and either by performing transformations in sequence, or designing a single procedure that achieves both types of decoupling simultaneously. In practice, the former strategy tends to be more efficient and less prone to truncation errors -an example is the VS-IMSRG decoupling of Hamiltonians that have been evolved to a low resolution scale by means of a prior SRG evolution (see Sections 2.2 and 2.3.3, as well as Ref. [76]). An added benefit of using low-momentum interactions is that the Shell Model wave functions will qualitatively resemble those obtained by a no-core method using the same Hamiltonian without valence decoupling. This facilitates qualitative comparisons and allow us to apply the same intuitive picture. For quantitative comparisons, the effects of all unitary transformations must be carefully taken into account [257]. Figure 13 illustrates the effect of the discussed transformations via the deviations between the computed and experimental energies of close to 400 levels in the sd-shell. Since the EM1.8/2.0 interaction used in these calculations has a low resolution scale, simply using the valence-space matrix elements of the input Hamiltonian without any further valence-space decoupling yields a root-mean-square (rms) deviation of "only" about 1.7 MeV, which is not outright disastrous. When we apply the VS-IMRSG to decouple the valence space, the newly evolved interaction yields a much Hergert A Guided Tour of Ab Initio Nuclear Many-Body Theory improved rms deviation of approximately 650 keV, which is better than for some of the older sd-shell interactions, albeit not as good as the USDB Hamiltonian, which is shown for comparison [258,259]. This is not really surprising: USDB essentially represents the best possible fit to experimental data under the model assumptions, i.e., the choice of a pure sd-shell valence space, the restriction to a two-body Hamiltonian, the omission of isospin-breaking effects from the Coulomb interaction and the nuclear interactions, and the empirical A-dependence multiplying the two-body matrix elements (TBMEs). The accuracy of the VS-IMSRG results, on the other hand, is affected by possible deficiencies in the input Hamiltonian and the use of the VS-IMSRG(2) truncation. Naturally, both of these aspects will be improved systematically in future calculations.
Phenomenological adjustments of effective Shell Model interactions like the A-dependent scaling factors in the USD Hamiltonians or Zuker's monopole shift [260] are typically attributed to the changes in the nuclear mean field away from the core, as well as missing three-body interactions. In Ref. [76], the VS-IMSRG is used to demonstrate that this is indeed the case. As described in Section 2.3.3, upon normal ordering, the three-body force gives contributions to operators of equal and lower particle rank, which in the Shell Model case amounts to the core energy, single-particle energy, and two-body matrix elements. All of these contributions become A-dependent in the VS-IMSRG, but one can shift the A-dependent parts completely into the TBMEs, like in phenomenological interactions, without changing the Hamiltonian matrix in the many-body Hilbert space or its eigenvalues.
Procedures like the VS-IMSRG decoupling also let us track in detail how operators besides the nuclear interactions evolve when they are subject to the valence-decoupling transformation. Recall from the discussion in Section 3.2.3 that this can even quantitatively explain the quenching of the Gamow-Teller strength in phenomenological Shell Model calculations, provided two-body current contributions to the initial transition operator are taken into account as well. For electromagnetic transitions, the renormalization of the one-body transition operator and the appearance of induced terms generate at least some part of the usual phenomenological effective charges, but a more complete treatment of nuclear collectivity (cf. Section 2.3.3) as well the inclusion of current contributions to these operators are developments that need to be undertaken in the coming years.
Emergence of Collectivity. Both NCCI and VS-IMSRG calculations with chiral NN+3N interactions have demonstrated that these interactions do indeed produce the telltale features of collective behavior in nuclear spectra [26,213,262,263,141]. Upon a bit of reflection, it is not surprising that reasonable results on rotational bands, for instance, should be found in these approaches: While they rely on particle-hole type expansions, the exact diagonalization is done in a complete model space of up to A v hA v p excitations, where A v is the number of valence nucleons. In contrast, euation-of-motion methods that typically employ 1p1h or 2p2h truncations struggle with the description of collectivity in low-lying states [122,83,203], but they do work reasonably well for giant resonances [241,242].
As argued in Sections 2.3.3 and 2.3.6, bases built on particle-hole type expansions are not ideally suited to the description of collective correlations. The SA-NCSM [139] instead uses irreducible representations of SU(3) or Sp(3,R), the dynamical symmetry groups of collective models [264], to achieve a much more efficient description of collective behavior in nuclei. This is illustrated for the case of 20 Ne in Fig. 14. The SA-NCSM calculations [140] based on the two-nucleon NNLO opt potential [261] describe the ground-state rotational band extremely well, all the way to the J = 8 +

Rethinking the Many-Body Expansion
A substantial part of the appeal of methods like CC, IMSRG and SCGF is their polynomial scaling. For the purposes of uncertainty quantification (UQ), we need to be able to evaluate at least two consecutive truncation levels to assess the convergence of the many-body expansion in nuclei for which exact calculations are not feasible. Efforts in that direction have been in progress for some time, and while some methods are at a more advanced stage than others, the improved truncations should be available for regular use within the next couple of years [12,100,244,116,86,265,10]. In part, this is owing to the development of computer tools that automate tasks like diagrammatic evaluation or angular momentum coupling [266,267]. The computational scaling of these approaches will be of order O(N 8 ) or O(N 9 ), which makes applications a task for leadership-class computing resources for the foreseeable future. It is clear that it will not be feasible to just push the calculations further, since we would then face a (naive) O(N 12 ) scaling.
Applications where we would expect to need high-order truncations involve nuclear states with strong collective correlations, provided we work from a spherical reference state. As explained Hergert A Guided Tour of Ab Initio Nuclear Many-Body Theory in Section 2.3, this issue can likely be addressed either by using mean-field reference states with spontaneously broken symmetries (cf. Section 2.3.4) or using correlated reference states in the first place (cf. Section 2.3.3), and the first applications of the IM-GCM give credence to that idea. Moreover, there is first evidence that the CC and IMSRG truncations converge much more rapidly for observables that are sensitive to collectivity [268], i.e., the current state-of-the-art truncations may be sufficiently precise.
The IMSRG framework also offers perspectives for the construction of further IMSRG hybrid methods (cf. Section 2.3.3). Based on the successes of both the IM-NCSM and IM-GCM it would be worthwhile to use IMSRG-evolved Hamiltonians in the SA-NCSM or techniques like the Density Matrix Renormalization Group, which is also capable of efficiently describing strong collective correlations under certain conditions [269,270].

Leveraging Computational and Algorithmic Advances
The progress in ab initio many-body calculations is not simply due to the availability of increasingly powerful computational resources, but also due to dedicated collaborations with computer scientists to ensure that the available high-performance computers are used efficiently. Such collaborations will only grow more important as hardware architectures change rapidly and a growing demand for computing time requires users to demonstrate sufficient efficiency to be granted access to supercomputers.
Measures to boost the numerical efficiency can also be taken at the many-body theory level. Efficient calculations rely on finding optimal representations of the relevant physical information that is encoded in the Hamiltonian. Algorithmic gains are possible whenever there is a mismatch, either because we made convenient choices, e.g., by expanding many-body states in terms of simple Slater determinants, or because we were not able to recognize simplifications beforehand, e.g., due to hidden or dynamical symmetries.
The SRG has played a key role in addressing the first points at the level of the nuclear interaction over the past two decades, and SRG and IMSRG can be applied in novel ways to explore dynamical symmetries [55]. In the construction of a configuration space, the selection of the single-particle basis leaves room for optimization. Indeed, the natural orbitals introduced in Ref. [131] lead to faster model-space convergence in NCSM and CC calculations, implying a more compact Hamiltonian matrix in natural orbital representation. The efficiency of this representation can be leveraged further by making robust importance truncations based on analytical measures, e.g., in MBPT, CC, or IT-NCSM [271,9].
The aforementioned steps make use of prior theoretical knowledge, e.g., to identify desired decoupling patterns in interactions, or define analytical measures for the importance of basis states. If such knowledge is not available, or we want to avoid bias, we can leverage a myriad of Principal Component Analysis (PCA) methods to factorize interactions or intermediate quantities in manybody calculations [272,271]. This can potentially even give us control over the computational scaling of nuclear many-body methods (see, e.g., [273,274,275,276,277]).
A very noteworthy development with origins in nuclear physics is Eigenvector Continuation (EVC) [279,280], a method for learning manifolds of eigenvector trajectories of parameter-dependent Hamiltonians. The method has been employed in several contexts, e.g., to stabilize high-order MBPT expansions [81] and to construct emulators for nuclear few-and many-body calculations   [278]. Figure reprinted with permission from the American Physical Society. [281,278]. As an example, Fig. 15 shows a global sensitivity analysis of CC results for 16 O under variations of the chiral LECs [278]. Eigenvector continuation was used to learn representations of the CCSD Hamiltonian and charge radius operators in a 64-dimensional subspace of the space of CCSD ground-state wave functions for interactions with 16 varying LECs. The subspace-projected Hamiltonian was then sampled more than a million times on laptop, while full CCSD calculations of the same ensemble would be completely unfeasible. The successful applications of EVC suggest that the method should be further explored as a tool for improvement, emulation and UQ in other many-body methods in the (near) future.

Uncertainty Quantification
In typical nuclear many-body calculations as discussed in Secs. 2 and 3 the main sources of theoretical uncertainties are the EFT truncation of the observables and the many-body wave function, either due to many-body expansion and/or model space truncations in configuration space approaches, lattice discretization effects in Lattice EFT, or the specific form of the wave function ansatz in QMC. If an SRG evolution is applied, there is an additional uncertainty associated with the truncation of induced operators (see Section 2.2).
The application of Bayesian methods has led to enormous progress in the quantification of the EFT uncertainties [282,283,284,34,36,35,285], and it would be highly desirable to apply the same approach to the many-body uncertainties as well. The most challenging amongst these are the truncation of the many-body expansion in methods like CC, IMSRG or SCGF, and the truncation of the free-space SRG flow of observables. In contrast, the infrared effects of finite-basis size truncations Hergert A Guided Tour of Ab Initio Nuclear Many-Body Theory in HO bases -or general orbitals that are at some point expanded in an HO basis -are well understood for the energy and other observables, and they can be systematically corrected for [286,287,288,289,290]. The situation is less clear for ultraviolet basis-size errors [291], but this error can be suppressed by working at appropriate values of the HO frequency.
An ideal uncertainty analysis would combine the exploration of EFT and many-body uncertainties for nuclear observables of interest, using EC or Machine Learning (ML) to construct emulators that allow an efficient sampling of the parameter space. In such an effort, the generation of sufficient training data poses a significant challenge, because it would require calculations at several truncation levels (see Section 4.1). A possible strategy for mitigating this issue is to combine non-perturbative methods with cheaper high-order MBPT in Bayesian mixed models (see Refs. [292,212,293] for applications in nuclear physics). The successful application of factorization methods to the nuclear many-body problem could likely resolve the issue once and for all by reducing the computational scaling of high-order truncations, at the cost of introducing an additional uncertainty from the factorization procedure.
On the road towards the destination represented by such a "complete" UQ framework, the intermediate milestones will already provide valuable insights into open issues in the EFTs of the strong interactions, and enable the design of better protocols for constraining and refining EFT-based interactions and operators (see, e.g., Refs. [294,295] and references therein.)

Strengthening and Employing the Hierarchy of Strong Interaction EFTs
Strong interaction physics is a multi-scale problem, and there are good reasons for making better use of the hierarchy of Effective (Field) Theories at our disposal. At the top level, we have QCD, followed by EFTs involving hyperons that can be eliminated progressively until we arrive at the "traditional" pionful and pionless chiral EFTs (see Refs. [296,297] and references therein). At even lower scales, one can formulate an EFT for nuclear halos (or clusters) [297] and make the connection to nuclear DFT and collective models, which can be understood as EFTs as well [298,299,300,301,302,303,304,305].
At least in principle, the different levels of this hierarchy can be connected either by computing observables with different theories and matching the LECs, or using RG flows to track in detail how the theories evolve from one into another. While matching procedures have been applied successfully to modern EFts in nuclear physics [306,307,308,309,310,311] as well as efforts to match more traditional models of nuclear structure to ab initio calculations [312,313,314], making the connection through RG methods is a more daunting task. While I must admit to bias in this regard, I still consider this an effort worth undertaking. The success of SRG techniques in nuclear physics demonstrate how these methods reveal the most effective degrees of freedom even in situations were the separation of scales is not perfectly clear. Moreover, RGs would also reveal unexpected features of the power counting schemes, like the enhancement or inadvertent omission of certain operators (see Ref. [51] and references therein).
Tackling Power Counting Issues. Throughout this work, I have alluded to shortcomings and issues of the current generation of chiral interactions, like the struggle to achieve a good simultaneous description of nuclear binding energies and radii (see Section 3.1). Recent efforts to construct new, accurate chiral interactions have revealed that this issue is connected to the use of local or nonlocal regulators, with the latter being favored for better descriptions [208,114]. In another exploration of nonlocally regularized chiral forces [209,210], a tension between the simultaneous description of Hergert A Guided Tour of Ab Initio Nuclear Many-Body Theory nuclear matter and finite was observed in the attempts to fit the chiral LECs. In QMC calculations, it was demonstrated that the use of local regulators breaks the equivalence of parameterizations of the interaction that are connected by Fierz identities, in certain cases with disastrous consequences [188]. Meanwhile, Epelbaum and co-workers have proposed the use of a more nuanced semilocal regularization scheme that applies local regulators to the long-range pion exchange and nonlocal regulators to the short-range contact terms [46,3]. While physical arguments can be made in favor of different regularization schemes, perhaps especially the semilocal one, the significant scheme dependence is at odds with the principles of EFT, which would predict regulator artifacts to be pushed beyond the order at which one currently works.
It has also been suggested that the scales of the chiral EFT interaction and the inherent scales of the many-body configuration space (e.g., IR and UV cutoffs inherited from a HO basis, see Section 4.3) or coordinate space wave functions should not be treated independently, and that by doing so, current many-body approaches might at least contribute to power counting issues. There have been a few efforts to explore this problem, but more work is clearly required [315,316,317,318,319,320].
Application Needs. Aside from the formal need to make progress on the power counting issue, there are also concrete application needs that call for a tighter coupling between QCD and the nuclear EFTs. For example, the chiral EFT transition operator for neutrinoless double-beta decay (see Section 3.2.3) contains counter terms whose LECs can only be determined from Lattice QCD [321,322,323,324].
The dawning of a new age in our understanding of neutron stars, heralded by the detection of gravitational waves from the neutron-star merger GW170817, has taken the demand for accurate neutron and nuclear matter equations of state to a new level (see, e.g., Ref. [159] and references therein). While ab initio calculations of infinite matter up to the saturation region based on chiral interactions are reasonably well controlled [285,325,190,159], the supranuclear densities probed in merger events are beyond the range of validity of regular pionful chiral EFT. To increase its validity, hyperons must be taken into account as dynamical degrees of freedom (see [296] and references therein), and the entire set of nuclear and hyperon LECs must be readjusted at the increased breakdown scale. For the NN sector, this is unproblematic due to the plethora of available scattering data. Since no direct experiments on three-neutron or three-proton systems are feasible, the only available experimental constraints come from finite nuclei, which implies that the corresponding channels of the 3N interaction are only constrained at sub-saturation densities. The world database of hyperon-nucleon scattering data is also quite limited. Thus, a high-precision interaction for describing the equation of state at high density can only be constructed with the help of Lattice QCD constraints on the 3N and hypernucleon LECs.

Interfacing with Reaction Theory
The final important research direction for the coming decade I want to discuss here are efforts to interface the advanced ab initio nuclear structure methods at our disposal with reaction theory [326].
As discussed in Sections 2.3.6 and 3.2.4, the NCSMC has been applied with great success to the reactions of light nuclei at low energies, but its computational complexity makes this approach unfeasible for nuclei beyond A ≈ 10 − 20. Work has begun on a similar approach that combines SA-NCSM with the RGM, leveraging the efficiency of the symmetry-adapted basis to reach mediummass nuclei [327] (cf. Sections 2.3.6 and 3.3). Since the RGM is just a special case of a Generator

A Guided Tour of Ab Initio Nuclear Many-Body Theory
Coordinate Method, the IM-GCM discussed in Sections 2.3.3 and 3.2.3 is a promising candidate for extending this type of reaction calculations to even heavier nuclei.
Methods that are similar in spirit to these combinations of structure approaches witht the RGM are the APM, which can provide an interface between structure and scattering in Lattice EFT (cf. Sections 2.3.8 and 3.2.4), as well as the GSM Coupled Channel (GSM-CC) approach, which was developed to describe reactions between light projectiles and targets that are treated in the GSM with a core [328,221,329]. Thus far, applications of the GSM-CC have been based on phenomenological valence-space interactions, but new efforts are underway to directly construct suitable Hamiltonians based on EFT principles [220,330], or derive the effective interactions from chiral forces with the techniques discussed in Section 2.3 (see [214,215]). Of course, the GSM-CC ideas could also be applied to the No-Core GSM [153,155,218], although the computational complexity would limit such an approach to light nuclei.
A complementary strategy for bridging nuclear structure and reactions for medium-mass nuclei is the construction of optical potentials for use in traditional reaction calculations. In SCGF theory, the optical potential for elastic nucleon-nucleus scattering is given by the one-body self energy, which is obtained as a byproduct of a nuclear structure calculation, and can be used with little effort in reaction codes [125]. Similarly, Rotureau et al. constructed optical potentials by extracting the self-energy from the Coupled Cluster Green's Function [123,124,331]. One can roughly view this procedure as performing a GF calculation with the similarity-transformed CC Hamiltonian, which does not require self-consistent iterations because of the CC decoupling (cf. Section 2.3.4). Optical potentials can also be constructed by folding scattering T -matrices with ab initio density matrices. This technique was applied for NCSM density matrices by two collaborations in Refs. [332,333] and [334,335], respectively, and more applications are underway.
While the published results from the optical-potential based approaches are promising, an important aspect of these calculations must be checked carefully in the near term: The optical potential depends on the resolution scale of the used chiral interactions, and the calculation scheme, which encompasses the truncations in the operators and many-body method, as well as the choice of regulator in the interaction [336,257]. To produce scale-and scheme-independent observables, these choices must be matched by the reaction theory. Matching the resolution scales is probably the easier of the two checks, but it will require the analysis of free-space SRG transformations on the reaction theory side. Once structure and theory are defined at a matching resolution scale, any residual scheme dependence of the observables will give rise to the remaining theoretical uncertainty of the combined calculation.

EPILOGUE
Thus concludes our little excursion through the landscape of state-of-the-art ab initio nuclear many-body theory, but of course, the road goes ever on. I hope that this guided tour has contributed to your appreciation of the immense progress the community has made in the last ten years, as well as the challenges that we are facing on the next stage of the road. None of the obstacles in our path are unsurmountable, and while we chip away at them, results from ab initio calculations can make meaningful contributions to the analysis and planning of nuclear physics and fundamental symmetry experiments. With new facilities launching in the next couple of years, the fun will begin in earnest, so here's looking forward to the next decade! Hergert A Guided Tour of Ab Initio Nuclear Many-Body Theory ACKNOWLEDGMENTS This work has been shaped by an enormous amount of discussions over the past decade, and naming all discussion partners would likely require multiple pages -it is safe to say that if your work is cited here, we have likely talked in person at some point. I am deeply grateful for these conversations, and the spirit of the ab initio nuclear structure and reactions community that fosters such exchanges.
I would like to thank the current (and former) members of my research group as well as colleagues at NSCL/FRIB, who had the most immediate impact on this work by virtue of being a short walk away. Particular thanks go out to S. Bogner, K. Fossez, M. Hjorth-Jensen, R. Wirth, and J. M. Yao. Special thanks are owed to R. Stroberg and S. Elhatisari, who helped out with sudden requests for figures.
I am also grateful to the Institute of Nuclear Theory for its hospitality, which was the venue for many of the aforementioned discussions, most recently during the program INT-19-2a, "Nuclear Structure at the Crossroads".  Next-to-Next-to-Leading Order (EFT) N 3 LO Next-to-Next-to-Next-to-Leading Order (EFT) N 4 LO Next-to-Next-to-Next-to-Next-to-Leading Order (