A Products and Commutators of Normal-Ordered Operators

# In-Medium Similarity Renormalization Group for Closed and Open-Shell Nuclei

## Abstract

We present a pedagogical introduction to the In-Medium Similarity Renormalization Group (IMSRG) framework for ab initio calculations of nuclei. The IMSRG performs continuous unitary transformations of the nuclear many-body Hamiltonian in second-quantized form, which can be implemented with polynomial computational effort. Through suitably chosen generators, it is possible to extract eigenvalues of the Hamiltonian in a given nucleus, or drive the Hamiltonian matrix in configuration space to specific structures, e.g., band- or block-diagonal form.

Exploiting this flexibility, we describe two complementary approaches for the description of closed- and open-shell nuclei: The first is the Multireference IMSRG (MR-IMSRG), which is designed for the efficient calculation of nuclear ground-state properties. The second is the derivation of nonempirical valence-space interactions that can be used as input for nuclear Shell model (i.e., configuration interaction (CI)) calculations. This IMSRG+Shell model approach provides immediate access to excitation spectra, transitions, etc., but is limited in applicability by the factorial cost of the CI calculations.

We review applications of the MR-IMSRG and IMSRG+Shell model approaches to the calculation of ground-state properties for the oxygen, calcium, and nickel isotopic chains or the spectroscopy of nuclei in the lower shell, respectively, and present selected new results, e.g., for the ground- and excited state properties of neon isotopes.

1

## 1 Introduction

More than 60 years have passed since Rainwater, Bohr, and Mottelson published the seminal works that led to them winning the 1975 Nobel Prize in Physics [[1], [2], [3], [4], [5]]. The collective model developed in these publications is an essential tool for nuclear physicists, and it qualifies as one of the most successful data-driven approaches to nuclear structure. A variety of approaches exist in this category, ranging from local to global applicability, and from microscopic to macroscopic views of nuclei: Collective Bohr-Mottelson Hamiltonians are fine-tuned to specific nuclei or limited regions of the nuclear chart, and do not treat the dynamics of all nucleons on a fully microscopic level. The traditional nuclear configuration interaction (CI) approach (see, e.g., [[6], [7]]) uses phenomenological interactions that are highly optimized, e.g., to -shell data [[8]], and treats only the dynamics of valence nucleons on top of an inert core in fully microscopic fashion. Finally, nuclear Density Functional Theory (DFT) takes a global perspective, and aims for a microscopic description of the entire nuclear chart based on energy density functionals (EDFs) that are optimized to data [[9], [10], [11], [12]].

The philosophy behind data-driven models is complementary to that of ab initio nuclear many-body theory, although the lines are somewhat blurry. An ab initio approach attempts to describe nuclear structure and dynamics based on fundamental degrees of freedom and their interactions. In the Standard Model, the fundamental theory of strong interactions is Quantum Chromodynamics (QCD), but a description of nuclear observables on the level of quarks and gluons is not feasible, except for the lightest few-nucleon systems (see, e.g., [[13]]). Instead, we start from nuclear interactions that describe low-energy QCD observables in the and systems, like scattering data or binding energies. Nowadays, such interactions are derived in Chiral Effective Field Theory (EFT), which provides a constructive framework and organizational hierarchy for , , and higher many-nucleon forces, as well as consistent electroweak operators (see, e.g., [[14], [15], [16], [17], [18], [19], [20], [21], [22], [23], [24], [25]]). Since Chiral EFT is a low-momentum expansion, high-momentum (short-range) physics is not explicitly resolved by the theory, but parametrized by the so-called low-energy constants (LECs).

In principle, the LECs can be determined by matching calculations of the same observables in chiral EFT and (Lattice) QCD in the overlap region of the two theories. Since such a calculation is currently not feasible, they are fit to experimental data, typically in the , , and sectors. Recently, Ekström et al. have developed an optimization protocol for chiral interactions that gives up on the reductionist approach of fixing the LECs in the few-nucleon system, and includes certain many-body data in the fit as well. The many-body data, e.g., selected radii, are chosen in order to improve the deficient saturation behavior of chiral interactions that are used as input for nuclear many-body calculations. The first interaction optimized with this protocol is [[26]], which is able to accurately describe the ground-state energies and charge radii of at the same time. Following the same philosophy but not the same approach, Shirokov et al. have produced Daejeon16, a softened chiral interaction that has been tuned for the description of light nuclei without explicit forces [[27], [28], [29], [30]].

Renormalization group (RG) methods are natural companions for EFTs, because they make it possible to smoothly connect theories with different resolution scales and degrees of freedom. Since they were introduced in low-energy nuclear physics around the start of the millennium [[31], [32], [33], [34]], they have provided a systematic framework for formalizing many ideas on the renormalization of nuclear interactions and many-body effects that had been discussed in the nuclear structure community since the 1950s. For instance, soft and hard-core interactions can reproduce scattering data equally well, but have significantly different saturation properties, which caused the community to move away from the former in the 1970s (see, e.g., [[35]]). What was missing at that time was the recognition of the intricate link between the off-shell interaction and forces that was formally demonstrated for the first time by Polyzou and Glöckle in 1990 [[36]]. From the modern RG perspective, soft- and hard-core interactions emerge as representations of low-energy QCD at different resolution scales, and the dialing of the resolution scale necessarily leads to induced forces, in such a way that observables (including saturation properties) remain invariant under the RG flow (see section 2 and [[33], [34]]). In conjunction, chiral EFT and nuclear RG applications demonstrate that one cannot treat the sectors in isolation from each other.

