# 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,

(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

(2) |

where the anti-Hermitian generator is related to by

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

(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

(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:

(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

(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:

(8) |

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

(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 space^{2}*antisymmetrized
product states*, or *Slater determinants*. Introducing
fermionic creation and annihilation operators and that
satisfy the canonical anticommutation relations

(10) |

we can write a generic -particle Slater determinant as

(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:

(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

(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

(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

(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

(16) |

and for future use, we also define

(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:

(18) |

The fractional occupation numbers are the eigenvalues of .

For , we denote full density matrices by

(19) |

and define

(20) | ||||

(21) |

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

(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:

(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:

(24) |

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

(25) |

For a two-body operator, we have the expansion

(26) |

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

(27) |

and see that

(28) |

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

(29) |

yielding normal-ordered operators of arbitrary rank that satisfy

(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)

(31) | ||||

(32) | ||||

(33) | ||||

(34) | ||||

(35) | ||||

(36) | ||||

(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

(38) |

where

(39) |

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

(40) |

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

(41) | ||||

(42) | ||||

(43) | ||||

(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,

(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

(46) | ||||

(47) | ||||

(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:

(49) | ||||

(50) | ||||

(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

(52) | |||

(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):

(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

(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

(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

(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,

(58) |

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

(59) |

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

(60) |

or, equivalently, to the Hamiltonian:

(61) |

The variation of the energy is

(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

(63) |

Expanding

(64) |

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

(65) | ||||

(66) | ||||

(67) | ||||

(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

(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:

(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,