During the 1960s, Kuo and Brown pioneered work on the ab initio derivation of effective interactions for nuclear valence-space CI calculations, culminating in the publication of Hamiltonians for the and shells [[37], [38]]. Their approach relied on Brueckner’s matrix to treat short-range correlations induced by the free-space interaction, and employed the so-called hole-line expansion to second order [[39], [40], [41]]. After some initial successes, Barret, Kirson, and others demonstrated a lack of order-by-order convergence of this expansion [[42], [43], [44], [45], [46]], and Vary, Sauer and Wong found a disturbingly strong model-space dependence in intermediate-state summations, with larger model spaces actually degrading the agreement with experimental data [[47]]. Bogner et al. revisited this issue with the help of the Similarity Renormalization Group (SRG) [[48], [49]], and demonstrated that the matrix retains significant coupling between low- and high-momentum nodes of the underlying interaction [[33]], so the convergence issues are not surprising from a modern perspective. In the SRG and other modern RG approaches, low- and high-momentum physics are decoupled properly, and the resulting low-momentum interactions are indeed perturbative [[50], [33]]. For such interactions, results from finite-order many-body perturbation theory (MBPT) are in good agreement with non-perturbative results if the expansion is based on a Hartree-Fock reference state [[51], [52]].

Of course, low-momentum interactions are well-suited inputs not just for MBPT, but for all methods that work in truncated configuration spaces. The decoupling of low- and high-momentum modes of the interaction leads to a greatly improved convergence behavior, which in turn extends the range of nuclei a many-body method can be applied to. With SRG-evolved interactions, the No-Core Shell Model (NCSM) and related large-scale diagonalization methods can be extended into the lower shell [[53], [54], [55], [56]], and methods with systematic many-body truncations like Coupled Cluster (CC) are nowadays applied to nuclei as heavy as tin [[57], [58], [59]].

Instead of merely using the SRG as a tool to “pre-process” the nuclear interactions that are used as inputs for other many-body methods, we can turn it into a method for solving the many-body Schrödinger equation itself. This leads us to the so-called In-Medium SRG (IMSRG), which is the main focus of the present work [[60], [61], [62]]. In a nutshell, we want to use SRG-like flow equations to decouple physics at different excitation energy scales of the nucleus, and render the Hamiltonian matrix in configuration space block or band diagonal in the process. This can also be viewed as a re-organization of the many-body expansion, in which correlations that are described explicitly by the configuration space are absorbed into an RG-improved Hamiltonian. With an appropriately chosen decoupling strategy, it is even possible to extract eigenvalues and eigenstates of the nuclear Hamiltonian, and therefore, the IMSRG qualifies as an ab initio method for solving quantum many-body problems.

In this work, we will discuss two distinct implementations of the IMSRG ideas. The first is the so-called Multireference IMSRG (MR-IMSRG), which is designed for calculations of the ground-state properties of closed- and open-shell nuclei. Like most many-body approaches, it relies on the organization of the many-body basis in terms of a reference state and its excitations. Contrary to approaches like MBPT or CI, which employ Slater determinant reference states, the MR-IMSRG is built for arbitrary correlated reference states. This gives us the greatest possible flexibility in the description of correlations: Static correlations, e.g., due to intrinsic deformation, can be built into the reference state, while dynamic correlations due to the excitation of nucleon pairs, triples, etc. are described by the MR-IMSRG transformation.

The second approach uses the IMSRG to construct RG-improved Hamiltonians for nuclear valence CI calculations. These interactions are subsequently used as input for existing Shell model codes like NuShellX [[63]]. Such a combined approach gives us immediate access to a much larger number of observables than the MR-IMSRG, but it is limited by the computational effort of the CI part of the calculation.

The idea of using flow equations to solve quantum many-body problems was already discussed in Wegner’s initial work on the SRG [[49]] (also see [[64]] and references therein). In the solid-state physics literature, the approach is also known as continuous unitary transformation (CUT) theory, see [[65], [66], [67], [68], [69]]. When we discuss our decoupling strategies for the nuclear many-body problem, it will become evident that the IMSRG is related to CC [[70], [58]], canonical transformation theory (CT) [[71], [72], [73]], and the Irreducible (or Anti-Hermitian) Contracted Schrödinger Equation (ICSE) approach [[74], [75], [76], [77], [78], [79], [80]], and there is even some overlap with purely variational methods (see section 4.3). What sets the IMSRG apart from these methods is that the Hamiltonian instead of the wave function is at the center of attention, in the spirit of RG methodology. This seems to be a trivial distinction, but there are practical advantages of this viewpoint, e.g., the simultaneous decoupling of ground and excited states (see section 4.2), the avoidance of -representability issues [[81]], and more. Inspired in part by our work on the IMSRG in nuclear physics, Evangelista and co-workers have recently presented the Driven SRG for ab initio calculations in quantum chemistry, which implements IMSRG transformations in terms of inhomogeneous nonlinear equations rather than flow equations [[82], [83], [84], [85]].

#### Organization of this work:

Let us conclude our introduction of the IMSRG with a look ahead at the main body of this work. In section 2, we briefly review the essential concepts of the free-space SRG, and how it is used to dial the resolution scale of (nuclear) interactions and operators. To set up the IMSRG formalism, we first discuss normal-ordering techniques and Wick’s theorem for arbitrary reference states in section 3. This will be followed by the presentation of the MR-IMSRG flow equations in section 4, which also discusses the choice of decoupling scheme and generators. The single-reference IMSRG equations are obtained as a limit of the more general MR-IMSRG framework. In section 5, general features of (MR-)IMSRG flows are discussed, and section 6 reviews applications of the MR-IMSRG to the ground-state properties of closed- and open-shell nuclei. In section 7, we describe the derivation of nonempirical interactions for the nuclear valence-space CI approach. Salient points are summarized in section 8, and we look ahead at new developments. Expressions for products and commutators of normal-ordered operators are collected in A, and B recapitulates elements of the Hartree-Fock Bogoliubov theory and particle-number projection, which are used in the construction of reference states for the MR-IMSRG.

## 2 The Similarity Renormalization Group

### 2.1 General Concept

The Similarity Renormalization Group (SRG) was first formulated independently by Wegner [[49]] and Glazek and Wilson [[48]], for applications in condensed matter physics and light-front quantum field theory, respectively. The general concept of the method is to simplify the structure of the Hamiltonian in a suitable representation through the use of a continuous unitary transformation,

 H(s)=U(s)H(0)U†(s). (1)

Here, is the starting Hamiltonian and denotes the so-called flow parameter, which parameterizes the unitary transformation. Taking the derivative of equation (1) with respect to , we obtain the operator flow equation

 ddsH(s)=[η(s),H(s)], (2)

where the anti-Hermitian generator is related to by

 η(s)=dU(s)dsU†(s)=−η†(s). (3)

By rearranging this relation, we obtain a differential equation for whose formal solution is given by the path-or S-ordered exponential [[86], [87]]

 U(s)=Sexp∫s0ds′η(s′)≡∑n1n!∫s0ds1∫s0ds2…∫s0dsnS{η(s1)…η(sn)}. (4)

Here, the -ordering operator ensures that the flow parameters appearing in the integrands are always in descending order, . Since our continuous unitary transformation preserves the spectrum of the Hamiltonian and any other observable of interest, it is an example of a so-called isospectral flow, a class of transformations which has been studied extensively in the mathematics literature (see, e.g., [[88], [89], [90]]).

With a suitable choice of generator , we can smoothly transform the Hamiltonian to almost arbitrary structures as we integrate the flow equation (2) for . Inspired by the work of Brockett [[88]] on the so-called double-bracket flow, Wegner [[49]] proposed the generator

 η(s)≡[Hd(s),Hod(s)], (5)

which is constructed by splitting the Hamiltonian into suitably chosen diagonal () and offdiagonal () parts. It can be shown analytically that the generator (5) will monotonically suppress as the Hamiltonian is evolved via equation (2) (see, e.g., [[49], [64], [62]]). Note that the label diagonal does not need to mean strict diagonality here, but rather refers to a desired structure that the Hamiltonian will assume in the limit . By working in bases that are ordered by momenta or energies, the capability to impose structure on the Hamiltonian allows us to make an explicit connection with renormalization group (RG) ideas.

### 2.2 SRG Evolution of Nuclear Interactions

In figure 1, we show schematic examples of RG evolutions that are applied to nucleon-nucleon interactions in momentum-space representation. Figure 1(a) implements the RG as a decimation: The interaction is evolved to decreasing cutoff scales , and we end up with a low-momentum interaction that only has non-zero matrix elements between states with initial and final relative momenta [[31], [33]]. In contrast, figure 1(b) results from a continuous unitary transformation via the flow equation (2), using a Wegner-type generator built from the relative kinetic energy in the two-nucleon system:

 η(λ)≡[→k22μ,v(λ)]. (6)

Here, , and is the reduced nucleon mass. We have parametrized the transformation with , which has the dimensions of a momentum (in natural units). As suggested by figure 1(b), is a measure for the “width” of the band-diagonal Hamiltonian in momentum space, i.e., it controls the scale of momentum transfers between nucleons. Because

 |→q|=|→k′−→k|≲λ (7)

low- and high-lying momenta are decoupled in a proper RG sense as is decreased.

The decoupling of low- and high-lying momenta significantly improves the convergence properties of configuration-space based many-body methods, because it prevents the Hamiltonian from scattering nucleon pairs from low to high momentum states. Methods like the NCSM or the IMSRG discussed below yield converged results in much smaller many-body Hilbert spaces, which in turn makes it possible to apply these methods to heavier nuclei [[91], [53], [54], [56], [55], [61], [92], [62], [93], [94], [95], [57], [96], [97], [98], [99]]. However, this improvement comes at a cost, which is best illustrated by considering the Hamiltonian in a second-quantized form, assuming only a two-nucleon interaction for simplicity:

 Hint=Trel+V=14∑pqrs⟨pq|→k2122μ+v12|rs⟩a†pa†qasar. (8)

If we plug and into the commutators in equations (6) and (2), we obtain

 [a†ia†jalak,a†pa†qasar]=δlpa†ia†ja†qasarak+{a†a†a†aaa}−δlpδkqa†ia†jasar+{a†a†aa}, (9)

where the bracketed terms with suppressed indices schematically stand for additional two- and three-body operators. Thus, even if we start from a pure two-body interaction, the SRG flow will induce operators of higher rank, i.e., three-, four-, and in general up to -nucleon interactions. Of course, these induced interactions are only probed if we study an -nucleon system. If we truncate the SRG flow equations at the two-body level, we preserve the properties of the two-nucleon system, in particular phase shifts and the deuteron binding energy. A truncation at the three-body level ensures the invariance of observables in nuclei, e.g.  and ground-state energies, and so on. Truncations in the SRG flow equation cause a violation of unitarity that manifests as a (residual) dependence of many-body results on . By varying this parameter, the size of the missing contributions can be assessed (see, e.g., [[33], [100], [101], [91], [61], [55], [57], [99], [102]]).

State-of-the-art SRG evolutions of nuclear interactions are nowadays performed in the three-body system, using relative (Jacobi) harmonic oscillator [[100], [103], [54]], relative momentum plane wave [[101]], or momentum-space hypherspherical harmonics representations [[104]]. Pioneering work on implementing the SRG evolution in the lowest partial waves of the four-body system has been carried out by A. Calci and co-workers [[105]], again working in Jacobi HO representation.

Figure 2 shows the evolution of NN and matrix elements of a chiral NNLO interaction by Epelbaum, Glöckle, and Meißner [[106], [107]], with cutoffs 550/600 MeV. As discussed for our schematic example, both the and interactions become band diagonal and the SRG decouples the low- and high-momentum regimes as we evolve to lower values of . In figure 3, the same family of SRG-evolved interactions is used to calculate the ground-state energy of the triton, as a function of . If only the part of the chiral interaction is used as input, and the SRG generator and flowing Hamiltonian are truncated at the two-body level (curve ’-only’), the SRG evolution is not unitary in the three-body system. The energy exhibits a significant dependence on , on the order of 5–6%. If the flow equations are truncated at the three-body level instead, induced interactions are properly included and the unitarity of the transformation is restored (’-induced’): The energy does not change as is varied. Finally, the curve ’-full’ shows the result for a calculation with initial and forces that are consistently SRG-evolved at the three-body level. The triton ground-state energy is again invariant under the SRG flow, and closely reproduces the experimental value that is used as a constraint in the adjustment of the force’s low-energy constants (see, e.g., [[14], [15], [108]]).

The SRG flow equations force us to manipulate large sections (or the entirety) of the Hamiltonian’s spectrum in order to avoid basis truncation artifacts (also cf. [[56], [57]]) We may ask, then, if it might be possible to avoid the use of matrix representations entirely by solving the operator flow equation (2) directly in the algebra of operators. This is the strategy that we will explore in the following, which will ultimately lead us to the formulation of the In-Medium SRG. First, we have to lay some groundwork on normal ordering techniques and Wick’s theorem.

## 3 Normal Ordering and Wick’s Theorem for Arbitrary Reference States

### 3.1 References States and Many-Body Bases

To describe the structure and dynamics of an atomic nucleus of mass , we need to work in an -body Hilbert space2, and choose a suitable -body basis. Since we are dealing with a system of fermions, a straightforward choice are antisymmetrized product states, or Slater determinants. Introducing fermionic creation and annihilation operators and that satisfy the canonical anticommutation relations

 {a†i,a†j}={ai,aj}=0,{a†i,aj}=δij, (10)

we can write a generic -particle Slater determinant as

 |Φ⟩=A∏k=1a†ik|vac⟩, (11)

where refers to the particle vacuum. Here, the indices run over a suitably chosen single-particle basis, e.g., spatially localized orbitals if we plan to describe a finite system like a nucleus. A complete basis for the many-body Hilbert space can be obtained by distributing nucleons over the available single-particle states in all possible ways.

Of course, not all of the states in this naively chosen basis are created equal. As alluded to in section 2, nuclear interactions and the nucleus itself have characteristic energy or momentum scales. The ground state and low-lying excitation spectrum of an -body nucleus is typically dominated by excitations of particles in the vicinity of its Fermi energy. Thus, we can find a Slater determinant that is a fair approximation to the nuclear ground state, and use it as a reference state for the construction and organization of our many-body basis. Slater determinants that are variationally optimized through a Hartree-Fock (HF) procedure have been shown to be reasonable reference states for interactions with low resolution scales around (see, e.g., Refs. [[33], [115], [53], [58], [62], [51]] and references therein), allowing post-HF methods like MBPT, CC, or the IMSRG discussed below to converge rapidly to the exact FCI result. Starting from such a HF reference state , we can obtain a basis consisting of the state itself and up to -particle, -hole () excitations:

 |ΦHF⟩,a†p1ah1|ΦHF⟩,…,a†p1…a†pAahA…ah1|ΦHF⟩. (12)

Here, indices and run over all one-body basis states with energies above (particle states) and below the Fermi level (hole states), respectively.

Many-body bases built from such a single Slater determinant and its particle-hole excitations work best for systems with large gaps in the single-particle spectrum, e.g., closed-shell nuclei. If the gap is small, particle-hole excited basis states can be near-degenerate with the reference state, which usually results in spontaneous symmetry breaking and strong configuration mixing. At best, these phenomena impede the convergence of a many-body calculation by forcing us to use model spaces that contain excitations with large , e.g., in a CI framework. At worst, the behavior of a truncated many-body method like IMSRG or CC may be completely uncontrolled. We want to overcome these problems by building correlations from configuration mixing into the reference state, and constructing a basis of generalized excitations on top of this state. A key element of such an approach are generalized normal ordering techniques.

### 3.2 Normal-Ordered Operators and Wick Contractions

In Ref. [[116]], Kutzelnigg and Mukherjee developed a generalized normal ordering for arbitrary reference states. Here, we present the essential elements of their discussion that we will need in the following, but use the slightly different notation of Kong et al. [[117]].

First, we introduce a pseudo-tensorial notation for strings of creation and annihilation operators, to facilitate book-keeping and make the formalism more compact. A particle-number conserving product of creators and annihilators each is written as

 Ai1…ikj1…jk≡a†i1…a†ikajk…aj1. (13)

We do not consider particle-number changing operators in the present work, because they cause ambiguities in the contraction and sign rules for the operators that are defined in the following. The anticommutation relations imply

 AP(i1…ik)P′(j1…jk)=(−1)π(P)+π(P′)Ai1…ikj1…jk, (14)

where indicates the parity (or signature) of a permutation . A general -body operator in second quantization can now be written in terms of the basis operators as

 O(k)=1(k!)2∑i1…ikj1…jkoi1…ikj1…jkAi1…ikj1…jk, (15)

where we assume that the coefficients are antisymmetrized, and therefore also obey equation (14) under index permutations.

Next, we introduce irreducible -body density matrices . In the one-body case, we have the usual density matrix

 λij≡⟨Φ|Aij|Φ⟩, (16)

and for future use, we also define

 ξij≡λij−δij. (17)

Up to a factor that unifies the sign rules for one-body contractions presented below, is simply the generalization of the hole density matrix for a correlated state. In the natural orbital basis, i.e., the eigenbasis of , both one-body density matrices are diagonal:

 λij=niδij,ξij=−¯niδij≡−(1−ni)δij. (18)

The fractional occupation numbers are the eigenvalues of .

For , we denote full density matrices by

 ρi1…ikj1…jk =⟨Φ|Ai1…ikj1…jk|Φ⟩, (19)

and define

 λijkl ≡ρijkl−A{λikλjl}, (20) λijklmn ≡ρijklmn−A{λilλjkmn}−A{λilλjmλkn}, (21)

etc., where fully antisymmetrizes the indices of the expression within the brackets, e.g.,

 A{λikλjl}=λikλjl−λilλjk. (22)

From equation (20), it is easy to see that encodes the two-nucleon correlation content of the reference state . If the reference state is a Slater determinant, i.e., an independent-particle state, the full two-body density matrix factorizes, and vanishes:

 λijkl=ρijkl−A{λikλjl}=λikλjl−λikλjl−(λikλjl−λikλjl)=0. (23)

Equation (21) shows that is constructed by subtracting contributions from three independent particles as well as two correlated nucleons in the presence of an independent spectator particle from the full three-body density matrix, and therefore encodes the genuine three-nucleon correlations. This construction and interpretation generalizes to irreducible density matrices of rank .

Now we consider the expansion of a (number-conserving) string of creation and annihilation operators in normal-ordered components. First, we define a normal-ordered one-body operator by subtracting from a given one-body operator its expectation value in the reference state:

 :Aab:≡Aab−⟨Φ|Aab|Φ⟩=Aab−λab. (24)

This implies that the expectation value of the normal-ordered operator in the reference state vanishes by construction:

 ⟨Φ|:Aab:|Φ⟩=0. (25)

For a two-body operator, we have the expansion

As a consequence of equation (14), the sign of each term is determined by the product of the parities of the permutations that map upper and lower indices to their ordering in the initial operator. Except for the last term, this expression looks like the result for the regular normal ordering, with pairwise contractions of indices giving rise to one-body density matrices. The last term, a contraction of four indices, appears because we are dealing with an arbitrary, correlated reference state here.

Taking the expectation value of equation (26) and using equation (25), we obtain

and see that

 ⟨Φ|:Aabcd:|Φ⟩=0. (28)

The normal ordering procedure can be extended in analogy to the one- and two-body cases, e.g.,

yielding normal-ordered operators of arbitrary rank that satisfy

 ⟨Φ|:Ai1…ikj1…jk:|Φ⟩=0. (30)

Finally, a generalized Wick’s theorem for arbitrary reference states can be formulated: Any product of two normal-ordered operators can be expanded in a sum of normal-ordered terms, with Wick contractions and operators containing at least one index from each of the original operators. For example, the basic contractions appearing in the expansion of a product of normal-ordered two-body operators are (notice the signs)

 :Aa\framebox$b$cd::Aij\framebox$k$l: =−λ\framebox$b$\framebox$k$:Aaijcdl:, (31) :Aab\framebox$c$d::Ai\framebox$j$kl: =−ξ\framebox$j$\framebox$c$:Abiadkl:, (32) :A\framebox$ab$cd::Aij\framebox$kl$: =+λ\framebox$ab$\framebox$kl$:Aijcd:, (33) :Aa\framebox$b$cd::A\framebox$i$j\framebox$kl$: =−λ\framebox$ib$\framebox$kl$:Aajcd:, (34) :Aab\framebox$c$d::A\framebox$ij$\framebox$k$l: =−λ\framebox$ij$\framebox$ck$:Aabdl:, (35) :A\framebox$ab$c\framebox$d$::A\framebox$i$j\framebox$kl$: =−λ\framebox$abi$\framebox$dkl$:Ajc:, (36) :A\framebox$ab$\framebox$cd$::A\framebox$ij$\framebox$kl$: =+λ\framebox$abij$\framebox$cdkl$. (37)

Only the first two contraction types, equations (31) and (32), appear in the regular Wick’s theorem for uncorrelated reference states. The additional contractions (33)–(37) increase the number of terms when we expand operator products using the generalized Wick’s theorem — examples are shown in appendix A. Fortunately, we will see in section 4 that the overall increase in complexity is manageable.

### 3.3 Normal-Ordered Hamiltonian and Normal-Ordered Two-Body Approximation

We conclude this section by applying the generalized normal ordering to an intrinsic nuclear -body Hamiltonian containing both and interactions, which will be relevant for the applications discussed later in this work. Let

 H=(1−1A)T(1)+1AT(2)+V(2)+V(3), (38)

where

 T(1)≡∑i→p2i2m,T(2)≡−∑i

(see, e.g., [[118]]). Choosing a generic correlated reference state , we rewrite the Hamiltonian as

 H =E+∑ijfij:Aij:+14∑ijklΓijkl:Aijkl: =+136∑ijklmnWijklmn:Aijklmn:, (40)

where the labels have been chosen for historical reasons. The individual normal-ordered contributions in equation (3.3) are given by

 E =(1−1A)∑abtabλab+14∑abcd(1Atabcd+vabcd)ρabcd =+136∑abcdefvabcdefρabcdef, (41) fij =(1−1A)tij+∑ab(1Atiajb+viajb)λab =+14∑abcdviabjcdρabcd, (42) Γijkl =1Atijkl+vijkl+∑abvijaklbλab, (43) Wijklmn =vijklmn. (44)

Here, we use the full density matrices for compactness, but it is easy to express equations (41)–(43) completely in terms of irreducible density matrices by using equations (20) and (21). Note that the normal-ordered zero-, one-, and two-body parts of the Hamiltonian all contain in-medium contributions from the free-space interaction. It has been shown empirically that the omission of the normal-ordered three-body piece of the Hamiltonian causes a deviation of merely 1–2% in ground-state and (absolute) excited state energies of light and medium-mass nuclei [[119], [91], [94], [120], [121]]. This normal-ordered two-body approximation (NO2B) to the Hamiltonian is widely used nowadays, because it provides an efficient means to account for force effects in nuclear many-body calculations without incurring the computational expense of explicitly treating three-body operators. In the next section, we will also see that the NO2B approximation meshes in a natural way with the framework of the MR-IMSRG, which makes it especially appealing for our purposes.

## 4 The Multireference In-Medium Similarity Renormalization Group

### 4.1 MR-IMSRG Flow Equations

We are now ready to formulate the MR-IMSRG flow equations by applying the tools discussed in the previous section to the operator differential equation (2). We express all operators in terms of normal-ordered strings of creation and annihilation operators. As discussed in section 2, each evaluation of the commutator on the right-hand side will induce operators of higher rank,

 [:Aabcd:,:Aijkl:]=δci:Aabjdkl:+…, (45)

and we would need to include operators up to rank if we want the MR-IMSRG flow to be unitary in an -body system, which is not feasible. However, in contrast to equation (9), we are now working with normal-ordered operators whose in-medium contributions have been absorbed into terms of lower rank. Consequently, we expect the induced operators to be much weaker than in the free-space SRG case. The empirical success of the NO2B approximation discussed at the end of section 3 certainly seems to justify this expectation in the case of nuclear + Hamiltonians with low resolution scales.

Following this line of reasoning further, we choose to truncate all flowing operators at a given particle rank in order to obtain a closed system of flow equations. For , we demand that for all values of the flow parameter

 η(s) ≈η(1)(s)+η(2)(s), (46) H(s) ≈E(s)+f(s)+Γ(s), (47) ddsH(s) ≈ddsE(s)+ddsf(s)+ddsΓ(s). (48)

This is the so-called MR-IMSRG(2) truncation, which has been our primary workhorse in past applications [[60], [122], [61], [55], [92], [123], [62]]. It is the basis for all results in this work. We note that the MR-IMSRG(2) at this specific level of truncation is a cousin to a variety of other truncated many-body schemes, e.g., Canonical Transformation theory with Singles and Doubles excitations (CTSD) [[71], [72], [73]], the two-body Antisymmetrized or Irreducible Contracted Schrödinger Equation approach (ICSE(2)) [[74], [75], [76], [77], [78], [79], [80]] and of course CCSD (Coupled Cluster with Singles and Doubles) (see, e.g., [[70]]), although the latter is based on non-unitary similarity transformations.

Plugging equations (46)–(48) into the operator flow equation (2) and organizing contributions by particle rank, we obtain the system of MR-IMSRG(2) flow equations:

 dEds Missing or unrecognized delimiter for \right =+14∑abcd(ddsΓabcd)λabcd+14∑abcdklm(ηabcdΓklam−Γabcdηklam)λbklcdm, (49) ddsfij =∑a(ηiafaj−fiaηaj)+∑ab(ηabΓbiaj−fabηbiaj)(na−nb) =+12∑abc(ηiabcΓbcja−Γiabcηbcja)(na¯nb¯nc+¯nanbnc) =+14∑abcde(ηiabcΓdeja−Γiabcηdeja)λdebc+∑abcde(ηiabcΓbejd−Γiabcηbejd)λaecd =−12∑abcde(ηiajbΓcdae−Γiajbηcdae)λcdbe+12∑abcde(ηiajbΓbcde−Γiajbηbcde)λacde, (50) ddsΓijkl =∑a(ηiaΓajkl+ηjaΓiakl−ηakΓijal−ηalΓijka−fiaηajkl−fjaηiakl+fakηijal+falηijka) =+12∑ab(ηijabΓabkl−Γijabηabkl)(1−na−nb) =+∑ab(na−nb)((ηiakbΓjbla−Γiakbηjbla)−(ηjakbΓibla−Γjakbηibla)). (51)

All single-particle indices and occupation numbers (cf. section 3) refer to natural orbitals, and the -dependence has been suppressed for brevity. The single-reference limit is readily obtained by setting the irreducible density matrices and to zero in the previous expressions.

We solve the flow equations (49)–(51) by integrating from to , using the components of the normal-ordered input Hamiltonian (equations (41)–(43)) as initial values. In this process, the flow equations will re-shuffle the correlations in the -body system, generating a highly nonperturbative resummation of the many-body expansion (see section 5.2 for numerical examples).

To interpret the multireference flow equations, we associate the fractional occupation numbers and with particle- and hole-like states, respectively (cf. equation (18)), and note that

 1−na−nb=¯na¯nb−nanb, (52) na−nb=na¯nb−¯nanb. (53)

For the typical ansätze that we use for (see section 4.3), the generator is proportional to the (offdiagonal) Hamiltonian, and we see the first two terms of the zero-body flow equation have the structure of second-order energy corrections, but evaluated for the flowing Hamiltonian . Furthermore, we recognize that the second and third lines of equation (51) have the structure of ladder (particle-particle / hole-hole) and ring (particle-hole) skeleton diagrams, respectively. They generate ladder and ring summations in the limit , but also ring-ladder interference diagrams with rich topologies that go far beyond traditional re-summation methods [[41], [124], [125]]. A detailed perturbative analysis is presented in [[62]].

For general reference states, the MR-IMSRG flow equations also include couplings to correlated pairs and triples of nucleons through the irreducible density matrices and . It is noteworthy that the MR-IMSRG(2) flow equations do not depend on or nonlinear powers of . While such contractions appear in the products of normal-ordered two-body operators, they cancel in the commutators (see A). This ensures that the MR-IMSRG only sums so-called connected many-body diagrams (i.e., diagrams which do not contain intermediate insertions of the reference state) [[124], [70]].

Let us conclude this section by briefly considering the numerical implementation of the MR-IMSRG(2) scheme. The computational effort is dominated by the two-body flow equation (51), which naively requires operations, where denotes the single-particle basis size. This puts the MR-IMSRG(2) in the same category as its aforementioned “cousins” CCSD [[70], [58]], CTSD [[71], [72], [73]], and ICSE(2) [[74], [75], [76], [77], [78], [79], [80]], as well as the Self-Consistent Green’s Function Approaches (SCGF) [[126], [127], [128], [97], [98]]. Fortunately, the flow equations can be expressed in terms of matrix products and traces, allowing us to use optimized linear algebra libraries provided by high-performance computing vendors.

Moreover, we can reduce the computational cost in the single-reference case by distinguishing particle and hole states, because the number of hole states is typically much smaller than the number of particle states . The best scaling we can achieve in the IMSRG(2) depends on the choice of generator (see section 4.3). If the one- and two-body parts of the generator only consist of ph and pphhh type matrix elements (and their Hermitian conjugates), respectively, the scaling is reduced to , which matches the cost of solving the CCSD amplitude equations.

### 4.2 Decoupling Strategy

Having set up the flow equations in the previous section, we now need to specify our decoupling strategy, i.e., how we split the Hamiltonian into diagonal parts that we want to keep, and offdiagonal parts that we want to suppress through the MR-IMSRG evolution (cf. section 2). To do this, we refer to the matrix representation of the Hamiltonian in a given -body basis, which is shown schematically for single- and multireference cases in figures 4 and 5, respectively. We stress that we do not actually need to construct the Hamiltonian matrix in this representation.

#### IMSRG Decoupling in the Single-Reference Case

Let us consider the simpler single-reference case first. We choose a Slater determinant reference and construct a basis by considering all possible particle-hole excitations (cf. section 3):

 |Φ⟩,:Aph:|Φ⟩,:App′hh′:|Φ⟩,…. (54)

Note that because contractions of particle and hole indices vanish by construction. Using Wick’s theorem, it is easy to see that the particle-hole excited Slater determinants are orthogonal to the reference state as well as each other. In the Hilbert space spanned by this basis, the matrix representation of our initial Hamiltonian in the NO2B approximation (or any two-body operator) has the structure shown in the left panel of figure 4, i.e., it is band-diagonal, and can at most couple and excitations.

We now have to split the Hamiltonian into appropriate diagonal and offdiagonal parts on the operator level, which is a non-trivial task (see, e.g., the extensive discussion in Refs. [[129], [130], [131]] ). Using a broad definition of diagonality is ill-advised because we must avoid to induce strong in-medium interactions to maintain the validity of the IMSRG(2) truncation. We choose what we call a minimal decoupling scheme that aims to decouple the one-dimensional block spanned by the reference state from all particle-hole excitations, as shown in the right panel of figure 4.

If we could implement this decoupling without approximations, we would extract a single eigenvalue and eigenstate of the many-body Hamiltonian for the nucleus of interest in the limit . The eigenvalue would simply be given by the zero-body piece of , while the eigenstate is obtained by applying the unitary IMSRG transformation to the reference state, . In practice, we end up with an approximate eigenvalue and mapping.

An important caveat is that we cannot guarantee a priori that we will target the true interacting ground state and its energy eigenvalue in this way. Empirically, the IMSRG flow is found to connect the reference state to the eigenstate with which it has the highest overlap. In single-reference scenarios, a HF Slater determinant will typically have the highest overlap with the exact ground state because it minimizes both the absolute energy and the correlation energy, the latter being due to admixtures from particle-hole excitations. In the multireference case, we have found examples where the MR-IMSRG flow targets excited states, as discussed in sections 6.3 and 7.4.

Analyzing the matrix elements between the reference state and its excitations with the help of Wick’s theorem, we first see that the Hamiltonian couples the ph block to ph excitations through the matrix elements

 ⟨Φ|H:Aph:|Φ⟩ =E⟨Φ|:Aph:|Φ⟩+∑ijfij⟨Φ|:Aij::Aph:|Φ⟩+14∑ijklΓijkl⟨Φ|:Aijkl::Aph:|Φ⟩ =∑ijfijδihδpjni¯nj=fhp (55)

and their Hermitian conjugates. The contributions from the zero-body and two-body pieces of the Hamiltonian vanish because they are expectation values of normal-ordered operators in the reference state (cf. equation (28)). Likewise, the and blocks are coupled by the matrix elements

 ⟨Φ|H:App′hh′:|Φ⟩=Γhh′pp′ (56)

and their conjugates. It is precisely these two-body matrix elements that couple and states and generate the outermost side diagonals of the Hamiltonian matrix. This suggests that we can transform the Hamiltonian to the shape shown in the top right panel of figure 4 by defining its offdiagonal part as

 Hod≡∑phfph:Aph:+14∑pp′hh′Γpp′hh′:App′hh′:+H.c.. (57)

In section 5, we will show that the IMSRG flow does indeed exponentially suppress the matrix elements of and achieve the desired decoupling in the limit .

#### Variational Derivation of Minimal Decoupling

Our minimal decoupling scheme is of course very reminiscent of the strategy followed in Coupled Cluster approaches [[70], [58]], except that we specifically use a unitary transformation instead of a general similarity transformation. It is also appealing for a different reason: As we will discuss now, it can be derived from a variational approach, tying the seemingly unrelated ideas of energy minimization and renormalization in the many-body system together. To this end, we consider the energy expectation value of the final IMSRG evolved Hamiltonian,

 ¯¯¯¯H≡H(∞), (58)

in the reference state (which is assumed to be normalized):

 E=⟨Φ|¯¯¯¯H|Φ⟩. (59)

Next, we introduce a unitary variation, which we can choose to apply either to the reference state ,

 |Φ⟩→eZ|Φ⟩,Z†=−Z, (60)

or, equivalently, to the Hamiltonian:

 Missing dimension or its units for \mskip (61)

The variation of the energy is

 Missing dimension or its units for \mskip (62)

with a suitable operator norm . The first term obviously vanishes, as does the commutator of with the energy. Thus, the energy is stationary if

 δE=⟨Φ|[¯¯¯¯H,Z]|Φ⟩=0. (63)

Expanding

 Z=∑phZph:Aph:+14∑pp′hh′Zpp′hh′:App′hh′:+H.c.+…, (64)

and using the independence of the expansion coefficients (save for the unitarity conditions), we obtain the system of equations

 ⟨Φ|[¯¯¯¯H,:Aph:]|Φ⟩ =0, (65) ⟨Φ|[¯¯¯¯H,:Ahp:]|Φ⟩ =0, (66) ⟨Φ|[¯¯¯¯H,:App′hh′:]|Φ⟩ =0, (67) ⟨Φ|[¯¯¯¯H,:Ahh′pp′:]|Φ⟩ =0, (68) …

which are special cases of the so-called irreducible Brillouin conditions (IBCs) [[76], [77], [78], [79]]. Writing out the commutator in the first equation, we obtain

 ⟨Φ|[¯¯¯¯H,:Aph:]|Φ⟩=⟨Φ|¯¯¯¯H:Aph:|Φ⟩−⟨Φ|:Aph:¯¯¯¯H|Φ⟩=⟨Φ|¯¯¯¯H:Aph:|Φ⟩=0, (69)

where the second term vanishes because it is proportional to . The remaining equations can be evaluated analogously, and we find that the stationarity conditions are satisfied if the IMSRG evolved Hamiltonian no longer couples the reference state and its particle-hole excitations, as discussed above. This connection between the decoupling conditions and the stationarity conditions of an energy functional will prove useful in the multireference case.

#### MR-IMSRG Decoupling for Correlated Reference States

In the multireference case, we choose a suitable correlated reference state, and construct its excitations by applying all possible one- and two-body operators:

 |Φ⟩,:Aij:|Φ⟩,:Aijkl:|Φ⟩,…. (70)

The properties of the normal ordering ensure that the excited states are orthogonal to the reference state, but they are in general not orthogonal to each other: for instance,

 ⟨Φ|